TailBound-FC-MNIST.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-20T20:13:59.673640Z",
     "start_time": "2019-01-20T20:13:57.997310Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Using TensorFlow backend.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initialized TensorFlow\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": 26,
   "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, 100]\n",
    "\n",
    "def get_p_arr(p):\n",
    "    \"\"\" p array with failure on the first layer \"\"\"\n",
    "    return [0, p]  + [0] * (len(N))\n",
    "\n",
    "# prob. of failure at inference\n",
    "p_inf = 1e-5\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": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "_________________________________________________________________\n",
      "Layer (type)                 Output Shape              Param #   \n",
      "=================================================================\n",
      "dense_1 (Dense)              (None, 100)               78500     \n",
      "_________________________________________________________________\n",
      "dense_2 (Dense)              (None, 100)               10100     \n",
      "_________________________________________________________________\n",
      "dense_3 (Dense)              (None, 10)                1010      \n",
      "=================================================================\n",
      "Total params: 89,610\n",
      "Trainable params: 89,610\n",
      "Non-trainable params: 0\n",
      "_________________________________________________________________\n",
      "Train on 60000 samples, validate on 10000 samples\n",
      "Epoch 1/10\n",
      "60000/60000 [==============================] - 1s 20us/step - loss: 0.2449 - categorical_accuracy: 0.8003 - mean_squared_error: 0.0473 - mean_absolute_error: 0.1386 - val_loss: 0.1097 - val_categorical_accuracy: 0.8862 - val_mean_squared_error: 0.0305 - val_mean_absolute_error: 0.1188\n",
      "Epoch 2/10\n",
      "60000/60000 [==============================] - 1s 15us/step - loss: 0.0777 - categorical_accuracy: 0.8873 - mean_squared_error: 0.0275 - mean_absolute_error: 0.1068 - val_loss: 0.0576 - val_categorical_accuracy: 0.9006 - val_mean_squared_error: 0.0232 - val_mean_absolute_error: 0.0916\n",
      "Epoch 3/10\n",
      "60000/60000 [==============================] - 1s 15us/step - loss: 0.0552 - categorical_accuracy: 0.8923 - mean_squared_error: 0.0259 - mean_absolute_error: 0.1014 - val_loss: 0.0504 - val_categorical_accuracy: 0.9021 - val_mean_squared_error: 0.0245 - val_mean_absolute_error: 0.0974\n",
      "Epoch 4/10\n",
      "60000/60000 [==============================] - 1s 15us/step - loss: 0.0495 - categorical_accuracy: 0.8955 - mean_squared_error: 0.0251 - mean_absolute_error: 0.0985 - val_loss: 0.0466 - val_categorical_accuracy: 0.9023 - val_mean_squared_error: 0.0233 - val_mean_absolute_error: 0.0931\n",
      "Epoch 5/10\n",
      "60000/60000 [==============================] - 1s 15us/step - loss: 0.0469 - categorical_accuracy: 0.8972 - mean_squared_error: 0.0243 - mean_absolute_error: 0.0960 - val_loss: 0.0435 - val_categorical_accuracy: 0.9064 - val_mean_squared_error: 0.0215 - val_mean_absolute_error: 0.0827\n",
      "Epoch 6/10\n",
      "60000/60000 [==============================] - 1s 15us/step - loss: 0.0453 - categorical_accuracy: 0.8985 - mean_squared_error: 0.0237 - mean_absolute_error: 0.0930 - val_loss: 0.0421 - val_categorical_accuracy: 0.9141 - val_mean_squared_error: 0.0210 - val_mean_absolute_error: 0.0862\n",
      "Epoch 7/10\n",
      "60000/60000 [==============================] - 1s 15us/step - loss: 0.0442 - categorical_accuracy: 0.8997 - mean_squared_error: 0.0232 - mean_absolute_error: 0.0914 - val_loss: 0.0422 - val_categorical_accuracy: 0.9129 - val_mean_squared_error: 0.0215 - val_mean_absolute_error: 0.0881\n",
      "Epoch 8/10\n",
      "60000/60000 [==============================] - 1s 15us/step - loss: 0.0435 - categorical_accuracy: 0.9003 - mean_squared_error: 0.0229 - mean_absolute_error: 0.0898 - val_loss: 0.0438 - val_categorical_accuracy: 0.9027 - val_mean_squared_error: 0.0233 - val_mean_absolute_error: 0.0973\n",
      "Epoch 9/10\n",
      "60000/60000 [==============================] - 1s 15us/step - loss: 0.0428 - categorical_accuracy: 0.9023 - mean_squared_error: 0.0225 - mean_absolute_error: 0.0884 - val_loss: 0.0430 - val_categorical_accuracy: 0.9040 - val_mean_squared_error: 0.0229 - val_mean_absolute_error: 0.0910\n",
      "Epoch 10/10\n",
      "60000/60000 [==============================] - 1s 15us/step - loss: 0.0424 - categorical_accuracy: 0.9021 - mean_squared_error: 0.0223 - mean_absolute_error: 0.0872 - val_loss: 0.0428 - val_categorical_accuracy: 0.9031 - val_mean_squared_error: 0.0228 - val_mean_absolute_error: 0.0919\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD9CAYAAABQvqc9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VPW9//HXJ5N9JZCwJggqyh6QVVHphaKoLdQFAbdiXa63gsjV24uVW9f+rlXrtVarxYqKUhB32+IOii1IEnZZBURmkggBspCELDPz/f1xJmESskySSSaZ+Twfj3nMmbPMfGbE9zn5nu/5HjHGoJRSKjSEBboApZRS7UdDXymlQoiGvlJKhRANfaWUCiEa+kopFUI09JVSKoT4FPoiMlVE9ojIPhFZWM/yM0TkcxHZJiJfiEiaZ/4IEVkvIjs8y2b6+wsopZTynTTVT19EbMBeYArgALKA2caYnV7rvAn83RjzqohMAm42xtwoIucAxhjzrYj0BjYCg4wxhW30fZRSSjXClyP9scA+Y8wBY0wlsAKYXmedwcBqz/Sa6uXGmL3GmG8907nAESDVH4UrpZRqPl9Cvw9g93rt8MzzthW4yjN9JZAgIt28VxCRsUAksL9lpSqllGqtcD+9z73AsyIyB1gL5ACu6oUi0gt4Dfi5McZdd2MRuR24HSAuLm7UwIED/VSWUkqFho0bNx41xjTZkuJL6OcA6V6v0zzzaniabq4CEJF44OrqdnsRSQT+AdxvjPm6vg8wxiwGFgOMHj3aZGdn+1CWUkqpaiLyvS/r+dK8kwUMEJH+IhIJzAI+qPNhKSJS/V73AUs88yOBd4Glxpi3fC1eKaVU22gy9I0xTmAu8DGwC1hpjNkhIg+LyDTPaj8C9ojIXqAH8FvP/GuBi4E5IrLF8xjh7y+hlFLKN0122Wxv2ryjlFLNJyIbjTGjm1pPr8hVSqkQoqGvlFIhRENfKaVCiIa+UkqFEH9dnKWUakL+iQr+sS2Xk1VuEmPCSYyOIDEmgsTocBKiI2rmRUfYAl2qCmIa+kq1Ibfb8NW+o6zIPMSnOw/jdDfdWy4yPMzaIUSHk+DZKVTvHLx3FIkxESTUmmftOGIibIhIO3w71Rlp6CvVBn4oKufNbDsrsuzkFJ4kOTaCmyf0Y+aYvvTpEkNxeRXFJ6us53KnZ9pZM+9EnXm5hSdrpiucp41kUostTLx2FN47hlM7iOTYCC4b1ouU+Kh2+kVUR6Ghr5SfOF1uvtybz/LMQ6zefQS3gQlnd2PhZQO5ZEgPosJPNdvERNrokRjdos+pcLpO2ynU3lFUUXzSyQmvHcqBoyU1y0srrWGx/vfD3dxyYX9uvehMkmIi/PIbqI5PQ1+pVnIUlLEyy87KbAc/FJeTEh/Fv088i1lj0jmjW5zfPy8q3EZUvK3FR+lOl5sDR0v5w+ff8sfV+1i6/nvumHgWcy7oR0yknk8IdnpFrlItUOVy8/muwyzPtLP223wALh6QyuyxfZk8qDsRts7RMe6bnCJ+/8ke1uzJJzUhirsmnc3MMX2JDO8c9atTfL0iV0NfqWY4eLSUFVl23tro4GhJBT0To7l2TDrXjk4jLTk20OW1WNbB4zzx0R4yDx4nLTmGBT8+h5+N7IMtTE8IdxYa+kr5SYXTxcc7DrMi8xDr9h/DFib827ndmT02nYnnpBLeSY7qm2KMYe23R3ni4918k1PMgO7x3HPJuVw6pIf2BuoENPSVaqV9R0pYkXmItzc5KCirIi05hllj0pkxOr3FJ2E7A2MMH37zA7//ZA/780sZnpbEf116LheenaLh34Fp6CvVAuVVLlZtz2N55iGyDhYQHiZcMqQHs8b05cKzUwgLoeYOp8vNu5tzePqzb8kpPMn4M7vyX5cOZNQZyYEuTdVDQ1+pZtiVV8yKzEO8uzmH4nIn/brFMmtsX64+L43UhNDuy17hdLEi084fV+/jaEkFPx7UnXsuOZdBvRIDXZryoqGvVBNKK5z8fVsuyzPtbLEXEmkLY+rQnswe25fxZ3bVpow6yiqdvPyvg/z5y/2cqHDy0+G9WTDlHPqn+L9bqmo+DX2lGrDdUcRfMw/xt625lFQ4Obt7PLPH9uWqkX1IjosMdHkdXlFZFYu/2s+Sfx6k0uXm2tFp3DV5AL2SYgJdWkjT0FfKS1FZFR9sy2VF5iF25BYTHRHGFcN6M3tsOqPOSNaj+hY4cqKcP63Zz7IN3yMi3DT+DP7jR2fRTYd2CAgNfRWyKpwuduWdYKu9kK32QrY4CjmQXwrAoF6JzB6bzvQRfXToAT9xFJTxh8++5e1NDmIibNxy0ZncelF/EqP1921PGvoqJLjdhgNHS62Ad1ghvzOvmCqX9e86JT6KEeldyEhL4uJzUhmelqRH9W1k35ES/u/Tvfxjex5dYiP4j4ln8fML+ulQ0e1EQ18FpcPF5WzxHMFvdRSyzV7EiQonAHGRNoalJZGR3oURaV3ISO9Cr6RoDfl29k1OEU98vIcv9+bTIzGKeZMGMHNMeqcZmqKz0tBXnV5xeRXbHUW1Qv5wcQUA4WHCwF4JZHjCfUR6F85KjddhAzqQDQeO8cTHe8j+voC+XWNZMGUA0zJ0aIe2oqGvOpUKp4vdeSfY6iisCfn9nnZ4gP4pcQxPS6oJ+SG9E7XZoBMwxvDF3nye+GgPO/OKObdHAvdccg5TBvt/aAdjDCerXBSdrKLoZBWFZVU108UnT02f9iirotLpJibSRlxUOLGRNuIiw4mNshEbaSM2Mpy4SBuxUZ7nSGsd79dxUafmV2/b3n/ZaOirDsvtNnx3rNTrRGsRu3KLqXRZNwdJiY/0tMNbAT88LYkusdqVsjNzuw2rvsnjqU/2cuBoKRnpXfjVpecy4eyUWuvVDe6istoBXR3ehSfrn199Lqc+IpAYHUFSjPXoEmvdUCYpJoKo8DBOVroorXRxstJJaYWLskonpZUuyio8z5XORt+/rkhbGLFR1k4gJtJ22g4iLspGTETt172SYpgyuEeLfmMNfdVhHKluh3cUstVexFZHISfKrXb42Egbw/okWSHvefTWdvig5XS5eWdTDk9/tpfconKG9Ukiwiae8LZu8lK9869P3eCufiR6QrzufO/lCVHh9Q+jYQy4KkHCQGzWhzTw76/S6fbsHJzWTqHCmq7eYdTsILx2FGWe55odSYX3fGv76hg+r28X3vnlhBb9tr6Gvt5EpZMyxlBS4aSgtIrjZZUUlFVSUFrJ8VLPdJl19OM2Bpfb4DbW0ZbLeE27DW5T/aDWa5fbWsdtrG2MZ7nLbTCeeS43XtOn1jn1Htb7VrOFCQN7JvDTjN41J1rP7q7t8KEk3BbGtWPSmTaiN8szD/H3bXnERFhHuIl1grpuiDca3A2pLIMTeXA4D078YE0X51nPJ36AE7nWs7O89nbVO4Awm9dzGJFhNiIljKRay8JqrXP6PK9l0WEQa6u1zISF4TZhOI3gTD4LaFno+8qn0BeRqcAfABvwF2PMY3WWnwEsAVKB48ANxhiHZ9nPgUWeVR81xrzqp9qDRnWAF5ZVcby00grxUiu4C2q9rqwJ+cKyygb/1LSFCcmxESRER2ALE2wihIUJYWItExFsXtPhYWGeaWtemFgPWxjWdPV7CF7Ttd+z1jae19Wf2yUmgoz0JIb0TtJ2eAVAdISNmyf05+YJ/Vv2Bq4qKDlcf4DXBPsPUFF0+rbhMZDYCxJ6QZ9R1nNMFzCAcYHbVefZDcZdzzLP/Oas73aCs6LWMnG7sRk3NuMiqh1uXtNk6IuIDXgOmAI4gCwR+cAYs9NrtSeBpcaYV0VkEvC/wI0i0hV4ABiN9ZNu9Gxb4O8v0pGUVTo5VuJ91F3J8dIqCstOzTteWlkT8gU+BHiX2Ei6xkbSLyWWkbFdSI6zXifHRZIcG3HqdWwkCdHNPBpSqqNwu6HsWD0BXifYS49iRYqXsHCI7wkJPSFlAJw50ZpO6OX16AnRSQ0234QCX470xwL7jDEHAERkBTAd8A79wcB/eqbXAO95pi8FPjXGHPds+ykwFVje+tI7Hpfb8IfPv+W5NftwuU8P8TCBZE9Qd42NpG/XWEaknwrwLrERdI2rDnJrnga46pDcbnBVWM0iTl+eG1hWWQolP1hBXpxnTbudp39eXOqpAO993qkAT+x9an5sitWsohrlS+j3Aexerx3AuDrrbAWuwmoCuhJIEJFuDWzbp8XVdmA/FJUzf8VmNnx3nOkjenPh2Sl0jYu0jtDjNMBVO3NVQXkxlBdCedGpR0WxZ7oYnCcbCerKxoPcXdX6Gm2REBFz6ui8/0X1H5nH94Bw7b3lL/46kXsv8KyIzAHWAjmAy9eNReR24HaAvn37+qmk9rNm9xHueXMr5VUunro2g6vOSwt0SaqzqyqvE9R1wru8uM5r70AvgqqyJj5AIDwawqMafo7t5nndyDrVzzbv9Rpa12sdW5QelQeIL6GfA6R7vU7zzKthjMnFOtJHROKBq40xhSKSA/yozrZf1P0AY8xiYDFYXTZ9Lz+wKp1unvxkD4vXHmBgzwSeu/48zkqND3RZqiNwu60QPllw+qOhoPYOdFdF4+8fFm61TUclWs/RSafaq+t7eK8XnQSR8Rq6IcqX0M8CBohIf6ywnwVc572CiKQAx40xbuA+rJ48AB8D/09Equ+vdolneadnP17G3OWb2Wov5MbxZ3D/FYO0Z0owcjmtIK4vvBt7lBdavTUaYouE6C61g7hL33qCus460Z7wjogN6ZORquWaDH1jjFNE5mIFuA1YYozZISIPA9nGmA+wjub/V0QMVvPOnZ5tj4vII1g7DoCHq0/qdmYfbs/jV29vA+D568/jsmG9AlyRapKrCsqONzO8C+vv8uctOskK5phk69Gl76npeh9drPUjgvfG6qpj0ytym6G8ysWj/9jJ618fIiO9C8/OHkl619hAlxWaXE6ra19pPpQdtbrwlR5t+HV5I+EtYU0EdQOPqESw6fWNqmPQK3L9bH9+CXP/upldecXcfvGZ3HvJuUS2w4UUIcPlhJPH6w/uMs+8Uq+QP9nApR5is05AxqVCXDfoPcLqyheXArFd6w/vyARt31YhQ0PfB+9scrDovW+ICg/j5Tlj+LeB3QNdUudgDBTZoeD7po/GTxZw2sU24DkK7+oJ8RToOdQT4p5Qj0v1ep1iNZ1ogCvVIA39RpRWOPnN+zt4e5ODsf278syskfRM0rbY0xgDRQ7I3wP5u+DIbus5fw9UltRZWayj6+qQ7j4QYi889ToupXaIxyRbY5QopfxCQ78Bu/KKmfvXTRw4Wspdkwdw16SzCQ/1O/8YA8W5tYP9yG5PuJ84tV5cdyvMR1xvPXc905oXl2IdtWs7uFIBo//31WGMYdmGQzz89510iYlg2a3juOCslKY3DCbGWJfF1xfu3r1Z4lIhdSBkzLLCPXUQdB9ktZ0rpTokDX0vxeVV3Pf2dv6xPY+Lz0nlqWszSImPCnRZbccYKDlST7jvqt3bJaarFebDrrGeuw+ywj4uxHaGSgUBDX2PrfZC5i7fRG5hOQsvG8jtF50ZXOPklOTXH+7evWBikq2j9aFXe47aPUfvcSl6IZBSQSLkQ98Yw0v//I7ffbSb7gnRrPz38xl1RnLTG3Zkxbnw3VfgyDwV7mXHTi2PTrLCfPD02uEe313DXakgF9Khf7y0knvf3Mrq3Ue4ZHAPHr9meOe8F2vJETj4lRX0362F4/ut+ZEJVlPMwCtqh3tCTw13pUJUyIZ+5nfHuWv5Zo6XVvLgTwfz8wv6dZ77spYdPxXyB7+C/N3W/MgEOOMCGP0La5jaHkO1u6NSqpaQC32X2/CnNfv4v8/20rdrLO/88gKG9kkKdFmNO1kI36+zjuIPfgWHv7HmR8RB3/FW75l+F0OvDO0OqZRqVEglxJET5Sx4Ywv/2neM6SN689srhxEf1QF/gooT8P16OLjWOpr/YZs1YmN4NKSPg0mLrJDvcx7YIgJdrVKqE+mAidc2vvo2nwVvbKGkwsnjVw9nxui0jtOcU1kKh74+1WSTu9m6cbItEtLGwMW/sppr0sZYN6FQSqkWCvrQd7rcPPXpXp7/cj8Dusfz19vGc06PhMAWVVVu9aypbpN3ZFu3nwsLhz6j4MIFnpAfC5E6iqdSyn+COvRzCk9y1/LNbPy+gNlj0/nNT4YQExmAE5vOSsjJPhXy9kzrzkgSBr1GwPm/tJpr+o6HKL3zllKq7QRt6H+68zD3vrkVl9vwzOyRTMvo3X4f7nJaTTTVbfKHvrZuQo1Az2Ew9jbodxGccb7VZ14ppdpJ0IV+hdPFYx/u5uV/HWRon0SenX0e/VLi2q+A/L3w0o9PDWPQfTCcd5PVXHPGBB2XRikVUEEV+gePljJ3+Sa+ySnm5gn9WHjZQKLC27k5Z99nVuBf+Wc4+8c6Po1SqkMJmtA/kF/CtGf/hS1MWHzjKC4Z0jMwhdg3QFK61XdeKaU6mKAJ/f4pcdw8oR+zxvalT5eYwBXiyLJOyCqlVAcUNKEvItxzybmBLaLIAcU5VldLpZTqgEL8VlB+Zs+0ntPHBLYOpZRqgIa+PzmyIDwGeg4PdCVKKVUvDX1/smdC75E6Ho5SqsPS0PeXqnLI26pNO0qpDk1D31/ytljj56SPC3QlSinVIJ9CX0SmisgeEdknIgvrWd5XRNaIyGYR2SYil3vmR4jIqyKyXUR2ich9/v4CHUb1SVztuaOU6sCaDH0RsQHPAZcBg4HZIjK4zmqLgJXGmJHALOBPnvkzgChjzDBgFPDvItLPP6V3MI5MSO4H8amBrkQppRrky5H+WGCfMeaAMaYSWAFMr7OOARI900lArtf8OBEJB2KASqC41VV3NMZYR/p6lK+U6uB8Cf0+gN3rtcMzz9uDwA0i4gBWAfM8898CSoE84BDwpDHmeGsK7pAKD0HJYUjX0FdKdWz+OpE7G3jFGJMGXA68JiJhWH8luIDeQH/gHhE5s+7GInK7iGSLSHZ+fr6fSmpHjizrWUNfKdXB+RL6OUC61+s0zzxvtwArAYwx64FoIAW4DvjIGFNljDkC/AsYXfcDjDGLjTGjjTGjU1M7YZu4fYN1k/LuQwJdiVJKNcqX0M8CBohIfxGJxDpR+0GddQ4BkwFEZBBW6Od75k/yzI8DxgO7/VN6B2LP9NykPGiGMlJKBakmQ98Y4wTmAh8Du7B66ewQkYdFZJpntXuA20RkK7AcmGOMMVi9fuJFZAfWzuNlY8y2tvgiAVNZBoe/0aYdpVSn4NOhqTFmFdYJWu95v/Ga3glMqGe7Eqxum8ErdzO4ndpzRynVKegVua1l32A9p+nwC0qpjk9Dv7UcWdDtbIjrFuhKlFKqSRr6raEXZSmlOhkN/dY4fgDKjurImkqpTkNDvzVqLsrSkTWVUp2Dhn5r2DMhMgFSBwa6EqWU8omGfms4MiFtFITZAl2JUkr5REO/pSpOwOEd2rSjlOpUNPRbKmcTGLf23FFKdSoa+i3lqL5T1qjA1qGUUs2god9S9kzrBG5McqArUUopn2not4QxVndNHXpBKdXJaOi3xLF9cLJAR9ZUSnU6GvotUTPImoa+Uqpz0dBvCXsmRCdByjmBrkQppZpFQ78lqtvzw/TnU0p1LppazVVeBEd2adOOUqpT0tBvLkc2YPQkrlKqU9LQby5HFiDQRy/KUkp1Phr6zWXPhO6DITox0JUopVSzaeg3h9ttNe9o045SqpPS0G+Oo3ugokhDXynVaWnoN4e9epA1DX2lVOekod8cjkyI6Qrdzgp0JUop1SIa+s1hz7SadkQCXYlSSrWIT6EvIlNFZI+I7BORhfUs7ysia0Rks4hsE5HLvZYNF5H1IrJDRLaLSLQ/v0C7KTsOR/fqyJpKqU4tvKkVRMQGPAdMARxAloh8YIzZ6bXaImClMeZ5ERkMrAL6iUg48DpwozFmq4h0A6r8/i3aQ85G61lP4iqlOjFfjvTHAvuMMQeMMZXACmB6nXUMUN1xPQnI9UxfAmwzxmwFMMYcM8a4Wl92ANg3gNig93mBrkQppVrMl9DvA9i9Xjs887w9CNwgIg6so/x5nvnnAEZEPhaRTSLyq1bWGzj2TOgxBKLiA12JUkq1mL9O5M4GXjHGpAGXA6+JSBhW89GFwPWe5ytFZHLdjUXkdhHJFpHs/Px8P5XkR26X1byjTTtKqU7Ol9DPAdK9Xqd55nm7BVgJYIxZD0QDKVh/Faw1xhw1xpRh/RVwWvuIMWaxMWa0MWZ0ampq879FWzuyEypLIH1coCtRSqlW8SX0s4ABItJfRCKBWcAHddY5BEwGEJFBWKGfD3wMDBORWM9J3YnATjqbmouytOeOUqpza7L3jjHGKSJzsQLcBiwxxuwQkYeBbGPMB8A9wIsisgDrpO4cY4wBCkTkKawdhwFWGWP+0VZfps04siAuFZL7BboSpZRqlSZDH8AYswqracZ73m+8pncCExrY9nWsbpudlz3TGnpBL8pSSnVyekVuU0qPwvH9ehJXKRUUNPSb4siynjX0lVJBQEO/KfZMCAuH3iMDXYlSSrWahn5T7JnQczhExAS6EqWUajUN/ca4nJC7SZt2lFJBQ0O/MYe/gaoy7Z+vlAoaGvqN0ZO4Sqkgo6HfGPsGSOgFSelNr6uUUp2Ahn5j7JlW045elKWUChIa+g05cRgKv9emHaVUUNHQb4jDM8iajqyplAoiGvoNsWeCLRJ6ZQS6EqWU8hsN/YY4sqzAD48KdCVKKeU3Gvr1cVZCziZt2lFKBR0N/fr8sB1cFXpRllIq6Gjo16fmJK723FFKBRcN/frYMyExDRJ7B7oSpZTyKw39+tgz9ShfKRWUNPTrKs6FYoeGvlIqKGno12X3tOenaegrpYKPhn5d9kwIj4aewwJdiVJK+Z2Gfl2OTOvWiOGRga5EKaX8TkPfm7MC8rZq/3ylVNDS0PeWtxVclXolrlIqaGnoe7NvsJ61545SKkhp6HuzZ0KXMyC+e6ArUUqpNuFT6IvIVBHZIyL7RGRhPcv7isgaEdksIttE5PJ6lpeIyL3+KtzvjLFG1tSjfKVUEGsy9EXEBjwHXAYMBmaLyOA6qy0CVhpjRgKzgD/VWf4U8GHry21DRXY4kaft+UqpoObLkf5YYJ8x5oAxphJYAUyvs44BEj3TSUBu9QIR+RnwHbCj9eW2oZqLsrTnjlIqePkS+n0Au9drh2eetweBG0TEAawC5gGISDzw38BDra60rTmyICIWegwNdCVKKdVm/HUidzbwijEmDbgceE1EwrB2Bv9njClpbGMRuV1EskUkOz8/308lNZM9E/qMAlt4YD5fKaXagS+hnwOke71O88zzdguwEsAYsx6IBlKAccDjInIQuBv4tYjMrfsBxpjFxpjRxpjRqampzf4SrVZ1En7Ypk07Sqmg58thbRYwQET6Y4X9LOC6OuscAiYDr4jIIKzQzzfGXFS9gog8CJQYY571R+F+lbsZ3E7tuaOUCnpNHukbY5zAXOBjYBdWL50dIvKwiEzzrHYPcJuIbAWWA3OMMaativY7HVlTKRUifGrANsaswjpB6z3vN17TO4EJTbzHgy2or33YM6HrWRDXLdCVKKVUm9Irco2xRtbUph2lVAjQ0C84CKX5ehJXKRUSNPSr2/P1SlylVAjQ0HdkQmQCdB8U6EqUUqrNaejbM6HPeRBmC3QlSinV5kI79CtK4PAObdpRSoWM0A793E1gXNpzRykVMkI79Gsuyhod2DqUUqqdhHboO7Ig5VyISQ50JUop1S5CN/SNsY7007V/vlIqdIRu6B/bDyeP63g7SqmQErqh76i+KEtDXykVOkI39O2ZEJVktekrpVSICO3QTxsNYaH7EyilQk9oJl55MRzZqU07SqmQE5qhn7MRMBr6SqmQE5qhb88EBProRVlKqdASmqHvyLRG1YxODHQlSinVrkIv9N1u60pcbdpRSoWg0Av9o3uhvEgvylJKhaTQC329KEspFcJCL/TtmdYAa93ODnQlSinV7kIv9B1ZVtOOSKArUUqpdhdaoX+yAPJ368iaSqmQFVqh79hoPetJXKVUiAqx0M8ECYM+owJdiVJKBYRPoS8iU0Vkj4jsE5GF9SzvKyJrRGSziGwTkcs986eIyEYR2e55nuTvL9As9g3QYwhExQe0DKWUCpQmQ19EbMBzwGXAYGC2iAyus9oiYKUxZiQwC/iTZ/5R4KfGmGHAz4HX/FV4s7ldVvOONu0opUKYL0f6Y4F9xpgDxphKYAUwvc46Bqge0yAJyAUwxmw2xuR65u8AYkQkqvVlt0D+bqg8of3zlVIhLdyHdfoAdq/XDmBcnXUeBD4RkXlAHPDjet7namCTMaaiBXW2nl0vylJKKX+dyJ0NvGKMSQMuB14TkZr3FpEhwO+Af69vYxG5XUSyRSQ7Pz/fTyXVYc+E2BRI7t8276+UUp2AL6GfA6R7vU7zzPN2C7ASwBizHogGUgBEJA14F7jJGLO/vg8wxiw2xow2xoxOTU1t3jfwlSPTOsrXi7KUUiHMl9DPAgaISH8RicQ6UftBnXUOAZMBRGQQVujni0gX4B/AQmPMv/xXdjOVHoNj+7RpRykV8poMfWOME5gLfAzswuqls0NEHhaRaZ7V7gFuE5GtwHJgjjHGeLY7G/iNiGzxPLq3yTdpjCPLetaeO0qpEOfLiVyMMauAVXXm/cZreicwoZ7tHgUebWWNrefIhLBw6D0y0JWoIFBVVYXD4aC8vDzQpagQFB0dTVpaGhERES3a3qfQ7/TsmdBzGETGBroSFQQcDgcJCQn069cP0XNEqh0ZYzh27BgOh4P+/VvWKSX4h2FwOa0boWvTjvKT8vJyunXrpoGv2p2I0K1bt1b9lRn8oX9kB1SV6Ulc5Vca+CpQWvtvL/hDv/qirDQdTlkppYI/9B1ZEN8TuvQNdCVKBUR8vP8GGHzvvffYuXOn396vMRdccEGLtnvwwQd58skn/VxN8Aj+0LdvsG6aon+OK9Vq7RH6TqcTgHXr1rXp57S16u/R0QR3752SfCg4CKNvCXQlKkg99Lcd7Mwt9ut7Du6dyAM/HdLg8oULF5Kens6dd94JWEekc2NCAAAROklEQVS24eHhrFmzhoKCAqqqqnj00UeZPr3uuIj1+93vfsfrr79OWFgYl112GY899hgvvvgiixcvprKykrPPPpvXXnuNLVu28MEHH/Dll1/y6KOP8vbbbwNw5513kp+fT2xsLC+++CIDBw5k//79XH/99ZSWljJ9+nSefvppSkpKMMbwq1/9ig8//BARYdGiRcycOZMvvviC//mf/yE5OZndu3ezd+9e4uPjKSkpaVaNsbFN99BraLvDhw9zxx13cODAAQCef/55LrjgApYuXcqTTz6JiDB8+HBee+015syZw09+8hOuueYagJpa6/seP/vZz7Db7ZSXlzN//nxuv/12AD766CN+/etf43K5SElJ4dNPP+Xcc89l3bp1pKam4na7Oeecc1i/fj3+HKkguEPfUT3IWt3x4ZTqvGbOnMndd99dE/orV67k448/5q677iIxMZGjR48yfvx4pk2b1uRJvw8//JD333+fDRs2EBsby/HjxwG46qqruO222wBYtGgRL730EvPmzWPatGm1wm7y5Mm88MILDBgwgA0bNvDLX/6S1atXM3/+fObPn8/s2bN54YUXaj7vnXfeYcuWLWzdupWjR48yZswYLr74YgA2bdrEN998c1pXxObW2JSGtrvrrruYOHEi7777Li6Xi5KSEnbs2MGjjz7KunXrSElJqfnsxtT9HkuWLKFr166cPHmSMWPGcPXVV+N2u7nttttYu3Yt/fv35/jx44SFhXHDDTewbNky7r77bj777DMyMjL8GvgQ7KFv3wBhEdArI9CVqCDV2BF5Wxk5ciRHjhwhNzeX/Px8kpOT6dmzJwsWLGDt2rWEhYWRk5PD4cOH6dmzZ6Pv9dlnn3HzzTfXHCF37doVgG+++YZFixZRWFhISUkJl1566WnblpSUsG7dOmbMmFEzr6LCGkR3/fr1vPfeewBcd9113HvvvQD885//ZPbs2dhsNnr06MHEiRPJysoiMTGRsWPH1tv3vDU11qeh7VavXs3SpUsBsNlsJCUlsXTpUmbMmEFKSkqtz25M3e/xzDPP8O677wJgt9v59ttvyc/P5+KLL65Zr/p9f/GLXzB9+nTuvvtulixZws033+zTd2qOIA/9LCvwI6IDXYlSfjVjxgzeeustfvjhB2bOnMmyZcvIz89n48aNRERE0K9fv1b15Z4zZw7vvfceGRkZvPLKK3zxxRenreN2u+nSpQtbtmxpxTc5JS4uzu81+nM7b+Hh4bjdbsD6HSorK2uWeX+PL774gs8++4z169cTGxvLj370o0b/u6Snp9OjRw9Wr15NZmYmy5Yta3ZtTQneE7muKsjdrE07KijNnDmTFStW8NZbbzFjxgyKioro3r07ERERrFmzhu+//96n95kyZQovv/wyZWVlADXNFydOnKBXr15UVVXVCp6EhAROnDgBQGJiIv379+fNN98ErKtFt27dCsD48eNr2vxXrFhRs/1FF13EG2+8gcvlIj8/n7Vr1zJ2bOPX0DS3xqY0tN3kyZN5/vnnAXC5XBQVFTFp0iTefPNNjh07Vuuz+/Xrx8aNGwH44IMPqKqqqvezioqKSE5OJjY2lt27d/P111/X/D5r167lu+++q/W+ALfeeis33HADM2bMwGaz+fy9fBW8of/DdnCetHruKBVkhgwZwokTJ+jTpw+9evXi+uuvJzs7m2HDhrF06VIGDhzo0/tMnTqVadOmMXr0aEaMGFHT1fGRRx5h3LhxTJgwodZ7zZo1iyeeeIKRI0eyf/9+li1bxksvvURGRgZDhgzh/fffB+Dpp5/mqaeeYvjw4ezbt4+kpCQArrzySoYPH05GRgaTJk3i8ccfb7IJqrk1NqWh7f7whz+wZs0ahg0bxqhRo9i5cydDhgzh/vvvZ+LEiWRkZPCf//mfANx22218+eWXZGRksH79+gb/Spk6dSpOp5NBgwaxcOFCxo8fD0BqaiqLFy/mqquuIiMjg5kzZ9ZsM23aNEpKStqkaQdArMEwO47Ro0eb7Ozs1r/R1y/AR/8NC3ZCUp/Wv59SHrt27WLQoEGBLqNDKysrIyYmBhFhxYoVLF++vGaHoBqXnZ3NggUL+Oqrrxpcp75/gyKy0Rgzuqn3D942fUcmJPbRwFcqADZu3MjcuXMxxtClSxeWLFkS6JI6hccee4znn3++TdryqwVv6NuzdLwdpTy2b9/OjTfeWGteVFQUGzZsaJPPu+iii2ra9wPlzjvv5F//qn3vpvnz57dZs4k/LFy4kIULF7bpZwRn6BfnQdEhGP8fga5EqQ5h2LBhfutl01k899xzgS6hQwrOE7k1F2Xpkb5SSnkLztC3Z4ItCnoOD3QlSinVoQRv6PceCeGRga5EKaU6lOALfWcF5G3R/vlKKVWP4Av9vG3gqtQrcZXCGgagvYYovvzyyyksLGz2dq+88gpz585tg4pUfYKv9071SVy9J65SfPHFF8THx7f4hiS+MMZgjGHVqlVt9hntofp7hIUF37Gwt+ALffsG6y5ZCT0CXYkKBR8utIb88Keew+Cyxxpdpe4Y79deey2PPvoolZWVdOvWjWXLlnHy5EleeOEFbDYbr7/+On/84x8ZOHAgd9xxB4cOHQKs4RImTJhAfn4+1113Hbm5uZx//vl8+umnbNy4kZSUFJ566qmai6tuvfVW7r77bg4ePMill17KuHHj2LhxI6tWrWLixIlkZ2eTkpJS7xj0f/vb306rsUePpv8/bWi7kpIS5s2bR3Z2NiLCAw88wNVXX33aOPWff/45Dz74IPHx8TWjfQ4dOpS///3vAKd9j8cee4ysrCxOnjzJNddcw0MPPQRAVlYW8+fPp7S0lKioKD7//HOuuOIKnnnmGUaMGAHAhRdeyHPPPUdGRscd2TcIQz8Lzmi7oxqlAq2+Md5FhK+//hoR4S9/+QuPP/44v//977njjjtqhd11113HggULuPDCCzl06BCXXnopu3bt4qGHHmLSpEncd999fPTRR7z00kuAdWXtyy+/zIYNGzDGMG7cOCZOnEhycjLffvstr776as14Mo3VB1Yg1ldjUxra7pFHHiEpKYnt262dbkFBAfn5+aeNU9+Uut/jt7/9LV27dsXlcjF58mS2bdvGwIEDmTlzJm+88QZjxoyhuLiYmJgYbrnlFl555RWefvpp9u7dS3l5eYcOfAi20C9ywIlcbc9X7aeJI/K2sHr16tPGeN++fTszZ84kLy+PysrKeselB2tseu/bHRYXF1NSUsI///nPmjHfp06dSnJyMmCNf3/llVfWDCh21VVX8dVXXzFt2jTOOOOM0wK/ofoAHA6HTzXW1dB2n332Wa0RPJOTk/nb3/5W7zj1jan7PVauXMnixYtxOp3k5eWxc+dORIRevXoxZozVQSQxMRGwhrh+5JFHeOKJJ1iyZAlz5szx6TsFUnA1Xtk9l5Rrzx0VYubNm8fcuXPZvn07f/7znxscs93tdvP111+zZcsWtmzZQk5OTotvnN7c8e99rdFf23nzHv8eqPUe3t/ju+++48knn+Tzzz9n27ZtXHHFFY1+XmxsLFOmTOH9999n5cqVXH/99c2urb35FPoiMlVE9ojIPhE5bWAIEekrImtEZLOIbBORy72W3efZbo+I+HZrm5ayZ0F4DPQY2qYfo1Qg1TfGe1FREX36WIMLvvrqqzXreo9/D3DJJZfwxz/+seZ19dAMEyZMYOXKlQB88sknFBQUANYYOu+99x5lZWWUlpby7rvvctFFFzW7PqDBGpvS0HZTpkypNdRCQUFBg+PU9+vXj02bNgHW7Qyrl9dVXFxMXFwcSUlJHD58mA8//BCAc889l7y8PLKysgBrTP7qG5/feuut3HXXXYwZM6bmL6SOrMnQFxEb8BxwGTAYmC0ig+ustghYaYwZCcwC/uTZdrDn9RBgKvAnz/u1DUcm9BkFtog2+wilAq2+Md4ffPBBZsyYwahRo2qaVQB++tOf8u677zJixAi++uornnnmGbKzsxk+fDiDBw+uuX/tAw88wCeffMLQoUN588036dmzJwkJCZx33nnMmTOHsWPHMm7cOG699VZGjhzZ7PqABmtsSkPbLVq0iIKCAoYOHUpGRgZr1qxpcJz6q6++muPHjzNkyBCeffZZzjnnnHo/KyMjg5EjRzJw4ECuu+46JkyYAEBkZCRvvPEG8+bNIyMjgylTptT8BTBq1CgSExM79EButVR3U2roAZwPfOz1+j7gvjrr/Bn4b6/119W3LvAxcH5jnzdq1CjTIpVlxjzUzZhPH2jZ9kr5aOfOnYEuwe/Ky8tNVVWVMcaYdevWmYyMjABX1Hnk5OSYAQMGGJfL1W6fWd+/QSDbNJHnxhifTuT2Aexerx1A3TOlDwKfiMg8IA74sde2X9fZtm0GuC8vhsHT4cx/a5O3VyqYHTp0iGuvvRa3201kZCQvvvhioEvqFJYuXcr999/PU0891Wn69/ur985s4BVjzO9F5HzgNRHxuWFdRG4Hbgfo27dvyypI6AHXvNSybZUKcQMGDGDz5s0BreG3v/1tzf12q82YMYP7778/QBU17aabbuKmm24KdBnN4kvo5wDpXq/TPPO83YLVZo8xZr2IRAMpPm6LMWYxsBis2yX6WrxSKnjcf//9HTrgg4Uvf49kAQNEpL+IRGKdmP2gzjqHgMkAIjIIiAbyPevNEpEoEekPDAAy/VW8UoFiOti9pVXoaO2/vSaP9I0xThGZi3US1gYsMcbsEJGHsU4cfADcA7woIgsAA8zxnFjYISIrgZ2AE7jTGONqVcVKBVh0dDTHjh2jW7duiEigy1EhxBjDsWPHiI6ObvF7SEc7Yhk9erTJzs4OdBlKNaiqqgqHw9Gii4SUaq3o6GjS0tKIiKjdNV1ENhpjRje1fXANw6BUO4iIiPB5CAGlOprO0cdIKaWUX2joK6VUCNHQV0qpENLhTuSKSD7wfSveIgU46qdyOjv9LWrT36M2/T1OCYbf4gxjTGpTK3W40G8tEcn25Qx2KNDfojb9PWrT3+OUUPottHlHKaVCiIa+UkqFkGAM/cWBLqAD0d+iNv09atPf45SQ+S2Crk1fKaVUw4LxSF8ppVQDgib0m7qPbygRkXTPPYt3isgOEZkf6JoCTURsnns4/z3QtQSaiHQRkbdEZLeI7PLcAyNkicgCz/8n34jIcs/Q8EErKELfx/v4hhIncI8xZjAwHrgzxH8PgPnArkAX0UH8AfjIGDMQyCCEfxcR6QPcBYw2xgzFGkl4VmCraltBEfrAWGCfMeaAMaYSWAFMD3BNAWOMyTPGbPJMn8D6n7ptblPZCYhIGnAF8JdA1xJoIpIEXAy8BGCMqTTGFAa2qoALB2JEJByIBXIDXE+bCpbQr+8+viEbct5EpB8wEtgQ2EoC6mngV4A70IV0AP2xbnD0sqe56y8iEhfoogLFGJMDPIl1I6g8oMgY80lgq2pbwRL6qh4iEg+8DdxtjCkOdD2BICI/AY4YYzYGupYOIhw4D3jeGDMSKAVC9hyYiCRjtQr0B3oDcSJyQ2CralvBEvo+3Ys3lIhIBFbgLzPGvBPoegJoAjBNRA5iNftNEpHXA1tSQDkAhzGm+i+/t7B2AqHqx8B3xph8Y0wV8A5wQYBralPBEvq+3Mc3ZIh1D7+XgF3GmKcCXU8gGWPuM8akGWP6Yf27WG2MCeojucYYY34A7CJyrmfWZKzbmYaqQ8B4EYn1/H8zmSA/sR0Ud85q6D6+AS4rkCYANwLbRWSLZ96vjTGrAliT6jjmAcs8B0gHgJsDXE/AGGM2iMhbwCasXm+bCfKrc/WKXKWUCiHB0ryjlFLKBxr6SikVQjT0lVIqhGjoK6VUCNHQV0qpEKKhr5RSIURDXymlQoiGvlJKhZD/Dz47gw+NC+uWAAAAAElFTkSuQmCC\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_4 (Dense)              (None, 100)               78500     \n",
      "_________________________________________________________________\n",
      "Crashes (Lambda)             (None, 100)               0         \n",
      "_________________________________________________________________\n",
      "dense_5 (Dense)              (None, 100)               10100     \n",
      "_________________________________________________________________\n",
      "dense_6 (Dense)              (None, 10)                1010      \n",
      "=================================================================\n",
      "Total params: 89,610\n",
      "Trainable params: 89,610\n",
      "Non-trainable params: 0\n",
      "_________________________________________________________________\n",
      "_________________________________________________________________\n",
      "Layer (type)                 Output Shape              Param #   \n",
      "=================================================================\n",
      "dense_7 (Dense)              (None, 100)               78500     \n",
      "_________________________________________________________________\n",
      "dense_8 (Dense)              (None, 100)               10100     \n",
      "_________________________________________________________________\n",
      "dense_9 (Dense)              (None, 10)                1010      \n",
      "=================================================================\n",
      "Total params: 89,610\n",
      "Trainable params: 89,610\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 = 'l2', reg_coeff = 1e-3, do_print = True,\n",
    "                                      name = 'experiment_weights', seed = None, batch_size = 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 478,
   "metadata": {},
   "outputs": [],
   "source": [
    "truth = exp.y_train"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 479,
   "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": 470,
   "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": 396,
   "metadata": {},
   "outputs": [],
   "source": [
    "correct = correct_loss()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 397,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 20/20 [00:19<00:00,  1.04it/s]\n"
     ]
    }
   ],
   "source": [
    "crashing = [crashing_loss() for _ in tqdm(range(20))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 398,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.640429319700057"
      ]
     },
     "execution_count": 398,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "correct * 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 480,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([3., 1., 2., 2., 2., 4., 2., 1., 1., 2.]),\n",
       " array([0.64042332, 0.64043843, 0.64045354, 0.64046865, 0.64048376,\n",
       "        0.64049887, 0.64051398, 0.6405291 , 0.64054421, 0.64055932,\n",
       "        0.64057443]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 480,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAErFJREFUeJzt3XuQnXV9x/H3xyTgteKY1WKSda3SacWi6Bq1tjMUxzaAQ2zFDlovWG1mbKna8Qa2QxH/0U4r1uLgpEABq4JF24kaa6mXoragSQzBgDhRUYK0xIBgtGKj3/7x/JDD6W72JHvO2cP0/Zo5w3P5nef5sJuzn30u52yqCkmSHrDUASRJk8FCkCQBFoIkqbEQJEmAhSBJaiwESRJgIUiSGgtBkgRYCJKkZvlS7XjlypU1MzOzVLuXpPulrVu3freqpkax7SUrhJmZGbZs2bJUu5ek+6Uk3xrVtj1lJEkCLARJUmMhSJIAC0GS1FgIkiTgIAohybIkX07ysTnWHZ7k8iS7klyTZGaYISVJo3cwRwivBW6YZ90rgTuq6gnAucA7FhtMkjReAxVCktXAScAF8wxZD1zSpq8AnpMki48nSRqXQY8Q3gW8CfjpPOtXATcDVNV+4E7gkYtOJ0kamwXfqZzkecBtVbU1yXGL2VmSDcAGgOnp6cVsShqpmTM+viT7ventJy3JfiUY7Ajh2cDJSW4CLgOOT/L3fWNuAdYAJFkOPBzY27+hqtpYVbNVNTs1NZKP4pAkHaIFC6Gqzqyq1VU1A5wKfLqqXtI3bBPw8jZ9ShtTQ00qSRqpQ/5wuyTnAFuqahNwIfC+JLuA2+mKQ5J0P3JQhVBVnwU+26bP6ln+I+CFwwwmSRov36ksSQIsBElSYyFIkgALQZLUWAiSJMBCkCQ1FoIkCbAQJEmNhSBJAiwESVJjIUiSAAtBktRYCJIkwEKQJDUWgiQJsBAkSY2FIEkCBiiEJA9M8sUk1ybZmeStc4w5LcmeJNvb41WjiStJGpVB/oTm3cDxVbUvyQrg80k+UVVX9427vKpOH35ESdI4LFgIVVXAvja7oj1qlKEkSeM30DWEJMuSbAduA66sqmvmGPaCJDuSXJFkzVBTSpJGbqBCqKqfVNVTgNXA2iRP6hvyUWCmqo4BrgQumWs7STYk2ZJky549exaTW5I0ZAd1l1FVfQ/4DLCub/neqrq7zV4APG2e52+sqtmqmp2amjqUvJKkERnkLqOpJEe06QcBzwW+2jfmyJ7Zk4EbhhlSkjR6g9xldCRwSZJldAXyoar6WJJzgC1VtQl4TZKTgf3A7cBpowosSRqNQe4y2gEcO8fys3qmzwTOHG40SdI4+U5lSRJgIUiSGgtBkgRYCJKkxkKQJAEWgiSpsRAkSYCFIElqLARJEmAhSJIaC0GSBFgIkqTGQpAkARaCJKmxECRJgIUgSWosBEkSMNjfVH5gki8muTbJziRvnWPM4UkuT7IryTVJZkYRVpI0OoMcIdwNHF9VTwaeAqxL8sy+Ma8E7qiqJwDnAu8YbkxJ0qgtWAjV2ddmV7RH9Q1bD1zSpq8AnpMkQ0spSRq55YMMSrIM2Ao8AXhPVV3TN2QVcDNAVe1PcifwSOC7fdvZAGwAmJ6ePuTQM2d8/JCfu1g3vf2kJdu3JI3SQBeVq+onVfUUYDWwNsmTDmVnVbWxqmaranZqaupQNiFJGpGDusuoqr4HfAZY17fqFmANQJLlwMOBvcMIKEkaj0HuMppKckSbfhDwXOCrfcM2AS9v06cAn66q/usMkqQJNsg1hCOBS9p1hAcAH6qqjyU5B9hSVZuAC4H3JdkF3A6cOrLEkqSRWLAQqmoHcOwcy8/qmf4R8MLhRpMkjZPvVJYkARaCJKmxECRJgIUgSWosBEkSYCFIkhoLQZIEWAiSpMZCkCQBFoIkqbEQJEmAhSBJaiwESRJgIUiSGgtBkgRYCJKkxkKQJAGD/U3lNUk+k+T6JDuTvHaOMccluTPJ9vY4a65tSZIm1yB/U3k/8Pqq2pbkYcDWJFdW1fV94z5XVc8bfkRJ0jgseIRQVbdW1bY2/X3gBmDVqINJksbroK4hJJkBjgWumWP1s5Jcm+QTSY6e5/kbkmxJsmXPnj0HHVaSNDoDF0KShwIfBl5XVXf1rd4GPLaqngz8DfBPc22jqjZW1WxVzU5NTR1qZknSCAxUCElW0JXB+6vqI/3rq+quqtrXpjcDK5KsHGpSSdJIDXKXUYALgRuq6p3zjPn5No4ka9t29w4zqCRptAa5y+jZwEuB65Jsb8veAkwDVNV7gVOAVyfZD/w3cGpV1QjySpJGZMFCqKrPA1lgzHnAecMKJUkaP9+pLEkCLARJUmMhSJIAC0GS1FgIkiTAQpAkNRaCJAmwECRJjYUgSQIsBElSYyFIkgALQZLUWAiSJMBCkCQ1FoIkCbAQJEmNhSBJAgb7m8prknwmyfVJdiZ57RxjkuTdSXYl2ZHkqaOJK0kalUH+pvJ+4PVVtS3Jw4CtSa6squt7xpwAHNUezwDOb/+VJN1PLHiEUFW3VtW2Nv194AZgVd+w9cCl1bkaOCLJkUNPK0kamUGOEH4myQxwLHBN36pVwM0987vbslv7nr8B2AAwPT19cEn/n5s54+NLHUFjsFTf55veftKS7HcpLeVralK/3gNfVE7yUODDwOuq6q5D2VlVbayq2aqanZqaOpRNSJJGZKBCSLKCrgzeX1UfmWPILcCanvnVbZkk6X5ikLuMAlwI3FBV75xn2CbgZe1uo2cCd1bVrfOMlSRNoEGuITwbeClwXZLtbdlbgGmAqnovsBk4EdgF/BB4xfCjSpJGacFCqKrPA1lgTAF/NKxQkqTx853KkiTAQpAkNRaCJAmwECRJjYUgSQIsBElSYyFIkgALQZLUWAiSJMBCkCQ1FoIkCbAQJEmNhSBJAiwESVJjIUiSAAtBktQM8ic0L0pyW5KvzLP+uCR3JtneHmcNP6YkadQG+ROaFwPnAZceYMznqup5Q0kkSVoSCx4hVNVVwO1jyCJJWkLDuobwrCTXJvlEkqOHtE1J0hgNcspoIduAx1bVviQnAv8EHDXXwCQbgA0A09PTQ9i1JGlYFn2EUFV3VdW+Nr0ZWJFk5TxjN1bVbFXNTk1NLXbXkqQhWnQhJPn5JGnTa9s29y52u5Kk8VrwlFGSDwLHASuT7Ab+HFgBUFXvBU4BXp1kP/DfwKlVVSNLLEkaiQULoapetMD68+huS5Uk3Y/5TmVJEmAhSJIaC0GSBFgIkqTGQpAkARaCJKmxECRJgIUgSWosBEkSYCFIkhoLQZIEWAiSpMZCkCQBFoIkqbEQJEmAhSBJaiwESRIwQCEkuSjJbUm+Ms/6JHl3kl1JdiR56vBjSpJGbZAjhIuBdQdYfwJwVHtsAM5ffCxJ0rgtWAhVdRVw+wGGrAcurc7VwBFJjhxWQEnSeAzjGsIq4Oae+d1tmSTpfmT5OHeWZAPdaSWmp6fHueuhmTnj40sdQRo6/10LhnOEcAuwpmd+dVv2f1TVxqqararZqampIexakjQswyiETcDL2t1GzwTurKpbh7BdSdIYLXjKKMkHgeOAlUl2A38OrACoqvcCm4ETgV3AD4FXjCqsJGl0FiyEqnrRAusL+KOhJZIkLQnfqSxJAiwESVJjIUiSAAtBktRYCJIkwEKQJDUWgiQJsBAkSY2FIEkCLARJUmMhSJIAC0GS1FgIkiTAQpAkNRaCJAmwECRJjYUgSQIGLIQk65LcmGRXkjPmWH9akj1JtrfHq4YfVZI0SoP8TeVlwHuA5wK7gS8l2VRV1/cNvbyqTh9BRknSGAxyhLAW2FVV36iqHwOXAetHG0uSNG6DFMIq4Oae+d1tWb8XJNmR5Ioka4aSTpI0NsO6qPxRYKaqjgGuBC6Za1CSDUm2JNmyZ8+eIe1akjQMgxTCLUDvb/yr27Kfqaq9VXV3m70AeNpcG6qqjVU1W1WzU1NTh5JXkjQigxTCl4CjkjwuyWHAqcCm3gFJjuyZPRm4YXgRJUnjsOBdRlW1P8npwCeBZcBFVbUzyTnAlqraBLwmycnAfuB24LQRZpYkjcCChQBQVZuBzX3LzuqZPhM4c7jRJEnj5DuVJUmAhSBJaiwESRJgIUiSGgtBkgRYCJKkxkKQJAEWgiSpsRAkSYCFIElqLARJEmAhSJIaC0GSBFgIkqTGQpAkARaCJKmxECRJwICFkGRdkhuT7EpyxhzrD09yeVt/TZKZYQeVJI3WgoWQZBnwHuAE4InAi5I8sW/YK4E7quoJwLnAO4YdVJI0WoMcIawFdlXVN6rqx8BlwPq+MeuBS9r0FcBzkmR4MSVJozZIIawCbu6Z392WzTmmqvYDdwKPHEZASdJ4LB/nzpJsADa02X1JblzgKSuB74421aKYb/EmPaP5FmfS88ESZMzBnVTvz/fYoYbpMUgh3AKs6Zlf3ZbNNWZ3kuXAw4G9/Ruqqo3AxkHDJdlSVbODjh838y3epGc03+JMej6Y/IzjzDfIKaMvAUcleVySw4BTgU19YzYBL2/TpwCfrqoaXkxJ0qgteIRQVfuTnA58ElgGXFRVO5OcA2ypqk3AhcD7kuwCbqcrDUnS/chA1xCqajOwuW/ZWT3TPwJeONxowEGcXloi5lu8Sc9ovsWZ9Hww+RnHli+e2ZEkgR9dIUm6R1UN9QGsA24EdgFnzDPmd4HrgZ3AB/rW/Rzdex3O61n2NOC6ts13045seta/HihgZZv/PWBHe86/A0+epHw9y58O7AdOmbR8wHHA9raPf5vA7/HDgY8C17Z9vGIp8gFn091lt709Tux5zplt/I3Ab01SPuC5wNb2nK3A8Uv1PT7Q17Ctnwb2AW+YtHzAMcB/tH1cBzxwUvIBK+jeMHwdcANw5lw57rPfhQYczIPuovPXgV8ADqN7sT6xb8xRwJeBR7T5R/Wt/2vgA31fiC8CzwQCfAI4oWfdGroL3t/i3h8Wv9qz/ROAayYpX0+WT9NdmzllkvIBR7R/qNP9+5igjG8B3tGmp+huZjhs3PnoXoxv6H1+W/7Etu/Dgce1TMsmKN+xwGPa9JOAW5bqezxfxp7nXQH8wz1jJiUf3TXYHbRfOOnejDtJ3+MXA5e16QcDNwEz832dq2rop4wG+ZiLPwDeU1V3AFTVbfesSPI04NHAv/QsOxL4uaq6urr/s0uB5/ds71zgTXS/PdK2+e/3bB+4mu69ExOTr/lj4MPAbT3LJiXfi4GPVNW3+/cxQRkLeFj7iJSH0hXC/iXKN5f1dC/Gu6vqm3S/1a2dlHxV9eWq+k6b3Qk8KMnhbX4iMrbnPR/4Zst4j0nJ95vAjqq6tu1jb1X9ZILyFfCQ9t6wBwE/Bu460BOGXQiDfMzFLwK/mOQLSa5Osg4gyQOAvwLeMMc2d8+1zSTr6X6zufYAmV5J16YTky/JKuC3gfPn2NaS52v7eESSzybZmuRlE5jxPOCXge/QHRK/tqp+Ou58zelJdiS5KMkjep4zV45JydfrBcC2qrp7gexjzZjkocCbgbfOsa0lz9f2UUk+mWRbkjdNWL4rgB8AtwLfBv6yqm7nAMb60RU9+zyK7hz1auCqJL8CvATYXFW7B/lcvCQPpjtt8JsHGPMbdIXwaxOW713Am6vqp4fwGYDjyLec7nzlc+h+s/iPJFdX1dcmKONv0Z0vPR54PHBlks+NM19zPvA2ut/G3kb3Yv79QZ+81PmSHE336cTzvo6WMOPZwLlVtW+pXicL5FtO97Pl6cAPgU8l2TpB+dYCPwEeAzwC+FySf62qbxwo1DAN8jEXu+nO6f8P8M0kX6P7wjwL+PUkf0h3CuCwJPvozqWtnmObj6c7N3tt+8KtBrYlWVtV/5nkGOACuvNs93yMxkTkA2aBy9rylcCJSfZPUL7dwN6q+gHwgyRXAU8GvjZBGV8BvL0dPu9K8k3gl8acj6r6r3sWJvlb4GNt9kA5JiEfSVYD/wi8rKq+3vP8SfkaPgM4Jclf0F3X+mmSH9FdBJ+EfLuBq6rqu23dZuCpwOcnJN+LgX9u+7gtyRfofvbMWwjDvqi8vO3scdx7MeXovjHrgEva9Eq6Q6tH9o05jQNfTDlxjn3fxL0XHKfpztn+6iTm61t+MfdeVJ6IfHSnYj7V8jwY+ArwpAnLeD5wdpt+NN2LY+W48wFH9oz5E+69iHc0972o/A26i42Tku+Itu/fWerX8XwZ+7Z1NvdeVJ6IfHS/dW+je40sB/4VOGmC8r0Z+Ls2/RC6G0WO6f/a3mefB1p5KA/gRLrfJL8O/Glbdg5wcpsO8M4W7jrg1Dm20f+FmKX7ofR1unPHmeM5N3HvD4sLgDu49zasLZOUr2/5xdz3ttOJyAe8se3jK8DrJvB7/Bi6i27Xtee9ZCnyAe9r29hB95levS/OP23jb+S+d00teT7gz+jOL2/veTxqkjL2bets7nvb6UTkozvFs7M97y8mKR/dEcY/tHzXA2/s30f/w3cqS5IA36ksSWosBEkSYCFIkhoLQZIEWAiSpMZCkCQBFoIkqbEQJEkA/C9/5KysQZI1vAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(np.array(crashing) * 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 481,
   "metadata": {},
   "outputs": [],
   "source": [
    "# computing probability delta of bad crash >= alpha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 482,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.1850873e-10 119.189575\n"
     ]
    }
   ],
   "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",
    "# balancedness\n",
    "q = q_factor(np.linalg.norm(exp.W[1], axis = 1, ord = 1) ** 2)\n",
    "\n",
    "# weight norm\n",
    "Wnorm = np.linalg.norm(exp.W[1].flatten(), ord = 1)\n",
    "\n",
    "# number of neurons\n",
    "n = N[1]\n",
    "\n",
    "# W matrix\n",
    "W = exp.model_correct.layers[1].weights[0]\n",
    "\n",
    "# loss function\n",
    "loss = exp.loss\n",
    "\n",
    "# probability of failure\n",
    "p = p_inf\n",
    "\n",
    "# get TF session\n",
    "sess = get_session()\n",
    "\n",
    "# q-factor\n",
    "print(q, Wnorm)\n",
    "\n",
    "# pseudo-weights\n",
    "W1 = exp.model_correct.layers[1].output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 483,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_feed():\n",
    "    \"\"\" Return feed dictionary for train dataset \"\"\"\n",
    "    return {exp.model_correct.input: exp.x_train, exp.output_tensor: exp.y_train}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 510,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.006751917012532552\n"
     ]
    }
   ],
   "source": [
    "# M1 ~ loss gradient\n",
    "M_1 = np.linalg.norm(sess.run(tf.multiply(W1, tf.gradients(exp.loss, W1)[0]), get_feed()), ord = 1) / len(exp.x_train)\n",
    "print(M_1)\n",
    "M_0 = 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 512,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "119.189575"
      ]
     },
     "execution_count": 512,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Wnorm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 511,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.5196012277701343e-05\n"
     ]
    }
   ],
   "source": [
    "# Maximal Hessian eigenvalue (handwavy)\n",
    "M_2 = 1\n",
    "\n",
    "# maximal alpha\n",
    "alpha_max = 1. / 2 * M_0 / M_2 / Wnorm ** 2\n",
    "\n",
    "# maximal alpha\n",
    "print(alpha_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 491,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.8324276688915037e-05\n"
     ]
    }
   ],
   "source": [
    "# maximal alpha\n",
    "alpha = alpha_max\n",
    "#alpha = max(p, alpha)\n",
    "if alpha > alpha_max:\n",
    "    alpha = alpha_max\n",
    "#alpha = min(alpha, alpha_max)\n",
    "if alpha < p:\n",
    "    alpha = p\n",
    "print(alpha)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 492,
   "metadata": {},
   "outputs": [],
   "source": [
    "# target probability to certify\n",
    "delta_0 = 1. / 8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 493,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.8324276688915037e-05 1e-05\n"
     ]
    }
   ],
   "source": [
    "# certifying low weight perturbation alpha\n",
    "print(alpha, p)\n",
    "assert alpha > p, \"Alpha must be greater than p\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 494,
   "metadata": {},
   "outputs": [],
   "source": [
    "# https://en.wikipedia.org/wiki/Binomial_distribution\n",
    "def kl_num(x, y):\n",
    "    \"\"\" KL divergence between numbers \"\"\"\n",
    "    return x * np.log(x / y) + (1 - x) * np.log((1 - x) / (1 - y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 501,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.116526371588494e-05"
      ]
     },
     "execution_count": 501,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kl_num(alpha, p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 498,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = 1\n",
    "delta = np.exp(-n * q * kl_num(alpha, p))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 509,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Certified 0.9988840967120627\n"
     ]
    }
   ],
   "source": [
    "print('Certified', delta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 497,
   "metadata": {},
   "outputs": [
    {
     "ename": "AssertionError",
     "evalue": "Delta must be small enough",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAssertionError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-497-87233110f836>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;32massert\u001b[0m \u001b[0mdelta\u001b[0m \u001b[0;34m<\u001b[0m \u001b[0;36m1.\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;36m4\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"Delta must be small enough\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mAssertionError\u001b[0m: Delta must be small enough"
     ]
    }
   ],
   "source": [
    "assert delta < 1. / 4, \"Delta must be small enough\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 433,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculating mean and variance using Taylor...\n",
    "err_mean_std = exp._get_bound_b3_loss(exp.x_train, exp.y_train, weights = exp.model_correct.layers[1].output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 434,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mean squared deviation\n",
    "mean_mean = np.mean(err_mean_std['mean'] ** 2)\n",
    "\n",
    "# mean variance\n",
    "mean_std = np.mean(err_mean_std['std'] ** 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 435,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1096.1976"
      ]
     },
     "execution_count": 435,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Wnorm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 436,
   "metadata": {},
   "outputs": [],
   "source": [
    "var_total = mean_mean + mean_std"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 437,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.114497952672306e-08\n"
     ]
    }
   ],
   "source": [
    "print(var_total)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 438,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-5.280081732462471e-08 0.00020284074530279778\n",
      "0.006404941863158373 4.257723123923872e-07\n"
     ]
    }
   ],
   "source": [
    "print(np.mean(err_mean_std['mean']), np.mean(err_mean_std['std'] ** 2) ** 0.5)\n",
    "print(np.mean(np.array(crashing)), np.std(crashing))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 449,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2 * variance =\n",
    "delta_L = (var_total) ** 0.5 * 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 450,
   "metadata": {},
   "outputs": [],
   "source": [
    "# total bad loss probability\n",
    "delta_total = var_total / delta_L ** 2 + delta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 454,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.28527163572124203\n"
     ]
    }
   ],
   "source": [
    "print(delta_total)\n",
    "assert delta_total < 1. / 3, \"Must have pfail < 1/3\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 455,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculating number of repetitions..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 456,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nout = exp.N[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 457,
   "metadata": {},
   "outputs": [],
   "source": [
    "# maximal probability of failure (<1/4)\n",
    "p_fail = 1. / 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 458,
   "metadata": {},
   "outputs": [],
   "source": [
    "# target probability\n",
    "delta_target = 1e-5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 459,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of repetitions 13.0 will achieve 1e-05\n"
     ]
    }
   ],
   "source": [
    "m = np.log(Nout / delta_target) / np.log(1 / p_fail)\n",
    "print('Number of repetitions', np.ceil(m), 'will achieve', delta_target)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 461,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "76.67171730593843"
      ]
     },
     "execution_count": 461,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alpha ** 2 * Wnorm ** 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 518,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'feed' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-518-7bb3db579f57>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0msess\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgradients\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mloss\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mexp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmodel_correct\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minput\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfeed\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'feed' is not defined"
     ]
    }
   ],
   "source": [
    "sess.run(tf.gradients(exp.loss, exp.model_correct.input), feed_dict = {exp.x_train[0]})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "sess = get_session()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 1000000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = tf.Variable(tf.random.uniform((N,)), tf.float64)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "loss = tf.reduce_sum(tf.exp(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def grad_and_mult(fcn, x):\n",
    "    return tf.reduce_sum(tf.multiply(tf.gradients(fcn, x)[0], tf.stop_gradient(x)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "loss = grad_and_mult(loss, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "sess.run(tf.global_variables_initializer())\n",
    "#sess.run([x, w])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "464207.6"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sess.run(loss)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "Wl = exp.model_correct.layers[-1].weights[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "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": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "feed = {exp.model_correct.input: x, exp.output_tensor: y}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "h = sess.run(tf.hessians(exp.loss, Wl), feed)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "h = np.abs(h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEACAYAAAB27puMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAD+lJREFUeJzt3H+onuV9x/H3Z8m02qHGH0iWhCWjoSMtjOpBMwpl1E2jlsU/pFhGDV3WMKptNwcz7h+hZTQdY9qw1hEaZwKlNriC2aoNwbaUwbQ5aUetZuLBanOCNtakcVugne67P87leMxOzjHPdU6enPh+wcO57+99Xff3fszBz7l/PE+qCkmSevzKqA9AkrTwGSaSpG6GiSSpm2EiSepmmEiSuhkmkqRuhokkqZthIknqZphIkrrNGiZJ7k9yOMmPBmoXJ9mb5Nn2c0mrJ8nWJBNJfpjkioE5G9r4Z5NsGKhfmeTJNmdrkgzbQ5I0Gpnt61SSfAD4T2BnVb231f4aOFJVW5JsBpZU1Z1JbgA+CdwAXA18oaquTnIxMA6MAQXsB66sqqNJvgd8CngCeATYWlWPnmqP2d7opZdeWitXrjzF/zyS9Pa2f//+n1XVZbONWzzbgKr6bpKVJ5TXA7/blncA3wHubPWdNZVQjye5KMnSNnZvVR0BSLIXWJfkO8AFVfV4q+8EbgIePdUeVfXiTO9j5cqVjI+Pz/Z2JUkDkrzwVsYNe8/k8oH/eb8EXN6WlwEHB8ZNttpM9clp6sP0kCSNSPcN+HaGMK9fPTxsjySbkownGX/55Zfn4cgkSTB8mPy0Xb6i/Tzc6oeAFQPjlrfaTPXl09SH6fH/VNW2qhqrqrHLLpv1kp8kaUjDhslu4I0nsjYADw/Ub21PXK0FjrVLVXuAa5MsaU9lXQvsadteTbK2PcV16wn7OpUekqQRmfUGfJKvMnUj/NIkk8DdwBZgV5KNwAvAh9vwR5h6ymoCOA58DKCqjiT5LLCvjfvMGzfjgU8ADwDnMXXj/dFWP6UekqTRmfXR4LPF2NhY+TSXJJ2aJPuramy2cX4CXpLUzTCRJHWb9Z7J2eLJQ8dYufkbQ819fsuNc3w0knR28cxEktTNMJEkdTNMJEndDBNJUjfDRJLUzTCRJHUzTCRJ3QwTSVI3w0SS1M0wkSR1M0wkSd0ME0lSN8NEktTNMJEkdTNMJEndDBNJUjfDRJLUzTCRJHUzTCRJ3QwTSVI3w0SS1M0wkSR1M0wkSd0ME0lSN8NEktTNMJEkdTNMJEndDBNJUjfDRJLUzTCRJHUzTCRJ3brCJMmfJXkqyY+SfDXJO5KsSvJEkokkX0tyTht7blufaNtXDuznrlZ/Jsl1A/V1rTaRZPNAfdoekqTRGDpMkiwDPgWMVdV7gUXALcDngXuq6l3AUWBjm7IRONrq97RxJFnT5r0HWAd8KcmiJIuALwLXA2uAj7SxzNBDkjQCvZe5FgPnJVkMnA+8CHwQeKht3wHc1JbXt3Xa9muSpNUfrKpfVNWPgQngqvaaqKrnquqXwIPA+jbnZD0kSSMwdJhU1SHgb4CfMBUix4D9wM+r6rU2bBJY1paXAQfb3Nfa+EsG6yfMOVn9khl6SJJGoOcy1xKmzipWAb8OvJOpy1RnjCSbkownGX/9+LFRH44knbV6LnP9HvDjqnq5qv4b+DrwfuCidtkLYDlwqC0fAlYAtO0XAq8M1k+Yc7L6KzP0eJOq2lZVY1U1tuj8CzveqiRpJj1h8hNgbZLz232Ma4CngW8DN7cxG4CH2/Lutk7b/q2qqla/pT3ttQpYDXwP2Aesbk9uncPUTfrdbc7JekiSRqDnnskTTN0E/z7wZNvXNuBO4I4kE0zd39jepmwHLmn1O4DNbT9PAbuYCqJvArdV1evtnsjtwB7gALCrjWWGHpKkEcjUH/pnv3OXrq6lG+4dau7zW26c46ORpIUhyf6qGpttnJ+AlyR1M0wkSd0ME0lSN8NEktTNMJEkdTNMJEndDBNJUjfDRJLUzTCRJHUzTCRJ3QwTSVI3w0SS1M0wkSR1M0wkSd0ME0lSN8NEktTNMJEkdTNMJEndDBNJUjfDRJLUzTCRJHUzTCRJ3QwTSVI3w0SS1M0wkSR1M0wkSd0ME0lSN8NEktTNMJEkdTNMJEndDBNJUjfDRJLUzTCRJHUzTCRJ3brCJMlFSR5K8u9JDiT5nSQXJ9mb5Nn2c0kbmyRbk0wk+WGSKwb2s6GNfzbJhoH6lUmebHO2JkmrT9tDkjQavWcmXwC+WVW/Bfw2cADYDDxWVauBx9o6wPXA6vbaBNwHU8EA3A1cDVwF3D0QDvcBHx+Yt67VT9ZDkjQCQ4dJkguBDwDbAarql1X1c2A9sKMN2wHc1JbXAztryuPARUmWAtcBe6vqSFUdBfYC69q2C6rq8aoqYOcJ+5quhyRpBHrOTFYBLwP/kOQHSb6c5J3A5VX1YhvzEnB5W14GHByYP9lqM9Unp6kzQ483SbIpyXiS8dePHxvmPUqS3oKeMFkMXAHcV1XvA/6LEy43tTOK6ugxq5l6VNW2qhqrqrFF5184n4chSW9rPWEyCUxW1RNt/SGmwuWn7RIV7efhtv0QsGJg/vJWm6m+fJo6M/SQJI3A0GFSVS8BB5O8u5WuAZ4GdgNvPJG1AXi4Le8Gbm1Pda0FjrVLVXuAa5MsaTferwX2tG2vJlnbnuK69YR9TddDkjQCizvnfxL4SpJzgOeAjzEVULuSbAReAD7cxj4C3ABMAMfbWKrqSJLPAvvauM9U1ZG2/AngAeA84NH2Athykh6SpBHI1C2Hs9+5S1fX0g33DjX3+S03zvHRSNLCkGR/VY3NNs5PwEuSuhkmkqRuhokkqZthIknqZphIkroZJpKkboaJJKmbYSJJ6maYSJK6GSaSpG6GiSSpm2EiSepmmEiSuhkmkqRuhokkqZthIknqZphIkroZJpKkboaJJKmbYSJJ6maYSJK6GSaSpG6GiSSpm2EiSepmmEiSuhkmkqRuhokkqZthIknqZphIkroZJpKkboaJJKmbYSJJ6maYSJK6dYdJkkVJfpDkn9v6qiRPJJlI8rUk57T6uW19om1fObCPu1r9mSTXDdTXtdpEks0D9Wl7SJJGYy7OTD4NHBhY/zxwT1W9CzgKbGz1jcDRVr+njSPJGuAW4D3AOuBLLaAWAV8ErgfWAB9pY2fqIUkaga4wSbIcuBH4clsP8EHgoTZkB3BTW17f1mnbr2nj1wMPVtUvqurHwARwVXtNVNVzVfVL4EFg/Sw9JEkj0Htmci/wF8D/tPVLgJ9X1WttfRJY1paXAQcB2vZjbfz/1U+Yc7L6TD0kSSMwdJgk+RBwuKr2z+HxzKkkm5KMJxl//fixUR+OJJ21FnfMfT/wB0luAN4BXAB8AbgoyeJ25rAcONTGHwJWAJNJFgMXAq8M1N8wOGe6+isz9HiTqtoGbAM4d+nq6nivkqQZDH1mUlV3VdXyqlrJ1A30b1XVHwLfBm5uwzYAD7fl3W2dtv1bVVWtfkt72msVsBr4HrAPWN2e3Dqn9djd5pyshyRpBObjcyZ3AnckmWDq/sb2Vt8OXNLqdwCbAarqKWAX8DTwTeC2qnq9nXXcDuxh6mmxXW3sTD0kSSOQqT/0z37nLl1dSzfcO9Tc57fcOMdHI0kLQ5L9VTU22zg/AS9J6maYSJK6GSaSpG6GiSSpm2EiSepmmEiSuhkmkqRuhokkqZthIknqZphIkroZJpKkboaJJKmbYSJJ6maYSJK6GSaSpG6GiSSpm2EiSepmmEiSuhkmkqRuhokkqZthIknqZphIkroZJpKkboaJJKmbYSJJ6maYSJK6GSaSpG6GiSSpm2EiSepmmEiSuhkmkqRuhokkqZthIknqZphIkroNHSZJViT5dpKnkzyV5NOtfnGSvUmebT+XtHqSbE0ykeSHSa4Y2NeGNv7ZJBsG6lcmebLN2ZokM/WQJI1Gz5nJa8CfV9UaYC1wW5I1wGbgsapaDTzW1gGuB1a31ybgPpgKBuBu4GrgKuDugXC4D/j4wLx1rX6yHpKkERg6TKrqxar6flv+D+AAsAxYD+xow3YAN7Xl9cDOmvI4cFGSpcB1wN6qOlJVR4G9wLq27YKqeryqCth5wr6m6yFJGoE5uWeSZCXwPuAJ4PKqerFtegm4vC0vAw4OTJtstZnqk9PUmaGHJGkEusMkya8B/wj8aVW9OritnVFUb4+ZzNQjyaYk40nGXz9+bD4PQ5Le1rrCJMmvMhUkX6mqr7fyT9slKtrPw61+CFgxMH15q81UXz5NfaYeb1JV26pqrKrGFp1/4XBvUpI0q56nuQJsBw5U1d8ObNoNvPFE1gbg4YH6re2prrXAsXapag9wbZIl7cb7tcCetu3VJGtbr1tP2Nd0PSRJI7C4Y+77gY8CTyb5t1b7S2ALsCvJRuAF4MNt2yPADcAEcBz4GEBVHUnyWWBfG/eZqjrSlj8BPACcBzzaXszQQ5I0AkOHSVX9C5CTbL5mmvEF3HaSfd0P3D9NfRx47zT1V6brIUkaDT8BL0nqZphIkroZJpKkboaJJKmbYSJJ6maYSJK6GSaSpG6GiSSpm2EiSepmmEiSuhkmkqRuhokkqZthIknqZphIkroZJpKkboaJJKmbYSJJ6maYSJK6GSaSpG6GiSSpm2EiSepmmEiSuhkmkqRuhokkqZthIknqZphIkroZJpKkboaJJKmbYSJJ6maYSJK6GSaSpG6GiSSpm2EiSepmmEiSui3YMEmyLskzSSaSbB718UjS29mCDJMki4AvAtcDa4CPJFkz2qOSpLevxaM+gCFdBUxU1XMASR4E1gNPz0ezlZu/MR+71Uk8v+XGUR+CpFO0UMNkGXBwYH0SuHpEx6I5Znir13z8QTIXv5e9xzXsMZyOP9BSVfPeZK4luRlYV1V/3NY/ClxdVbefMG4TsKmtvht4ZsiWFwLHhpx7JjtT39cojut09JzrHnO1v579DDv3UuBnQ/bUqev5N/6NqrpstkEL9czkELBiYH15q71JVW0DtvU2S7KtqjbNPnJhOVPf1yiO63T0nOsec7W/nv0MOzfJeFWNDdNTp+50/H4vyBvwwD5gdZJVSc4BbgF2z2O/f5rHfY/Smfq+RnFcp6PnXPeYq/317OdM/R3Sm837v9OCvMwFkOQG4F5gEXB/Vf3ViA9J0lvkmcnZZ8GGiaSFK8mmdhlaZwnDRJLUbaHeM5EknUEME0lSN8NEktTNMJF0Rknym0m2J3lo1Meit84wkTRnktyf5HCSH51Qf8vf8l1Vz1XVxvk9Us21hfoJeElnpgeAvwN2vlEY+Jbv32fqe/T2JdnN1GfEPnfC/D+qqsOn51A1lwwTSXOmqr6bZOUJ5Wm/5buqPgd86PQeoeaLl7kkzbfpvuV72ckGJ7kkyd8D70ty13wfnOaGZyaSzihV9QrwJ6M+Dp0az0wkzbe39C3fWtgME0nz7XR/y7dGwDCRNGeSfBX4V+DdSSaTbKyq14DbgT3AAWBXVT01yuPU3POLHiVJ3TwzkSR1M0wkSd0ME0lSN8NEktTNMJEkdTNMJEndDBNJUjfDRJLUzTCRJHX7X3Ba1hfsKxRBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(np.abs(h.flatten()))\n",
    "plt.xscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0031329184"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(np.abs(h))"
   ]
  },
  {
   "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
}