decoding.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Copyright (c) Meta Platforms, Inc. and affiliates. All rights reserved."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stable Signature - Decoding\n",
    "\n",
    "Each model has its own key, which can be found in the `keys.txt` file when the fine-tuning is done. \n",
    "The key is a string of 48 bits, which can be converted to a boolean array of 48 elements. \n",
    "The `msg_extractor` is a TorchScript model that extracts the message from the image.\n",
    "\n",
    "Based on the number of matching bits between the key and the message, the image can be classified as genuine or generated by our model."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Imports and setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from PIL import Image\n",
    "import torch\n",
    "import torchvision.transforms as transforms\n",
    "\n",
    "def msg2str(msg):\n",
    "    return \"\".join([('1' if el else '0') for el in msg])\n",
    "\n",
    "def str2msg(str):\n",
    "    return [True if el=='1' else False for el in str]\n",
    "\n",
    "msg_extractor = torch.jit.load(\"path/to/dec_48b_whit.torchscript.pt\").to(\"cuda\")\n",
    "transform_imnet = transforms.Compose([\n",
    "    transforms.ToTensor(),\n",
    "    transforms.Normalize(mean=[0.485, 0.456, 0.406],std=[0.229, 0.224, 0.225])\n",
    "])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Decode images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Extracted message:  111010110101000001010111010011010100010000100111\n"
     ]
    }
   ],
   "source": [
    "img = Image.open(\"path/to/generated/img.png\")\n",
    "img = transform_imnet(img).unsqueeze(0).to(\"cuda\")\n",
    "msg = msg_extractor(img) # b c h w -> b k\n",
    "bool_msg = (msg>0).squeeze().cpu().numpy().tolist()\n",
    "print(\"Extracted message: \", msg2str(bool_msg))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute bit accuracies and run statistical test\n",
    "\n",
    "Metrics are:\n",
    "- **Bit accuracy**: number of matching bits between the key and the message, divided by the total number of bits.\n",
    "- **$p$-value**: probability of observing a bit accuracy as high as the one observed, assuming the null hypothesis that the image is genuine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bit accuracy:  1.0\n",
      "p-value of statistical test:  BinomTestResult(k=48, n=48, alternative='greater', statistic=1.0, pvalue=3.552713678800501e-15)\n"
     ]
    }
   ],
   "source": [
    "key = '111010110101000001010111010011010100010000100111' # model key\n",
    "bool_key = str2msg(key)\n",
    "\n",
    "# compute difference between model key and message extracted from image\n",
    "diff = [bool_msg[i] != bool_key[i] for i in range(len(bool_msg))]\n",
    "bit_acc = 1 - sum(diff)/len(diff)\n",
    "print(\"Bit accuracy: \", bit_acc)\n",
    "\n",
    "# compute p-value\n",
    "from scipy.stats import binomtest\n",
    "pval = binomtest(len(diff)-sum(diff), len(diff), 0.5, alternative='greater')\n",
    "print(\"p-value of statistical test: \", pval)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Typically, one can set a threshold $\\tau$ on the bit accuracy to classify an image as genuine or generated.\n",
    "Here is the table of FPRs for different thresholds on bit accuracy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk0AAAGwCAYAAAC0HlECAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABLjklEQVR4nO3de1yUZcI+8OuZE4zAcBABB1A8HwMUFUkrKcqoNHMre7c2D5u7ubJrkbuLu7+l3HdbayvXDrbs2hrZ265mpW5pLoWn8hBCjmfJA4qKnFQYGIaBmXl+fyCThIcBBu4Z5vp+PnxknnnmnmtuUS+foyTLsgwiIiIiuiGF6ABEREREnoCliYiIiMgJLE1ERERETmBpIiIiInICSxMRERGRE1iaiIiIiJzA0kRERETkBJXoAJ7KbrejpKQEAQEBkCRJdBwiIiJygizLqKmpgV6vh0LRtm1HLE3tVFJSgujoaNExiIiIqB3Onj2LqKioNr2GpamdAgICADRNuk6ng81mAwAolcoOjWs2mwEAWq22YwEBl2XyhrE472LG4ryLGYvzLiYT511Mph/Ou9FoRHR0tOPf8bZgaWqn5l1yOp3OpaVJrVYD4B+qrh6L8y5mLM67mLE472Iycd7FZLrevLfn0BoeCE5ERETkBJYmIiIiIiewNBERERE5gaWJiIiIyAksTUREREROYGkiIiIicgJLExEREZETWJqIiIiInMDSREREROQEry5Nn332GYYMGYJBgwbhnXfeER2HiIiI3JjX3kbFarUiPT0dW7duRWBgIBISEvDQQw+hZ8+eoqMRERGRG/LaLU15eXkYMWIEIiMj4e/vj9TUVOTk5IiORURERG7KY7c07dixA6+88goKCgpw4cIFrFu3DtOmTWuxzvLly/HKK6+gtLQUcXFxePPNNzFu3DgAQElJCSIjIx3rRkZG4vz58135Ea6p2tyIWosVvvVyh8ey2+wAAIWy4924u49VX18PAE7N+81u8mi/cqNJRTtvNHn16M03rVSplJCuPNP89tIPXiBBgiQ1PZSkprUVVxZIUtNcKSRArWp6rFRIUEgSFFL7blxJRORtPLY0mUwmxMXFYc6cOZg+fXqr59esWYP09HRkZWUhMTERy5Ytw+TJk1FYWIiwsLA2v5/FYoHFYnE8NhqNAACz2Qy1Wu2yOzK/ve0kVuWVdGgMovZQXilQCoUEpUKC6qpfm79XKiSolAooJUCjUsBHpYBGpYBG2fL75l+1agV6aJTw81Ghh0bZ4stPo4RStsFPo0CvIDvUHSy+7nindncdq/k/Ca7gjp/PlWO5MhPnXUymH8672Wxu91geW5pSU1ORmpp63eeXLl2KuXPnYvbs2QCArKwsbNy4EStXrkRGRgb0en2LLUvnz593bIW6liVLlmDx4sWu+wDXoVRI8FF57V5TryI7tTFRhtz0C676xfHaK89Clr9/rr1ssgybDMDe8a2c7eHvo0SgVo0grRpBPdQI1KqavteqEahVI9RfgwidD8J1PugV4AONC7ZUEhG1hSTLzv3V7c4kSWqxe66hoQE9evTARx991GKX3cyZM1FVVYUNGzbAarVi2LBh2LZtm+NA8F27dl33QPBrbWmKjo5GdXU1dDqdy1pxcwPWarUdGgdwz/89uOtY3WneZVl2lChZlmGXm8qVzWqDDEBSKGCzX1l+5XmbXXZ8b7XbYbcDjXY7bHYZVpsMm11u8bjBakWD1Q6rHbBY7Wiw2mGx2mCx2q96bEd9ow11DVaYLDaYGqyoa/61wQaTxYpaS9P37RHqr0G4zhcROl+EB/oizF8DfZAWA8L80S/UH8E91O3e7Sj697Czx+pOP++dPZYrM3HexWT64bwbjUYEBgY6/v1uC4/d0nQjlZWVsNlsCA8Pb7E8PDwcx44dAwCoVCq89tprSE5Oht1ux29+85sbnjnn4+MDHx+fTs1N5AqSJEH64QFPAGxXvnWnvxjNZnNTIYMKl+sacLmuEVV1Daiqa8TlH/xaXlOPUmM9yqotaLDZUVnbgMraBhwuMV5zbJ2vCv1C/dAv1A8xV35t/l7nq+5QbiLyTt2yNDlr6tSpmDp1qugYRF5NqZDgr9Ug2E/j1PqyLOOSqaGpQBnrUVptQamxHqVVZpy9XIczF+tQUl0PY70V+89VY/+56lZjRAZpMTJSh5H6QIy48muYztfVH42IupluWZpCQ0OhVCpRVlbWYnlZWRkiIiIEpSIiV5AkCT39fdDT3wcj9IGO5Vdv/TI32HDmkgmnK00oqqxDUWUtTlfW4VSlCZW1FpyvMuN8lRn/Pfz93xGh/j4YGanDCL0OwyMCEBsViOie/l3++YjIfXXL0qTRaJCQkIDc3FzHMU12ux25ublIS0sTG46IOp1Wo8TQCB2GRrQ+XqHa3IgjJUYcLqnG4RIjDp2vxsmKWlTWWrCtsALbCisc60aHaDG+X08k9u+JxH4hiA7p0ZUfg4jcjMeWptraWpw4ccLxuKioCAaDASEhIejTpw/S09Mxc+ZMjBkzBuPGjcOyZctgMpkcZ9MRkXcK1KqRNKAnkgZ8fwyjucGGo6VGHC4x4vD5ahw8X41jpTU4e8mMs5fOYW3BOQBNu/US+4dgfP+eGN+vJ6JDtLzGFZEX8djSlJ+fj+TkZMfj9PR0AE1nyGVnZ2PGjBmoqKhAZmYmSktLER8fj82bN7c6OJyISKtRYnSfYIzuEwygaVdfrcWKfWer8U3RJew5dREHz1XjfJUZn3x7Hp9823S5ksggLZKH9sLdwyMwvn8IfFQdP8ieiNxXt7jkgAg/PGWRpwJ79licdzFjedK8myxWFJy5jG+KLuKbU5ew/1wVGm3f//Xp76PCHYN7IWV4GJKHhCGoh6ZLcrWHJ8276LF4yQExY/GSA0REHszPR4XbB/fC7YN7AWjapbfn1EV8cbQMXx4pQ3mNBRsPXsDGgxegVEgYGxOMu4dH4M4hoejDY6GIugWWJiKidtBqlEgeGobkoWH404MjceB8Nb48UoYvjpShsKwGe05dwp5Tl/C/nwEj9To8nBCFqfGRCHHy0gpE5H5YmoiIOkihkBAfHYT46CAsnDwExRfr8MXRMnxxpBR7T1/GoRIjDpUcwYubjiJ5SBgeTojCpCFh0PCWSUQehaWJiMjF+vTsgZ9O7IefTuyHcqMZn+4vwTpDCQ6dNyLnSBlyjpQhxE+DqXF6PJwQhRF6Hc/CI/IALE1ERJ2op58Gs26NwU9vG4BjpUZ8XHAO6/aVoLLWguxdp5G96zSGhAfgkTFReGRMNAK1vMULkbvitmEioi4yNEKH398/HHsW3Yl3Z43F/bG9oVEpUFhWgz9tPIpbl+Tihf8cxpmLJtFRiegauKWJiKiLqZQKx0Hk1XWN+PRACVbtPo3vymqRves03tt9GinDwvHTif2Q2C+Eu+6I3ARLExGRQIE91HhifF88ntgHX5+oxD+/LsK2wgp8ceVMvBF6HX46sR8eiNVDye5EJBRLExGRG5AkCbcN6oXbBvXCifIarNx5Gp98ew6HS4xI/3A/Xvr8GJ5I7IOfjO+DYH9eeZxIBB7TRETkZgaGBeDPD92C3Rl34deThyAswAflNRYs/fI4bn91O97MPY5ai1V0TCKvw9JEROSmgv00mJ88EF//9k78dUYcBof7o6beite++A63vbwF/9hxEuYGm+iYRF6DpYmIyM1pVAo8NCoKG9MmYNmMOPQP9cPlukb8edMx3P7KVry36zQsVpYnos7G0kRE5CEUCglTYnsj59nb8crDsYgK1qKixoLn/3MYya9sw+q8YjTa7KJjEnVbLE1ERB5GpVTgkTHR2PLcJPxp2kiE63xQUl2PjE8O4q7XtmOD4TxkWRYdk6jbYWkiIvJQGpUCT4zvi+2/TsYfHhiOUH8Nii/VYcFqAx79+24cKTGKjkjUrbA0ERF5OF+1Ej+d2A87fpOMhfcMhlatxN7Tl/HAm18hc8MhVNc1io5I1C2wNBERdRM9NCqk3TkIuc/dgftje8MuA6t2n0Hya03HO9nt3GVH1BEsTURE3Yw+SIvlPx6Nfz2ViEFh/rhkakDGJwfx0Ns7YThbJToekcdiaSIi6qZuHRiKTQtuwx8eGI4AHxX2n6vGtOU78fsNR3HR1CA6HpHHYWkiIurG1EoFfjqxH3IX3oEfjY4CAHxiuIDUN/fgw/yzPMuOqA1YmoiIvEBYgC9eezQOH89LwrAIf9RYrPjNRwfw1Hv5KDfWi45H5BFYmoiIvEhC3xB8OHcM0u/qD41Sgdxj5bj7rzt4bSciJ7A0ERF5GZVCgbkTY/DpLydiZKQO1eZGLFhtwLz/+xaVtRbR8YjcFksTEZGXGhIRgHW/mIBnUwZDpZCw+XAp7vnrDnx+8ILoaERuiaWJiMiLqZUKLEgZhPXzJ2BoRAAumRow74Nv8at/78NlnmFH1AJLExERYWRkIDakTUBa8kAoFRL+s78E9yzbga2F5aKjEbkNliYiIgIA+KiUWDh5CD6ZdysGhvmjosaC2e/uxV82H4PVZhcdj0g4liYiImohLjoIn/1yIp5M6gsAeHvbSTz+zje8NAF5PZYmIiJqxVetxB8fHIk3/2cU/DRKfFN0Cfe98TV2nawUHY1IGJYmIiK6rilxevznlxMxNCIAlbUWPPHON3hry3He/Je8EksTERHd0IBe/lj3iwl4JCEKdhl4Nec7zM7ei0s8u468DEsTERHdlFajxCuPxOEvD8fCR6XA9u8qcP8bX+Hb4suioxF1Ga8tTWfPnsWkSZMwfPhwxMbGYu3ataIjERG5vUfHRGP9/AnoH+qHC9X1+J8VeXh352negoW8gteWJpVKhWXLluHIkSPIycnBM888A5PJJDoWEZHbG9Zbhw1pE3B/bG9Y7TL+tOkYMj4+iAYrL0tA3ZvXlqbevXsjPj4eABAREYHQ0FBcunRJbCgiIg8R4KvGW/8zCv/v/qFQSMCa/LOYuTIP1XWNoqMRdRq3LU07duzAlClToNfrIUkS1q9f32qd5cuXIyYmBr6+vkhMTEReXl673qugoAA2mw3R0dEdTE1E5D0kScLsW2Pwj58kwE+jxO5TF/HQ33bidCW32lP3pBId4HpMJhPi4uIwZ84cTJ8+vdXza9asQXp6OrKyspCYmIhly5Zh8uTJKCwsRFhYGAAgPj4eVqu11WtzcnKg1+sBAJcuXcKTTz6JFStW3DCPxWKBxfL93b+NRiMAwGw2Q61Ww2azAQCUSmX7PvAV9fWuu3icqzJ5w1icdzFjcd7FjOXqeR/fxx8fzEnAvH/tx6kKE6Yt34k3Z9yCMX2D2jwW4F5z5cpM/HkXk+mH8242m9s9liR7wNF7kiRh3bp1mDZtmmNZYmIixo4di7feegsAYLfbER0djV/+8pfIyMhwalyLxYK7774bc+fOxU9+8pMbrvvCCy9g8eLFrZaXlpZCp9O5/C8zX1/fDo0DuOcfBHcdi/MuZizOu5ixOmvey2ssSFt9AAdLaqBWSvjT1GGYGhvRrrFcmcsdxgH48y4q0w/n3Wg0IiIiAtXV1dDpdG0ay223NN1IQ0MDCgoKsGjRIscyhUKBlJQU7N6926kxZFnGrFmzcOedd960MAHAokWLkJ6e7nhsNBoRHR0NrVYLrVbr0t9gANBqtR0ewx3/ILjzWADnXcRYAOddxFiA6+e9r1aLD5+egPQPDfj8UCl+u+4Izlc34Nm7B0OSpDaN5cpc7jDO1fjz3rWZmjXPe2Nj+4+7c9tjmm6ksrISNpsN4eHhLZaHh4ejtLTUqTF27tyJNWvWYP369YiPj0d8fDwOHjx43fV9fHyg0+lafBERUUtajRLLfzwa8yYNAAC8seUEfrXagPpGm+BkRB3nkVuaXGHixImw23l6LBGRqykUEn5771D0C/XD7z45iE/3l+Dc5TqseHIMQv19RMcjajeP3NIUGhoKpVKJsrKyFsvLysoQEeH8/nMiIuo8j46JxqqfjkOgVo19xVV49O+7caG6/QfhEonmkaVJo9EgISEBubm5jmV2ux25ublISkoSmIyIiK5264BQfPKLW6EP9MWpChMe/ttunLnISxKQZ3Lb0lRbWwuDwQCDwQAAKCoqgsFgQHFxMQAgPT0dK1aswHvvvYejR49i3rx5MJlMmD17tsDURET0QwN6+WPtvFvRL9QP56vMeCRrN74rqxEdi6jN3LY05efnY9SoURg1ahSAppI0atQoZGZmAgBmzJiBV199FZmZmYiPj4fBYMDmzZtbHRxORETiRQZp8eHPkzA0IgDlNRY8+vfdOHCuSnQsojbxiOs0uSOj0YjAwEDHdR5cdXpk80W3eEpq147FeRczFuddzFgi572qrgEz392L/Wer4O+jwj9njkFi/57tGsuVuTp7HIA/76Iy/XDef/jvd1u47ZYmIiLqfoJ6aPDBU4kY3z8EtRYrnlyZh22F5aJjETmFpYmIiLqUv48K2bPH4c6hYbBY7Zi7Kh+fH7wgOhbRTbE0ERFRl/NVK5H1RALuj+2NRpuM+f/6Fh9/e150LKIbYmkiIiIhNCoF3nhsFB4dEwW7DPzm44N4f88Z0bGIrouliYiIhFEqJLw0PRazJ8QAAF749Cj+nVcsNhTRdbA0ERGRUAqFhMwHhmPubf0AAL9bdxDr9p0TnIqoNZYmIiISTpIk/HbyYPwksQ9kGXjuw/08OJzcDksTERG5BUmSkPnAMDyc0HSM069W78PWY7wcAbkPliYiInIbCoWEl38UiweunFX38/8rwK4TlaJjEQFgaSIiIjejVEj464x43D08HA1WO55alY/805dExyJiaSIiIvejVirw1o9H4bZBoahrsGH2u3t5rzoSjqWJiIjcko9KiX/8ZAzG9QtBzZVbrhwrNYqORV6MpYmIiNyWVqPEylljERcdhKq6RjzxTh5OVdSKjkVeiqWJiIjcmr+PCqtmj8Pw3jpU1lrw+Dvf4HyVWXQs8kIsTURE5PYCe6jx/k/HYWCYPy5U12P2u3moNjeKjkVehqWJiIg8Qk9/H6yaMw7hOh98V1aLef9XgAarXXQs8iIsTURE5DH0QVqsnDUWfholdp28iIxPDkCWZdGxyEuwNBERkUcZoQ/E8sdHQ6mQ8Mm357Hsy+OiI5GXYGkiIiKPM2lIGP73wZEAgNdzj2Nt/lnBicgbsDQREZFH+nFiH8ybNAAAsOiTg9jJ261QJ2NpIiIij/Xre4ZgSpweVruMp98vQGFpjehI1I2xNBERkcdSKCS8+kgsxsU0XTV89rt5KDPWi45F3RRLExEReTQflRL/eDIB/Xv5oaS6HnOy98JksYqORd0QSxMREXm8oB4aZM8ah55+GhwuMeJXq/fDauM1nMi1WJqIiKhb6NOzB/45ayx81Qps+64Cf9p4THQk6mZYmoiIqNuIjw7C64+NgiQB739TzEsRkEuxNBERUbcyeUQEFtw5EADw+/WHcOBcldhA1G2wNBERUbczf9IApAwNQ4PVjqffL0BlrUV0JOoGWJqIiKjbaboUwS3oH9p0Rl3av77lgeHUYSxNRETULQX4qvGPJxPgp1Fiz6lLeOlzHhhOHcPSRERE3dbAsAC89mgcAOCdr4uwwXBecCLyZCxNRETUrd07sjd+ceUedb/9+ACOlBgFJyJPxdJERETd3nP3DMHtg3uhvtGOn/9fPqrqGkRHIg/k9aWprq4Offv2xcKFC0VHISKiTqJUSHjjsXj0CemBs5fM+NVqA2x2WXQs8jBeX5pefPFFjB8/XnQMIiLqZEE9NMh6IgG+agV2fFeBpV8Uio5EHsarS9Px48dx7NgxpKamio5CRERdYLheh5d/FAsAWL71JDYfKhWciDyJ25amHTt2YMqUKdDr9ZAkCevXr2+1zvLlyxETEwNfX18kJiYiLy+vTe+xcOFCLFmyxEWJiYjIEzwYH4mnJvYDADz3oQGnK02CE5GnUIkOcD0mkwlxcXGYM2cOpk+f3ur5NWvWID09HVlZWUhMTMSyZcswefJkFBYWIiwsDAAQHx8Pq9Xa6rU5OTnYu3cvBg8ejMGDB2PXrl03zWOxWGCxfH9FWaOx6ewLs9kMtVoNm80GAFAqle36vM3q6+s79PqruSqTN4zFeRczFuddzFicd2BBcl/sP3sZe89UIe1fBfhgTgI0ytbbEVyZifMuJtMP591sNrd7LLctTampqTfcbbZ06VLMnTsXs2fPBgBkZWVh48aNWLlyJTIyMgAABoPhuq/fs2cPVq9ejbVr16K2thaNjY3Q6XTIzMy85vpLlizB4sWL2/+BiIjIbagUCvxl+nBMy8rDoZIavL7lFH5990DRscjNSbIsu/3pA5IkYd26dZg2bRoAoKGhAT169MBHH33kWAYAM2fORFVVFTZs2NCm8bOzs3Ho0CG8+uqr113nWluaoqOjUV1dDZ1O57JW3NyAtVpth8YB3PN/D+46FuddzFicdzFjcd6/l3O4FD97vwAAsGrOONw+uFenZeK8i8n0w3k3Go0IDAx0/PvdFm57TNONVFZWwmazITw8vMXy8PBwlJZ2zkF9Pj4+0Ol0Lb6IiMiz3TMiAk8m9QUApH+4HxU1vLEvXZ/b7p7rSrNmzRIdgYiIBPndfcOQV3QJx0prsHDtfrw7aywUCkl0LHJDHrmlKTQ0FEqlEmVlZS2Wl5WVISIiQlAqIiLyRL5qJd78n1HwUSmw/bsKrNxZJDoSuSmPLE0ajQYJCQnIzc11LLPb7cjNzUVSUpLAZERE5IkGhQcgc8pwAMDLm4/h4LlqwYnIHbltaaqtrYXBYHCcAVdUVASDwYDi4mIAQHp6OlasWIH33nsPR48exbx582AymRxn0xEREbXFj8f1wb0jItBok/Gr1ftgsrS+ZA15N7c9pik/Px/JycmOx+np6QCazpDLzs7GjBkzUFFRgczMTJSWliI+Ph6bN29udXA4ERGRMyRJwks/ugUHzlWhqNKE5/9zGC9PHyk6FrkRj7jkgDv64SmLPBXYs8fivIsZi/MuZizO+43lFV3CY//YDbsM/PXRWEyN0/OSA108Fi85QERE5AHG9QvBL+8cBAD4w4bDKL5UJzgRuQuWJiIioh/45Z0DMTYmGLUWG55Zsx+NNrvoSOQGWJqIiIh+QKVUYNljo6DzVWH/uWq8kXtcdCRyAyxNRERE1xAZpMWL05oOBH9720kcOs/LEHg7liYiIqLruO+WCNx3SwRsdhkL1+5Hg5W76bwZSxMREdENvDBlOEL8NDhWWoO3tp4QHYcEYmkiIiK6gZ5+Gvzvg1d20209wd10XoyliYiI6Cbuj+2N+26JgNUu49cfHeBuOi/F0kREROSEPz44EsE91Dh6wYi3t3E3nTdiaSIiInJCqL8P/nhlN91bW07gcAl303kbliYiIiInPRDbG/eOaNpNt3DtAV700suwNBERETlJkiT877SRCGreTbf1pOhI1IVYmoiIiNqgV4APFk8dAQB4c8txHCkxCk5EXYWliYiIqI2mxukxeUT4lbPpeG86b8HSRERE1EZX76Y7XGLE37ZxN503YGkiIiJqh7AA3xa76Y5e4G667o6liYiIqJ2mxulx9/BwNNpk/OajA7DZZdGRqBOxNBEREbWTJEl4cdpIBPiqcPB8NT745ozoSNSJWJqIiIg6IEzni99MHgIAeOW/hSivqReciDoLSxMREVEH/TixL26JDERNvRVLNh0THYc6CUsTERFRBykVEv40bSQkCVi37zx2n7woOhJ1ApYmIiIiF4iLDsLjiX0AAH/YcAgNVl67qbthaSIiInKRX98zFD39NDhRXot/fl0kOg65GEsTERGRiwT2UON39w0DALyRexznLtcJTkSuxNJERETkQtNHR2JcTAjMjTb88dMjouOQC7E0ERERuVDzLVZUCgk5R8qw5ViZ6EjkIixNRERELjYkIgBzJvYDADz/n8Oob7QJTkSuwNJERETUCRbcNQgROl+cvWTG21tPiI5DLsDSRERE1An8fFR4fspwAEDW9lMousiDwj0dSxMREVEnuXdkBO4Y3AsNNjv+d1MhZJk39PVkLE1ERESdRJIkLJ46AhqVArtPXcbmI+WiI1EHsDQRERF1ophQP/xi0gAAwJLNx1FT3yg4EbWXV5emoqIiJCcnY/jw4bjllltgMplERyIiom7o6TsGoE+wFhW1DfjbtpOi41A7eXVpmjVrFv74xz/iyJEj2L59O3x8fERHIiKibshXrcRv7hkIAPjn10U4X2UWnIjaw2tL0+HDh6FWq3HbbbcBAEJCQqBSqQSnIiKi7urOIaEY0zcIFqsdr/63UHQcage3LU07duzAlClToNfrIUkS1q9f32qd5cuXIyYmBr6+vkhMTEReXp7T4x8/fhz+/v6YMmUKRo8ejT//+c8uTE9ERNSSJEn47ZWtTev2nceBc1ViA1Gbue2mFZPJhLi4OMyZMwfTp09v9fyaNWuQnp6OrKwsJCYmYtmyZZg8eTIKCwsRFhYGAIiPj4fVam312pycHFitVnz11VcwGAwICwvDvffei7Fjx+Luu+++Zh6LxQKLxeJ4bDQaAQBmsxlqtRo2W9PVXpVKZYc+d319fYdefzVXZfKGsTjvYsbivIsZi/MuJlN9fT0GhmgwJTYcnx4owx8/PYxVM0dBkiShudxxLFfP+9XM5vbvGnXb0pSamorU1NTrPr906VLMnTsXs2fPBgBkZWVh48aNWLlyJTIyMgAABoPhuq+PjIzEmDFjEB0dDQC47777YDAYrlualixZgsWLF7fz0xARETV59s4ByDlSgfwzVcgtrETK0F6iI5GT3LY03UhDQwMKCgqwaNEixzKFQoGUlBTs3r3bqTHGjh2L8vJyXL58GYGBgdixYwd+/vOfX3f9RYsWIT093fHYaDQiOjoaWq0WWq3Wpa0YALRabYfHcMf/PbjzWADnXcRYAOddxFgA572rMzXrH6zFU7f1w/KtJ7E09xQm3xIFjaptR8u441y5cqzOmPfmn/fGxvZf8sFtj2m6kcrKSthsNoSHh7dYHh4ejtLSUqfGUKlU+POf/4zbb78dsbGxGDRoEB544IHrru/j4wOdTtfii4iIqD3mTRqIUH8NiipN+OCbM6LjkJM8sjS5SmpqKg4ePIhDhw5h6dKlouMQEZGX8PdR4dm7BwMAXs89jmozL3jpCTyyNIWGhkKpVKKsrKzF8rKyMkRERAhKRURE5LwZY6IxKMwfVXWNWL71hOg45ASPLE0ajQYJCQnIzc11LLPb7cjNzUVSUpLAZERERM5RKRX43X3DAADZO0/j7KU6wYnoZty2NNXW1sJgMDjOgCsqKoLBYEBxcTEAID09HStWrMB7772Ho0ePYt68eTCZTI6z6YiIiNzdpCG9MHFgKBpsdry0+ZjoOHQTblua8vPzMWrUKIwaNQpAU0kaNWoUMjMzAQAzZszAq6++iszMTMTHx8NgMGDz5s2tDg4nIiJyV5Ik4Xf3DYMkARsPXEDBmcuiI9ENuO0lByZNmgRZlm+4TlpaGtLS0rooERERkesN1+vwSEIUPsw/hz9tPIJP5t3argteUudz2y1NRERE3uK5e4ZAq1ZiX3EVNh68IDoOXQdLExERkWDhOl/87Pb+AICXNx+DxWoTnIiuhaWJiIjIDfz8jv4IC/DB2UtmrNrFC166I5YmIiIiN9BDo8LCe4YAAN7edgK1ltY3nCexWJqIiIjcxPTRkegf6ofLdY3I3lkkOg79AEsTERGRm1ApFViQMggA8I8dp3h7FTfD0kRERORGHojVY1CYP4z1Vqz8mlub3AlLExERkRtRKiQ8k9J0M9+VXxehqq5BcCJqxtJERETkZlJHRmBoRABqLFb8Y8cp0XHoCpYmIiIiN6NQSEi/u2lrU/au07hYaxGciACWJiIiIrd09/Bw3BIZiLoGG/7OrU1uweWl6aOPPnL1kERERF5Hkr7f2rRq92mU19QLTkRtLk1WqxWHDh3Cd99912L5hg0bEBcXh8cff9xl4YiIiLzZpCG9MKpPEOob7fjbtpOi43i9NpWmQ4cOYeDAgYiLi8OwYcMwffp0lJWV4Y477sCcOXOQmpqKkyf5m0pEROQKkiThububrhL+wTfFuFBtFpzIu7WpNP32t7/FwIEDsWHDBjz22GNYv349Jk2ahClTpuDcuXN46aWXEBUV1VlZiYiIvM6EgT0xLiYEDVY7lm89ITqOV2tTadq7dy9effVVPPDAA3j77bcBAL/73e+wcOFCaLXaTglIRETkzSRJQvo9Tcc2rdl7Fucu1wlO5L3aVJoqKyuh1+sBAIGBgfDz88P48eM7JRgRERE1Gd+/JyYM7IlGm4y3tnBrkyhtKk2SJKGmpgZGoxHV1dWQJAlmsxlGo7HFFxEREblW85l0awvO4cxFbm0SoU2lSZZlDB48GMHBwQgJCUFtbS1GjRqF4OBgBAcHIygoCMHBwZ2VlYiIyGsl9A3BHYN7wWaX8RaPbRJC1ZaVt27d2lk5iIiI6CbS7x6M7d9VYL2hBPPu6I9BEYGiI3mVNpWmO+64o7NyEBER0U3ERQchZVgYvjxajje2nMSbPx4tOpJXadPuObvdjpdffhkTJkzA2LFjkZGRAbOZ14wgIiLqKs9eObbps4MXcLKiVnAa79Km0vTiiy/id7/7Hfz9/REZGYnXX38d8+fP76xsRERE9AMj9IFIGRoGWQb+sZ33pOtKbSpNq1atwttvv43//ve/WL9+PT799FN88MEHsNvtnZWPiIiIfuBnt/cDAKzbdx5lRt6Trqu0qTQVFxfjvvvuczxOSUmBJEkoKSlxeTAiIiK6toS+wRjTNxgNNjtWfl0kOo7XaFNpslqt8PX1bbFMrVajsbHRpaGIiIjoxn5+ZWvTB98Uo9rMf4e7QpvOnpNlGbNmzYKPj49jWX19PZ5++mn4+fk5ln3yySeuS0hEREStTBrcC4PD/fFdWS0++OYMfjFpoOhI3V6btjQ9+eSTCAsLQ2BgoOPriSeegF6vb7GMiIiIOpdCIeHpOwYAAFZ+fRr1jTbBibq/Nm1pyszMRExMDBSKNnUtIiIi6gRT4vR4Lec7nK8y4+Nvz+HxxL6iI3VrbWo/gwYNQmVlpePxjBkzUFZW5vJQREREdHNqpQI/ndh0bNOKHadgs8uCE3Vvbb733NU2bdoEk8nk0kBERETkvMfGRSOohxqnL9Zh86FS0XG6Ne5nIyIi8mA9NCo8mRQDAMjafrLVBg5ynTaVJkmSIElSq2VEREQkzqxbY+CrVuDg+WrsOnlRdJxuq0OXHLjW5QYAz7nkwF//+le88847kGUZKSkpeP3111kCiYjI44T4afDY2D7I3nUaf9t2EhMGhoqO1C21qTTNnDmzxeMnnnjCpWG6UkVFBd566y0cPnwYarUat99+O/bs2YOkpCTR0YiIiNrspxP74f09Z/D1iUocPFeNW6J4CSBXa1NpevfddzsrhxBWqxX19U337GlsbERYWJjgRERERO0THdIDU2J7Y72hBFk7TmL5j0eLjtTtuO2B4Dt27MCUKVOg1+shSRLWr1/fap3ly5cjJiYGvr6+SExMRF5entPj9+rVCwsXLkSfPn2g1+uRkpKCAQMGuPATEBERda2fX7nY5ecHL+DMRZ7d7mpt2tLUlUwmE+Li4jBnzhxMnz691fNr1qxBeno6srKykJiYiGXLlmHy5MkoLCx0bDGKj4+H1Wpt9dqcnBxotVp89tlnOH36NLRaLVJTU7Fjxw7cfvvt18xjsVhgsVgcj41GIwDAbDZDrVbDZmu6EqtSqezQ527e8uUKrsrkDWNx3sWMxXkXMxbnXUymrpj3mCA1bhvYE1+duIi3txzHCw8MafdYrswlahyg9bybzeZ2j+W2pSk1NRWpqanXfX7p0qWYO3cuZs+eDQDIysrCxo0bsXLlSmRkZAAADAbDdV+/du1aDBw4ECEhIQCA+++/H3v27LluaVqyZAkWL17czk9DRETUNZ6a0AdfnbiIdYYLmD8pBr38fW7+InKK25amG2loaEBBQQEWLVrkWKZQKJCSkoLdu3c7NUZ0dDR27dqF+vp6qNVqbNu2DT/72c+uu/6iRYuQnp7ueGw0GhEdHQ2tVgutVuvSVgwAWq22w2O44/8e3HksgPMuYiyA8y5iLIDz3tWZmnX2vN8+1Bej+hRhX3EVVheU4jf3Dm33WK7MJWKcqzXPe2NjY7vHcNtjmm6ksrISNpsN4eHhLZaHh4ejtNS5q6GOHz8e9913H0aNGoXY2FgMGDAAU6dOve76Pj4+0Ol0Lb6IiIjcjSR9fyPf9/ecQU19+0sCteSRpclVXnzxRRw9ehSHDx/GG2+8wWs0ERFRt3D3sHAM6OWHmnor1uw9KzpOt+GRpSk0NBRKpbLVzYLLysoQEREhKBUREZF7UCgkPHVbfwDAqt1neCNfF/HI0qTRaJCQkIDc3FzHMrvdjtzcXF6ckoiICMCD8XrofFUovlSH7d+Vi47TLbhtaaqtrYXBYHCcAVdUVASDwYDi4mIAQHp6OlasWIH33nsPR48exbx582AymRxn0xEREXmzHhoVHh0TDQB4b9cZwWm6B7c9ey4/Px/JycmOx81nrs2cORPZ2dmYMWMGKioqkJmZidLSUsTHx2Pz5s2tDg4nIiLyVj9J6ot/7izC9u8qUFRpQr9Qv5u/iK7LbUvTpEmTIMs33geblpaGtLS0LkpERETkWfr29EPykDBsOVaO93efQeaU4aIjeTS33T1HREREHfdkUl8AwNqCszBZWt8lg5zH0kRERNSN3T6oF2J69kBNvRXrDedFx/FoLE1ERETdmEIh4SdJMQCA93advumhL3R9LE1ERETd3MMJUeihUeK7slrsOXVJdByPxdJERETUzQVq1XhoVCQAYNXu02LDeDCWJiIiIi/w5JVddDlHylBSZRYbxkOxNBEREXmBIREBGN8/BDa7jH99Uyw6jkdiaSIiIvISM69sbfp3XjEsVpvYMB6IpYmIiMhL3D08HL0DfXHR1IBNBy+IjuNxWJqIiIi8hEqpwBPjmy52mc370bUZSxMREZEXmTE2GhqlAvvPVsFwtkp0HI/C0kRERORFQv198EBsbwC8/EBbsTQRERF5mSdvjQEAfLb/Ai6aGsSG8SAsTURERF4mPjoIcVGBaLDZ8WH+WdFxPAZLExERkRdqvtjlB9+chdVmFxvGQ7A0EREReaH7Y3ujp58GF6rrkXusQnQcj8DSRERE5IV81Uo8Ni4aAPD+Hl5+wBksTURERF7qx4l9IUnA7lOXUHyxTnQct8fSRERE5KUig7SYOKAnAOCjAh4QfjMsTURERF7s4YQoAMBHBedgs8uC07g3liYiIiIvdvewMOh8VSiprseuk5Wi47g1liYiIiIv5qNWYmqcHgCwNv+c4DTujaWJiIjIyz2cEAkA2Hy4FNV1jYLTuC+WJiIiIi83Uq/D0IgANFjt+M+BEtFx3BZLExERkZeTJAmPjGm6ZtNa3lbluliaiIiICNPi9VApJBw4V43C0hrRcdwSSxMRERGhp78P7hoWBoBbm66HpYmIiIgAAI9e2UW3bt95NPImvq2wNBEREREA4I7BvdArwAcXTQ3YcqxcdBy3w9JEREREAACVUoHpo5suP8BddK2xNBEREZHDIwlNu+i2FlagvKZecBr3wtJEREREDgPD/DGqTxBsdhnr950XHceteEVpeuihhxAcHIyHH3641XOfffYZhgwZgkGDBuGdd94RkI6IiMi9NB8Q/mH+Ocgyb+LbzCtK04IFC7Bq1apWy61WK9LT07Flyxbs27cPr7zyCi5evCggIRERkft4ILY3fNUKnCivheFsleg4bsMrStOkSZMQEBDQanleXh5GjBiByMhI+Pv7IzU1FTk5OQISEhERuY8AXzXuG9kbQNPWJmoivDTt2LEDU6ZMgV6vhyRJWL9+fat1li9fjpiYGPj6+iIxMRF5eXkuee+SkhJERkY6HkdGRuL8ee6/JSIienhMFADgs/0lMDfYBKdxDyrRAUwmE+Li4jBnzhxMnz691fNr1qxBeno6srKykJiYiGXLlmHy5MkoLCxEWFjTlUvj4+NhtVpbvTYnJwd6vd4lOS0WCywWi+Ox0WgEAJjNZqjVathsTT9QSqWyQ+9TX++6MxVclckbxuK8ixmL8y5mLM67mEyeNu9xvXsgMsgX56vq8em+YkyJjeiyXJ0572azud1jCS9NqampSE1Nve7zS5cuxdy5czF79mwAQFZWFjZu3IiVK1ciIyMDAGAwGNr13nq9vsWWpfPnz2PcuHHXXHfJkiVYvHhxu96HiIjI0ygkCQ/F98Zb24rwieGCU6WpuxNemm6koaEBBQUFWLRokWOZQqFASkoKdu/e3eHxx40bh0OHDuH8+fMIDAzE559/jj/84Q/XXHfRokVIT093PDYajYiOjoZWq4VWq3VpKwYArVbb4THc8X9t7jwWwHkXMRbAeRcxFsB57+pMzTxp3h9LjMHy7UXYU3QZlWYZ0SE9uiRXZ857Y2Nju8cQfkzTjVRWVsJmsyE8PLzF8vDwcJSWljo9TkpKCh555BFs2rQJUVFRjsKlUqnw2muvITk5GfHx8XjuuefQs2fPa47h4+MDnU7X4ouIiKg7iwrugQkDQgEAHxXwgHC33tLkKl9++eV1n5s6dSqmTp3ahWmIiIg8xyNjovD1iUp8VHAOC+4aBIVCEh1JGLfe0hQaGgqlUomysrIWy8vKyhARwX2rREREnW3yiAgE+KpwvsqMPUXefS1Dty5NGo0GCQkJyM3NdSyz2+3Izc1FUlKSwGRERETewVetdFyzaeOBC4LTiCW8NNXW1sJgMDjOgCsqKoLBYEBxcTEAID09HStWrMB7772Ho0ePYt68eTCZTI6z6YiIiKhzPRDXVJo2HyqF1WYXnEYc4cc05efnIzk52fG4+Qy1mTNnIjs7GzNmzEBFRQUyMzNRWlqK+Ph4bN68udXB4URERNQ5kvr3RIifBhdNDdhz6hImDgoVHUkI4aVp0qRJN70ZYFpaGtLS0rooEREREV1NpVTg3pER+Nc3xfjsQInXlibhu+eIiIjI/T1wy5VddIdL0eilu+hYmoiIiOimxvULQai/BlV1jdh5olJ0HCFYmoiIiOimVEoFUr38LDqWJiIiInLK/bFNpem/h0vRYPW+XXQsTUREROSUsTEhCAvwgbHeiq9PVIiO0+VYmoiIiMgpSoWE+64cEP7Zfu/bRcfSRERERE574Mouui+OlKG+0SY4TddiaSIiIiKnje4TjAidL2osVuz4zrt20bE0ERERkdMUCslxQPjGg961i46liYiIiNqkuTR96WW76FiaiIiIqE1GRQchMkgLU4MN2wrLRcfpMixNRERE1CaS9P0uuk+96EKXLE1ERETUZs1n0W05Wo66BqvgNF2DpYmIiIja7JbIQPQJ6QFzow1bjnnHLjqWJiIiImqzq3fRecu96FiaiIiIqF3uv3J18C3HylFr6f676FiaiIiIqF1G6HXoF+oHi9WO3KNlouN0OpYmIiIiahdJkhxbmz7zgl10LE1ERETUbg/ENZWm7YUVqKnv3rvoWJqIiIio3YaEB2BALz802OzI7eZn0bE0ERERUbs1nUWnB9D9z6JjaSIiIqIOab7Q5VcnKmE0NwpO03lYmoiIiKhDBocHYHC4PxptMnK68Vl0LE1ERETUYQ9c2UX330MsTURERETXlTIsHACw69RF1DfaBKfpHCxNRERE1GHDegcgItAX9Y127D51UXScTsHSRERERB0mSRImDQ4FAGztppceYGkiIiIil0geEgag6V50siwLTuN6LE1ERETkErcOCIFGpcC5y2acKK8VHcflWJqIiIjIJXpoVEjsFwKgaWtTd8PSRERERC6TPKQXAJYmIiIiohtqLk35Zy7DWN+9rg7O0kREREQu0yekBwb08oPNLuOr7ypFx3EpryhNDz30EIKDg/Hwww+3WH727FlMmjQJw4cPR2xsLNauXSsoIRERUfdx59Dvz6LrTryiNC1YsACrVq1qtVylUmHZsmU4cuQIcnJy8Mwzz8BkMglISERE1H0kXylN278rh93efS494BWladKkSQgICGi1vHfv3oiPjwcAREREIDQ0FJcuXeridERERN3L2JgQ+PuoUFnbgAPnq0XHcRnhpWnHjh2YMmUK9Ho9JEnC+vXrW62zfPlyxMTEwNfXF4mJicjLy3N5joKCAthsNkRHR7t8bCIiIm+iVipw26Cmq4N3p110KtEBTCYT4uLiMGfOHEyfPr3V82vWrEF6ejqysrKQmJiIZcuWYfLkySgsLERYWNPmv/j4eFit1lavzcnJgV6vv2mGS5cu4cknn8SKFSuuu47FYoHFYnE8NhqNAACz2Qy1Wg2brenmhEql8qbvdyP19fUdev3VXJXJG8bivIsZi/MuZizOu5hM3jbvE/sH4fNDpcg9Uop5E9u2QaIz591sNrd7LOGlKTU1Fampqdd9funSpZg7dy5mz54NAMjKysLGjRuxcuVKZGRkAAAMBkO7399isWDatGnIyMjArbfeet31lixZgsWLF7f7fYiIiLzJbYN6AgAOX6hBeY0FYQE+ghN1nPDSdCMNDQ0oKCjAokWLHMsUCgVSUlKwe/fuDo8vyzJmzZqFO++8Ez/5yU9uuO6iRYuQnp7ueGw0GhEdHQ2tVgutVuvSVgwAWq22w2N4wv9E3GksgPMuYiyA8y5iLIDz3tWZmnnLvPfRahEbFYgD56rxzZkaPDo2SEimZs3z3tjY/mtHCT+m6UYqKyths9kQHh7eYnl4eDhKS0udHiclJQWPPPIINm3ahKioKEfh2rlzJ9asWYP169cjPj4e8fHxOHjw4DXH8PHxgU6na/FFRERE13f1DXy7A7fe0uQqX3755TWXT5w4EXa7vYvTEBEReYc7h4bh9dzj+PpEJRqsdmhUbr2t5qbcOn1oaCiUSiXKyspaLC8rK0NERISgVEREROSMWyIDEervg1qLFfmnPf+SPm5dmjQaDRISEpCbm+tYZrfbkZubi6SkJIHJiIiI6GYUCgmTutENfIWXptraWhgMBscZcEVFRTAYDCguLgYApKenY8WKFXjvvfdw9OhRzJs3DyaTyXE2HREREbkvxy1VCj2/NAk/pik/Px/JycmOx81nqM2cORPZ2dmYMWMGKioqkJmZidLSUsTHx2Pz5s2tDg4nIiIi9zNxUChUCgmnKkw4c9GEvj39REdqN+GladKkSZDlG9+XJi0tDWlpaV2UiIiIiFxF56vG2JgQ7D51EVuOlWP2hH6iI7Wb8N1zRERE1L05dtF5+HFNLE1ERETUqZKvlKZvTl2CydL6tmeegqWJiIiIOtWAXn7oE9IDDTY7dp6oFB2n3ViaiIiIqFNJkoTkK5ce2OrBZ9GxNBEREVGna95Ft/VYxU1PAHNXLE1ERETU6cb37wmtWolSYz2OXqgRHaddWJqIiIio0/mqlZgwsCcAz91Fx9JEREREXSLZwy89wNJEREREXaL5ek37ii+juq5RcJq2Y2kiIiKiLtE7UIv+oX6wy0BB8SXRcdqMpYmIiIi6zNiYEABAXtFlwUnajqWJiIiIuszYfs2l6aLgJG3H0kRERERdZtyVLU0Hz1ejvtEmOE3bsDQRERFRl4kO0SJc54NGm4x9xVWi47QJSxMRERF1GUmSMK5f0/Wa9p72rIPBWZqIiIioS42LCQbA0kRERER0Q80HgxecuQyrzS44jfNYmoiIiKhLDQ4LQKBWjboGGw6XGEXHcRpLExEREXUphULCmL6et4uOpYmIiIi63PfXa2JpIiIiIrqucVdK097TlyDLsuA0zmFpIiIioi43Uh8IX7UCl+sacaK8VnQcp7A0ERERUZfTqBQYFd10XFOehxzXxNJEREREQjQf17TXQ45rYmkiIiIiIZrvQ7f39GXBSZzD0kRERERCjO4bBJVCwvkqM85drhMd56ZYmoiIiEiIHhoVRkQGAvCM6zWxNBEREZEwzfehyyty/110LE1EREQkzNiY76/X5O5YmoiIiEiY5tJ0orwWl0wNgtPcGEsTERERCRPsp8GgMH8A7r+1iaWJiIiIhBrnIfeh84rS9NBDDyE4OBgPP/zwNZ+vq6tD3759sXDhwi5ORkRERFffh86deUVpWrBgAVatWnXd51988UWMHz++CxMRERFRs+bjmg6XGGGyWAWnuT6vKE2TJk1CQEDANZ87fvw4jh07htTU1C5ORURERACgD9IiMkgLm13Gt8Xue+kB4aVpx44dmDJlCvR6PSRJwvr161uts3z5csTExMDX1xeJiYnIy8tz2fsvXLgQS5Yscdl4RERE1HaecFyTSnQAk8mEuLg4zJkzB9OnT2/1/Jo1a5Ceno6srCwkJiZi2bJlmDx5MgoLCxEWFgYAiI+Ph9XaenNeTk4O9Hr9dd97w4YNGDx4MAYPHoxdu3bdMKfFYoHFYnE8NhqNAACz2Qy1Wg2bzQYAUCqVN//QN1BfX9+h11/NVZm8YSzOu5ixOO9ixuK8i8nEeb+x+Eh/rNsH7DlZCfOtkS7L9MN5N5vN7R5LeGlKTU294a6xpUuXYu7cuZg9ezYAICsrCxs3bsTKlSuRkZEBADAYDO167z179mD16tVYu3Ytamtr0djYCJ1Oh8zMzFbrLlmyBIsXL27X+xAREdGNjekbBAA4cN6IBqsdGpXwnWGtCC9NN9LQ0ICCggIsWrTIsUyhUCAlJQW7d+/u8PhLlixx7JrLzs7GoUOHrlmYAGDRokVIT093PDYajYiOjoZWq4VWq3VpUwcArVbb4THc5X8PnjIWwHkXMRbAeRcxFsB57+pMzTjv1zY8yhc9/TS4aGrAiUsNSOgb3Cnz3tjY2O4x3K/GXaWyshI2mw3h4eEtloeHh6O0tNTpcVJSUvDII49g06ZNiIqKalfh8vHxgU6na/FFREREriFJEsZcuQ9d/hn3PBjcrbc0ucqXX35503VmzZrV+UGIiIjousbGhOC/h8uQV3QZP79ddJrW3HpLU2hoKJRKJcrKylosLysrQ0REhKBURERE1Bmaz6ArKL4Mm10WnKY1ty5NGo0GCQkJyM3NdSyz2+3Izc1FUlKSwGRERETkasN76+CnUaKm3orvympEx2lF+O652tpanDhxwvG4qKgIBoMBISEh6NOnD9LT0zFz5kyMGTMG48aNw7Jly2AymRxn0xEREVH3oFIqMLpvML46Xom9py9jZFSw6EgtCC9N+fn5SE5OdjxuPkNt5syZyM7OxowZM1BRUYHMzEyUlpYiPj4emzdvbnVwOBEREXm+cTEhTaXpzGXMnig6TUvCS9OkSZMgyzfeb5mWloa0tLQuSkRERESijG2+eW/RZciyDEmSBCf6nlsf00RERETeJT46CGqlhIpaC85crBMdpwWWJiIiInIbvmolYqMCAQB5p93rPnQsTURERORWxvYNRg+NElV1DaKjtCD8mCYiIiKiqz19xwA8kzIIvhq16CgtsDQRERGRWwnwdc96wt1zRERERE5gaSIiIiJyAksTERERkRNYmoiIiIicwNJERERE5ASWJiIiIiInsDQREREROYGliYiIiMgJLE1ERERETmBpIiIiInICSxMRERGRE1iaiIiIiJzA0kRERETkBPe8jbAHkGUZAGA0GgEANpsNAKBUKjs0rtlsBgA0NjZ2aBxXZvKGsTjvYsbivIsZi/MuJhPnXUymH85787/bzf+OtwVLUzvV1NQAAKKjowUnISIioraqqalBYGBgm14jye2pWgS73Y6SkhIEBARAkiQAwNixY7F3794OjWs0GhEdHY2zZ89Cp9N1OKcrMnnDWJx3MWNx3sWMxXkXMw7nXcw4P5x3WZZRU1MDvV4PhaJtRylxS1M7KRQKREVFtVimVCpd8gcBAHQ6nUvGcmUmbxiL8y5mLM67mLE47107TjPOe9eO0+zqeW/rFqZmPBDchebPny86QiuuzOQNY7mKu34+dx3LVdz187nrWK7irp/PVWO545wD7jlXrhzLHeedu+fcjNFoRGBgIKqrq13asOnGOO9icN7F4LyLwXkXw5Xzzi1NbsbHxwfPP/88fHx8REfxKpx3MTjvYnDexeC8i+HKeeeWJiIiIiIncEsTERERkRNYmoiIiIicwNJERERE5ASWJiIiIiInsDS5meXLlyMmJga+vr5ITExEXl6e6Ejdyo4dOzBlyhTo9XpIkoT169e3eF6WZWRmZqJ3797QarVISUnB8ePHxYTtRpYsWYKxY8ciICAAYWFhmDZtGgoLC1usU19fj/nz56Nnz57w9/fHj370I5SVlQlK3D387W9/Q2xsrOOifklJSfj8888dz3POO99LL70ESZLwzDPPOJZx3jvHCy+8AEmSWnwNHTrU8bwr5p2lyY2sWbMG6enpeP755/Htt98iLi4OkydPRnl5ueho3YbJZEJcXByWL19+zef/8pe/4I033kBWVha++eYb+Pn5YfLkyaivr+/ipN3L9u3bMX/+fOzZswdffPEFGhsbcc8998BkMjnWefbZZ/Hpp59i7dq12L59O0pKSjB9+nSBqT1fVFQUXnrpJRQUFCA/Px933nknHnzwQRw+fBgA57yz7d27F3//+98RGxvbYjnnvfOMGDECFy5ccHx9/fXXjudcMu8yuY1x48bJ8+fPdzy22WyyXq+XlyxZIjBV9wVAXrduneOx3W6XIyIi5FdeecWxrKqqSvbx8ZH//e9/C0jYfZWXl8sA5O3bt8uy3DTParVaXrt2rWOdo0ePygDk3bt3i4rZLQUHB8vvvPMO57yT1dTUyIMGDZK/+OIL+Y477pAXLFggyzJ/1jvT888/L8fFxV3zOVfNO7c0uYmGhgYUFBQgJSXFsUyhUCAlJQW7d+8WmMx7FBUVobS0tMXvQWBgIBITE/l74GLV1dUAgJCQEABAQUEBGhsbW8z90KFD0adPH869i9hsNqxevRomkwlJSUmc8042f/583H///S3mF+DPemc7fvw49Ho9+vfvj8cffxzFxcUAXDfvvGGvm6isrITNZkN4eHiL5eHh4Th27JigVN6ltLQUAK75e9D8HHWc3W7HM888gwkTJmDkyJEAmuZeo9EgKCioxbqc+447ePAgkpKSUF9fD39/f6xbtw7Dhw+HwWDgnHeS1atX49tvv8XevXtbPcef9c6TmJiI7OxsDBkyBBcuXMDixYtx22234dChQy6bd5YmIupS8+fPx6FDh1oca0CdZ8iQITAYDKiursZHH32EmTNnYvv27aJjdVtnz57FggUL8MUXX8DX11d0HK+Smprq+D42NhaJiYno27cvPvzwQ2i1Wpe8B3fPuYnQ0FAolcpWR/KXlZUhIiJCUCrv0jzP/D3oPGlpafjss8+wdetWREVFOZZHRESgoaEBVVVVLdbn3HecRqPBwIEDkZCQgCVLliAuLg6vv/4657yTFBQUoLy8HKNHj4ZKpYJKpcL27dvxxhtvQKVSITw8nPPeRYKCgjB48GCcOHHCZT/vLE1uQqPRICEhAbm5uY5ldrsdubm5SEpKEpjMe/Tr1w8REREtfg+MRiO++eYb/h50kCzLSEtLw7p167Blyxb069evxfMJCQlQq9Ut5r6wsBDFxcWcexez2+2wWCyc805y11134eDBgzAYDI6vMWPG4PHHH3d8z3nvGrW1tTh58iR69+7tup/3Dh6sTi60evVq2cfHR87OzpaPHDki/+xnP5ODgoLk0tJS0dG6jZqaGnnfvn3yvn37ZADy0qVL5X379slnzpyRZVmWX3rpJTkoKEjesGGDfODAAfnBBx+U+/XrJ5vNZsHJPdu8efPkwMBAedu2bfKFCxccX3V1dY51nn76ablPnz7yli1b5Pz8fDkpKUlOSkoSmNrzZWRkyNu3b5eLiorkAwcOyBkZGbIkSXJOTo4sy5zzrnL12XOyzHnvLM8995y8bds2uaioSN65c6eckpIih4aGyuXl5bIsu2beWZrczJtvvin36dNH1mg08rhx4+Q9e/aIjtStbN26VQbQ6mvmzJmyLDddduAPf/iDHB4eLvv4+Mh33XWXXFhYKDZ0N3CtOQcgv/vuu451zGaz/Itf/EIODg6We/ToIT/00EPyhQsXxIXuBubMmSP37dtX1mg0cq9eveS77rrLUZhkmXPeVX5YmjjvnWPGjBly7969ZY1GI0dGRsozZsyQT5w44XjeFfMuybIsu2hLGBEREVG3xWOaiIiIiJzA0kRERETkBJYmIiIiIiewNBERERE5gaWJiIiIyAksTUREREROYGkiIiIicgJLExEREZETWJqIyONs27YNkiS1uvlmZ8vOzkZQUFCHxjh9+jQkSYLBYLjuOqI+HxHdGEsTEbkVSZJu+PXCCy+IjkhEXkolOgAR0dUuXLjg+H7NmjXIzMxEYWGhY5m/vz/y8/PbPG5DQwM0Go1LMhKRd+KWJiJyKxEREY6vwMBASJLUYpm/v79j3YKCAowZMwY9evTArbfe2qJcvfDCC4iPj8c777yDfv36wdfXFwBQVVWFp556Cr169YJOp8Odd96J/fv3O163f/9+JCcnIyAgADqdDgkJCa1K2n//+18MGzYM/v7+uPfee1sUPbvdjj/+8Y+IioqCj48P4uPjsXnz5ht+5k2bNmHw4MHQarVITk7G6dOnOzKFRNRJWJqIyGP9/ve/x2uvvYb8/HyoVCrMmTOnxfMnTpzAxx9/jE8++cRxDNEjjzyC8vJyfP755ygoKMDo0aNx11134dKlSwCAxx9/HFFRUdi7dy8KCgqQkZEBtVrtGLOurg6vvvoq3n//fezYsQPFxcVYuHCh4/nXX38dr732Gl599VUcOHAAkydPxtSpU3H8+PFrfoazZ89i+vTpmDJlCgwGA5566ilkZGS4eKaIyCVkIiI39e6778qBgYGtlm/dulUGIH/55ZeOZRs3bpQByGazWZZlWX7++edltVotl5eXO9b56quvZJ1OJ9fX17cYb8CAAfLf//53WZZlOSAgQM7Ozr5uHgDyiRMnHMuWL18uh4eHOx7r9Xr5xRdfbPG6sWPHyr/4xS9kWZbloqIiGYC8b98+WZZledGiRfLw4cNbrP/b3/5WBiBfvnz5mjmISAxuaSIijxUbG+v4vnfv3gCA8vJyx7K+ffuiV69ejsf79+9HbW0tevbsCX9/f8dXUVERTp48CQBIT0/HU089hZSUFLz00kuO5c169OiBAQMGtHjf5vc0Go0oKSnBhAkTWrxmwoQJOHr06DU/w9GjR5GYmNhiWVJSktNzQERdhweCE5HHunq3mSRJAJqOKWrm5+fXYv3a2lr07t0b27ZtazVW86UEXnjhBfz4xz/Gxo0b8fnnn+P555/H6tWr8dBDD7V6z+b3lWXZFR+HiNwctzQRkdcYPXo0SktLoVKpMHDgwBZfoaGhjvUGDx6MZ599Fjk5OZg+fTreffddp8bX6XTQ6/XYuXNni+U7d+7E8OHDr/maYcOGIS8vr8WyPXv2tPGTEVFXYGkiIq+RkpKCpKQkTJs2DTk5OTh9+jR27dqF3//+98jPz4fZbEZaWhq2bduGM2fOYOfOndi7dy+GDRvm9Hv8+te/xssvv4w1a9agsLAQGRkZMBgMWLBgwTXXf/rpp3H8+HH8+te/RmFhIf71r38hOzvbRZ+YiFyJu+eIyGtIkoRNmzbh97//PWbPno2KigpERETg9ttvR3h4OJRKJS5evIgnn3wSZWVlCA0NxfTp07F48WKn3+NXv/oVqqur8dxzz6G8vBzDhw/Hf/7zHwwaNOia6/fp0wcff/wxnn32Wbz55psYN24c/vznP7c6E5CIxJNk7ownIiIiuinuniMiIiJyAksTERERkRNYmoiIiIicwNJERERE5ASWJiIiIiInsDQREREROYGliYiIiMgJLE1ERERETmBpIiIiInICSxMRERGRE1iaiIiIiJzw/wHVc4ZtZFmq2AAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.stats import binom\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def fpr(threshold, n):\n",
    "    return binom.cdf(n-threshold, n, 0.5)\n",
    "\n",
    "thresholds = np.linspace(0, 48, 49)\n",
    "fprs = [fpr(threshold, 48) for threshold in thresholds]\n",
    "\n",
    "plt.plot(thresholds, fprs)\n",
    "plt.yscale(\"log\")\n",
    "plt.xlabel(\"Threshold\")\n",
    "plt.ylabel(\"FPR\")\n",
    "\n",
    "plt.minorticks_on()\n",
    "plt.grid(which='major', linestyle='-', linewidth=0.2, alpha=0.8)\n",
    "plt.grid(axis='both', which='minor', linestyle='-', linewidth=0.1, alpha=0.5)\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ldm",
   "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.10.10"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "df2bc3fed92e5aa2486479205013f8a4acaa6b99dec94d8ee399d56842a5d582"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}