notebooks/certified_adversarial_training.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7f6145ce-083e-4213-af7a-9f51804b2d4a",
   "metadata": {},
   "source": [
    "# Certified Training using Zonotopes with DeepZ"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6248072b",
   "metadata": {},
   "source": [
    "In this notebook we will look at how different training strategies impact the certified performance, and see that by using certified training very high levels of certification can be obtained. \n",
    "\n",
    "The drawback is that using certified training is computatianally expensive, even when compared to PGD adversarial training.\n",
    "\n",
    "The zonotopes abstraction used here is defined by:\n",
    "\\begin{equation}\n",
    "    \\hat{x} = \\eta_0 + \\sum_{i=1}^{i=N} \\eta_i \\epsilon_i \n",
    "\\end{equation}\n",
    "\n",
    "where $\\eta_0$ is the central vector, $\\epsilon_i$ are noise symbols, $\\eta_i$ are coefficients representing deviations around $\\eta_0$.\n",
    "\n",
    "By pushing zonotopes through a neural network we can obatin an overapproximation of the worst case loss and optimise the weights for certifiable robustness."
   ]
  },
  {
   "attachments": {
    "ART_Cert_Training.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwEAAAD0CAYAAADDqCoiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAOxAAADsQBlSsOGwAAABl0RVh0U29mdHdhcmUAd3d3Lmlua3NjYXBlLm9yZ5vuPBoAACAASURBVHic7N13eFTF+sDx7+ymQcAQkA4qRbh4FUQFC1h+CIpioSWgoKLSIkUQUeSKIiKioigEMRTxoihSlHJtNEEUbKiIXbwignDpEkLa7r6/P+ZssoRNg4RNeT/Ps0/YU99zNmxm5rwzA0qVcTVr1nwZkFC/3G53JnBlcV+vUkoppVR+XKEOQKnilpGRccWmTZsQkZC+hg0bJsDFob4fSil1SokYROKc140ncZzuzjFuLsLo1MkQqRfw2V4Q6nBU4YSFOgCllFJKnSCRBsA0IDzHmhTgIPA18CXwKcb4TnF0fm7gTSAVyABiT/A4bwJp2KerlYomtCIiMgzoFGRNBvAr8AOwGGP2n9K4it+lwFzAAPOAu0MbjioMrQQopZRSpVcsNs2wYi7re2ILoumIPA4kYoycquACCDbGzJM8TkXg6MmHU+SuANrnsT4FeA6R/hjz+imK6VTJAE7DVgRUKaLpQEoppVTplhnwM815pQPJzrLTgOrABGAVIjmfGpwKLmxh0XMSxxDnGN4iiah4+Mj+DNKwFZbD2MpLNDATkfNDF55S2fRJgFJKKQUtgWqhDgKbOvLHCe6bBpwZ8L4W0Bp4wFleCbgEeAG4J9ejiLiAGODvQqUQiVQBjmDMsQV9YzyIVHXe5X48kQggEmOSc9miun/LfOKoBPgwpnBPDEQigYg8zl8QPuBp4HnnfTjQAngNqAFEAUOBu/KIoxLgxZjUk4ij4ERigJTjPrfg2xogBmMOFfIcEUAFjPn7BOKLxn6ep+Z+lCNaCVBKKVXetahUqdKG8847LyPUgXz99dcmLS0thvwKusEJxhwMeH8Q+BGR14A5QFdsa3QfRCZizPbsPcUF3AGMBJpgW+zDENkMjMaYD4KfUToD/wLOx7bQhyFyCHgbm3q0GRE38DG2ALyDwFHSbJ+GSUA7oDLgdWL5C1gCPBsQ5zon/j3YXPTAOOoDTwI3A5HOsiPA68Cjx+Tii9QA1jrbbcH2NRiPrSh5nf2eByZgzIk8dUjN8TmsROQp7JOYSKDZcXuInAVMBG4AIgCDyGFs5WHsMccTaQ28ii3DLcaYB3IcKw57LwzwFMbMcJY/BvTGfk5jgC5AZ2cvFyJfAAMxZkuQ+BoDicDVZH9GK7D38Xi2YtEFuA24CPu5eRAB2Ij9ndqYY58+wMPOu+eBc4Fbsb83BpE/sf8vjgJtMOZwkPO+gL2HXux9K2upV0qpwoiNjd26adMmCbURI0ZkYFvklFIlS7tWrVodDPV3hIgItiW54LnVIhcgcsgZhiz3VlaRSEQOOtulIjIix7qViCTnMrzZEURG5zieC5E38tjHh8giZ9sw5704FQT/Mc5G5G9EMnM5xlGnQOvf3ussT8nlHmQEOUYqIrudyoZ/+3rOsQURj3N9wa755UJ8Dm85+2UiMirI+scDrvOtHOsuzuM+pCKyw6nk+Ldvj8hhZ/2yIOcahEias35swPLXnWXJzrlyns/nHPfsHMdr4Sz35Nje4xzjb+f9ywH79Ay4x8F+N1IQic9xnlEBx/QFXIP/d8H/PgWRhCDXXcW5X/7PoXbeH5rSPgFKKaVUWWdMOrZlXbAtqx0D1iYCbbCttSnAJ8Bs7MhCR5zl/0IksPX9MeBGbIpRCvAjMANYDPzJ8Xn7wZ5sPOUc2wvsBf4NTAHewT7FiMj3ukQqAu9i05cA/naOs8CJKwKbRrQcm8ri509LcmNTdrYBb2FH8fFfcw9Ejm+1z19jp6DeHpHOiEwAhmNb7o9gnzz44z8N+A+234YPOAS8AizCtnhHYtO6luQ4R35PKPJK44rAXncG8BGwHNt/xIftuzA5IL5wYCn2c/YCB4Ak7EhAac5xKuRynjDsiFAbsZ/JfGA39nOpCMx2rj+QYMum/grx99jP9yj2niQ7+94X5Hy3OTEKsAZjduVxDxSaDqSUUkqVF18CPbCFtjMAfxpKL2dZCtARYz7O2sOOKDQMWyh+FOiISCy2EFbR2ecRjHnumDOJtAca5RNPK2wh0gM0y5GyE4ZN7dmWzzFudWIDW0D8J8bsdo7RDNgUcL0dsCksgXzAHIwZ6OxjgPexI/24nH1+zCeGQG4gznn52XQWWwF4F1gYsO4OsgvRh4BzMWavE0tz4FNn/dmIXIExHxUiltxEYAvV52HMf51znY9N2YoG/i9g25uAqk78R4F/BMTXElvAjwxyjr3YlKOZGHMga6n9XFdgU8IE+xm/mmNfF/Z3oivGvO/sF4Xts9MemzZWG5FWGPNFwH73OfEnA8f+Pqqg9EmAUkopVT4cJnt0nsrOzy5kpx8tPaYCYD2OP0cdLncKydcEHGfHcRUAAGNWYUxSPvH4U4MygZFOhcS/vwdjFmPMpnyO0QPbSn0UeDyrAmCP8SO2wO11tukSZH8vth+Efx/BFkpTsIX3JvmcP5iIgFcF56e/dbsm2ZUWf/z+JzCPZhWwbSzfYlvpfc42RTVJmg94MasCYM/1DbaSAuAOaKHvhv1dSQOeyxHf18CHQc9gzGqMeQrbD+BaRO5C5EFgBLaCcITsjurBfJVVAbDHS8OYndi+AkfJ7mBt2adUpzvvUoGVed4BBWglQCmllCovorEt1WALnWDTgKKwhfpbj8vftkON+tNy3NgCYXOyKxHrTiKe+dgCXSVsge4HRPYg8orzJKEg/uH8zAC+CrJ+hXMOAwSb0fZAkNGA/kd2us2ZFI4X21l2rPN6FNspO5XsQu+8gO39+feZecSfii2vXVjIWHKTAnwTZLm/gJ+BTUEC+1mDrQR8EmSf1QSb+8HOEj0eey/nA1OxHaMnAvFkzyuQW97+xlyWz8LeCzfQDRH/7+G92CdTqcDUEE6MV6poJUAppZQqH84lO/XE32LuHxa1IOnBLmxFojq2AOcBfjuJeJ7G5sOnkJ1bXh2b2/0WIh86Of958a93A8FywHeTXaDPmX8OttCYU2AB8kTKSesw5inn9STG3IUt/KdhU2euRsSfKuWPP5zg8e8iu5AdE2T9ifARfMI1/3X78/Ihe2ZmF3ZUppx2Y68rJ38aWRS2Eunl2Dki0p2fwX7vMoGdQSM3Zh+2z4jPOWYv7EhENzsxGmxFQRWA9glQSimlyjo7pGNXbCEpg+x0Cf8wiynYFuz1+RwpGds5VLBliMK2lGczJhPb+bYFtnX4ZmzLeAb2ScPFwDPAoDyO4i+AerEViF9zrK9B9tOPkxn//+QY8z0iR7EVnUzgPGwFKg1b0PZg4/8zx541yC6r+T+rwE7WwSopFSi68p2/suAjO90mUHVyduC2cwLcR3Y/kwnYtKz/YcxhRG4BXiJ4X4KCeA64Fnvf7sNeq79D8NpjUsJUnrQSoJRSSpV9D5KdwuPBjrQCtuNsR2wLf5qTx503kS3YnO7KQNuTjsyYzcBm7AhElbGpQf54O5N3JeAXoA62Jb01sCHH+suxhVEBvj3pWE+UrYT5+wIYbAs52IrA6diKSiuOTwm6Avu0wIe9R2BTbPz9OGoGOduF2PtRFL7HzmsQ4Rx3bY71bTm+MN8QW9GpAPyFMRNyrL+I7CcMhWfMBkT2OMeoA4xDOwSfEE0HUkoppcoiOz5/K0QWYif0qoRteV6FMT84W71FdprG9Yi0C3KcSETuQOQhZ8kHZDciNkKkd5B9mjstvnnF9woilx+zzObnzw44frAUnkCLsK3N0cBDTmqI//j1saPvhGELiDmH2Tw17Mg2E8m+z26yC/T+oUArAY9gZwv279cA23HYhb1G/5wAf5Cd1tUMkWoB+/yD4B2gT9TbZA/L+UBADj7OfALXB9knmeynLzHOiED+fRoC/Tn58ucksjtv+ytUR7F9FFQB6ZMApZRSqmyojIg/FSISWzBOxRbg/GO2/4XNubeM+Q6Rd7Bj/lcA/oPIl9ihKQX4J7Y1PRJY5eyzD5GXgAHOOZIQ6eesr4btbHwBtoLxRh7x9gTinJjXY1v1KwO3Y1u6fcBn+VzzXLJbgk8DfnImraoA3OnELcA+4L18jlUUXNgnGsOc9xHYa0p3YswANjsjF4Gt8Phnya3qxD8HWym4y4nfh51p2RZwjUlGZBt25CI38B12YrZ62CFN/Tn4/oL4yXgL27peGXt/f0BklnMt/Z347OzSfsbsRMTfvyAa+Bw7QVp155qij9un8F7DVgT8fUlSgSnaIbhwtBKglFJKlW7+VlXD8ekh4dgWUw92VJphGJOSY5s+2D4CLbAFtMvJTvPxp514OHaCqgexee2XOvtc4ewXOCFXYGtvbrMgV8Smj/hTSPxpLB5si3JeqUD+AnFXbGfRCthRbfxPLIxzzMPATRgTGH9+LdEn2lLtwhbgc6a7+J9GbCNwqE9jDiLSHfuUogJQFxgdcH5//J2d4Uv9Bjn7RGOvebCz3BOwX26VgIJfuzFpzv1dgf2s6mHH/w+ML4Pjr/dhbF+AaKAl9nfFHbBPGsGf8hRstmz7uc/HVhjDnP1mF2hflUXTgZRSSqnSKx2bn78v4LUX+B2bcvIqdvjEuhjTL0gFAIw5ii3E3+/sl+kcMxVb8P8vMB47xrt/n0xsX4KRwHYnjhRsSkYGtoA6ydlasE8g9nHs5F/3YXPMU53905yfqdjUl5YY83PA9jucY/yRI/512P4AqwPux1HneEuAFhjzfcAeHuyMxPuca8sp1dl/H8d31M2NP7acrz3YJxxvYZ9MtHRGuAmMfxVwGXa41XTn3ClOHIuB5hjzS5B94p34Up1rzsCmaj2ArXDsc67Tz/8Z+I+f0+/O+mQCh/00ZgP29+ObgPgysf1JxjjH2kfg6EbGTMHOkryX7M/Ugx36czDZ93dHwPkPYTudHyR7Dom8fO9cs3+G4P8VYB+lVHkSGxu7ddOmTRJqI0aMyMB+OSulSpZ2rVq1Ohjq7wgREWzqR8FaQ4uLSCVEGiNSB5GCpZSIxDr7BOuoWpD9qzj7x57Q/tnHiUSkISJnIFJUnWNPHZEoJ/76BY5fpBYijRCpkP/GJyn7c8qvr0bgPjHONeU33Gth4nAhstOZz+Jw0L4sKl+aDqSUUkqpbMYcAbYWcp+DHNvqXNhzHqJgrb/5HSed4K37pYMxaRQ2/lM5JOaJfE7G/A38XSTnF7kSm5Z0M9npRLswZk2RHL+c0UqAUkoppZQqDf6DTf+JwHZKTsFOSqZOgFYClFJKKaVUaeElu7/KBIw5FaM+lUnFXQmoClTBdgwJ3Ux9SimllFKqtOuNfQqwF9iCMftDHE+pdiKVgHBsj2+DHUd4XR7bPoAdRqwfMOsEzqWsjtih2z4BPg5xLEoppZRSp54xS0MdQllyIkOEdgKews5+N7lowzkhK7DTblcNdSD5WIeN80Smyu6Cvd/tizSiU+M97HXXCHUgSimllFLKOpFKwF3Oz7+xE0CcX3ThnJD62ElGimJmvOLkj7O8zc1Qj9Lx+SillFJKlRuFLZDWxKam/AU85izrU5QBKaWUUkoppYpXYfsE3I7tEzAPeA2bFtQbm/ufkc++ZzjbXYod1mkT8DR2xrecqmKHfGrl7JcJ/A/43Dn3T0AT7OyFtZx9JmFnB/QbhJ2d7mbgeuB14Ednn8uBakBfYL1zvuuAa7Ct1jWxM9l9CiSS95i9rYC7gQuAGGA38DUw39n/HOxsjac727/AsfdqQB7Hzk9HbKrQYudc9zvXEIOdofAFINjYuZOBKCABuBHoDzTCjv27HHge2/M+UBfnfHOxfRNyGoh9MvQMdnzpxtiZJOs465/GzhDoNxQ7i6BSSimllCrhfsCOz3qu836J8757LttPdNY/BezHTgf9EXbacMEWNDvk2KcWdkpwcfbZCHyBnXpbgLHOdm2d43mc5Qed9/5XpLPd4876J7FPMAQ7tfVBbAUBZ51gKxE/A185xxDsqEZX5nJ9E5xtxIn5c+wU2OIcA6CdcyxvLnEWRFKOa/cb5SyfiK0YeYFfsfdNnPddgxzvIPa+PeRs9zO2MnTEef85UDmXa707lxgXO+svd95fwrGfj386cP+r6GYOzIfOGKyUyofOGKyUKncKkw50KdAMW7j9zlk21/l5Zz77jgBWYfPDrwDOAgZjC+rzgMBpwodhW/+nYSsEl2Jb22s4P1c7232MbcH/1XnfxHnvf+VsZR4BfAbUdl7VAo71JdAZ22m3KbZVvyZ2FKRo4BWOz2nvjy1E78VWZM4EWpOdAz/d2W6NE88fzvszc8RZFO4DvgWqA2djnzqMwX6+kwj+B8UFPAJ0w17z5dj7vhZ7nyeeZEyfYq/vR+d9M4697qO57KeUUkoppYpZYSoB/oL+3IBl/8G2Ol9LdtpHMCnYQnNgwW8a8AG24NonYHlj5+e/sWlAgb7EtlifiL+x6Uz+6bV92JZvsK3YS7Gt1n6Z2LSYJdhKy6UB68LJbpXvha3gBPodmHmCcZ6IvdgO2/4nC4Jtuf8daADUDbKPAV4F3gpYdsA5jgfb4l+tmOJVSimlCk/EhUj1UIehVFlQ0EpARaAHtmD8RsDyDGzuuxtbwM7NQmwhPCf/3AGBQ1/68+9HYlvsi8oi8p+wLBLbin+5E1N74LCzrlnAdhc6sW0FVhZhjCfqP2RXaPx82KcDYJ8+BBNs7obfsU9IIslO7VFKKaVCS8QFvAy8GepQlCoLCloJ6AacBryPzc0PFJgSlFse42+5LN/q/GwYsOwFbM5+HDa//ktsSstVnNwwk7/nsa4KtuX+gBPrR9jC/UrgDmebmIDt/U8rfj6JeIrStlyWH3R+xuayPr/PpcGJBqSUUkoVsUuxffT6InITIvVCHZBSpVlBRwfyzw1wJrAgyHovNif/MoKPHJNbC7y/lT0qYNlO4J/AEOAm7IgzF2Jz+r8GbsV2gi2s3GIwwDJsq/dG7KhH27EF6DRsJWAIx1aY/PEGe7oRCp5clks+++V2T/zLo3JZr5RSSp1axnyCSDNs2WAucAu2sVApdQIKUgloSPboOM2dV27uJHglIFhOOthOtGCH4wx0EBjnvGKBq7GddC/DFtIvyjfqgrsYWwH4CnudOfsh9Amyjz/es4owjlCoS/AnJP7PK/Bz8Vc0cvudKekzNiullCrtjElDJBxb3vgGkVXALRizN8SRKVXqFCQdqA+2tfzfHDu6S+CrGTYHPR47mk5ObXM5tj/n/Os8zn8Qm89/NXaYyQux6Tt+/rkBIvK+jFw1cn5u4PgKANgRf3L6DNvK3hybJlUQJxtncQiW8+8G2jj/Dvxcdjo/g3UAj8A+vQmmJF63Ukqp0sqYecDb2EE5dgGHEakZ2qCUKn3yqwS4yM6J/ze2QB7s9ROwDju2fLA5A9oC/5djWSx2mFCww4T65ZaHbpx4vM7Lz/8osNFxexSMf//zOb5PQzzBKwG7sJNqVSL3ITjDcznPicZZHEZyfMpPL+zTn5+xE7r5+YdijeP4vhn3YUd5CqYkXrdSSqnSzBh/um5fbBnilZDGo1QplF860NXYseN3YAv5eZmHLejfha0wBPoBW2sfj52Iqj52jP3a2OE5A4/9BrbVeAG24LnX2f4ebKv7Ao7NZV8P3IBNE1pC9jCkD5N7rnygL51ztCV7JuR07AzCQ4HNQIsg+w3GphL1w3YUnoOtHNTGtqSfg50TITDODthRDZYFxDmqADEWBx/281+NnUF4L/bzHumsG86xfQo+wnYkbgaswF6vGzsZWjfsfACBIyj5rcfOwfBv7OeT4ix/hPxnmVZKKaWCM+ZjRK7CTmB6NSI9gI8xZmfeOyqlCuINbEHwmQJsWwU7A7CP7NFz/DMG9wemkD1rrv/1JsfPHPsStqAoOV4ebE2/Uo7tI4FEbAE8cPucMwYn5BH75UH2TwXud15C8JleG2JHEMot1kAVsBOI/S/HtgWR34zBI3PZ72Vn/Q05lvtnDD4b2JIjnr+B3rkc73yyZ3v2v3Zjh1LNOWOwXwR2xKe/cuwXLG2sWOiMwUqpfOiMwaWZSAQidyGSisi1oQ5HqbIiBpu2U9B87tNybF/Bee8vkNcHumJTTprkcZwI7Ky1NwO3YVvUCzI5SJRzvsAhMXPGkNe+bbAF4OvJHhLUf8y8Rso5C+iCnSuhA3a24cLGmZeKzrYVcjlObrFFO+tzpib5KwFgU6wuxd7nG7ApXXmJxD4xuA24huz7Wsk5V15PlyLJvu5T9kdOKwFKqXxoJaC0E+ntDBv6D0Q+QkQnu1QqH/mlAxV2CMzDOd6nOi+/P51XfjKALwp5brCdUNNyLMsZQ177fsLxoxsFO2ZO28h9rP7czpXfMQMd5djZlgt6nBSy029y48MOjbqxgLGkY1OIcso5WVlu+6YX8DxKKaVUwRjzGiLR2P5sy4EURGphzO4QR6ZUiVXQeQKUUkoppUouY1IQ6Q58A7yFbXTqEtqglCq5CjpjsFJKKaVUyWbMp9i+aXWBAYj0QSS3uYqUKtf0SUD59AI6br9SSqmyyJiViFyIHVXwSWw/Px0xSKkctBJQPo0NdQBKKaVUsTHGg8gB7AAjfyOyEbgOYw6FODKlSgytBCillFKq7LGdhU/Dzjn0byBTOwsrlU0rAUoppZQqm4w5jMj1wE/AO9g5gW4JbVBKlQzaMVgppZRSZZcxm4ArsXPqJCCSgEidEEelVMhpJUAppZRSZZsx72InBB0EPAXUR0QnZVPlmlYClFJKKVX2GeMFtgMdseWfzxGpFNqglAodrQQopZRSqnww5lVs/4B3gWVAGCK1QxuUUqGhHYOVUkopVX4YcwCRK7FPBVYA3wN3hjYopU49fRKglFJKqfLFmG+xnYVTgCGIjNAnAqq80UqAUkoppcofY5YCVwMPAw8BNRDRcpEqN/SXXSmllFLlkzE+4DugPVAN+AKRCqENSqlTQysBSimllCq/jHkNO4nYUmAuUElTg1R5oB2DlVJKKVW+GfM/RFoD+4DVwMfAPaENSqnipU8ClFJKKaWM+RG4AtgBjEDkX/pEQJVlWglQSimllAIwZjHQCXgGGAzEIOIObVBKFQ+tBCillFJK+RkjwEagHdAY+BKRiNAGpVTR00qAUkoppVQgY+YBh4E3gUSgOiK1QhuUUkVLKwFKKaWUUjkZsxNoAbwPfAiMDG1AShUtrQQopZRSSgVjzFagLbAZ+Bci4/WJgCortBKglFJKKZUbY97EmDjgJeA2oAIiOsS6KvW0EqCUUkoplb8VwFVAa+zMwloRUKWaVgKUUkoppfJjzOuAAHOAiUB9RGqGNiilTpxWApRSSimlCsKYbcA5wBfAWuxcAkqVSloJUEoppZQqKFsRuBRYA4xHZJI+EVClkVYClFJKKaUKw84jcBcwD7gJcOuEYqq00UqAUo7bbrsNYwxffvllqENRSpVQXbp0wRjDTz/9FOpQVKjZmYUXYTsLXw9sRETLVarU0F9WpRxhYXagh8aNG4c4EqVUSRUWFoYxhkaNGoU6FFUSGDMfqARMAR4FmmhqkCottBKglKN58+YA/PrrryGORClVUjVv3hwRYevWraEORZUUxvwCNAF+w3YWvjuk8ShVQFoJUMrRu3dvqlevzvPPP4+IhDocpVQJdOedd1KlShWef/75UIeiShJjdgCtsOlBkxGZqk8EVEmnE10o5YiMjGTq1KkMGDCANm3aMGTIEJo0aUJsbGzQ7WvXrk2FChVOcZRKqVCqXLkyL7zwAvfccw8///wz99xzD2effTYxMTFBt69bty6RkZGnOEoVEsbMReRV4B2gHjAOkSiMSQtxZKoEktk3VebPC1LM2LG+UMWglQClHPPnz2fAgAEAbNy4kY0bN+a5/erVq2nXrt2pCE0pVULMnj2bESNGALBu3TrWrVuX5/afffYZrVu3PhWhqZLAGEHkFeBDoBdwKyIXO52IVQHJjLiFiLkKfNvMgEWtspYn9bgAfO+ACQOmmAELHs9eFz8WGAS+NCI855k7lxwqtvhm3lITr7cPAEbeMAMWbi/YfnHtEBmLmEvwItT53siMuF2ILMEVPsH0e+N/xRVzMFoJUMpRvXp1LrzwwgJvX7ly5WKMRilVEtWsWbNQ3xPR0dHFGI0qkYxZgMi5wOPAzcB5iOzGmD0hjqz0ELaCdEFMtEyPq2sSFu4EwMXV+CQWiAR6YO+xf6eewOnAjuKsAADg89XFyFhA8MlnQL6VAHkpfgI+GQomGkhDSENwYaiLcA8+z53AacUadw5aCVDK0aVLF7p06RLqMJRSJVivXr3o1atXqMNQJZ0x3yFyNlAb+1RgPDA5tEGVIiKrcJkEEBfGXAm8bpf7bgQTCRzFcLbMvS3a3P5qisy+qTIeGoAIuFZkHWbsWBe1v68HQCw7TfxCb66nnHJdJOExtfGmJ5vBb+8HkAVxERykFj68VGOviV+YEbBLOoYCPeGRGT16IjIUiETMQYz3cYzrK4yJQrgWF/0RjmkxkLFXhVG3ai18kWlmwBv7cj32tLhKGF8NjO8IA9/ea4LEJHNvi+bI0ZpEuZLNXQv3+peHicgZaGVAFa3dxpijoQ5CKaWUChljdiNyHTADmInITOBf+kSgAEz6RoiqACYc5DrgdVkQ5+YAF2E4CkQicpSjaZcBK8mMbIMhFTFpuPhABMNL3Ufi+v5h/IPgHMQnM3tMoO+bTxmDyEvdzrYFcQywER9tMD4frrAwmRVXBw+zOGhuAMnAheGguOWluI24zBOIbxlgECO4eF+S4nyAD5/vcpOw+OvAS5GxY134vp+MIRokGcP5ZsDiwCcHH8isuBfxynJwUoZ8jAcuxGcyEY9LkuK8iJlKVRnjr8hIUreLwfUycDa40hGXmxlxGZJkppkBC/5lt4k7A2EWqalX4HZlkmnCZEbc7xhuM/0WbgoD/ijOz1GVS/2AWaEO4mRlZmZy8OBB0tLSOOOMM0IdjlKqBMrIyODgwYNkZGRQv379UIejShpj5jgTiK0BKgIZiFREG8ryZAYsPypJcT8D5wFXwjUgEwAAIABJREFUA/A352PwAGEIezGmBsjVwEpcvg6IqxJGPHg9a5kV9zguMwwkCjEpGAGkIj7fw8yMj4UFD2JwOy35BrgaI6kYfIDBw0KQK8CkI7hABGMyMeYKhBfs0wjxgNi9wY2YSNxhxzeq1/7+nxiiEXzgmmMGLDgudcj0XbhVZsW1AcAr1+Eyl+AjBZe4AA9CBYzcy0EM8JAk3XI64vkAQwzCYYxEYsgAUxHjuwH4l8y8pSY+zyYwVYGjgAuRDKAZwlp5qdsFOkSoKg4m1AGcjA0bNtCpUydiYmKoWbMmbdu2zVqXnp5Or169GDNmTAgjVEqF2tq1a7n22muJiYmhVq1aXHvttVnrjhw5Qs+ePXniiSdCGKEqMYzxAYnANcBQYCUipfrv5Cnyji1om1iZ3bMOPnMVEAUkY8wbGAzGdAJAXNcBbnCesvi4DyMVgd8x0hbjuwzMfzFSEZEhktS1tt3PiJOGcxRjeiO+1gjPYGhrC9VMI4IzCQ87G2P6gKzHI58idAZSMRzFZ/rg8TXD6z2btMPfHncVRv7pnOcoRr4+br1/s74LD9h/mHSQJbi4gUwaga81yEYgHGSoCAbJbOOUtNJw0wVvRG3c5gKMGQ2ub+w9yJwEVMFIGi73MPA2xmXigSNARXBNCAP48ssv+emnn3QYM3VSMjIyaNSoEZdcckmoQzlhr7zyCn379sXr9RITE0NUVNQx6yMjIzl69ChPP/00999/f67DAiqlyq4XX3yRIUOG4PP5qFKlynF/OytVqsSBAweYOHEiI0aMOO57RJVDxixC5ELgPqAj0AqR3zFmbz57ll/+fgFiDJm+K3HJDQjhiKzFsBqkL8Y0keldamBohCAgKzHhF4NkIGTiYoDpt3ALgMyI64/PLAPAFXYZ4v3enshkgjxt+i98G0CS+j8GB0diTBhIV7xmOz553wxYsBhYbLfpsRN/6r34dpl7Fv839wtxVUZ8YbhMBj5253vdu/aOpXqt2oT7+hMmw8BVGxexCOngMrwUVwe3eJ3Te/AygjDPy2SErzH3vD7J3joMM11dEAkD8zU+z2FwtUV8IOZzoB1waRjAhx9+yPbt22nRokVhPyKlsvzwww9s3bq11FYCtm3bxsCBA4mNjWXmzJl07tyZvn37smLFimO269mzJ0uWLGHNmjXakVipcubHH3/k3nvvpUaNGrz88stcd9119OzZk2+/PbYBsEePHqxcuZL169fToUOHEEWrShRjNiHSGDu78Argfmx/ARVMVr8ACQfTCaE1kIrLvIfwMWIiwaThCrsfIRVIA97H+BqBiQKBzIjslnd3+ldIZBQG43TadioBchTYknXaATMyJan7a4irB0hjRJ4B39OSFLcHZIAZsOjYQkF+fHIY4/LYArnUznf7WqffivFOR4hEcGPwIsbYSocvGbe7Dkd8q4kmGaiB4Xp8tMWdESlJcR/j9t7NbHcyIm5AMJwPrnnZ9xWcGVFrhQFUrFiRDh06cNNNNxXqupQKtHLlyuP+EJYm8+bNIz09nddee43OnTsDYMzxT2ybNGkC2EpDeSQiI7B/xJQqKgI8YIw5HOpA8jN37lw8Hg9z5syhY8eOgH5PqEIwZh8iXYCJwEJE5gL3YUyuo7+UV8f2C/B1BzLAuPH4PjL3LP5bZsT/iU8agAzCODn6Pu86XGHxgCBGiMrMztF3udx29CADQkZW4rIRwedKP+bku87tR+0fPwe5B8M/nL4BZ4HrLZnR7eKCjQmUZQuIC6iIy3URMCfYRjItrpYZtHA3mKnYtKd9GJkFZhdIU5DeGAPicZn7FqfKi7c2x53xCNADQxXE+DBchde9lnBa4xUDxgeSjI/9QU6ZrqMCKeX46aefcLvddOrUKc/tatWqBcDu3fk/1Suj4oCLQx1EeebxePjuu+94/PHHueWWW+jWrVvQgmgpMwEo8ZWAn376iQoVKuTbuq/fEypXxsxEJBz4BNthMxWRyhiTHOLISqJ3EJrhcgkiUUByVuqNz7cSYwaAqWg3de00CYv2yEtxWzGkY4zBI5cCNgUoM6wtmDQEcLE1t4K8TIurRNjPUWbAm0lAkkyLq4Tb3IiRl4EIMDeA713EPlPA4M7zCgYs/JGkuIMYKiFyu8yKm2z6Ltx6zDmn9zgXl2+5zOx6OYIbMckYGWj6L7IpSjPjx+EjGkix28fVxVVpjxkwYzgwXKb1rI/b+zDG9EKoitecjvjSMUQi8qsZuOiYFA0RDLPiztJKgFIOY4zzhCxvf/75J6CTAKlTLz09na+++oq5c+ciIowcOZLExERcLhddunQpCxWBEk+/J1SRMCYTkSexaUFjgYuA/wtpTCVRVr8AiQTCEdYHrF0J9MYQgUg4hpUAZERsICozDJFohJnyUnf7xShmBhADHMUTvh53es2g5zTUBc9mmd5tAW7XIgw7ESrYPge+TIT9uNmJVyLBeHGZgZIU1wiAcN4OHIcfwBhEknxDMa7XgAp45StJin8O4/vSHleux/h6AJFkuA8RTjhIGGLulJk9duGjJT7fKCA866AuboWDD0pS/AuIbwNuScFLBOIziAnHk7kf4/43xvRDzLmSFDcP4RVEMnGZc5nBcES26OhASjkaN26Mz+djw4YNWcuCFarWrVsHwHnnnXfKYlPlW3p6OuvWrWPYsGG8+eabJCQkMH36dC655BKmT5/Oe++9x/z58/H5fKEOtcxr3LgxaWlpfPnll1nL9HtCnRBj3gZaAH2AYYhcjkj10AZVwpj0jUAFbHrMUVzmnax1EWY9SDhCBGJS8PE+gBk67zBehmNMCkaqY8xcMHMxVAdSMDLS3PP6QXt8cSPB5soyLoyrF8Jr+FiNTxLtsKDmCOFmqem78AAu8xkQhpiuwIuITMPjOivoZQxYvATMEwjpGFMJkX8hrteA2RjXHYABs9cMWngE4QObISnXIr73ML7nATeGFDjmqUMsyBiMeQuf730MvcC4MawyCW/vwWtGI2wD4wbiMSzEZZYCzyI0BBOulQClHN26dcPlcjFs2LBcH+Fv2rSJJ554gqpVq2pnP1Xs0tPTWbVqFYMHD+bdd99l8ODBTJ48mebNm2cVPCtXrsykSZP45JNPePXVV/F6c50QUxWBuLg4jDEMHjyYffuCp3F/8sknTJo0idq1a3P55Zef4ghVqWLMBmwfq+rA+8C1ee9QvpgBy48CyxD+i/AnLvkwa91dC/cifITwXww7iGBd1rqEBTPx+Xohrp+AaAzR+PgFTG/Tf9GLAPjC0xD+BPMXwt9ZJ/VE7AImY/jZVjCIwUg64luEi0uyWvrDMm/GZ54F+QzYBuY3fL6UXK+l/4InMa6r8PEORo6AxNhUJvNfYCTpyWcCkB7RG1yvAukgMYjsB5mOjx34+AtfWAriWgHyOsg+7PwTlTDsw7gmEBEdD2AGLTxChaiLwPcchj1AZUQqAgcwZhZuHkFEJDExUZYuXSpKnYwVK1bIpEmTRET6nfx//aITGxu7ddOmTQW6huHDhwsg1apVk8GDB0ubNm2kWrVqMmXKFOnRo4eEh4cLILNmzSr0/RkxYkQG8ECIb8dJE5FPC33xqlBSU1Nl+fLl0qdPHxkzZoz8/PPP4vP58tzn0KFDMnz4cJk1a5Z4PJ5TFGmRCuWMfO1atWp1sKCB9u/fXwCpWbOmDB06VFq1aiV169aVyZMnS7du3SQsLEwAmT9/fqFvAjiTFanyRWQgIkMQqYHIfESqhTokVfZpnwClAjzzzDNERUUxadIkEhMTs5YPHToUgAoVKvDcc89x9913hypEVYalp6fzzjvvsHjxYpo1a8aYMWNo0KBBgXL9Y2JieOyxxxg7diyzZ8/m7rvvxu3Ou7+aOjGJiYlER0czZcoUpkyZkrV8+PDhgJ0n4IUXXqBHjx6hClGVNsa85OS9fwHsBI4iEoMxf+ezp1InTCsBSgVwu91MmDCBhIQE3nrrLTZv3syRI0eoVKkSLVu2JC4uLmvUD6WKSmZmJm+//Tbz58/n4osvZuLEidStWxeXq3AZm5UrV2bcuHE88sgjJCUlMXDgwEIfQ+UvPDyc5557jiFDhrBkyRI2b97M0aNHqVy5MhdeeCHx8fGcfvrpoQ5TlTbGpCPyL2AV8BzQEE0PUsVIKwFKBVG/fn3uvffeUIehyjh/4X/27Nlcc801JCUlcfrpp5/UKD/R0dFMmDCB0aNHk5iYyJAhQ3TUoGLSoEGDrNZ/dXJEpANlIF2yKCyBR/vBue/BV7+LbLoMkutCZqjjKgOeMMasDXUQJYlWApRy/Pjjj0RFRdGgQYNQh6LKuOTkZJYtW8brr79Ou3btWLx4MZUqVSqy40dGRvL0008zatQoJk+ezL333qupQUVky5YtVKlShfr164c6lLLmLKB9qIMoCToD7YC1cOntwGtAt1McQ2pqKlu2bKFVq1ZlqRHh1VAHUNLoc2KlHOvXr6dx48Z06NCBhQsXkpmpDS+q6IgI+/fvZ86cOdx9993s37+fefPmMWLEiCKtAPi53W6eeuop9u3bx7PPPqu/z0Vk5cqVnHXWWVnfEx6PJ9QhqTLoNGzHgGnYmRm7AgeL8XwiQkpKCu+99x4DBw6ke/fujB8/njFjxpCSkuuAN6qU00qAUo6mTZtSu3ZtVq1aRXx8PGeddRZjxoxh+/btoQ5NlXJ79uwhKSmJoUOHkpaWxksvvcTQoUOpUqVKsZ7X5XIxbtw40tLSeOaZZ0hPTy/W85UH55xzDjVq1Mj6nmjYsCHjxo1j586doQ5NlTEJwC3YpwIu7DiQh4rw+F6vlwMHDrB06VIGDx5Mz5492bJlC8OHD+c///kPy5Yto2XLlsTHx/Phhx9qhbcM0kqAUo4rr7yS7du3s3LlSuLi4tizZw/jx4+nQYMG2uqnTsiuXbuYOnUqI0aMICIigsTERBISEqhateopiyEsLIxRo0ZlPRlITU09Zecuizp27MiOHTuyvif++usvHn30Uc444wz9nlBFrgLwPDAP6AvcfpLH83g87Ny5k0WLFjF8+HD69+/Ptm3bGDJkCMuWLeOBBx6gadOmWSlA3bp1Y9asWbz11ls88sgj/PXXXycZQSF8lASv9oPF2lWkuGifAKUCuFwu2rdvT/v27fnrr7949dVXmT59OqtWrWLVqlXUrl2b22+/nYSEBM4888xQh6tKqJ07d/L666/z3Xff0aFDB6ZOnVrsrf55iYiIYPjw4SQmJvLkk0/y0EMPUaFChZDFU9q53e6s74kdO3Ywb948pk2blvU9UbduXXr37s2gQYO074A6adcDK4CPgQ+B5cClQEHHn/J4PGzfvp2NGzeyceNGjhw5QqtWrUhISKBp06b5jiBWu3ZtJk+ezLvvvsvIkSOJj4/nuuuuIyIi4mQuK3/JeyD1b/BqKmNx0ScBSuWiTp06PPjgg/z2228sXbqUTp06sWfPHp566ikaNWrEzTffzAcffICd30cpW/h/4okneOihh2jQoAFTpkyhd+/eIa0A+EVERDBkyBCqVauWlSKkTl69evWyvicWLVrENddcw65du7K+J7p3786aNWtCHaYq5a4BvgM2AHHA5/ls7/F4+OWXX3j55ZcZNGhQ1lPAIUOGMHv2bAYNGkSzZs0KPIRwWFgYN910E5MmTeLzzz9n1KhR/PHHHyd3UUXJ54HkvZB2OI+NvEZm96wj0+Pqypw+UcG2kKS4GJkZ10Bm31S5mCItUfRJgFL5cLvdXH/99Xi9Xnbt2sVXX32F1+tl2bJlLFu2jBYtWpCYmEjbtm1DHaoKkd27d/Piiy/yyy+/cOuttzJs2DCio6NDHdZxwsPDGTRoEElJSYwePZqJEycWf2teOREeHs6NN96I1+vlr7/+4rvvviMzM5PFixezePFiLrroIl588UVatWoV6lBVKRWN7Rw8F2gK3IgdOSjGWe/1evnll19Yu3YtX375JVFRUVxyySU8+OCDnHHGGYSFnXyRr3bt2jz22GOsXbuWUaNGceONNxIfH18kxz4hnnT4ZDZs/QRcbhAfVIiBtv2g/vl2G58XPp0L36+Yic+TictAZkqYJMVvxiU9TL+Fv8vMHl3x+RKBqnjJxBcZKTPi/ofPPGAGLngjNBdX/LQSoFQetm3bxqxZs3j55ZfZtWsXABdeeCEJCQnExsYybdo01qxZQ7t27diwYQMXXXRRiCNWp4qIsHv3bhITE/n+++9JSEjg4YcfLvGF6rCwMBISEpg1axZDhw4lMTExdH/Ay4itW7cyc+ZMXnnlFfbs2QPAJZdcQkJCApGRkUybNo3169dzxRVXsGnTJs4555wQR6xKq0FAMnAutsNwhAif/vILn73/Pp999hmnn346bdu25dFHH6VOnTrF8n87LCyM9u3bc8EFFzB9+nQGDhzI6NGjadiwYZGfK08i8O4TsPc3wBnG1BVmnwisfBauuR/qtYDPX4efVgMCLiMIBiEdkZZ4TB15qVs1fDIXiEJIwWUi8JEK1AGuALQSoFR54fF4WL58OTNmzGDFihX4fD4qVKhAnz59SEhIoHXr1lnbdu3alXnz5tG7d29eeeUVrQSUAz6fj23btjFz5kx++eUXEhISGD9+fKkaS9vlctG/f3/+/e9/M3DgQKZOnap9BAopIyODpUuXMmPGDFavXo2IEB0dTb9+/UhISKBly5ZZ2/bo0SNrBue5c+cyceLEEEauSrtKIozdsYPdixbR8NJLqRQZyaSzzuLWW28lNjb2lFXqq1atykMPPcTnn3/Ogw8+SMeOHbn99tsJDw8/Jedn+1ew/w8wBuq3gItvh7//gpWTwZMG62fALdPgtw3gyYCYmt+T/NeNCF5MWHOQe3AbLz5zE0gFkBRwX0FGxg7crjMw5jaMK6/8olJPKwFKOXbv3s20adN4+eWXs0ZAOPvssxk4cCB9+vTJdUSXXr16MWjQIHbs2HEqw1WnmP9R+9y5c9mxYwd33nkn48ePL9WTcN1xxx243W6GDx/O008/zWmnnRbqkEq8HTt2kJiYyJw5c7Ja/Zs1a0ZCQgK33347MTExQffr168fQ4YM0e8JdUI8Hg+//vorS5YsYf369Zx55pnUuPtuKl9wAWuAjS1b4uLUF+pcLheXXHIJc+bMYfr06dxxxx2MGzeOxo0bF//Jd34LmangjoArB0FERTitJpx3PXy7DI4egrRkW0kASDnYEJ+7N8a33PRf8AHwAYAkxXfEiCAmDCODCXPPJyPlYzP0vfuK/yJCSysBSjmWLVvG+PHjCQsLo0uXLiQkJNC+ffsCtfC2bt361D8KVaeE1+vl+++/Z+7cuRw8eJDbbruNtm3blpkUml69euFyuXjwwQd54oknTunwpaXRggULeOqppwgPDyc+Pp6EhASuuuqqfPdzuVxcfPHFnHXWWcUeoyr9RITU1FR+/fVX3nnnHT799FPq1atH586dGT58OFFRtl/rKOw0uPcC7wL/F6J4K1WqxMiRI9myZQujRo2iXbt29OnTh4oVKxbfSQ84c/hUiLEVAL+aZ4M7EhA4vBua3wBfzIfMtNNAxiKuf8lL8Wm4XCNN//mv4PXOJcx1H5hokLswJo7ISpGS1P11wjKGmbuXJRffRYRW2fgrplQRqF69Oo888gj9+/enbt26hdp3xYoVxRSVChWPx8PXX3/Na6+9RmZmJrfddhutWrUqM4V/P2MMPXv2JDw8nNGjRzN+/HhOP72ggw+WP7Vr1+bxxx+nb9++1KpVq1D7rl+/vpiiUmWBiHDo0CF++uknVq5cyQ8//ECdOnXo1KkT9913X1bBP1AFIANYBNQEOgILgVANbXPeeefx2muvMWvWLPr27cvIkSM5//zziydd0uU8hRXfsct9PkAAA8YF514PlarB+tm/kXqoPoIPI6cj3qkyI+6A6b9wmbzU7UKMGQemE4IbF27E9MITVRXoUvTBlwxl66+ZUiehS5cudOlSZv+vqwLyeDx89tlnzJs3j8jISHr37s0FF1xQqtN+8uNyuejWrRthYWE89NBDTJgwgerVq4c6rBLplltuCXUIqozZv38/3377LatXr2b79u3Ur1+fDh06MGLEiAKNMjYYO5PwP7DDh4YBB4BQPdOLiopi8ODBdOjQgYkTJ3LOOefQt29fYmNji/ZEVc+AnVvssKCpf9snAgC7foDMDJsGFFMLjuyDsy6Gsy4exyt3LiDzyP8hvAmmMj5uk6S4b8yAhb8Ct4hgmNXtPLyuRAxtQa6XsWNdZuxYX56xlFJaCVBKKWyH348//phXXnmFatWq0b9/f84999wy1/KfG5fLxc0330x4eDj3338/zz77rD4RUKqY7N27l6+++orVq1ezd+9eGjRowPXXX895551H5cqFb8evAswHLgY6A9Wxw4eGUtOmTZk+fTpvvPEG9957L3379qVNmzaFb1DxZsLn845fLgJhkXYI0JXPQquecGgn/LgSxAtV6kJENMztZ0cJOq3mRaQn/4DLuO0jAsC49mJ8z8qMuPPxmWeZ6fsWr8vgFoMAyNGyWgEArQQodRyPx8OHH37I5s2b2bdvX67b9evXj0aNGp3CyFRxEBHWr1/PtGnTOPPMM3nggQdo3LhxuSn8B3K5XHTq1ImwsDAGDRrE9OnTtY9ALjIzM1m1ahVbtmzhwIEDuW6nswYrv3379rFhwwZWr15NcnIyTZs2JS4ujmbNmlGpUqWTPv5VwFvA/4B52PkEOgHVTvrIJy4qKoo+ffrQpk0bEhMTWbNmDYMHDy5cA4M3E75ZevzyZldD3eawYzP871dY8Qx4veDJtJWDqwbZ7cRnRxIypj8u120I4UBFIBmX5yV87vEgjXExGTHpuMQgRGM4iuH5orgPJVX5+yunVB6+/vpr4uLi+O233/Ld9pprrtFKQCnm9XpZv349U6ZMoWHDhjzzzDPUq1evwDNollXGGDp27EhUVBR33HEHc+fOLfrH+KXcp59+So8ePdi+fXu+23bt2lUrAeXYgQMHWLt2Le+99x7p6em0bNmSvn370qBBgyIp+OfUFVvwnwY8AqwALivysxSOMYYmTZrw9NNPs3z5cu655x769+/P1VdfnXdfgdh6cHhP7usr14Q2feGbt+H7D7JHAqrVFC67E05vYLe7oDts+xwObHchxACpwHp8jDQDFn8rM7pPx2fCwXcJmNMQ8WHMdjBP0nfhLPoV5d0oWbQSoJQjOTmZTp06sWvXLjp37swFF1zAI488Qu/evTnnnHPYsmULb7/9NlWqVGH06NH84x//CHXI6gSkpqby0UcfMWfOHM4880ymTp1a6I7g5cFVV11FZGQkd9xxB0lJSdSuXTvUIZUI+/bt44YbbmD//v3Ex8fTpEkTxo8fz913303Dhg3ZvHkzS5YsoWbNmjz44IM6GlA5dPDgQVasWMG7775LWloabdu2ZdSoUdStWzdo596iFglEAMuBKOBqYBl2xuFQioqKIi4ujjZt2jBx4kSWL1/OuHHjch1Wl8vutK/8XNDdvnze7M7CgS7sbl/Ql4Xx80z8Qm/gatN/0XvAewCyIM5t4hdlr+9f0KsrpUREEhMTZenSpaLUyVixYoVMmjRJRKRE1ZtjY2O3btq0Kd/4k5KSBJChQ4eKiMjOnTsFkKSkpKxtfvjhB6lbt6507ty50PdnxIgRGcADIb4dJ01EPi30xZcAycnJsmzZMundu7eMGTNG/vjjj1CHVCp88sknEh8fL9u2bSvuU50Rwl/rdq1atTpYkCCfffZZAWT06NEiIvLzzz8LIPPmzcva5ptvvpHq1atL7969C30TAB9Z05+WHyLSr9A3q4Twer2yf/9+mTdvnvTq1Uu6du0q06dPl+3bt0tmZmbI4torIlVF5CEROSwi+0MWyfG8Xq988MEHcs0118i7774rXq/3VJz29lD/npc0+iRAKccXX3wBwP3333/Mco/Hk/XvZs2aMWnSJG655RY2bNjAZZeF+kGrys/hw4dZuXIlS5cupUmTJkyYMEHTMwrhsssuQ0QYPXo0jz322KmZBKgEK8j3RIsWLZgwYQL9+vVj5MiRNG/e/JTGqIpfRkYGe/fuZfXq1axevZojR47QsWNHnn76aerUqRPq8AA4HfsEoDlwHXb0oFkhjSiby+XimmuuoXXr1jz++OO89957PPDAA9StW7dUzb5e2mklQClHcnIyERER1KtXDyArX/PQoUPHbHf11VcD8NFHH2kloAQ7dOgQ77zzDitWrOCf//wnTz75pKb9nKA2bdrgcrl47LHHePjhh2natGmoQwqZ5ORkYmNjs/pJFOR7QisBZUN6ejp//vkna9eu5bPPPiMtLY127drxxBNPZP3dKGnaAG8CHuBZYAbQndANH5pTlSpVeOaZZ1izZg0jRoyge/fu3HjjjackbUppJUCpLDVq1CAjI4OjR48SHR3NaaedRuXKlY/rJHzkyBHAtjCrkufQoUO89dZbfPjhh7Rq1YoJEyZo4b8IXHrppbjdbsaNG8fDDz9Ms2bNQh1SSNSoUYMjR47g9Xpxu93UrFmTsLAw/Z4oo9LS0vj9999Zs2YN33zzDcYYrrzySh555BHq1atXKlqte2A7DD8BPA9cRMmpBIB9KtC+fXtatmzJ1KlTWbduHffeey9nn312qEMr87QSoJTDn+awdetWWrRo8f/s3XdcVfX/wPHX594LiOBMUstJ5dYcaTlypJaZTUUqU3OB5jZxI7hFyz1Ay5VagdrOLBxZP78tS0tcaZojZ25k3PH5/XGQ1CBBgXOB9/Px8FGfy/mcz/te9NzzPp8FQJUqVfj44485ceJE6sTIyMhIAPz9/c0JVKTp8uXLrFq1is2bN9OiRQveeOMN7r777lzxJZ1bNGjQgKFDhxIWFsb48ePz5eT4+++/H7vdzuHDh7nvvvuwWq088MADrF27ltDQ0NSlD+U6kXslJiayd+9eNm/eTFzfxXHWAAAgAElEQVRcHAUKFKB58+aEhYXl2uEqHhh7B2zAWBqnecr/e5kY083uuusuxo4dy7Zt2xg/fjwtWrSgU6dOeHm5U5R5S/5eC0+I6zz++OMAREdHp74WFBTE2bNnqVu3Li+++CJNmjRh2rRpFCxYkCeffNKsUMV1Ll++zPz583n55Zfx8vIiMjKS3r17U7JkyVz5Ze3u6tSpQ3h4OMOGDWP//v1mh5Pj0rtOHD9+nNq1a/PSSy/xyCOPsGDBAgoXLkyrVq3MClVkgt1u58cff2Tq1KkEBwezcuVK7rvvPsaNG8fcuXPp0KFDrnnyn56+QDmM+QFPAFeB86ZG9G8Wi4UmTZowb948zp07R58+fdi9e7fZYeVZ0hMgRIoqVaowc+ZMihYtmvpa9+7d2bVrF3PnzuX9998HoHjx4ixZskSGmJjs4sWLLFu2jA0bNhAYGEh0dDTe3t5mh5UvVKtWjVmzZvHaa68xa9asfNUjUK9ePaZPn37Dkqn9+/dnz549LF68mPfeew8APz8/Vq5cKbsuuzG73c4vv/zC+vXriYuLw9/fn8cee4zu3btTokSJPLlnSGlgI3A/0ApjD4G5pkaUtiJFivD666/zyy+/MHXqVB566CF69eol1/gsJkmAENcZNGjQDWWLxcKsWbMICQlh586d+Pn5UaNGDbkQmURrzenTp1m1ahVbtmwhMDCQDz74QLqLTeDv78+iRYvo06cPEydOpE6dOmaHlCOUUv9aGchqtRIVFcWoUaOIi4ujZMmSVK9eXSY3uiGHw8H27dv5+OOP2bVrF9WqVeOJJ55g4MCBFC5cOE/e+N+sPvAOUASYCswGugDutiWgUoq6deuyYMEC3nnnHbp168aoUaNkon0WkiRAiAy49957U5/8Hzx4kEOHDlG3bl2KF3en6VV5l9aa48ePs2rVKrZv30779u2Jjo6WmyyTlStXjvnz5xMSEsLQoUN5+OGHzQ7JVOXLl6d8+fIA7N+/nyNHjtCgQQMKFy5scmT5l9Yau93Ozz//zAcffMCOHTuoW7cuzzzzDKGhoXh5eeXqIT63qzPwMsbGNSswegWK4p4bVPj6+tK7d29atmxJeHg4tWvXpn///vIwLgvk/ZRXiCy2fPlyWrduzY4dO8wOJc/TWnPo0CGmTJnC8OHDqVSpEsuXLycwMFASADdRoUIFpk2bxpw5c/jmm2/QWpsdkluIjIykdevW7N271+xQ8h2Xy8WFCxfYvHkzI0eOpH379nzyyScEBASwfv16pkyZQsOGDSlQoEC+TACusQL3YQwP+gtoBthNjSh9SikqVarEihUr8PPz4+WXX+aHH37A5XKZHVquJj0BQgi3o7Vm3759rF69mmPHjvH8888zZMgQufF3UxUrVmTSpEmEhYVht9tp0aJFvr65EjnP5XJx5swZfvnlFzZt2sSRI0eoXr067du3Z+LEidhscruTlteAI0BDYBpwDvDE/YYGXWOz2ejWrRuPPfYY48aNw9/fnz59+nDXXXeZHVqa9OKAimiiVFDM42bHkhb5VyGEcBtaa/bs2cPy5cs5f/48HTp04NFHH5Vu31ygQoUKTJgwgTFjxuByuWjZsqUkAiJbuVwuTpw4wfbt29myZQtnz56lRo0aBAYGUrNmTTw9Pc0OMVcoB2wDSgGPYawe9IapEd1a+fLlWbRoEdHR0QwYMIBXX32VFi1auF+yZ8cTK+XNDiM9bvZpCSFM1AZjY8mNQI6P6di9ezcLFy4kKSmJTp068fDDD8uT/1ymXLlyTJ06lZCQEABZHjPv+xxjBMlsYDM5cN3QWnP06FG+++47vvnmGxISEqhRowZdu3alatWqcuN/mx4ElgAVgIlABBCE+/YIgNEr8PLLL9OoUSPmzp3Lli1b6Nu3L/fcc4/ZoeUakgSI3KI4xvXpZ5PjyMsuYnypHwW+x/geOJDdje7evZvZs2fjcDgIDg6mdu3a8kWei91zzz3Mnj2bPn36oLWmdevWZoRRGzgGnDWj8XxkNhCN8QD5D+A7YDrZcN04cuQIW7duZcuWLVitVurUqUOfPn3w9/eXhwVZpDvwKtAH+BRjl2F3nSx8vQoVKjBlyhS++OILQkJCCAgI4Omnn8ZqtZodGliL/YHr0hNmh5EeSQJEbpEArAXOAN9gDF88ZWpEec//gF1AE6Am8CywDyMxWAhcyKqGXC5X6v4LDoeDAQMGUKtWLfe4aIs7VqJECZYsWUK3bt1wOBy0adMmp4cGlQXWAzuAD4FlQFJOBpCNLBjDtt3B1xhzSqsAtVL+BACHgC+AmdxBIvbXX3/x8ccf8+233+Lt7U3Dhg0ZPXo099xzjywLnE0sGBn0ECAOIzHYgjGJ2J15enryzDPPUL9+febNm8eGDRsYPHgwlSpVMjcwy/kyaBZi9LS7HUkCRHbyxditPKv8H9AJY5njTsCvwGdAFJCYhe2YrRUQalLb1+/65JfypwnwOkYvzJg7Obndbmfnzp0sWrQIrTX9+/enZs2aMnY8DypUqBDLli0jODiYxMREnn322XTXYNdas3z58iJk3eiDbzEeErQBHsf4+7sDWAx8lUVtmOUVYLnZQfyHYil/6mKMKPkI6Hk7J/r999/ZunUrb7zxBn5+fnh4ZOXXiUhPH+Ag8BLG058TGF/mRf+rkjtwOSntOMPEx+5m2+4ENm3caH4SYMwJqGhuEOmTJEBkp8+BR7Pp3CWB1hg3zH0xuqOPZeYEP/74I6NGjcp0wwcPHsx0nUw6h3nDnu4Fbt7iNAkjph0YXf6ZlpSUxPbt21m+fDk2my315l/kbb6+vixYsIDXX3+d5ORkOnTokGZvz/nz5+nWrduv2RSGBXgg5c9TwCbg6YxW3rp1KxMmTMh0o/v37890nQz6P6Bbdp38NozH6Hm53kWMYYWLMIaa3xYfHx8eeOABGeNtgvuAn4DCQHOMhCDz/wqymXbB2T/gwLdw7ghcOA7xF1AFi3LSrykFk51mR+j2JAkQ2aklWbsXxWKMPU4A/gR+wRiPGg1k+l/733//TWxsbNZFl3V+xpwkQGE89Qfjxv9AShxvAjtv54SJiYls27aN6OhovL29GTBgANWrV8+aaEWuUKRIEd544w1Gjx6N0+kkMDDwX4lA8eLF2b59+wP16tU7moVN/w+ogzFxdS/wIzAP47qRYadOnXK368TBlD/uoBTGPFIwhgseBGIx5gqcMCsokTWqY2Rx9TC6gMcBA3GTHgGnHfZshL0bIekyXL1gJAXFysJzE3Gs+5jbuC3IeiU4xEX3HAoEkgSI7JWV+474YvQqbAU2YHyZX7qTE1avXp2oqKjbrl+lSpVbH5S7dALKYwyXeBtYw21eRRMSEtiyZQvr1q3Dz8+PAQMGULVqVRn2k08VLVqUyZMnExoait1u55VXXvlXIlC3bt1ksm7cfgBQGvgEY+jMB8Bt7SpUr169O7pOVKhQ4bbr5gJvAsnAaoxVJTOVYAn3F4QxlqszxjKi3TB6B0zfadbqAQ80gUsn4I/vweYBRcpAuzHg4UZLSv9tvQeLcwHGyqtuR5IAkVtYMYb8HMqqE5YtW5agoKCsOl1e8B1wP3eQXCUlJbFhwwbWrl2Lv78/Q4cO5YEHHkh3LLjIP4oUKcLEiRMZO3YsDoeDbt26Zeffi18whv9cudMT+fv7y3UifVOBPRhLC4s8ygI0BSZjPIVbhDFZ2LSrussJ+zbDL2vh/ibwwKNw6HtoOwo8fcyKKm3a6QX4mx1GeuSbWeQWF8nCBECk6QC3mQAkJyfz4Ycf8sorr7Br1y7GjRtHaGgolStXlgRApCpcuDCTJk1i165dd/R0PQMOkAUJgLil35AEIF8IBq5i9AoMAA6ThcvFZcbx32DtMPhzOzw9Hhp0AmcyPDseChQyI6JcTb6dhRC3zeFwsGzZMgICAti7dy9z585l5MiRVKhQQW7+RZp8fHyYOnUq+/btY8qUKWaHI4TIoKoYS/LVBpoBC3Ky8YsnIHYmfL8SHu0JbYZDIT/jZ/VfBG+3mKnwbyU4BJa2ZoeRHvmWFkLcNq0127dvp0GDBgwdOpRSpUrJuH/xnxITE/nggw+Ii4vDbrejdY5vTi2EuE2VgC8x1t0dAozC6KbPNvZE2B4Dn4RDqcrw/BQoVfXGY9xpDsDN/rbeA67ZZoeRHkkChBC3zcPDg5kzZxIfH8+MGTNwOGRkgPg3rTUXLlxgxYoVBAYGcuLECdatW8fYsWMlaRQil3kNY6m+QIwZ9wlkw5gwrWH/1xA9GJKuQMeZUKMtqFx222rMCXjA7DDSk8s+TSGEu7HZbEyYMIErV64wbdo07PasXBRK5HZnzpxhyZIlBAcHc+nSJZYvX87gwYMpVEjG7wqRW1mAdhibbnyIsWFPlvXpnf4dPhoDv2+FtqOhUTfwLJhVZxfXkSRACHHHrFYr4eHhAEyePJnExLy0gbO4HSdPnmTevHkMHjwYrTWRkZH069ePokXddOyuECJTemFsyz0YGA7sI5NDg1w39R/En4PN82DLfKjbAZ4KhWJlsihak5TgENrVzuww0iNJgBAiS1gsFoYNG4avry+TJ08mISHB7JCECY4fP86MGTMYNmwYRYoUYd68efTs2ZNixYqZHZoQIovVBnZjbBndDFia0YqOZPgpJuX/k2DHR/DRaGOzrw5vQLk62RFuzvvbeg/K8qbZYaRHkgAhRJax2WwMGDCAEiVKMGHCBK5evWp2SCKHHD16lClTpjBq1Cjuvfde5s6dS+fOneXJvxB5XEWMYUEvYmwuNoRb9Ai4HPD5REi8ZCz1uWYonD8Kz0+F2s+CJQ9tYWXMCahsdhjpyUOftBDCHXh4eNCnTx8WL15MWFgY48aNo2BBGc+ZVx0/fpzFixdz8OBB2rdvT79+/WS8vxD5zACM7eXbAGeA0YAd8Lj5QJcTvoiAMwfBngCXTsLjQ6F4+ZwNWACSBIh8LC4ujjlz5tx2/SFDhlC5stsm+Kby8PAgODiYt956i5EjRxIREUGBAgXMDktkoRMnTjBnzhwOHDhA9+7dGT58ON7ebrxU323avn07ixYtuu361/bNECKvswIvAc9iDAtaD2y8/gDtgg1T4NhvgDZW/fEuAmf/hMKlweaZ80FnN8eVw+D7jNlhpEeSAJFvHT169I6+3AMDAyUJ+A9Wq5WgoCCWLl3KkCFDmDFjhiQCuZzWmqNHj7JgwQL2799P3759mThxIlar1ezQss0ff/xxR9eJHj16SBIg8o3uwA9AOPAZxuZi5YEiWhsr/pw+aNzsF/ID3xJQpjbcWyNvJgAABb1L4WQaRm7kdiQJEPlW9erViYqKyvDxV69eZeHChezfvz8bo8pblFJ0794dm81Gv379mDt3bp58WpzXOZ1ODh48yJIlSzh8+DDBwcFMnjw5X+wKXa9evUxdJ65cucL8+fP5448/sjEqIdxXA2APcBpoDkxBE7zjA0BB/UCo1Ax87jIzxJyTbC2AlSpmh5EeSQJEvlW2bFmCgoJueZzdbmfJkiVMnz6dv/76C5vNRteuXaldu3YORJk3dOnSBU9PT/r168fMmTMpXLiw2SGJDHA6nezZs4fly5dz5swZunXrRpMmTfL0k/+b+fv7Z+g6kZycTGRkJBEREZw+fRoPDw969uxJlSpu+/0vRLYpC8RgbCz2Mop+dV5gUp0XKGJyXOJGkgQIkQ6tNWvWrGH06NH8/vvvALRq1Yo333yTWrVqmRxd7vPiiy9itVoZMmQI06ZNo3jx4maHJNLhcDj49ddfWblyJVeuXKFLly40bNgwX938Z5TL5WLt2rWMHDmSgwcPopQiICCASZMm8cADbrtRqBDZbgiQDLQAXIACkgAvM4PKaY4rh/Eo/JzZYaRHkgAh0hAbG8uIESPYvn07AI0bN2bq1Kk0adLE5Mhytw4dOmCz2Rg+fDiTJ0/Gz8/P7JDEdRwOBz/99BOrV6/G5XLRqVMnGjRoIDf/6YiNjSUkJIQdO3YAxkOCiIgI6tata3JkQrgHT4xlQ58H3sSYL7De1IhyWEHvUjhdkzE+ArcjSYAQ1/n+++8ZNWoUmzZtAqB+/fqEhoby9NNPmxxZ3qCU4tlnn8VqtTJ8+HCmTp3K3XffbXZY+Z7T6WTbtm2sWrUKHx8funbtyoMPPojNJl8Radm2bRsjR45k69atADz88MNMmTKFFi1amByZEO6nK7AVmAN8CfwEVALyxaBQY05ANbPDSI9c4YUA9uzZQ1hYGGvWrEFrTeXKlZkwYQIdOnRAKWV2eHmKxWKhXbt22Gw2hg4dyvTp0ylZsqTZYeVLWmu2bt3K22+/TalSpejbty9Vq1aVm/90xMXFMW7cOGJijJ1Oq1WrRnh4uFwnhLiFphg7Cx8A2gILgM6mRiRAkgCRzx09epSJEyfy9ttv43Q6KVOmDKGhoakr2ojsYbFYePLJJ/H09GTAgAHMmzdPhgbloGs3/3PnzuW+++5j7NixVKxYUYb9pOPPP/9k8uTJqdeJcuXKMXr0aHr06CGfmRAZVBpYCYwBnsIYJvQGebxHwHHlMLaCL5gdRnrkLkfkWz/88AOPPvooycnJ3HXXXYwcOZK+ffvKWvY5RClFq1at8PT0pFevXixevFgSgWxmt9v5+uuviYyMxN/fn5kzZ1KmTBl5iv0fNm/ezBNPPIHdbqdkyZKMHj2a4OBgPD3z6LrmQmSjECARaAQUw7gJTQDy7MLRBb1L4bSMB9qbHUpaJAkQ+da5c+dITk4GjI2tFixYwIIFCzJcf+XKlTRs2DC7wss3mjZtiqenJ927d2fRokWULl3a7JDynKtXr/L111+zYsUKKlasyKxZsyhTpozZYeUKZ8+exW63A0biOmvWLGbNmpXh+uvWrePBBx/MrvCEyHUKYKwc9AIwGvgD+MjUiLKRMSeghtlhpEeSACGA06dPc/r06UzVSUhIyKZo8p9HHnmE0NBQ+vbty4wZM2SH1Sxy5coVYmNjWbNmDffffz/Tpk2jbNmyZoeVa508eTLTdZKSkrIhEiFyt1cwJgmvAjYC3wIPAoXMDCofkiRA5FtNmzbl4MGDt11fnlhnrQYNGjBy5EiGDx/OhAkTqFSpktkh5VoXLlzgiy++4PPPP6d69epERERw7733mh1WrtS2bds7uk7I5y5E2h4H4oCfMdbPXAF0MDWibODt8yeJCW45FAgkCRD5WMGCBfH39zc7DHGd+vXrExISQnh4OKGhoVStWtXskHKV8+fP89FHH7F582Zq167NlClT5Cb0Dvn4+Mh1Qohs4oeRBLwBNAFeBeaSh3oEkuPvxkoYEGB2KGmRJEAI4VYeeughQkJCmDBhAqNGjaJGDbcdTuk2Lly4QHR0NFu3bk3d2K5UqVIy4VcI4fZGAvFAfeB+wAO4ChQ0M6is4sQbK7XMDiM9FrMDEEKIm9WpU4cRI0YwYcIEdu3aZXY4buvSpUvMnz+fnj17YrVamTVrFn369KF06dKSAAghcg0fIAyIBvoDXcwNJ9+QngCRb+3YsYPJkyffdv2xY8fKU+psVKtWLcaNG8eIESOYNGkS1atXNzskt3Hp0iWWLl3KV199RceOHVm6dCmFCuWZDnS38t133zFjxozbrj958mTuv//+LIwo/0hOTpaJ1flIIPAp8DmwBWPCcAMyPzRIa8358+fdY68fb58/cVx2y6FAIEmAyMdOnjyZuvPn7ejdu3cWRiPSUqVKFWbNmkXfvn2JiIjI10mX1ppz586xfPlyNm7cyMsvv8zatWvx8vIyO7Q87ejRo3d0nRg6dGgWRpO/FCpUiF27drFq1SpefPFF2ZgtH2gH/IaRBLwMrEl57Va01hw6dIgvvviCr7/+Gh8fH7p0cYP+hOT4u9FqNEaO43YkCRD5VvXq1YmKisrw8ZcvX2bevHkcPnw4+4IS/1KhQgWioqLo3bs348aNo169emaHlKO01pw8eZJ33nmH7777jsDAQNatWyc3/zmkXr16mbpOXLhwgblz53Ls2LFsjCp/qFmzJqtWrWLOnDl06dKFadOmyUT3fKA4sBdYBNQCOgFRgG8axx4+fJgPP/yQTZs2UbJkSZ544gnmzZtHsWLF3KMnwIk3VlXb7DDS4wafkBDmKFu2LEFBQbc8LikpicjISCIiIjhz5gyenp4EBQVRt27dHIhSAJQpU4aFCxcyePBgBg8eTOPGjc0OKdtprTly5AirV69mx44ddOzYkQEDBsiO1jnM398/Q9eJhIQE5s6dS0REBOfOnaNAgQK89tprVKtWLQeizLuKFClCaGgoO3bsICgoiA4dOtCpUyfZsTmPGwVcBOpiTBj2wJg87JGczJ9//slnn33Gxo0b8fPz49lnn2X16tX4+PjIXKhMkiRAiHS4XC7Wrl3LiBEj+OOPP7BYLAQEBMgYX5OULVuWGTNmMGrUKJKTk2nRooXZIWULl8vFwYMHeffddzlw4ADt27dnyJAh8uTfTTkcDlavXk1YWBiHDx9OvU5ERERQsWJFs8PLM2rXrs26deuYPXs2PXv2ZMSIEVSpUgWLRdY3yauKANOB5vHxPGe3c+rECR4ID6do0aI89dRT9OrVCx8fH7PD/G/ePn+SePlFs8NIjyQBQqQhNjaWoUOHsnPnTgBatWrFtGnTqFOnjsmR5W/lypVj8uTJjB07FrvdTuvWrfPMkx+tNXv37mXVqlWcPHmSF154gZCQELy9vc0OTaRBa82aNWsYM2YM+/fvB4zrxJtvvkmtWm67ImCu5uXlRUhICDt37mTKlCk88sgjdO7cmcKFC5sdmshily5dYv/+/ez98ktm+/qyo1Mn5pw8ifeSJTzp45N79hFIjr8bi3UY8JLZoaRFkgAhrrNt2zZGjhzJ1q1bAXjkkUeYMmUKzZs3NzcwkapcuXJMnDiRsWPHorXm8ccfz/WJQFxcHEuWLOHKlSt07NiRRo0ayc2/G4uNjWXEiBFs374dgMaNGzNlyhQeffRRkyPL+5RS1K5dm8jISJYuXcqgQYPo3bs39evXz/XXgfzu/PnzxMXFsXHjRg4fPoyfnx+PP/44rz30EKpoUT5u0YIuwGdAS7ODzSgn3li1244dliRACGDXrl2MHz8+dRWQatWqER4eTkCA267sla+VKVOGSZMmMWLECFwuF23atMmVNwBxcXEsWLAAu91O165deeihh2TYjxv74YcfGDlyJJs2bQKMiauhoaFynTCBj48P/fr1Y9euXcydO5eNGzfSt29f6RXIZS5cuMD27duJjY3l1KlTlClThpYtW1KrVi2KFSt2w7FHgNVAOYyldpaSRzYUM5EkASJfO3z4MFOmTOGtt97C5XJRrlw5Ro8eTY8ePWQ5OjdXunRp3nzzTQYMGADAk08+aXJEGRcXF8esWbPQWtO7d28efPBBPDw8zA5LpGPPnj2EhYWxZs0atNZUrlyZCRMm0KFDh1yZfOYlNWrUYNasWURHR9OrVy/69+9PkyZNzA5L/IeLFy/y3Xff8cUXX/D3339TqVIlnnvuOSpXrkzRokXTrTca+Bt4EGiDsdvtFdJeNchtJHIEH4tbDgUCSQJEPva///2Ppk2b4nA4KFWqFGPGjKFXr16y6kQuUqJECaKioujZsycOh4N27dq57U2Z0+nkt99+Y968eWitGThwIDVq1JCJjW4uNjaWNm3a4HQ6KVOmDGFhYbz66qvusfygAMDb25uuXbvStGlTpk+fzvr16wkJCfnPG0qRsy5fvszXX3/N+vXrOXfuHLVr16Zr167cd999mdro8C5gIdAa6AiUAJZkT8hZw8taAuUYirHtgduRq5jIty5evIjD4QCMC1RoaCihoaEZrv/RRx/JGGA3UKhQIZYuXUqvXr1ITk7m+eefd6sb6+TkZHbu3MnixYtRSjFw4EBq1qxpdlgig86fP4/T6QSMoQvDhg1j2LBhGa4fGxsrywnnkIoVKzJr1iw++eQTunTpwpAhQ2jWrJnbPhjI6xISEvjyyy/56KOPOHfuHE2bNqVfv35UrFjxjpY6fhpjWNABjL0E1gGP46Y9AspZEK3cdnMbSQKEAOLj44mPj89UHbvdnk3RiMwqWLAgUVFRDBkyhMTERLfYXTQxMZHvv/+e1atX4+XlxcCBA6levbqpMYk7c+XKlUzXufagQeQMT09P2rdvT9OmTQkLC+OTTz5hxIgR+Pn5mR1anud0Orlw4QKbN29m/fr1nDp1itatWxMWFka5cuWyNBl7GXgGWAEMBb4C8v7uMVlPkgCRbzVt2pSDBw/edv3SpUtnYTTiTvn6+vLGG28wevRo7HY7nTt3NiURSEhIYOvWraxdu5ZixYoxcOBA2TAqF2vbtu0dXSdkh1tz+Pn5MX/+fNavX09QUBDdu3enTZs2Mvcmi9ntdk6fPs3XX3/N1q1bOXfuHM2aNWPMmDFUqFAhW3thfIFzwAdAYaA9sApwq+0UEzlCQVcns8NIjyQBIt8qWLAg/v7+ZochslDhwoWZNGkS48aNY+nSpTk6djsxMZGvvvqKDz/8kDJlyjBkyBAqVarkVkOTROb5+PjIdSKXUkrRtm1bGjRoQEREBJs2bWLQoEGUL1/e7NByteTkZI4dO8Y333zDd999R3x8PI8++ighISH4+/vn6PCrMcApoCbQCXDhZpOFvawlUHog0NnsUNIiSYAQIk8pXLgw4eHhjB8/nrfffpuePXtma49AcnIyn376KWvXrqVatWqMGjWKihUrys2/EG6iRIkSTJs2jY0bNzJq1CieeuopXnjhhTsal57fJCcnc/DgQb755ht+/vlnXC4XjRs3JiQkhIoVK5o676IkxrCgphjzBSphTB52C8pZEFQDs8NIjyQBQog8p1ChQoSHhxMWFkZkZCR9+vTJ8pvypKQkPvnkE959993UJ4333nuvTEIUwg0ppWjVqhW1a9dOnT80cOBAKleubHZobsvhcLBnz2Q5v+kAACAASURBVB42b97Mjh078PHxoXHjxowYMYJy5cq51YOONsAyjOFBEzEmDj+DG/UIuClJAoQQeZKPjw+TJk1i9OjRzJkzh0GDBmXJeZOTk1m3bh2rVq2iadOmLFy4ED8/P7n5FyIXKFGiBCNGjOD7779n0qRJPProo3Tp0kU26UvhcDiIi4vjq6++YseOHfj5+dGsWTPat29PqVKlTF9w4b+8irFs6DxgAkaPwENmBgTGnABv9xwKBJIECCHyMC8vL6ZNm8bIkSOZPn06r7/++m0/vYqPj+eDDz7gvffeo2XLlqxcuZIiRYpkccRCiOxmtVpp1KgR1atXZ8mSJXTr1o2wsLB82yvgcrn49ddf+eyzz/j555+pUKECLVu2pEuXLpQoUcKtnvjfSkHADnwOaOBZIAYwbfcfL2sJrK6+wA9mhfBfJAkQQuRpFouFqVOnMnbsWKZOnUpISEimVgi5cOECH3zwAR999BEtW7Zk1apVcvMvRB5QpEgRBg0axO7duxkxYgQtW7akd+/e+WIjOK01v/32G+vWrWPbtm3UqFGDtm3b0r9/f3x9fXPVjf/NRgPHMCYL9wcSMRIDHzOCUc6CaB4xo+mMyL2/ZSGEyCClFGFhYWitiYiIICkp6ZZ1Tp8+zaJFiwgODiYpKYkVK1bQv39/SQCEyEOUUlSvXp01a9YAEBgYyM6dO3G5XCZHlrW01ly9epUff/yRsLAw2rZty9tvv02LFi1Yv349M2bMoFWrVhQuXDhXJwDXlAHWYuwh0AYjMRD/lvfTXSGEAGw2G8OGDWPmzJlMnTqVYcOG4e3t/a/j/vrrL9auXcuPP/5Is2bNiIqKomjRoiZELITIKVarlX79+tG6dWsmTZrEgw8+SI8ePXL1v32tNRcvXiQuLo4NGzYQFxfHfffdR7t27Rg9ejSenqYNkskRjwGLU/5/PNAVmqH1WpTK3M6gd8KaeBQ8u+RYe5kkSYAQIt/w8PBg0KBBzJ8/n8mTJzNy5EgKFiwIwPHjx3n33XfZtWsXrVq1Yvbs2RQrVszkiIUQOaly5cosXryYlStX0r9/f4KDg2nYsKFbT4i9ntaas2fP8ttvv7Fx40YOHTqEv78/bdq0YdSoUfluWdReQBeMFYPegxeBWcBvORaAxaM4LtUH+D7H2swEG/yzvb3W2ux4RC62Y8eOPNGNKPI2T09P+vXrR2RkJOPGjSMoKIj33nuPvXv38swzz9CrVy8Z8iNEPubl5UWPHj1o1qwZM2bMYNOmTfTt25e77rrL7NDSpLXm1KlT/Pzzz2zZsoUzZ86kPvGvXbt2mj2e+YkX4A3Mhml9oBhafwo8h1KObG/cbvXBSsNsb+c22QDq1avHjh07OHHihNnxiFysaNGiVKtWzewwhLglDw8P+vTpw5w5c+jXrx99+vRh4MCB+PrKqtJCCMP999/PzJkz+fDDD+nbty89evSgZcuWbvOw68SJE/zwww/ExsZy9epVqlatSmBgIFWrVk3t4RSGUcB6YzPhz4FxQGG0TsrRoUFuyAbMa968uWfz5s3NjkXkHb+aHYAQt2Kz2RgwYAB9+vTJ90/KhBBp8/LyIjAwkEaNGjFz5ky++uorXn/9dUqWLGlKPCdPnmTr1q1s3LgRh8NBrVq16N27NxUrVpQb/1t4Es5gzBPeA2wENgAjsrVRa+JRlFe3bG3jDtiUUv3NDkIIkascAfzMDiIr2Gy2fLEcYC5x6yWbhDBJ2bJliYiIYMOGDfTt25euXbvStm3bHJkrcObMGTZs2MCXX36JUoqGDRsydOhQypQpIw8wMkupb9G6L3AWGI/WrwHLUOpqtrRn8SiOy9IT2JYt579D8u0nhMgUpVRHs2MQQoic5uHhQbt27XjkkUeYNm0an332GePHj+fuu+/O8rbOnDnDZ599xqefforVaqVNmzZMnToVPz+/TO1zItKg1Hy0XgxMx5gsHAvsz5a27FYfrLpxtpw7C0gSIIQQQgiRQSVKlGDatGl88803dOvWjVdffZXnnnvujm7Ok5OTOXXqFBs2bGDjxo04nU6eeuopoqKi3HZCcq6mVDJanwRaAveh9SygHUrlrQ0ibkGSACGEEEKITHr00UepXbs2ERERxMbGMnz4cCpWrIhSKkP1ExISOH78OLGxsWzbZowWeeKJJ5g9e3a29C6Imyg1Ba0rA+uAwYAfWl/O0qFB1sSjUKBHlp0vi0kSIIQQQghxGwoVKsTEiRPZtm0bo0eP5vHHHycwMDDdSbpXr17l4MGDbN68mZ07d2K1WmnVqhXTp083bbJxvqbUPrRujjHXbQvwHsbqQVnD4lEczavAt1l2ziwkSYAQQgghxB1o1KgR1atXZ9GiRQwePJi+fftSs2ZNlFIkJCQQFxfHpk2b2LdvH76+vrRo0YKOHTtSsmTJDPcciGyi1PdoPRhjXkAEWg8BIrOkR8DYJ+DROz5PNpEkQAghhBDiDhUpUoTXX3+d7du3M3v2bKpUqUJSUhIHDhzAz8+Pxx57jC5dusiNvztSaiZazwMigScwhggdNjWmHCBJgBBCCCFEFrBYLNSvX58qVaqwePFiGjRoQHBwMH5+eWJV5bxNKTta7wemAPXQegHwFErp2z6nb4FjJCX0zKoQs5okAUIIIYQQWahQoUIMGTLE7DBEZikVgdY1gVVAT+BetD5320OD4hOLougMfJOFUWYZSQKEEEIIYbaTwHazgxB52ukMHaXUb2jdEDgHbAUWAVNvq0WNLxbV7Lbq5gBJAoQQQghhKqXUJ8AnZschBABK/YLWIcB3wEy0HgnMzradhU0iSYAQQgghhBDXU2o6WnsAK4GHgBVonZCpOQK+BY6RmBiUXSHeKYvZAQghhBBCCOF2lLIDPwHNgVbAp5mqH59YFK1eyvrAsoYkAUIIIYQQQqRFqelASYzlQxejdUW09s5QXY0v6BbZGd6dkCRACCGEEEKI9Cj1E8aQoN0Yu/+67RCfzJAkQAghhBBCiP+iVBzQAfgMiELr8Fv2CPgWOIaFPjkR3u2QJEAIIYQQQohbUWoy0A9YAwQAvmid/r10fGJRXHTIoegyTZIAIYQQQgghMkKpZGAT8BjQEfgw3WM1vkDLnAks8yQJEEIIIYQQIqOUmgHcD0wH5qJ15QxPFnYjkgQIIYQQQgiRGUr9H1AHOIUxWbjzv47RHEfr13I4sgyTJEAIIYQQQojMUmof8DywAliF1pNv6BEoYC0MlhfMCu9WZMdgIYQQQgghbodS49DaC/gCKAxMR+tklHKS5CyElVYmR5gu6QkQQgghhBDidimVhDFBuBUQDMSYG1DGSBIghBBCCCHEnVBqNsYcgdHADLSuMbfGk+fQ9DM5snRJEiCEEEIIIcSdUmoT8CCQDHwbXanRSyieMTmqdEkSIIQQQgghRFZQ6g+gHTBrwvfvbTvtXaSTuy4fKkmAEEIIIYQQWUWpsUDE4EZd30q0ehQAbGjtYXZYN5PVgYQQQgghhMhKSiWcufL3cqU4BIzA2Fws0OSobiBJgBB5hF4YWAOrrkyv6HVKoc2ORwiRd+jwcAv37mmA1kVUUPQGs+MRIjc49n7f938tXGYB0BRojdb1gDiUSjQ5NECSACFyLT0/wBcPSzvQz6N1VXCVwmltIAmAECIr6KiXSoCrPcrVDnZXwqWtQD2z4xIi10hyFqp17s+qQG3AH/gUeBVYY2ZY10gSIEQuot8KuB+X6gWuQLQujnb5YMztOYNWX2J1DNdRAVOAl1G6PAB222Q8HK8A5VLLNldnlKssAMo1CZety3+XrV1RugwATjURC6/+Z1npbli4FwCrbQJ2R/fUsnaNB1tPlOueDJXtahweBIEunaGy8gjHZe+NolSGyk5nGBaP11CukmmWbWosyfTDou/OUFl7hIJ9AAq/DJWTnWPw9BgErhIAWNRonAxB6bsyVHZ4jsKaPBRF8QyVPRwjSfIchsVVLM2yViOAEShdNEPlRK/hFEgaBRTJUNmaPAxHgTEoV+F/yl6hKF0o5fcZgo2x/1m26jAs+AJQwHuo6vJOPCJL6IXt62C1BuPSz4KjAFAYjQVcZ9CWN1HO0npRh6nGwZY94PwSpQYatdVeHHyBTQ8yfs4+lPVzcA5Ou2zZj3Z+ikUNSbPsUr9j42Nc+vWU+gfA+iHKOTTdslN/gE2FGPFYDuJwrruhrPRatB5mHK/+wMKa9MscwmWLxuoYbsRjOYzN8R4uy4g0y1r9CaxG6ZFpluEIdttKPByjUt7vUSyOFWjL6LTL6hgulmHVY9IsuziOh20JTkdoSv2/wPEWyjI2zTLqBHYW4aHD0ixrTmLxiETbw9Mtu5wLsVrHpZz/FB56Pg49PuX3dRpP5t1QVra5KPuElPpnwGNOumUsZ0m2z8LTOjHl/f6NlRm49KS0y5zD6fkGtuTJaZZdlvN4JU/DbpuSUv8CMBWlp6ZZhoskek2mQFJEmmVtuYQtcSJOz2kp9S/jYDweenqaZRdXwLUILKDUEbTuDYQBW9B6NjDc7B4BSQKEyE2sPsdwXfkGrRqhABcOFMWML2l9Fq32ge0q2nEA1EUAbCoBl/odxYXUstb7QZ0HQHsloh3/lJN1Ah6u39HXl/V+tDoHgNORiMX6T9kD4/ya68qW/Wj9t3F+ZyKW68pWaxIu9qHV2dSyQ++HlLKzQCIW+z9l5ZGEy74Ppc6klknej04pF7QnYrf9U/ZISMJu24dWp9Msa0ci2rIftFG22JLQ7AN1Ks2ynWRQe9GcBOCKZxJeyfvQKu1ykiUZL7UXzQnj845Pwum5D63SLns7k3Fa96KVcVOskuwozz3plpNIxqr2AsZNsbbbjZsxnXbZ15pMotqLTrlpthayY0ncg1Zpl70LOEhM3H1D+WrCHlA+xl/CRDva65+y1yUHLq/dqOvK2nsP6IIAOJKvxWOUL3o58GE3Wl1XVnFoCqa8Pye268oFfRzY43ejlfd15Tg0RtlpdSKyjqf9AC7r1yjqoXVxLAo0RUGBcl3Aoq/gsvwKgOYY2nIZCylldQylr6DVP2WXKx5LSllxHKu6guNaWR/H5hGPw5lyvOsvo37K+dF/4XTFgzXl5/ovlP0q+rqy0xWPLeV8qBNofRWdEg+cwMMzHpfjWnsnU66P1+I5icMRj8V27XyncDhvLJN0FZ1SxnUKu05IfT8WTmOxXsWZ0p5Sp7HpBOzX3j9nUNe1pzljXH+v/VyfweWVCI5r5z9Lkk7A49r70WdxOhNT41Hqbzx0As5r7au/0c7E1PMp/sZiTUqNR+tzuAokYrVfa/8cyiMJbb/2eZynoCORZI9/yraEpBvKLkcipP4+L2BRSamfr+ZCyvXx2ud/gSseSXglX4vnYsr18Nrv48ay0pewxifh8Eopuy7h7UzGYfvn5yrZjva61v5lknRy6u9bq8vG9S71/V/B15pMwrXPR1/BWshOcuI/Py9QwEFC4rXPJx6VZEcX+KfsdcmRWkZdNa5fqfFfTbl+/fPzJO1MjUeRQEEfB8nx1/5+JODhe5Skq+Epv79RaF0I+AY4DSi09kSpZEyizGpYiJxSrFixA7GxsffVrVvX1DiGDh1qf/PNN8cA07LqnHpxQEU0nUE1R+vaeFjbqe7vbcuq8wuRTzxWv379tT/88ENRswNRSmnACu4xrE9HBdQE1Q30I2juw0Zj1TPmgNlxCZFraR0EvINxL1AUpTqbFYr0BAiRi6leMYeA8cB4vfTVAiRfflBrlMwLEEJkBRUc8xswBEC//Uwh7B6VTQ5JiNxNqUVo/QIQADyG1o2An80YGiRJgBB5hOq2LBH4nmCzIxFC5EWqx8eXgZ/MjkOIXE+pdWj9P6AuxiThDsBnOR2GbBYmhBBCCCFETlLqBMbSoQOBX9B6AVoXyMkQpCdACCGEEEKInKbUcLQuCnwH7AFAay+USsqJ5iUJEEIIIYQQwhwXgQjgXWAh4AJ65ETDMhxICCGEEEIIMyilUWop8DzQCpiM1s1zYmiQJAFCCCGEEEKY6z2MHbnrAl8ATbK7QUkChBBCCCGEMJPRI3AaqI8xHOgAWi9Ca6/salLmBAgh3IpeHPgITt0s9QXl0ih1AsW3KfsiZH2b4c1tlPazgxqogqPnZEcbQgghxC0pNQytSwA/Av8HZNtkYekJEEK4F+1sgdJTUfpllA4A1RnNAlz8rqMCRpsdnhBCCJHN/gZCMXoEVgKzsqMRSQKEEO7JZn1SBcc8pIJjahLP3cBHwAT91gvlzQ5NCCGEyDbG0KCVwEtAHWASWj+R1UODZDiQEMLtqSExCTqywyqUegGHtRLwJ4Ce/2JZPFyBaBqg9V0o9Sear1Tv6HdvPoeObP8AFktvNA8CdlB70Xqx6h2zO712dWRABxT1gIUqOOaIXtSxF2g7uH5EW19H6/Io9QdO5xT12to/UutFdXwOxf0o6zu4nCFoXQeLWqWCopfo+QG+2OgHNEPjiVI/YbPMVj3e+yu1/qKAZ3Cpqjj1cqyMAGqi1GlceqbqE/NDFn60Qggh3Ndy4GPgWWAB0Br4NqtOLj0BQojcQakWgB2rc0/qazZnF7R+ATiCYhPa5YPSK3Vkx8nXV9WLAx5DWXageRatf0HzA+jKWNRL6TWnozqGoliF1ntUcMwRAFz6RTRD0ZbNoONR/A90K6yW7/XCFytcV/txtB6My/F/oB8A/TMukvWKzj7Y+AYYAxxAqZ9Av4LD+bNe0N7/n+q6JUoPwca3oEpg4WsjXr7RkR2bZsXHKYQQws0ZPQLngOpAIHAGrZegtWdWnF56AoQQ7snujNJRHRPQ2oaiEuADvKp6rTuWeoynz5uq27JJ11fTUR1HoXSojnp6ogr+5KqODvDkPMuBvThopvquuZJ6bHjzf10DjUnCdy8E3RFNO9V7zVc3HVId1PMqOPpDAL34pbm4HHtRronAK9cddw/owSp4TepYzpQ5DbXRqpnqHb3VqP/CXFzW3VitU4GO19W/G63nqt4xEwH0jIDp+PA9Ss8BamfikxRCCJGbKTUUrUsD2zF6BixZMVlYkgAhhHtS+gJwFaVtoE4AjVDqER3ePFqFb3EAqG7LEnV0gJUL1ARdBo03yqXRqgCugv7ALi7QECgDesj1CQDAtfNcx5fSfh+DfhBtaap6v78zjcgOExT9EcEp5+j17im9qMN7GN2110ummFpw02vNUfyggo0EwKi/7piO7PAeihduOlaT7JW6UpEaEpOgozpGgZ6n5z1/l+r3wd//9fEJIe6MnhHgja9qfMOLTnUSl+uw6htzJZ1qd97uooCxuCiuescMyq42RK50EhiMMT9uHbAfuKO/I5IECCHck802/IZx8pEdWgKxlL77EDATQEd1bMF5vRgoA+pP4DIabwCUsySwC5cuj1KgdLpj//+hQwFPXK6HVJ+YtBIA0PyhFPqm1w4ApfSMAG81JCYh5dXjqmNM8k21/XGp//3rnBbLfrS+S0cFFFHBMRdTXj2jBqy6dMNxigNowGapiLF6hBAiuxRwlUZbbuwJtGiwkKCjAuYRFDP8X9eCrKCphYVSWX5ekbsppYH30ToYKAe8itbPABtut0dA5gQIIXIF1XvNRuAsSj/5z6t6KYp9JHnerYJjKqvgmIdQuhcAFosy/otxI+60+GSgmbnAASyW5XrxSyXTDiQlybiBpSBgpwzX3/Tb06gdj9Jp1NcFASfFSLjuxfSOA6stPu3whRDZYAxOj+LGH9d9wDIghEUdu5gcl8ifFgGNMeYIvAvUvN0TSU+AECJX0FFPFwQK4+IKgH4roDhOygOhNzwx1zS8sab1Z3CBRTcHbrGyjjqGQzfDRiwuxxa9MKCV6hNz/KaDqun5Ab43DAfQugGKg6pjjPMWb2M/0ECHN7fdMBTJRSMUf97Uc1BILw6sqnq9/89EaPTDoOzYvLNl0zQhRBo0Ceq11edTSud1VNBAON8T7WqMsXqLcdhbL5THZWuJS5cHLqL1z/Re83VavQU6sv1DYG2J0oVAHwW1PnUBgrRCWNHZh4SEtqBOqOCYb/WC9v5YVT1OnP2A0iUeQ6vGKK4A719/Hr0kwA+7bo7TMxaPpLJo9RTakqiCo1N6UwPKgXoOdGngENjWqeB3z6bWj3qpBNhbYFUbcVoqgH4K7dIo9ZEKjvntTj9acRuMHoGLaF0eYxiqHa2XAz1RKq2HT+mSngAhhHty6XJ6QXt/vaC9v44MaAgF3gM8UfpzAApzEbgC+slrE3yNL1Y17PrTqKD396P5FNRIvTjgsWuv66iXSuiF7evc3KzqG3MSbM2BBCxsvXHVn5SWbUzT4eEWAL0o4HkUbYEVt3xPWi0HylLq7vFao4yYO76C4nFQy24+Gu2arVd09gHQiwPqoVUw6PdVt2WJt2xLCJFNzpcGrs1VAkAvaO+P03oIrV/HwkMoXsKiNrOo4+pr/9YBdHi4RUcGLEVZfkTpbkA9UEOAzem1phe/VJKExC2gwgHjBt9qaQUqmtJ+80EtQ1EfeB2I01EBTVIrJ+nqoKKxOCbisvwfWj0P+lkAHdnhBWAf6OHAQ8BkcOw3rrcpXM7KoKJxqghwfQW6CUp1A37RkR2kJ8RMSg3FeLC0ETgNeGZ21SBJAoQQ7snl+h9Wy0GsloMotgENUAwiaM3bAMZTdz0UrQIo7XdERwXsQ1k2kTJf4AYOx6tofsbFRh0V8JeOCvgdHH9htTydVtMq+N2zeNofA85gcX6tI9s/cN2PN6KoSum40zoq4CCadSg+oRgzbvWWVO/oT0CNQ+lhLAo4q6MCjqH0Oyj9HsV0xE2H/4lWf5KQeFJHBezDxQ/AIfAYnMFPUAiRFZR6REcFBOmogD46KmAisAXYgdX2z8R/L/sZlKWKCo6proJinlLBMfVBPQ/6RRZ3fDz1uNK7B6J4FXhNBcdUUcExT6rgmMo4aZtW03pBQGVcjv8B8Tg9mqTRW/AgtqTKKjimLVAZ46ZwmY4OsN74HvQTQHUVHNNABcc011EvlUCpZcC3xHO/Co5pDbYqwHEU7+qoII+bImmDw1pbBcc8wYkzlYEPUWqBXvj83Zn9OEWWOgr0BMKBz4CJmaksw4GEEO7FrpbiyZepZYfLhctywnhCDwT9c6gKXhOlFwVuBtfDoO0oj81cdlyiEF+R4PF76nHGSjotU3oKaoLFjnbGqT5rfwEgbIuTtwIewmE/mlqn24cX9IrOLUlKrAKWf4buaBIpxnOc082xWPxRrr30XPPtTV3+U1Gu+Wm9PRUcHa4XvrgM5WqCRdtwWX5RwWmuQoQKju6lo9q/hbbUwKL/oqjaqDq+e/NkYyFEttKPAfVSCj5AYbRei83r2iR+VI+PLwOXAfTSVwuQmFgKp/1XrJaTaP0QsCHlXEHAVhUcs/D6FtRrMfv+1erCjo2x6I+AWJKudFUD1v978qdW41PaRgXHXNSLAiaiWcc5Vx3gp+uOnHZDAqGdbVAUAstwNeT9BKP+u2f1oo4T0Pp9XBcaAP93Xf05qu97R8FYVU0vChyFdrXHamuLMUdCmMEYGvQhWg8AfIEpaP3/7d1fbFN1FMDx77kdf7KxOJb4B8OIxslwb5o0AY0PEo380SjZuqFEIg+0JRgxJAhqIouJxCcfMCAt8f+fjHV7gQej20SNUQFJkEiUmTkSJUBsMpawsX+9x4fbdbtdXathdMzzSZqsd+d3T7uH9f5uf+f8GvAmaXk/K2wSYIyZUdIX+xcLjg8f6sK7+zXRyZyx0bYf8X8wescFhcSkMbLxo/5c50qv/e9MP2Bz1u+nWNsLIFuazwHnpooZP1fbMeBYIbHGmGmg7JFoIvNNn8Ybg+C2MzywGHgKQFtCAXplJ2iU4f4qb51FerGF6GLI7EtSAxwpIGcNjnaAfMyF2og0Nbm549R/A0E4hQKOVOP/X+fvjibuXSApuMm/rl/lJCgE9G58kwB/Hgkf6tJYaAClOu97MdfDW8C7wDZgJ/ArcDrfIFsOZIwxxhhTIAkfOgGSAG0cq9mhV18A3Q36OqOBJTBYJpGEAN24eEtzFi0VQFDJ3p8kl4vAz6CPsOiXO/8xykn5JwcuXnMC18m6yesM+N+EUwK4LOzNGi/p8fiXEym5mh6kvH1cTNF5OwtfASqA1UApqhvyDbNJgDHGFMqRZhzapj2PSCcqB6Y9jzHmP3IdQLiSShdiykMoRyXSGpOtzX9I5MiA7gstwOvn7kVE4iNAD6L5d/wW+oCHEf4E92uN1S3LHVjiP67Upn/onjqBdgNz6HOW+vO6XrtJzRovco9v9L71VUA5OHnymOtKZAfeBHId8Gm+cJsEGGNMgSTcclDCifenP0/isERbsguFjTEzgMbr7wcJAWcmtA4dRFiS6VSmCCX6GpBdYPsBsErj9Y/5ju5/emF2Hokk+hhhFdAFzlcaC03uB6+6c6yIV/eunofwInCBytxLIjNStAPDuPrqWBGxxh4vRXQX8BdXJXsJ4nP63pMV4+/NfQWvNeXnU+YxxfA7sCtdLzAl+xrHzHqBQCC5adOm28rKyvL1cJ9WPT09c/HaeBljjLlRCC9rLLQ1/awSpQI4i7rPZmJU30bkCxbdfELj9T8Ql3uBAN4F2bih/jeYu2AFyGGNhb4E6QG9A0aW4NUL+FNvTVzRN0NrKaMNOKoHQ4/K5gn1S8IQ9J7UWMMx0BVADUJdjt3K/efdkjivBxq2IbqfXpZprP4nkAeAKkTrZXvr1awhZxmec0rjoQ7i1IKuQGSHhFumrH8yRVDAxf8YmwSYWS+ZTK5JJpP3Fft1AMPAN8V+EcYYYwowOj9JYDjiO6Y6iDrdXKr9fmKxrkRbOzXeuBzVRlQqUfcTZOgdmLcWnEuZuOc/G9KmpjXcfuYJlJXglqPyHbgvjecgNnFnctmeuKp7V69jfvkGUlqNv1nBRkTqcPVBRNpxU89kup4BzC3pYjQVAWfSxbpEWw7owdAJ3NVrMgAAALZJREFUXBqBWxBpxk19KJG237Jjwd2DOLeirAU5jepuibS0/7s/qJlpJH+IMcYYM6utDAaDbcePH68o9gsR7y5eACbvMmvMGI2FwkCMUcp9u5df6zxem9JvwV2e7lRmZhGrCTDGGGOMMeZ/xiYBxhhjjDE3lvNAB0M5W3deO4HUZaAD6MsXaowxxhhzo1kZDAZ7dQYAXGyprjHGGGOMMdOuqrS09DLeOvyiPiorK/Pu8mmMMdfC3+73wX+/neAOAAAAAElFTkSuQmCC"
    }
   },
   "cell_type": "markdown",
   "id": "4dd4686c",
   "metadata": {},
   "source": [
    "![ART_Cert_Training.png](attachment:ART_Cert_Training.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8229684f-9136-4e08-8743-a76600a3880f",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/giulio/Documents/Projects/certified_training_art/venv/lib64/python3.8/site-packages/tqdm/auto.py:22: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html\n",
      "  from .autonotebook import tqdm as notebook_tqdm\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "\n",
    "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 deep_z\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",
    "\n",
    "if not os.path.isdir('notebook_models/mnist/certified/'):\n",
    "    os.makedirs('notebook_models/mnist/certified/') \n",
    "\n",
    "if not os.path.isdir('notebook_models/mnist/pgd/'):\n",
    "    os.makedirs('notebook_models/mnist/pgd/')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f92b485c-7010-4315-a5d8-c9a59b08d7c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# For all of the demonstrations we will use this example MNIST model.\n",
    "\n",
    "class MNISTModel(nn.Module):\n",
    "    def __init__(self):\n",
    "        super(MNISTModel, self).__init__()\n",
    "        self.conv1 = nn.Conv2d(in_channels=1,\n",
    "                               out_channels=16,\n",
    "                               kernel_size=(4, 4),\n",
    "                               dilation=(1, 1),\n",
    "                               padding=(0, 0),\n",
    "                               stride=(2, 2))\n",
    "        self.conv2 = nn.Conv2d(in_channels=16,\n",
    "                               out_channels=32,\n",
    "                               dilation=(1, 1),\n",
    "                               padding=(0, 0),\n",
    "                               kernel_size=(4, 4),\n",
    "                               stride=(2, 2))\n",
    "        self.fc1 = nn.Linear(in_features=800,\n",
    "                             out_features=1000)\n",
    "        self.fc2 = nn.Linear(in_features=1000,\n",
    "                             out_features=10)\n",
    "        self.relu = nn.ReLU()\n",
    "\n",
    "\n",
    "    def forward(self, x):\n",
    "        if isinstance(x, np.ndarray):\n",
    "            x = torch.from_numpy(x).float().to(device)\n",
    "        x = self.relu(self.conv1(x))\n",
    "        x = self.relu(self.conv2(x))\n",
    "        x = torch.flatten(x, 1)\n",
    "        x = self.relu(self.fc1(x))\n",
    "        x = self.fc2(x)\n",
    "        return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c743dd9b-7800-4777-9bdd-f0298db2942f",
   "metadata": {},
   "outputs": [],
   "source": [
    "model = MNISTModel()\n",
    "model = 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": "ab5146d2-8b04-4968-9d9f-409f2bf3cba9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "End of epoch 0 loss 0.4505019783973694\n",
      "End of epoch 1 loss 0.14715342223644257\n",
      "End of epoch 2 loss 0.08925174176692963\n",
      "End of epoch 3 loss 0.06519503891468048\n",
      "End of epoch 4 loss 0.052043616771698\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().detach().numpy())\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": "0ac5fd81-fb23-4d00-a550-caa38fd1ebbb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test acc:  98.36\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": "2d4a4cc5-5ca0-41eb-982a-3e9feb055661",
   "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.linear.Linear'>\n",
      "registered <class 'torch.nn.modules.activation.ReLU'>\n",
      "registered <class 'torch.nn.modules.linear.Linear'>\n",
      "Inferred reshape on op num 4\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/giulio/Documents/Projects/certified_training_art/adversarial-robustness-toolbox/art/estimators/certification/deep_z/pytorch.py:239: UserWarning: \n",
      "This estimator does not support networks which have dense layers before convolutional. We currently infer a reshape when a neural network goes from convolutional layers to dense layers. If your use case does not fall into this pattern then consider directly building a certifier network with the custom layers found in art.estimators.certification.deepz.deep_z.py\n",
      "\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "# But how robust are these predictions? \n",
    "# We can now examine this neural network's certified robustness. \n",
    "# We pass it into PytorchDeepZ. 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 PytorchDeepZ 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.deepz.deep_z.py \n",
    "# can be used to directly build a certified model structure.\n",
    "\n",
    "zonotope_model = deep_z.PytorchDeepZ(model=model, \n",
    "                                     clip_values=(0, 1),\n",
    "                                     optimizer = optim.Adam(model.parameters(), lr=1e-4),\n",
    "                                     loss=nn.CrossEntropyLoss(), \n",
    "                                     input_shape=(1, 28, 28), \n",
    "                                     nb_classes=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "8b7b0b1a-76e8-40ef-b919-7a686eef8ef5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lets now see how robust our model is!\n",
    "# First we need to define what bound we need to check. \n",
    "# Here let's check for L infinity robustness with small bound of 0.15\n",
    "\n",
    "# lets now loop over the data to check its certified robustness:\n",
    "# we need to consider a single sample at a time as due to memory and compute footprints batching is not supported.\n",
    "# In this demo we will look at the first 50 samples of the MNIST test data.\n",
    "\n",
    "original_x = np.copy(x_test)\n",
    "def certification_loop(model, x, y, preds, bound):\n",
    "    num_certified = 0\n",
    "    num_correct = 0\n",
    "    for i, (sample, pred, label) in enumerate(zip(x[:50], preds[:50], y[:50])):\n",
    "\n",
    "        # we make the matrix representing the allowable perturbations. \n",
    "        # we have 28*28 features and each one can be manipulated independently requiring a different row.\n",
    "        # hence a 784*784 matrix.\n",
    "        eps_bound = np.eye(784) * bound\n",
    "\n",
    "        # we then need to adjust the raw data with the eps bounds to take into account\n",
    "        # the allowable range of 0 - 1 for pixel data.\n",
    "        # We provide a simple function to do this preprocessing for image data.\n",
    "        # However if your use case is not supported then a custom pre-processor function will need to be written.\n",
    "        sample, eps_bound = model.pre_process(cent=sample, \n",
    "                                              eps=eps_bound)\n",
    "        sample = np.expand_dims(sample, axis=0)\n",
    "\n",
    "        # We pass the data sample and the eps bound to the certifier along with the prediction that was made\n",
    "        # for the datapoint. \n",
    "        # A boolean is returned signifying if it can have its class changed under the given bound.\n",
    "        is_certified = model.certify(cent=sample,\n",
    "                                     eps=eps_bound,\n",
    "                                     prediction=pred)\n",
    "\n",
    "        if pred == label:\n",
    "            num_correct +=1\n",
    "            if is_certified:\n",
    "                num_certified +=1 \n",
    "\n",
    "        print('Classified Correct {}/{} and also certified {}/{}'.format(num_correct, i+1, num_certified, i+1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "26d0ea51",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Classified Correct 1/1 and also certified 0/1\n",
      "Classified Correct 2/2 and also certified 0/2\n",
      "Classified Correct 3/3 and also certified 0/3\n",
      "Classified Correct 4/4 and also certified 0/4\n",
      "Classified Correct 5/5 and also certified 0/5\n",
      "Classified Correct 6/6 and also certified 0/6\n",
      "Classified Correct 7/7 and also certified 0/7\n",
      "Classified Correct 8/8 and also certified 0/8\n",
      "Classified Correct 9/9 and also certified 0/9\n",
      "Classified Correct 10/10 and also certified 0/10\n",
      "Classified Correct 11/11 and also certified 0/11\n",
      "Classified Correct 12/12 and also certified 0/12\n",
      "Classified Correct 13/13 and also certified 0/13\n",
      "Classified Correct 14/14 and also certified 0/14\n",
      "Classified Correct 15/15 and also certified 0/15\n",
      "Classified Correct 16/16 and also certified 0/16\n",
      "Classified Correct 17/17 and also certified 0/17\n",
      "Classified Correct 18/18 and also certified 0/18\n",
      "Classified Correct 19/19 and also certified 0/19\n",
      "Classified Correct 20/20 and also certified 0/20\n",
      "Classified Correct 21/21 and also certified 0/21\n",
      "Classified Correct 22/22 and also certified 0/22\n",
      "Classified Correct 23/23 and also certified 0/23\n",
      "Classified Correct 24/24 and also certified 0/24\n",
      "Classified Correct 25/25 and also certified 0/25\n",
      "Classified Correct 26/26 and also certified 0/26\n",
      "Classified Correct 27/27 and also certified 0/27\n",
      "Classified Correct 28/28 and also certified 0/28\n",
      "Classified Correct 29/29 and also certified 0/29\n",
      "Classified Correct 30/30 and also certified 0/30\n",
      "Classified Correct 31/31 and also certified 0/31\n",
      "Classified Correct 32/32 and also certified 0/32\n",
      "Classified Correct 33/33 and also certified 0/33\n",
      "Classified Correct 34/34 and also certified 0/34\n",
      "Classified Correct 35/35 and also certified 0/35\n",
      "Classified Correct 36/36 and also certified 0/36\n",
      "Classified Correct 37/37 and also certified 0/37\n",
      "Classified Correct 38/38 and also certified 0/38\n",
      "Classified Correct 39/39 and also certified 0/39\n",
      "Classified Correct 40/40 and also certified 0/40\n",
      "Classified Correct 41/41 and also certified 0/41\n",
      "Classified Correct 42/42 and also certified 0/42\n",
      "Classified Correct 43/43 and also certified 0/43\n",
      "Classified Correct 44/44 and also certified 0/44\n",
      "Classified Correct 45/45 and also certified 0/45\n",
      "Classified Correct 46/46 and also certified 0/46\n",
      "Classified Correct 47/47 and also certified 0/47\n",
      "Classified Correct 48/48 and also certified 0/48\n",
      "Classified Correct 49/49 and also certified 0/49\n",
      "Classified Correct 50/50 and also certified 0/50\n"
     ]
    }
   ],
   "source": [
    "# We can toggle how zonotope_model will process data though the method set_forward_mode.\n",
    "# with 'abstract' it will use abstract operations on the input (which is expected to be a zonotope).\n",
    "# with 'concrete' the nerual network will be run normally and expects regular data as input.\n",
    "\n",
    "# As we want to do certification analysis now, lets set the model to run in abstract mode\n",
    "\n",
    "zonotope_model.model.set_forward_mode('abstract')\n",
    "certification_loop(model=zonotope_model,\n",
    "                   x=np.copy(x_test),\n",
    "                   y=y_test,\n",
    "                   preds=test_preds,\n",
    "                   bound=0.15)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92f68aea",
   "metadata": {},
   "source": [
    "We can see that our NN has a very low certified robustness. You can modify the `bound` parameter in the function `certification_loop` to see how the certfied robustness varies. We can now try and improve this through different robust training stratgies. We will look at:\n",
    "\n",
    "+ Using certified adversarial training\n",
    "+ Using PGD adversarial training."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "9e25ac62-b575-468d-99d6-55d069e16731",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<All keys matched successfully>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from art.defences.trainer import AdversarialTrainerCertifiedPytorch\n",
    "\n",
    "# We will now train the model to improve its certified accuracy. \n",
    "# Regular PGD training will boost certified performance, however even higher certification scores can \n",
    "# be obtained by training the nerual network with the objective of certified performance. Let's compare the \n",
    "# two methods.\n",
    "\n",
    "# NB! Certified Adversarial training takes about 9 hours on an NVIDIA V100 with the following parameters.\n",
    "\n",
    "pgd_params = {\"eps\": 0.3,\n",
    "              \"eps_step\": 0.05,\n",
    "              \"max_iter\": 20,\n",
    "              \"num_random_init\": 1,\n",
    "              \"batch_size\": 32,}\n",
    "\n",
    "trainer = AdversarialTrainerCertifiedPytorch(zonotope_model,\n",
    "                                             pgd_params=pgd_params,\n",
    "                                             batch_size=10,\n",
    "                                             bound=0.15)\n",
    "\n",
    "# Uncomment if you wish to train your own model, but it will take several hours!\n",
    "\n",
    "#trainer.fit(x_train,\n",
    "#            y_train,\n",
    "#            nb_epochs=30)\n",
    "# torch.save(trainer._classifier.model.state_dict(), 'notebook_models/mnist/certified/model.pt')\n",
    "\n",
    "# Here we will load a model which was trained using the above function.\n",
    "trainer._classifier.model.load_state_dict(torch.load('notebook_models/mnist/certified/model.pt', map_location=torch.device(device)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "e7855a8c-0ab0-4ee6-a7b9-380d6c2360cc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test acc:  98.7\n"
     ]
    }
   ],
   "source": [
    "# Like before we obtain the model's test time accuracy\n",
    "# Make sure to set the model in concrete forward mode\n",
    "\n",
    "with torch.no_grad():\n",
    "    trainer._classifier.model.set_forward_mode('concrete')\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": 11,
   "id": "65c3e186",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Classified Correct 1/1 and also certified 1/1\n",
      "Classified Correct 2/2 and also certified 2/2\n",
      "Classified Correct 3/3 and also certified 3/3\n",
      "Classified Correct 4/4 and also certified 4/4\n",
      "Classified Correct 5/5 and also certified 5/5\n",
      "Classified Correct 6/6 and also certified 6/6\n",
      "Classified Correct 7/7 and also certified 7/7\n",
      "Classified Correct 8/8 and also certified 8/8\n",
      "Classified Correct 9/9 and also certified 8/9\n",
      "Classified Correct 10/10 and also certified 9/10\n",
      "Classified Correct 11/11 and also certified 10/11\n",
      "Classified Correct 12/12 and also certified 11/12\n",
      "Classified Correct 13/13 and also certified 12/13\n",
      "Classified Correct 14/14 and also certified 13/14\n",
      "Classified Correct 15/15 and also certified 14/15\n",
      "Classified Correct 16/16 and also certified 15/16\n",
      "Classified Correct 17/17 and also certified 16/17\n",
      "Classified Correct 18/18 and also certified 17/18\n",
      "Classified Correct 19/19 and also certified 17/19\n",
      "Classified Correct 20/20 and also certified 18/20\n",
      "Classified Correct 21/21 and also certified 19/21\n",
      "Classified Correct 22/22 and also certified 20/22\n",
      "Classified Correct 23/23 and also certified 21/23\n",
      "Classified Correct 24/24 and also certified 22/24\n",
      "Classified Correct 25/25 and also certified 22/25\n",
      "Classified Correct 26/26 and also certified 23/26\n",
      "Classified Correct 27/27 and also certified 24/27\n",
      "Classified Correct 28/28 and also certified 25/28\n",
      "Classified Correct 29/29 and also certified 26/29\n",
      "Classified Correct 30/30 and also certified 27/30\n",
      "Classified Correct 31/31 and also certified 28/31\n",
      "Classified Correct 32/32 and also certified 29/32\n",
      "Classified Correct 33/33 and also certified 30/33\n",
      "Classified Correct 34/34 and also certified 30/34\n",
      "Classified Correct 35/35 and also certified 31/35\n",
      "Classified Correct 36/36 and also certified 32/36\n",
      "Classified Correct 37/37 and also certified 33/37\n",
      "Classified Correct 38/38 and also certified 34/38\n",
      "Classified Correct 39/39 and also certified 35/39\n",
      "Classified Correct 40/40 and also certified 36/40\n",
      "Classified Correct 41/41 and also certified 37/41\n",
      "Classified Correct 42/42 and also certified 38/42\n",
      "Classified Correct 43/43 and also certified 39/43\n",
      "Classified Correct 44/44 and also certified 40/44\n",
      "Classified Correct 45/45 and also certified 41/45\n",
      "Classified Correct 46/46 and also certified 42/46\n",
      "Classified Correct 47/47 and also certified 43/47\n",
      "Classified Correct 48/48 and also certified 44/48\n",
      "Classified Correct 49/49 and also certified 45/49\n",
      "Classified Correct 50/50 and also certified 46/50\n"
     ]
    }
   ],
   "source": [
    "# then let's compare the certified robustness using certified training against the \n",
    "# model that was trained regularly. We will see that the model now has extremly high certified robustness.\n",
    "\n",
    "trainer._classifier.model.set_forward_mode('abstract')\n",
    "certification_loop(model=trainer._classifier,\n",
    "                   x=np.copy(x_test),\n",
    "                   y=y_test,\n",
    "                   preds=test_preds,\n",
    "                   bound=0.15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "b48fd5da",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<All keys matched successfully>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# For the final comparison, we will look at how PGD trainng compares. \n",
    "# Let's make a new model and train it using art's AdversarialTrainer \n",
    "\n",
    "from art.attacks.evasion.projected_gradient_descent.projected_gradient_descent import ProjectedGradientDescent\n",
    "from art.estimators.classification import PyTorchClassifier\n",
    "from art.defences.trainer import AdversarialTrainer\n",
    "\n",
    "pgd_params = {\"eps\": 0.3,\n",
    "              \"eps_step\": 0.01,\n",
    "              \"max_iter\": 30,\n",
    "              \"batch_size\": 32,\n",
    "              \"num_random_init\": 1}\n",
    "\n",
    "model = MNISTModel().to(device)\n",
    "optimizer = optim.Adam(model.parameters(), lr=1e-4)\n",
    "criterion = torch.nn.CrossEntropyLoss()\n",
    "\n",
    "classifier = PyTorchClassifier(\n",
    "    model=model,\n",
    "    clip_values=(0, 1),\n",
    "    loss=criterion,\n",
    "    optimizer=optimizer,\n",
    "    input_shape=(1, 28, 28),\n",
    "    nb_classes=10,\n",
    ")\n",
    "\n",
    "attack = ProjectedGradientDescent(\n",
    "    estimator=classifier,\n",
    "    eps=pgd_params[\"eps\"],\n",
    "    eps_step=pgd_params[\"eps_step\"],\n",
    "    max_iter=pgd_params[\"max_iter\"],\n",
    "    num_random_init=pgd_params[\"num_random_init\"],\n",
    ")\n",
    "\n",
    "trainer = AdversarialTrainer(classifier, attack, ratio=1.0)\n",
    "\n",
    "# Uncomment if you wish to train your own pdg model, but it may take 30min to ~2 hours depending \n",
    "# on your hardware.\n",
    "\n",
    "# trainer.fit(x_train, y_train, nb_epochs=20, batch_size=32)\n",
    "# torch.save(trainer._classifier.model.state_dict(), 'notebook_models/mnist/certified/pgd_model.pt')\n",
    "\n",
    "# Here we will load a model which was trained using the above function.\n",
    "trainer._classifier.model.load_state_dict(torch.load('notebook_models/mnist/pgd/pgd_model.pt', map_location=torch.device('cpu')))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "fddea1c9",
   "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.linear.Linear'>\n",
      "registered <class 'torch.nn.modules.activation.ReLU'>\n",
      "registered <class 'torch.nn.modules.linear.Linear'>\n",
      "Inferred reshape on op num 4\n",
      "Test acc:  98.67\n"
     ]
    }
   ],
   "source": [
    "# Let's get the model we trained with PGD and pass it into PytorchDeepZ so we can use \n",
    "# the certification methods on the model\n",
    "\n",
    "zonotope_model = deep_z.PytorchDeepZ(model=trainer._classifier.model, \n",
    "                                     clip_values=(0, 1),\n",
    "                                     optimizer = optim.Adam(model.parameters(), lr=1e-4),\n",
    "                                     loss=nn.CrossEntropyLoss(), \n",
    "                                     input_shape=(1, 28, 28), \n",
    "                                     nb_classes=10)\n",
    "\n",
    "# Get the test time predictions. As always, set the forward_mode to the correct type\n",
    "with torch.no_grad():\n",
    "    zonotope_model.model.set_forward_mode('concrete')\n",
    "    test_preds = zonotope_model.model.forward(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)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "f1edd03d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Classified Correct 1/1 and also certified 1/1\n",
      "Classified Correct 2/2 and also certified 2/2\n",
      "Classified Correct 3/3 and also certified 3/3\n",
      "Classified Correct 4/4 and also certified 3/4\n",
      "Classified Correct 5/5 and also certified 3/5\n",
      "Classified Correct 6/6 and also certified 4/6\n",
      "Classified Correct 7/7 and also certified 4/7\n",
      "Classified Correct 8/8 and also certified 4/8\n",
      "Classified Correct 9/9 and also certified 4/9\n",
      "Classified Correct 10/10 and also certified 4/10\n",
      "Classified Correct 11/11 and also certified 5/11\n",
      "Classified Correct 12/12 and also certified 5/12\n",
      "Classified Correct 13/13 and also certified 5/13\n",
      "Classified Correct 14/14 and also certified 6/14\n",
      "Classified Correct 15/15 and also certified 7/15\n",
      "Classified Correct 16/16 and also certified 7/16\n",
      "Classified Correct 17/17 and also certified 7/17\n",
      "Classified Correct 18/18 and also certified 8/18\n",
      "Classified Correct 19/19 and also certified 8/19\n",
      "Classified Correct 20/20 and also certified 9/20\n",
      "Classified Correct 21/21 and also certified 9/21\n",
      "Classified Correct 22/22 and also certified 9/22\n",
      "Classified Correct 23/23 and also certified 9/23\n",
      "Classified Correct 24/24 and also certified 10/24\n",
      "Classified Correct 25/25 and also certified 10/25\n",
      "Classified Correct 26/26 and also certified 10/26\n",
      "Classified Correct 27/27 and also certified 11/27\n",
      "Classified Correct 28/28 and also certified 12/28\n",
      "Classified Correct 29/29 and also certified 13/29\n",
      "Classified Correct 30/30 and also certified 14/30\n",
      "Classified Correct 31/31 and also certified 15/31\n",
      "Classified Correct 32/32 and also certified 16/32\n",
      "Classified Correct 33/33 and also certified 17/33\n",
      "Classified Correct 34/34 and also certified 17/34\n",
      "Classified Correct 35/35 and also certified 18/35\n",
      "Classified Correct 36/36 and also certified 18/36\n",
      "Classified Correct 37/37 and also certified 18/37\n",
      "Classified Correct 38/38 and also certified 19/38\n",
      "Classified Correct 39/39 and also certified 19/39\n",
      "Classified Correct 40/40 and also certified 20/40\n",
      "Classified Correct 41/41 and also certified 21/41\n",
      "Classified Correct 42/42 and also certified 21/42\n",
      "Classified Correct 43/43 and also certified 22/43\n",
      "Classified Correct 44/44 and also certified 22/44\n",
      "Classified Correct 45/45 and also certified 23/45\n",
      "Classified Correct 46/46 and also certified 23/46\n",
      "Classified Correct 47/47 and also certified 24/47\n",
      "Classified Correct 48/48 and also certified 24/48\n",
      "Classified Correct 49/49 and also certified 25/49\n",
      "Classified Correct 50/50 and also certified 26/50\n"
     ]
    }
   ],
   "source": [
    "# And finally, lets see the PGD model's certified robustness.\n",
    "\n",
    "# We set the appropriate forward mode\n",
    "zonotope_model.model.set_forward_mode('abstract')\n",
    "\n",
    "# And then run the certification. We see that although PGD trainig boosts the certified robustness\n",
    "# in comparison to normal training, it's performance is much weaker than certified training.\n",
    "certification_loop(model=zonotope_model,\n",
    "                   x=np.copy(x_test),\n",
    "                   y=y_test,\n",
    "                   preds=test_preds,\n",
    "                   bound=0.15)"
   ]
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}