{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import arviz as az\n",
    "import numpy as np\n",
    "\n",
    "import os, sys\n",
    "sys.path.append(os.path.join(\"../../\"))\n",
    "\n",
    "import privugger as pv\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "# Governor case study\n",
    "\n",
    "In this tutorial, we explore the famous case study presented by [Sweeny 2002](https://dl.acm.org/doi/10.1142/S0218488502001648) where she re-identified the the medical records of the governor of Masachussets. The problem the dataset was that records were \"anonymized\" simply by removing identifiers. Nowadays, it is well known that this is a naïve and insufficient method for anonymization data. Yet many real applications still use this type of \"anonymization\". \n",
    "\n",
    "This tutorial shows how privugger can be used analyze the uniqueness of records after identifiers have been removed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Naïve anonymization program (removing identifiers)\n",
    "\n",
    "We consider a program that takes a dataset with attributes: `name`, `zip`, `age`, `gender` and `diagnosis`. In this dataset, `name` is considered an identifier (this is an unrealistic assumption, since in reality there can be many people with the same name, but it serves the purpose of the tutorial). \n",
    "\n",
    "As you can see below, the program simply returns a dataset without names."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def remove_names(names, zips, ages, genders, diagnoses):\n",
    "    output = []\n",
    "    output.append(zips)\n",
    "    output.append(ages)\n",
    "    output.append(genders)\n",
    "    output.append(diagnoses)\n",
    "    return output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Input specification\n",
    "\n",
    "Note that this specification is _compositional_, we define different parts of the dataset and concatenate them as the bottom.\n",
    "\n",
    "We consider a dataset with 100 records. At this point, we fix the attributes of the victim's record (from now on called Alice). Alice is 25 years old, lives in ZIP code 50 (this is just a nominal placeholder for a real zip code), she is female and she is ill. These values are defined as constants of the form `ALICE_XXX` below.\n",
    "\n",
    "We define a constant record (using point distributions) that we may add to the specification to ensure that Alice's record is in the dataset. This step is optional. In fact, when the analysis is performed without this record, we model an adversary who doesn't know whether Alice is in the dataset.\n",
    "\n",
    "The variables `xxx_others` contain the spec for the records in the dataset. We have 5 different names uniformly distributed, ages are normally distributed with the most likely age being around 55. We consider 100 different ZIP codes uniformly distributed. Two genders (though acknowledging that this is a simplification and more genders exist) which are uniformly distributed. As for the diagnosis, there is a 20% chance of being ill.\n",
    "\n",
    "The following lines show how to merge Alice's record with the rest of the dataset. This also exemplifies how to compositionally define input specifications.\n",
    "\n",
    "Finally, we create the `Dataset` object with the complete input spec."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "N = 100\n",
    "\n",
    "# Victim's record\n",
    "\n",
    "## Constant values for Alice's record\n",
    "ALICE_NAME      = 0\n",
    "ALICE_AGE       = 25.0\n",
    "ALICE_ZIP       = 50\n",
    "ALICE_GENDER    = 1\n",
    "ALICE_DIAGNOSIS = 1\n",
    "\n",
    "## Spec for Alice's record (all point distributions with the values above)\n",
    "alice_name      = pv.Constant(\"alice_name\",ALICE_NAME, num_elements=1)\n",
    "alice_age       = pv.Constant(\"alice_age\", ALICE_AGE, num_elements=1)\n",
    "alice_zip       = pv.Constant(\"alice_zip\",ALICE_ZIP, num_elements=1)\n",
    "alice_gender    = pv.Constant(\"alice_gender\",ALICE_GENDER, num_elements=1) # 1: female, 0: male\n",
    "alice_diagnosis = pv.Constant(\"alice_diagnosis\", ALICE_DIAGNOSIS, num_elements=1) # 0: healthy, 1: ill\n",
    "\n",
    "# Spec for the records of others    \n",
    "names_others     = pv.DiscreteUniform(\"names_others\", 0, 5, num_elements=N)\n",
    "ages_others      = pv.Normal(\"ages_others\", mu=55.2, std=3.5, num_elements=N)\n",
    "zips_others      = pv.DiscreteUniform(\"zips_others\", 0, 100, num_elements=N)\n",
    "genders_others   = pv.Bernoulli(\"genders_others\",p=.5,num_elements=N)   \n",
    "diagnoses_others = pv.Bernoulli(\"diagnoses_others\", p=.2,num_elements=N)\n",
    "\n",
    "# Merging all in a single dataset spec\n",
    "names     = pv.concatenate(names_others, alice_name,  \"discrete\")\n",
    "zips      = pv.concatenate(zips_others, alice_zip, \"discrete\")\n",
    "ages      = pv.concatenate(ages_others, alice_age,  \"continuous\")\n",
    "genders   = pv.concatenate(genders_others, alice_gender, \"discrete\")\n",
    "diagnoses = pv.concatenate(diagnoses_others, alice_diagnosis, \"discrete\")\n",
    "\n",
    "# Dataset spec\n",
    "ds = pv.Dataset(input_specs = [names, zips, ages, genders, diagnoses])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Program specification\n",
    "\n",
    "The program specification takes the input specification above, and the program to analyze `remove_names`. We give the name `'output'` to the distribution of the output of the program. In this example, it is important to remark that the output of the program is a matrix of floats (i.e., a numeric dataset). In this matrix, each row models a row in the dataset, and each column models an attribute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "program   = pv.Program('output', dataset=ds, output_type=pv.Matrix(pv.Float), function=remove_names)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inference\n",
    "\n",
    "We use the pymc3 backend to perform the inference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (2 chains in 4 jobs)\n",
      "CompoundStep\n",
      ">CompoundStep\n",
      ">>Metropolis: [alice_diagnosis]\n",
      ">>Metropolis: [alice_gender]\n",
      ">>Metropolis: [alice_age]\n",
      ">>Metropolis: [alice_zip]\n",
      ">>Metropolis: [zips_others]\n",
      ">>Metropolis: [alice_name]\n",
      ">>Metropolis: [names_others]\n",
      ">NUTS: [ages_others]\n",
      ">BinaryGibbsMetropolis: [genders_others, diagnoses_others]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "        <style>\n",
       "            /* Turns off some styling */\n",
       "            progress {\n",
       "                /* gets rid of default border in Firefox and Opera. */\n",
       "                border: none;\n",
       "                /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "                background-size: auto;\n",
       "            }\n",
       "            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "                background: #F44336;\n",
       "            }\n",
       "        </style>\n",
       "      <progress value='22000' class='' max='22000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [22000/22000 03:21<00:00 Sampling 2 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 2 chains for 1_000 tune and 10_000 draw iterations (2_000 + 20_000 draws total) took 202 seconds.\n",
      "/home/pardo/.local/lib/python3.8/site-packages/arviz/stats/diagnostics.py:561: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  (between_chain_variance / within_chain_variance + num_samples - 1) / (num_samples)\n",
      "/home/pardo/.local/lib/python3.8/site-packages/xarray/core/nputils.py:227: RuntimeWarning: All-NaN slice encountered\n",
      "  result = getattr(npmodule, name)(values, axis=axis, **kwargs)\n",
      "/home/pardo/.local/lib/python3.8/site-packages/xarray/core/nputils.py:227: RuntimeWarning: All-NaN slice encountered\n",
      "  result = getattr(npmodule, name)(values, axis=axis, **kwargs)\n",
      "/home/pardo/.local/lib/python3.8/site-packages/xarray/core/nputils.py:227: RuntimeWarning: All-NaN slice encountered\n",
      "  result = getattr(npmodule, name)(values, axis=axis, **kwargs)\n",
      "/home/pardo/.local/lib/python3.8/site-packages/xarray/core/nputils.py:227: RuntimeWarning: All-NaN slice encountered\n",
      "  result = getattr(npmodule, name)(values, axis=axis, **kwargs)\n",
      "/home/pardo/.local/lib/python3.8/site-packages/xarray/core/nputils.py:227: RuntimeWarning: All-NaN slice encountered\n",
      "  result = getattr(npmodule, name)(values, axis=axis, **kwargs)\n",
      "The rhat statistic is larger than 1.4 for some parameters. The sampler did not converge.\n",
      "The estimated number of effective samples is smaller than 200 for some parameters.\n"
     ]
    }
   ],
   "source": [
    "trace = pv.infer(program, cores=4, draws=10000, method='pymc3', return_model=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Privacy Risk Analysis\n",
    "\n",
    "In this tutorial, we focus on uniqueness queries. That is, we quantify how unique Alice's record is based on her attribute values.\n",
    "\n",
    "### How many records share the Alice's attribute values?\n",
    "\n",
    "First we compute the average of the number of records that have the same values as Alice's record.\n",
    "\n",
    "Currentely, privugger does not have built-in functions to compute these type of query so the information in the trace directly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Avg. number of rows with gov's name:  16.20235\n",
      "Avg. number of rows with gov's zip:  1.95685\n",
      "Avg. number of rows with gov's age:  1.0\n",
      "Avg. number of rows with gov's gender:  50.98105\n",
      "Avg. number of rows with gov's diagnosis:  21.01775\n"
     ]
    }
   ],
   "source": [
    "trace_attr   = lambda attr : np.concatenate([trace.posterior['output'][i][:,attr,:] for i in [0,1]])\n",
    "\n",
    "names_db     = np.concatenate([trace.posterior['names_others'].values[i] for i in [0,1]])\n",
    "zips_db      = trace_attr(0)\n",
    "ages_db      = trace_attr(1)\n",
    "genders_db   = trace_attr(2)\n",
    "diagnoses_db = trace_attr(3)\n",
    "\n",
    "\n",
    "print(\"Avg. number of rows with gov's name: \",sum([np.count_nonzero(names_db[i]==ALICE_NAME) for i in range(0,len(names_db))])/len(names_db))\n",
    "print(\"Avg. number of rows with gov's zip: \",sum([np.count_nonzero(zips_db[i]==ALICE_ZIP) for i in range(0,len(zips_db))])/len(zips_db))\n",
    "print(\"Avg. number of rows with gov's age: \",sum([np.count_nonzero(ages_db[i]==ALICE_AGE) for i in range(0,len(ages_db))])/len(ages_db)) # 1.0 is expected\n",
    "print(\"Avg. number of rows with gov's gender: \",sum([np.count_nonzero(genders_db[i]==ALICE_GENDER) for i in range(0,len(genders_db))])/len(genders_db))\n",
    "print(\"Avg. number of rows with gov's diagnosis: \",sum([np.count_nonzero(diagnoses_db[i]==ALICE_DIAGNOSIS) for i in range(0,len(diagnoses_db))])/len(diagnoses_db))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Attributed based uniqueness histograms\n",
    "\n",
    "A more interesting perspective is to plot histograms showing the probability of having $n$ records equal to Alice's attributes.\n",
    "\n",
    "Here we focus on the attributes `gender` and `zip`. The plots below show the probability of having $n \\in \\mathbb{N}$ number of records with: 1) the same gender as Alice, 2) the same zip code as Alice, 3) the same gender _and_ zip code as Alice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEWCAYAAABollyxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAem0lEQVR4nO3df5wdZWHv8c/XxAQqCjZsrSbRjTdRGxQDjfyoP2u0hlKJvQYFvSW2WOTa9LbF1sbelovUewvWglVRSw01jdZgY3/sq8RyW+OvKsYsiGDEtAumEsCyJAEFiRD59o951h6OZ7NnN2ez2Sff9+t1XjvzPM/MPHNm93tmZ+bMyDYREVGvx0x1ByIiYnIl6CMiKpegj4ioXII+IqJyCfqIiMol6CMiKpegj4ioXIK+IpLeJOndE5z2C5JO6HGXIsZF0vGSvthW9ieS/udU9akGCfppRNIOSQ9Kul/Sf0j6sKSjSt0s4PeBP+4w3VGSbi/DvyLpsrb6VwLftf2VlrJFkjZIGpb0HUn/Jum9kuaV+pdIeqT0pfV1aqn/jKS9kua3zPNlknZ0WJ/vSrpX0hclnS/pMS1tPizpobZlfLXU9UtyS/kOSWt68V53Ut7H+yV9skPdDkkv61D+Ekk728peIelzZb2HJX1W0hml7g2SftDhfX1KF/07S9IWSQ9IursMv1mSSn037+Wmtnl+RNJFLevSus13Svq4pOe1TePSh9blvLXUXSTp4VI2ss1PHZnW9k3AveV3csS7gN8rv+MxAQn66eeVto8CTgSW0oQ7wArgG7bv6DDNCcBIiP80cENb/fnA+pERSQuBLcCdwAm2nwA8H7gVeEHLdHfaPqrtdV1L/QPAH3SxPo8HngZcAvwusLatzTvblvHctvpjynuyEvgDSS8fY5kT9Wrg+8DLJf3kRGYgaSXw18BfAvOAJwEXAq3Bdl2H9/XOMeb7FuBPaT7of7LM93ya7dYakGO9lydL+pn9LOrO8l4/HjgF+AbweUnL2to9t20572ypu7rM41jg0+X9aPVR4E0jI7bvKss5Yz/9iv1I0E9TJdA/CTy7FJ0GfHaU5kuB61uGfxj0ZS/ppW3TXgR8wfYFtneW5d1t+922N4yjm+8Bzpb038ZqaPs+2wPAa4FVkp491jQd5jEIbAOWjHfaLq0CPgjcBPyP8U5c9qwvA/7Q9ofKOj9i+7O2f3WinZJ0NHAx8GbbG21/142v2H697e+PY3bvBP7vWI3K/HfavhD4EHDpePttex9NqM+V1NdS9RlgmaTZbWWnj3cZ0UjQT1PlkMjP81976s8Btre1WSvpXpo/wt8pw0uBL0raVpotAh4ZCfTiZcAnetDNO4A/B97e7QS2vwzsBF443oVJOoXmg29oP21uKocMOr3ev5/pnga8hCaYPgqcM97+Ac8E5gMbJzDt/pwKzAb+vgfzej/wjE6Hofbjb4ATJT1uPAsqOxnnALuAPSPlZSfmYZr3a8QtQPt/H9GlBP3083clsP+FZi/8/5XyY4Dvtja0fS6wANhB82/ym4EP2j7G9nGjTVfafntkRNLqEoT3S/rzlnZP6RCW7X/sfwS8UtJxdO9O4Mdbxn+7bRnr2trfI+lB4DqaoPq70WZs+/iy/p1eb95Pn34JuMn214ENwHEa/8nrOeXnXWO0O6VtfW8do/2xwD1lDxmAcuz7XjXnQF7U0nas9/JBmj36d3SzQsWdgGh+l0bc0LacV7TUvab8Dj8I/CqwsrXvxXfb5tc+HuOQoJ9+XlVC6Wm232z7wVK+h+a4KQCSzih/TDtpjn9/G1gHnFP+8JZ2mq7YBTx5ZMT2+2wfA7wbeGxLuzs7hOUDrTOyPQy8j+bQQrfmArtbxt/VtoxVbe2PBY4C3kKz1/1Yeu8cmj35kT3Oz9IcyhmPXeXnk/fbCr7Utr5jHfraBRwraeZIge2fKdtsF4/+Ox/rvYTmUMyT2k6I7s9cwMC9LWUnti3n2pa6j5e+PQn4Gs15o3aPb5tf+3iMQ4K+HjcBzxgZsT1Q/pjWA28ow7uBvvKHN1iaDtEcPp7bMq9PAf+9h337Y+Bn6fwH/SjlCo65NP+xdM32D2xfBuyl+c9ltPlv63BFy8jrg6NM8zM0h7jeJunbkr4NnAy8rjVcu7AduJ3mpG4vXUdzknhFL2Zm+yGaw21/SLOnPpZfBG5o/5DvYjn3AOcBF0n64Ydf+V2cxaMPRf4U8NXxzD/+S4K+HpuAF3co/2maf6MXAHfZ3ttaWf6o/7lt2ouAF0q6bOQDQNKxNH9s42b7XuBPgLeO1kbSEyT9As1hkY/Yvnkiy6K5cuetko4YpS/HdbiiZeR1/ijzXAX8E7CY5kTvEppzAUfSnAQf8VhJR7S8HvUhYNvABTRXBv1yWefHSHqBpCsnuL4j7+/bgfdLWinp8WW+S4BxHTdvsR44AljeqVKNuZL+D/BG4PcmshDb24FrefTvxouBzW0nkV9Mc/FBTITtvKbJi+ZY+8tGqXss8C3gKW1lwzR7Za8G1o4y7enAJ9vKngl8HLiH5vjoduC9wPxS/xLgEeD+tterS/1ngDe2zO8o4G5gR9v6PFjmfx/NnumvATNa2nwYeKhtGfeUun6aQwYzW9qL5sqbX+/Re34EzeGtV3aoez+wsWVd3PZ6R3mfdrZNtxz4fFmX4fJenV7q3gD8oMP7+rwu+vp64MvA98p8t9DsMc+a4Hv5mlJ2UYdt/gDNsfmNwClt/XCpb13Ou0vdRTQf5K3tTy7tf6KMXwOc0VL/ZJpDkLOm+m9wur5U3siogKTzgMW2f3MC034BWO2WL01FHGySjgf+zPapLWV/Atxqe9SromL/EvQREZUbz4mkiJgikp4KfH2U6sW2v3Uw+xPTS/boIyIqd8jt0R977LHu7++f6m5EREwr119//T22+zrVHXJB39/fz+Dg4NgNIyLihyT9+2h1uY4+IqJyCfqIiMol6CMiKpegj4ioXII+IqJyCfqIiMol6CMiKpegj4ioXII+IqJyXQW9pOWStksakrSmQ/1sSVeX+i2S+kv5YyWtk3SzpFskva3H/Y+YMv1rrvnhK+JQNmbQS5oBXEHzJJ3FwNmSFrc1OxfYY3shcDlwaSk/E5ht+zk0Tzp608iHQEREHBzd7NGfBAzZvs3NY+c28KPPplxB8+BpaJ44s0ySaJ4087jySLUjaZ5u852e9DwiIrrSTdDPpXmg8YidpaxjG9v7aB4LN4cm9B8A7qJ5zN27bO9uX4Ck8yQNShocHh4e90pEHCw5VBPT0WSfjD2J5vmXTwEWAG+R9PT2RravtL3U9tK+vo532YyIiAnqJujvAOa3jM8rZR3blMM0RwO7gNcB/2j7Ydt3A18Alh5opyMionvdBP1WYJGkBZJmAWcBA21tBoBVZXglsNnNo6u+BbwUQNLjgFOAb/Si4xER0Z0xg74cc18NXAvcAnzc9jZJF0s6ozRbC8yRNARcAIxcgnkFcJSkbTQfGH9h+6Zer0RERIyuqydM2d4EbGoru7BleC/NpZTt093fqTwiIg6efDM2IqJyCfqIiModcg8Hj5iOWq+t33HJ6VPYk4gflT36iIjKJegjIiqXoI+IqFyCPiKicgn6iIjKJegjIiqXoI+IqFyuo4/osVxTH4ea7NFHjCEPG4npLkEfEVG5BH1EROUS9BERlUvQR0RUrqugl7Rc0nZJQ5LWdKifLenqUr9FUn8pf72kG1tej0ha0ttViIiI/Rkz6CXNoHkk4GnAYuBsSYvbmp0L7LG9ELgcuBTA9kdtL7G9BPgl4Ju2b+xd9yMiYizd7NGfBAzZvs32Q8AGYEVbmxXAujK8EVgmSW1tzi7TRkTEQdRN0M8Fbm8Z31nKOrYpDxO/D5jT1ua1wMc6LUDSeZIGJQ0ODw930++IiOjSQTkZK+lk4Hu2v9ap3vaVtpfaXtrX13cwuhQRcdjoJujvAOa3jM8rZR3bSJoJHA3saqk/i1H25iMiYnJ1E/RbgUWSFkiaRRPaA21tBoBVZXglsNm2ASQ9BngNOT4fETElxrypme19klYD1wIzgKtsb5N0MTBoewBYC6yXNATspvkwGPEi4Hbbt/W++xERMZau7l5pexOwqa3swpbhvcCZo0z7GeCUiXcxIiIORL4ZGxFRuQR9RETlEvQREZVL0EdEVC5BHxFRuQR9RETlEvQRHeQ5sVGTBH1EROUS9BERlevqm7ERMTGth392XHL6FPYkDmfZo4+IqFyCPiKicgn6iIjKJegjIiqXoI+IqFyCPiKicgn6iIjKdRX0kpZL2i5pSNKaDvWzJV1d6rdI6m+pO17SdZK2SbpZ0hE97H9ERIxhzKCXNAO4AjgNWAycLWlxW7NzgT22FwKXA5eWaWcCHwHOt30c8BLg4Z71PiIixtTNHv1JwJDt22w/BGwAVrS1WQGsK8MbgWWSBPwccJPtrwLY3mX7B73pekREdKOboJ8L3N4yvrOUdWxjex9wHzAHeAZgSddKukHSWzstQNJ5kgYlDQ4PD493HSIiYj8m+2TsTOAFwOvLz1+UtKy9ke0rbS+1vbSvr2+SuxQRcXjpJujvAOa3jM8rZR3blOPyRwO7aPb+P2f7HtvfAzYBJx5opyMionvdBP1WYJGkBZJmAWcBA21tBoBVZXglsNm2gWuB50j6sfIB8GLg673pesT0MvIwkzzQJA62MW9TbHufpNU0oT0DuMr2NkkXA4O2B4C1wHpJQ8Bumg8DbO+RdBnNh4WBTbbzWx6HpJEAzu2EozZd3Y/e9iaawy6tZRe2DO8Fzhxl2o/QXGIZERFTIN+MjYioXII+IqJyCfqIiMol6CMiKpegj4ioXII+IqJyCfqIiMol6CMiKpegj4ioXII+IqJyCfqIiMol6CMiKpegj4ioXII+IqJyCfqIiMol6CMiKtdV0EtaLmm7pCFJazrUz5Z0danfIqm/lPdLelDSjeX1wR73PyIixjDmE6YkzQCuAF5O87DvrZIGbLc++/VcYI/thZLOAi4FXlvqbrW9pLfdjoiIbnWzR38SMGT7NtsPARuAFW1tVgDryvBGYJkk9a6bERExUd0E/Vzg9pbxnaWsYxvb+4D7gDmlboGkr0j6rKQXdlqApPMkDUoaHB4eHtcKRETE/k32ydi7gKfaPgG4APgrSU9ob2T7SttLbS/t6+ub5C5FRBxeugn6O4D5LePzSlnHNpJmAkcDu2x/3/YuANvXA7cCzzjQTkdERPfGPBkLbAUWSVpAE+hnAa9razMArAKuA1YCm21bUh+w2/YPJD0dWATc1rPeR0xT/Wuu+eHwjktOn8KexOFgzKC3vU/SauBaYAZwle1tki4GBm0PAGuB9ZKGgN00HwYALwIulvQw8Ahwvu3dk7EiERHRWTd79NjeBGxqK7uwZXgvcGaH6T4BfOIA+xgREQcg34yNw1r/mmsedRglokYJ+oiIyiXoIyIql6CPiKhcgj4ionIJ+oiIyiXoIyIql6CPiKhcgj4ionIJ+oiIyiXoIyIql6CPiKhcgj4ionIJ+oiIyiXoIyIql6CPiKhcgj4ionJdBb2k5ZK2SxqStKZD/WxJV5f6LZL62+qfKul+Sb/do35HRESXxgx6STOAK4DTgMXA2ZIWtzU7F9hjeyFwOXBpW/1lwCcPvLsRETFe3ezRnwQM2b7N9kPABmBFW5sVwLoyvBFYJkkAkl4FfBPY1pMeR0TEuHQT9HOB21vGd5ayjm1s7wPuA+ZIOgr4XeDt+1uApPMkDUoaHB4e7rbvERHRhck+GXsRcLnt+/fXyPaVtpfaXtrX1zfJXYqIOLzM7KLNHcD8lvF5paxTm52SZgJHA7uAk4GVkt4JHAM8Immv7fcdaMcjIqI73QT9VmCRpAU0gX4W8Lq2NgPAKuA6YCWw2baBF440kHQRcH9CPiLi4Boz6G3vk7QauBaYAVxle5uki4FB2wPAWmC9pCFgN82HQcQhqX/NNQDsuOT0Ke5JY6Q/cOj0KerSzR49tjcBm9rKLmwZ3gucOcY8LppA/yIi4gDlm7EREZVL0EdEVC5BHxFRuQR9RETlEvQREZVL0EdEVC5BHxFRuQR9RETlEvQREZVL0EdEVC5BHxFRuQR9RETlEvQREZVL0EdEVC5BHxFRuQR9RETlugp6ScslbZc0JGlNh/rZkq4u9Vsk9ZfykyTdWF5flfSLPe5/RFf611zzqCc5RRxOxgx6STOAK4DTgMXA2ZIWtzU7F9hjeyFwOXBpKf8asNT2EmA58Gfl4eEREXGQdLNHfxIwZPs22w8BG4AVbW1WAOvK8EZgmSTZ/p7tfaX8CMC96HRERHSvm6CfC9zeMr6zlHVsU4L9PmAOgKSTJW0DbgbObwn+H5J0nqRBSYPDw8PjX4uIiBjVpJ+Mtb3F9nHA84C3STqiQ5srbS+1vbSvr2+yuxQRcVjp5nj5HcD8lvF5paxTm53lGPzRwK7WBrZvkXQ/8GxgcMI9jqhY6wnjHZecPoU9iZp0s0e/FVgkaYGkWcBZwEBbmwFgVRleCWy27TLNTABJTwOeBezoSc8jIqIrY+7R294naTVwLTADuMr2NkkXA4O2B4C1wHpJQ8Bumg8DgBcAayQ9DDwCvNn2PZOxIhER0VlXlzra3gRsaiu7sGV4L3Bmh+nWA+sPsI8REXEA8s3YiIjKJegjIiqXoI+IqFyCPiKicgn6iIjKJegjIiqXoI+IqFyCPiKicgn6iIjKJegjIiqXoI+IqFyCPiKicnl+a8QhKvemj17JHn1Uq3/NNY8Ky4jDVYI+IqJyCfqIiMol6CMiKtdV0EtaLmm7pCFJazrUz5Z0danfIqm/lL9c0vWSbi4/X9rj/kdExBjGDHpJM4ArgNOAxcDZkha3NTsX2GN7IXA5cGkpvwd4pe3n0Dw8PI8VjIg4yLrZoz8JGLJ9m+2HgA3AirY2K4B1ZXgjsEySbH/F9p2lfBtwpKTZveh4RER0p5ugnwvc3jK+s5R1bGN7H3AfMKetzauBG2x/v30Bks6TNChpcHh4uNu+R0REFw7KyVhJx9EcznlTp3rbV9peantpX1/fwehSRMRho5ugvwOY3zI+r5R1bCNpJnA0sKuMzwP+FjjH9q0H2uGIiBifboJ+K7BI0gJJs4CzgIG2NgM0J1sBVgKbbVvSMcA1wBrbX+hRnyMiYhzGDPpyzH01cC1wC/Bx29skXSzpjNJsLTBH0hBwATByCeZqYCFwoaQby+sner4WERExqq5uamZ7E7CprezCluG9wJkdpnsH8I4D7GNERByAfDM2IqJyCfqIiMol6CMiKpegj4ioXII+IqJyCfqIiMol6KMqeXxgxI9K0EdEVC5BHzFN5L+VmKgEfURE5RL0ERGVS9BHRFQuQR8RUbkEfURE5RL0ERGVS9BHRFSuq6CXtFzSdklDktZ0qJ8t6epSv0VSfymfI+nTku6X9L4e9z0CyPXlEWMZM+glzQCuAE4DFgNnS1rc1uxcYI/thcDlwKWlfC/wB8Bv96zHERExLt3s0Z8EDNm+zfZDwAZgRVubFcC6MrwRWCZJth+w/S80gR8RPZL/YmI8ugn6ucDtLeM7S1nHNuVh4vcBc7rthKTzJA1KGhweHu52soiI6MIhcTLW9pW2l9pe2tfXN9XdiYioSjdBfwcwv2V8Xinr2EbSTOBoYFcvOhgREQemm6DfCiyStEDSLOAsYKCtzQCwqgyvBDbbdu+6GREREzVzrAa290laDVwLzACusr1N0sXAoO0BYC2wXtIQsJvmwwAASTuAJwCzJL0K+DnbX+/5mkREREdjBj2A7U3AprayC1uG9wJnjjJt/wH0LyIiDtAhcTI2IiImT4I+YprLNfUxlgR9TEsJt4juJegjIiqXoI+IqFyCPiKicl1dXhkx1VqPx++45PQp7EnE9JM9+oiIyiXoIyIql6CPqEguO41OEvQREZVL0EdEVC5X3cQhK4cgInoje/QRlcrx+hiRoI+IqFyCPuIwkL37w1uO0cchJWEU0Xtd7dFLWi5pu6QhSWs61M+WdHWp3yKpv6XubaV8u6RX9LDvETEB2bs//Iy5Ry9pBnAF8HJgJ7BV0kDbc1/PBfbYXijpLOBS4LWSFtM8P/Y44CnAP0t6hu0f9HpFYnpJ0EQcPN0cujkJGLJ9G4CkDcAKoDXoVwAXleGNwPskqZRvsP194Jvl4eEnAdf1pvsxnSTcDz0j22THJac/ajjqItv7byCtBJbbfmMZ/yXgZNurW9p8rbTZWcZvBU6mCf8v2f5IKV8LfNL2xrZlnAecV0afCWw/wPU6FrjnAOcx3WSdDw9Z58PDRNb5abb7OlUcEidjbV8JXNmr+UkatL20V/ObDrLOh4es8+Gh1+vczcnYO4D5LePzSlnHNpJmAkcDu7qcNiIiJlE3Qb8VWCRpgaRZNCdXB9raDACryvBKYLObY0IDwFnlqpwFwCLgy73pekREdGPMQze290laDVwLzACusr1N0sXAoO0BYC2wvpxs3U3zYUBp93GaE7f7gF87SFfc9Oww0DSSdT48ZJ0PDz1d5zFPxkZExPSWWyBERFQuQR8RUbmqgn6sWzXUQNJ8SZ+W9HVJ2yT9Rin/cUn/JOnfys8nTnVfe03SDElfkfQPZXxBueXGULkFx6yp7mMvSTpG0kZJ35B0i6RTa9/Okn6r/F5/TdLHJB1R23aWdJWku8v3j0bKOm5XNd5T1v0mSSdOZJnVBH3LrRpOAxYDZ5dbMNRmH/AW24uBU4BfK+u5BviU7UXAp8p4bX4DuKVl/FLgctsLgT00t+KoyZ8C/2j7WcBzada92u0saS7wv4Cltp9Nc/HHyC1VatrOHwaWt5WNtl1Po7lacRHNl0o/MJEFVhP0tNyqwfZDwMitGqpi+y7bN5Th79L88c+lWdd1pdk64FVT0sFJImkecDrwoTIu4KU0t9yAytZZ0tHAi2iuaMP2Q7bvpfLtTHMl4JHl+zg/BtxFZdvZ9udork5sNdp2XQH8pRtfAo6R9OTxLrOmoJ8L3N4yvrOUVavcJfQEYAvwJNt3lapvA0+aqn5NkncDbwUeKeNzgHtt7yvjtW3vBcAw8BflcNWHJD2Oirez7TuAdwHfogn4+4DrqXs7jxhtu/Yk12oK+sOKpKOATwC/afs7rXXly2rVXDcr6ReAu21fP9V9OYhmAicCH7B9AvAAbYdpKtzOT6TZg11Ac7fbx/GjhziqNxnbtaagP2xutyDpsTQh/1Hbf1OK/2PkX7ry8+6p6t8keD5whqQdNIfkXkpz/PqY8i8+1Le9dwI7bW8p4xtpgr/m7fwy4Ju2h20/DPwNzbaveTuPGG279iTXagr6bm7VMO2VY9NrgVtsX9ZS1XobilXA3x/svk0W22+zPc92P8123Wz79cCnaW65AfWt87eB2yU9sxQto/mGebXbmeaQzSmSfqz8no+sc7XbucVo23UAOKdcfXMKcF/LIZ7u2a7mBfw88K/ArcD/nur+TNI6voDm37qbgBvL6+dpjll/Cvg34J+BH5/qvk7S+r8E+Icy/HSaeycNAX8NzJ7q/vV4XZcAg2Vb/x3wxNq3M/B24BvA14D1wOzatjPwMZpzEA/T/Od27mjbFRDN1YS3AjfTXJE07mXmFggREZWr6dBNRER0kKCPiKhcgj4ionIJ+oiIyiXoIyIql6CPaU3SH0n6WUmvkvS2cU7bV+6K+BVJL5ysPo6jP5+RdFg9BDsOjgR9THcnA18CXgx8bpzTLgNutn2C7c+P1qjcGbWnWr7pGTHpEvQxLUn6Y0k3Ac8DrgPeCHxA0oUd2vZL2lzu5/0pSU+VtAR4J7BC0o2SjmybZoekSyXdAJwp6eckXSfpBkl/Xe41hKTnSfqipK9K+rKkx5d7qP+FpJvLfws/W9q+QdKApM3ApyQdKWlDudf83wJHlnYzJH243JP9Zkm/NXnvZBwOslcR05Lt31Hz4PlzgAuAz9h+/ijN3wuss71O0q8A77H9qvKhsNT26lGm22X7REnH0tx35WW2H5D0u8AFki4BrgZea3urpCcAD9LcN9+2nyPpWcD/l/SMMs8TgeNt75Z0AfA92z8l6XjghtJmCTDXzT3ZkXTMBN+mCCB79DG9nQh8FXgWj34gSbtTgb8qw+tpbiPRjavLz1NoHmbzBUk30tyL5GnAM4G7bG8FsP0dN7fTfQHwkVL2DeDfgZGg/yfbI/cif1FLu5tobnUAcBvwdEnvlbQceNTdSSPGK3v0Me2Uwy4fprmT3z00D6hQCeFTbT/Yo0U9MLJImoA+u60fzzmAeY7K9h5JzwVeAZwPvAb4lQksKwLIHn1MQ7ZvtL2E5gZ2i4HNwCtsLxkl5L9Ic9dLgNcDo554HcWXgOdLWggg6XHlUMx24MmSnlfKH19Osn6+LIfS7qmlbbvPAa8r7Z4NHF+GjwUeY/sTwO/T/OcSMWHZo49pSVIfsMf2I5KeZfvr+2n+6zRPavodmqc2/fJ4lmV7WNIbgI9Jml2Kf9/2v0p6LfDecjL3QZp7qr+f5sTwzTTP+H2D7e83d959lA+Uft1Cc+hp5MEqc0v5yI7YuC4bjWiXu1dGRFQuh24iIiqXoI+IqFyCPiKicgn6iIjKJegjIiqXoI+IqFyCPiKicv8JRSJUKAW8BsUAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEWCAYAAABv+EDhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZFUlEQVR4nO3df5BdZYHm8e8zwaADOqD0TmlCTNCgRsGATdRC0VF+BFkJO/4KjjNBncpgkV1dalzjauFUXHcBd6jZtVDJllHGESPK6LRjKAZB/IVoGohggtEmk5FkUCKJIopgzLN/3DdyuN5O39t9O53O+3yqbvU573nfc99zcnKee37cc2WbiIiozx9MdQciImJqJAAiIiqVAIiIqFQCICKiUgmAiIhKJQAiIiqVAIiIqFQCIPpG0l9J+rtxtv2mpBP63KXokaRDJX1f0kCj7NWSPjOV/YrJkQCIrknaKukhSQ9K+omkT0g6vEybCbwX+GCHdodLuqcMv0XSZW3TXw38wvbtZfzBDq/fSNpSpp8n6Rvd9KvPyy9JWyRt6jDtJkl/2aF8riRLOqRRtkjSOkk/k7RT0nckvblMe7mkPR2W/8X76NdLR1lneyStKXU+Iel/tPVpb72tklYC2H4YWAOs3Dt/218Enivp+HGvvDggJQCiV6+2fThwIjBIa6cPsAT4vu3tHdqcANxehl8A3NY2/Xzgk3tHbB/efAHHAjuB94+jX/10CvAfgGMknTSeGZQd+Y3AV4FnAk8B3gac2aj27+3rwPa3Rpun7a93WGd/CvwSuGy0dsARpe65wEWSFpfyq4Blkg5t1P00sLzHxY0DXAIgxqXs6K8FnleKzqS1U+tkELi1Mfy7AChHDq8YrW355Hw18EXbHx9Hv/ppGfBPwLoyPB4fBK60fYntn7rlVtuv71cnJR0NfAp4m+3vjVW/hMtGyjqzvQ3YBbyoUe0m4Kx+9TEODAmAGJeyk3kVj36yPw7Y3FbnY5J+BlwCvLMMDwI3S9pYqs0H9pSdTieXAocBK8bZr/bp/1xOvXR6/fM+5vuHwGtp7Vg/BSwt4dW1Mo8XA5/rpV2P7zET+CzwOdv/0EV9SToZeC6PXWd3Ac9vG58r6Un97G9MrUPGrhLxGF+QtBv4OfAl4H+W8iOAXzQr2n6rpL8Gvg0cT+u0xMm2L2hU+712e0l6DfBm4AW2fz3Ofj2G7f84xnxG86fAw8C/0Pp/8zhan4g/38M8jqT1oeveMeo9rYRl0yzbv+ziPf4WmAG8o4u6PwUM/BhYafuGxrRf0Pq3aY5Tyh7oYt4xDSQAolfn2P5yh/JdwBP3jkg6G/h7WjvKQ2jtZA4Dfi3pz4BTbQ+3t2u0Pxb4GPAXtrdMoF/9sgy42vZuYLeka0pZLwGwC9gDPBX4/j7q/bvt2b12UNJS4I3AieVi7liOKsvTyROBn7WN01YW01xOAUW/3EHrYi0AtodsH0Hr4u55ZXgnMGD7iLLzBxihdSZi1t625VTJNcBHbQ/1s5OSrh3ljpkHJV07SpvZtK5TvEnSjyX9mNbpoFdJOqrb97b9K+BbwGv6sSxtfXwOsBr4c9v/1odZPgf4btv4Vtv59H8QSQBEv6wDXtah/AXAbZLmAfe2n8qx/Qjw5ba2HwXuB97T707aPrPDHTZ7X2eO0uzPgR8AzwIWltexwDZad9DsdYikxzdej+swr/8GnCfpnZKeAiDp+ZLWjneZJB1GKzD/j+11451PY36zgCcDtzSKX0br4nocRBIA0S9fBJ4t6Wl7C8oOcC6tneeJPHonULsraO1kkTSnDL8I+Hn7p/RJ7P++LAM+bPvHzRetoGreDfQR4KHG6/fuWrJ9M62jiVcAWyTtpPXJvbnjflqHo5N9HTW8htYn9Au7PaoZwxtp3anUPI10Lq1/pziIKL8IFv0iaTmwwPY7xtH2m8CKvV8Gi6lR7v3/LnCK7ftK2atpnVrq262qcWBIAEREVCp3AUVME/s4BXam7a/v187EQSFHABERlTrgjgCOOuooz507d6q7ERExrdx6660/tT0wds1HHXABMHfuXIaHh8euGBERvyOp5+9/5DbQiIhKJQAiIiqVAIiIqFQCICKiUgmAiIhKJQAiIiqVAIiIqFQCICKiUgmAiIhKHXDfBJ4Mc1d+qat6Wy8+a5J7EhFx4MgRQEREpRIAERGVquIU0HjktFFEHOxyBBARUakEQEREpRIAERGVSgBERFQqARARUakEQEREpRIAERGVSgBERFSqqwCQtFjSZkkjklZ2mH6+pDslbZD0DUkLSvlcSQ+V8g2SPtrvBYiIiPEZ85vAkmYAlwOnAduA9ZKGbG9qVLvK9kdL/bOBy4DFZdrdthf2tdcRETFh3RwBLAJGbG+x/QiwFljSrGD7gcboYYD718WIiJgM3QTALOCexvi2UvYYki6QdDdwKfBfGpPmSbpd0lclvbTTG0haLmlY0vCOHTt66H5ERIxX3y4C277c9jOAdwHvLcX3AnNsnwBcCFwl6Ukd2q62PWh7cGBgoF9dioiIfegmALYDRzfGZ5ey0awFzgGw/bDt+8vwrcDdwLHj6mlERPRVNwGwHpgvaZ6kmcBSYKhZQdL8xuhZwA9L+UC5iIykY4D5wJZ+dDwiIiZmzLuAbO+WtAK4DpgBrLG9UdIqYNj2ELBC0qnAb4BdwLLS/BRglaTfAHuA823vnIwFiYiI3nT1gzC21wHr2souagy/fZR21wDXTKSDERExOfJN4IiISiUAIiIqlQCIiKhUAiAiolIJgIiISiUAIiIqlQCIiKhUAiAiolIJgIiISiUAIiIqlQCIiKhUAiAiolIJgIiISiUAIiIqlQCIiKhUAiAiolIJgIiISiUAIiIq1VUASFosabOkEUkrO0w/X9KdkjZI+oakBY1p7y7tNks6o5+dj4iI8RszACTNAC4HzgQWAOc2d/DFVbaPs70QuBS4rLRdACwFngssBj5c5hcREVOsmyOARcCI7S22HwHWAkuaFWw/0Bg9DHAZXgKstf2w7X8FRsr8IiJiih3SRZ1ZwD2N8W3AC9srSboAuBCYCbyi0faWtrazOrRdDiwHmDNnTjf9joiICerbRWDbl9t+BvAu4L09tl1te9D24MDAQL+6FBER+9BNAGwHjm6Mzy5lo1kLnDPOthERsZ90EwDrgfmS5kmaSeui7lCzgqT5jdGzgB+W4SFgqaRDJc0D5gPfmXi3IyJiosa8BmB7t6QVwHXADGCN7Y2SVgHDtoeAFZJOBX4D7AKWlbYbJV0NbAJ2AxfY/u0kLUtERPSgm4vA2F4HrGsru6gx/PZ9tP0A8IHxdjAiIiZHvgkcEVGpBEBERKUSABERlUoARERUKgEQEVGpBEBERKUSABERlUoARERUKgEQEVGpBEBERKUSABERlUoARERUKgEQEVGpBEBERKUSABERlUoARERUKgEQEVGpBEBERKUSABERleoqACQtlrRZ0oiklR2mXyhpk6Q7JN0g6emNab+VtKG8hvrZ+YiIGL8xfxRe0gzgcuA0YBuwXtKQ7U2NarcDg7Z/JeltwKXAG8q0h2wv7G+3IyJioro5AlgEjNjeYvsRYC2wpFnB9lds/6qM3gLM7m83IyKi37oJgFnAPY3xbaVsNG8Frm2MP17SsKRbJJ3TqYGk5aXO8I4dO7roUkRETNSYp4B6IelNwCDwskbx021vl3QMcKOkO23f3WxnezWwGmBwcND97FNERHTWzRHAduDoxvjsUvYYkk4F3gOcbfvhveW2t5e/W4CbgBMm0N+IiOiTbgJgPTBf0jxJM4GlwGPu5pF0AnAFrZ3/fY3yIyUdWoaPAk4GmhePIyJiiox5Csj2bkkrgOuAGcAa2xslrQKGbQ8BHwQOBz4rCeBHts8GngNcIWkPrbC5uO3uoYiImCJdXQOwvQ5Y11Z2UWP41FHa3QwcN5EORkTE5Mg3gSMiKpUAiIioVAIgIqJSCYCIiEolACIiKpUAiIioVAIgIqJSCYCIiEolACIiKpUAiIioVAIgIqJSCYCIiEolACIiKpUAiIioVAIgIqJSCYCIiEolACIiKpUAiIioVFcBIGmxpM2SRiSt7DD9QkmbJN0h6QZJT29MWybph+W1rJ+dj4iI8RszACTNAC4HzgQWAOdKWtBW7XZg0PbxwOeAS0vbJwPvA14ILALeJ+nI/nU/IiLGq5sjgEXAiO0tth8B1gJLmhVsf8X2r8roLcDsMnwGcL3tnbZ3AdcDi/vT9YiImIhuAmAWcE9jfFspG81bgWt7aStpuaRhScM7duzooksRETFRfb0ILOlNwCDwwV7a2V5te9D24MDAQD+7FBERo+gmALYDRzfGZ5eyx5B0KvAe4GzbD/fSNiIi9r9uAmA9MF/SPEkzgaXAULOCpBOAK2jt/O9rTLoOOF3SkeXi7+mlLCIiptghY1WwvVvSClo77hnAGtsbJa0Chm0P0TrlczjwWUkAP7J9tu2dkt5PK0QAVtneOSlLEhERPRkzAABsrwPWtZVd1Bg+dR9t1wBrxtvBiIiYHPkmcEREpbo6AojuzF35pa7qbb34rEnuSUTE2HIEEBFRqQRARESlEgAREZVKAEREVCoBEBFRqQRARESlEgAREZVKAEREVCoBEBFRqQRARESlEgAREZVKAEREVCoBEBFRqQRARESlEgAREZVKAEREVKqrAJC0WNJmSSOSVnaYfoqk2yTtlvTatmm/lbShvIba20ZExNQY8xfBJM0ALgdOA7YB6yUN2d7UqPYj4DzgrzvM4iHbCyfe1YiI6KdufhJyETBiewuApLXAEuB3AWB7a5m2ZxL6GBERk6CbU0CzgHsa49tKWbceL2lY0i2SzumlcxERMXn2x4/CP932dknHADdKutP23c0KkpYDywHmzJmzH7oUERHdHAFsB45ujM8uZV2xvb383QLcBJzQoc5q24O2BwcGBrqddURETEA3AbAemC9pnqSZwFKgq7t5JB0p6dAyfBRwMo1rBxERMXXGDADbu4EVwHXAXcDVtjdKWiXpbABJJ0naBrwOuELSxtL8OcCwpO8CXwEubrt7KCIipkhX1wBsrwPWtZVd1BheT+vUUHu7m4HjJtjHiIiYBPkmcEREpRIAERGVSgBERFQqARARUakEQEREpRIAERGVSgBERFQqARARUakEQEREpRIAERGVSgBERFQqARARUakEQEREpRIAERGVSgBERFQqARARUakEQEREpRIAERGVSgBERFSqqwCQtFjSZkkjklZ2mH6KpNsk7Zb02rZpyyT9sLyW9avjERExMWMGgKQZwOXAmcAC4FxJC9qq/Qg4D7iqre2TgfcBLwQWAe+TdOTEux0RERPVzRHAImDE9hbbjwBrgSXNCra32r4D2NPW9gzgets7be8CrgcW96HfERExQd0EwCzgnsb4tlLWja7aSlouaVjS8I4dO7qcdURETMQBcRHY9mrbg7YHBwYGpro7ERFV6CYAtgNHN8Znl7JuTKRtRERMom4CYD0wX9I8STOBpcBQl/O/Djhd0pHl4u/ppSwiIqbYmAFgezewgtaO+y7gatsbJa2SdDaApJMkbQNeB1whaWNpuxN4P60QWQ+sKmURETHFDummku11wLq2sosaw+tpnd7p1HYNsGYCfYyIiElwQFwEjoiI/S8BEBFRqQRARESlEgAREZVKAEREVCoBEBFRqQRARESlEgAREZVKAEREVCoBEBFRqQRARESlEgAREZVKAEREVCoBEBFRqQRARESlEgAREZVKAEREVCoBEBFRqa4CQNJiSZsljUha2WH6oZI+U6Z/W9LcUj5X0kOSNpTXR/vc/4iIGKcxfxNY0gzgcuA0YBuwXtKQ7U2Nam8Fdtl+pqSlwCXAG8q0u20v7G+3IyJioro5AlgEjNjeYvsRYC2wpK3OEuDKMvw54JWS1L9uRkREv3UTALOAexrj20pZxzq2dwM/B55Sps2TdLukr0p6aac3kLRc0rCk4R07dvS0ABERMT5jngKaoHuBObbvl/QC4AuSnmv7gWYl26uB1QCDg4Oe5D4dUOau/FJX9bZefNYk9yQiatPNEcB24OjG+OxS1rGOpEOAPwLut/2w7fsBbN8K3A0cO9FOR0TExHUTAOuB+ZLmSZoJLAWG2uoMAcvK8GuBG21b0kC5iIykY4D5wJb+dD0iIiZizFNAtndLWgFcB8wA1tjeKGkVMGx7CPgY8ElJI8BOWiEBcAqwStJvgD3A+bZ3TsaCREREb7q6BmB7HbCureyixvCvgdd1aHcNcM0E+xgREZMg3wSOiKhUAiAiolIJgIiISiUAIiIqlQCIiKhUAiAiolIJgIiISiUAIiIqlQCIiKhUAiAiolIJgIiISiUAIiIqlQCIiKhUAiAiolIJgIiISiUAIiIqNdk/Ch+TID8kHxH9kCOAiIhKJQAiIirVVQBIWixps6QRSSs7TD9U0mfK9G9LmtuY9u5SvlnSGX3se0RETMCY1wAkzQAuB04DtgHrJQ3Z3tSo9lZgl+1nSloKXAK8QdICYCnwXOBpwJclHWv7t/1ekBhdt9cMINcNImrSzRHAImDE9hbbjwBrgSVtdZYAV5bhzwGvlKRSvtb2w7b/FRgp84uIiCnWzV1As4B7GuPbgBeOVsf2bkk/B55Sym9pazur/Q0kLQeWl9EHJW3uqvePOgr4aY9tfo8umdZtfm8d7K++HSD6sg1MY7UvP2QdPKvXBgfEbaC2VwOrx9te0rDtwT52adqpfR1k+etefsg6kDTca5tuTgFtB45ujM8uZR3rSDoE+CPg/i7bRkTEFOgmANYD8yXNkzST1kXdobY6Q8CyMvxa4EbbLuVLy11C84D5wHf60/WIiJiIMU8BlXP6K4DrgBnAGtsbJa0Chm0PAR8DPilpBNhJKyQo9a4GNgG7gQsm6Q6gcZ8+OojUvg6y/FH7Ouh5+dX6oB4REbXJN4EjIiqVAIiIqNS0D4CxHlNxsJO0VdKdkjaM5zaw6UjSGkn3Sfpeo+zJkq6X9MPy98ip7ONkGmX5/0bS9rIdbJD0qqns42SSdLSkr0jaJGmjpLeX8pq2gdHWQU/bwbS+BlAeU/EDGo+pAM5te0zFQU3SVmDQdjVfgJF0CvAg8Pe2n1fKLgV22r64fBA40va7prKfk2WU5f8b4EHb/3sq+7Y/SHoq8FTbt0l6InArcA5wHvVsA6Otg9fTw3Yw3Y8AunlMRRxkbH+N1t1mTc3HkVxJ6z/DQWmU5a+G7Xtt31aGfwHcResJAzVtA6Otg55M9wDo9JiKnlfCNGfgXyTdWh6pUas/tn1vGf4x8MdT2ZkpskLSHeUU0UF7+qOpPHn4BODbVLoNtK0D6GE7mO4BEPAS2ycCZwIXlNMDVStfQpy+5zbH5yPAM4CFwL3A305pb/YDSYcD1wDvsP1Ac1ot20CHddDTdjDdA6D6R03Y3l7+3gd8nnqftvqTcl507/nR+6a4P/uV7Z/Y/q3tPcD/4yDfDiQ9jtaO71O2/7EUV7UNdFoHvW4H0z0AunlMxUFL0mHlAhCSDgNOB76371YHrebjSJYB/zSFfdnv9u74iv/EQbwdlEfNfwy4y/ZljUnVbAOjrYNet4NpfRcQQLnN6e949DEVH5jaHu0/ko6h9akfWo/1uKqG5Zf0aeDltB7/+xPgfcAXgKuBOcC/Aa+3fVBeKB1l+V9O67DfwFbgrxrnww8qkl4CfB24E9hTiv87rXPgtWwDo62Dc+lhO5j2ARAREeMz3U8BRUTEOCUAIiIqlQCIiKhUAiAiolIJgIiISiUAYlqT9L8k/YmkcyS9u8e2A5K+Lel2SS+drD720J+bJFX7o+ax/yUAYrp7IXAL8DLgaz22fSVwp+0TbH99tErlqbN9JWnMn2ONmGwJgJiWJH1Q0h3AScC3gL8EPiLpog5150q6sTwg6wZJcyQtBC4FlpTnpj+hrc1WSZdIug14naTTJX1L0m2SPluewYKkkyTdLOm7kr4j6YmSHi/p4+V3Gm6X9Cel7nmShiTdCNwg6QmS1kq6S9LngSeUejMkfULS98o8/uvkrcmoWT6FxLRk+52Srgb+ArgQuMn2yaNU/xBwpe0rJb0F+L+2zylhMWh7xSjt7rd9oqSjgH8ETrX9S0nvAi6UdDHwGeANttdLehLwEPD2Vhd9nKRn03pa67FlnicCx9veKelC4Fe2nyPpeOC2UmchMKvxrP8jxrmaIvYpRwAxnZ0IfBd4Nq3noY/mxcBVZfiTwEu6nP9nyt8XAQuAb0raQOs5M08HngXca3s9gO0HbO8u8/+HUvZ9Wo8l2BsA1zceT3BKo94dwB2lfAtwjKQPSVoMPOZJlxH9kiOAmHbK6ZtP0Hr660+BP2wVawPwYtsP9emtfrn3LWntuM9t68dxE5jnqGzvkvR84AzgfFq/8vSWcbxXxD7lCCCmHdsbbC+k9XOgC4AbgTNsLxxl538zrSfFAvwZrYdo9eIW4GRJz4TfPYX1WGAz8FRJJ5XyJ5aLu18v70OpN6fUbfc14I2l3vOA48vwUcAf2L4GeC+tI52IvssRQExLkgaAXbb3SHr2GL8D/Z+Bj0t6J7ADeHMv72V7h6TzgE9LOrQUv9f2DyS9AfhQuYj8EHAq8GFaF6TvBHYD59l+uPUE38f4SOnXXbROYd1aymeV8r0f0Hq6vTWiW3kaaEREpXIKKCKiUgmAiIhKJQAiIiqVAIiIqFQCICKiUgmAiIhKJQAiIir1/wFEDjmoEiwDEgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEWCAYAAAB2X2wCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAcUklEQVR4nO3de5gcZZ328e9NIiAHAcnoahJIhHCInB0QxAVXQRIOCZcgggrEw2ZdN57igmF10WXfXTns8roqi8sqwgtoQBDNaiCwIKuiYAbklMRozEYTDjKcFBCFyO/943kGik6fZtKTnnm4P9fV13RVPVX16+6au6ufqq5WRGBmZqPfRt0uwMzMOsOBbmZWCAe6mVkhHOhmZoVwoJuZFcKBbmZWCAe6mVkhigx0SX8l6XNDnPdmSXt3uCQze5GQtImkn0nqqYw7StLlw77yiBh1N2AV8BTwBPAb4CJgizxtY2A1ML7OfFsAq/P99wLn1kw/Cri2MvxEndszwMo8fRbww3bq6vDjF7ASWFpn2k3A++uMnwQEMLYybj9gIfAY8AjwE+A9edqbgGfrPP4DmtT15w2es2eBC3Obi4D/U1PTQLtVwLwWj/0E4HsDr2OLtlsC5+blPgn8GrgSeH2lTeRp1XpPzdM+k6cfV2k/No+bVHk8TwOP59s9wGeBrSrzzAL+VOd5eXWd7eaB4dpuaraRR4FNasY/99rUmSeAHSvDOwHfAB4CfgvcBcwFxtR5XQdu72ijtl7gO7m+x4ClwD8B2wziuXwQ2LyyzPcDNzV4zR8GbqitLT9Hf6hZx3/V+d94HFhO/r+pzH8q8K814+4B9hiu1zUiRvUe+lERsQWwD2kj+FQePxP4WUTcW2eevYGf5vuvA26vmf4B4JKBgYjYonojbcSPAP84hLo66SDgFcBrJO07lAVIOgC4EfgfYEdgW+CvgemVZvfVPgcR8eNGy4yIH9R5zt5G+uc5t0k5W+e2JwCnS5rWpO0RwFeARyTt2eTxbZIf3+7AkcDLgF2B+TWPEWDPmrrPrkx7BPgHSWOa1HR2RGwJ9ADvAfYHbpa0eaXNj+s8l/dVpg9sN3uRttPTmqxvyCRNIr3xBjBjiMvYAbiVtOO0e0RsBbydtL1vWWm6dc3jbbqHKukNpCC9GdglIrYGpgFrgepr3eq5HAN8pMXD2DM/3zuT3si+KOnTNW3m1KzjqMq0+/L8LwM+BvynpJ0r078GnJy3wwFfB2a3qGu9jOZAByAH9zXAbnnUdFJI1dML3Fa5/1ygS9oYeHOjeSWNBa4gvUt/dQh1ddLJwLdJe9cnD3EZ5wAXR8RZEfFQJLdFxHGdKlLSROAy4K8j4p5W7fObxRIaPGeSNgIOBa4lPfbDmyzuRGACcHRE3BMRf4qIJyPiyoj4zCAexrWkPfB3t1H/HyJiMSkotyWF+6BExAPAIlKwD4eTgFtIITbUbecfgB9FxNyIuB8gIpZHxDsj4rH1qO1s4KsR8dmI+E1e7q8j4tMRcdMglnMO8LeStm7VMG/7l5B2Zk6TtO1gCs7/NwtJb/x7VMavIX3K2L/S/CbSDsmwGfWBnkPjcJ7f896d9BGo2uYrkh4DzgJOyfd7gR9JWpKbTQGezS9EPWcDmwNzhlhX7fTvSHqswe07TZa7GXAsKSgvA47Pb0Zty8s4gNT9MCxyTd8AroyIS9toL0kHAq+lwXNG6iJaGREPAd+leaAfAiyKiCcHV/k6Avh74NOSXtLWDBGPA9eT9oQHRdIE0k7JiiZt/r3JtnNXi1WcxPPbzmGSXjnYGknPbUe3nfxp5gDgqg4sro8Unn87iHm+TepO228wK5K0kaQZwDjWfc2W8cJPFsuASZJeNph1DMZoDvRv5WD+IWmv+p/z+K1J/VrPiYj3AZNJ/WvjgA8CX4qIrSPitY3mGyDpGNLe1jER8Ych1vUCEXFkXn+925FNlv824I/AdaRQewmDf9ffhvTa39+i3avrBMbmLeYZ8K+kj74fbaPtQ6Q9nC+T+tBvaNDuCNKeOcCPgZ0lbdOg7ThSfzQAkvbK9f9O0vKatrfXPMbDqhMjYgHQT+qLbdd9wMsrw/vXrOOXNe2/JelxUjfGg0Dtx/9qPR9ssu3s0Wg+SW8EtgeuiIjbgF8C7xzEYxqwLa23HYCHah7zrk3aDmyT1dfs7Dzfk5KqXZetnkuA04EPVQ9MNhMRz5C2w+pr9vma9VS7Wl+d/8+fAq4G5kZE7Y7I46RcqQ5TM66jRnOgH5034O3zBv5UHv8olX48STPyE7+GtDE/AFwMnJRfpN5681Xm34nUZ3tyRKxcj7o65WTSP+Ta/OZyFYP/6Pwo6aDOq1q0u69OYLTc45V0PCkojo2IP7ZRz7iI2CYido2Izzdpdzg50CPiT6Q+8sMatH2YyuOLiDtyn+zbgE1q2u5T8xgX1Vnep4BPApu28XgAxpPepAbcUrOOHWraHx2pH/5NwC6kN6ROOxm4Ln/CgdzPO4TlvOC5bWJczWNe1qTtOttkRJyaX7OrSXvPA1o9l+Quvu8A89qok/zpq4cXvmYfrlnP31em3ZdrexnweVJ3ba0tSQd2q8PUjOuo0RzojdxFOngJpL2r/MRfAszK9x8BevKL1JebriB98h8/MG/umriKtDe/oJNFSrpG0hMNbtc0mGcCacN5t6QHJD1A6n45XFLbARARvyft4R7TicdSU+OuwAXAiRHxqw4u989I/+zVA9nNul1uAN46iE8UTUXE9aRt5IOt2kragtQt8YMhrOd/SP3b/9Jk+V9qsu0saTDPS4HjgIMr287HgD3V5OByA/9Nh7edvKNwK+kNt1M+Dfwl6c21lZmkg68/GcwK8g7LJ4DdJR1dM3lX4M6a4VUR8bvBrGMwSgz0hcDBdca/jvTRejJwf23XSUQ8TdpQq/N+ibQ38slOFxkR02PdI/UDt9qzMAacCPycdGR+r3zbifTp44RKu7GSNq3c6vX9ngrMknTKwIEgSXtKmj/Ux5TD8yrg3/KBok6aTjqltHoB/2tJoV1vO/5/pG6BqyXtJmmMpE1Jx06G6pOk560upfOPXwd8i7TH2fLgeQOfAw5tFLQR8YEm285r680DHE063W8qz287u5LedE6qtBtTs+3UOz7zaeANks7Jb7RI2lHSpWrjQGQTpwLvlTRP0ivycieQuksHLSJWAJcDH27URtLLJb0LOA84KyIeHsJ6niZ1MZ5eWe54UvfNLZWmB5NOlBg+MYznRA7XjdQXfkiDaS8hnW/86ppx/aTzt48BvtJg3iOAa/L97UgHxGrPRX0CeCKePye29jz0unV16HH/DPhQnfGnAn3x/PmzUXO7lMbnoV9DOo/4EdIe0kmx7rm21dsxTeo7ifrndT9ReV4vYt3z0Me28divJHXh1I6/lcp55TXTtiKF469yTb8iveHsV2lTr97P5WmfAS6tWeZCGp+H/gTpLJ2zSKfsDcwzi/rnTu/baLsBzgeu6uC2cy0150Xn8ceRuiHH5sdSu+38sPI8Vc9D35l00PvhvP3cSTpe0uw89Llt1Pl6nv9uxGOkc7f/Cdh2KM8lMJH0P3xTg9f8EdL3Gt5ZU8dNrPu/f1vlf2NNTfvNSH3wR+XhU1j3ey53k06XHLZsVF5RUSTNBqZGxEeHMO/NpPNPG51pYRuY0imjDwCviZqPq0rnDm8UEQ0PIpptSErnnt8JHBQRD+ZxR5G6ITt2WnDddZcY6FaW/PH7mIg4v860HUjfXm15aqRZ6RzoNmiSnmgwaXpEDPpAoL145IO229eZ9FcRcdmGrqc0DnQzs0KMbd1keIwbNy4mTZrUrdWbmY1Kt91220MRUfcLU10L9EmTJtHX19e6oZmZPUdSw+93lHgeupnZi5ID3cysEA50M7NCONDNzArhQDczK4QD3cysEA50M7NCONDNzArhQDczK0TXvik6Ukya9931mn/VmcP6I95mZm3zHrqZWSHaCnRJ0yQtl7RCUt0fXZV0nKSlkpZI+lpnyzQzs1ZadrlIGkP6vb1DSb9duVjSgohYWmkzBTgNODAiHh34PUAzM9tw2tlD3w9YERErI/0Y6nzSL2RX/SVwXkQ8CjDws0tmZrbhtBPo44HVleE1eVzVTsBOkm6WdIukafUWJGm2pD5Jff39/UOr2MzM6urUQdGxwBTSr2GfAPynpK1rG0XEBRHRGxG9PT11r89uZmZD1E6g3wtMrAxPyOOq1gALIuKZiPhf4OekgDczsw2knUBfDEyRNFnSxsDxwIKaNt8i7Z0jaRypC2Zl58o0M7NWWgZ6RKwF5gCLgGXAFRGxRNIZkmbkZouAhyUtBb4HnBIRDw9X0WZmtq62vikaEQuBhTXjTq/cD2BuvpmZWRf4m6JmZoVwoJuZFcKBbmZWCAe6mVkhHOhmZoVwoJuZFcKBbmZWCAe6mVkhHOhmZoVwoJuZFcKBbmZWCAe6mVkhHOhmZoVwoJuZFcKBbmZWCAe6mVkhHOhmZoVwoJuZFcKBbmZWCAe6mVkhHOhmZoVwoJuZFcKBbmZWCAe6mVkh2gp0SdMkLZe0QtK8OtNnSeqXdEe+vb/zpZqZWTNjWzWQNAY4DzgUWAMslrQgIpbWNL08IuYMQ41mZtaGdvbQ9wNWRMTKiHgamA/MHN6yzMxssNoJ9PHA6srwmjyu1jGS7pJ0paSJ9RYkabakPkl9/f39QyjXzMwa6dRB0f8CJkXEHsD1wMX1GkXEBRHRGxG9PT09HVq1mZlBe4F+L1Dd456Qxz0nIh6OiD/mwS8Dr+tMeWZm1q52An0xMEXSZEkbA8cDC6oNJL2qMjgDWNa5Es3MrB0tz3KJiLWS5gCLgDHAhRGxRNIZQF9ELAA+LGkGsBZ4BJg1jDWbmVkdLQMdICIWAgtrxp1euX8acFpnSzMzs8HwN0XNzArhQDczK4QD3cysEA50M7NCONDNzArhQDczK4QD3cysEA50M7NCONDNzArhQDczK4QD3cysEA50M7NCONDNzArhQDczK4QD3cysEA50M7NCONDNzArhQDczK4QD3cysEA50M7NCONDNzArhQDczK4QD3cysEA50M7NCONDNzArRVqBLmiZpuaQVkuY1aXeMpJDU27kSzcysHS0DXdIY4DxgOjAVOEHS1DrttgQ+Atza6SLNzKy1dvbQ9wNWRMTKiHgamA/MrNPuH4GzgD90sD4zM2tTO4E+HlhdGV6Txz1H0j7AxIj4brMFSZotqU9SX39//6CLNTOzxtb7oKikjYBzgY+3ahsRF0REb0T09vT0rO+qzcysop1AvxeYWBmekMcN2BLYDbhJ0ipgf2CBD4yamW1Y7QT6YmCKpMmSNgaOBxYMTIyI30bEuIiYFBGTgFuAGRHRNywVm5lZXS0DPSLWAnOARcAy4IqIWCLpDEkzhrtAMzNrz9h2GkXEQmBhzbjTG7R90/qXZWZmg+VvipqZFcKBbmZWCAe6mVkhHOhmZoVwoJuZFcKBbmZWCAe6mVkhHOhmZoVwoJuZFcKBbmZWCAe6mVkhHOhmZoVwoJuZFcKBbmZWCAe6mVkhHOhmZoVwoJuZFcKBbmZWCAe6mVkhHOhmZoVwoJuZFcKBbmZWCAe6mVkhHOhmZoVwoJuZFaKtQJc0TdJySSskzasz/QOS7pZ0h6QfSpra+VLNzKyZloEuaQxwHjAdmAqcUCewvxYRu0fEXsDZwLmdLtTMzJob20ab/YAVEbESQNJ8YCawdKBBRPyu0n5zIDpZ5Ggyad5313sZq848ogOVmNmLTTuBPh5YXRleA7y+tpGkvwHmAhsDb663IEmzgdkA22233WBrNTOzJjp2UDQizouIHYBPAJ9q0OaCiOiNiN6enp5OrdrMzGgv0O8FJlaGJ+RxjcwHjl6PmszMbAjaCfTFwBRJkyVtDBwPLKg2kDSlMngE8IvOlWhmZu1o2YceEWslzQEWAWOACyNiiaQzgL6IWADMkXQI8AzwKHDycBZtZmbrauegKBGxEFhYM+70yv2PdLguMzMbJH9T1MysEA50M7NCONDNzArhQDczK4QD3cysEA50M7NCONDNzArhQDczK4QD3cysEA50M7NCONDNzArhQDczK4QD3cysEA50M7NCONDNzArhQDczK4QD3cysEA50M7NCONDNzArhQDczK4QD3cysEA50M7NCONDNzArhQDczK0RbgS5pmqTlklZImldn+lxJSyXdJekGSdt3vlQzM2umZaBLGgOcB0wHpgInSJpa0+ynQG9E7AFcCZzd6ULNzKy5dvbQ9wNWRMTKiHgamA/MrDaIiO9FxO/z4C3AhM6WaWZmrbQT6OOB1ZXhNXlcI+8Drqk3QdJsSX2S+vr7+9uv0szMWuroQVFJ7wZ6gXPqTY+ICyKiNyJ6e3p6OrlqM7MXvbFttLkXmFgZnpDHvYCkQ4BPAgdHxB87U56ZmbWrnT30xcAUSZMlbQwcDyyoNpC0N/AfwIyIeLDzZZqZWSstAz0i1gJzgEXAMuCKiFgi6QxJM3Kzc4AtgG9IukPSggaLMzOzYdJOlwsRsRBYWDPu9Mr9Qzpcl5mZDZK/KWpmVggHuplZIRzoZmaFcKCbmRXCgW5mVggHuplZIRzoZmaFcKCbmRXCgW5mVggHuplZIRzoZmaFcKCbmRXCgW5mVggHuplZIRzoZmaFcKCbmRXCgW5mVggHuplZIRzoZmaFcKCbmRXCgW5mVggHuplZIRzoZmaFcKCbmRXCgW5mVoi2Al3SNEnLJa2QNK/O9IMk3S5praRjO1+mmZm10jLQJY0BzgOmA1OBEyRNrWn2a2AW8LVOF2hmZu0Z20ab/YAVEbESQNJ8YCawdKBBRKzK054dhhrNzKwN7XS5jAdWV4bX5HGDJmm2pD5Jff39/UNZhJmZNbBBD4pGxAUR0RsRvT09PRty1WZmxWsn0O8FJlaGJ+RxZmY2grQT6IuBKZImS9oYOB5YMLxlmZnZYLUM9IhYC8wBFgHLgCsiYomkMyTNAJC0r6Q1wNuB/5C0ZDiLNjOzdbVzlgsRsRBYWDPu9Mr9xaSuGDMz6xJ/U9TMrBAOdDOzQjjQzcwK4UA3MytEWwdFrbsmzfvuei9j1ZlHdKASMxvJvIduZlYIB7qZWSEc6GZmhXCgm5kVwoFuZlYIB7qZWSEc6GZmhXCgm5kVwoFuZlYIB7qZWSEc6GZmhXCgm5kVwoFuZlYIB7qZWSEc6GZmhXCgm5kVwoFuZlYIB7qZWSEc6GZmhfBvir5I+XdKzcrT1h66pGmSlktaIWlenembSLo8T79V0qSOV2pmZk21DHRJY4DzgOnAVOAESVNrmr0PeDQidgT+L3BWpws1M7Pm2uly2Q9YERErASTNB2YCSyttZgKfyfevBL4oSRERHazVRjh345h1l1plrqRjgWkR8f48fCLw+oiYU2lzT26zJg//Mrd5qGZZs4HZeXBnYHmnHkg2DnioZavuc52dNRrqHA01guvstOGoc/uI6Kk3YYMeFI2IC4ALhmv5kvoione4lt8prrOzRkOdo6FGcJ2dtqHrbOeg6L3AxMrwhDyubhtJY4GtgIc7UaCZmbWnnUBfDEyRNFnSxsDxwIKaNguAk/P9Y4Eb3X9uZrZhtexyiYi1kuYAi4AxwIURsUTSGUBfRCwAvgJcImkF8Agp9Lth2LpzOsx1dtZoqHM01Aius9M2aJ0tD4qamdno4K/+m5kVwoFuZlaIIgK91aUJRgJJEyV9T9JSSUskfaTbNTUjaYykn0r6TrdraUTS1pKulPQzScskHdDtmuqR9LH8mt8j6euSNu12TQCSLpT0YP4eycC4l0u6XtIv8t9tulljrqlenefk1/0uSVdL2rqLJQ7UtE6dlWkflxSSxg1nDaM+0Nu8NMFIsBb4eERMBfYH/maE1jngI8CybhfRwr8B10bELsCejMB6JY0HPgz0RsRupBMLunXSQK2LgGk14+YBN0TEFOCGPNxtF7FundcDu0XEHsDPgdM2dFF1XMS6dSJpIvBW4NfDXcCoD3QqlyaIiKeBgUsTjCgRcX9E3J7vP04Kn/Hdrao+SROAI4Avd7uWRiRtBRxEOsOKiHg6Ih7ralGNjQVemr+jsRlwX5frASAivk86K61qJnBxvn8xcPSGrKmeenVGxHURsTYP3kL6fkxXNXg+IV3f6lRg2M9AKSHQxwOrK8NrGKFBOSBfjXJv4NYul9LI50gb4LNdrqOZyUA/8NXcNfRlSZt3u6haEXEv8C+kvbP7gd9GxHXdraqpV0bE/fn+A8Aru1lMm94LXNPtIuqRNBO4NyLu3BDrKyHQRxVJWwBXAR+NiN91u55ako4EHoyI27pdSwtjgX2A8yNib+BJRkb3wAvkPuiZpDegVwObS3p3d6tqT/5y4Ig+r1nSJ0ndmZd1u5ZakjYD/g44fUOts4RAb+fSBCOCpJeQwvyyiPhmt+tp4EBghqRVpO6rN0u6tLsl1bUGWBMRA59yriQF/EhzCPC/EdEfEc8A3wTe0OWamvmNpFcB5L8PdrmehiTNAo4E3jVCv5m+A+mN/M78/zQBuF3Snw3XCksI9HYuTdB1kkTq710WEed2u55GIuK0iJgQEZNIz+WNETHi9igj4gFgtaSd86i38MJLOo8Uvwb2l7RZ3gbewgg8eFtRvYzHycC3u1hLQ5KmkboFZ0TE77tdTz0RcXdEvCIiJuX/pzXAPnnbHRajPtDzgZGBSxMsA66IiCXdraquA4ETSXu8d+Tb4d0uapT7EHCZpLuAvYB/7m4568qfIK4EbgfuJv3PjYivrUv6OvBjYGdJayS9DzgTOFTSL0ifLs7sZo3QsM4vAlsC1+f/pS91tUga1rlhaxiZn1TMzGywRv0eupmZJQ50M7NCONDNzArhQDczK4QD3cysEA50GxUkfVbSX0g6WtKgLsQkqUfSrfkSAX8+XDUOop6bJI34Hzi20ceBbqPF60kXYToY+P4g530LcHdE7B0RP2jUKF+5s6PyBbnMNggHuo1o+brXdwH7kr608X7gfEnrXB9D0iRJN+ZrZN8gaTtJewFnAzPzF1BeWjPPKklnSbodeLukt0r6saTbJX0jX3sHSftK+pGkOyX9RNKWkjaV9FVJd+e9/7/IbWdJWiDpRuAGSS+VNF/pmu1XAy/N7cZIukjpOul3S/rY8D2T9mLgvQcb0SLiFElXACcBc4GbIuLABs2/AFwcERdLei/w+Yg4Ood/b0TMaTDfwxGxT/7xgW8Ch0TEk5I+AcyVdCZwOfCOiFgs6WXAU6RrxkdE7C5pF+A6STvlZe4D7BERj0iaC/w+InaVtAfpW6OQvt06Pl8nnZHwIw02unkP3UaDfYA7gV1ofh2UA4Cv5fuXAG9sc/mX57/7k34k5WZJd5CuZbI9sDNwf0QsBoiI3+VLTrwRuDSP+xnwK2Ag0K+PiIFrYx9UaXcXcFcevxJ4jaQv5GuTjLirb9ro4j10G7Fyd8lFpKvUPUT6cQjlsD0gIp7q0KqeHFglKYhPqKlj9/VYZkMR8aikPYHDgA8Ax5Gu7W02JN5DtxErIu6IiL1IPzE2FbgROCwi9moQ5j/i+Z93exfQ8ABoA7cAB0raEUDS5rkLZTnwKkn75vFb5oOdP8jrIbfbLret9X3gnbndbsAe+f44YKOIuAr4FCPz8r82ingP3UY0ST3AoxHxrKRdIqLZJXI/RPoFo1NIv2b0nsGsKyL68zW2vy5pkzz6UxHxc0nvAL6QD6o+RboS4b+TDtDeTfqRhVkR8cd0ldwXOD/XtYzUZTTw4yHj8/iBHauR8LuYNor5aotmZoVwl4uZWSEc6GZmhXCgm5kVwoFuZlYIB7qZWSEc6GZmhXCgm5kV4v8DdtlhBCdcsVkAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "total_samples=len(genders_db)\n",
    "total_N=N\n",
    "\n",
    "count_gov_genders=np.array([np.count_nonzero(\n",
    "    (genders_db[i]==ALICE_GENDER)\n",
    ") for i in range(0,total_samples)])\n",
    "\n",
    "count_gov_zips=np.array([np.count_nonzero(\n",
    "    (zips_db[i]==ALICE_ZIP)\n",
    ") for i in range(0,total_samples)])\n",
    "\n",
    "count_gov_genders_zips=np.array([np.count_nonzero(\n",
    "    (zips_db[i]==ALICE_ZIP) &\n",
    "    (genders_db[i]==ALICE_GENDER)\n",
    ") for i in range(0,total_samples)])\n",
    "\n",
    "\n",
    "y1=[np.count_nonzero(count_gov_genders==i)/len(count_gov_genders) for i in range(0,total_N)]\n",
    "y2=[np.count_nonzero(count_gov_zips==i)/len(count_gov_zips) for i in range(0,total_N)]\n",
    "y3=[np.count_nonzero(count_gov_genders_zips==i)/len(count_gov_genders_zips) for i in range(0,total_N)]\n",
    "plt.bar(range(0,total_N),y1)\n",
    "plt.title(\"P(#(GENDER = ALICE_GENDER))\")\n",
    "plt.xlabel(\"# of records\")\n",
    "plt.show()\n",
    "plt.bar(range(0,total_N),y2)\n",
    "plt.title(\"P(#(ZIP = ALICE_ZIP))\")\n",
    "plt.xlim((-1,25))\n",
    "plt.xlabel(\"# of records\")\n",
    "plt.show()\n",
    "plt.bar(range(0,total_N),y3)\n",
    "plt.title(\"P(#(ZIP = ALICE_ZIP /\\ GENDER = ALICE_GENDER))\")\n",
    "plt.xlim((-1,15))\n",
    "plt.xlabel(\"# of records\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Discussion of results\n",
    "\n",
    "The results show that gender or zip code in isolation have a low probability of uniquely identifying Alice in the dataset. However for both of them combined we observe that the probability of having only one record with those attribute values is close to 60%. This level of risk may be unacceptable for some types of sensitive data."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "argv": [
    "python",
    "-m",
    "ipykernel_launcher",
    "-f",
    "{connection_file}"
   ],
   "display_name": "Python 3",
   "env": null,
   "interrupt_mode": "signal",
   "language": "python",
   "metadata": null,
   "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.10"
  },
  "name": "Governor.ipynb"
 },
 "nbformat": 4,
 "nbformat_minor": 4
}