{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "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",
    "# move to previouse directory to access the privugger code\n",
    "import os, sys\n",
    "sys.path.append(os.path.join(\"../../\"))\n",
    "\n",
    "import privugger as pv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example of using Privugger on OpenDP "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial shows the use of privugger on a program using the differential privacy libray [OpenDP](https://github.com/opendp).\n",
    "\n",
    "## OpenDP program\n",
    "\n",
    "We consider a program that takes as input a dataset with attributes: age, sex, education, race, income and marriage status. The program outputs the mean of the incomes and adds Laplacian noise to protects the individuals privacy.\n",
    "\n",
    "For each attribute, the program takes a parameter of type array (int or float) with the attribute value for each individual in the dataset. For example, to model a dataset of size 2 where the first individual is 20 and the second is 40, we set the `age` input parameter as `age=[20,40]`. The remaining parameters are defined in the same way. The last parameter `N` indicates the number of records in the dataset.\n",
    "\n",
    "This way of defining the input may seem unnatural, but, as we will see below, it allows for a structured manner to specify the prior of the program.\n",
    "\n",
    "Furthermore, note that the first lines of `dp_program` simply defined a pandas dataframe. This snippet of code can be adapted to other programs. The part of the code after the comment `## After here the...` can contain arbitrary code working on a pandas dataframe with the attributes defined in the parameters of the program."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dp_program(age, sex, educ, race, income, married, N):\n",
    "    import opendp.smartnoise.core as sn\n",
    "    import pandas as pd\n",
    "\n",
    "    # assert that all vectors are have the same size\n",
    "    assert age.size == sex.size == educ.size == race.size == income.size == married.size\n",
    "    \n",
    "    ## Dataframe definition (can be automatized)\n",
    "    temp_file='temp.csv'\n",
    "    var_names = [\"age\", \"sex\", \"educ\", \"race\", \"income\", \"married\"]\n",
    "    data = {\n",
    "        \"age\":     age,\n",
    "        \"sex\":     sex,\n",
    "        \"educ\":    educ,\n",
    "        \"race\":    race,\n",
    "        \"income\":  income,\n",
    "        \"married\": married\n",
    "    }\n",
    "    df = pd.DataFrame(data,columns=var_names)\n",
    "    \n",
    "    ## After here the program works on a pandas dataframe\n",
    "    df.to_csv(temp_file)\n",
    "    with sn.Analysis() as analysis:\n",
    "        # load data\n",
    "        data = sn.Dataset(path=temp_file,column_names=var_names)\n",
    "\n",
    "        # get mean income with laplacian noise (epsilon=.1 arbitrarily chosen)\n",
    "        age_mean = sn.dp_mean(data = sn.to_float(data['income']),\n",
    "                              privacy_usage = {'epsilon': .1},\n",
    "                              data_lower = 0., # min income\n",
    "                              data_upper = 200., # max income                   \n",
    "                              data_rows = N\n",
    "                             )\n",
    "    analysis.release()\n",
    "    return np.float64(age_mean.value)   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Input specification\n",
    "\n",
    "The next step is to specify the prior knowledge of the attacker using probability distributions (i.e., they are defined as _random variables_). In this example, we show specify each of the attributes that conform the dataset.\n",
    "\n",
    "The variable `N` defines the size of the dataset (`N_rv` is a point distribution with all probability mass concentrated at `N`, this is necessary because the input specification must be composed by random variables). In this example, we consider a dataset of size 150.\n",
    "\n",
    "For non-numeric attributes such as `sex`, `educ` or `race` we use distribution over natural numbers with each number denoting a category. We treat them as nominal values (i.e., we assume there is not order relation, $\\leq$ among them). For these attributes (and `married`) we specify a uniform distribution over all possible categories.\n",
    "\n",
    "For `age`, we set a binomial distribution prior with support 0 to 120; this distribution gives highest probability to ages close to 60 years old. We remak here that this prior may be refined by using statistical data about age data.\n",
    "\n",
    "Finally, `income` is distributed according to a Normal distribution with mean 100 and standard deviation 5. This gives high probability to values close to 100 (i.e., 100k DKK)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "N         = 150\n",
    "N_rv      = pv.Constant('N', N)\n",
    "age       = pv.Binomial('age', p=0.5, n=120, num_elements=N)\n",
    "sex       = pv.DiscreteUniform('sex', 0,2,num_elements=N)\n",
    "educ      = pv.DiscreteUniform('educ', 0,10, num_elements=N)\n",
    "race      = pv.DiscreteUniform('race', 0,50, num_elements=N)\n",
    "income    = pv.Normal('income', mu=100, std=5, num_elements=N)\n",
    "married   = pv.DiscreteUniform('married', 0,1,num_elements=N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dataset & Program specification  \n",
    "\n",
    "The input specification above is always wrapped into a `Dataset` object. This object is used as the input for the program."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds       = pv.Dataset(input_specs = [age, sex, educ, race, income, married, N_rv])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Program specification\n",
    "\n",
    "The program specification takes the `Dataset` above and the program to analyze. To this end, we define a `Program` object including the `Dataset`, a python function corresponding to the program to analyze (the input parameters of the function must match those of the `Dataset`). Finally, it is necessary to specify the type of the output of the program. In this case, since we are analyzing a program compute the mean income it is a float. The first parameter of the `Program` constructor is the name of the output distribution (i.e., the distribution of the output of the program under analysis).\n",
    "\n",
    " output type specifies the output type of the program, in this case it is a floating point number as the program calulates the mean. The function is the program specified above. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "program = pv.Program('output',dataset=ds, output_type=pv.Float, function=dp_program)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inference\n",
    "\n",
    "Lastly we use the privug interface to perform the inference. This is done by calling `infer` and specifying the `program`, number of cores, number of chains, number of draws, and the backend (which is pymc3 in this example)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (2 chains in 4 jobs)\n",
      "CompoundStep\n",
      ">CompoundStep\n",
      ">>Metropolis: [N]\n",
      ">>Metropolis: [married]\n",
      ">>Metropolis: [race]\n",
      ">>Metropolis: [educ]\n",
      ">>Metropolis: [sex]\n",
      ">>Metropolis: [age]\n",
      ">NUTS: [income]\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='42000' class='' max='42000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [42000/42000 02:59<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 20_000 draw iterations (2_000 + 40_000 draws total) took 180 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",
      "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, chains=2, draws=20000, method='pymc3')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Privacy Risk Analysis\n",
    "\n",
    "### Mutual Information\n",
    "\n",
    "To quantify privacy risks, in this tutorial we are using mutual information (though we remark here that privugger can be used to compute many more leakage measures).\n",
    "\n",
    "We study the risks for the individual in the first record. We compute the mutual information between the output of the program (mean income + laplace noise) and each of the other attributes in the dataset.\n",
    "\n",
    "Since the mutual informaiton estimator we use is not exact, we compute 100 estimates and use box plots to get an impresison of the accuracy of the estimates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEWCAYAAABbgYH9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA2XklEQVR4nO3dd5wU9f3H8df7juOOrsIpAipEsYAiiYg10SD2KCkW1FhRYq/52WIiGruxl6ioUdSIiKJYUaOgUZGiKF1RMIANEZB2B3f3+f0x39Pl3LvbPXb2yn6ej8c9bnfmO9/5TNn5zMx3iswM55xzLtPy6jsA55xzTZMnGOecc7HwBOOccy4WnmCcc87FwhOMc865WHiCcc45FwtPMFkmqaskk9Ssmv7TJe2dYl3bSJoiabmkszMZZyZIukfSX2Ood4ikRzNdb0MR13xrCiSNlXRyhupqlOtR4vohaW9JC+o7puo06QQjaZ6kNZI6VOn+QdjId02xHpO0VSxBVmFmPc1sbIrFLwTeMLM2ZnZ7jGHVStIJkv6b2M3MTjWzv9dXTMk0ho1KQ5xvjV1D2xAn29FM9htKJpPrR9zbtiadYIK5wFGVXyTtALSsv3Ayagtgel0GrO4IyjmXvmz9niTlZ2M8GWNmTfYPmAdcBkxM6PYP4C+AAV1Dt7HAyQllTgD+Gz6/GcquBFYARyb2TxjGgK3C54OBD4DvgfnAkIRyXUPZZjXE3D98HgKMAIYBy4mSSZ/Q73WgHCgJcW0NtAtlFwGfh2nPS5imt4FbgMXAVcBDwN3AS6GOt4GOwK3AEmAW8POE2C4GPg2xzAB+F7pvF+IoD/UsDd0fAq5KGP4UYA7wHTAa6FRl/p0KfAIsBe4CVM08GgKMBJ4IsbwP7JjQvxPwVJgPc4GzQ/cDgDXA2hDnh8CvgakJw77KuuvLW8Bva6o39MtLmD+Lw3LbqMoyPx74H/At8Jca1tsf5huwN7AAuAD4BvgSODGhbAvgprC8lwH/BVqEfocSrTNLidbx7aqsZ/8HfES0bj8AbBLWheXAa8CGCeV3Bd4JdX0I7J3G7/BJ4KsQ35tAzyrTehfwQhjve8CWCf33JVoPlwF3AuNI+K1WGU8h0br7Rfi7NXRrBawGKsJyXxGW5RCq+X2lsLyHEK2DjxL9zn8SEzVvB/4X1onKeHaj+t/QP4EXw3LqX836cSnRejUPOCZhPGNJY9sWuv8GmBKW9TtAr4ThLwIWhvk1G9inxmUf5wa+vv/CzO4fZsR2QH5YGFuQYoIJ339IHsn6Vy0TFvoORBudXsDX/LiR6kp6CaYEOCjEfi0wvoaVZxjwLNAmjOdjYFBCzGXAWUAzog3TQ2Gl3AkoIkpac4HjwviuIjoFV1n/4UQ/ujyiRLsS2LSGefIQP/4Q+oVx/YLoR38H8GaV+fc8sAGwOdGP+oBq5tEQoiRxGFAA/DnEXRBimwz8DWgO/Az4DNg/YdhHE+pqEeZxhzD810Q/oDah32qgfQr1ngOMB7qE6bsXeLzKMh8a6twRKCVhg1/DfNs7LLcrQ3wHAasIG3+ijfNYoHNYZruH8W8dls++YbgLiZJ784T1bDxRUulMlLzeB36esC5cHsp2JkqaB4X5sG/4Xpzi7/CkMD8rE8CUKtO6GOhLtF4+BgwP/ToQbcgql/N5YV5Ul2CuDNO0MVBMtHH8e8J8XJBkPUr6+0pheQ8hWgd/G8q2SBLP3qSxHaD639AyYI9QTxHJ14+bw/zdKyz3beq4bft5WBd2CfPkeKJ1pRDYhihRdkqYhi2TLYvKv1w4RQbwCNFGc19gJtEGJDZmNtbMpppZhZl9BDxOtODr4r9m9qKZlRNNx47JCoVD54HAJWa23MzmEe3ZHptQ7Aszu8PMysxsdeg2yswmm1kJMAooMbNhYXxPEK1wldP1pJl9EabrCaKjjb4pTscxwINm9r6ZlQKXALtVaQe7zsyWmtn/gDeA3jXUN9nMRprZWqIfVxHRXvbORBu+K81sjZl9RrRhH5iskjAfJgK/Ikq0HxIdye0R6vvEzBanUO+pREclC8L0DQEOq3Lq5AozW21mH4bxJF2WSawFrjSztWb2ItHe5jaS8og23ueY2UIzKzezd8L4jwReMLNXwzz6B1Fy2z2h3jvM7GszW0h0pPaemX2QsC5ULvs/Ai+G9bDCzF4FJhFtmGtlZg+GdbJyvuwoqV1CkVFmNsHMyogSTO/Q/SBgesJyvpXoSKg6x4T59I2ZLQKuYN31P5nqfl+prEfvmtkzYZ6srlJvJrcDz5rZ26GekmrK/NXMSs1sHNHR4BF1GA/AYOBeM3svrE8PE+0M7Up0dFUI9JBUYGbzzOzTmirLlfPwjxAdDnYj2suPlaRdgOuA7Yn2fgqJThPUReIPahVQJKlZ+DEmqtwD/zyh2+dEe5+V5iep/+uEz6uTfG9d+UXSccD5RHsuhH7rXEBRg05Ee8gAmNkKSYtDfPNC56rT2prq/TAtZlYRGnA7Ee2RdZK0NKFsPtEGtDrj+PFUwzii04N7Ef2wxoUyW9RS7xbAKEkVCf3LiY4QKqUzfYkWV1nelcN2IEqsyX7knUhYF8I8ms+660Oqy34L4HBJhyT0LyDaCahR2PG5mujot5joNBUh9mXhc3XzpRPrLmcL01CddaY5fO5US4hJf1/Uvrwh+e/pBxncDtQ4HmCJma1M+J7KdFdnC+B4SWcldGtOdNQyTtK5RDsJPSWNAc43sy+qqywnjmDM7HOiUygHAU8nKbKSdRv+O9ZS5TrlJVUt/2+iNobNzKwdcA+gNMNO17dEe7pbJHTbnHWP1qyulUvagmgP7kygvZltAEzjx+mqre4vEmOT1Iro1FNdjyY3S6grj+jU1BdEP8a5ZrZBwl8bM6vc204WZ2WC+VX4PI4owezFjwmmtnrnAwdW6V8Ujg7i8i3RKZ4tk/SrOr9FNM/qEs984JEq09bKzK5LYdijgQFEp6rb8ePOSSq/hy9Zdzkr8XsS60wz0fpfufFLd92vbXmnUmdN24Fkw1ZXX23j2TD8niolTne627b5wNVVprulmT0OYGb/NrM9+bGZ4fqaKsuJBBMMAvpVyfSVpgC/l9QyXLI3qEr/r4nOwVb6kCiD95ZURJTRE7UBvjOzEkl9iX5ksQqH+COAqyW1CQnhfKJGyExoRbRCLQKQdCLRnlmlr4EukppXM/zjwIlhnhUC1xCdkplXx3h2kvT7sLd5LtHRxnhgArBc0kWSWkjKl7S9pJ0T4uwaklKld4jOL/cFJpjZdKIf0C5ER76kUO89RPN+CwBJxZIG1HHaUmJmFcCDwM2SOoWYdgvzdwRwsKR9JBUQXSRQGqY1XY8Ch0jaP4yjKFz22wV+uPR7bDXDtgnjXUy0obsmjfG+QPQ7q1zOZ1PzBvJx4LIw7zsQtZ9Urv9fA+2rnJqrSW3LOxU1bQcWER3NJW5XavsN1eQKSc0l/ZKokb7ySGkK6W3bhgKnStpFkVaSDg7blG0k9QvrVwk/XjhRrZxJMGb2qZlNqqb3LURXF30NPEx0HjjREOBhSUslHWFmHxM1KL5G1A5R9dr104ErJS0nWslHZGYqanUW0R7LZyGmfxNtgNabmc0gatN5l2g+7UDUVlHpdaKrcL6S9G2S4V8D/kp0Vc6XRHvdSdtFUvQsUTvDEqLz7L8PbRTlRD+w3kRHrd8C9xPtPcOPP7zFkt4Psa0kOn033czWhP7vAp+b2TehTG313ka0t/pKWO7jiRJU3P4MTCVqR/qOaI8yz8xmE7Wd3BFiPQQ4JGH6UmZm84mOQi4l2jDOJ7oCrXL7sRnrrguJhhGdsllIdOXh+DTG+y3RqbXriBJU9xrGA9FFKZOIroybSrRMrwp1zSJKQJ+F33GNp5BSWN6pqHY7YGariE4dvh3i2ZVafkM1+Irod/AF0bbr1DC9kP62bRLR1Z53hjrnEF0YANEpvuuI5sVXRBdTXFJTYDKr81kT55xD0hSiy1UX13csrmHxBOOccy4WOXOKzDnnXHZ5gnHOORcLTzDOOedikSs3Wv6gQ4cO1rVr1/oOwznnGpXJkyd/a2bF6QyTcwmma9euTJpU3dXKzjnnkpH0ee2l1uWnyJxzzsXCE4xzzrlYeIJxzjkXC08wzjnnYuEJxjnnXCw8wTjnnIuFJxjnnHOx8ATjnHMuFjl3o6VzmRS9ZDF9/hRzlws8wTi3HmpKFJI8kbic5qfInHPOxcKPYBoBPw3jnGuMPME0An4axjnXGPkpMuecc7HwBOOccy4WnmCcc87FwhOMc865WHiCcc45FwtPMM4552LhCcY551wsPME455yLhScY55xzsfAE45xzLhaeYJxzzsXCE4xzzrlYeIJxzjkXi6wlGEkHSJotaY6ki5P0L5T0ROj/nqSuCf0uCd1nS9o/ofs8SVMlTZE0KUuT4pxzLgVZeVy/pHzgLmBfYAEwUdJoM5uRUGwQsMTMtpI0ELgeOFJSD2Ag0BPoBLwmaWszKw/D/drMvs3GdDjnnEtdto5g+gJzzOwzM1sDDAcGVCkzAHg4fB4J7KPoTVsDgOFmVmpmc4E5oT7nnHMNWLYSTGdgfsL3BaFb0jJmVgYsA9rXMqwBr0iaLGlwdSOXNFjSJEmTFi1atF4T4pxzLjWNvZF/TzP7BXAgcIakXyUrZGb3mVkfM+tTXFyc3Qidcy5HZSvBLAQ2S/jeJXRLWkZSM6AdsLimYc2s8v83wCj81JlzzjUY2UowE4HukrpJak7UaD+6SpnRwPHh82HA6xa9bH40MDBcZdYN6A5MkNRKUhsASa2A/YBpWZgW55xzKcjKVWRmVibpTGAMkA88aGbTJV0JTDKz0cADwCOS5gDfESUhQrkRwAygDDjDzMolbQKMiq4DoBnwbzN7ORvT45xzrnaKDhJyR58+fWzSpKZzy4wkcm0ZNha+bFxTImmymfVJZ5jG3sjvnHOugfIE45xzLhaeYJxzzsXCE4xzzrlYeIJxzjkXC08wzjnnYuEJxjnnXCw8wTjnnIuFJxjnnHOx8ATjnHMuFp5gnHPOxcITjHPOuVh4gnHOORcLTzDOOedi4QnGOedcLDzBOOeci4UnGOecc7HwBOOccy4WnmCcc87FwhOMc865WHiCcc45FwtPMM4552LhCcY551wsPME455yLhScY55xzsfAE45xzLhZZTTCSDpA0W9IcSRcn6V8o6YnQ/z1JXRP6XRK6z5a0f5Xh8iV9IOn5LEyGc865FGQtwUjKB+4CDgR6AEdJ6lGl2CBgiZltBdwCXB+G7QEMBHoCBwB3h/oqnQPMjHcKnHPOpSObRzB9gTlm9pmZrQGGAwOqlBkAPBw+jwT2kaTQfbiZlZrZXGBOqA9JXYCDgfuzMA3OOedSlM0E0xmYn/B9QeiWtIyZlQHLgPa1DHsrcCFQUd2IJQ2WNEnSpEWLFq3HJDjnnEtVo27kl/Qb4Bszm1xTOTO7z8z6mFmf4uLiLEXnnHO5LZsJZiGwWcL3LqFb0jKSmgHtgMU1DLsHcKikeUSn3PpJejSO4J1zzqUnmwlmItBdUjdJzYka7UdXKTMaOD58Pgx43cwsdB8YrjLrBnQHJpjZJWbWxcy6hvpeN7M/ZmNinHPO1axZtkZkZmWSzgTGAPnAg2Y2XdKVwCQzGw08ADwiaQ7wHVHSIJQbAcwAyoAzzKw8W7E755xLn6IDhNzRp08fmzRpUn2HkTGSyLVl2Fj4snFNiaTJZtYnnWGydgTjnGv4orsC0ueJ1CWTdoKR1Aoo8VNUzjU91SUKPxpzdVFrI7+kPElHS3pB0jfALOBLSTMk3Shpq/jDdM4519ikchXZG8CWwCVARzPbzMw2BvYExgPXS/Irt5xzzq0jlVNk/c1sbdWOZvYd8BTwlKSCjEfmnHOuUav1CKYyuUjqLulBSXdVV8Y555yrlM6Nlo8ATwK/BJC0vaRhsUTlnHOu0UsnweSZ2UtAOYCZTQO2jyUq55xzjV46CeaL8JgWAwiP0W8RS1TOOecavXTugzmX6J0rHSWdSPTir2lxBOWcc67xSznBmNk8SQcAvwN6AeOAB+MKzDnnXOOW7p38MrMniRr7nXPOuWqlnGAkDQV+I6kM+AL4CPjIzO6IKzjnnHONVzpHML8CuphZuaTOwI5Ep8qcc865n0gnwbwHtCd6RfFCojdKvhhLVM455xq9dC5TvhcYJ+nPkn4pqV1cQTnnnGv80kkwjwLDiI56TgfekfRpLFE555xr9NI5RbbAzK5N7CCpMMPxOOecayLSOYKZIumcxA5mVprheJxzzjUR6RzBbAL0l3QR8D7wITAl3BfjnHPOrSOdO/mPgB9Oi/UEdgB2wW+6dM45l0Q6N1puBJwHbAzMAIaZ2cNxBeacc65xS6cNZjiwHHgOaAn8V1LfWKJyzjnX6KXTBlNsZjeEz89LegL4N7Br5sNyzjnX2KVzBPOdpB0qv5jZZ0RHMs4559xPpHMEczrwlKS3gKlEDf1+o6Vzzrmk0jmCaQPsDrxB1NA/BTgqhpicc841AekkmGHAGjMbYWZDgFFA/3RGJukASbMlzZF0cZL+hZKeCP3fk9Q1od8loftsSfuHbkWSJkj6UNJ0SVekE49zzrn4pJNgSsyspPKLmX0LXJnqwJLygbuAA4EewFGSelQpNghYYmZbAbcA14dhewADiU7LHQDcHeorBfqZ2Y5Ab+AASX7RgXPONQDpJJjPJB1YpVvzNIbvC8wxs8/MbA3RZc8DqpQZAFTeWzMS2EeSQvfhZlZqZnOBOUBfi6wI5QvCn6URk3POuZik08h/FvCSpGOB8aTfyN8ZmJ/wfQHRkwCSljGzMknLiN5B0zmMM3HYzvDDkdFkYCvgLjN7r+qIJQ0GBgNsvvnmaYTsnHOurlI+gjGzL4GdgKeAYqJXJh8dU1wpM7NyM+sNdAH6Sto+SZn7zKyPmfUpLi7OeozOOZeL0jmCwczKiRLMU3UY10Jgs4TvXUK3ZGUWSGoGtAMWpzKsmS2V9AZRG820OsTnnHMug9J5Flk/4BhgKdEG/CNgWhqP7J8IdJfUjSg5DOSnR0CjgeOBd4HDgNfNzCSNBv4t6WagE9AdmCCpGFgbkksLYF/ChQHOZUpFRQVr166t07Clpam/0SIvL4+CgoI6jce5hiidI5gHgXOJGtJ7Ab8laofZKpWBQ5vKmcAYIB940MymS7oSmGRmo4EHgEckzQG+I0pChHIjiB6yWQacYWblkjYFHg7tMHnACDN7Po1pcq5Wg08/k38NvQ/lpXNNDKA8WrVuk3LxVq3b8MWC/9GqVas0I3SuYZJZahddSRpnZnvFHE/s+vTpY5MmTarvMDJGEqkuQ1c3+x50CB8W9qLlNrvHOp5v7j6G+XM/pX379rGOpy58PXOSJptZn3SGSWeX7E1J54XLhp1zzrkapXOKrAfRS8YukjSZ6FEx/kZL55xzSdV6BCNpSPh4C9Flyt2AvwGf8NP7WJxzzjkgtSOYMeH/2USN+s2JGts/An5yU6NzzjkHKSQYM3s3/D8CogdSEiWaHYge/+KnyJxzzv1EWjdaAoT7Xt4Pf84551xSaV7Y75xzzqXGE4xzzrlYpJxgJHWs8n3T0B7jnHPO/UQ6bTAPAAcnfH8E2FLSU2b258yGlXsuuexyht5/f9rDqaCIDh07p1y+Y8eOTH7vHQoLfd/ANW11vSfcn1iQOSknGDM7uMr3/uGu/qpvpXR1MH7CBOznh9GiW++0huuU5ng+eehMSkpKPMG4Jq+mROGPvsmOtK8iS2TREpqeoVhyXn7LdjRru3Gs45C82c05lx1pJxhJRwOHAuWAgOfM7PFMB+acc65xq8sRzF5mNrDyi6S7AE8wzjnn1lGXBFMo6WBgPtGbJVtkNiTX1Hnjq3O5IaUT8pLyJF0avp4ObAgcFP6fGVNsrokys6R/NfXz5OJc45PSEYyZVUj6DXCNma0CHo03LOecc41dOpcUfSTpcvllSM4551KQThvMRsBewGmS3iN6XP9H/sIx55xzyaSTYC40s3n+uH7nnHOpSCfBPA38IvFx/ZJ2jScs55xzjV0qr0w+QtJ1QBtJ21Vpg7kvvtCcc841ZqkcwbwNFAEnAzcD20haCnwBrI4vNOecc41ZKq9MXggMkzTXzN4CkNQe6ArMijc855xzjVU6lxzfVvnBzBab2WSihn7nnHPuJ7wNxjnnXCxSOYJ5G5hB9FiYm4E5kt6X9DxptMFIOkDSbElzJF2cpH+hpCdC//ckdU3od0noPlvS/qHbZpLekDRD0nRJ56Qai3POufil0wbzqZm9Dem3wUjKB+4C9gUWABMljTazGQnFBgFLzGwrSQOB64EjJfUABhLde9MJeE3S1kAZcIGZvS+pDTBZ0qtV6nTOOVdP0mmDmSXpNEknAlsCM8xsZYrD9gXmmNlnZrYGGA4MqFJmAPBw+DwS2Ce8MXMAMNzMSs1sLjAH6GtmX5rZ+wBmthyYCaT+7mDnnHOxSifBjAKKgWuAG4FlklK9iqwz0eP9Ky3gp8nghzJmVgYsA9qnMmw4nfZz4L1kI5c0WNIkSZMWLVqUYsjOOefWRzoJpo2ZXQl8bWZ7AUcBI+IJK3WSWgNPAeea2ffJypjZfWbWx8z6FBcXZzdA55zLUekkmJLwv1RSCzN7CtgvxWEXApslfO8SuiUtI6kZ0A5YXNOwkgqIkstjZvZ06pPinHMubukkmH9I2gh4AnhQ0lnABikOOxHoLqmbpOZEjfajq5QZDRwfPh8GvG7RW6ZGAwPDVWbdgO7AhNA+8wAw08xuTmM6nHPOZUHKD7sMRywAN0s6lugmy9+nOGyZpDOBMUA+8KCZTZd0JTDJzEYTJYtHJM0BviNKQoRyI4gulS4DzjCzckl7AscCUyVNCaO61MxeTHWanHPOxafWBCNJVuV9tWb2SG1lqgob/herdPtbwucS4PBqhr0auLpKt/8CdXu5u3POudilcorsDUlnSdo8saOk5pL6SXqYH09tOeecc0Bqp8gOAE4CHpf0M2AJ0IIoOb0C3GpmH8QXonPOucYolQRzC9HrkS8iunO/AFhtZktjjMs551wjl0qC+YCoQb/ycS0rgY8kTQU+MrPhMcbnnHOukUrlWWTrPDFZUheihNMLOJjosS/OOefcOlK+TLmSmS0gelzLS5kPxznnXFORdoJxLte0adWK8uVJn0KUMRVrSihbs4bmzZvHOh6AtWvXUstdBUmtWbMmrfIFBQVE90O7XJXOnfzO5aQTjj0aPh4b6zhWzXqTPX+1F23atIl1PAsXLqRtuw1o2ap1Wn8oL63yRS1acNU118Y6La7h8yMY52px0EEHocGnsuabuTTfuFss47CZr3H+HdfHUneiJUuW0GKjjmzyx9tjHc/y959n3uf/i3UcruHzIxjnatGsWTNOG3wKa6a9Ekv9a77+lLySpRx44IGx1O9cffEE41wK/jT4FFbNHEfFmpTfEp6ytdNe4YxT/0R+fn7G63auPnmCcS4FXbp0Yfc992TlzDczWm9F6SpWzHyLwaecnNF6nWsIPME4l6ILzj4Tm/lqRutcNXMce+29F506dcpovc41BJ5gnEvRfvvtR8GaFZR+NScj9ZkZFTNe5fyzz8xIfc41NJ5gnEtRfn4+Z55+KmXTx2SkvjVffkyhldK/f/+M1OdcQ+MJpoEoKCigYs2qWMdhZWspL1vrjcnr4ZSTB7Fi1ttUlK7/slo7/RXOPv1U8vL8Z+iaJl+zG4gTjhmIzXo91nGsnPEGu+3xS1q3bh3reJqyjh070q9fP1bOGLte9VSUrGDl7HcYNOikzATmXAPkCaaBOPzww8lb8TWlX34SS/1mFZRNGc3lf7k4lvpzyXlnnUHFjFfq9LiVSiunv8G+++3PxhtvnMHInGtYPME0EAUFBVx4wfmUTRkdS/2rP5vMxhu0pl+/frHUn0v69etHy7xy1nwxu07DmxkVM1/lvLNOz3BkzjUsnmAakD8NPoWSz6dQtuzrjNddPuVZLv/Lxf7wwQzIy8vj7NNPZW0dG/tLF86gdYHYe++9MxuYcw2MJ5gGpG3btpxy8iBKPnguo/WWfvkxeSu+4YgjjshovbnspJNOZOXH4ykvWZH2sOXTX+WcM0/zZO+aPE8wDcwF553Lqulv1GnDVZ21U0Zz4QXnU1BQkLE6c11xcTH7H3Agq6b9J63hyld/z4pP3uPEE06IJzDnGhB/mnID07lzZ35zyG94/cOXab3LYetdX9myrymZN4U/DX42A9G5RIccuB/PzD8feCL1gYqAPwvat48rLOcaDE8wDdBlF1/Ii3v3p9VOA1Cz9TvqKHl/NKecPIi2bdtmKDpX6YFhj3FR21No3fPXKQ9TsbaUb4cOYvrxc+nWLZ5H/zvXUPgpsgaoV69e/HzHXut9r0X56uWsnPEGfz7/vMwE5n4wZ84cpkyZQqtt9khruLyCQlr2+DV333NvTJE513D4EUwDdflfLuaw4wdjO/Svc2Pw6o9e5pBDDqFz584Zjq5mZsYLL7zAihXptyMNHz48rfL77LMPxcXFaY9nfd19z7207NEPNUv/FceF2+/H0Pv/xtV/vzIrr0h2rr5kLcFIOgC4DcgH7jez66r0LwSGATsBi4EjzWxe6HcJMAgoB842szGh+4PAb4BvzGz7LE1KVuyzzz5s3K4Vyz6bRIstd057eCtbS8mHL3LZba/FEF3Npk2bxmFHHk27rfukNVzRz3binOvvS7n8qm/mc/wfxnPn7bemGeH6KS0t5YEH/0Xrw66p0/AFHTaj2UadefbZZzn88MMzHF1uef/993n55ZfrNOw116S+/IqLizn55JP9yr80ZSXBSMoH7gL2BRYAEyWNNrMZCcUGAUvMbCtJA4HrgSMl9QAGAj2BTsBrkrY2s3LgIeBOosTUpEji8r9czDlX3AR1SDArZ7zBz3v3olevXjFEV7Py8nJadehIi/0vSGu4FmmOZ+2kZykrK09zqPU3atQoCoq3oGCjuh8Zart9uen2uzzBrKfLrriKN2d9RfMOm6U1XNtd/sA/Xvgw5fLLJ43m0EMPZZNNNkk3xJyWrSOYvsAcM/sMQNJwYACQmGAGAEPC55HAnYp2FwYAw82sFJgraU6o710ze1NS1+xMQvYdccQRnPvnCyn9ag6FHbdKebgfHgszbGiM0eWuf9x2J9puv/Wqo+XWuzN16P188skndO/ePUOR5SCDFtv+ilbb/TLW0ZROz/6ZgKYgW438nYH5Cd8XhG5Jy5hZGbAMaJ/isE1SQUEB5599FmVTX0pruNL/TWPD1kXss88+MUWWu2bNmsXMmbNo2X2X9apHzQpo0XMf7vznPRmKzLmGJyca+SUNBgYDbL755vUcTXo+mj6DRYe8C7yb+kDbRH/LV6ygTZs2cYWWk+68+x5a9OyP8tf/ptXCHfbjXw9dwvXXXE1RUVEGonOuYclWglkIJJ4k7RK6JSuzQFIzoB1RY38qw9bIzO4D7gPo06dP3R+Bm2VffPEFzz47ms1Oupf8FuklihUv3sil7e7n/PP8EuVMWb16NQ8PG0bbgTdmpL6CDTvRvLgbTz/9NEcffXRG6nSuIcnWKbKJQHdJ3SQ1J2q0r/rY4NHA8eHzYcDrFj0PfTQwUFKhpG5Ad2BCluKuVzffehste+yddnIBKOh9KNfdeBNlZWUxRJabRo4cSeGm3SnYoGPG6lSPffnHbXdmrD7nGpKsJJjQpnImMAaYCYwws+mSrpR0aCj2ANA+NOKfD1wchp0OjCC6IOBl4IxwBRmSHic6d7SNpAWSBmVjerJh+fLl3HvfUAp/fkidhi/stA3lLTswcuTIDEeWu/5x253kbbdvRutsudUufPzxJ8yYMaP2ws41Mlm7k9/MXjSzrc1sSzO7OnT7m5mNDp9LzOxwM9vKzPpWXnEW+l0dhtvGzF5K6H6UmW1qZgVm1sXMHsjW9MRt6P33U7h5r/XaW87f8VCGXHXter0Yy0WmTp3KZ3Pn0WKrvhmtV/nNKOy5D7ff9c+M1utcQ+CPimmAysrKuO7Gmyj4+YD1qqfFVjvz1eJljB07NjOB5bDb77qb5j36o7z8jNfdYof9efTRR1m1alXG63auPnmCaYBGjhxJeasOFG669XrVI+XRrPchXHHNdbUXdtVauXIljz8+nBY7ZPb0WKVm7TamqPO2jBgxIpb6nasvnmAaGDNjyFXXkr/jobUXTkGrnv2YNHEy06dPz0h9uWj48OEUdelBs7bxPfMsL9zZH7fWrVtT8v13GX3fUFJLFrBBO3+Cd67zBNPAjB07lq+++75Ozx9LRs2aU7jjgVx13Q0ZqS8XjXz2efjZbrGOo8WWffh45nSWLl0a63i6du3KsX88htX/uTu2trnV86bA/yZxyUUXxlK/azw8wTQwV1xzHc16H4KUuUXTsvdBPPvMM3z55ZcZqzOXlJWVkdesMNZxKC+f/ILmlJfH/2y1227+B+3WLEr7bZypKF/9PStfvZ3hjw6jQ4cOGa8/l0mq0199yok7+RuL6dOnM3HiZDqcfGpG681v0ZaWPfbi5ltv48brvT0m1xUVFfHsUyPYdY9f0bxLDwo27JSRes2MVa/dxYnHHsO++8bTXtWUVVRUMHToUJYtW5a0//XXX5+0+0UXXVRtP4Abbvjp2YsBAwawzTbb1C3QNHiCaUCuvv5GWvQ+qE7vGKlNYe9DuefeC7n8r5fRunXrjNfvGpftt9+eq64cwhU330LbI65F+eu/KVj10Rg68D3/uP7a9Q8wB3366aecemr1O5cb7vaHpN3b9BnAdc/89N7zkgWzWD0/edvr7E8+5YGh8b/0zhNMA/Hll1/yzKhnaH9SPPdDFGy4KYWb78DQ++/nvHPPjWUcrnE55+yzeOa5F5j27uO03vPY9apr7eL5rH73MZ55920KC+M9ndhUtWvXrsb+3+33ah1qbYuu+P4nXVu2TPflGHXjCaaBuOW222nZYy/yW8R35U1B7wFcd+NNnH3WWeTnZ/5+Dte4SOKJx4axbc8dKNm8N0Wb71Cneqx8LavG3MIN11zNdtttl+Eoc0dhYSG7/XJvllRzimy7UcmHmz17do2nu7at8kooM2PnnTNzEVFtPME0EBPfn0Jel/TeAJmuwk7b8NWSJaxYsaLWvSWXGzbZZBMeG/YQRx53EgV/vI38ovRPn658+zF27tmd007LbNthrmnXrh3vvPlG2sNJYuaH78cQ0frzq8galCxc8eGvfHVVHHTQQRx1+B9YVYdLl1fPm4LNeYvHhv2r3q9Ycg2PJxjnHDdcezXLPx5P2ZIv0hqu7P1R/PXSiykuju8mVBep6TJkv0zZOddgXfrXy2nTfRcWdvq/9AY8Glh6OYsXH0v79u1jic1FGuNDaz3BOJfjxowZwyP/foIN/ngrXUvSf/fQinEPsvPkk3hx9DP1vsfsGhY/ReZcDlu0aBFHH3s8rfY7p04vtgNotcexvDtlBkOH3p/h6Fxj5wnGuRxlZhz1x+NR919RtEWv2geohpoV0HL/8zn/wouYPXt2BiN0jZ0nGOdy1J133c2kmZ/Ravej17uu5h02p8WuR/G7wweyZs2aDETnmgJPMM7loBkzZnDJZX+l1f7nofyCjNTZcscD+bqsBRf/5bKM1OcaP08wzuWY0tJSfnvYkbTY/Y8UtO+SsXol0bL/mdz3wEO88Ub6Nwy6pscTjHM55oL/u4jFeRvScof9Ml53fst2tNr3LA4/6hgWL16c8fpd4+IJxrkcMm/ePB586GFa7nNabJcUt+j2C6zLL7j2en/JXa7zBONcDlmxYgVFbTeK9aGqAGy0GUuWJn9oo8sdnmCcc87FwhOMyzhJrF21EitfG+t4KlZ/T16e3zme0wRWFu9l0VZRTkV5WazjaKo8wbiM69mzJ7vvshMrnruWirUlsYxj5dRXyfv4Dc447U+x1O8ah5OPP5aStx+mdOHMWOqvWFvK8uevY9dddvUHetaBJxiXcc2aNeO5UU/x6x1/xopn/k5F6aqM1r9y8mjyPhjJ+LffomfPnhmt2zUuv//973jy8UdZ/tw1rP5sckbrrihdyYpnrmTvnpsx5sXnyMvzzWW6fI65WBQUFDDi8cf4bb++LH/6b5Sv/ulrW9NlZqwYP5yiOa8xcfw7bL311hmI1DV2Bx54IGNeeI6S125n1aw3M1Jn+cqlfP/UX/lD/90Y+cTjNG/ePCP15pqsJhhJB0iaLWmOpIuT9C+U9ETo/56krgn9LgndZ0vaP9U6Xf3Jy8vjgfvu5YTDf8P3I/9C2Yrv6lyXmbHqrYfY8Kv3mfju22y++eYZjNQ1dnvssQdvjX2d8nceZuWUl9arrrJlX7PsyUs4/bgjuO+eu/3IZT1kbc5JygfuAg4EegBHSepRpdggYImZbQXcAlwfhu0BDAR6AgcAd0vKT7FOV48kcfONN3DBqSfx/ZOXUrbsm7TrsIpyVr7+Tzqunsd777xFx44dY4jUNXa9evViwrtvUzDjeVa+N6JO709Z8+3/WDriUoZcdAHXXPV3f/3Aesrm+2D6AnPM7DMAScOBAcCMhDIDgCHh80jgTkVLeAAw3MxKgbmS5oT6SKHORkESpVNfxhZMiXU8ZWtKY60fYPXq1Vxy2V9ZvnzFOt2Xf7OA5fecxEa9+5PXvEXK9X074TkA9jnqaP7v4kvX6Td40EnsskvfZINljCRKpr2CfTE11vGsWZ3ZtqrqrF76Lflv3BvrOEq/nou67hHrOABeeeVVnhj51DrdftG7Ny88NwzGDqND30NSrmvt8sUsm/kOrdu0YfrMmQwafOoP/TbZeGP+fsXl5OfnZyz2XJDNBNMZmJ/wfQGwS3VlzKxM0jKgfeg+vsqwncPn2upE0mBgMNBgT63cftMNjBs3Lmm/M888s0513nnnnT/p1vHC42nXrl2d6kvV7NmzGfrgMIp2/sMP3RL3JhcPmJBehQdGNwW2e6mEZu1+vEGw5POPKGj2SOwJ5tYbr6v22VqZXDab/N9xsb8VskePHtx79x0sX748af9MTs/BBx9cp7rSceud/+TNecsp7LjVD93WfL0aALu8LZD8N1W9tuRfV86ouaxz9LLikZs456wz2GSTTTIQde7IiTdamtl9wH0Affr0aZDvHd1+++3Zfvvtk/Y744wzshzN+itquwFtfvEbAKy8jBWv3kGfXffgP2NehLbp30V+0823UPHpDRT9/goKNor2LayiPKMxV6dHjx706JH8zGtjWzZ5eXkcd9xx1fZvbNMD0KLbTrTa7pcAlMyfRumEJxkx4kk4/LC061q0aBHbjtqfr5fOp9Wv/4TyoiOWkgkjMhpzrshm69VCYLOE711Ct6RlJDUD2gGLaxg2lTpdPbKyNax48QZ6dchn3H9eoW0dkgvABeefx83X/Z1lIy9jzTefZThK1xSsmvMeK1+8gVFPPs7hdUguAMXFxYz/7zi2KlzBipdvjv1m4aYumwlmItBdUjdJzYka7UdXKTMaOD58Pgx43aJzK6OBgeEqs25Ad2BCinW6elKxpoTlo69m9607MubF52jZsuV61XfKySdz/9138P3TQyhdOCtDUbqmYOX011nzxj/5z5iX2HfffderrjZt2jD2tTHsvFlblo++moo1qzMUZe7J2imy0KZyJjAGyAceNLPpkq4EJpnZaOAB4JHQiP8dUcIglBtB1HhfBpxhZuUAyerM1jS56pWVrGb5M0M4cI9f8MhDD2ascXTgwCNp06Y1Rx5zHBWbbANb75SRel3jtXLqK7Ra/TXv/vdNtttuu4zUWVRUxPPPPs3xJ53MC09fTvlaf0tnXWT1Am8ze9HMtjazLc3s6tDtbyG5YGYlZna4mW1lZn0rrw4L/a4Ow21jZi/VVKerf99/s4AjD9ibRx/+V8avvDn44IN54dmnqZj/UUbrdY1P8+YFdLDvmTT+nYwll0rNmjXj0Yf/xXG/O4A1q1f6/TB1oLpcK96Y9enTxyZNmlTfYTRpFRUVTJw4kb59+8Z6H8GsWbPYeOON2WijjWIbh2vYSkpKWL16NRtuuGGs4/nqq69y/v4rSZPNrE86w+TEVWQuu/Ly8thll59cLZ5x2267bezjcA1bUVERRUVFsY8n15NLXfkxn3POuVh4gnHOOReLnGuDkbQcmF3fcWRQB+Db+g4iQ5rStEDTmp6mNC3Q9KYnG7Yws7ReipOLbTCz022oasgkTWoq09OUpgWa1vQ0pWmBpjc9DZWfInPOORcLTzDOOedikYsJ5r76DiDDmtL0NKVpgaY1PU1pWqDpTU+DlHON/M4557IjF49gnHPOZYEnGOecc7HwBONcBkk6QdJPX+/YxEjqJGlkmsM8JKluL2r5sY531md4l12eYFxWKeLrXSMSXv63zncz+8LM1itZ1IWZ7Z7tcbq6a3I/dEnPSJosabqkwaHbIEkfS5ogaWjlHqakYklPSZoY/vao3+h/SlIrSS9I+lDSNElHStpJ0rgwnWMkbSqpnaTZkrYJwz0u6ZT6jh9AUtcQ2zBgGvCApElhGV2RUG5nSe+EaZ0gqY2kfEk3huXzkaQ/1d+UgKQ/htimSLo3xHdi5foF7JFQdp09dkkrEj5fJGlqmNbrYoq1q6RZIY6PJT0mqb+ktyV9Iqlv+HtX0gdh3leuPydIGi3pdeA/Sb53lTQtlE26jMLOxJ1h2b8GbJyBaVoR/u8taaykkWEaH1N4dHc161GRpH+Fef6BpF8nTOczkl6VNE/SmZLOD2XGS9oolNtS0svhN/eWJH/SairMrEn9ARuF/y2INmadgXnARkAB8BZwZyjzb2DP8HlzYGZ9x59kev4ADE343g54BygO348ketEawL7Au0Qvanu5vmNPiLkrUAHsWmUZ5QNjgV5Ac+AzYOfQry3RkyYGA5eFboXAJKBbPU3HdsBzQEH4fjfRG1j/BxSHaXg7Yf16CDgsYfgV4f+BYRm2TJwfMc33MmAHop3JycCDgIABwDOV8zmU7w88FT6fACxIWFZVv3cFpoXPSZcR8Hvg1bCcOwFLE+dHHaepch7uDSwjek16Xljv96xhPbog4XeybVhmRWG65gBtwjJcBpwayt0CnBs+/wfoHj7vQvS23Xr/bTX0v6b4qJizJf0ufN4MOBYYZ2bfAUh6Etg69O8P9NCP7yxpK6m1ma2g4ZgK3CTpeuB5YAmwPfBqiDsf+BLAzF6VdDhwF7Bj/YRbrc/NbHz4fISio8tmwKZAD8CAL81sIoCZfQ8gaT+gV8KRQDuiV2bPzWbwwT7ATsDEMO9bALsDY81sEYCkJ/hx/apOf+BfZrYKoHLdjMlcM5saYpsO/MfMTNJUoiTRDnhYUneiZVCQMOyrVWKr+r1SdcvoV8DjFr199otw9JNJE8xsAYCkKUTTs4zk69GewB2h2yxJn/PjcnrDzJYDyyUtI9qJgOi310tSa6Ll/GTCtqIww9PSJDWpBCNpb6If725mtkrSWGAW0Z5nMnlEe9UlWQmwDszsY0m/AA4CrgJeB6ab2W5Vyypq29gOWAVsSLTH2VCsBJDUDfgz0R7mEkkPEe1JVkfAWWY2Jv4QayXgYTO75IcO0m+J9tSTKSOchg7LpnncASZRmvC5IuF7BdHv/+9EG9jfSepKdERZaWWVuqp+r5R0GUk6qI4xpypx2sqp+/astnmUByw1s951rD9nNbU2mHbAkpBctgV2BVoBe0naUFFj5R8Syr8CnFX5RVLvbAabCkmdgFVm9ihwI9HhebGk3UL/Akk9Q/HzgJnA0cC/JBUkq7OetSXaUC2TtAnR6SKInnC9qaSdAcJ582bAGOC0ymmRtLWkVvUQN0SnSQ6TtHGIZSPgA6L1q32I8fCE8vOIjngADuXHo4NXgRMltUyop760AxaGzyfUsY7qltGbwJGhjWZT4NfrG2wKqluP3gKOqYyP6JR4Sk9VD0dBc8PZgcq2pYZ2hqBBalJHMMDLwKmSZhKtPOOJfjzXABOA74iOaJaF8mcDd0n6iGhevAmcmu2ga7EDcKOkCmAtcBrRnvHtktoRxX2rpDLgZKCvmS2X9CZwGXB5PcWdlJl9KOkDouUwn6jNAjNbI+lI4A5JLYDVREej9xOd+ng/NOIuAn5bD6FjZjMkXQa8Eo5I1gJnAEOI2gCWAlMSBhkKPCvpQ6J1c2Wo5+WwMzNJ0hrgReDS7EzFT9xAdIrsMuCFOtZR3TIaBfQDZhC1eby7vsHWpob16G7gn+HUYBlwgpmVKvVXeh8Thr+MaEdhOPBhxiegicmJR8VUtquEPZlRRI19o+o7Lueca8qa2imy6gwJjYDTiBqHn6nXaJxzLgfkxBGMc8657MuVIxjnnHNZ5gnGOedcLDzBOOeci4UnGOfqmaRLEz5vIOn0+ozHuUzxRn7n6pmkFWbWOnzuCjxvZtsnKdfMzMqyHZ9zddXUbrR0rkGT9AzRM/KKgNuAnwEtwmX004meLbdl+P4q0c2Pfyd6Bt221P6cM+caDD+CcS6LJG1kZt+Fu8wnAnsRPQg06RFMeL7eC8D2ZlYfD/h0rs78CMa57Kr6tO/uKQwzwZOLa4w8wTiXJdU87bumJ0lXqu4pxs41aH4VmXPZk+xp3wBrE558vZzo5VfONXqeYJzLnpeBZuFp39cRPe0b4D7gI0mPmdli4G1Fr8e+sb4CdS4TvJHfOedcLPwIxjnnXCw8wTjnnIuFJxjnnHOx8ATjnHMuFp5gnHPOxcITjHPOuVh4gnHOOReL/wdtaqR6rFGaKgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "trace_length=20000\n",
    "attrs=['age','sex','race','educ','married','income']\n",
    "\n",
    "trace_attr = lambda attr : np.concatenate(trace.posterior[attr],axis=0)\n",
    "\n",
    "y=[[pv.mi_sklearn([trace_attr(attr)[:trace_length,0], trace_attr('output')[:trace_length]],\n",
    "                  n_neigh=40,input_inferencedata=False)[0]\n",
    "   for attr in attrs] for i in range(0,100)]\n",
    "\n",
    "plt.boxplot(np.array(y),attrs,\n",
    "            showmeans=False, showfliers=False, patch_artist=True, vert=True)\n",
    "plt.xticks(range(0,len(attrs)), attrs)\n",
    "plt.xlabel('attr')\n",
    "plt.ylabel('$I(attr_0;income)$')\n",
    "plt.title(\"Mutual information between income, and other attributes\")\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the figure above, we observe that the mutual information between the output and any of the attributes is very low $I(\\mathit{attr};\\mathit{output})\\leq 0.005$ for all attributes ($\\mathit{attr}$)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}