notebooks/certification_interval_domain.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f9bf98ec",
   "metadata": {},
   "source": [
    "# Certification of Robustness using Intervals"
   ]
  },
  {
   "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)"
   ]
  },
  {
   "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 one using different abstract domains\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=32,\n",
    "                                      kernel_size=(3, 3),\n",
    "                                      stride=(1, 1))\n",
    "\n",
    "        self.conv_2 = torch.nn.Conv2d(in_channels=32,\n",
    "                                      out_channels=32,\n",
    "                                      kernel_size=(4, 4),\n",
    "                                      stride=(2, 2))\n",
    "\n",
    "        self.conv_3 = torch.nn.Conv2d(in_channels=32,\n",
    "                                      out_channels=64,\n",
    "                                      kernel_size=(3, 3),\n",
    "                                      stride=(1, 1))\n",
    "\n",
    "        self.conv_4 = torch.nn.Conv2d(in_channels=64,\n",
    "                                      out_channels=64,\n",
    "                                      kernel_size=(4, 4),\n",
    "                                      stride=(2, 2))\n",
    "\n",
    "        self.fc1 = torch.nn.Linear(in_features=1024, out_features=512)\n",
    "        self.fc2 = torch.nn.Linear(in_features=512, out_features=512)\n",
    "        self.fc_out = torch.nn.Linear(in_features=512, out_features=number_of_classes)\n",
    "\n",
    "        self.relu = torch.nn.ReLU()\n",
    "\n",
    "    def forward(self, x: \"torch.Tensor\") -> \"torch.Tensor\":\n",
    "\n",
    "        x = self.relu(self.conv_1(x))\n",
    "        x = self.relu(self.conv_2(x))\n",
    "        x = self.relu(self.conv_3(x))\n",
    "        x = self.relu(self.conv_4(x))\n",
    "\n",
    "        x = torch.flatten(x, 1)\n",
    "\n",
    "        x = self.relu(self.fc1(x))\n",
    "        x = self.relu(self.fc2(x))\n",
    "\n",
    "        return self.fc_out(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b239e7da",
   "metadata": {},
   "outputs": [],
   "source": [
    "model = MNISTModel(number_of_classes=10)\n",
    "model.to(device)\n",
    "opt = optim.Adam(model.parameters(), lr=1e-4)\n",
    "criterion = nn.CrossEntropyLoss()\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": 4,
   "id": "bd22bbf7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "End of epoch 0 loss 0.2816831171512604\n",
      "End of epoch 1 loss 0.07683608680963516\n",
      "End of epoch 2 loss 0.051124896854162216\n",
      "End of epoch 3 loss 0.03650381416082382\n",
      "End of epoch 4 loss 0.027725202962756157\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": 5,
   "id": "bb5195fa",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test acc:  98.72\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": 6,
   "id": "b14657ad",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "registered <class 'torch.nn.modules.conv.Conv2d'>\n",
      "registered <class 'torch.nn.modules.activation.ReLU'>\n",
      "registered <class 'torch.nn.modules.conv.Conv2d'>\n",
      "registered <class 'torch.nn.modules.activation.ReLU'>\n",
      "registered <class 'torch.nn.modules.conv.Conv2d'>\n",
      "registered <class 'torch.nn.modules.activation.ReLU'>\n",
      "registered <class 'torch.nn.modules.conv.Conv2d'>\n",
      "registered <class 'torch.nn.modules.activation.ReLU'>\n",
      "registered <class 'torch.nn.modules.linear.Linear'>\n",
      "registered <class 'torch.nn.modules.activation.ReLU'>\n",
      "registered <class 'torch.nn.modules.linear.Linear'>\n",
      "registered <class 'torch.nn.modules.activation.ReLU'>\n",
      "registered <class 'torch.nn.modules.linear.Linear'>\n",
      "Inferred reshape on op num 8\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_model = 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": 7,
   "id": "06b46793",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test acc:  98.72\n",
      "Certified score  0.0\n"
     ]
    }
   ],
   "source": [
    "bound = 0.01\n",
    "num_certified = 0\n",
    "num_correct = 0\n",
    "\n",
    "# Use the test data to check its certified robustness.\n",
    "original_x = np.copy(x_test)\n",
    "\n",
    "# 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))\n",
    "\n",
    "# Here we will manually convert the data into its interval representation\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",
    "interval_x = np.concatenate([lower_bounds, upper_bounds], axis=1)\n",
    "\n",
    "with torch.no_grad():\n",
    "    interval_preds = interval_model.predict_intervals(x=interval_x,\n",
    "                                                      is_interval=True,\n",
    "                                                      batch_size=32)\n",
    "    cert_results = interval_model.certify(preds=interval_preds, labels=y_test)\n",
    "    print('Certified score ', np.mean(cert_results))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ec41143",
   "metadata": {},
   "source": [
    "We can see that this 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 (or also for lower dimensional inputs). Let's now use it to certify against a pixel brightening attack"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "20bd7f5b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Certified score  0.4766\n"
     ]
    }
   ],
   "source": [
    "bound = 0.005\n",
    "int_x = np.expand_dims(x_test, axis=1)\n",
    "\n",
    "# pixels of a certain brighness can be raised to the maximum value.\n",
    "upper_bounds = np.where(int_x > 1 - bound, 1, int_x)\n",
    "interval_x = np.concatenate([int_x, upper_bounds], axis=1)\n",
    "\n",
    "with torch.no_grad():\n",
    "    interval_preds = interval_model.predict_intervals(x=interval_x,\n",
    "                                                      is_interval=True,\n",
    "                                                      batch_size=32)\n",
    "    cert_results = interval_model.certify(preds=interval_preds, labels=y_test)\n",
    "    print('Certified score ', np.mean(cert_results))"
   ]
  }
 ],
 "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
}