TailBound-exploration-FC-MNIST.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-20T20:13:59.673640Z",
     "start_time": "2019-01-20T20:13:57.997310Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline\n",
    "from experiment_train import *\n",
    "from experiment_datasets import *\n",
    "from helpers import *\n",
    "from tqdm import tqdm\n",
    "import matplotlib"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-20T20:13:59.823133Z",
     "start_time": "2019-01-20T20:13:59.812310Z"
    }
   },
   "source": [
    "Problem: Taylor is WAY OFF for $np>1$...\n",
    "\n",
    "Hypothesis: need to add second term to the mean as only MEAN is off and variance is fine.\n",
    "Need to calculate $L''_{ii}$, diagonal of the Hessian\n",
    "\n",
    "A bit better (for second layer) with it...."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 262,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-20T20:13:59.823133Z",
     "start_time": "2019-01-20T20:13:59.812310Z"
    }
   },
   "outputs": [],
   "source": [
    "# tf._hessian_vector_product -- use for VarDelta regularization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-20T20:13:59.823133Z",
     "start_time": "2019-01-20T20:13:59.812310Z"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 363,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-20T20:13:59.823133Z",
     "start_time": "2019-01-20T20:13:59.812310Z"
    }
   },
   "outputs": [],
   "source": [
    "# network configuration\n",
    "N = [100, 200]\n",
    "\n",
    "def get_p_arr(p):\n",
    "    \"\"\" p array with failure on the first layer \"\"\"\n",
    "    return [0, 0, p]  + [0] * (len(N) - 1)\n",
    "\n",
    "# prob. of failure at inference\n",
    "p_inf = 1e-4\n",
    "p_inference = get_p_arr(p_inf)\n",
    "\n",
    "# prob. of failure at train\n",
    "p_trains = np.linspace(0, 1.2 * p_inf, 10)\n",
    "\n",
    "# Lips. coeff\n",
    "KLips = 1.\n",
    "\n",
    "# activation function\n",
    "activation = 'sigmoid'\n",
    "\n",
    "# number of repeats for error estimation\n",
    "# see below for explanation\n",
    "repetitions = 6\n",
    "\n",
    "# number of runs for each p\n",
    "runs = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 364,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "_________________________________________________________________\n",
      "Layer (type)                 Output Shape              Param #   \n",
      "=================================================================\n",
      "dense_43 (Dense)             (None, 100)               78500     \n",
      "_________________________________________________________________\n",
      "dense_44 (Dense)             (None, 200)               20200     \n",
      "_________________________________________________________________\n",
      "dense_45 (Dense)             (None, 10)                2010      \n",
      "=================================================================\n",
      "Total params: 100,710\n",
      "Trainable params: 100,710\n",
      "Non-trainable params: 0\n",
      "_________________________________________________________________\n",
      "Train on 60000 samples, validate on 10000 samples\n",
      "Epoch 1/10\n",
      "60000/60000 [==============================] - 5s 87us/step - loss: 97.9709 - categorical_accuracy: 0.7881 - mean_squared_error: 0.0492 - mean_absolute_error: 0.1479 - val_loss: 29.3187 - val_categorical_accuracy: 0.8936 - val_mean_squared_error: 0.0325 - val_mean_absolute_error: 0.1341\n",
      "Epoch 2/10\n",
      "60000/60000 [==============================] - 3s 55us/step - loss: 27.6892 - categorical_accuracy: 0.9074 - mean_squared_error: 0.0259 - mean_absolute_error: 0.1119 - val_loss: 26.5174 - val_categorical_accuracy: 0.9175 - val_mean_squared_error: 0.0210 - val_mean_absolute_error: 0.0981\n",
      "Epoch 3/10\n",
      "60000/60000 [==============================] - 3s 51us/step - loss: 25.8253 - categorical_accuracy: 0.9216 - mean_squared_error: 0.0241 - mean_absolute_error: 0.1118 - val_loss: 25.2323 - val_categorical_accuracy: 0.9289 - val_mean_squared_error: 0.0220 - val_mean_absolute_error: 0.1031\n",
      "Epoch 4/10\n",
      "60000/60000 [==============================] - 3s 50us/step - loss: 24.7458 - categorical_accuracy: 0.9305 - mean_squared_error: 0.0226 - mean_absolute_error: 0.1095 - val_loss: 24.3196 - val_categorical_accuracy: 0.9314 - val_mean_squared_error: 0.0263 - val_mean_absolute_error: 0.1239\n",
      "Epoch 5/10\n",
      "60000/60000 [==============================] - 3s 50us/step - loss: 24.0036 - categorical_accuracy: 0.9347 - mean_squared_error: 0.0228 - mean_absolute_error: 0.1130 - val_loss: 23.7626 - val_categorical_accuracy: 0.9341 - val_mean_squared_error: 0.0214 - val_mean_absolute_error: 0.1040\n",
      "Epoch 6/10\n",
      "60000/60000 [==============================] - 3s 54us/step - loss: 23.4680 - categorical_accuracy: 0.9403 - mean_squared_error: 0.0230 - mean_absolute_error: 0.1157 - val_loss: 23.3776 - val_categorical_accuracy: 0.9410 - val_mean_squared_error: 0.0219 - val_mean_absolute_error: 0.1117\n",
      "Epoch 7/10\n",
      "60000/60000 [==============================] - 3s 55us/step - loss: 23.2128 - categorical_accuracy: 0.9418 - mean_squared_error: 0.0265 - mean_absolute_error: 0.1296 - val_loss: 23.0739 - val_categorical_accuracy: 0.9466 - val_mean_squared_error: 0.0230 - val_mean_absolute_error: 0.1184\n",
      "Epoch 8/10\n",
      "60000/60000 [==============================] - 3s 55us/step - loss: 23.0150 - categorical_accuracy: 0.9444 - mean_squared_error: 0.0268 - mean_absolute_error: 0.1300 - val_loss: 22.7022 - val_categorical_accuracy: 0.9435 - val_mean_squared_error: 0.0271 - val_mean_absolute_error: 0.1263\n",
      "Epoch 9/10\n",
      "60000/60000 [==============================] - 3s 54us/step - loss: 22.9175 - categorical_accuracy: 0.9464 - mean_squared_error: 0.0285 - mean_absolute_error: 0.1362 - val_loss: 22.8396 - val_categorical_accuracy: 0.9477 - val_mean_squared_error: 0.0299 - val_mean_absolute_error: 0.1464\n",
      "Epoch 10/10\n",
      "60000/60000 [==============================] - 3s 55us/step - loss: 22.7825 - categorical_accuracy: 0.9505 - mean_squared_error: 0.0235 - mean_absolute_error: 0.1213 - val_loss: 22.7722 - val_categorical_accuracy: 0.9473 - val_mean_squared_error: 0.0239 - val_mean_absolute_error: 0.1205\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VPW9//HXJztZCUkIS9gUkD0gYVFUvFAUtWJdKC70XqxLbQWXto/+sHpb69La1nrt4rXFiopSEW1V7NW6orY1iQnKDrJKVmBCyEr2+fz+OJMwxEAGMmGSzOf5eOSRmTPnnPmcEd9z8j3f8/2KqmKMMSY4hAS6AGOMMaePhb4xxgQRC31jjAkiFvrGGBNELPSNMSaIWOgbY0wQsdA3xpggYqFvjDFBxELfGGOCSFigC2gtOTlZhw4dGugyjDGmW1m3bl2Jqqa0t16XC/2hQ4eSm5sb6DKMMaZbEZF9vqxnzTvGGBNELPSNMSaIWOgbY0wQsdA3xpggYqFvjDFBxELfGGOCiIW+McYEkS7XT98YY4JGUyOU7YNDu+HQLgjvBRk3dupbWugbY0xnUoWqg06oH9rp+e0J+dK94G44um7aVAt9Y4zpFuoqjw30Eq+Ar688ul5oJPQ5A1LOglGXQdJwz88IiO7T6WVa6BtjglqTW6moaSAkRAgLEUI9P2Ehgogcu3Jjvac5plWoH9oFVfu9VhToPcgJ80FTnUBPOtN5npAGIaGn9Ri9WegbY4LC4ep69pRUsdtVzR5XNXtLqtjjqmbfoSPUN7m91lT6UcqwkP2cIcWcGbKfM0KKGUYxA+UgYRxd9zDx5IcMpDBkLIWRF1EcOpDisIEcDB+IOzSS0GohtEYIKwrxfKkcJDTERagIoaGeLxnxfMmECoP7xPDdC8/s1M/BQt+YIKCq7Dt0hE2F5WwqLGfngUrie4XTP6EX/ROiPD+96N87iqSYiK+e4XYTtQ1N5JUeYY/LCfe9JdXscVWxp6SasiP1xFFDglSRHFLNyPhGrolrYMiIOvpJGfFH9hFfvY+EI/sId9e27LMhJJLSqMEcihzDJ5EXczBiMAcj0jgYlkZlSBxNbjdNCk1uN41NiluVWLfS5PlpdCs1DU00uhW353mT233M683LxwyIt9A3xpwcVSW/tIZNheVsLCxjc2E5mwrKqahtBCAiNIQzUmLY5apif3kxDU16zPYRoSH0a/kiiKJ/715HvxQ8y/oE8ItBmxo54DpAYVERBw8Uc6jkIFVlLmorSqDmMAlUkSDVjKaKC8KOkBRSTYJU06tXJSHadHRHNZ4fAAmFxCHQfzgkzT7aFJM0gvC4/qSGhJAaiIPtBBb6xnRjqkrB4Ro2F5azsbDc+V1QTnmN0yMkPFQY1S+eyyYMYEJaAuMHJjAyNY6IMOcWHbdbOVRdT3F5DcXltRSX1VBcUUtxWS3F5TXk7jvMgU1tfDGEhRz7F0Krvxb6J/QiMTr8xF8M7iaodoK65ae2rOVxfdUhjpSVUFd5CD1SSmh9OVGN5cRpNf2Afm3tMwzqw+PRqN6ExfQhNCYNonpDr8QT/0T3gdBw//xH6eIs9I3pJlSVovJaNhWUs6mwjE2FFWwqKOPwESfgw0KEs/rFccm4foxPS2DCwN6M7BdLZNjxLxqGhAgpcZGkxEUyIa3tddxupaS6zvNF4HwZ7C+vpai8lv3lNXy6t5QDFbU0up0vhl7UkiLlpIWWMyKmmmGR1aSFl5MaUk4fPUx8Yym96koIrT2EqLvN92wihEqNoUxjKSeGco2lIWIExCYSEZtMTO9kEvr0JTmlH4lJfZHoPk54R8YTEWqxdiL26RjTBakq+yuaA97zU1DOoep6AEJDhJGpccwZk8r4tN5MGJjAWf3iiAr3f6+QkBChb0wEfUOqSA8vgaj9EHMQ4vZD74OQtB+tPEBTxX6k+iChDVVHN65zfhoJxaUJHNTebNXeuHQcLhJxaQJlGksZsZRpDPRKJDEplf4pyQzrG8+w5BjOTIlhbFL0Cb+8jO98Cn0RmQv8FggF/qyqj7R6fQiwHEgBSoGFqlrg9Xo8sBV4TVUX+6l2Y3qMA56A926iKamqAyBEYGRqHP8xqm9LE83o/vH+CfiGWqg64Nw8VLUfKvcffVx18Ojz6oPgbvzq9hFxEJeKxKYSNnAixKZCXCrE9oPYvhDXD2JTCevVh74IIVV1uMtqaCyvpaa8lqaaBib2iW4J997RER0/JnNC7Ya+iIQCTwBzgAIgR0TWqOpWr9UeBVao6nMiMgv4BfAtr9cfBD72X9nGdF8HK2tbgr3598HKowE/vG8sM0emMH5gPOPTejOmfzy9IloFvLsJasqgvsq5KaiuyrkBqM7zvN77dxvLaiucIK8t/2qBEgIxKU5ox/aD1HGeIPf8xHkCPTYVImJ8Pu5QIDU+itT4qA58eqajfDnTnwrsUtU9ACKyCrgC58y92Rjg+57Ha4HXml8QkclAKvAPIMMPNRvTbTR3lczac4jsvaVs2p1PXaWLWGqJkxqGJ8DivsIZw5UhsU2kRjUS0VjtBHNxFXzZOsA9wd5Y0/6bA4SEQ2QsRMY5Z+WRsc6FzYRBEPsfR4O95Qw9FaKTwdrFeyxf/ssOBPK9nhcA01qtswG4CqcJ6EogTkSSgMPAb4CFwNc6XK0xXZyqsrekmuy9pXy+M4+KvetIq/mCCSF7uDt0L4PZD5FeG9QAhZ6fZuHREOEJ6shYJ6zjB3iCO/bospYwjz32Ne+AD4vEGG/++jr/IfAHEVmE04xTCDQB3wPeVNWCE3XdEpFbgVsBBg8e7KeSjOl8qsqekmpyduRTvP1TtPBzhjXsYKrsYUHIfkJQCIeG2IGEDZoC/SdC/EBPcLcK7eZlAbxF3/R8voR+ITDI63kax56XoKpFOGf6iEgscLWqlonIOcD5IvI9IBaIEJEqVV3aavtlwDKAjIyMYzsEG9OFqCp7il3s2JhF5e4cepVsZGTTLuZLIaHi/NM9Et0XHTAJGfptGHA2DJhIeExyYAs3xsOX0M8BRojIMJywvxa43nsFEUkGSlXVDdyD05MHVb3Ba51FQEbrwDemK9OGWvK257B/ayZNhZ+RXLGVMzSfM8XpX14ZmkhF8ngqhn6T3sOnIgMmER3X5m1DxnQJ7Ya+qjaKyGLgbZwL8MtVdYuIPADkquoa4ELgFyKiOM07t3dizcZ0jsZ63Ae2cPCLbMr3fEov10b61+1lCI0MAcqIozhmNF/0v4iUkdNJOWs6cfEDieum49SY4CSqXas1JSMjQ3NzcwNdhunpmhrBtR134WeU7/6UxoLP6F2xg3Ccu1vLNZodoSOoSBxH1NAMho6bwYAhI5AQm2HUdE0isk5V2+0haf2yTM9XfwRKd8P+zWjRZ9TsW0eEazNh7jpCgDDtxXb3MPZGXkpTv4kkjZjG+HETmZLkex90Y7oLC33TM7jdUFEAJTtpcu2k7sAXNLl2En54N1FHilpWqyGSTe5hbHLPoih6FJGDMxg+agLTzkzmnMToAB6AMaeHhb7pNmobmigpcVFZuJX6Azvh0E4iy/cQX/UlSfUFRKpzV2so4NYo9mp/9uhQ9rjPZY/2pzJhBP3PmMC04SlcOiyJAb17BfaAjAkAC30TUKpKRW0jrspaDlbU4aqqw1VWTf2hvYQd3kV0xV561+yjX0M+Q7SINDk6bECTCvmksjdkAJ9GjqciZii1CWegSSOI6TOAvvFRDI2PYkpcJCmxkS3DCRsTzCz0TadoaHJTWl2Pq9IT5BV1HKys5WBlHa7KOg5W1FJfeZC4yi9J00JnWjopZpwUMVgOEi5HJ7uoDEmgNGYIJXEzKe59JiEpI4jqN4qEgSMZFB/L0BDrPWOMryz0jc/qGps4VFVPSVWd81NZj6v5cVU9JZV1La81j/EeShNnShFnSDFnSBGTw/czImQ/g7WIWK1q+RfYFBJOXdxQ3H0m0pQ6krDUkUjySEgaTlx0H+ICeNzG9CQW+kGutqHpmNB2VdV5hbdXqFfWtUy311pMRCjJcZEkx0ZyRkoM5w6J4ezG9Yyt+JjBrg+JqPcayTFuACQPh6TzIGkEJI+ApOGE9h5MtA0/YEyns9Dvoeob3WwuKsdVefSsvOUM3SvkK+vaDvK4qDBSYp0gH9UvjuThySR7nifHRpDsaSdPjo10hv2tq4Sd78C2N+CLd50RIaMSYNQlMHw2eM7aiYw9zZ+EMcabhX4PU3D4CC9+msdLOfmUVNUf81rv6PCW0B47IJ7kWGeavOTYiKOBHhdJUkyEbxN0HCmFLa86Qb/7A2iqc8ZhHz8fRl8OQ8+HMJsUw5iuxEK/B2hyKx/vcPFC1j4++OIgAswalcpVZw9kUGI0yXERJMX4qfdKRTFs/7sT9F/+C7TJGZt9yk0weh4MmmqjRBrThVnod2MlVXWszs3nL9l5FByuISUuksX/MZxrpw5moD/7oJfudYJ+6xoo+NRZljQCzrvLOaPvPxFs/BljugUL/W5GVcn58jAvZO3jrc3FNDQp55yRxD2XjOaisamEh/rhbF4VXNuds/lta2D/Jmd5vwkw6z7njD7lrI6/jzHmtLPQ7yYqaxt49fNCXsjax44DVcRFhbFw+hBumDaE4X39cHFUFYo+d0J+2xtwaBcgMGgaXPQwjP46JA7t+PsYYwLKQr+L21JUzgtZeby+vpAj9U1MSEvgV1dP4PL0AV+dLPtkuZsgL8tzRv+GM3aNhMKw82H692DUZc4k2MaYHsNCvwuqbWji/zYW80L2Pj7PKyMqPITLJwxg4fQhpA/q3bGdN9bDlx877fPb/w+OlEBopNOtcta9MHIuRPfxz4EYY7ocC/0u5MuSalZm7+PldQWUHWngjJQY/vvrY7jm7DQSosNPfcf1R2D3+54+9P+AunJnLtYRF8GYeTB8jvWfNyZI+BT6IjIX+C3OAIZ/VtVHWr0+BGeKxBSgFFjomQx9IvAkEI8zUfrDqvqSH+vv9hqb3Ly37SArs/fxz50lhIUIF41NZeG0IZxzZhInmlD+hOqPeLpWroGd70FjDfRKdHrbjL4czrgQwqP8eSjGmG6g3dAXkVDgCWAOUADkiMgaVd3qtdqjwApVfU5EZgG/AL4FHAH+U1V3isgAYJ2IvK2qZX4/km7mQEUtL36ax6pP89lfUUv/hCi+P2ck104ZRN/4DoRxQy2sexb+9RhUHYDYfjBpoRP0Q2ZAqP1xZ0ww8yUBpgK7VHUPgIisAq4AvEN/DPB9z+O1wGsAqrqjeQVVLRKRgzh/DQRl6Ksqn+w+xPOZ+3h32wGa3MoFI1N44IqxzBrVl7COdLdsrIPPVsA/H4PKIhhyHlz1lHNXrE3xZ4zx8CX0BwL5Xs8LgGmt1tkAXIXTBHQlECciSap6qHkFEZkKRAC7W7+BiNwK3AowePDgk6m/Wyg7Us8r6wr4S3Yee0qqSYwO5+bzhnH9tMEM6eiUfE0NsH4lfPwolOfDoOlw1Z9g2AX+Kd4Y06P462/9HwJ/EJFFwMdAIU4bPgAi0h94HvgvVXW33lhVlwHLwJkY3U81BZSqsqGgnBey9vHGhiLqGt1MHpLI/8weziXj+vs2ts2JNDXCxlXw0a+gbB8MzIDLfwtnzrK7Y40xx+VL6BcCg7yep3mWtVDVIpwzfUQkFri6ud1eROKB/wPuVdUsfxTdlTU2ufnrZwU8n7WPzYUVxESEcs3kNBZOH8Lo/vEdfwN3E2x6BT56BEr3OEMgXPoojJhjYW+MaZcvoZ8DjBCRYThhfy1wvfcKIpIMlHrO4u/B6cmDiEQAr+Jc5H3Fn4V3RVV1jSz5y2es/cLFqH5xPPiNcVw5aSCxkX74g8rthi1/g49+CSU7IHU8XPsXOOtSC3tjjM/aTSNVbRSRxcDbOF02l6vqFhF5AMhV1TXAhcAvRERxmndu92z+TeACIMnT9AOwSFXX+/cwAq+4vIZvP5vLjgOVPPSNcdwwbfCpd7f05nY73S4/fARc2yBlNHxzBYy63C7QGmNOmqh2rSb0jIwMzc3NDXQZJ2VLUTnffjaH6romnrjhbGaOTOn4TlXhizdh7S/gwCZnVMsLl8LYqyzsjTFfISLrVDWjvfWs03YHfbD9AIv/8jkJvcJ5+bZzOt5ur+rMQLX251C8HvqcAVcug/HX2Dj1xpgOs9DvgBWZX3L/mi2MGRDP0/81hdSO3FSl6sw+tfbnUJgLvYfAFf8LExbYDVXGGL+xNDkFTW7l529u4+l/7eVro/vyu+smER3RgY9yz0dO2OdnQXya0/Vy4g0Q2oHxdowxpg0W+ifpSH0jd65az7tbD3DjjKHcd9kYQkNO8YLtvk+csP/ynxA3AC77DUz6FoRF+rdoY4zxsNA/CQcrarnpuVy2FJVz/+VjWDRj2KntKP9TWPsw7PkQYlNh7i9h8iIbAM0Y0+ks9H20fX8F334mh7KaBp76zwxmj049+Z0UrIMPfw673oPoZGdGqoxvQ0S0/ws2xpg2WOj74KMdLm5f+RkxkaGs/s45jBuYcHI7KN7gdL3c8ZYzvPHX7oept0JEB8fdMcaYk2Sh346/ZOfx369vZkTfWJ65cQr9E3r5vvGBLU6b/fa/Q1SCM6n41O9AlB+GYzDGmFNgoX8cbrfyy7e386eP9nDhWSn84fqzfR9O4eB2Z2ycLa9CZDxceA9M/64T/MYYE0AW+m2obWji7pfW89bm/SycPpj7Lx/r+1j3+zfDny6A8F5w/g/hnNttzlljTJdhod9KSVUdNz+Xy4aCMu67bDQ3nTfs5MbQ2fkOaBPcng0JaZ1XqDHGnAILfS87D1Ry47M5lFTV8eQNk5k7rt/J7yQvC5LPssA3xnRJFvoen+wq4TsvrCMyLJSXbj2H9EG9T34nbrdzV+2Yb/i/QGOM8QMLfWB1bj4//tsmzkiJYfmiKaQlnmK/+ZIvoLYcBk/3b4HGGOMnQR36qspv3tnBH9bu4rzhyfzvwrOJj+rAeDd5mc5vC31jTBflU5cUEZkrIl+IyC4RWdrG60NE5H0R2SgiH4pImtdr/yUiOz0//+XP4juitqGJO1et5w9rd7EgYxDP3DilY4EPTnt+bCoknuLwDMYY08naPdMXkVDgCWAOUADkiMgaVd3qtdqjOFMiPicis4BfAN8SkT7AT4EMQIF1nm0P+/tATkZpdT23rsgld99h/t/cUdw28wz/zHKVlwWDptn0hcaYLsuXM/2pwC5V3aOq9cAq4IpW64wBPvA8Xuv1+sXAu6pa6gn6d4G5HS/71O1xVXHl//6bjYXlPHH92Xz3wjP9E/gVxVC2Dwaf0/F9GWNMJ/El9AcC+V7PCzzLvG0ArvI8vhKIE5EkH7c9bbL3HOKqJz+hsraRF2+ZzmUT+vtv5/lZzu/B0/y3T2OM8TN/Tbb6Q2CmiHwOzAQKgSZfNxaRW0UkV0RyXS6Xn0o61qufF7Dw6Wz6xETw2vdmMHlIon/fIC8LwqOh3wT/7tcYY/zIl9AvBAZ5PU/zLGuhqkWqepWqTgLu9Swr82Vbz7rLVDVDVTNSUvwwqfix++bx93Zw90sbmDwkkVe/O4PBSZ0wlHFeJqRl2GxXxpguzZfQzwFGiMgwEYkArgXWeK8gIski0ryve4DlnsdvAxeJSKKIJAIXeZadFnWNTfxg9QYef28nV5+dxopvTyMhuhNCua4S9m+CQdZV0xjTtbXbe0dVG0VkMU5YhwLLVXWLiDwA5KrqGuBC4BciosDHwO2ebUtF5EGcLw6AB1S1tBOO4yvKjtTznefXkb23lO/PGcmSWcP9c8G2LQW5oG7rn2+M6fJ8ujlLVd8E3my17Cdej18BXjnOtss5euZ/Wuw7VM2Nz+ZQUFrD4wsm8o1JnXztOD8bJATSpnTu+xhjTAf1uDty1+0r5ZYV63Cr8sLN05g67DQMa5yXCaljbXIUY0yX56/eO13CGxuKuO6pbOKjwvjbd889PYHf1Aj5OdY/3xjTLfSYM/1dB6u4Y9XnZAxJ5E/fyqBPTMTpeeMDm6Gh2rkT1xhjurgeE/rD+8byp4WTuWBkClHhoafvjfOab8qyM31jTNfXY0If4KKxpzDpSUflZULCYEgI2I3Gxhjjsx7Vpn/aqTo9d2zoBWNMN2Gh3xFl+6Cy2PrnG2O6DQv9jsjLdn7bnbjGmG7CQr8j8jIhMgH6jg50JcYY4xML/Y7Iy4JBUyHkNPYWMsaYDrDQP1U1h8G1zS7iGmO6FQv9U5X/qfPb+ucbY7oRC/1TlZcJIWEw4OxAV2KMMT6z0D9VednQfyJEdMKELMYY00ks9E9FYx0UrrP++caYbsen0BeRuSLyhYjsEpGlbbw+WETWisjnIrJRRC71LA8XkedEZJOIbBORe/x9AAFRvAGa6iz0jTHdTruhLyKhwBPAJcAY4DoRGdNqtfuA1Z45cq8F/tezfD4QqarjgcnAd0RkqH9KD6C8TOe33ZRljOlmfDnTnwrsUtU9qloPrAKuaLWOAs0ziCQARV7LY0QkDOgF1AMVHa460PKyIGk4xPp3EndjjOlsvoT+QCDf63mBZ5m3+4GFIlKAM63iEs/yV4BqoBjIAx49XXPkdhpVz01ZdpZvjOl+/HUh9zrgWVVNAy4FnheREJy/EpqAAcAw4AcickbrjUXkVhHJFZFcl8vlp5I6SclOqCm19nxjTLfkS+gXAoO8nqd5lnm7CVgNoKqZQBSQDFwP/ENVG1T1IPBvIKP1G6jqMlXNUNWMlJQu3mSS3zxpioW+Mab78SX0c4ARIjJMRCJwLtSuabVOHjAbQERG44S+y7N8lmd5DDAd2O6f0gMkLwuik5w2fWOM6WbaDX1VbQQWA28D23B66WwRkQdEZJ5ntR8At4jIBuBFYJGqKk6vn1gR2YLz5fGMqm7sjAM5bfIynaEXRAJdiTHGnDSfpktU1TdxLtB6L/uJ1+OtwIw2tqvC6bbZM1QdhNI9MPnGQFdijDGnxO7IPRk2Cboxppuz0D8ZeVkQFgX90wNdiTHGnBIL/ZORnwUDJ0NYRKArMcaYU2Kh76v6I86YO9ZV0xjTjVno+6pwHbgb7U5cY0y3ZqHvq7wsQGDQlEBXYowxp8xC31d5mdB3DPRKDHQlxhhzyiz0feFugoIcmwTdGNPtWej74uBWqKuw/vnGmG7PQt8XzTdlDbIzfWNM92ah74u8LIgbAL0HB7oSY4zpEAt9X+RlOf3zbZA1Y0w3Z6HfnrJ8qCiwm7KMMT2ChX578rOd3xb6xpgewEK/PXmZEBELfccGuhJjjOkwC/325GVD2hQI9WnqAWOM6dJ8Cn0RmSsiX4jILhFZ2sbrg0VkrYh8LiIbReRSr9cmiEimiGwRkU0iEuXPA+hUteVwYLP1zzfG9Bjtnr6KSCjOtIdzgAIgR0TWeGbLanYfzjSKT4rIGJxZtoaKSBjwAvAtVd0gIklAg9+PorMU5ABqd+IaY3oMX870pwK7VHWPqtYDq4ArWq2jQLzncQJQ5Hl8EbBRVTcAqOohVW3qeNmnSV4WSCgMzAh0JcYY4xe+hP5AIN/reYFnmbf7gYUiUoBzlr/Es3wkoCLytoh8JiI/ausNRORWEckVkVyXy3VSB9Cp8rKg33iIjA10JcYY4xf+upB7HfCsqqYBlwLPi0gITvPRecANnt9Xisjs1hur6jJVzVDVjJSUFD+V1EFNDVCQa+35xpgexZfQLwQGeT1P8yzzdhOwGkBVM4EoIBnnr4KPVbVEVY/g/BVwdkeLPi2KN0JjjfXPN8b0KL6Efg4wQkSGiUgEcC2wptU6ecBsABEZjRP6LuBtYLyIRHsu6s4EttId5GU6vy30jTE9SLu9d1S1UUQW4wR4KLBcVbeIyANArqquAX4APCUid+Nc1F2kqgocFpHHcL44FHhTVf+vsw7Gr/KzIHEoxPULdCXGGOM3Pt1xpKpv4jTNeC/7idfjrcCM42z7Ak63ze5D1bmIO/xrga7EGGP8yu7IbUvpHqh22fj5xpgex0K/Lc2TpljPHWNMD2Oh35a8TIjqDckjA12JMcb4lYV+W/KznV47IfbxGGN6Fku11qoPQckO66ppjOmRLPRba540ZZCFvjGm57HQby0vE0IjYMCkQFdijDF+Z6HfWl6WE/jh3WfYf2OM8ZWFvreGGij63NrzjTE9loW+t6LPwd1g/fONMT2Whb635kHW7E5cY0wPZaHvLS8bks+C6D6BrsQYYzqFhX4zt9sZWdPa840xPZiFfrOSL6C23ELfGNOjWeg3s0lTjDFBwKfQF5G5IvKFiOwSkaVtvD5YRNaKyOcislFELm3j9SoR+aG/Cve7vCyI6QuJwwJdiTHGdJp2Q19EQoEngEuAMcB1IjKm1Wr3AatVdRLOdIr/2+r1x4C3Ol5uJ8rztOeLBLoSY4zpNL6c6U8FdqnqHlWtB1YBV7RaR4F4z+MEoKj5BRH5BrAX2NLxcjtJRRGU7bP++caYHs+X0B8I5Hs9L/As83Y/sFBECnCmVVwCICKxwP8DftbhSjtTy6Qp1j/fGNOz+etC7nXAs6qaBlwKPC8iIThfBv+jqlUn2lhEbhWRXBHJdblcfirpJORnQ3g09Jtw+t/bGGNOI18mRi8EBnk9T/Ms83YTMBdAVTNFJApIBqYB14jIr4DegFtEalX1D94bq+oyYBlARkaGnsqBdEheJgycDKHhp/2tjTHmdPLlTD8HGCEiw0QkAudC7ZpW6+QBswFEZDQQBbhU9XxVHaqqQ4HHgZ+3DvyAq6uE/ZusPd8YExTaDX1VbQQWA28D23B66WwRkQdEZJ5ntR8At4jIBuBFYJGqnv4z9lNRkAvqtv75xpig4EvzDqr6Js4FWu9lP/F6vBWY0c4+7j+F+jpfXhZICKRNCXQlxhjT6eyO3PwsSB0LUfHtr2uMMd1ccId+UyPk59h8uMaYoBHcoX9gMzRUW3u+MSZoBHfot9yUZT13jDHBIchDPxMSBkFC6xuMjTGmZwre0Fd17sS1ph1jTBAJ3tAv2weVxRb6xpigEryh39yebz13jDFBJLhDPzIB+o4OdCXGGHPaBHfoD5oCIaGBrsQYY06b4Az9msPg2mbt+caYoBOcoZ//qfPb+ucbY4JMcIZ+XiaEhMEezjBlAAAV30lEQVSAswNdiTHGnFZBGvrZ0H8iREQHuhJjjDmtgi/0G+ugcJ215xtjgpJPoS8ic0XkCxHZJSJL23h9sIisFZHPRWSjiFzqWT5HRNaJyCbP71n+PoCTVrQemuos9I0xQandSVREJBR4ApgDFAA5IrLGM3FKs/twZtR6UkTG4Ey4MhQoAS5X1SIRGYcz+1ZgB7rJt5uyjDHBy5cz/anALlXdo6r1wCrgilbrKNA8C0kCUASgqp+rapFn+Ragl4hEdrzsDsjLgj5nQmxKQMswxphA8CX0BwL5Xs8L+OrZ+v3AQhEpwDnLX9LGfq4GPlPVulOo0z9UndC3rprGmCDlrwu51wHPqmoacCnwvIi07FtExgK/BL7T1sYicquI5IpIrsvl8lNJbSjZCTWl1p5vjAlavoR+ITDI63maZ5m3m4DVAKqaCUQByQAikga8Cvynqu5u6w1UdZmqZqhqRkpKJza75GU6vy30jTFBypfQzwFGiMgwEYkArgXWtFonD5gNICKjcULfJSK9gf8Dlqrqv/1X9inKz4boJEgaHuhKjDEmINoNfVVtBBbj9LzZhtNLZ4uIPCAi8zyr/QC4RUQ2AC8Ci1RVPdsNB34iIus9P3075Uh8kZfp9NoRCVgJxhgTSO122QRQ1TdxLtB6L/uJ1+OtwIw2tnsIeKiDNfpH1UEo3QOTbwx0JcYYEzDBc0euTYJujDFBFvphUdA/PdCVGGNMwARP6OdnwcDJEBYR6EqMMSZggiP066uheAMMmhboSowxJqB8upDb7RWuA3ejtecbv2hoaKCgoIDa2tpAl2KCUFRUFGlpaYSHh5/S9sER+nnZgDhz4hrTQQUFBcTFxTF06FDEuv+a00hVOXToEAUFBQwbNuyU9hEczTt5mdB3NPRKDHQlpgeora0lKSnJAt+cdiJCUlJSh/7K7Pmh726CghwbesH4lQW+CZSO/tvr+aF/cCvUVVh7vjHGEAyh33xTlvXcMUEqNjbWb/t67bXX2Lp1a/sr+sG55557Stvdf//9PProo36upucIjtCPGwC9Bwe6EmO6vdMR+o2NjQB88sknnfo+na35OLqant97Jy8LBk+zQdZMp/jZG1vYWlTh132OGRDPTy8fe9zXly5dyqBBg7j99tsB58w2LCyMtWvXcvjwYRoaGnjooYe44orWE9y17Ze//CUvvPACISEhXHLJJTzyyCM89dRTLFu2jPr6eoYPH87zzz/P+vXrWbNmDR999BEPPfQQf/3rXwG4/fbbcblcREdH89RTTzFq1Ch2797NDTfcQHV1NVdccQWPP/44VVVVqCo/+tGPeOuttxAR7rvvPhYsWMCHH37If//3f5OYmMj27dvZsWMHsbGxVFVVnVSN0dHR7R7v8bY7cOAAt912G3v27AHgySef5Nxzz2XFihU8+uijiAgTJkzg+eefZ9GiRXz961/nmmuuAWipta3j+MY3vkF+fj61tbXceeed3HrrrQD84x//4Mc//jFNTU0kJyfz7rvvctZZZ/HJJ5+QkpKC2+1m5MiRZGZm4s8h53t26JflQ0UBDL4j0JUY4zcLFizgrrvuagn91atX8/bbb3PHHXcQHx9PSUkJ06dPZ968ee1e9Hvrrbd4/fXXyc7OJjo6mtLSUgCuuuoqbrnlFgDuu+8+nn76aZYsWcK8efOOCbvZs2fzxz/+kREjRpCdnc33vvc9PvjgA+68807uvPNOrrvuOv74xz+2vN/f/vY31q9fz4YNGygpKWHKlClccMEFAHz22Wds3rz5K10RT7bG9hxvuzvuuIOZM2fy6quv0tTURFVVFVu2bOGhhx7ik08+ITk5ueW9T6T1cSxfvpw+ffpQU1PDlClTuPrqq3G73dxyyy18/PHHDBs2jNLSUkJCQli4cCErV67krrvu4r333iM9Pd2vgQ89PfTzs53f1nPHdJITnZF3lkmTJnHw4EGKiopwuVwkJibSr18/7r77bj7++GNCQkIoLCzkwIED9OvX74T7eu+997jxxhtbzpD79OkDwObNm7nvvvsoKyujqqqKiy+++CvbVlVV8cknnzB//vyWZXV1zmyomZmZvPbaawBcf/31/PCHPwTgX//6F9dddx2hoaGkpqYyc+ZMcnJyiI+PZ+rUqW32Pe9IjW053nYffPABK1asACA0NJSEhARWrFjB/PnzSU5OPua9T6T1cfzud7/j1VdfBSA/P5+dO3ficrm44IILWtZr3u+3v/1trrjiCu666y6WL1/OjTf6f1Tgnh36eZkQEQt9T///mMZ0pvnz5/PKK6+wf/9+FixYwMqVK3G5XKxbt47w8HCGDh3aob7cixYt4rXXXiM9PZ1nn32WDz/88CvruN1uevfuzfr16ztwJEfFxMT4vUZ/buctLCwMt9sNOJ9DfX19y2vex/Hhhx/y3nvvkZmZSXR0NBdeeOEJ/7sMGjSI1NRUPvjgAz799FNWrlx50rW1p2dfyM3LhrQpENqzv9tM8FmwYAGrVq3ilVdeYf78+ZSXl9O3b1/Cw8NZu3Yt+/bt82k/c+bM4ZlnnuHIkSMALc0XlZWV9O/fn4aGhmOCJy4ujsrKSgDi4+MZNmwYL7/8MuDcLbphwwYApk+f3tLmv2rVqpbtzz//fF566SWamppwuVx8/PHHTJ061a81tud4282ePZsnn3wSgKamJsrLy5k1axYvv/wyhw4dOua9hw4dyrp16wBYs2YNDQ0Nbb5XeXk5iYmJREdHs337drKyslo+n48//pi9e/ces1+Am2++mYULFzJ//nxCQ0N9Pi5f+RT6IjJXRL4QkV0isrSN1weLyFoR+VxENorIpV6v3ePZ7gsR8e3vL3+oLYcDm61px/RIY8eOpbKykoEDB9K/f39uuOEGcnNzGT9+PCtWrGDUqFE+7Wfu3LnMmzePjIwMJk6c2NLV8cEHH2TatGnMmDHjmH1de+21/PrXv2bSpEns3r2blStX8vTTT5Oens7YsWN5/fXXAXj88cd57LHHmDBhArt27SIhIQGAK6+8kgkTJpCens6sWbP41a9+1W4T1MnW2J7jbffb3/6WtWvXMn78eCZPnszWrVsZO3Ys9957LzNnziQ9PZ3vf//7ANxyyy189NFHpKenk5mZedy/UubOnUtjYyOjR49m6dKlTJ/u5FFKSgrLli3jqquuIj09nQULFrRsM2/ePKqqqjqlaQdwvp1P9AOEAruBM4AIYAMwptU6y4Dveh6PAb70erwBiASGefYTeqL3mzx5svrFjndVfxqvunutf/ZnjMfWrVsDXUKXV11drW63W1VVX3zxRZ03b16AK+o+cnJy9LzzzjvhOm39GwRytZ08V1Wf2vSnArtUdQ+AiKwCrgC8O+sqEO95nAAUeR5fAaxS1Tpgr4js8uwv8+S+mk5BfhZIKAzM6PS3MsYca926dSxevBhVpXfv3ixfvjzQJXULjzzyCE8++WSntOU38yX0BwL5Xs8LgNa3t94PvCMiS4AY4Gte22a12nZg6zcQkVuBWwEGD/bTTVR5WdBvPET6725EY7qrTZs28a1vfeuYZZGRkWRnZ3fK+51//vkt7fuBcvvtt/Pvf//7mGV33nln5zWb+MHSpUtZuvQrLeh+5a8rnNcBz6rqb0TkHOB5ERnn68aqugyniYiMjAztcDVNDVCQC5MXdXhXxvQE48eP91svm+7iiSeeCHQJXZIvoV8IDPJ6nuZZ5u0mYC6AqmaKSBSQ7OO2/le8ERprnDtxjTHGtPCl904OMEJEholIBHAtsKbVOnnAbAARGQ1EAS7PeteKSKSIDANGAJ/6q/jjyvNcMhhkPXeMMcZbu2f6qtooIouBt3F68ixX1S0i8gDO1eI1wA+Ap0TkbpyLuos8V5O3iMhqnIu+jcDtqtrUWQfTIj8LEodCfP9OfytjjOlOfGrTV9U3gTdbLfuJ1+OtwIzjbPsw8HAHajw5qs5F3DNnn7a3NMaY7qLn3ZFbugeqXXZTljE4wwCcriGKL730UsrKyk56u2effZbFixd3QkWmLT1vfILmSVNspixj+PDDD4mNjT3lCUl80XzTz5tvvtn+yl1Y83GEhPS8c2FvPTD0MyGqNySPDHQlJhi8tRT2b/LvPvuNh0seOeEqrcd4/+Y3v8lDDz1EfX09SUlJrFy5kpqaGv74xz8SGhrKCy+8wO9//3tGjRrFbbfdRl5eHuAMlzBjxgxcLhfXX389RUVFnHPOObz77rusW7eO5ORkHnvssZabq26++WbuuusuvvzySy6++GKmTZvGunXrePPNN5k5cya5ubkkJye3OQb9G2+88ZUaU1NT2/04jrddVVUVS5YsITc3FxHhpz/9KVdfffVXxql///33uf/++4mNjW0Z7XPcuHH8/e9/B/jKcTzyyCPk5ORQU1PDNddcw89+9jMAcnJyuPPOO6muriYyMpL333+fyy67jN/97ndMnDgRgPPOO48nnniC9PT0U/tvfxr0vNDPz3aadnr4t7UJXm2N8S4iZGVlISL8+c9/5le/+hW/+c1vuO22244Ju+uvv567776b8847j7y8PC6++GK2bdvGz372M2bNmsU999zDP/7xD55++mnAubP2mWeeITs7G1Vl2rRpzJw5k8TERHbu3Mlzzz3XMp7MieoDJxDbqrE9x9vuwQcfJCEhgU2bnC/dw4cP43K5vjJOfXtaH8fDDz9Mnz59aGpqYvbs2WzcuJFRo0axYMECXnrpJaZMmUJFRQW9evXipptu4tlnn+Xxxx9nx44d1NbWdunAh54W+tUlULID0q8LdCUmWLRzRt4ZPvjgg6+M8b5p0yYWLFhAcXEx9fX1bY5LD87Y9N7THVZUVFBVVcW//vWvljHf586dS2JiIuCMf3/llVe2DCh21VVX8c9//pN58+YxZMiQrwT+8eoDKCgo8KnG1o633XvvvXfMCJ6JiYm88cYbbY5TfyKtj2P16tUsW7aMxsZGiouL2bp1KyJC//79mTJlCuCMMArOENcPPvggv/71r1m+fDmLFi3y6ZgCqWedDrdMmmLt+Sa4LFmyhMWLF7Np0yb+9Kc/HXfMdrfbTVZWFuvXr2f9+vUUFhae8sTpJzv+va81+ms7b97j3wPH7MP7OPbu3cujjz7K+++/z8aNG7nssstO+H7R0dHMmTOH119/ndWrV3PDDTecdG2nW88K/bwsCI2AAZMCXYkxnaatMd7Ly8sZONAZ1uq5555rWdd7/HuAiy66iN///vctz5uHZpgxYwarV68G4J133uHw4cOAM4bOa6+9xpEjR6iurubVV1/l/PPPP+n6gOPW2J7jbTdnzpxjhlo4fPjwccepHzp0KJ999hngTGfY/HprFRUVxMTEkJCQwIEDB3jrrbcAOOussyguLiYnJwdwxuRvnvj85ptv5o477mDKlCktfyF1ZT0v9AdMgvCoQFdiTKdpa4z3+++/n/nz5zN58uSWZhWAyy+/nFdffZWJEyfyz3/+k9/97nfk5uYyYcIExowZ0zJ/7U9/+lPeeecdxo0bx8svv0y/fv2Ii4vj7LPPZtGiRUydOpVp06Zx8803M2nSiU+qjjcG/fFqbM/xtrvvvvs4fPgw48aNIz09nbVr1x53nPqrr76a0tJSxo4dyx/+8AdGjmy7o0d6ejqTJk1i1KhRXH/99cyY4dx+FBERwUsvvcSSJUtIT09nzpw5LX8BTJ48mfj4+C49kJs3cW6c7ToyMjI0Nzf35DdsqIFfDIJzvgdzHvB/YcZ4bNu2jdGjRwe6DL+qq6sjNDSUsLAwMjMz+e53vxt0A7SdqqKiIi688EK2b99+2rp7tvVvUETWqWq7Y8n3nAu5dZUw9htw5qxAV2JMt5OXl8c3v/lN3G43ERERPPXUU4EuqVtYsWIF9957L4899li36d/fc870jTlNeuKZflfw8MMPt8y322z+/Pnce++9Aaqo67IzfWNMt3fvvfdawJ8G3ePvEWO6mK72F7IJHh39t2ehb8xJioqK4tChQxb85rRTVQ4dOkRU1Kn3ULTmHWNOUlpaGgUFBbhcrkCXYoJQVFQUaWlpp7y9T6EvInOB3+JMovJnVX2k1ev/A/yH52k00FdVe3te+xVwGc5fFe8Cd6qdIpluLDw83OchBIzpatoNfREJBZ4A5gAFQI6IrPFMnAKAqt7ttf4SYJLn8bk4k6tM8Lz8L2Am8KGf6jfGGHMSfGnTnwrsUtU9qloPrAKuOMH61wEveh4rzny5EUAkEA4cOPVyjTHGdIQvoT8QyPd6XuBZ9hUiMgQYBnwAoKqZwFqg2PPztqpu60jBxhhjTp2/L+ReC7zSPPm5iAwHRgPNVx3eFZHzVfWf3huJyK3ArZ6nVSLyRQdqSAZKOrB9T2KfxbHs8ziWfR5H9YTPYogvK/kS+oXAIK/naZ5lbbkWuN3r+ZVAlqpWAYjIW8A5wDGhr6rLgGW+FNweEcn15a60YGCfxbHs8ziWfR5HBdNn4UvzTg4wQkSGiUgETrCvab2SiIwCEoFMr8V5wEwRCRORcJyLuNa8Y4wxAdJu6KtqI7AYeBsnsFer6hYReUBE5nmtei2wqlV3zFeA3cAmYAOwQVXf8Fv1xhhjTopPbfqq+ibwZqtlP2n1/P42tmsCvtOB+k6FX5qJegj7LI5ln8ex7PM4Kmg+iy43yqYxxpjOY2PvGGNMEOkxoS8ic0XkCxHZJSJLA11PIInIIBFZKyJbRWSLiNwZ6JoCTURCReRzEfl7oGsJNBHpLSKviMh2EdkmIucEuqZAEpG7Pf+fbBaRF0WkR8+32iNC32uoiEuAMcB1IjImsFUFVCPwA1UdA0wHbg/yzwPgTqznWLPfAv9Q1VFAOkH8uYjIQOAOIENVx+GML3ZtYKvqXD0i9Dn5oSJ6NFUtVtXPPI8rcf6nbvMu6mAgImk4g/79OdC1BJqIJAAXAE8DqGq9qpYFtqqACwN6iUgYzoCRRQGup1P1lND3eaiIYCMiQ3EGwMsObCUB9TjwI8Ad6EK6gGGAC3jG09z1ZxGJCXRRgaKqhcCjOPcUFQPlqvpOYKvqXD0l9E0bRCQW+Ctwl6pWBLqeQBCRrwMHVXVdoGvpIsKAs4EnVXUSUA0E7TUwEUnEaRUYBgwAYkRkYWCr6lw9JfRPZqiIoOC5A/qvwEpV/Vug6wmgGcA8EfkSp9lvloi8ENiSAqoAKFDV5r/8XsH5EghWXwP2qqpLVRuAvwHnBrimTtVTQt+noSKChYgITpvtNlV9LND1BJKq3qOqaao6FOffxQeq2qPP5E5EVfcD+SJylmfRbGDrCTbp6fKA6SIS7fn/ZjY9/MJ2j5guUVUbRaR5qIhQYLmqbglwWYE0A/gWsElE1nuW/dhzZ7UxS4CVnhOkPcCNAa4nYFQ1W0ReAT7D6fX2OT387ly7I9cYY4JIT2neMcYY4wMLfWOMCSIW+sYYE0Qs9I0xJohY6BtjTBCx0DfGmCBioW+MMUHEQt8YY4LI/wdZbV4ez6w6yQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Creating network for 784-dimensional input and 10-dimensional output, with 2 hidden layers\n",
      "_________________________________________________________________\n",
      "Layer (type)                 Output Shape              Param #   \n",
      "=================================================================\n",
      "dense_46 (Dense)             (None, 100)               78500     \n",
      "_________________________________________________________________\n",
      "dense_47 (Dense)             (None, 200)               20200     \n",
      "_________________________________________________________________\n",
      "Crashes (Lambda)             (None, 200)               0         \n",
      "_________________________________________________________________\n",
      "dense_48 (Dense)             (None, 10)                2010      \n",
      "=================================================================\n",
      "Total params: 100,710\n",
      "Trainable params: 100,710\n",
      "Non-trainable params: 0\n",
      "_________________________________________________________________\n",
      "_________________________________________________________________\n",
      "Layer (type)                 Output Shape              Param #   \n",
      "=================================================================\n",
      "dense_49 (Dense)             (None, 100)               78500     \n",
      "_________________________________________________________________\n",
      "dense_50 (Dense)             (None, 200)               20200     \n",
      "_________________________________________________________________\n",
      "dense_51 (Dense)             (None, 10)                2010      \n",
      "=================================================================\n",
      "Total params: 100,710\n",
      "Trainable params: 100,710\n",
      "Non-trainable params: 0\n",
      "_________________________________________________________________\n"
     ]
    }
   ],
   "source": [
    "exp = MNISTExperiment(N = N, p_inference = p_inference, p_train = get_p_arr(0), KLips = KLips, epochs = 10,\n",
    "                                      activation = activation, reg_type = 'balanced', reg_coeff = 1e-5, do_print = True,\n",
    "                                      name = 'experiment_weights', seed = None, batch_size = 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 265,
   "metadata": {},
   "outputs": [],
   "source": [
    "truth = exp.y_train"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {},
   "outputs": [],
   "source": [
    "def crashing_loss():\n",
    "    prediction = exp.predict_crashing(exp.x_train, repetitions = 1)\n",
    "    prediction = prediction.reshape(60000, 10)\n",
    "    loss = np.mean(np.square(prediction - truth))\n",
    "    return loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "metadata": {},
   "outputs": [],
   "source": [
    "def correct_loss():\n",
    "    prediction = exp.predict_correct(exp.x_train)\n",
    "    #prediction = prediction.reshape(60000, 10)\n",
    "    loss = np.mean(np.square(prediction - truth))\n",
    "    return loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {},
   "outputs": [],
   "source": [
    "correct = correct_loss()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 20/20 [00:15<00:00,  1.28it/s]\n"
     ]
    }
   ],
   "source": [
    "crashing = [crashing_loss() for _ in tqdm(range(20))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.7836008603728715"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "correct * 1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([2., 2., 1., 2., 2., 5., 3., 0., 2., 1.]),\n",
       " array([7.79040083, 7.79100509, 7.79160934, 7.7922136 , 7.79281786,\n",
       "        7.79342211, 7.79402637, 7.79463062, 7.79523488, 7.79583914,\n",
       "        7.79644339]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAC8hJREFUeJzt3H+M5Hddx/HX2x6oVASxiyL0OIJAICRCXIkJakKJpFJ+RI0JBE1QkwsmIBiEVEkMakyKxh//GM0pRCMIQaEJUn4UAwUxUmxriy0tBskRaJCCaICYYNq+/WOHuL3u3c6V+e72vTweyeR2Zr478/7czD7ve9+Z2eruADDHtxz2AACcH+EGGEa4AYYRboBhhBtgGOEGGEa4AYYRboBhhBtgmGNL3OhFF13UJ06cWOKmAY6k66+//ovdvbXOtouE+8SJE7nuuuuWuGmAI6mqPr3utg6VAAwj3ADDCDfAMMINMIxwAwyz1rtKqup0kq8kuSvJnd29veRQAJzd+bwd8Bnd/cXFJgFgLQ6VAAyzbrg7ydVVdX1VnVxyIADObd1DJT/S3bdX1cOTvK+qbuvuD+3eYBX0k0ly/PjxDY8Jm3Hi8qsO5X5PX3HZodwvR9Nae9zdffvqzzuSXJnkaXtsc6q7t7t7e2trrY/bA3Af7Bvuqrqwqh789a+TPCvJzUsPBsDe1jlU8j1Jrqyqr2//1939nkWnAuCs9g13d38qyQ8cwCwArMHbAQGGEW6AYYQbYBjhBhhGuAGGEW6AYYQbYBjhBhhGuAGGEW6AYYQbYBjhBhhGuAGGEW6AYYQbYBjhBhhGuAGGEW6AYYQbYBjhBhhGuAGGEW6AYYQbYBjhBhhGuAGGEW6AYYQbYBjhBhhGuAGGEW6AYYQbYBjhBhhm7XBX1QVV9S9V9c4lBwLg3M5nj/vlSW5dahAA1rNWuKvqUUkuS/Lny44DwH7W3eP+oySvTnL3grMAsIZ9w11Vz0lyR3dfv892J6vquqq67gtf+MLGBgTgntbZ4356kudV1ekkb0lySVW98cyNuvtUd2939/bW1taGxwTg6/YNd3f/Wnc/qrtPJHlBkvd3988uPhkAe/I+boBhjp3Pxt19TZJrFpkEgLXY4wYYRrgBhhFugGGEG2AY4QYYRrgBhhFugGGEG2AY4QYYRrgBhhFugGGEG2AY4QYYRrgBhhFugGGEG2AY4QYYRrgBhhFugGGEG2AY4QYYRrgBhhFugGGEG2AY4QYYRrgBhhFugGGEG2AY4QYYRrgBhhFugGGEG2AY4QYYZt9wV9W3VdVHq+qmqrqlqn7zIAYDYG/H1tjma0ku6e6vVtUDkny4qt7d3R9ZeDYA9rBvuLu7k3x1dfYBq1MvORQAZ7fWMe6quqCqbkxyR5L3dfe1y44FwNmsc6gk3X1XkqdU1UOTXFlVT+7um3dvU1Unk5xMkuPHj298UJjsxOVXHdp9n77iskO7b5ZxXu8q6e7/TvKBJJfucd2p7t7u7u2tra1NzQfAGdZ5V8nWak87VfXtSX48yW1LDwbA3tY5VPKIJH9ZVRdkJ/Rv7e53LjsWAGezzrtKPpbkqQcwCwBr8MlJgGGEG2AY4QYYRrgBhhFugGGEG2AY4QYYRrgBhhFugGGEG2AY4QYYRrgBhhFugGGEG2AY4QYYRrgBhhFugGGEG2AY4QYYRrgBhhFugGGEG2AY4QYYRrgBhhFugGGEG2AY4QYYRrgBhhFugGGEG2AY4QYYRrgBhtk33FV1cVV9oKo+XlW3VNXLD2IwAPZ2bI1t7kzyyu6+oaoenOT6qnpfd3984dkA2MO+e9zd/bnuvmH19VeS3JrkkUsPBsDezusYd1WdSPLUJNcuMQwA+1vnUEmSpKq+I8nbkryiu7+8x/Unk5xMkuPHj9/ngU5cftV9/t6pTl9x2aHc7zfj3zUH5zCfX4f1M3VQ1trjrqoHZCfab+rut++1TXef6u7t7t7e2tra5IwA7LLOu0oqyeuT3Nrdf7D8SACcyzp73E9P8nNJLqmqG1enZy88FwBnse8x7u7+cJI6gFkAWINPTgIMI9wAwwg3wDDCDTCMcAMMI9wAwwg3wDDCDTCMcAMMI9wAwwg3wDDCDTCMcAMMI9wAwwg3wDDCDTCMcAMMI9wAwwg3wDDCDTCMcAMMI9wAwwg3wDDCDTCMcAMMI9wAwwg3wDDCDTCMcAMMI9wAwwg3wDDCDTCMcAMMs2+4q+oNVXVHVd18EAMBcG7r7HH/RZJLF54DgDXtG+7u/lCSLx3ALACs4dimbqiqTiY5mSTHjx/f1M1+Uzhx+VWHPQIcKYf1M3X6issO5H429uJkd5/q7u3u3t7a2trUzQJwBu8qARhGuAGGWeftgG9O8k9JnlBVn62qX1x+LADOZt8XJ7v7hQcxCADrcagEYBjhBhhGuAGGEW6AYYQbYBjhBhhGuAGGEW6AYYQbYBjhBhhGuAGGEW6AYYQbYBjhBhhGuAGGEW6AYYQbYBjhBhhGuAGGEW6AYYQbYBjhBhhGuAGGEW6AYYQbYBjhBhhGuAGGEW6AYYQbYBjhBhhGuAGGEW6AYdYKd1VdWlWfqKpPVtXlSw8FwNntG+6quiDJHyf5iSRPSvLCqnrS0oMBsLd19rifluST3f2p7v7fJG9J8vxlxwLgbNYJ9yOTfGbX+c+uLgPgEBzb1A1V1ckkJ1dnv1pVn9jUbR+Qi5J88bCH2ICjso7k6KzlUNdRr9vYTR2VxyNZaC3f4N/1o9fdcJ1w357k4l3nH7W67B66+1SSU+ve8f1NVV3X3duHPcc36qisIzk6a7GO+5/pa1nnUMk/J3lcVT2mqh6Y5AVJ3rHsWACczb573N19Z1W9NMl7k1yQ5A3dfcvikwGwp7WOcXf3u5K8a+FZDtvYwzxnOCrrSI7OWqzj/mf0Wqq7D3sGAM6Dj7wDDHPkwl1VT6iqG3edvlxVrzhjm1ftuv7mqrqrqh62uu7lq8tu2f19VfUzq8vurqoDeTV6wbX8XlXdVlUfq6orq+qhQ9fx26s13FhVV1fV9y25jiXXsut7X1lVXVUXTVxHVb22qm7f9X3PnriO1XUvW/2c3FJVv7vkOs5bdx/ZU3ZeTP2PJI8+xzbPTfL+1ddPTnJzkgdl5/j/3yf5/tV1T0zyhCTXJNkevpZnJTm2+vp1SV43dB3fuet7fjnJn059TFbXX5ydNwF8OslFE9eR5LVJfvUgH4eF1vGM1flvXZ1/+GGs6WynI7fHfYZnJvn37v70ObZ5YZI3r75+YpJru/t/uvvOJB9M8lNJ0t23dvdhfqhok2u5enVZknwkO+/NPyibXMeXd33PhUkO+gWbja1l5Q+TvDrz13FYNrmOX0pyRXd/LUm6+46FZr5Pjnq4X5D/f5DupaoelOTSJG9bXXRzkh+tqu9eXffs3PPDR4dpqbX8QpJ3b3jWc9noOqrqd6rqM0lelOQ3Fpt6bxtbS1U9P8nt3X3TsiPvadPPrZeuDmG9oaq+a6mh97DJdTx+dd21VfXBqvqhBec+b0c23LXzYaHnJfmbc2z23CT/2N1fSnb2qrNz6ODqJO9JcmOSuxYedV9LraWqXpPkziRvWmDse1liHd39mu6+ODtreOlCo9/LJteyisav5+D/4VniMfmTJI9N8pQkn0vy+8tMfk8LrONYkocl+eEkr0ry1qqqZaY/f0c23Nn5NbQ3dPfnz7HNvf6F7u7Xd/cPdvePJfmvJP+24Izr2vhaqurFSZ6T5EW9Ooh3AJZ8TN6U5Kc3Nun+NrmWxyZ5TJKbqup0dg5d3VBV37vI5Pe00cekuz/f3Xd1991J/iw7v130IGz6ufXZJG/vHR9Ncnd2fr/J/cNhH2Rf6pSdXz/78+e4/iFJvpTkwjMuf/jqz+NJbkvy0DOuvyYH/OLkpteSnf8ufjzJ1vB1PG7XNi9L8rdT13LGNqdzQC9OLvCYPGLXNr+S5C1D1/GSJL+1+vrx2fkNqXVQz69913vYAyz0IF6Y5D+TPGTXZS9J8pJd51+815MqyT+sonZTkmfuuvwns/Ov8NeSfD7Jewev5ZOrJ+KNq9Pi78ZYaB1vy85xyo8l+bskj5z6mJyxzYGEe6HH5K+S/OvqMXlHdoV82DoemOSNq+fXDUkuOYjn1ronn5wEGOYoH+MGOJKEG2AY4QYYRrgBhhFugGGEG2AY4QYYRrgBhvk/dNOXw67IVb8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(np.array(crashing) * 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1. Experimental evaluation -> quite fast, doesn't ruin things much..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# computing probability delta of bad crash >= alpha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 181,
   "metadata": {},
   "outputs": [],
   "source": [
    "def q_factor(array):\n",
    "    \"\"\" q-balanced factor \"\"\"\n",
    "    array = np.abs(array)\n",
    "    assert np.max(array) > 0, \"Cannot compute q-factor for zero array\"\n",
    "        \n",
    "    return np.min(array) / np.max(array)\n",
    "\n",
    "q = q_factor(np.linalg.norm(exp.W[1], axis = 1) ** 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1189.2539"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(exp.W[1].flatten(), ord = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0001"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_inf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = N[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.02"
      ]
     },
     "execution_count": 178,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_inf * n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "metadata": {},
   "outputs": [],
   "source": [
    "W = np.linalg.norm(exp.W[1], ord = 1, axis = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 182,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.38271666"
      ]
     },
     "execution_count": 182,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "# target probability\n",
    "delta_0 = 1e-1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.log(2 / delta_0) / (2 * q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = (a/n)**0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = (a / n) ** 0.5\n",
    "alpha = t + p_inf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3.6787729679967103, 0.1918012765337267, 0.19190127653372668)"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a, t, alpha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "W = exp.model_correct.layers[1].weights[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 351,
   "metadata": {},
   "outputs": [],
   "source": [
    "sess = get_session()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sess.run(tf.hessians(exp.loss, W), feed_dict = )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 344,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.1181999267578125"
      ]
     },
     "execution_count": 344,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_inf * np.linalg.norm(exp.W[1].flatten(), ord = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 322,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.10010000000000001\n"
     ]
    }
   ],
   "source": [
    "t = p_inf * 1000\n",
    "alpha = (p_inf + t)\n",
    "print(alpha)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 312,
   "metadata": {},
   "outputs": [],
   "source": [
    "# target probability of big error\n",
    "delta_target = 1e-5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 317,
   "metadata": {},
   "outputs": [],
   "source": [
    "#t = np.log(1/delta_target) / (2 * n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 318,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 321,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.01"
      ]
     },
     "execution_count": 321,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n * p_inf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 310,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.061153622438558e-09"
      ]
     },
     "execution_count": 310,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Chernoff bound\n",
    "delta = np.exp(-2*t*n)\n",
    "delta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 193,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.10999999999999999"
      ]
     },
     "execution_count": 193,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alpha ** 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.009816558831670864\n",
      "1.228175672145885e-05\n"
     ]
    }
   ],
   "source": [
    "# for CHEBYSHEV\n",
    "print(np.mean(crashing))\n",
    "print(np.std(crashing))\n",
    "a = np.std(crashing) * 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 195,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.24644639434772372\n"
     ]
    }
   ],
   "source": [
    "print(np.std(crashing) ** 2 / a ** 2  + delta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {},
   "outputs": [],
   "source": [
    "def _get_bound_b5_loss(self, data, outputs, weights = None, add_diag_hess = False):\n",
    "    \"\"\" Exact error up to O(alpha^2) \"\"\"\n",
    "\n",
    "    # default value: first hidden layer\n",
    "    if weights is None:\n",
    "        weights = self.model_correct.layers[0].output\n",
    "        self.check_p_layer0()\n",
    "\n",
    "    def get_graph():\n",
    "        # resulting gradient w.r.t. first layer output\n",
    "        grad = []\n",
    "        grad_sq = []\n",
    "\n",
    "        # list of layers\n",
    "        layers = self.model_correct.layers\n",
    "\n",
    "        # loss tensor\n",
    "        loss = self.loss\n",
    "\n",
    "        # w.r.t. first layer output\n",
    "        if self.check_shape:\n",
    "            grad += [tf.reduce_sum(tf.multiply(tf.gradients([loss], [weights])[0], weights), axis=1)]\n",
    "            grad_sq += [tf.reduce_sum(tf.square(tf.multiply(tf.gradients(\n",
    "                [loss], [weights])[0], weights)), axis=1)]\n",
    "        else:\n",
    "            grad += [tf.reduce_sum(tf.multiply(tf.gradients([loss], [weights])[0], weights))]\n",
    "            grad_sq += [tf.reduce_sum(tf.square(tf.multiply(tf.gradients([loss], [weights])[0], weights)))]\n",
    "\n",
    "        # compute the result\n",
    "        p = np.max(p_inference)\n",
    "            \n",
    "        if add_diag_hess:\n",
    "            # squared weights\n",
    "            Wsq = tf.multiply(weights, weights)\n",
    "\n",
    "            wshape = weights.shape\n",
    "            if wshape[0].value is None:\n",
    "                wshape = wshape[1:]\n",
    "            \n",
    "            # diagonal of the hessian. TODO: compute in O(N) instead of O(N^2)\n",
    "            hess_diag = tf.linalg.tensor_diag_part(tf.reshape(tf.hessians([loss], weights)[0], (*wshape, *wshape)))\n",
    "\n",
    "            # quadratic term\n",
    "            S = -0.5 * tf.reduce_sum(tf.multiply(Wsq, hess_diag))\n",
    "            grad[0] += S\n",
    "        \n",
    "        return {'mean': tf.transpose(tf.multiply(-p, grad)), 'std': tf.transpose(tf.sqrt(tf.multiply(p, grad_sq)))}\n",
    "\n",
    "    return self.run_on_input_output(get_graph(), data, outputs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 196,
   "metadata": {},
   "outputs": [],
   "source": [
    "W1 = exp.model_correct.layers[1].output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 197,
   "metadata": {},
   "outputs": [],
   "source": [
    "# why is it negative and not equal to the experiment???\n",
    "theory = _get_bound_b5_loss(exp, exp.x_train, exp.y_train, weights = W1, add_diag_hess = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 200,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-3.4397557936293496e-05, 0.002791661591393622)"
      ]
     },
     "execution_count": 200,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(theory['mean']), np.mean(theory['std'] ** 2) ** 0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Trying for a single object...\n",
    "idx = 10\n",
    "x = exp.x_train[idx:idx+1]\n",
    "y = exp.y_train[idx:idx+1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 156,
   "metadata": {},
   "outputs": [],
   "source": [
    "crashing = exp.predict_crashing(x, repetitions = 1000)#.reshape(1000, 10)\n",
    "correct = exp.predict_correct(x)\n",
    "\n",
    "crashing = np.mean(np.square(crashing - y), axis = 2).reshape(-1)\n",
    "\n",
    "correct = np.mean(np.square(correct - y))\n",
    "\n",
    "delta = crashing - correct\n",
    "\n",
    "np.mean(delta), np.std(delta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'mean': array([[0.00183269]]), 'std': array([[0.000824]])}\n",
      "{'mean': array([[1.0128133e-07]]), 'std': array([[0.000824]])}\n"
     ]
    }
   ],
   "source": [
    "# WTF? here crashing error is less than correct...\n",
    "print(_get_bound_b5_loss(exp, x, y, weights = W1, add_diag_hess = True))\n",
    "print(_get_bound_b5_loss(exp, x, y, weights = W1, add_diag_hess = False))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exp.get_bound_b3(x)['mean']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tf.gradients(loss, W0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tensorflow_probability as tfp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g0 = tf.gradients(loss, W0)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tf.gradients(loss, W0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = tf.placeholder(tf.float32, shape = (10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loss = tf.reduce_sum(tf.square(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = tf.gradients(loss, x)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1 = tf.unstack(x)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = tf.square(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tf.gradients(x, x1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sess.run(tfp.math.diag_jacobian(g, x), feed_dict = {x: np.random.randn(10, 11)})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sess.run(tfp.math.diag_jacobian(x, x)[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sess.run(tfp.math.diag_jacobian(g0, W0, sample_shape = [1,1])[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g1 = tf.gradients(g0, W0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = sess.run(g1, feed_dict = {exp.model_correct.input: x, exp.output_tensor: y})[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Tail bound test\n",
    "# https://en.wikipedia.org/wiki/Binomial_distribution\n",
    "def kl_num(x, y):\n",
    "    return x * np.log(x / y) + (1 - x) * np.log((1 - x) / (1 - y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.45260976348657705"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kl_num(0.001,1e-200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see https://github.com/tensorflow/tensorflow/issues/25040\n",
    "\n",
    "import tensorflow as tf\n",
    "import numpy as np\n",
    "from tensorflow.python.ops import array_ops\n",
    "from tensorflow.python.ops import math_ops\n",
    "from tensorflow import gradients\n",
    "\n",
    "def hessian_vector_product(ys, xs, v):\n",
    "  \"\"\"Multiply the Hessian of `ys` wrt `xs` by `v`.\n",
    "  This is an efficient construction that uses a backprop-like approach\n",
    "  to compute the product between the Hessian and another vector. The\n",
    "  Hessian is usually too large to be explicitly computed or even\n",
    "  represented, but this method allows us to at least multiply by it\n",
    "  for the same big-O cost as backprop.\n",
    "  Implicit Hessian-vector products are the main practical, scalable way\n",
    "  of using second derivatives with neural networks. They allow us to\n",
    "  do things like construct Krylov subspaces and approximate conjugate\n",
    "  gradient descent.\n",
    "  Example: if `y` = 1/2 `x`^T A `x`, then `hessian_vector_product(y,\n",
    "  x, v)` will return an expression that evaluates to the same values\n",
    "  as (A + A.T) `v`.\n",
    "  Args:\n",
    "    ys: A scalar value, or a tensor or list of tensors to be summed to\n",
    "        yield a scalar.\n",
    "    xs: A list of tensors that we should construct the Hessian over.\n",
    "    v: A list of tensors, with the same shapes as xs, that we want to\n",
    "       multiply by the Hessian.\n",
    "  Returns:\n",
    "    A list of tensors (or if the list would be length 1, a single tensor)\n",
    "    containing the product between the Hessian and `v`.\n",
    "  Raises:\n",
    "    ValueError: `xs` and `v` have different length.\n",
    "  \"\"\" \n",
    "\n",
    "  # Validate the input\n",
    "  length = len(xs)\n",
    "  if len(v) != length:\n",
    "    raise ValueError(\"xs and v must have the same length.\")\n",
    "\n",
    "  # First backprop\n",
    "  grads = gradients(ys, xs)\n",
    "\n",
    "  # grads = xs\n",
    "\n",
    "  assert len(grads) == length\n",
    "\n",
    "  elemwise_products = [\n",
    "      math_ops.multiply(grad_elem, array_ops.stop_gradient(v_elem))\n",
    "      for grad_elem, v_elem in zip(grads, v) if grad_elem is not None\n",
    "  ]\n",
    "\n",
    "  # Second backprop  \n",
    "  grads_with_none = gradients(elemwise_products, xs)\n",
    "  return_grads = [\n",
    "      grad_elem if grad_elem is not None \\\n",
    "      else tf.zeros_like(x) \\\n",
    "      for x, grad_elem in zip(xs, grads_with_none)]\n",
    "  \n",
    "  return return_grads\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ws = exp.model_correct.weights[2:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<tf.Variable 'dense_9/kernel:0' shape=(784, 100) dtype=float32_ref>,\n",
       " <tf.Variable 'dense_9/bias:0' shape=(100,) dtype=float32_ref>,\n",
       " <tf.Variable 'dense_10/kernel:0' shape=(100, 100) dtype=float32_ref>,\n",
       " <tf.Variable 'dense_10/bias:0' shape=(100,) dtype=float32_ref>,\n",
       " <tf.Variable 'dense_11/kernel:0' shape=(100, 100) dtype=float32_ref>,\n",
       " <tf.Variable 'dense_11/bias:0' shape=(100,) dtype=float32_ref>,\n",
       " <tf.Variable 'dense_12/kernel:0' shape=(100, 10) dtype=float32_ref>,\n",
       " <tf.Variable 'dense_12/bias:0' shape=(10,) dtype=float32_ref>]"
      ]
     },
     "execution_count": 105,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Ws = "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [],
   "source": [
    "sess = get_session()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 300,
   "metadata": {},
   "outputs": [],
   "source": [
    "W = exp.model_correct.layers[1].weights[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 278,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Variable 'dense_41/kernel:0' shape=(100, 200) dtype=float32_ref>"
      ]
     },
     "execution_count": 278,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "W"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 279,
   "metadata": {},
   "outputs": [],
   "source": [
    "hvp = hessian_vector_product(exp.loss, [W], [tf.stop_gradient(W)])[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 280,
   "metadata": {},
   "outputs": [],
   "source": [
    "wHw = tf.reduce_sum(tf.multiply(tf.abs(hvp), tf.stop_gradient(tf.abs(W))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 281,
   "metadata": {},
   "outputs": [],
   "source": [
    "feed = {exp.model_correct.input: x, exp.output_tensor: y}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 188,
   "metadata": {},
   "outputs": [],
   "source": [
    "cert_den = sess.run(wHw, feed_dict = feed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "metadata": {},
   "outputs": [],
   "source": [
    "cert_enu = sess.run(tf.reduce_sum(tf.multiply(tf.abs(tf.gradients(exp.loss, W)[0]), tf.abs(W))), feed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.02951110303401947"
      ]
     },
     "execution_count": 174,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# This is alpha!\n",
    "cert_enu / cert_den / 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 183,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "200 0.38271666\n"
     ]
    }
   ],
   "source": [
    "print(n, q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 184,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.03 0.0301\n"
     ]
    }
   ],
   "source": [
    "t = 0.03\n",
    "print(t, t + p_inf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.8691927046718584e-05"
      ]
     },
     "execution_count": 185,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.exp(-n*q*kl_num(t + p_inf, p_inf))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.1894852870413295e-05, 0.00018028125517313703)"
      ]
     },
     "execution_count": 159,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "crashing = exp.predict_crashing(x, repetitions = 10000)#.reshape(1000, 10)\n",
    "correct = exp.predict_correct(x)\n",
    "\n",
    "crashing = np.mean(np.square(crashing - y), axis = 2).reshape(-1)\n",
    "\n",
    "correct = np.mean(np.square(correct - y))\n",
    "\n",
    "delta = crashing - correct\n",
    "\n",
    "np.mean(delta), np.std(delta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 187,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0001"
      ]
     },
     "execution_count": 187,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_inf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 214,
   "metadata": {},
   "outputs": [],
   "source": [
    "wHwg = tf.gradients(wHw, W)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 213,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(100, 200)"
      ]
     },
     "execution_count": 213,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sess.run(wHwg, feed)[0].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 284,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.21682757,  0.2660623 , -0.13696536, ...,  0.04763535,\n",
       "        -0.17063715,  0.00432912],\n",
       "       [-0.3852396 , -0.15810953, -0.01746978, ...,  0.09028391,\n",
       "        -0.04298741, -0.01647916],\n",
       "       [ 0.03244729, -0.11132674, -0.1511244 , ..., -0.08711679,\n",
       "         0.02930702, -0.07788351],\n",
       "       ...,\n",
       "       [ 0.00601823,  0.127277  , -0.150358  , ...,  0.07315686,\n",
       "         0.12544541, -0.0302185 ],\n",
       "       [ 0.08320521,  0.01221226,  0.00548319, ..., -0.09786451,\n",
       "        -0.09717263,  0.10645077],\n",
       "       [-0.32533476, -0.1721567 , -0.06282311, ...,  0.38471898,\n",
       "        -0.32876042, -0.28661048]], dtype=float32)"
      ]
     },
     "execution_count": 284,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "W1 = tf.Variable(tf.zeros_like(W))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 331,
   "metadata": {},
   "outputs": [],
   "source": [
    "hvp = hessian_vector_product(exp.loss, [W], [W1])[0]\n",
    "wHw = tf.reduce_sum(tf.multiply(tf.abs(hvp), tf.stop_gradient(tf.abs(W1))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 332,
   "metadata": {},
   "outputs": [],
   "source": [
    "a_op = tf.train.AdamOptimizer(0.1).minimize(-wHw)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 333,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.21682757,  0.2660623 , -0.13696536, ...,  0.04763535,\n",
       "        -0.17063715,  0.00432912],\n",
       "       [-0.3852396 , -0.15810953, -0.01746978, ...,  0.09028391,\n",
       "        -0.04298741, -0.01647916],\n",
       "       [ 0.03244729, -0.11132674, -0.1511244 , ..., -0.08711679,\n",
       "         0.02930702, -0.07788351],\n",
       "       ...,\n",
       "       [ 0.00601823,  0.127277  , -0.150358  , ...,  0.07315686,\n",
       "         0.12544541, -0.0302185 ],\n",
       "       [ 0.08320521,  0.01221226,  0.00548319, ..., -0.09786451,\n",
       "        -0.09717263,  0.10645077],\n",
       "       [-0.32533476, -0.1721567 , -0.06282311, ...,  0.38471898,\n",
       "        -0.32876042, -0.28661048]], dtype=float32)"
      ]
     },
     "execution_count": 333,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sess.run(tf.global_variables_initializer())\n",
    "sess.run(W1.assign(W))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 356,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Variable 'dense_41/kernel:0' shape=(100, 200) dtype=float32_ref>"
      ]
     },
     "execution_count": 356,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "W"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 355,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor 'sub_4:0' shape=(100, 200) dtype=float32>"
      ]
     },
     "execution_count": 355,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "W + (W - W1) / tf.norm(W - W1) * eps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 342,
   "metadata": {},
   "outputs": [],
   "source": [
    "nrm = tf.norm(W)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 353,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.8591837e-07, None]"
      ]
     },
     "execution_count": 353,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sess.run([wHw / nrm ** 2, a_op], feed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 259,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15.261283"
      ]
     },
     "execution_count": 259,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(exp.W[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tf.zeros(W.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 305,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<tf.Tensor 'gradients_76/dense_41/MatMul_grad/MatMul_1:0' shape=(100, 200) dtype=float32>]"
      ]
     },
     "execution_count": 305,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hvp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 283,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.21682757,  0.2660623 , -0.13696536, ...,  0.04763535,\n",
       "        -0.17063715,  0.00432912],\n",
       "       [-0.3852396 , -0.15810953, -0.01746978, ...,  0.09028391,\n",
       "        -0.04298741, -0.01647916],\n",
       "       [ 0.03244729, -0.11132674, -0.1511244 , ..., -0.08711679,\n",
       "         0.02930702, -0.07788351],\n",
       "       ...,\n",
       "       [ 0.00601823,  0.127277  , -0.150358  , ...,  0.07315686,\n",
       "         0.12544541, -0.0302185 ],\n",
       "       [ 0.08320521,  0.01221226,  0.00548319, ..., -0.09786451,\n",
       "        -0.09717263,  0.10645077],\n",
       "       [-0.32533476, -0.1721567 , -0.06282311, ...,  0.38471898,\n",
       "        -0.32876042, -0.28661048]], dtype=float32)"
      ]
     },
     "execution_count": 283,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 365,
   "metadata": {},
   "outputs": [],
   "source": [
    "feed = {exp.model_correct.input: x, exp.output_tensor: y}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 367,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sess.run(tf.gradients(exp.loss, exp.model_correct.input)[0], feed_dict = feed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 370,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = np.abs(g.flatten())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 371,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([237., 191., 144.,  94.,  57.,  35.,  15.,   8.,   1.,   2.]),\n",
       " array([1.0964959e-07, 1.2013664e-04, 2.4016363e-04, 3.6019064e-04,\n",
       "        4.8021763e-04, 6.0024462e-04, 7.2027161e-04, 8.4029860e-04,\n",
       "        9.6032559e-04, 1.0803526e-03, 1.2003796e-03], dtype=float32),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 371,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAADixJREFUeJzt3X2MXNV5x/HvE5zQ9/LijeUa6ELiSjVSQ9CKoraqaJESMFVN/ggFVYlLiFwloL5K1aZIpYqEZIKaSlQNyFFcTEUBtwnCkp0XarWiUQthiQgxEMKGGGHXYAMtDYqa1uTpH3tMx2Z3Z3ZePLuPvh9pNHfOnHvOeTy7P9+9d3Y2MhNJUl1vG/cCJEmjZdBLUnEGvSQVZ9BLUnEGvSQVZ9BLUnEGvSQVZ9BLUnEGvSQVt2rcCwBYvXp1Tk5OjnsZkrSiPPbYYy9n5kS3fssi6CcnJ5mZmRn3MiRpRYmI53vp56kbSSrOoJek4gx6SSrOoJek4gx6SSrOoJek4gx6SSrOoJek4gx6SSpuWfxm7CAmp3ePbe79W68Y29yS1CuP6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpuK5BHxFnR8Q/RcRTEfFkRPx+az8jIh6MiGfb/emtPSLitoiYjYgnIuLCURchSVpYL0f0R4E/zswNwMXA9RGxAZgG9mbmemBvewxwObC+3bYAtw991ZKknnUN+sw8lJlfb9vfA54G1gGbgB2t2w7gyra9Cbgr5zwMnBYRa4e+cklST5Z0jj4iJoH3Ao8AazLzUHvqRWBN214HvNCx24HWJkkag56DPiJ+Avg88AeZ+V+dz2VmArmUiSNiS0TMRMTMkSNHlrKrJGkJegr6iHg7cyF/d2Z+oTW/dOyUTLs/3NoPAmd37H5WaztOZm7LzKnMnJqYmOh3/ZKkLnp5100AnwOezsxPdzy1C9jctjcDD3S0f7i9++Zi4LWOUzySpJNsVQ99fhn4EPDNiHi8tf0psBXYGRHXAc8DV7Xn9gAbgVng+8C1Q12xJGlJugZ9Zn4ViAWevnSe/glcP+C6JElD4m/GSlJxBr0kFWfQS1JxvVyM1QImp3ePZd79W68Yy7ySViaP6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpOINekooz6CWpuK5BHxHbI+JwROzraPvziDgYEY+328aO5z4REbMR8UxEvH9UC5ck9aaXI/o7gcvmaf/LzLyg3fYARMQG4Grg/LbPZyLilGEtVpK0dF2DPjMfAl7tcbxNwL2Z+YPM/C4wC1w0wPokSQMa5Bz9DRHxRDu1c3prWwe80NHnQGuTJI1Jv0F/O/Au4ALgEPAXSx0gIrZExExEzBw5cqTPZUiSuukr6DPzpcx8IzN/CHyW/z89cxA4u6PrWa1tvjG2ZeZUZk5NTEz0swxJUg/6CvqIWNvx8APAsXfk7AKujohTI+JcYD3wtcGWKEkaxKpuHSLiHuASYHVEHABuAi6JiAuABPYDvwuQmU9GxE7gKeAocH1mvjGapUuSetE16DPzmnmaP7dI/5uBmwdZlCRpeLoGvZafyendY5t7/9Yrxja3pP74EQiSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFGfSSVJxBL0nFdQ36iNgeEYcjYl9H2xkR8WBEPNvuT2/tERG3RcRsRDwREReOcvGSpO56OaK/E7jshLZpYG9mrgf2tscAlwPr220LcPtwlilJ6lfXoM/Mh4BXT2jeBOxo2zuAKzva78o5DwOnRcTaYS1WkrR0/Z6jX5OZh9r2i8Catr0OeKGj34HW9hYRsSUiZiJi5siRI30uQ5LUzcAXYzMzgexjv22ZOZWZUxMTE4MuQ5K0gH6D/qVjp2Ta/eHWfhA4u6PfWa1NkjQm/Qb9LmBz294MPNDR/uH27puLgdc6TvFIksZgVbcOEXEPcAmwOiIOADcBW4GdEXEd8DxwVeu+B9gIzALfB64dwZolSUvQNegz85oFnrp0nr4JXD/ooiRJw9M16KVOk9O7xzLv/q1XjGVeqQI/AkGSijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16SijPoJak4g16Sils17gVIvZic3j22ufdvvWJsc0vD4BG9JBU30BF9ROwHvge8ARzNzKmIOAO4D5gE9gNXZeZ/DLZMSVK/hnFE/2uZeUFmTrXH08DezFwP7G2PJUljMopTN5uAHW17B3DlCOaQJPVo0KBP4CsR8VhEbGltazLzUNt+EVgz4BySpAEM+q6bX8nMgxHxTuDBiPhW55OZmRGR8+3Y/mPYAnDOOecMuAxJ0kIGOqLPzIPt/jBwP3AR8FJErAVo94cX2HdbZk5l5tTExMQgy5AkLaLvoI+IH4+Inzy2DbwP2AfsAja3bpuBBwZdpCSpf4OculkD3B8Rx8b5u8z8UkQ8CuyMiOuA54GrBl+mJKlffQd9Zj4HvGee9leASwdZlCRpePzNWEkqzqCXpOIMekkqzqCXpOIMekkqzqCXpOIMekkqzqCXpOL8U4JSF+P6M4b+CUMNi0f0klScQS9JxRn0klScQS9JxRn0klScQS9JxRn0klScQS9JxRn0klScQS9JxRn0klScQS9JxRn0klScQS9JxRn0klScQS9JxRn0klScf2FKWqb8y1YaFo/oJak4g16SijPoJak4g16SijPoJak433Uj6TjjercP+I6fUfGIXpKKM+glqbiRBX1EXBYRz0TEbERMj2oeSdLiRhL0EXEK8NfA5cAG4JqI2DCKuSRJixvVxdiLgNnMfA4gIu4FNgFPjWg+Sepb9QvQowr6dcALHY8PAL84orkkFTHOwK1sbG+vjIgtwJb28PWIeKbPoVYDLw9nVWNnLctTlVqq1AGFaolbBqrlZ3vpNKqgPwic3fH4rNb2pszcBmwbdKKImMnMqUHHWQ6sZXmqUkuVOsBalmpU77p5FFgfEedGxDuAq4FdI5pLkrSIkRzRZ+bRiLgB+DJwCrA9M58cxVySpMWN7Bx9Zu4B9oxq/A4Dn/5ZRqxleapSS5U6wFqWJDJz1HNIksbIj0CQpOKWRdB3+7iEiDg1Iu5rzz8SEZMdz32itT8TEe/vNma7QPxIa7+vXSxeiXXc3dr3RcT2iHj7sOo42bV0PH9bRLw+zDpOdi0x5+aI+HZEPB0Rv7eCa7k0Ir4eEY9HxFcj4t3LvI7tEXE4IvadMNYZEfFgRDzb7k8fVh1jqOXWiPhWRDwREfdHxGk9LTIzx3pj7mLtd4DzgHcA3wA2nNDn48Adbftq4L62vaH1PxU4t41zymJjAjuBq9v2HcDHVmgdG4Fot3uGVcc4amn7TQF/C7y+wr++rgXuAt7WHr9zBdfybeDnO8a9c7nW0Z77VeBCYN8JY30KmG7b08Aty/k16VLL+4BVbfuWXmtZDkf0b35cQmb+D3Ds4xI6bQJ2tO1/AC6NiGjt92bmDzLzu8BsG2/eMds+v97GoI155UqrA+YudmcDfI2531UYlpNaS8x9NtKtwJ8MsYax1AJ8DPhkZv4QIDMPr+BaEviptv3TwL8v4zrIzIeAV+eZr3OsYX7Pn/RaMvMrmXm0PXyYHr/vl0PQz/dxCesW6tOKfA04c5F9F2o/E/jPjn+o+ebq18ms403tlM2HgC8NXME861xoXoZbyw3Arsw8NKT1z7vOeeZ9S58h1PIu4LciYiYivhgR64dUx3HrnGfet/QZQi0fBfZExAHmvsa2DqWK0dSxmDUdX1svAmv6W/a8TnYtnT4CfLGXjssh6DWYzwAPZea/jHsh/YiInwE+CPzVuNcyJKcC/51zv+n4WWD7mNcziD8ENmbmWcDfAJ8e83oG1n4CXvFvNYyIG4GjwN299F8OQd/14xI6+0TEKuZ+jHxlkX0Xan8FOK2NsdBc/TqZddDGuAmYAP5oKBXMs8755j2xz4C1vBd4NzAbEfuBH4uI2WEVssi88/YZwutyAPhC274f+IWBK5hnnfPM+5Y+g9QSERPAezLzkdZ+H/BLwyljJHUs5qWIWNvGWgsM83Taya6FiPgd4DeA327/cXU3rIsSA1zMWAU8x9zFiGMXM84/oc/1HH8xY2fbPp/jL2Y8x9zFkQXHBP6e4y/GfnyF1vFR4F+BH13pr8kJ4w77YuzJfl22Ah9p25cAj67EWlr7y8DPtf2vAz6/XOvo2G+St17AvJXjL8Z+ajm/Jl1quYy5j3ufWNI6h/lNNcA/1kbmrvB/B7ixtX0S+M22/SPMBfQscxcez+vY98a23zPA5YuN2drPa2PMtjFPXaF1HG1tj7fbn63U1+SEeYca9GN4XU4DdgPfBP6NuaPilVrLB1od3wD+uXOsZVrHPcAh4H+Z+8nqutZ+JrAXeBb4R+CMFfCaLFTLLHPn9Y9939/Ryxr9zVhJKm45nKOXJI2QQS9JxRn0klScQS9JxRn0klScQS9JxRn0klScQS9Jxf0f8iNFnd4aOvgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 372,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.579237"
      ]
     },
     "execution_count": 372,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.max(g) / np.mean(g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.8"
  },
  "notify_time": "5",
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "position": {
    "height": "407px",
    "left": "1561px",
    "right": "20px",
    "top": "120px",
    "width": "339px"
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}