notebooks/certification_IBP_training.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f9bf98ec",
   "metadata": {},
   "source": [
    "# Certification and Training using Interval Bound Propagation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a04e2990",
   "metadata": {},
   "source": [
    "In this notebook we will look at using interval bound propagation (IBP) to assess a neural network's certified robustness. \n",
    "\n",
    "To do so we will use a datapoint's interval, also called box, representation. This captures the minimum and maximum values that a feature can take. \n",
    "\n",
    "Then we propagate the interval through the neural network and, by using interval arithmetic on the neural network components we can determine if a datapoint could have its class changed. \n",
    "\n",
    "The interval domain has the great advantage that it is *fast*. However, this speed comes at the expense of precision - we aggressively over-approximate in the forward pass and so we can often only certify a small subset of the data which is safe. More formally, this technique is sound but incomplete.\n",
    "\n",
    "We can see an example of how imprecision arises by looking at a neural network layer that causes a rotation: "
   ]
  },
  {
   "attachments": {
    "box_domain.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxoAAADjCAYAAAACNWIeAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdd3gU5cL+8e9mk5BAAiGVXgSMICBdEQ8CCkEgSokCIoSA0o6A+Io/FbGAEIqAKP3QqwRIQgdF1COggCivEGnSSxqEhPS2+/uD477mAIKSzaTcn+vyutzZ2Zl7YUn2npnnGZPVarUiIiIiIiKSf7Y4GJ1ARERERESKHxUNERERERHJdyoaIiIiIiKS71Q0REREREQk36loiIiIiIhIvlPREBERERGRfKeiISIiIiIi+U5FQ0RERERE8p2KhoiIiIiI5DsVDRERERERyXcqGiIiIiIiku9UNEREREREJN+paIiIiIiISL5T0RARERERkXynoiEiIiIiIvlORUNERERERPKdioaIiIiIiOQ7FQ0REREREcl3KhoiIiIiIpLvVDRERERERCTfqWiIiIiIiEi+U9EQEREREZF8p6IhIiIiIiL5TkVDRERERETynYqGiIiIiIjkOxUNERERERHJdyoaIiIiIiKS71Q0REREREQk36loiIiIiIhIvlPREBERERGRfKeiISIiIiIi+U5FQ0RERERE8p2KhoiIiIiI5DsVDRERERERyXcqGiIiIiIiku9UNEREREREJN+paIiIiIiISL5T0RARERERkXynoiEiIiIiIvlORUNERERERPKdioaIiIiIiOQ7FQ0REREREcl3KhoiIiIiIpLvVDRERERERCTfqWiIiIiIiEi+U9EQEREREZF8p6IhIiIiIiL5TkVDRERERETynYqGiIiIiIjkOxUNERERERHJdyoaIiIiUmRZrVajI4gUekb9O3E0ZK+3YbVamTp1Kg4ODpjNZqPjiBR6OTk5mM1mXn/9daOjiIgUuO+++44vvvgCLy8vlQ2Ru3BwcCAuLo5+/frh7+9fYPstNEUjJycHgNdee83gJCJFg8Vi4ZNPPsFiseDgoJOTIlKyxMfH06VLFx599FHbdwgRuT1HR0c2btzItWvXCna/Bbq3u3BycsLRsVBFEim0LBYLTk5ORscQETGEg4OD7WegvjuI3J2Tk1OBH5jUYVAREREREcl3Ra5oZGRkkJWVZXQMEbvKzMwkMzPT6BgiIiJFSlJSktER5A/sVjSys7OZO3cus2bNyrN8zZo1fPjhh0RFReVZ7uTkhMlk+tNtrl+/nhYtWvDYY4+xffv2fM8sUhhERkbSokULHn30UTZt2mR0HBERu4iNjWXs2LF89dVXtmXJycl88sknTJs2jeTkZA4fPswHH3xAeHj4La93dHTUIHCxuXT9Os/27k2dRo0YMHw4N3JzjY4k2HGMhqOjIx07dmTRokW2ZadOneL8+fO89NJLrFq1iocffhiAK1euEB0dTWpq6h23FxcXx4gRI4iOjgYgJCSEkJAQ/ZCRYsNkMmG1Wlm4cKFtsNaIESNo3bo1Hh4eBqcTEclffn5+tGzZkgsXLtiWffnll9SsWRMnJye2bdvG6dOnefnll1m8eDEtWrSgSpUqpKenc/nyZc6ePUvlypX/fCfnzsH8+aDvCsVbqVIkHjrEI1u30hzInTWLhPh4ytauDSVhooCcHGjeHHr2NDrJLexWNEwmE35+fri6utqW3bhxA29vb2rVqpVnhojo6GiOHj36p0UjOTmZ5ORk2+O0tDS8vb1xc3NT2ZBiwcHBgcTERNLS0mzLbty4QWpqqoqGiBRLPj4+JCQk2B4nJibSpEkTnJ2d2b17NwBVqlShfPny3LhxA4D09HSioqK4ePHi3Qe2njwJmzfDqFGQnW239yEGMZnAxYXUq1cJO3CAS4AzkA1ccnSkRs2axf/v3dERfv4ZVq0qWUUD4Pjx45w9e5YrV64QHx+Pt7c3ly5dYvPmzXh5ednWa9q0KU2bNmXmzJl33FatWrUYOnQoU6dOxdHRkTfffJP/+Z//sWd8EUMkJSUxceJEzGYzr7766t2P2ImIFEEWi4Xjx49z5swZrl69SlxcHLVq1WLfvn04ODhQt25djhw5wo4dO4iLi6NixYoAeHp68txzz+Hk5HRv09o2bgwDB9r53YhR1m3dypywMDwDAvjhhx+4cuYMAM/07g2dOxucroB8/z3MmGF0ituyW9GwWCxER0fTtGlTLl68CEDlypXp2rUrx44do1+/fnnWz87OxmKx/Ok2p0yZQlRUFC+++CJ9+vSxV3QRQ02YMIGjR4/StWtXQkJCjI4jImIXv1+lUKVKFWJiYsjMzOTJJ58kPT2d3Nxc2rVrxyOPPMLOnTvp2bMn5cuXz/P6rKysu47txGQqGZfOlEB79uxh0qRJmEwmPnrnHVq1akV8fDy//fYb165dY9yHH1L/gQeoW7eu0VHtLyPD6AR3ZLei4eDgQOfbNMlGjRrRqFGjv73d0qVLU7NmzfuJJlLoubq66nMuIsVauXLl6Nu37y3LO3bsaPt/b29vHViUPE6ePEloaCinT59m5MiR9OjRw/acj48PPj4+wM1L7ENCQggPD6dSpUpGxS3xitwdbqxWq+4AKsWePuciIiL/5+rVq8yYMYOvvvqK3r17M2fOnDzjgP/bCy+8QFxcHH379iUyMhJ3d/cCTCu/K3L30RARERGRkiEzM5PZs2fTuXNnsrOz2bhxIyNHjvzTkvG7V199lVatWtG/f38dvDOIioaIiIiIFDqRkZF07NiR/fv3s2zZMqZMmYKfn99f2sa4cePw8fFh6NChdkopf0ZFQ0REREQKjYMHD9K9e3fmz5/Pe++9x/Lly3nooYf+9vZmz55NYmIi7777bj6mlHuhoiEiIiIihjt79iyDBw9m1KhR9OjRg23bttG2bdv73q7ZbGbJkiXs37+fTz75JB+Syr1S0RARERERwyQmJjJu3Dh69uxJrVq12LFjB3369Ln79MV/gZubG8uXLycsLIw1a9bk23blzxW5WadEREREpOjLyclh2bJlLFy4kBYtWhAREWHXm9RWrFiRpUuX0qtXL7y9vWnfvr3d9iU3qWiIiIiISIHavn0706ZNw9vbm3nz5vHII48UyH4ffPBB5syZw9ChQ1m0aBFNmjQpkP2WVCoaIiIiIlIgDh8+TGhoKAkJCYwePZqAgIACz/DYY48xbtw4Bg8ezLp166hRo0aBZygpVDRERERExK4uXbrE1KlTOXjwIC+//DLBwcGYzWbD8gQGBhIfH0+/fv0IDw/H29vbsCzFmYqGiIiIiNhFSkoKs2fPJjw8nC5durB9+3bKlStndCwABgwYQGxsLMHBwWzYsAEXFxejIxU7mnVKREREDHH69GkiIyO5fv06AFlZWezYsYN169bx3Xffce3aNSIiIti+fTvp6el5XpufMxJJ/rNaraxcuZKAgADOnTtHWFgYY8eOLTQl43dvv/02Dz30EAMHDjQ6SrGkoiEiIiIFLjk5maVLl2KxWFi0aBFw834HtWvXJi0tjT179nDy5Em2b9+Ou7s7pUqVyvN6i8ViRGy5B7t376ZTp05ERETwySefMHfuXKpXr250rDuaNm0ajo6OvPbaa0ZHKXZUNERERKTARUdH4+3tTffu3UlOTiYrK8tWNJKTk+nUqRO1atWiTZs2fP311/z2228AxMbGMmfOHLZs2YKjo64AL0yOHj1Knz59+Oijjxg2bBgbNmygefPmRse6JwsWLOD06dNMmDDB6CjFioqGiIiIFLgyZcpw48YN4uPjsVqtpKWlAXDt2jUuX75M/fr18fLy4sUXX6ROnTr88ssvAPj6+jJs2DC6dOlCTk6OkW9B/iMmJobRo0fzyiuv8MQTT7Bjxw4CAwONjvWXlCpVimXLlrFz507bGTa5fzoUICIiIgWucuXK1K5dm08//ZTOnTvz7bff0qZNG86fP0+rVq0wm838+uuvbN68GScnJ0JCQoD/G5vh4KBjpUbLyMhg3rx5rFmzhvbt27Nlyxa8vLyMjvW3eXp6smLFCoKCgvDx8eHZZ581OlKRp6IhIiIihujduzdWqzXPwO4mTZrYbqJWr149/P39bzsNqtVqLbCccquwsDBmzZpFnTp1WLVqFbVr1zY6Ur6oXr06//rXvwgJCcHLy4tWrVoZHalIU9EQERERw9xt9igj77Ugt9qzZw+TJ0/GZDIRGhpaLL+IN2rUiI8//pjhw4ezatUq6tata3SkIktFQ0RERET+1MmTJwkNDeXMmTMMHz6coKAgoyPZ1VNPPcVbb71FSEgI4eHhVKpUyehIRZKKhoiIiIjc1rVr15g+fTpfffUVvXr1Ys6cObi6uhodq0C88MILxMbG0q9fPyIiInB3dzc6UpGjkVQiIiIikkdWVhZz5syhc+fOZGdnExkZyWuvvVZiSsbvhg8fTsuWLQkJCdEsZ3+DXYtGWloaMTExeZalpqZy+fJle+5WRERERP6mjRs3EhAQwPfff8/ixYuZMmUKFSpUMDqWYcaPH4+XlxfDhg0zOkqRY7dLp1JSUpg5cyY5OTk88cQTPPXUUyQlJfHZZ59RtmxZPD09eemll2zrm83muw4IExERERH7OHDgAKGhoWRkZPDee+/Rtm1boyMVGnPmzKFnz56MGTNGN/X7C+x2RiMqKooKFSowcuRI9u3bZ1uekZGB1WrNc/opMTGRc+fOkZGRYa84IiIiInIbZ8+eZfDgwbz++usEBQWxbds2lYz/YjabWbp0KQcPHmTGjBlGxyky7FY0LBYLjo6OODo62ua6vnHjBm5ubjz88MNcvXrVtm5UVBRbt27l+vXr9oojIiIiIn+QmJjIhx9+SK9evahduzY7duygT58+usLkDtzc3Fi+fDnr1q1jzZo1RscpEuxWNOrWrcvZs2eZNWsWjzzyCNu2bSMrK4u0tDQuXLiAi4uLbd1WrVoxfPhwTR0mIiIiYmc5OTksXryYZ555hoSEBDZs2MDo0aNxc3MzOlqhV6FCBZYsWcLUqVP56quvjI5T6NltjIaHhwcjRozg2rVr1KlTh9jYWPz8/BgxYgRxcXHUqVMnz/rZ2dlYLBZ7xREREREp8Xbs2MG0adPw8vJi7ty5NGrUyOhIRY6/vz+zZ89m2LBhLF68mMaNGxsdqdCy6300PD098fT0BMDPzw8Ab29vvL297blbEREREfmDw4cPExoaSkJCAm+88QYBAQFGRyrSWrZsybhx4xg0aBDr1q2jRo0aRkcqlHTDPhERETFETk4O165dsx2MhJvjBm7cuIGnpydubm4kJSUBUK5cOaNiFmkXL15k6tSpHDx4kEGDBtGvXz/MZrPRsYqFwMBA4uPj6devH+Hh4TqQfhsqGiIiIlLgLBYLs2bNIikpifr169OjRw8yMzMZO3Ys9evXp0OHDly/fp2FCxcC0L9/f2rWrGl7vbOzs22yGblVSkoKs2bNIjw8nMDAQLZv346Hh4fRsYqdAQMGEBsbS3BwMBs2bMgzBll0Z3ARERExwPnz58nMzGTMmDEcPnyYnJwcTCYTnp6eJCcnk5mZybfffkunTp3o1KkT33zzDQDXr19n8+bNfPvttzg66njpf7NaraxatYqOHTty7tw5wsLCGDt2rEqGHb399tv4+/vz8ssvGx2l0FHREBERkQKXm5trmwbfZDJhMplwdnbmww8/pGfPnoSFhZGZmYmLiwulSpUiNzcXABcXF/z9/alUqZImkfkvu3fvplOnTmzYsIEZM2Ywb948jR0oINOnT8fBwYHXXnvN6CiFig4FiIiISIGrVq0aGRkZzJgxgwceeIBdu3ZRv359Dh06xKVLl6hevTotWrRg3bp1WK1WunXrBoCrqysPPvggJ0+eVNH4j6ioKEJDQ7l8+TKjRo3i2WefNTpSibRgwQKef/55Jk6cyDvvvGN0nEJBRUNEREQKnLOzM8OHD+f8+fPUq1ePq1evUq5cOfz9/alTpw7+/v44ODjQr18/gFuOzP9+qVVJFhMTw7Rp0/juu+/o378/AwcOxMnJyehYJZaLiwvLli2ja9eu+Pn5MXDgQKMjGU5FQ0TEjrZv386ZM2du+9xzzz1HlSpVCjhR0XP9+nVWr1592+c8PT3p3bt3ASeS/FK2bFkaNGgA/N80+P7+/nnW0aU/t0pPT2f+/PmsWbOG9u3bs3XrVry8vIyOJdz8mbRixQqCgoLw9vbmueeeMzqSoVQ0RETsaP78+WzcuPG2zz344IMqGvcgJiaGV1999bbPPfTQQyoaUqKEhYUxa9Ys6tSpw6pVq6hdu7bRkeS/VK9enYULF9K/f3+8vb1p1aqV0ZEMo6IhImJnLVq0YMuWLbcs1yww9+bBBx8kLi7uluUjR47k559/NiCRSMHbu3cvkydPxmq1EhoaWqK/vBYFjzzyCNOmTWPEiBGsWrWKhx56yOhIhlDREBGxM0dHR3x8fIyOUWSZzebb/vlpvnopCU6dOsWkSZM4deoUw4cP5/nnnzc6ktyjdu3aMXr0aPr3709ERAQVK1Y0OlKBU9EQERERKWSuXr3KjBkz2LVrF71792b27Nkq10VQr169iI+Pp2/fvkRERODu7m50pAKloiEiUkicOHGCBQsWcPjwYaxWKw899BDdu3fn6aefNjra37Z06VKuXLmiqR5F7lFWVhaLFi1i2bJltG7dmk2bNtkGykvRNHz4cOLi4ggJCWHt2rWYzWajIxUY3bBPRKQQWLJkCcOGDaNTp06sWbOGGTNmkJycTPv27enSpQvJyclGR7wnN27c4OjRoyxfvpynnnqKkJAQDh48aHQskSJh48aNBAQEsG/fPpYsWcKUKVNUMoqJ8ePH4+3tzdChQ42OUqB0RkNExGD79u3j888/Z/v27Tg7OwPg6+vLihUrSE1NJSIigr59+xIZGWlw0rsbO3YssbGxNGrUiGeffZbdu3cbHUmk0Dtw4AChoaFkZGTw3nvv0bZtW6MjiR3Mnj2bnj17MmbMGCZMmGB0nAKhMxoiIgabMGECzz33HA4Ot/5Ifu+994CbRzp/+OGHgo72l82cOZPPP/+ct956i+rVqxsdR6RQO3fuHEOHDmXUqFEEBQWxbds2lYxizGw2s3TpUg4cOMDMmTONjlMgVDRERAy2Z88e/vnPf/LPf/7zlucaNWpkGzz45ZdfFnQ0EbGDpKQkxo0bR8+ePalRowY7d+6kT58+Jf5O5yWBm5sby5cvZ+3atXz++edGx7E7XTolImIwDw8Pbty4cdt7RQB4e3uTnJxMQkJCAScTkfyUm5vL0qVLWbRoEc2bN2fDhg26aWcJVLFiRZYsWULv3r3x9vYu0hN+3I2KhoiIwbZu3coXX3xx2ztc5+TkcOnSJQCqVq1a0NFEJJ/s2LGDadOm4eXlxZw5c2jUqJHRkcRA/v7+zJ49m2HDhrF48WIaN25sdCS7UNEQETFY/fr1qV+//m2f++qrr8jOzsbR0ZFu3boVcDIR+/riiy84dOgQzz33HPXq1SMzM5O1a9eSkJBAYGAgDg4OrF69mqpVq9K9e3fc3Nxsry0qU4QePnyY0NBQEhISeOONNwgICDA6khQSLVu2ZNy4cbzyyiusW7eOmjVrGh0p32mMhohIITZr1iwA3njjjWL5S0hKrri4OPbs2UP37t0JCwsDbpaHgIAAWrduTXh4OHFxcVy8eJFq1apRpkwZ4OZ9JmJjY4mNjS3UYxouXbrEyJEjGTZsGB06dGDHjh0qGXKLwMBAhg4dSnBwMNeuXTM6Tr5T0RARKaS2bt3Kli1bCAoK4qOPPjI6jki+un79Ot7e3vj7+wM3LxN0dHTEz8+PAwcO8Pjjj/PII48wYsQIDh48yN69e4Gb92r56quv+OWXX247U5vRUlNTmTx5Mj169MDb25tt27YxcODAInMGRgrewIED6dSpE8HBwWRmZhodJ1/Z9V/or7/+ys6dO/P8ocXGxrJ9+3YuXLhgz12LiBRply5dYuDAgXTu3JlVq1bpS4oUOz4+PsTHx7Nv3z6cnZ05c+YMaWlpfPzxx5hMJpo2bUp6ejqenp6UKVOGlJQU4ObkCC+++CIBAQHk5OQY/C7+j9VqZdWqVQQEBHDu3DnWrl3L2LFj8fDwMDqaFAFvvfUWDz74IC+//LLRUfKV3YrGlStXWLt2LTExMaxZswaAlJQUZs2aRWxs7C3rOzk5FcojEyIiBS0lJYWuXbvy1FNPER4ebruJn0hx4unpSWBgIL/88gvBwcHEx8eTlJSEk5MTJpOJQ4cOkZyczKZNm6hSpQrt27fP8/rs7OxCc+nU7t276dSpExs2bGDGjBnMnTuXGjVqGB1Lipjp06djMpl47bXXjI6Sb+w2GPz06dP4+/vTo0cPJk+ebFt29OhRqlWrxnfffUefPn0AOHr0KCdPniyW16aJiPwV2dnZPP/88zz++ON88sknOgAjxVqLFi1o0aIFAJUrVwZg5MiRedYZNGhQgee6V7/++isTJ07kypUrvPbaazz77LNGR5IibsGCBTz//PNMnDiRd955x+g4981uv8Hc3d1JSEggJiYGR0dHcnNzcXZ2pl69evTs2ZOoqCjbuqVLl6Z8+fK6NEBESjSr1crAgQNp2bIln376aZ6ScfXqVRYuXGhgOhH5XUxMDKNHj2bAgAE88cQT7NixQyVD8oWLiwvLly9nx44dLFq0yOg4981uRaN+/fo4ODiwePFiOnXqRHh4ONWrV6dOnTp88skndOzY0bbuAw88QNu2bXUdo4iUaG+//Tb169fnvffeu+W5/fv3Y7FYDEglIr9LT09n5syZPPfcczg7O7NlyxaGDBmiyxslX5UvX55ly5Yxb948Nm3aZHSc+2K3S6ccHR0ZNmwYFosFBwcHHnnkEUwmE/3797fNLPFH2dnZ+iUqIiXWp59+Snp6OkFBQZw5cybPc0lJSfzrX/9i1KhRBqW7P1ar1egIIvdt/fr1fPrpp9SuXZuVK1dSp04doyNJMVazZk3+9a9/ERISgre3N48//rjRkf4Wu9+w7/dT/38csPXfJUNEpCQLCwtj1KhRWCwWPv300zuuN2/evAJMdX9yc3P56aefADhz5gxJSUmUK1fO4FQif93evXuZNGkSAKGhobRq1crgRFJSNGrUiGnTpvHqq6+yatUq6tata3Skv0zf+EVEDDZ69Oi7fgl3d3enQoUKBZTo73vnnXfYuHEjABkZGTzwwAOkpqbSrFkz3NzcSE9P5/jx4wanFLm7kydPMnnyZE6dOsWIESMICgoyOpKUQO3atePNN98kJCSEiIgIKlasaHSkv0RFQ0TEYOfPnzc6Qr6ZOHEiEydONDqGyN929epVZsyYwa5du+jduzezZs3C1dXV6FhSgvXq1Yv4+Hj69u1LREQE7u7uRke6Z5o3UUREREq8rKws5s2bR5cuXcjMzGTjxo289tprKhlSKAwfPpxHH32UkJAQcnNzjY5zz1Q0REREpETbtGkTHTt2ZM+ePSxevJiPP/64SFyqKCXLhAkT8PLyYujQoUZHuWcqGiIiIlIiHTx4kO7duzN37lzeffddVq5cSb169YyOJXJHs2fPJiEhgXfffdfoKPdERUNERERKlLNnzzJ06FBGjRpFjx492LZtG+3atTM6lshdOTo6snTpUvbv38/MmTONjnNXGgwuIiIixZOj483//iMxMZHPPvuMLVu20L17d6ZOnYqbm5uBAUX+Ojc3N1asWEG3bt3wrVqV3lWqgNlsdKzbUtGQP5WamsrmzZsxm80EBgbi4uJidCQRESkmrly5wvHjx3n00UcpU6YMANHR0Rw7dsy27KeffsJqtdK0adM8r/3j/bluxwKcOX6crCNHSPz5Z347coS5c+fSokULwsPDqVy5sr3eVol28OBBfvrpJ5588kkeeugho+MUWxUqVGDp0qUEv/IKro0b8/Dx41w7fJjHGjUyOloeKhryp4YMGcLKlSsBGDRoEPPnzzc4kYj9xMXFMWvWLJ5++mlatWqFuZAeIRIpCMuXLyctLY2goCC8vb3zfftpaWnMnz+fGjVqEBUVxfDhw0lPT2f+/PlUr16dEydO0KJFC7788ktMJhPZ2dk89thjttffrWisXreO5cOG8RIQ3KQJnTp1Yv78+TRs2DDf34vctGvXLtusXX5+fnzzzTcqG3bk7+9Pt+7dGT9qFK8DfR9/nB+++YYWLVoYHc1GYzTkjqKjo9m6davt8caNG0lISDAwkYh9+fr68o9//INnn30Wf39/2rdvz4wZM4iNjTU6mkiBe+GFF5g7dy7+/v40b96ckJAQ9u3bh8ViyZftX7p0ifLlyxMSEkJ8fDwWi4WYmBjc3d0JCQkhISGB77//nsDAQDp16sSRI0cAiI2NZd68eWzatAlHxzsfL12/YQOpgPN/Hrdp00Ylw842bdpEZmYmJpOJ2NhYdu/ebXSkYu/o0aM4cfNzbk1PZ/OWLUZHykNnNOSOfHx8aNy4se0HRbNmzfD09DQ4ldyrPn36cPHiRaNjFElVqlTh5MmTnD59ml27dvHWW2/h4eHB4MGDefPNN3VNdxETGRnJ9OnTjY5RJLm7u5OZmcmPP/7Ijz/+yIoVK3B2dqZ169Z89NFHNGvW7G9v28XFhbS0NJKTk7FarWRlZeHo6Eh6erptWdmyZbl69aptfbj5u2nQoEF4enqSk5Nzx+3/44kniFi7lqz/PG7evPnfzir3pmXLlnz22WdYrVacnJxo3Lix0ZGKLYvFwqpVq/j222+pDrbP+eMtWxoZ6xYqGnJHv89sMH/+fMxmM0OGDDE6kvwFr732Gunp6UbHKHKys7N56623yM3NpVy5cvj4+ODn58cLL7xASEiISkYR9Oijj/LRRx8ZHaNI+vrrrzlx4gRZWVn4+Pjg6elJw4YNGTt27H1fElOtWjX8/PyYPHky7dq1Y8eOHTzxxBNUrVqVyZMn06ZNGxo3bsycOXOwWq0MGjQIAAeHmxdjlCpV6k+3P3ToUJpGR+P6+ec0qViREydO0KZNm/vKLH+uV69eWK1W9u7dS5cuXWhZyL70Fhe7du1i6tSplApYMvIAACAASURBVC1blrVr15K9eTNOCxeybNo0nnnmGaPj5aGiIX+qatWq+gVdROno3V+XlZVFUFAQFouFMWPG0K9fP2rXrm10LLlPFStWpGLFikbHKHI2bdrE0qVLad++PS+99BJPPfXUXb/c/1UDBw4kIyMDFxcXcnNzMZvNBAcH25YBvP322wC3XCZltVr/dNulzWbadOgAly6xacYMuj/zDD4+PnTv3j1f34P8H5PJxIsvvsiLL75odJRi6ciRI4SGhhIbG8vrr79O586dbz6RkQHHjtHyhReMDXgbKhoiIv+RmprK2rVrcXV1NTqKiOGaN2/O2bNn7b6f3wvFHydf+OMMh382DuOusrIgK4vKnp4sXLiQfv364e3tTevWrf/+NkUKWHR0NNOmTWPv3r3079+fAQMG4OTk9H8rZGZCPo2dym8aDC4i8h/ly5dXyRD5j+J2FqhBgwbMmDGDkSNHEhUVZXQckbtKS0tj+vTpdO3aFVdXV7Zu3crgwYPzloxCTkVDRKSQSUlJISQkhH/9619GRykwu3btolmzZne9HEbkfrRp04axY8cyYMAATZYhhdrnn39OQEAAv/76K6tXr2b8+PFFckIeXTolIlII/PTTT5w/f569e/eyZs0arly5QvXq1Y2OZTexsbGcOHGCY8eOERkZyRdffIHFYiE3N/f+LpURuYvu3bsTGxtLcHAwkZGRlC1b1uhIIjb//ve/mTJlCmazmSlTphT5AfX6aS4iUggsWbKEsmXL0rZtW1xcXJgwYYLRkezqyJEjbNy4kdq1azN79mz8/f3z7f4MInczdOhQ4uLi6NevH+vWrStSl6JI8XTixAlCQ0M5d+4cI0aMKDaTFqhoiIgUAp999pnt/w8ePGhgkoLx9NNP8/TTTxsdQ0qw999/n1dffZXBgwezePFio+NICRUfH8/06dP5+uuv6dOnD/Pnz8/32d2MpDEaIiIiUiJ9+umnpKWl8dZbbxkdRUqYzMxMZs+eTWBgIBaLhU2bNjF8+PBiVTJARUNERERKKAcHBxYvXszPP//MtGnTjI4jJURkZCQdO3bkwIEDLFmyhMmTJ+Pr62t0LLuw66VTly5dIiYmhiZNmtju5Alw6tQp/Pz8NABLREREDFW6dGlWrFhBt27d8PX1pW/fvkZHkmJq//79hIaGkp2dzfvvv18i7lRvtzMa8fHxLFiwgH379rF+/Xrb8t9++43Ro0dz7NixPOs7OTlhMpnsFUdERETktnx9fVm2bBmffPIJO3bsMDqOFDNnzpxh0KBBvPHGG/Tq1YutW7eWiJIBdjyjcfLkSR588EF69OjB1KlTAUhMTOTrr7+mQ4cOeeZK//rrrzlw4ABXr161VxwRERGRO6pduzbz5s1j0KBBeHl50bx5c6MjSRGXmJjIjBkz2LFjB88//zwzZsygTJkyRscqUHY7o1GqVClSU1NJTk62XTZ19uxZTp48yTfffMO3335rW/exxx5jwIABxfb6NBEREblVSkpKnse5ubncuHHDdjAyMzOT69evk5GRUSB5mjdvzoQJExg6dChnzpwpkH1K8ZOTk8PChQt55plnSEpKIjw8nDfeeKPElQyw4xmN+vXr88033zBr1iwCAgIICwsjMDCQqVOnsnbtWmrWrGlb19XVFWdnZ81jLSJSCKWkpLB9+/a//fpy5crRoUOHfEwkxcGCBQu4cOECLVu2pHPnzgB89913/Pjjj5jNZl599VXee+89KlSoQPv27alfv36e1zs5OdnlTvKdOnUiPj6efv36ERERgY+PT77vQ4qv7du3M23aNHx8fJg3bx6PPPKI0ZEMZbei4eLiwqhRo0hLS8Pd3Z3MzEycnZ0BeP75528Zj2GxWOzyA0NERO7PtWvXiIiIICcn52+9vmLFiioaJZzFYiEzMxMAs9nMxYsXuXr1Kh988AHvv/8+7du3x9nZmSeeeII2bdowZswYEhISMJvN5OTk5PnOkJSUxE8//cSBAwfo1q2bXfIGBwcTFxdH3759CQ8Pp3Tp0nbZjxQfP//8M6GhoSQmJvLmm2/qZ95/2HXWKbPZjLu7O0CeeYH/OAOViIgUbtWrV2f16tVGx5AiLDo6mrCwMHJzczGZTFSvXh13d3ccHR0xmUzk5uYC4OjoSGRkJE2aNMHPz4/x48dz6dIl5s6dy6RJk4Cb3y3Kli1r9y//o0eP5urVqwwcOJDVq1drwhq5rQsXLjBlyhR++uknBg0aRL9+/fQ99w90Z3ARERGxq8qVKzNq1Cjb48zMTCZPnsyCBQuoXLky165dIy4ujgsXLhAeHs6QIUNISkri0KFDXL58Oc8YTjc3N5o2bUpMTIytoNjL5MmTCQkJYcSIEXz22Wd23ZcULcnJycyaNYuIiAi6du3KhAkTKFeunNGxCh0VDRERESlQpUqVYtiwYZw4cYJmzZqRmZmJ2WymXLlyeHp6AjfHYHh5eVG6dGleeOGFW7aRnZ1dIGcZFixYQFBQEB9++CHvv/++3fcnhZvFYmHlypXMmzePxo0bs379eqpVq2Z0rEJLRUNERAz3+xg9jdUrOby9vfH29gZuFo/fb+Jbq1Yt2zqFYSCtk5MTy5Yto1u3blSoUIHBgwcbHUkM8tVXXzF16lTKlCnDp59+SrNmzYyOVOipaIiIFBJJSUkkJCSwd+9eAH744QfOnTuHh4cHHh4eBqfLf5mZmaSmpvLvf//bdgnMpk2baNu2LWXKlMkztk/ESB4eHixfvpwePXrg6+trt0HoUjgdPXqU0NBQoqOjef311+nSpYvRkYoMFQ0RkUIgODiYqKgoypcvD8DTTz9Nbm4ur7zyCtevX6devXosX77c4JT55/vvv2f48OG4uLjg6urK008/DcC8efP45JNPSE9PZ9q0aTz55JMGJxW5qWrVqixatIjg4GC8vLxo3bq10ZHEzqKjo/n444/Zt28f/fv3Z8CAAboVw1+koiEiUggsW7bM6AgFqmXLlvz4449GxxD5Sxo0aMD06dMZOXIkK1eu5OGHHzY6kthBWloa8+bNY+3atQQEBLB161bb2CH5azT/loiIiMg9atOmDWPHjmXAgAFcvHjR6DiSz34vF8eOHWP16tWMGzdOJeM+6IyGiIiIyF/QvXt3YmJiCA4OJjIy0jaQXYquPXv2MHnyZEwmE5MnT+bxxx83OlKxoKIhIiIi8hcNGzaMuLg4goODCQsL07X7RdSJEyeYNGkSZ86cYcSIEfTo0cPoSMWKLp0SERER+Rs++OADKlasyJAhQ4yOIn9RfHw877zzDsHBwTRu3JidO3eqZNiBioaIiIjI3zRr1ixSUlJ48803jY4i9yAzM5PZs2cTGBhIbm4umzZtYsSIEbi4uBgdrVhS0RARERH5mxwcHFiyZAm//PIL06ZNMzqO/ImIiAg6duzI/v37Wbp0KZMnT8bX19foWMWaxmiIiNhZdnY2cXFxtywvX768ruu+B7m5uVy7du2W5RkZGQakEblV6dKlWb58Od26dcPX15e+ffsaHUn+4IcffmDSpElkZ2fz/vvv06ZNG6MjlRgqGiIidnbw4EH8/PxuWf7FF1/Qvn17AxIVLSdPnqRevXq3fe6hhx4q4DQit+fr68uyZcvo2bMnvr6+BAQEGB2pxDtz5gyTJ08mKiqKYcOG8eKLLxodqcRR0RARsaMRI0bQtWvX2z5Xv379Ak5TNFWqVIklS5bc9jkPD48CTiP55fvvv+fAgQMEBgbywAMPAHDo0CH27NlDzZo1efbZZ22Pn3nmGR588ME8r3dwKHxXf9euXZu5c+cyePBgvLy8aNasmdGRSqTExERmzpzJtm3bCAoKYvr06ZQpU8boWCWSioaIiB21a9fO6AhFXrly5ejfv7/RMSQfJSUlsW3bNrp168aKFSt4//33Afjll1+4cuUKTz31FKmpqURERNCjRw9WrlzJhx9+iMlkIicnh4yMDJKSkqhSpYrB7+RWLVq0YMKECQwZMoSwsDBbiRL7y8nJYcmSJSxevJhHH32UiIgIKlWqZHSsEk1FQ0REROwqNjaWrVu3YrFYMJlMeHp6UqFCBZo0aUJkZCQpKSm4ubnRvn17atWqxfr162nXrh1eXl40btyYjRs3kpSUhIeHBwkJCWzZsoWDBw/SsGFDo9/abXXq1In4+Hj69etHREQEPj4+Rkcq9rZt28b06dPx8fFh/vz5hfazUdKoaIiIiIhdlStXjrZt22K1WjGZTAAcPnyYX3/9FbPZTG5uLjExMbi7u9O0aVP+/e9/4+DgwI0bNzh27Bi5ubm2u2/7+voyYMAAfH19ycnJMfJt/ang4GDi4uLo27cvERERuLq6Gh2pWPr555+ZOHEiSUlJvPnmm3To0MHoSPIHKhoiIiJiVy4uLtSsWTPPsnbt2rF161b69OlDUlIS0dHRmM1mvvvuOxo1asQTTzyBg4MDmzdvpmfPnreMycjNzbWVlsJq9OjRxMfHM2DAAFavXl3o8xYlFy9eZMqUKRw6dIhXXnmFfv36YTabjY4l/0VFQ0RERArck08+yZNPPml7XK1aNYA8A6gff/xxHn/88du+3mq12jdgPpkyZQohISGMHDmSTz/91Og4RV5KSgqzZs0iIiKCwMBAtm/fTrly5YyOJXdQ+KZsEBERESlG5s+fz/nz5xk/frzRUYosi8XCihUrCAgI4MKFC4SFhfHuu++qZBRydj2jsWvXLs6cOUOPHj3w8vIiJSWFjRs3kpubS2BgIOXLl7et6+DgoFOKIiIiUuw4OzuzbNkyunbtiq+vL4MHDzY6UpHy1VdfMXXqVNzc3Jg5c6amDS5C7HZG47fffmP//v3UqlWLVatW2ZY/9thjlCtXjg0bNtiWpaenc+3aNbKzs+0VR0RERMQwHh4erFixgkWLFuX5DiR3dvToUfr06cOECRN49dVXWb9+vUpGEWO3MxpXrlyhVq1aPPHEE+zduxcANzc33Nzc2Lp1Ky1btrStu3//fg4cOEB8fLy94oiIiIgYqmrVqixatIjg4GB8fHxo3bq10ZEKpejoaD7++GP27dtH//79GTBgAE5OTkbHkr/Bbmc0fHx8uHDhAkeOHKF06dLExsaSlZXF9OnTqVixIs2bN7et26ZNG958800qV65srzgiIiIihmvQoAHTp09n5MiRREVFGR2nUElPT2fGjBl07dqV0qVLs2XLFgYPHqySUYTZrWjUrVuXWrVq8e2339KzZ09++uknrly5QnR0NNHR0Rw4cCDP+tnZ2UVmBgkRERGRv6tNmzaMGTOGAQMGcOnSJaPjFAphYWEEBAQQFRXFqlWrGD9+PF5eXkbHkvtk18HgPXr0sP1/1apVAZg6dao9dykiIiJS6AUFBREXF0e/fv2IjIy03ZCwpNmzZw+TJk3CwcGBSZMm3XE6YymadB8NEREREQMMGzaMuLg4goODCQsLK1GXCJ04cYLQ0FDOnj3LiBEj8hycluJD99EQERERMcgHH3xAxYoVGTJkiNFRCkR8fDzvvPMOwcHBNG3alJ07d6pkFGMqGiIiIiIGmjVrFsnJybz11ltGR7GbrKws5syZQ5cuXcjJyWHjxo0MHz4cFxcXo6OJHaloiIiIiBjIwcGBpUuXcvjwYaZNm2Z0nHwXGRlJQEAA33//PUuXLmXKlCn4+fkZHUsKgMZoiIiISIFLTEzk9OnTNGjQAGdnZwAuXrzIxYsXKVWqFPXq1ePy5cvExcVRo0YNKlWqlOf1JpPJiNh2U7p0aZYtW0b37t3x8/PjpZdeMjrSfTtw4AChoaFkZmby3nvv0bZtW6MjSQHTGQ0REREpULm5uXz22Wfs27ePxYsX25ZnZGRw/fp1pk2bxrVr15g6dSpRUVHk5OTcsg2z2VzspsX38/Nj2bJlzJgxg507dxod5287c+YMgwcP5vXXX+eFF15g27ZtKhkllIqGiIiI2FVqaipHjhzhl19+4ejRoxw8eBA3NzeGDx/OpUuXSEtLA6BOnTo0a9aMhx9+mCpVqtCkSROuX7/OL7/8YttWXFwcixYtYsuWLTg6Fr8LM2rXrs3cuXN56623+PHHH42O85dcv36dDz74gN69e1OrVi127txJ7969jY4lBip+/0JFRESkUElLS+PXX3/FYrEA4OPjYztLYbVacXBwwGq1YjKZ2LJlCy1atABg8ODBJCUlMW7cOLp06QKAp6cnvXr1olSpUuTm5hrzhuysRYsWTJw4kcGDBxMWFkatWrWMjvSncnJyWLp0KYsWLeKxxx4jIiLilkvdpGRS0RARERG78vHxoWfPnnmWnThxgvfee4/HHnuMy5cvc/bsWZ5++mni4uLo3bs3mZmZrF27litXrtCyZUvb6xwdHXF0dCz2N7h75plniI2NJTg4mIiICHx8fIyOdFvbt29n+vTpeHl5MX/+fBo2bGh0JClEVDRERESkwA0dOpTExEQ8PT3JycmhSpUqWK1WXn/9dUqVKgXAs88+i8ViwdPT85bX/352pDjr37+/7e7h4eHhuLq6Gh3J5ueff2bSpEkkJCTwxhtvEBAQYHQkKYRUNERERKTAOTg42ArE72cpAFvJAPDw8DAkW2Hy5ptvcvXqVQYMGMDq1asNn23r4sWLTJkyhUOHDvHKK6/Qr18/zGazoZmk8NJgcBEREZFCbMqUKbi4uDBixAjDMqSkpDBp0iSCgoLw8/Njx44dhISEqGTIn1LREBERESnk5s+fz4ULFxg/fnyB7tdisbB8+XICAgK4cOECYWFhvPvuu8V+jIzkD106JSIiIlLIOTs7s2zZMrp27Yqvry+DBw+2+z537drFxx9/jJubGzNnzqRZs2Z236cULyoaIiIiIkWAh4cHy5cvJygoCB8fH7p3726X/URFRREaGsqVK1cYNWoUgYGBdtmPFH8qGiIiIiJFRLVq1Vi4cCHBwcH4+Pjwj3/8I9+2HRMTw8cff8yePXvo378/AwcOxMnJKd+2LyWPxmiIiIiIFCENGzZkxowZjBgxgqioqPveXnp6OtOnT+e5557D1dWVrVu3MmTIEJUMuW8qGiIiIiJFTJs2bRgzZgwDBgzg8uXLf3s7a9euJSAggF9//ZVVq1Yxfvx4vLy88jGplGS6dEpERESkCAoKCiI+Pp6+ffsSERFBuXLl7vm13333HZMnT8ZsNjNp0iQef/xxOyaVkkpFQ0RERKSIGjp0KDExMYSEhBAWFma78eGdnDx5ktDQUM6cOcPw4cMJCgoqoKRSEhXJS6ccHIpkbJG/RJ9zERG5Fx9++CEVKlRgyJAhd1zn6tWrjBkzhn79+tGoUSN27typkiF2Z/dvMrm5ufe07F5ZrVZSU1PvJ5JIkaDPuYgUZ1ar9bbLLRZLnsc5OTkFEafImzVrFjdu3ODtd97BYrHYxm1kZWUxZ84cOnfuTHZODhs3bmTkyJG4uLgYnFhKArtdOpWTk8P8+fOJjY0lKCiIhg0bkpOTw4IFC4iOjqZHjx40atTItr6Tk9Ndj+DOnz+fnTt38v333zN16lT69Oljr/gihlm0aBHbtm3j22+/ZdKkSfTv39/oSCIi+W7lypX87//+Lx9//LFt2enTp1mxYgXu7u4MHz6cbdu2cfDgQR577LFb7uXg5OR0x7JiY7XCXS4lKi4cHBxYGhZG165dqV63LtevX+fRRx/FwcGBihUrsjQsjLrVqxsdU+yhEJdGu/3rO3LkCFarlZdffpmVK1fSsGFDoqKiyMnJ4ZVXXmHFihW2onH69GnOnTtHYmLiHbd3+fJlxowZQ2pqKqmpqfzzn/9k165d9oovYgir1cr69etJTU0lJSWFd955h86dO+Pj42N0NBGR+/LHMxUODg506dKFY8eOkZOTYxtXsHXrVoKCgjhy5AiRkZGcPHmSd999l48++oi2bdvi5uZGSkoKx48f58iRI1SuXPnPd+rgALt3w8svw31cTVFUlHZ2ZmZMDF+fPIkzkLllC/Vr1uRJPz+YMAGys42OKPnNbIbTp8HDw+gkt2W3opGSkoKnpyd+fn62055paWl4eHhQoUIFsv/wYc/MzCQ5OflPL6nKycnJc/rU0dGRDh064ObmdstpVpGiyGw2c/36dSIiImzLcnJy7utSQxGRwuDy5ct8/vnn5ObmkpOTQ2BgIA0aNMDV1TXP4OX09HQqVapEbGwsx48fx8XFBVdXV0qVKkVaWhpubm7k5uaSkpJCRkbG3XfcpAnMn18iSgYmEzg58UNcHBuA0kA6UO/ZZ6FjR8jMNDig2I3VCjVqGJ3ituxWNGrVqsWuXbtYs2YNVatWZf/+/VSqVImdO3eyevVqqv/h9F29evWoV68e58+fv+P2qlevzttvv83EiRMxm82MGzeO3r172yu+iGHi4uIYP348JpOJd955hwoVKhgdSUTkvlSqVImRI0faHjs6OvK///u/REVFceTIEfz8/IiOjqZevXqEhYVx7do1nnvuOTZv3kxYWBgODg74+voCUK5cOdq0aUNqaurdD8R4ekKXLvZ8a4VOs9q1mXLuHMd//ZV/PPkk9T78EP7CtLci+cluRaNSpUr06tWLixcv0rp1a06dOkW1atXo3bs358+f58knn8yzfnZ29l3PTPy///f/6NWrF2azmSpVqtgruoih/ud//oegoCBMJhPVqlUzOo6IyH0zmUy3TLualZVFUFAQ6enpODs74+rqSmBgILt378bT05P69evj6+vLgQMHbjubUnZ2NiaTqaDeQpHRoE4dDu3dy9mzZ/H397/rdLci9mTXT1/dunWpW7cuAA0aNADA398ff3//v73N6hrIJCWAPuciUtw1b96c5s2b2x57/Oca83bt2tmW+fr60qWEnZHID6VLl+bhhx82OoZI0byPhoiIiIiIFG6FqmjoBmUi987BwUGXDYhIiWUymXBycjI6hkiRYcRldIXmwj2TyURMTAw//PDDHQuH1WrFbDYDN2/6V1i+ZJlMJhwcHMjJySk0meDmLEaFbcaiwpbJarXi6OiIxWK5+3zsBeReP+c5OTnEx8cXqs+ciEhByc7O5uDBg2RnZ//pTf2M/L2jfWvfhWXfzs7OHDlyhJYtWxZgKjBZC8u3K+Df//43iYmJdywaTk5O7N+/n1KlStGoUaM8U+QaxWw2c+HCBU6dOkWHDh3IysoyOhImk4m0tDR27dpF165dC80XewcHByIiInjmmWcoVapUofhi7+TkxI4dO6hfvz6VKlUqFH9WTk5OHDx4ELPZTJMmTe74ObdYLHh5edGqVasCTigiYrxLly7x448/3vEorclkIjMzk+3bt9OtW7cCnQpf+95O165dC/T3vMlkIisri23btmnfd2CxWPjHP/5B+fLlCyrelkJzRgOgdevWd13H1dUVV1fXPAPIjHb+/Hl8fHzo0KGD0VFssrOzSUxM5JlnnjE6Sh4XLlygW7duheoofGJiIq1bt6ZixYpGR7Fxc3PDwcGhwI88iIgUFVWqVLnrDJRWq5X4+Hg6depUQKkKz77j4uIM3Xfnzp0LfN8AsbGx2nchUqjOaNyL5ORkHBwcKFOmjNFRbDIzM0lLSyvIhnhXFouFa9euFbo7SsfHx+Pl5VWoxuMkJCTg5uaGs7Oz0VFskpOTMZlMuLm5GR1FRKTIMvJ3ofZd8Pu2Wq1cvXpV+y48thS5oiEiIiIiIoXeFvMHH3zwgdEp/kxGRgbLly/n0KFD1KtXD0dHRzIzM1mxYgU//vgjdevWNWTWiZMnT7JkyRKysrKo8Z/bvn/xxRds3ryZ3NxcQ260du3aNRYtWsSFCxfyzJ/922+/sXbtWpo1a1bglyxlZmaycuVKDhw4kOfv6ssvv2T37t1UqFCBsmXLFmgmgE2bNrFt2zaqV6+Ou7s7FouFTZs28eWXX1KmTBnbHWgL0s6dO9myZQtNmza1DQb/6aefWLFiBaVKlaJSpUoFnklEpKj54YcfWL16Ne7u7vj5+QGwf/9+1qxZg5ubGxUqVLDbvq9fv87ixYs5ffo0Dz/8MCaTicOHD7N8+XIuXrxI3bp17XZG/4svvmDTpk15fof8/PPPrFixAicnJypXrmyX/QJcvnyZOXPmULZsWduf+fHjx1m0aBHnz5+3240DLRYL69evZ9euXbi7u9t+dyclJbF48WJOnjzJww8/bJc/8zvt+5dffmHZsmVcuHDBrn/fZ86c4fPPP+e3337jwQcfxNHRkRs3brBo0SK7vu+/6KThCe7myy+/xMXFBXd3d7Zv3w7Arl27cHJyoly5crZlBW316tUEBASwe/du4uLigJs3KOzZsycRERGkp6cXeKawsDDq16/PlStXOHToEIBtgNDx48cLdEDY777++mtMJhPe3t5s3boVgB9//JGNGzfi4eFhyKVBJ0+eJCoqipYtW7JmzRoA4uLi2L9/P7Vr12bHjh0FngmgYcOGnD9/nrS0NOBmyQ4PDycwMJCIiAgyMjIMySUiUlSkpKSwefNmunTpwrp168jJySEtLY3NmzfTuXNn1q1bZ9dJW8LDw6lTpw7Xr1/nhx9+AODUqVPExsZSq1YtWwGwh4YNG3LhwgVSU1OBmwf6wsPD6dKlCxs3brT9brEHHx8ffHx8OHHihG3Z2bNnuXTpEg888AClSpWyy36tViuPPvooHTp0sP0+B4iIiKBGjRokJyezb98+u+wbuO2+T506RUxMDLVq1bLrdLK+vr506tSJY8eO8dtvvwEQGRlJ9erVSUlJYe/evXbb919R6ItGbGws9evXp0GDBrYv9HFxcdSrV48GDRoQHx9f4JmysrLIzc2lQYMGeHt7c/36dQCqVq3KsWPHqFu3Lq6urgWeKzExkSZNmlCrVi3bn9WyZcsoXbo0JpOJhISEAs8UExPDww8/TMOGDW1/V8eOHaN06dLk5ubaykdBio2NpUaNGjRp0oSUlBT+f3v3/9LUHsdx/FlzJ78w2ZcGbrEtQ2TZwwAABKNJREFUt0a1X4b9UEE/2Ez9IQKFgogI+sf6I8Ir9UNKRqaLakcslzfvRMUvGya66VLOpveH7j1wL3Eh2Jl2eT1+Gvvl/T5ng/fnfT5fDoDH46GlpYXXr187+rTrv4RCIfx+vz3rtLe3R2trq/1/+jtXERH5sZ2dHTweD6lUCrfbbe+hbG9vJ5VKYRiGowPura0tenp6SCaTdh3OZDI8ePCA0dFRPn/+7Fjsrq4uAoGAXUOq1SqGYZBKpWhvb3e0hhiGwfnz5/8xsL5+/TqPHz9mbGwM0zQdietyuYhGo+RyOa5du2Z/v7m5STqd5uLFixSLRUdinz59+oexe3t7efjwIc+ePWNubs6R2PD90Jh8Pk+5XMbr9QLNue6fdeIbjXA4jGma5HI5fD4fpVKJcDjM7OwsuVzuWAaFhmHgdrt5//49W1tbHBwcUKlUeP78OePj4wwPDzc9J4BAIMD09DQLCwu0tbVRLBbtabuFhQW2t7ebntO5c+eYmZnhw4cP+Hw+e5AfiUTo6uqiUqk0PadQKMTi4iLT09N0dnayurrK4uIibrebR48ekc/nm54TQLlcplgs8vXrV1ZWVjg8PMSyLEzTZH9/H4/Hcyx5iYj8KrxeL7u7u5imyeHhIevr6xwdHbG/v08ul6NWqzl6mEwwGCSbzTI/P09HRwcbGxu4XC4SiQStra2Uy2XHYlcqFTY2NuwaUq/XqdVqmKbJt2/fHK8hq6urlEoldnd3WVpa4tSpUyQSCTo6OtjZ2XEs7pMnT1hfX2dgYMC+B6FQiHfv3jE3N+fosuN/x15eXsblchGPx2lra3P0ujc3N7l69Srd3d32eLRZ1/0zTvwejVgsZg/8+vv7+fjxI319feTzeer1OsPDw8fypsNoNMqLFy/IZDL2U5N8Pm/PHCSTSUenSH8kkUgwOTlJd3c3ly9fZmlpiRs3btDT00M4HObKlStN36MRjUb58uULlmUxODjI7OwsN2/eZG1tjbW1Ne7du9f02R+/30+1WmV+fp779+/z6dMn0uk01WqVmZkZhoaGjuXkhsnJSYrFon0OeTAYJB6PMz4+zp07dxxdXysi8n9gGAZer5dXr15x9+5dSqUSwWCQSCTCxMQEQ0NDjj6gjMfjZLNZQqEQ6XSalZUVDg4OGBkZ4dKlS/T29jpWh9+8eWPXEMuyCAQCXLhwgbGxMW7fvk0kEnEkLnxfKTA1NYVlWXR2drK9vc3R0RFPnz4lFovR39/vyH6BWq3GxMQEZ86cYW9vD6/XS6FQ4NatW2SzWc6ePcvAwIAj97xer/Py5Us7ts/nY3l5GcuyGBkZIZlMkslkHPu9C4UCo6OjxGIx0uk0hUKBvr4+3r59i9/vZ3Bw8CS8SuB3nTolIiIiIiKN9tuJXzolIiIiIiK/HjUaIiIiIiLScGo0RERERESk4dRoiIiIiIhIw6nREBERERGRhlOjISIiIiIiDadGQ0REREREGk6NhoiIiIiINJwaDRERERERaTg1GiIiIiIi0nBqNEREREREpOFagKm/PieAduDo+NIRkRPoFLALLB53IiJyYniBKHB43ImIyIn0B5D/E+Pnl+47PnslAAAAAElFTkSuQmCC"
    }
   },
   "cell_type": "markdown",
   "id": "35f7152d",
   "metadata": {},
   "source": [
    "![box_domain.png](attachment:box_domain.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd0f95db",
   "metadata": {},
   "source": [
    "The exact operation when multiplied with the weight matrix should lead to the rotated rectangle. However, in the interval domain we only consider the maximums and minimums of each feature, thus resulting in the larger red rectangle which contains many excess regions.\n",
    "\n",
    "Let's see how this does in practice!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c6a8846d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using device  cuda:0\n"
     ]
    }
   ],
   "source": [
    "import torch\n",
    "import torch.optim as optim\n",
    "import numpy as np\n",
    "\n",
    "from torch import nn\n",
    "from sklearn.utils import shuffle\n",
    "\n",
    "from art.estimators.certification import interval\n",
    "from art.utils import load_mnist, preprocess, to_categorical\n",
    "\n",
    "device = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\n",
    "print('Using device ', device)\n",
    "\n",
    "bound = 0.05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "00ac4fad",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We make an example pytorch classifier\n",
    "\n",
    "class MNISTModel(torch.nn.Module):\n",
    "    \"\"\"\n",
    "    The base model which we will then convert into a certified classifer.\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, number_of_classes: int):\n",
    "        super(MNISTModel, self).__init__()\n",
    "\n",
    "        self.device = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\n",
    "\n",
    "        self.conv_1 = torch.nn.Conv2d(in_channels=1,\n",
    "                                      out_channels=16,\n",
    "                                      kernel_size=4,\n",
    "                                      stride=2)\n",
    "\n",
    "        self.conv_2 = torch.nn.Conv2d(in_channels=16,\n",
    "                                      out_channels=32,\n",
    "                                      kernel_size=4,\n",
    "                                      stride=1)\n",
    "\n",
    "        self.fc1 = torch.nn.Linear(in_features=3200, out_features=number_of_classes)\n",
    "\n",
    "        self.relu = torch.nn.ReLU()\n",
    "\n",
    "    def forward(self, x: \"torch.Tensor\") -> \"torch.Tensor\":\n",
    "        \"\"\"\n",
    "        Computes the forward pass though the neural network\n",
    "        :param x: input data of shape (batch size, N features)\n",
    "        :return: model prediction\n",
    "        \"\"\"\n",
    "        x = self.relu(self.conv_1(x))\n",
    "        x = self.relu(self.conv_2(x))\n",
    "        x = torch.flatten(x, 1)\n",
    "        return self.fc1(x)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b239e7da",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup the MNIST Model\n",
    "\n",
    "model = MNISTModel(number_of_classes=10)\n",
    "model.to(device)\n",
    "opt = optim.Adam(model.parameters(), lr=1e-4)\n",
    "criterion = nn.CrossEntropyLoss()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "6427dcff",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the MNIST data\n",
    "\n",
    "(x_train, y_train), (x_test, y_test), min_, max_ = load_mnist()\n",
    "\n",
    "x_test = np.squeeze(x_test)\n",
    "x_test = np.expand_dims(x_test, axis=1)\n",
    "y_test = np.argmax(y_test, axis=1)\n",
    "\n",
    "x_train = np.squeeze(x_train)\n",
    "x_train = np.expand_dims(x_train, axis=1)\n",
    "y_train = np.argmax(y_train, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "bd22bbf7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "End of epoch 0 loss 0.4864298701286316\n",
      "End of epoch 1 loss 0.2043914496898651\n",
      "End of epoch 2 loss 0.12787242233753204\n",
      "End of epoch 3 loss 0.09274934977293015\n",
      "End of epoch 4 loss 0.0746498852968216\n"
     ]
    }
   ],
   "source": [
    "# Train the model normally\n",
    "\n",
    "def standard_train(model, opt, criterion, x, y, bsize=32, epochs=5):\n",
    "    num_of_batches = int(len(x) / bsize)\n",
    "    for epoch in range(epochs):\n",
    "        x, y = shuffle(x, y)\n",
    "        loss_list = []\n",
    "        for bnum in range(num_of_batches):\n",
    "            x_batch = np.copy(x[bnum * bsize:(bnum + 1) * bsize])\n",
    "            y_batch = np.copy(y[bnum * bsize:(bnum + 1) * bsize])\n",
    "\n",
    "            x_batch = torch.from_numpy(x_batch).float().to(device)\n",
    "            y_batch = torch.from_numpy(y_batch).type(torch.LongTensor).to(device)\n",
    "\n",
    "            # zero the parameter gradients\n",
    "            opt.zero_grad()\n",
    "            outputs = model(x_batch)\n",
    "            loss = criterion(outputs, y_batch)\n",
    "            loss_list.append(loss.data.cpu())\n",
    "            loss.backward()\n",
    "            opt.step()\n",
    "        print('End of epoch {} loss {}'.format(epoch, np.mean(loss_list)))\n",
    "    return model\n",
    "\n",
    "model = standard_train(model=model,\n",
    "                       opt=opt,\n",
    "                       criterion=criterion,\n",
    "                       x=x_train,\n",
    "                       y=y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "bb5195fa",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test acc:  97.96000000000001\n"
     ]
    }
   ],
   "source": [
    "# Lets now get the predicions for the MNIST test set and see how well our model is doing.\n",
    "with torch.no_grad():\n",
    "    test_preds = model(torch.from_numpy(x_test).float().to(device))\n",
    "\n",
    "test_preds = np.argmax(test_preds.cpu().detach().numpy(), axis=1)\n",
    "print('Test acc: ', np.mean(test_preds == y_test) * 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b14657ad",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Inferred reshape on op num 4\n"
     ]
    }
   ],
   "source": [
    "# But how robust are these predictions? \n",
    "# We can now examine this neural network's certified robustness. \n",
    "# We pass it into PyTorchIBPClassifier. We will get a print out showing which \n",
    "# neural network layers have been registered. There will also be a \n",
    "# warning to tell us that PytorchInterval currently infers a reshape when \n",
    "# a neural network goes from using convolutional to dense layers. \n",
    "# This will cover the majority of use cases, however, if not then the \n",
    "# certification layers in art.estimators.certification.interval.interval.py \n",
    "# can be used to directly build a certified model structure.\n",
    "\n",
    "interval_classifier = interval.PyTorchIBPClassifier(model=model, \n",
    "                                                    clip_values=(0, 1), \n",
    "                                                    loss=nn.CrossEntropyLoss(), \n",
    "                                                    input_shape=(1, 28, 28), \n",
    "                                                    nb_classes=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "06b46793",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test acc:  0.9796\n"
     ]
    }
   ],
   "source": [
    "# Regular accuracy on normal data\n",
    "with torch.no_grad():\n",
    "    test_preds = model(torch.from_numpy(x_test).float().to(device))\n",
    "\n",
    "test_preds = np.argmax(test_preds.cpu().detach().numpy(), axis=1)\n",
    "print('Test acc: ', np.mean(test_preds == y_test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "9fdff88c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Certified score  0.0\n"
     ]
    }
   ],
   "source": [
    "# Now to get the certified performance. \n",
    "\n",
    "# Here we will manually convert the data into its interval representation\n",
    "# upper_bounds represents the maximum any feature could take given the bound and the pixel limits\n",
    "# lower_bounds represents the minium any feature could take given the bound and the pixel limits\n",
    "upper_bounds = np.clip(np.expand_dims(x_test, axis=1) + bound, 0, 1)\n",
    "lower_bounds = np.clip(np.expand_dims(x_test, axis=1) - bound, 0, 1)\n",
    "\n",
    "# We stack this to make a representation of the data defined by its upper and lower bounds\n",
    "interval_x = np.concatenate([lower_bounds, upper_bounds], axis=1)\n",
    "\n",
    "with torch.no_grad():\n",
    "    interval_preds = interval_classifier.predict_intervals(x=interval_x,\n",
    "                                                      is_interval=True,\n",
    "                                                      batch_size=32)\n",
    "    cert_results = interval_classifier.certify(preds=interval_preds, labels=y_test)\n",
    "    print('Certified score ', np.mean(cert_results))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "1ca42f95",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "For sample 0:\n",
      "The lower bound pred: -35.76298141479492 and upper pred: 18.96987533569336\n",
      "The lower bound pred: -45.13648223876953 and upper pred: 17.85682487487793\n",
      "The lower bound pred: -29.591373443603516 and upper pred: 26.661205291748047\n",
      "The lower bound pred: -28.94809341430664 and upper pred: 26.747478485107422\n",
      "The lower bound pred: -47.43501281738281 and upper pred: 14.983299255371094\n",
      "The lower bound pred: -35.376808166503906 and upper pred: 24.173919677734375\n",
      "The lower bound pred: -56.462249755859375 and upper pred: 8.382795333862305\n",
      "The lower bound pred: -25.337738037109375 and upper pred: 31.648338317871094\n",
      "The lower bound pred: -35.181827545166016 and upper pred: 19.416484832763672\n",
      "The lower bound pred: -36.235626220703125 and upper pred: 23.697036743164062\n",
      "------------------------------------------\n",
      "For sample 1:\n",
      "The lower bound pred: -30.613969802856445 and upper pred: 18.29610252380371\n",
      "The lower bound pred: -32.045406341552734 and upper pred: 24.777889251708984\n",
      "The lower bound pred: -19.47723960876465 and upper pred: 33.95309829711914\n",
      "The lower bound pred: -29.63457489013672 and upper pred: 21.600341796875\n",
      "The lower bound pred: -51.245262145996094 and upper pred: 9.134675979614258\n",
      "The lower bound pred: -35.606964111328125 and upper pred: 18.399505615234375\n",
      "The lower bound pred: -33.94319152832031 and upper pred: 22.73233985900879\n",
      "The lower bound pred: -45.692466735839844 and upper pred: 12.092912673950195\n",
      "The lower bound pred: -26.79229736328125 and upper pred: 21.609710693359375\n",
      "The lower bound pred: -54.669654846191406 and upper pred: 7.804758071899414\n",
      "------------------------------------------\n",
      "For sample 2:\n",
      "The lower bound pred: -37.00175476074219 and upper pred: 17.81367301940918\n",
      "The lower bound pred: -30.86052894592285 and upper pred: 25.852563858032227\n",
      "The lower bound pred: -29.69383430480957 and upper pred: 26.391416549682617\n",
      "The lower bound pred: -31.832534790039062 and upper pred: 22.93268585205078\n",
      "The lower bound pred: -37.91720199584961 and upper pred: 20.41739273071289\n",
      "The lower bound pred: -35.551265716552734 and upper pred: 20.825023651123047\n",
      "The lower bound pred: -39.76785659790039 and upper pred: 21.17462158203125\n",
      "The lower bound pred: -32.9035758972168 and upper pred: 22.67422866821289\n",
      "The lower bound pred: -27.38610076904297 and upper pred: 22.835006713867188\n",
      "The lower bound pred: -40.48732376098633 and upper pred: 21.314517974853516\n",
      "------------------------------------------\n",
      "For sample 3:\n",
      "The lower bound pred: -20.57539176940918 and upper pred: 24.196565628051758\n",
      "The lower bound pred: -44.418601989746094 and upper pred: 12.289190292358398\n",
      "The lower bound pred: -27.15709686279297 and upper pred: 20.966163635253906\n",
      "The lower bound pred: -35.07857131958008 and upper pred: 12.747844696044922\n",
      "The lower bound pred: -37.69193649291992 and upper pred: 15.888842582702637\n",
      "The lower bound pred: -32.40180969238281 and upper pred: 18.97934341430664\n",
      "The lower bound pred: -34.47942352294922 and upper pred: 18.297475814819336\n",
      "The lower bound pred: -32.931392669677734 and upper pred: 19.342220306396484\n",
      "The lower bound pred: -34.2726936340332 and upper pred: 12.594486236572266\n",
      "The lower bound pred: -33.94245147705078 and upper pred: 18.43587303161621\n",
      "------------------------------------------\n",
      "For sample 4:\n",
      "The lower bound pred: -36.567726135253906 and upper pred: 15.803500175476074\n",
      "The lower bound pred: -42.30445861816406 and upper pred: 20.1270751953125\n",
      "The lower bound pred: -32.78966522216797 and upper pred: 20.093400955200195\n",
      "The lower bound pred: -36.702457427978516 and upper pred: 19.05740737915039\n",
      "The lower bound pred: -27.400218963623047 and upper pred: 29.82215118408203\n",
      "The lower bound pred: -38.59861755371094 and upper pred: 20.023099899291992\n",
      "The lower bound pred: -36.92227554321289 and upper pred: 20.892738342285156\n",
      "The lower bound pred: -32.97300720214844 and upper pred: 22.323326110839844\n",
      "The lower bound pred: -31.611032485961914 and upper pred: 21.05921745300293\n",
      "The lower bound pred: -31.900264739990234 and upper pred: 23.774620056152344\n",
      "------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "# We can also have a closer look at what the predict_intervals method returns. \n",
    "# Lets examine the last batch of predictions we obtained.\n",
    "# It is of shape (batch_size, 2, featues) \n",
    "# The axis 2 represents the lower and upper bound predictions for the datapoint.\n",
    "# We can print some of them out to manually see:\n",
    "\n",
    "for i, pred in enumerate(interval_preds[0:5]):\n",
    "    print(f'For sample {i}:')\n",
    "    for cnum, (class_pred_lower, class_pred_higher) in enumerate(zip(pred[0], pred[1])):\n",
    "        print(f'The lower bound pred: {class_pred_lower} and upper pred: {class_pred_higher}')\n",
    "    print('------------------------------------------')\n",
    "\n",
    "# Note that intervals for different classes overlap substantially.\n",
    "# For a sample to be certified correct, a class's lower bound prediction must be higher than \n",
    "# all other class's upper bound. In that way there is no perturbation in the input that could \n",
    "# cause the correct class to have a final precition smaller than an incorrect class."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ec41143",
   "metadata": {},
   "source": [
    "We can see that the certified performance is very low! Much lower than the certification you could get with using Zonotopes for example (see the deepz notebook for comparison). So why use intervals?\n",
    "\n",
    "- Computationally it is very fast: only x2 overhead compared to normal classification as each datapoint is represented by upper and lower bounds. \n",
    "- By comparison Zonotopes can grow (particularly in terms of memory) hundreds of times larger. \n",
    "- We can improve the performance by orders of magnitude if we combine it with methods like certified adversarial training.\n",
    "\n",
    "None the less, we can use the interval domain to certify for smaller regions of inputs."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26df9573",
   "metadata": {},
   "source": [
    "We can compare against the empirical results with PGD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "c693af9e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "PGD - Batches:   0%|          | 0/313 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test acc:  94.66\n"
     ]
    }
   ],
   "source": [
    "from art.attacks.evasion.projected_gradient_descent.projected_gradient_descent import ProjectedGradientDescent\n",
    "from art.estimators.classification import PyTorchClassifier\n",
    "\n",
    "interval_classifier.model.set_forward_mode(\"attack\")\n",
    "attack = ProjectedGradientDescent(\n",
    "    estimator=interval_classifier,\n",
    "    eps=bound,\n",
    "    eps_step=0.001,\n",
    "    max_iter=20,\n",
    "    num_random_init=1,\n",
    ")\n",
    "\n",
    "i_batch = attack.generate(x_test.astype('float32'), y_test.astype('float32'))\n",
    "adv_preds = interval_classifier.predict(i_batch)\n",
    "adv_preds = np.argmax(adv_preds, axis=1)\n",
    "print('Test acc: ', np.mean(adv_preds == y_test) * 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b9bea9f",
   "metadata": {},
   "source": [
    "## Training with IBP"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0d473a8",
   "metadata": {},
   "source": [
    "Now lets look at the effects of training with certification and we will be able to see orders of magnitude improvement on the certified performance!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "17c64da3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Inferred reshape on op num 4\n"
     ]
    }
   ],
   "source": [
    "from torch.optim.lr_scheduler import MultiStepLR\n",
    "from art.defences.trainer import AdversarialTrainerCertifiedIBPPyTorch\n",
    "\n",
    "# Starting from a fresh model\n",
    "\n",
    "model = MNISTModel(number_of_classes=10)\n",
    "model.to(device)\n",
    "criterion = nn.CrossEntropyLoss()\n",
    "\n",
    "\n",
    "interval_classifier = interval.PyTorchIBPClassifier(model=model, \n",
    "                                                    clip_values=(0, 1), \n",
    "                                                    loss=nn.CrossEntropyLoss(), \n",
    "                                                    input_shape=(1, 28, 28),\n",
    "                                                    optimizer=torch.optim.Adam(model.parameters(), lr=0.001),\n",
    "                                                    nb_classes=10)\n",
    "\n",
    "\n",
    "# Uncomment if you wish to train your own model, but can take around 90 min with a GPU.\n",
    "\n",
    "# scheduler = MultiStepLR(interval_model._optimizer, milestones=[50, 80], gamma=0.1)\n",
    "# trainer = AdversarialTrainerCertifiedIBPPyTorch(classifier=interval_model, bound=0.4,)\n",
    "# trainer.fit(x=x_train, y=y_train, scheduler=scheduler, batch_size=100, nb_epochs=200, limits=[0, 1])\n",
    "# torch.save(trainer._classifier.model.state_dict(), 'notebook_models/mnist/certified/certified_ibp_model.pt')\n",
    "\n",
    "# Here we will load a model which was trained using the above function.\n",
    "interval_classifier.model.load_state_dict(torch.load('notebook_models/mnist/certified/certified_ibp_model.pt', \n",
    "                                                map_location=torch.device(device)))\n",
    "interval_classifier.model.re_convert() # this is needed because the dense equivalent layer has the initial random values\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "08413cf2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The clean acc is 0.9617\n"
     ]
    }
   ],
   "source": [
    "# Get the clean accuracy using the trained interval_classifier\n",
    "\n",
    "# as we are working with normal data, and just want clean performance, we set the forward mode to concrete\n",
    "interval_classifier.model.set_forward_mode(\"concrete\")\n",
    "preds = interval_classifier.model.forward(x=x_test)\n",
    "\n",
    "clean_acc = interval_classifier.get_accuracy(preds, y_test)\n",
    "print(f'The clean acc is {clean_acc}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "8820fbc9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "PGD - Batches:   0%|          | 0/313 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Adv acc: 92.2\n"
     ]
    }
   ],
   "source": [
    "# Empirical PGD accuracy\n",
    "\n",
    "# We are interfacing the classifier with one of the ART attacks, \n",
    "# So we need to set the forward mode to attack.\n",
    "\n",
    "interval_classifier.model.set_forward_mode(\"attack\")\n",
    "attack = ProjectedGradientDescent(\n",
    "    estimator=interval_classifier,\n",
    "    eps=0.4,\n",
    "    eps_step=0.01,\n",
    "    max_iter=30,\n",
    "    num_random_init=1,\n",
    ")\n",
    "\n",
    "i_batch = attack.generate(x_test.astype('float32'), y_test.astype('float32'))\n",
    "adv_preds = interval_classifier.predict(i_batch)\n",
    "\n",
    "print(f'Adv acc: {np.mean(np.argmax(adv_preds, axis=1) == y_test) * 100}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "4b94ace4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Certified acc 0.7509\n"
     ]
    }
   ],
   "source": [
    "# Certified Accuracy\n",
    "\n",
    "# To get the certified accuracy we set the forward mode to abstract\n",
    "interval_classifier.model.set_forward_mode(\"abstract\")\n",
    "interval_preds = interval_classifier.predict_intervals(x=x_test,\n",
    "                                                  bounds=0.4,\n",
    "                                                  limits=[0.0, 1.0])\n",
    "\n",
    "samples_certified = interval_classifier.certify(preds=interval_preds, labels=y_test)\n",
    "print(f'Certified acc {np.sum(samples_certified) / len(y_test)}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "d009a548",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "For sample 0:\n",
      "The lower bound pred: -36.567726135253906 and upper pred: 15.803500175476074 adv sample pred: -13.662408828735352\n",
      "The lower bound pred: -42.30445861816406 and upper pred: 20.1270751953125 adv sample pred: -10.511754989624023\n",
      "The lower bound pred: -32.78966522216797 and upper pred: 20.093400955200195 adv sample pred: -4.524026393890381\n",
      "The lower bound pred: -36.702457427978516 and upper pred: 19.05740737915039 adv sample pred: -0.1533353626728058\n",
      "The lower bound pred: -27.400218963623047 and upper pred: 29.82215118408203 adv sample pred: -9.642894744873047\n",
      "The lower bound pred: -38.59861755371094 and upper pred: 20.023099899291992 adv sample pred: -5.101047515869141\n",
      "The lower bound pred: -36.92227554321289 and upper pred: 20.892738342285156 adv sample pred: -20.015060424804688\n",
      "The lower bound pred: -32.97300720214844 and upper pred: 22.323326110839844 adv sample pred: 6.100724697113037\n",
      "The lower bound pred: -31.611032485961914 and upper pred: 21.05921745300293 adv sample pred: -11.406320571899414\n",
      "The lower bound pred: -31.900264739990234 and upper pred: 23.774620056152344 adv sample pred: -3.169931411743164\n",
      "------------------------------------------\n",
      "For sample 1:\n",
      "The lower bound pred: -36.567726135253906 and upper pred: 15.803500175476074 adv sample pred: -4.840512275695801\n",
      "The lower bound pred: -42.30445861816406 and upper pred: 20.1270751953125 adv sample pred: -0.9820030927658081\n",
      "The lower bound pred: -32.78966522216797 and upper pred: 20.093400955200195 adv sample pred: 4.475952625274658\n",
      "The lower bound pred: -36.702457427978516 and upper pred: 19.05740737915039 adv sample pred: -2.523003578186035\n",
      "The lower bound pred: -27.400218963623047 and upper pred: 29.82215118408203 adv sample pred: -6.6240386962890625\n",
      "The lower bound pred: -38.59861755371094 and upper pred: 20.023099899291992 adv sample pred: 0.5403081178665161\n",
      "The lower bound pred: -36.92227554321289 and upper pred: 20.892738342285156 adv sample pred: 2.613499879837036\n",
      "The lower bound pred: -32.97300720214844 and upper pred: 22.323326110839844 adv sample pred: -10.114892959594727\n",
      "The lower bound pred: -31.611032485961914 and upper pred: 21.05921745300293 adv sample pred: -2.541588306427002\n",
      "The lower bound pred: -31.900264739990234 and upper pred: 23.774620056152344 adv sample pred: -7.538417816162109\n",
      "------------------------------------------\n",
      "For sample 2:\n",
      "The lower bound pred: -36.567726135253906 and upper pred: 15.803500175476074 adv sample pred: -16.553003311157227\n",
      "The lower bound pred: -42.30445861816406 and upper pred: 20.1270751953125 adv sample pred: 3.684479236602783\n",
      "The lower bound pred: -32.78966522216797 and upper pred: 20.093400955200195 adv sample pred: -2.3804402351379395\n",
      "The lower bound pred: -36.702457427978516 and upper pred: 19.05740737915039 adv sample pred: -2.618450403213501\n",
      "The lower bound pred: -27.400218963623047 and upper pred: 29.82215118408203 adv sample pred: -6.235487461090088\n",
      "The lower bound pred: -38.59861755371094 and upper pred: 20.023099899291992 adv sample pred: -3.8312270641326904\n",
      "The lower bound pred: -36.92227554321289 and upper pred: 20.892738342285156 adv sample pred: -6.061868667602539\n",
      "The lower bound pred: -32.97300720214844 and upper pred: 22.323326110839844 adv sample pred: -6.5338287353515625\n",
      "The lower bound pred: -31.611032485961914 and upper pred: 21.05921745300293 adv sample pred: -3.47023344039917\n",
      "The lower bound pred: -31.900264739990234 and upper pred: 23.774620056152344 adv sample pred: -9.340883255004883\n",
      "------------------------------------------\n",
      "For sample 3:\n",
      "The lower bound pred: -36.567726135253906 and upper pred: 15.803500175476074 adv sample pred: 3.856935977935791\n",
      "The lower bound pred: -42.30445861816406 and upper pred: 20.1270751953125 adv sample pred: -8.496208190917969\n",
      "The lower bound pred: -32.78966522216797 and upper pred: 20.093400955200195 adv sample pred: -1.1915643215179443\n",
      "The lower bound pred: -36.702457427978516 and upper pred: 19.05740737915039 adv sample pred: -1.7413063049316406\n",
      "The lower bound pred: -27.400218963623047 and upper pred: 29.82215118408203 adv sample pred: -1.4715278148651123\n",
      "The lower bound pred: -38.59861755371094 and upper pred: 20.023099899291992 adv sample pred: -0.8217955827713013\n",
      "The lower bound pred: -36.92227554321289 and upper pred: 20.892738342285156 adv sample pred: 1.4374966621398926\n",
      "The lower bound pred: -32.97300720214844 and upper pred: 22.323326110839844 adv sample pred: -4.385831832885742\n",
      "The lower bound pred: -31.611032485961914 and upper pred: 21.05921745300293 adv sample pred: -1.8728493452072144\n",
      "The lower bound pred: -31.900264739990234 and upper pred: 23.774620056152344 adv sample pred: -1.5020307302474976\n",
      "------------------------------------------\n",
      "For sample 4:\n",
      "The lower bound pred: -36.567726135253906 and upper pred: 15.803500175476074 adv sample pred: -6.092507839202881\n",
      "The lower bound pred: -42.30445861816406 and upper pred: 20.1270751953125 adv sample pred: -8.040410995483398\n",
      "The lower bound pred: -32.78966522216797 and upper pred: 20.093400955200195 adv sample pred: -6.087466239929199\n",
      "The lower bound pred: -36.702457427978516 and upper pred: 19.05740737915039 adv sample pred: -6.674464702606201\n",
      "The lower bound pred: -27.400218963623047 and upper pred: 29.82215118408203 adv sample pred: 4.123129367828369\n",
      "The lower bound pred: -38.59861755371094 and upper pred: 20.023099899291992 adv sample pred: -5.928511142730713\n",
      "The lower bound pred: -36.92227554321289 and upper pred: 20.892738342285156 adv sample pred: -6.545170783996582\n",
      "The lower bound pred: -32.97300720214844 and upper pred: 22.323326110839844 adv sample pred: -4.635064125061035\n",
      "The lower bound pred: -31.611032485961914 and upper pred: 21.05921745300293 adv sample pred: -2.049896717071533\n",
      "The lower bound pred: -31.900264739990234 and upper pred: 23.774620056152344 adv sample pred: -3.035524368286133\n",
      "------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "# We can also convince ourselves that the upper and lower bounds contain the \n",
    "# predictions for the adversarial data by having a look at some of the predictions. \n",
    "# Lets do so for the first 10 samples\n",
    "\n",
    "for i, (adv_pred, cert_pred) in enumerate(zip(adv_preds[0:5], interval_preds[0:5])):\n",
    "    print(f'For sample {i}:')\n",
    "    for cnum, (adv_p, class_pred_lower, class_pred_higher) in enumerate(zip(adv_pred, pred[0], pred[1])):\n",
    "        print(f'The lower bound pred: {class_pred_lower} and upper pred: {class_pred_higher} adv sample pred: {adv_p}')\n",
    "    print('------------------------------------------')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69dbd57c",
   "metadata": {},
   "source": [
    "## Analysis\n",
    "\n",
    "We can see that the certified accuracy is very high, with the trained model reaching 75%\n",
    "\n",
    "Also, note the difference between the empirical adversarial accuracy and the certified accuray. \n",
    "The empirical adversarial accuracy is higher (92% in this case) and the gap with the certified accuracy shows that\n",
    "there could be several adversarial examples that exist which PGD was not able to find. \n",
    "\n",
    "Note! That IBP overapproxaimates the impact of aritmetic operations on the inputs, and will give us a more pessimistic view on robustness than the model is likely capable of, but it gives us two important peicies of information: \n",
    "\n",
    "+ The certified performance with IBP of the model is a lower bound on the performnace. Thus there may be some datapoints that IBP says *could* be misclassified, but in fact cannot be.\n",
    "\n",
    "+ If IBP says a datapoint cannot be misclassified then we can be absolutly certain that the model is robust to the given input and allowable perturbation.\n",
    "\n",
    "The true performance of the model therefore will lie somewhere between the certified and empircial robustness."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}