analysis/unknown_dataset_size.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Working with Unknown Dataset Sizes\n",
    "\n",
    "This notebook will demonstrate the features built into SmartNoise to handle unknown or private dataset sizes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up libraries and load exemplar dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-15T16:41:40.471981Z",
     "iopub.status.busy": "2021-04-15T16:41:40.458021Z",
     "iopub.status.idle": "2021-04-15T16:41:40.634785Z",
     "shell.execute_reply": "2021-04-15T16:41:40.635247Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dimension of dataset: (1000,)\n",
      "Names of variables:   ('age', 'sex', 'educ', 'race', 'income', 'married')\n"
     ]
    }
   ],
   "source": [
    "# load libraries\n",
    "import os\n",
    "import opendp.smartnoise.core as sn\n",
    "import numpy as np\n",
    "import math\n",
    "import statistics\n",
    "\n",
    "# establish data information\n",
    "data_path = os.path.join('.', 'data', 'PUMS_california_demographics_1000', 'data.csv')\n",
    "var_names = [\"age\", \"sex\", \"educ\", \"race\", \"income\", \"married\"]\n",
    "\n",
    "data = np.genfromtxt(data_path, delimiter=',', names=True)\n",
    "age = list(data[:]['age'])\n",
    "\n",
    "print(\"Dimension of dataset: \" + str(data.shape))\n",
    "print(\"Names of variables:   \" + str(data.dtype.names))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see above this dataset has 1000 observations (rows).  Oftentimes the number of observations is public information.  For example, a researcher might run a random poll of 1000 respondents and publically announce the sample size.\n",
    "\n",
    "However, there are cases where simply the number of observations itself can leak private information.  For example, if a dataset contained all the individuals with a rare disease in a community, then knowing the size of the dataset would reveal how many people in the community had that condition.  In general, a dataset maybe composed of some defined subset of a population, and the dataset size is then equivalent to a count query on that subset, and so we should protect it like any other query we want to provide privacy guarantees for.\n",
    "\n",
    "SmartNoise assumes the sample size is private information.  If it happens that you actually know the dataset size, then that information will be valuable if you add it into your analysis graph.  However, SmartNoise will not assume you truthfully or correctly know the size of the dataset.  (Moreover, it can not directly send you an error message if you get this value incorrect, or this would permit an attack whereby an analyst keeps guessing different dataset sizes until the error message goes away, thereby leaking the exact dataset size.)\n",
    "\n",
    "If we know the dataset size, we can incorporate it into the analysis as below, where we provide `data_n` as an argument to the release of a mean on age:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-15T16:41:40.645482Z",
     "iopub.status.busy": "2021-04-15T16:41:40.644181Z",
     "iopub.status.idle": "2021-04-15T16:41:40.686094Z",
     "shell.execute_reply": "2021-04-15T16:41:40.685529Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DP mean of age: 45.2\n",
      "Privacy usage: approximate {\n",
      "  epsilon: 0.1\n",
      "}\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# get mean of age, using correct dataset size\n",
    "with sn.Analysis() as analysis:\n",
    "    # load data\n",
    "    data = sn.Dataset(path = data_path, column_names = var_names)\n",
    "\n",
    "    # get mean of age\n",
    "    age_mean = sn.dp_mean(data = sn.to_float(data['age']),\n",
    "                          privacy_usage = {'epsilon': .1},\n",
    "                          data_lower = 0.,\n",
    "                          data_upper = 100.,\n",
    "                          data_rows = 1000       # Here we're providing the correct value of n.\n",
    "                         )\n",
    "analysis.release()\n",
    "\n",
    "print(\"DP mean of age: {0}\".format(age_mean.value))\n",
    "print(\"Privacy usage: {0}\".format(analysis.privacy_usage))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Providing incorrect dataset size values\n",
    "\n",
    "However, if we provide an incorrect value of `data_n` we still receive an answer as we see below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-15T16:41:40.694235Z",
     "iopub.status.busy": "2021-04-15T16:41:40.693539Z",
     "iopub.status.idle": "2021-04-15T16:41:40.711013Z",
     "shell.execute_reply": "2021-04-15T16:41:40.711551Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DP mean of age: 44.666666666666664\n",
      "DP mean of age: 44.18181818181818\n",
      "Privacy usage: approximate {\n",
      "  epsilon: 0.2\n",
      "}\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# get mean of age, using correct dataset size\n",
    "with sn.Analysis() as analysis:\n",
    "    # load data\n",
    "    data = sn.Dataset(path = data_path, column_names = var_names)\n",
    "\n",
    "    # get mean of age using too small n\n",
    "    age_mean_low_n = sn.dp_mean(data = sn.to_float(data['age']),\n",
    "                          privacy_usage = {'epsilon': .1},\n",
    "                          data_lower = 0.,\n",
    "                          data_upper = 100.,\n",
    "                          data_rows = 900       # Here we're setting n below the actual value in the dataset.\n",
    "                         )\n",
    "\n",
    "    # get mean of age using too large n\n",
    "    age_mean_high_n = sn.dp_mean(data = sn.to_float(data['age']),\n",
    "                          privacy_usage = {'epsilon': .1},\n",
    "                          data_lower = 0.,\n",
    "                          data_upper = 100.,\n",
    "                          data_rows = 1100       # Here we're setting n above the actual value in the dataset.\n",
    "                         )\n",
    "analysis.release()\n",
    "\n",
    "print(\"DP mean of age: {0}\".format(age_mean_low_n.value))\n",
    "print(\"DP mean of age: {0}\".format(age_mean_high_n.value))\n",
    "print(\"Privacy usage: {0}\".format(analysis.privacy_usage))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's examine what is actually happening when these values are provided.  When we provide all of the metadata arguments (`data_lower`, `data_upper`, `data_n`) to the function `sn.dp_mean`, it works as a convenience method that knits together a number of library components to provide a mean.  A clamping, imputation and resize step are run on the dataset, in order for the validator to certify the analysis is privacy preserving (for more detail see the notebook \"data_analysis_tutorial\")."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-15T16:41:40.719385Z",
     "iopub.status.busy": "2021-04-15T16:41:40.718844Z",
     "iopub.status.idle": "2021-04-15T16:41:40.723981Z",
     "shell.execute_reply": "2021-04-15T16:41:40.724432Z"
    }
   },
   "outputs": [],
   "source": [
    "with sn.Analysis() as analysis:\n",
    "    # load data\n",
    "    data = sn.Dataset(path = data_path, column_names = var_names)\n",
    "    \n",
    "    # establish data \n",
    "    age_dt = sn.to_float(data['age'])\n",
    "    \n",
    "    # clamp data to range and impute missing values\n",
    "    age_dt = sn.clamp(data = age_dt, lower = 0., upper = 100.)\n",
    "    age_dt = sn.impute(data = age_dt, distribution = 'Gaussian',\n",
    "                       lower = 0., upper = 100., shift = 45., scale = 10.)\n",
    "    \n",
    "    # ensure data are consistent with proposed n\n",
    "    age_dt = sn.resize(data = age_dt, number_rows = 1000, distribution = 'Gaussian',\n",
    "                       lower = 0., upper = 1000., shift = 45., scale = 10.)\n",
    "    \n",
    "    # calculate differentially private mean of age\n",
    "    age_mean = sn.dp_mean(data = age_dt, privacy_usage = {'epsilon': .1})\n",
    "    \n",
    "analysis.release()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `resize()` step takes `data` and `n` as arguments and builds a new dataset from the original data, with the number of rows equal to the supplied `n`.  This is done by two rules that guarantee the privacy of the analysis:\n",
    "1. If `n` less than or equal to the size of the supplied data, sample `n` rows (without replacement) from the supplied dataset.\n",
    "2. If `n` greater than the size of the supplied data, create a new dataset by using the entire supplied dataset, and adding the remaining required rows by imputing new observations as if they were missing values.\n",
    "\n",
    "Thus, in general, we get the first `n` rows from the supplied dataset, and any additional rows are treated as if they were in the original data, but missing values.  \n",
    "\n",
    "(Note, we would get a better answer if we resampled any extra rows we needed from the dataset, instead of imputing them.  However, then an individual could appear in the dataset twice or more, and this would change all of our worst-case sensitivity calculations in the differentially private mechanisms.  We plan to add this as an optional parameter to allow a user to inflate the sensitivity by a maximum value, and pay the corresponding privacy loss in their budget, in a future release.)\n",
    "\n",
    "We'll demonstrate the effects that resizing can have on differentially private releases at the end of this notebook, but first we show how to deal with the case where `data_n` is simply unknown."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analysis with no provided dataset size\n",
    "If we do not believe we have an accurate estimate for `data_n` we can instead pay for a query on the dataset to release a differentially private value of the dataset size.  Then we can use that estimate in the rest of the analysis.  Here is an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-15T16:41:40.731918Z",
     "iopub.status.busy": "2021-04-15T16:41:40.731318Z",
     "iopub.status.idle": "2021-04-15T16:41:40.740106Z",
     "shell.execute_reply": "2021-04-15T16:41:40.739600Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DP number of records: 1021\n",
      "DP mean of age: 40.59745347698335\n",
      "Privacy usage: approximate {\n",
      "  epsilon: 0.15000000000000002\n",
      "}\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# get mean of age, using DP release of dataset size as estimate\n",
    "with sn.Analysis(dynamic=True) as analysis:\n",
    "    # load data\n",
    "    data = sn.Dataset(path = data_path, column_names = var_names)\n",
    "\n",
    "    age_dt = sn.to_float(data['age'])\n",
    "    \n",
    "    # get dp release of the size of the dataset\n",
    "    dp_num_records = sn.dp_count(data= age_dt,\n",
    "                                 privacy_usage={'epsilon': .05},\n",
    "                                 lower=0,\n",
    "                                 upper=10000\n",
    "                                )\n",
    "    # get mean of age\n",
    "    age_mean = sn.dp_mean(data = age_dt,\n",
    "                          privacy_usage = {'epsilon': .1},\n",
    "                          data_lower = 0.,\n",
    "                          data_upper = 100.,\n",
    "                          data_rows = dp_num_records   # Here we're providing the DP release of n.\n",
    "                         )\n",
    "analysis.release()\n",
    "print(\"DP number of records: {0}\".format(dp_num_records.value))\n",
    "print(\"DP mean of age: {0}\".format(age_mean.value))\n",
    "print(\"Privacy usage: {0}\".format(analysis.privacy_usage))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that our privacy usage has increased because we apportioned some epsilon for both the release count of the dataset, and the mean of the dataset.\n",
    "\n",
    "The DP released estimate of the dataset size is a noisy answer.  As we saw above (in section [Providing incorrect dataset size values](#Providing-incorrect-dataset-size-values)), if the DP released estimate happens to be too small, we will create a new dataset by subsampling from the `age` variable.  If the DP released estimate happens to be too large, we will add missing values to the column and impute.  Likely, the mean of the imputations is not the mean of the variable, so in this case, the answer will have some bias.  We can see this in simulation by plotting the number of records against the DP released mean through many simulations.  \n",
    "\n",
    "(Note, we're adjusting the epsilon for the dp count downwards to intentionally increase the range of released datasets sizes for dramatic visual effect in this simulation.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-15T16:41:40.747506Z",
     "iopub.status.busy": "2021-04-15T16:41:40.746944Z",
     "iopub.status.idle": "2021-04-15T16:41:48.342734Z",
     "shell.execute_reply": "2021-04-15T16:41:48.343165Z"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "n_sims = 1000\n",
    "history = np.zeros(shape=(n_sims, 2), dtype=float)\n",
    "\n",
    "# loop over the previous release to accumulate many simulations\n",
    "for i in range(n_sims):\n",
    "    with sn.Analysis(dynamic=True) as analysis:\n",
    "        data = sn.Dataset(path = data_path, column_names = var_names)\n",
    "        age_dt = sn.to_float(data['age'])\n",
    "\n",
    "        # get dp release of the size of the dataset\n",
    "        dp_num_records = sn.dp_count(data= age_dt,\n",
    "                                 privacy_usage={'epsilon': .05},\n",
    "                                 lower=0,\n",
    "                                 upper=10000\n",
    "                                )\n",
    "        # get mean of age\n",
    "        age_mean = sn.dp_mean(data = age_dt,\n",
    "                          privacy_usage = {'epsilon': 1},\n",
    "                          data_lower = 0.,\n",
    "                          data_upper = 100.,\n",
    "                          data_rows = dp_num_records   # Here we're providing the DP release of n.\n",
    "                         )\n",
    "    analysis.release()\n",
    "    history[i,0] = dp_num_records.value\n",
    "    history[i,1] = age_mean.value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-15T16:41:48.358006Z",
     "iopub.status.busy": "2021-04-15T16:41:48.350871Z",
     "iopub.status.idle": "2021-04-15T16:41:48.680883Z",
     "shell.execute_reply": "2021-04-15T16:41:48.681384Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABd2klEQVR4nO2deXxU5b3/398zM5msQALIKpuAgGJBIovIqqJFxa22aq1L9Vrbequ3tbZeb217b71L77X6s4sttQouqHVBoeKCGDaBIFGUVYKyKIQ1AZLMmpnn98eZM5mZzGRmksn+vF8vXkxmzjnznEnm83yf7/NdRCmFRqPRaDovRlsPQKPRaDQtixZ6jUaj6eRooddoNJpOjhZ6jUaj6eRooddoNJpOjr2tBxBLr1691JAhQ9p6GBqNRtOhKCsrO6aU6h3vtXYn9EOGDGHTpk1tPQyNRqPpUIjIvkSvadeNRqPRdHK00Gs0Gk0nRwu9RqPRdHK00Gs0Gk0nRwu9RqPRdHLaXdSNRqPRdDVKy70sK/NQURWgX6GNuROymTTCmbHra6HXaDSaNqS03MvrpW5umZXH8L52dh+qY2FJLUDGxF67bjQajaYNWVbm4ZZZeYwa4MBuE0YNcHDLrDyWlXky9h4pC72I2ETkYxH5R+jnBSKyR0Q2h/6NS3Deb0Vkm4jsEJHHRUQyNHaNRqPp8FRUBRjeN9q5MryvnYqqQMbeIx2L/h5gR8xzP1VKjQv92xx7goicD0wFzgHOBs4DZjRxrBqNRtPp6FdoY/ehuqjndh+qo1+hLWPvkZLQi8hA4DLgyTSvr4BsIAtwAg7gcJrX0Gg0mk7L3AnZLCypZecBP3UBxc4DfhaW1DJ3QnbG3iPVzdjHgPuBgpjnHxaRh4AVwM+VUt7IF5VS60WkBKgABPiDUip2VaDRaDRdFmvD9YU1rnDUzVWTclo36kZELgeOKKXKRGRmxEsPAIcwrfX5wM+Af485dzgwGhgYemq5iExTSq2JOe5O4E6AQYMGNelGNBqNpqMyaYQzo8IeSyqum6nAPBHZC7wIzBaR55RSFcrECzwNTIxz7tXABqVUjVKqBngLmBJ7kFJqvlKqWClV3Lt33CqbGo1Go2kiSYVeKfWAUmqgUmoIcD3wvlLqJhHpBxCKorkK2Brn9P3ADBGxi4gDcyNWu240Go2mFWlOHP3zIrIF2AL0An4DICLFImJt2r4CfB465hPgE6XU0ma8p0aj0WjSRJRSbT2GKIqLi5VuPKLRaDTpISJlSqnieK/pzFiNRqPp5Gih12g0mk6OFnqNRqPp5Gih12g0mk6OFnqNRqPp5Gih12g0mk6OFnqNRqPp5Gih12g0mk6OFnqNRqPp5Gih12g0mk6OFnqNRqPp5Gih12g0mk6OFnqNRqPp5Gih12g0mk6OFnqNRqPp5Gih12g0mk6OFnqNRqPp5KQs9CJiE5GPReQfoZ8XiMgeEdkc+jcuwXmDRORdEdkhIttFZEhmhq7RaDSaVLCncew9mI29u0U891Ol1CtJznsGeFgptVxE8oFgmmPUaDQaTTNIyaIXkYHAZcCTyY6NOW8MYFdKLQdQStUopVxpj1Kj0Wg0TSZV181jwP00tMYfFpFPReRREXHGOW8kcEJEXgu5ff5XRGyxB4nInSKySUQ2HT16NK0b0Gg0Gk3jJBV6EbkcOKKUKot56QFgFHAeUAT8LM7pdmAacF/ouGHArbEHKaXmK6WKlVLFvXv3TusGNBpN+6O03MsvXzzJnU9U8ssXT1Ja7m3rIXVpUvHRTwXmichcIBvoJiLPKaVuCr3uFZGnMcU8lq+AzUqpLwBE5HVgMvC3Zo9co9G0S0rLvbxe6uaWWXkM72tn96E6FpbUAjBpRLyFv6alSSr0SqkHMK13RGQmcJ9S6iYR6aeUqhARAa4CtsY5/UOgh4j0VkodBWYDmzI0do1Gk4TSci/LyjxUVAXoV2hj7oTsFhfbZWUebpmVx6gBDgBGDXBwy6w8Xljj0kLfRqQTdRPL8yLSGxBgM3AXgIgUA3cppe5QSgVE5D5gRWhCKAP+2swxazSaJJSWe3llnYsTtYqeBQa3zc6lMN/WKpZ1RVWA4X2jpWV4XzsVVYEWe09N46Ql9EqplcDK0OPZCY7ZBNwR8fNy4Jwmj1Cj0aSF5TqxGcK9l+dhswkLS2q5alJOq1jW/Qpt7D5UF7boAXYfqqNfYYM4DE0roTNjNZpOhuU6qawJcuYAR9h1sqzM0yqW9dwJ2SwsqWXnAT91AcXOA34WltQyd0J2i76vJjHNcd1oNJp2iOU6ibSsLYFvDcvaWi28sMYV3hu4alJO0lVEW+wndBW00Gs0nQxL4C3L+pZZeQSCiqJ8I+zCaWkmjXCmJdI6Uqdl0a4bjaaTYQl891yDeeflsOD9Wh5bWkMwqFKyrNuCyEgdu02i3E2a5qMteo2mkxHPdXLHxXntUuAtdKROy6KFXqPphKTrOmlrdKROy6JdNxqNps3RkToti7boNRpNm9PUSB1Namih12SU1gqR06F4nY+O5m7qSGih12SM1gqR06F4Gk16aB+9JmO0VoicDsXTaNJDC70mY7RWiJwOxdNo0kO7bjQZo7VC5HQoXjR6v0KTDG3RazJGa4XI6VC8eqz9ihum5fKnOwu5YVour5e6dUcnTRTaotdkjNYKkdOhePV0hiYfekXS8mih12SU1gqR06F4Ji2xX9GawqsjqFqHlF03ImITkY9F5B+hnxeIyB4R2Rz6N66Rc7uJyFci8ocMjFmj0YSw9isiac5+RWu7gnQEVeuQjo/+HmBHzHM/VUqNC/3b3Mi5/wGsTndwGo2mcTK9X9HawqsjqFqHlFw3IjIQuAx4GPhxOm8gIhOAPsDbQHG6A9RoNInJdJOP1hZeHUHVOqTqo38MuB8oiHn+YRF5CFgB/FwpFbW+ExEDeAS4Cbgo0cVF5E7gToBBgwalOCSNRgOZbfLR2sIb2Rwlciyt0RylK5HUdSMilwNHlFJlMS89AIwCzgOKgJ/FOf0HwDKl1FeNvYdSar5SqlgpVdy7d+/URq7RaNKmtNzLgvdrOX4qyAtrXJR94Ytyz7R26OqkEU6umpTDC2tc/GB+FS+scXXZCKqWJBWLfiowT0TmAtlANxF5Til1U+h1r4g8DdwX59wpwDQR+QGQD2SJSI1S6ueZGLxGo0kdy5KvC8Af7ujBnqOBsCU/YVgWFVWBNgld7fIRVIe+hM93wNQ5LfYWSYVeKfUApvWOiMwE7lNK3SQi/ZRSFSIiwFXA1jjnftt6LCK3AsVa5DUak5YOY1y0upY1O7zUBcBuA6dduOvSfF5Y42LP0UBUzH33XCPsnklVeHX8ewZ4/o+w8h/gzIEJF0B2bou8TXPi6J8Xkd6AAJuBuwBEpBi4Syl1R/OHp9F0Tlo6fnzR6lpWbfNy7ZQcZoxxsmq7l5fXudm02xvlFx/a20ZFZSBtv7iOf28iwQB8vB7OmQiOLBg8HC67AWZd0WIiDyBKqRa7eFMoLi5WmzZtauthaDQtyi9fPMkN03KjNj13HvDzwhoXv76+e7Ov//2/VHL1pBzmjKsX7588XUWNV/GXu4rC1vjBygB2G9w6O72esi09/k6HzwsfLId3X4WjFXDH/TB5dkbfQkTKlFJxIxt1ZqxG0wa0dBhjXQBmjIkW7msm57CgxMXOA34mDMuie64RtuQTiXx7CcPssNTVwdt/hxVvQPVJGHomfON2GD+lVYehhV6jaQOaG8aYzD9ut8Gq7d4oi77WqzAk9Y3W9hSG2eHwuExXjM1mumqGngmXXgcjzgaRVh+OFnqNpg1oTvz4otW1rN7uJRiEvoUG44Y6eL3UDdT7x6eNdvLqevM5y0f/6no3M85ycuP0vJTG2FjBtObGv3fajdwvv4C3X4YtH8J/PQ15BXD//4KzbSuraqHXaNqApoYxLlpdS8lWMy+xX6HB+GFZbCz3MXFEFsvKPOHzLTFfXOrm5XVu7DaSinys+FZUJnbPNCcMs9Nt5CoFn30Kb/0dtpWZETTTv24+D20u8qCFXqNpM5qS0bp6uynykXHwE0dk8fEXPg6dCEYdf+P0vJSt93ji+9jSapZucnP1pPpokEj3TFPj3ztDaeUoDn8F//cz6FYI19wKMy4zLfl2hBZ6jaaDsKzMQyBoWvKRcfCLVtdyqCpIv6Km+8fjie8l47N552MPowc6Muqe6fAbuVYEzfEj8I3vQt/T4Ue/htHjzZDJdogWeo2mndGYQFrumqg4+KogNoNmlSmIJ75XFOewrMyTcfdMh93IrTkFJUvh/SVmBM0ZYyAQMDdcz5nU1qNrFC30Gk07IplAjhvqCPvkF62upaLKdNdMH9O8bNZE4tu/yNakuPiW3MhtEz5eB3/9H9OaH3sefP2bbRZB0xS00Gs0bUhpuZdX1rk4UWtu3BkCl56bnVAgXy91h33yliU/fUxqkTSNTSKZFt/G3DMdphXkoa+gzg8Dh8KQkVA8DeZca/7cwdBCr9G0EaXlXl5a68IQuPfyfBB4bGkNq7d56V9kCwtfrEAuK/Nw6ESQ/kXphSU2ZmVbVnumxDeZe6ZdFzLb/zksewnK1sBZE+De30BhL/huvLqNHQMt9BpNG2CVC64LQM8CgxqvYtIIJz0LDPx1KipUMlMCmWwTNJ1rJ4uD75Dumc93wNJFsPVDyMk13TMXXtXWo8oIWug1mlbGcqEEAqaL9zszc3lulQuAKydm89QKF9XuAHUBlZZAJhPfTG2CphIH3xT3TJskUSll/jMMKN8K+3aZIZIzr4Dc1EJTOwJa6DWaVsZyobywxoXXr7AZEv75hmm55GeDxw8/mF+VUpmCZWUeKioDGAZcMj6bX1zXLa74ZsrKTjUOPt0VQqsmUQUDUPaB6aKZczVMucisIDnrinaR4JRptNBrujRtYUVaLpS5E7J5aa2LJ5fXcMvMXA5WBvjr8hpEhFtn56aUJWuVQrAZMOEMMyLH8u/Him+mNkFbIg6+1ZKo6vywfoVZpuDwAeg7ELJDlnsnFHgLLfSaLktLWpGNTSCWC8X6+ZV1Lh5fZr6vTeDaqamJfGQphIqqILsPBZg0sr4UQjzxzcQmaEvEwbdaEtXvf2WWKRh0Btz1IJx7PhjtPH4/A2ih13RZWsqKTDaBRLpQUi0XbF3XqiFvYZVCeGxpNcP72qJKIbRUElJzXEDpxvE3e/y11bDqTdMlk5NnhkdefLUZTdNBYuAzQcpCLyI2YBNwQCl1uYgsAGYAJ0OH3KqU2hxzzjjgCaAbEAAeVkq91Pxha1qbzlhtsKWsyGQTSFM3Kq3Jw0qUMgSWfezh6km54XIFVokEq6l3S0S5NNUF1Jpx/FSfgOWL4f2lZsng3v3hvOlw1rlNu14HJx2L/h5gB6ZoW/xUKfVKI+e4gJuVUuUi0h8oE5F3lFIn0h+qpq3odNUGQ7SUFZnKBJKKCyWy5yvAuCEORg1wcOhEkH6FBoN628K1aOaOz2ZZmSf0/kFeWONKuRpmZF/ZaaNTS75qiguoVeL4AwF49SmzD6vfBxOmwdxvma6aJtBZDJyUhF5EBgKXAQ8DP0714kqpXRGPD4rIEaA3cCK9YWrakk5XbTBES8V6Z2ICie35eveTJ/hkr59Fq2vDpRBKd/kIBIkqhTDr7NSEetHqWlZu85qRhQLTx2TRp4ctXMM+1aqX6ZDJOP4GRDb6OLAXzp0Kc6+H/oOaPN7OZOAYKR73GHA/EIx5/mER+VREHhWRRu9cRCYCWcDncV67U0Q2icimo0ePpjgkTWvR4asNJmDSCCdXTcrhhTUufjC/KmUrOBnWBLLzgJ+6gAq7UdIpOrZmhynyc8bl4Mwy6F9kY/pZWazZYTb33ljuY3g/GzaDcCmEdER+1TYvAlw9ycm1U3JYu8PHsVNBrp2Sw5od3mbcfWKsCTCSZq+gjh+GZ38P990ElSHt+NG/mz1ZmyHyEG3g2G0SNnCslVNHIqlFLyKXA0eUUmUiMjPipQeAQ5jiPR/4GfDvCa7RD3gWuEUpFTtZoJSaH7oGxcXF7atbuabjVhtMgZZIxc9EGGNsz9e5E7JZvMFNXQAmDMviYGWAdz72EAySdikEaxJ5eZ2bC8eaEwmYTUoeu60HL69zp3G3Jqm4ODK6gjp8EN560QyVROCCOWbSE5hWfQboTAZOKq6bqcA8EZkLZAPdROQ5pdRNode9IvI0ELcQhIh0A94EHlRKbcjEoDWtS3tPZ29JP2pr+rAjie35OmmEky17fRyvDoYTqW67MK9JPnhrEllc6g6/x4wxTl5eZ/5sT1MnU3VxZKyYWc0p+OX3zMczLjN7sRb1Tu8aKdCZDJykQq+UegDTeidk0d+nlLpJRPoppSpERICrgK2x54pIFrAYeCbJpq2mHdMeqw1GhhrakmSENpVYP7nVdxWa78NONjnF6/n64W5/Su6ZR5ecYvtX9S6S07oLD13XPTx+EXMSiXyPuoDZONzqK9sYsZOH0y7cdWl+Sns4TZ4AD+yFTz+Er18H+d3gth/DmedAj57pXytF2ruBkw7NiaN/XkR6AwJsBu4CEJFi4C6l1B3AN4HpQE8RuTV0XoMwTE37pz1VG4wNNbxgtJOSrd6EGaFNub41iXTLEbrnGTizjLB1vbjU3WShLy338sJqF7Ve00PZI0/iNvduSs9XgHv/VklthIs92wFHTir+9HYN/zLPDJh7ZZ2bV9e7uXZKDlNHO3h1vZugMr/IMxuZSErLvTxTUosvNId0z4FRpzso3eVn025vlOWbMRfH/t1mobGP15mbredfBN0LYdKs5l87Ce3RwGkqaQm9UmolsDL0eHaCYzYBd4QePwc816wRajQxRG6SHToRZPbYbAb1tofFvTkiEzmJPPJGNbfMzuWFNfUibLk4mnrtl9a6cPsU107O5vRedp5+v5Y1271MG+OMqlgJ6fV8BfjViyfCIt+3h3DohMLjh8JcCVv41vhnnmW6biyrfEYCl1S8JK3cLJh+tpN1O31s/7KOnCxYs8PHTTPqz2u2i+PYIVj0J/h0o5nodMW3zUqS+a3bi7U9GTjNQWfGdhE6cjzwr148wYHK6D18a5PM8qNGinu6IhP52dgMmDMuOxRpAYeqglErhHR92LGx8E6HWSzxoq/lYLcJd1wsLHi/Nm5z71SuvWqbl2BM+MLvbu3GgSrFI29U07ubcPRU/QHW+FOZRKxaOoHQsLJs4AvAtDFONu32c8FoJ6W7fBiicPtg5wF/810cXo9ZcyYnDw7uh6tuhtlXdqpKkm2BFvouQEeOB75/YRVVoe5Lp3UDZ5bBl8eC/GLRCf7rO4VhP+qss5307ZF+RmjsZ/P9v1RRusssDGb5sK+elM3BygDvbnan5MO2iBcL7/VDrpPwJt/wvnYqa4IoZUbPpEqkH74wD9w+s+IlwPz3avnJvO6MG+Jg815/+Jx0xh9ZS8dmQCBoirwA6z/z8U9z8lm0upbj1eYsYEgzXRx7y2Hpc1B1DH7xB8grgP/8W5eoQ9MaaKHvAnTUhKfSci9VtYphfQx+elX38ATVM184Vm3Gp1uhhq9tcBMMknYsfOxn07/IxtRRZmEwK1vz9Y1m3PTiUndKfvJIK94QOHYqGI6FP1ETwOMjvMkXCCgKsoVar0opzj52o3VQLxvXTc3lkTeqw8/t/CrAc6tquHlmDpsX1At9svHH7h8ATBuTxdrtPgyBLDsEFZxyK/Ye8VNRFSQ3C1y+1PYP4rJ/Nyx5DjZvgNx8mHMNBOrA7tAin0G00HcBOmo8sJWYcvfX86MSVp4tMUUt0oJMFmqYyHUV+9mY8eouKqsVdQHFuWdksWW/P+XJI9KKf3mdmysnZvNGaKKYOyGbF9e4CCpF8XAHT6+oobJGIQIzz0ruC47nwtp/LMBbZW56FhicqAkSCGn0qm0+Vm3zAZDnhMduL2r02qXlXp5bWYvXDwXZUO0xa36t2e4j1wljB5ubrmC6nxZvMO/J5YMxA+2NinzCENWtm+CxfzMF/srvmD547aJpEbTQdwE6ajywNREtWOnin+eaESPD+9o5csp83bK4k9GY6yr2s5k0whly03hSavxhXd+aRJSCSSMczBmXw+JSN3abcO2UnHCkzpa9PkrL/bz9kekW6ZEnfOP85GWJS8u9HKgMcs4QO1v21WEIBIMwsJfB9q/q+O6FeTy1orbBednbyzn5j1Ku+KeT9BrUnZsfnsWsb49tcNyyMg++AHzj/BxeWe82a+n0srFxtx+3F8oP1q8ivKFFQp4Tbpje+AQb674qW72TTesOsYhJ3Hj+1+C6O2Da17XAtzBa6LsAHTUeuF+hDZ8/yKd76/j9slPcOjOXP7xVA8CAotSqd0T2Zn1hjStsyVuuq3ifzcZyH7fOTi0ZydoMNQyzXsyqbT52HwpQWl4fp37leWaS0rub3SnHwsdirW5un53H/7xew2ndDT7Z4+fLY6aFX37QtN6twrv9imycfuAL3nq6BK/LVOaj+07yhzvfBGgg9tYkNWOMkw92+kK1dMxaOACVNfXunHQyccOlHHpWwJPPc/5HH3BO4UB+uv1r5mdwyTfS+hw0TUMLfRegvcYDp9Jg+vVSN70K4NO9dfx4gWnKF+YJv7q+R0rXt3qzWnXbLUt+wrAsKqoCTfpsYl0oPfNh9jk54WSkswfZo338H5rWeyo+/kRuDmt181qpO/y5jB1s59N9pqW9Zocfpx2+M6t+grptyJqwyFt4XX6eebCkgdD3K7Rx6ESAVdu94esP72enqsZPUJkbstPHpD9B9ao9wIUfvAEfrzUbbl/xbRwz5lG3qH27DTsbWui7CO0tHjidBtPLyjyIpB4WGhn73bPAoEeesOdoIGoTunuuEXZdpfLZxIsnz3NCrReO18C6nabl+so6N6u3+VCA1xekVzcjaSKSReQEYhMY3NvGqm1WFykbDhthv/vc8Vm89IFp5WfZ4eZZDVcgx/afJB5H953k3r9VRrld5k7I5rmVtbyyzs01k7OZMMzOu5/4UMqMEroxgYsm4WStFIhwWt0xVHkZXH4jXHQ15BewarMbu61puQgtQUcOPU4VLfSaNqElGkxD9ATyuyXVfGdmLn97r5b579Zw55x8hva2cbAy0KQQzIkjzAiffoVCRZXigtFO3vvUS1G+cKAyGE5GssIa737yRNKM1thYeEPg5pm51HoVr653M7SPjTU7vNw6O4/XS90h8feFzx9QZCRc3fQa1J2j+xqKfUG/Alw+eG5l/cQaXtmsdvFqaKPVkOSZsrGT9ZtvlnPGksX0GtIPrv0uPSedx305jzN3SC9mZDlZlWKIZ2uJb0cOPU4HLfSaNqE1ujv1K7RhM4Q75+Tzt/dqTPdMZQC7jZQ7IkWuDNbvNK3rn15ZwE+fOcX6z3xcMj6btz8yhbFkqwfDgK8qA9xxcXIff2yoJJjhiy+udXHTzDyunZLDa+vdBFTs6oaUxG/4d87n2G+Xo3z17yFZds6+/QLOm5LDaxvcUdm46U6qUZP10UOMeud5zvxwBT7DAWd+C4AbZ+SDSFqlHFpTfDtq6HG6aKHXtAnNjQRKJVzS2mi9aXouJ2oVt1+Um1JvVmjoh891wpfHTLN7wUoXF4w2N15zsghb469u8JDnlKTXj7XivzbYzichX/vpvQwOVgZZVubhX68pCIsjpC7E4Qmq3zD63DKLI6+sJ1hVQ+9B3Tn3e1PZ3nMo94VWH+lOrLGZvpt2exn1+Xvw/B9BBHXhVfzrkYt55PKh4XPSLeXQmuLbUUOP00ULvaZNaG6D6VTCJS1ReHaVC0g9merRJafCoYy3z87jgedO8uWxIHabGdL46d46hvUxyLbDq+vrm1Aki6aJZ8EDfLqvDgH6FUk4iuZgZYD575kRRtNGN8119cgb1fzowfH8ftxIsh3ww7kFDO1t4+4nT7BquxebAX16NLHr1eke/vvlSlZvL6Sw/2Aum3YpXHY9u1zdyV/jSvma8WhN8e2oocfpooVe0yY0JxKoMYsvdgLpnmsgQqOulNjVgbXh+v1LCrDbhBtn5PFsKJkIoFsOfHG43tovzBN+e0thwvHGCnz3XDgZoYWjBtrYeSBARWV00ZpP99YxoMhosjVs1eoZN9TB5j1+/vx2NeOGml/5V9a7cdpJmo0bmyMwbYiXOXuWwF9f54fDJvKbHnextKIvl/zTD9tVK8ZU6aihx+mihV7TZjQ1Eqgxiy/dCcQq3BUMQt9Cg3FDHWZ54myiVgY+X4BnVtVng5oWeON+8sh6MZGcdJmRMvk5QmW1YucBU0T7FwoHq0yxN6RxX3YqrqvIWj2b9/hxeeGDneaEk5uVPNkpauXUw8c//vtFLvnkTQi44LwZ9Jp3E9edyGVBiSvl5LJUaE3xba+hx5lGC72mw5HM4kvHl716u5drJucwe2x2lAvI7SdKbDaHfOipbLJCdDG2WLIdoQJkoZetpKSDVapBLHy8MS9YUUtdRCUEh6HC9ewjPxtrkrBq9djS6JAFMSunV57hqoqX+XJgMUsHXMcP/ulrANQeNvcQnvhe4yUW0qG1xbe9hR63BFroNe2SxsLrMmXxLSvzEAjC7LHZUbV0fv9mNb466NVNeG5lNYdDEYoDioyUsmVjrfjzz7Sz7rN6141VZbJXN4PKGtNNlOcUbpjeeCmE2A3igUVmcbd9x4IM7mWwrMzT4LNJp1ZP5Gc+sLtiRPl7DB8/BhgLc67lHWMCrx0eRNBv5ghYHatSreaZDunkNnTm+PdMkbLQi4gN2AQcUEpdLiILgBmAFagbt3OUiNwC/Fvox98opRY2a8SaTk+y8LqmWHyxfnKnHXx10K/QiFodDO9rx19nulZ2flW/+TdmoD3coSneeF9Z5+JEAgt+75EghblQFbNHuetgIKVrxyZqAUwb42DNdj+De9s4ejLAvmNBRNKzhuNVqxw/2OChAZsIvPEcWSeO8PmbX+eMM8dCtx5ccs04dr9VzeY9qeUItBSRn3fPAoPbZudSmG/rlPHvmSIdi/4eYAcQ+Rf508Z6wYpIEfBLoBhzoVomIkuUUlVNGayma5BKeF06y+1IN4ohZvclyxfep7stygJ+f4sZCx8v0zQei1bXsnJbfU2YUQMMdh6ot7pzHHCwKsj5ZzpY91l0OYJU3DSRETQWPXKFiSOyMUSikqfSzfS1qlVeOzmb97d4GHb4I65a9jw2bwW2ISP4YPL3ePboKO6NaCjy1fHkOQItaWlbn4nNEO69PA+bTcKrufYa/94eVh4pCb2IDAQuAx4GfpzG9S8BliulKkPXWQ5cCryQ5jg1XYhMhdeVlntZtLoWV8iTMmd8Ft1zbLy63h3e+Nyy388l47NZtLqWiqpguKZLKkIZ2T9VMC2Zo6HKmoaY8fXukLZHinyyssGxLprHlpoi36c7HD4JLq9iYUkt3zw/m1XbfDjt4K1LHkFjjTtyhTBppINLx+fwWqmHqwZXor6y8efh/8JdP5vDpCAs/EtVWhvbVn6A9TkWD8/KqKVtGQG/W1LNmQMc2G0SFvhfXNet3cW/t5fM21Qt+seA+4HYho0Pi8hDwArg50qp2BCDAcCXET9/FXouChG5E7gTYNCgQSkOSdNZyUR4XWwbPICy3X6unuwwywZvMDcvA0HYvMfPoRPBlKoyJoqFz84yuzw57YqskFsoHo25aeL59616OgA2w8BuBPEFoGeB8OR7pi/IWxd//6Dk+S0882AJx/abJYqn/vM0Pu8zNCyUg2q+YMbbL/L5qQvpV3g+xyddzuPHZxPEABF2H/LTr8jWaDnoeK6lonxhRH87a7Z7MYSMWtqWERD5N2IZAe0x/r29ZN4mFXoRuRw4opQqE5GZES89ABwCsoD5wM+Af2/KIJRS80PXoLi4OL6jU9MpiBQGuw0CgYZhis3dbI2MprEaeQvmqiAy4xTMkrup1LVPZMFbnDPEwUef+8MuIav9nkVjAg8NRf78UQ7W7fRT64V8J9R4TRfQ6b3MVoqfHUi8f1Dy/Bbm3/MO1cfrC4cd3XeSN/71ba5/5OuMsvXkni+f5KwjpXic3Vi218/cK7J5aqWLbKcNn1+l1JIxtgbQOUPsfLq3jkkjsvjwcz/njXCwZoeXb07NzZilbQl85N9IIKgoyjfaZfx7e8m8TcWinwrME5G5QDbQTUSeU0rdFHrdKyJPA/fFOfcAMDPi54HAyqYPV9ORibWGB/eyccJlJvRY4YFN3WyNbfBtRdMsLnXTM8/g8KkgG8v9KOCRpaZ/xWFL7u5IZMErzMSnbrmm8H5Y7o9q0h0IJo+FT3RtgO/MyCcQrKF0l5+aCCPfypx1bdyFd9kGvMeqqRvUnZJqs6FIyfNb+MOdbzYoTwygfHXkvPIUfLqbUfYslg24hrd7zcVty+VHDvD6VTgpLJUsYstafSGUCXvnRfn88qVTbNzt59bZeTy/qoa6QGaTnSIFft55OSx43+xbW5gnXJtCA5fWpr1k3iYVeqXUA5jWOyGL/j6l1E0i0k8pVSEiAlwFbI1z+jvAf4qIlTY4x7qWpn0RG4GRTKTSZdHqWrZ/VceMs7LYdbCOkf3trNnu45zBZtZmUzdbYzdDc53g8prjX7rJzbTRTlZt83JaN+HIKfOgPYeDpkshSXORSF+5wwB/dBc/vH6o9ZjXjBT5VEIlv/dEZdQ5kbg37uLW/yzl5MFTSI98CuZNIXfiSOw2uHV2Hq6Nu/jDyyvjNhR55sGSBiKfn+WnLmjgqbOxfa9BQdbpLFx7GqrQjePrX+E4dySPL6vFkMbLOMSLLqqqCVBRFcBuwKrtXq6cmMNTK2oJBBWHTihsBhm1tOMZAanmNrQF7SXztjlx9M+LSG/MVexm4C4AESkG7lJK3aGUqhSR/wA+DJ3z79bGrKb9EBuBMaDIxt9WuFgZciVkQuzX7DCvdf0FefxgfhUPfTOP07rbzKbeqmnLWUvkUXD1RCeGIbxW6kGAoX1svPOxh3uvKCCoYM128/1TmcBKy728us4VlfD0w7n5PPaPmiiXjccPBTHf12Qumnh9X3MdYGmza+MuTi4qAb9p6auqGk4uKsEQuPtfxzNphJMbprwTt6HIEz9ZgetwfXSO0x5g3pgKrj3nK17f2p8XNw9iw76ebNgXOsB9Cnm2hCtHO7njJ+MTjtn6TF5a68IQuPfyfBB4/B9mRdAeucLg0+y8ut7NtDFZdMuB+e+YdXqcjtQqhaZDR0pwai+Zt2kJvVJqJSHXi1JqdoJjNgF3RPz8FPBUk0eoaVEifcOGmC3jLj3XyV2XGvz13RrW7PA2SehjQ8rqAvUx69Zy1qrf3r/IltZyNraCIsCc8bnYbYJhCK+ud7P3SIBAkKhomlkprFDu/VtlePPTQoBdFXUM7m1j39HoyehoaJVgNfiOd/144h6Jy282GgkoqF6yPizyYfx1BN8tZdLTkyl5fkuU7z3qOoersfcsQFWeZM6Zh7lh/JcU5frZsK+I9ft6NtxYwHTnfPD7NUmFflmZhyy7cOvs+o3FS881SzTbDPjyWIChfaJr5Seb9LoK7WFi0pmxXRgrHM7iyonZvBFKl//m1FyqPSrsEkmHeCFlj7xRTa7TjHmeOCKLBe/XMqCn2fd13FBHo8tZK3rk6P6T2AsLyLtiMn1nnMkpd/3g/vJuDT/8ekF48kCZrfhSjaYpLffy1Hu1BBXYDaJKDIjAOx97+OfL8nlsaU3Uecli4RsrhWBdWylT5AGCVTVxjztVYe4rPPNgScJr2fKdXPermfR55VEuGn6YbYe68V8rRrHjSDcKeuaAij9BxHaiitfO0FptRW4sXlGcw7IyDzUeRY1HcdKVWg0gTeujhb4LYzVu/mCnj0MnAthtYoYelro594wsCrIFly99pY8XUmZVUJw0wsHHX/g4Xh3keLWpppv3JE7RL3l+C49+9x8EQ+EudZXVnHqhhLnF2aw0BpHjgMF9bGze42fnAT97D/sxxBTQdBt8W9w9N5/fL6uhMM8sMRBUgII//KNehBvr7ATwT39KzUMZO5Eahflxxb7XIDMyKFF7wK/1P8FX1XnsPVxHj8u+xX88vZbSHdkYhQXcs+BC1NdG8McL/0ygMvG1IaYU8RhnuMxBThbkZEVnEe8+VEfPAgOnQ1KKXNK0HVrouzB1AZgxxkn3PCPcL3RecRZ1Afjz2zW4vIqZZ2emuuT35uTz/b9UUfaFv0Hj60SUlnt57EfvhUXeQvnqePXXK7lw4e1s3utn32HT2pz/TjXVofLws85qfLkcrziYxcptXi4Zn807H3sYUCQcCJUPrguJcjKXRKoiH7tyABh121Q+//P7UT54Z66Dmx+eBTRsDzi0qJZbz9vLhIEnWLylPwv+tIafHbmXN/Inc0EfG58fCvD3WrCXupn4g2l89Lv3El67tNwbduN9sNNH9zyDOePMVdar693YDMWTy2u4bXYeCDy1ohal4OrJjW8stofM0K6OFvoujN1mRkpYX+YXVrt4Y5PpY3X7VEoNreORKKQs1Zh1qN9ojWeBAviOVfO9S8zJwxOaB6o9jTeytojnTsl1Ci6vwm6YK4yZZzt5qwwOheLie+QJ32gkfC9RTZrGiBR5uwFP3FUETKTk3JyoRKfe105m0ckBLPpTJb6ZE5EXSujpqOWmc/cze8QRan12niwdwj+29yMQrGHHlz4KsoU9hwMxUTRFlIxyhq+dc1oB2XPNa7/4RCVZITX4f9/tzv7jwXAGp+UO+9YFubyyzsVjoZVNYZ5w7dTGo4vaS2ZoV0cLfTugrSweq145mF/mkxOyw9UImxNp09yQsqgN4gSuDKMwnx1f+uiWI9R4kk9KiWLWzxpksG1/ELuhEOrFd+H7tShM33meE/63kcYi8TZwXRt3UbNkPYGqGozC+hDJeMSGNM769lhmfXts1LhnnJXFNZNyeO2ss3kLuP7zJ5g+7CiLtwzg758MpNZnD38uVqhkvN+jde1YF82/Pn+Sao/5Gazd6WPOuPraMSdrze5aTdlUbC+ZoV0dLfRtTFtaPJYIpNO4ORWalfBUGQgHhohAwbwpnHqhJKrBNQ47BfOmNCpokXz/z5VxXTQA+dk2Tu9lJiJJxPOWxd9Y96i751fijZPvFBsiGQyFSAJhsXdt3IX7zQ34j1ezfFB3+j08Kyzu0DBSZ88BN7lr3+Om4Wfx4QUjWcJN/H3JWo5URXyFHXZ6XTslpW5aBysDdMsRuucZOLMMqj2KayZn89oGT3jyv2CUmfHanFLE7SUztKujhb6NyaTF05SVQbqNm1MlneYfz5bUNhBMmwE5WULPC89kP1CzdD2BynrruNuUkSlttt6/sIq6oBkNcvvF2Ty53BP1eukuP1dMyOLLY77wBJPnTN59qTE/fKIQyeol6ym+Zgzjqvfxh5dX4otIeHrkO6+z44MvmfAvF4b3DpwO8PoUV9tKOW/tC+A7Cpdexzem3Mgzvgn0Pa0Hx/+2Nvy5jLptKr99ZGLCcUV20wKYfKYjqmFJ/0IbSsHMs53hyR+aN/m3l8zQro4W+jYmk5UaO5ovdNHqWlZu9aIwy/lOP8vJO5tNH0ivfIOj1UFshnD67JEcj3B7JOvRCvWTnmWZK2DCsByexIPDBv6Ij3dpWX3sd7JomsZWBxaJQiSDJ2r4l3nduG1Iw+xVFCx7oowPXEVhq3+2czfjtyxgqOsLjnUbxH+7b+Kz+45zdP//YfTIx3njVPr85pZGJ/V4ewf9Cg0OnwjywQ4fM87ODjcseWpFLTbDDK0994yssLutOX8/TXHj6c3bzKOFvo3JlMXT0XyhsUW8Zox1smm3n7xs8Pvh8ClTTUXgeEg3k5UWaFBLp3f9ZyiYJREGFBkNkpdcG3dRvWQ9wRM1Zt2YQLQbJVksfCyJ9hV6JwmRBKhdup7ciSPpVQC28k/pUXeCpwbfxZe1ffnw+Wh3UMXT73P3/BxmXT+2wXXiNRUBs7Lk+GFZvP+pB5cX1m73csqtOFkbxOUze8lmsv9rum68jmiwdAS00LcxmaqF0VF8ofHKCxTlC+t2+rhgtJN1OzzUhoZsN+BkrQrXNm+sBovl/snNMqOJhvQxKylCfY/Wt8o8fP/SXJ542xV207g27qL6hZJwCGdk3ZhZ3x7LrbdvoHLxeoIpbKpaFMybEuWjh8ZDJC0Kc3x8e/RmdnzYg+3nXUnevG/w0MbL8ZDF4UcXNnAHeV1+nnmwJGpSgujcgIJsCAK1HnPSHNHPzuY9fn44t4An3q4OJ50tLnUnzOxtLuls4nY0g6WjoIW+jclULYxMrAxaeslsWWuGYebiF+ULlTWKwnyDY6eCbPjMywkXYddKXZCkWa3WNb11Znu9iSOyeeSNaorPyKIwz2DVNh8eP/jLdlG5eD2/+GG9YP/owfH8+bcbG8Tpe11+HvvRezy5vDbppmo8uk0eyR0X55nZvPtOYtgkLMoA/lkTYcHy8PFOe4Brxh7g2rEHsBmKis0bcanRbD19FB6ycDoSu4NiVwdWieb8bOGUW3H7RWaNnrxsCAZh0+d+gkEY2tuGy2vuhdx2YfK9jtZyp3QUg6WjoYW+HZCJWhiZqOHeUkvmSD9xzwKD49VB+hUajB+WxfLNHj4/FKD4DAebPjf91v5A8qSkeAlP5Qf83DyzgJ4FBm9s9PDQdQWs2uZDPtnFsWejBdvz8kpcF+fFtawBApU1jW6qJhL6+nGbHaQiSwYf3XeSR2563TwwNJtNGXycu6Z8Qc88H2u+6MnCTYM5VJ0Dm5fz/pL1FMybAhNHktWrAN+x6gbv12tQ97i1f065Ff0KDWw2oWeBgccXxB3airAZcPeTJ4DUu2m1ljtFb962DKKaUsykBSkuLlabNm1q62F0SJpjdf3yxZPcMC036gu284CfF9a4mpXeHtv39N7LzfICYwc7+Op4gOLhDt7+qN5X31hxsMhrPv1ebbg+zJ0X5zB/uRkhMqDI4JLxZqlcqxGG+z+f5cSBUw2uU9Azh+pKd4NCX0B9L8AEWH74SHfOX38Q3R7wtiGPJ5xIzDcVZp5xhCvGVPBk6VB2HGk4sRlOO//yt8sBGtSZd2TbKbxpFrbxI+mWI1wzOZsFJebnkO2A2edks7HcF/6MIztfNeYOi/07cnmD3H5Rfsb/NuKRaFJpi4qPHQ0RKVNKFcd7TVv0nYjmrAxaaskc6XPtX2TDZpNweYFLxmez4TNT5JP54aFhBMmgXgZ1QaEg1x6uLHmgMki2Q+GwEfbRxxN5IGEVSKBRkYd6V0qwqoZTL5Rwx8V5WFa8RTyR79/Nza3n7WXH4W4s3jqAVZ/3ZtXnvVFRUfwR7+Ot45kHS3h6748AwsXdbD3yyZ03Bdt4c2UxrI+df5R5KcwT+hYa7PgqwNodXs4flcUHoRLNlsg35g5LVJCuqiZ6A7ul3CnNbTqjo3Tio4VeA7TckjlyArHcSzfNyOWtMijd5aOyRiUtLwANY8ABatxBvjbUaV5zem44Nf9Pb5sdjyw3ym2/jb/5mSmUr67Bpuj3nqiMir7Jy6rjhvFfctnoCuqCwtZDpiWcSOAjsfzws749lorTh4VDUqeOtvPBDlO9N+/1M26Ig71H6zhysr55jLVaSqU0RGm5lwXv11IXMIXWEkzTFeZmypn157akOyUdg0VH6aSGFnoN0Dwff2MWVeQEYj337Eoz6sXpkKTdgSKjdLrlmBPCGxvdHK8OYrMJnx2s46pJOTy7yhU+J/aaNz88K2F7vUxxdN9JrjD+w1RUBFUbKrwDTB1yjB9O/Zx8Zx3Ld/Xh2bJBnHBnpXxtq7pkg2zZQ2bf3TGn29m6r44t+812hhOH2zleHeSkS6VUohnqBTMQgD/c0YM9RwNhwbxyYjZPrXCx84C/TbskxUNH6aRGykIvIjZgE3BAKXV5xPOPA99VSuXHOccBPAmcG3qvZ5RS/9XsUWsyTlOjf5JZVLETSPdcA5GGYhzvuotW1+KKqCEzdVQWSz50c85gByu3eUNNPwI4DMUplymAA4qMuNfNyrG3qNBDqORwuOiNwu5xU4dBpSuLz4/n8dTGoeypTDN8USDv8skNMnGddjhYpchzwtZ9deFSyt1yhA93+xttCQgN2wIaYjYS2bzHz56jgSjBvGFaLoV50uZdkuKho3RSIx2L/h5gBxDeMQq1DWwsRfE6wKmUGisiucB2EXlBKbW3KYPVtCwtUbSqqT7Xl9a68PjM1oavbTDLFqzc5mXWWFOMbpudy1MrTCv+iXfM/+NltTbWLDsWozAf22k98H/2VVqfQSxDi2q4Y9JevjyRw5/Xn8GOI934xdtnN3pOQc8cfO66uBmz+xd9QMG8IH1mjKQ6tK0wbqiDj77wU+uFM/qa5YgBajyq0ZIFkSskQ+DS8U5GDXTw2NIaVm/zct7wrPDEPLS3jYOVpmXfHhtvg47SSZWUhF5EBgKXAQ8DPw49ZwP+F7gRuDrBqQrIExE7kAP4gPg7Y5oOSSoWVSoTSGyIYLbDtJAv+loO63f5Ob2XwYflfjZ85qWqVlHtVinFgD/xkxUpW/LhRKcmUpTr5aZz93PRyCNUe+2s3dMz5XNrKt38+NmrwrH3kUTG7w+68EyOVStKy/18bbDwyT7FliU7wpm9vUMF0iB+tqyVx1CUDxeek03JVi8De9npWWDgr1NhV9gLa1xUVJquofZivcejvTTfbu+katE/BtwPFEQ8dzewRClVIZJwQ+kV4EqgAsgF/qVFm4Nv3QQvzW/4/J0/h9OHQdlaeP2Zhq/f/Svo0x/WvwfL/t7w9Z/8F/ToCavehPfeaPj6A49Cbh68+yqseafh67/8E9jtsPR52Lgq6qXagI3fnvPfVFQF+M6xlyiuKSMny6g/IDcfHvid+XjRn2DH5uhrF/aCH/+n+XjB7+DzndGv9x0IP3zIfPyX/4Kv9kS/PugM+KefmY8f/yUcrYh+fcRZcPM9AJz8j/vxHjsebhxSkC3clnsOuw9937So/vNecLsI1Cl+4wrCL2ww/ny45tbQ53AXUbupAJNn82hgLrv2u/nFzn+NasZhCBxf/HXmTriCt9cc5Zfbfhk+TbbBJKfQfcy1MOLrUHkUHn0w6tKHTwQ4Ly+PVfRmYHcX/3phzGcDLNw0mNL9PRk+NMiPXf8LV0SP7y8bhvHJwR6c1eckP5z6eYPzH187nJ1HunHugCr+7aIdiMDrW/vz0ubTw6WDU6HXoO5UnD4M+/0DMP5tYcMEqVD8/onJZpSNCHyyT+HauItTi0pQ/viZvY3VyZ89NptBve28sMbFlRPNkNRqd4AJw7Lonmt0iLDG9tJ8u72T9C9RRC4HjiilykRkZui5/phumZlJTp8IBID+mC6eNSLynlLqi5j3uBO4E2DQoEHp3UEkObkwYHDD5x2hja/c/ASvh5Z9ed3iv24LfUz53eO/boSEuVth/NetibBHTxgwmOPVQSqqAnh8iqBhY9xQB7+4rhvH3uhL+YcDGdDNRs+C0DWzIyyTot4Nr1/Qo/5xzz7giQkZ7Hla/ePe/UDFCG3vvvWP+/SHrJhNwiLz/EWraxng6k2eI4fsPKFHnkF5TRBnn95hi2pkn9OpOeVi75EAAwbaoMAwJyKLAYMbCP2Gw3ls99RhQ6jsdjo98gwOnwggoTD23ftyGDfFYNxwJwd3DAQFDjsM7GnHW1nN4z9Zw/KyTQwbns2D1/bgtMHd2bbfj9unwAm1WaY17wsY7KvKjf3NmGLssGOffjb79u9q8LrbZ7oAXH5b3PM9/vrX3yvvw2tb+psJT2lgc9iQiyeGa/8kLIpWVUNdEE7rBn2LQjkCb24Ii7yFlYVbcfowVm/3Egia4auBIBTlg7cO3F6z9s8VxTlUVAUozDfIywavP7O1blqD9tB8u72TNGFKRP4L+A5QB2Rj+ui9oX9WzddBwBdKqeEx5/4R2KCUejb081PA20qpOGazSWdPmIrcvFy0upYLRjsp2eoNf6laKhGlOZSWe3l6RS3XTM5h9tjs8PJ44ogsNu/xM3dCdtpxzLE1b26ekUPvHnYWltQysKeNT/b4UZg+5A93+zleHYyqPZ/I9y552XT7xrSouu+xdWcwBMlxmpExocQoyXOC24dKEj9v4ThzIP4vDjXMnE1CQU9zErBi+GPHC3A4nkWPuYfQ5ze3hH8eUGRQduPvEzZw7/fHH5rnCZw3wsGm3X6KhzvYdaAOj1/h9cPVk7JZuc2HP2A2gv/WBe3TF69JTrMSppRSDwAPhC40E7gvMuom9HxNrMiH2A/MBp4VkTxgMqYbqMsSuXl56EQwavk8aYQz7YiB1kgWWVbmIRg0l/p2m4Q3XBetruXQiWDKFlW8RtxgCtGrGzzcMD03fF3rkLc/Nq3cwjyJ2hB85sE4pX4BVeuJqkeTO3EkDjtUvmYWJrOEXSlVb+YCKrZFVAIihfnEi6twr9ma0nlgFja78/9dwqxvj220Ima8omhWs5XY0hC3/Tx+joCtyAyC+9HcPB5fVsu2/X7GDjbFPhCEW2flsKDEzauhje5U2gJqOi4Zj6MXkXlAsVLqIeCPwNMisg2zUuzTSqlPM/2eHYnIzUsrYiBS3NOJGGitZJGKqgB9C42o6AZzzEH6F6U2Vqv2fK5TcDrgrNPtrNnhR4DTe9nYfzTA39e6+I8bulFRZYpvYyGCjZX6tfzZBZNHEgiC49yRFNQRLZ6u1IQ9lr6/vT382Lt1b8rnFfTM4c7/dwmrbKezKCZUskcunKhPAwhb99VLoqtmXvbdsQ0+j5sfnsXjd/wDv6d+UpAsO/lXTKFfoYHDYYQ3Wg+fCBAImvsrC0vMjmKpNG/JNDqTtfVJS+iVUiuBlXGez494vARYEnpcg+nL14SIDAezIgZmne2kbw+DnQf8aUUMZDpZJNEXsF+huY8QGd3w/hYPNsOMekhGZO15j19R64Ubp+fzyd4TVLsVXx0PcM4QO5/sreOep0wBHzPQ3mgceKJSvxbBqhrLWAcSdH1KE6MwOlUkkS89kt6DzebeVWcMZ9FJMIsGR3PCBVk28EUs5PJCqxGIbBzekNyJI+lz62xqlm7g5MFTSA9zUug/60zGDzNDJc8bYda6OeU2QyrnjDNr4LSFD15nsrYNOjO2lYkMB5swzOzJ+doGN8GgGTmQzpcvk8kijX0B507I5vVSNxNHZJnumqoghpG88mG8pKerJmXz6noPf3m3hmsm57CgxEUgCF8eM8ecSs0bgOHfOZ+j/7M8oXg3RZQbJeQ6iX2PRNd15jq4e/5lvOEZQFUKiwdfzK/Mcuo47fCHOxuKfFRF0PPP5LYfj2PSCCe/fPEkNe4gdQHYWO5j4oiscD0hMDe4N+/xt9lGq85kbRu00Lcy8cLBmloPPJPJIo19Aa2N4WVlHg6dCNIvhbT6cJ14v5n0tHijB6XA7YNJIx2U7vLTPdeMRjIMUqp5ExVr328Y3W+cxcmXVzd0w8SIst2A7v27cfJgGikcERu2iRqOFMybEtW0JPx8zxwufOBCFp0ckPr7RZBI3C0iJ+XfLanmOzNzeS5UAmLuhGxeXFNLjQcuGGOKfGWNSloVNHZzPJX6Q01BZ7K2DVro24B0w8ESWdtjBzkyliyS7AuY7pitieORN6qZdXY2lTWKkq1e3vrIw11zcind5WfVNrNAevccSZp5GdkmsG8P4dAJRe7EkXSfMpKT63c18GfnThwZ1Vu2pGB2gygdsQlkO8N1aUSkUWEPnwfc0P0Az6w0m5YYNiEYUPQe3J2bH57FW3UDKY3YaO2WA6fiFMrs0x0Ox3if8pzw2O3xRT62rv9JV5B+hTZshkRNygcrA7z9kSdc0CyZaFuZyCJw7xX5oODp92t5aa05eWTaDagzWVsfLfQZoKU3lxqztsNZjM1MFsn0F9CaOOw2WLXdy43T8wgqxaptvnDJgmQ9YC1Ky71s/6qOGWdlcf0FeeHSuXnZZou8vjNGNhDl2I1cq7LkMw+WcGz/SWxFBeRdPjl8Xq98OJaCd2dAkcGMwJdRk0YwoDCy7PhmTgxZ8dHRNPFEHqJFPpHAR4q7zYBLxmdzsDLAd2bk8txqV3iyv2l6LgcrA+w84GdjuY/vXpT6JuuyMg9ZduHW2fV/Y3dcLCx4v5ZlZZ60/56SuQF1Jmvro4W+mbTG5lJj1namkkUy/QW0Jo5po528ut5UunOGOPj4C3/SeiwWsVmdQ0+zhcM7xw11sHmPKbSRnpt49W4sNhcMxn7/TfSN81oykY/cEL1tyDMNQjuDvsY7TzVGom5aVmnmQNB8/wnDHWws91GYJ9hsEjXZW9U7093nAfPvSymi/saG97VTWdNw4zgVUnED6kzW1kULfTNpjc2l1ljuZjqVPHLiCCp4bYObQJCopKfGiO1MZTPg7+vc2O1mdcrvzcnn+3+pAsySCXYbTBvd8Lolz2+p792aYnPvWGKFOFFoZ7obvokmpdhN7B55wolaxe6KOiaNdLJupzdc17+iKpByRdBFq2tZuc0bTrASYObZpjvF61dRf2O7D9VRlG/gdCSvlx9Lpt2Amuajhb6ZtMbmUmstdzP5BYw3cTTm0opNphLg6xOyw52pTutu8MleP39f62LcYAcvrzct2MZi7Uue38Kjt/+DoDe95t7QuK88UWinFeljtj6nwWPr5/k/SLzRGltz/vRewkmXOUEO72fn4y98nHQpbr8o16zrr1Kz4q08BgVcPdEJIrxe6qFkq5cxA+18eTzIX5fX8N0L88I++qCCqyen/zem/fDtDy30zaQjWtutRaoTR+RGa2EeDO/v4MNyPys+9YQbZ7xe6mbsYLO+i9XYOlms/e/vXREW+TChZKrC80divTS4N+w7Wn9IY+6fR5ecwjdzIiTIXIVoYY98nEqoZCQi8OUxxeDeNoLBYDirtV+hkbIVb7FmhxcR+MaUHOaMM8XbbhNeXe9mV0Udt87O49V1Lh5baq5KeuQJ32pipqz2w7c/tNA3k45obbcnrI3W3Cy4c04+NpuwsKSWbIfZ43RZmSfs1128wfT1J+uaZLk+fMeq474erKohUv8jRb6xFYI1IeXFZK7aCvPJb8QlZAj85fuJI2leL3WT54x2kTgM8AdNsd93NHoCOFQVTNsXb5V/njGm/vgZY5y8vM5NMNDyq7mOYJh0ZrTQNxP9R900Yq1Ylw/OHODAHtpkfHJ5NR4/VFQGqAuotDpTWU1LEiU0xSZTQephjblOweVV4To6NoGAgix7ffPtSCJDPCOJdVUdCz3fpzu4fUJQKWwB8MSU80nWOSpyvJEuM7vNLBy6ars3bNGv2u7FkPriq+les7HfQ2c1TDoqWugzQFv/UXe02iGxG62GmC6OpZvcXD0pl+F97Zx0mdaszUivbK4VKlitFAXzpkTVagfiZrg2FvVi+bUtLjjTwbuf+nDYwB8wRR6iRT5ZwlOkqwpgwnA7ZbvNn22GQVAFqfHARWMdvLfFVPrGNrFj/frmdcwM4+JQx6iR/ezs+KqOl9e5qasLhn30CpgxunntIjXtHy30HZyW+hK2xOQRL+HH2mjdvMfPWx95OKOPjc17zEQqpx1umpncgo8cZ0VlAMT0Y2fNHc0OGhYHy504MqG4W9w9vzLKvZPjALcf3tviY3Avg33HguHnImnsupHZpyJwyXgnH5b7KT9Y75o5WBXk9F4GtZ5gWOTtBtyaIHvaEvme+VAZKs4ZUGaOwpqQxW5Fgc0828nKbV4WbzTDeYTUVgi6bEHHRwt9B6clvoQtMXkkStsfO8jBlv1+vhYqavb7Zeb7pCryseN8bGk1OVnC+GFZbCz3MXruaPYlSaaK5d6/RYt87wI4Wm0KY1DBkZOm9RzpVmnsmvEsbqXgw3I/5wx28NEXvqjOWl8eqz82eSJVEKcdqmrhG+fn8Mp6N2MH2diyP8B15+ewuNTNN6eaYZi/vr57UlGPhy5b0PHRQt/BydSXMNIythlmhcNMTB6R9WlsBnz0ua9B2v5Vk3JS3mhNdM1RAxyMGuDgkvHZvPWRhw92eDn/TCel5ab1mix+P9adYtGzwMAbUJzWDY6cMv00lhWvaNyCj6zaCXD26Ta2fmn+XkRg9lgnH+z0ce2UXJ5aUdvg/GSJVFazrgtGO1ixxU+eU+hXaGPqqCy27HeHN1ubGwWmwyU7PlroOziZ+BLGWsbf/0sVpbt89C+yhQW3KQ1RXlhdS60XuuUIp9yKqyZms7jUw6gB9qi0/XRCBa1NzGun5PDyOnf4mgA3Ts/jiuIclpV5EIG3Po5oqtFIrRdL5GeclcXOr3zh0gQ2gUkjs/hghxd/6NZFzEkjWSG6eFb89q8C5DrNyanaDW9ucuP2QXmF6arKdjS+ioncwDUEpp+VxaptPlZv92OImZT2zQty+Vto0ijZapaSbm4UmA6X7Phooe/gZOJLGOv+6V9kWoWRdU6a0hDF44drp2Qz5DQHjy2txu2Ha6eY7oRbZ+elnLYfudpQCiaNcDBnXA4f7PQxpI+Da6cIi0vd3DjdrIPTv8iWVitGS+RvmpHPnU+YTUGs8MZVW72cf6aD5Z+aYoyC6Wcl3nyPteIBhvQW9h5V9O1hcOhEkB6h1rNun7kqWLPdn9RVZYl8frY5aV49OZvFGzzkOaHWawr/KbfizU0u/AFzQnp1g4dcJ82OAtORZR0fLfQdnEx8CWPdP3MnZLN4g4vKakVdQKU9eURVrjwrG2eWwSXjs3nnYw8/+HoedQHSCpeMXW3sPhSgtNwbnuRumJZDXYC0G7dEck3onH6FNk65AtSEuiHXelW9yAO3JxhvbHkBgCsmOFha5mfvUcVp3YWDoc5ZeTlCZa1CkXo54DU7zFXMK+vd9Cs0GHJa/QTXv1A4WGW+cUXof5sRvyREU2nryDJN80hZ6EXEBmwCDkT2jBWRx4HvRnaZijnvHOAvmE3Fg8B5SilPvGM1TaO5X8JY98+kEU4OVgZ4d7MnrdBGi9jKlXPG5XBFcQ5vlXlY+H7qVvwr61ycCJX7feq9Gq49Pxe7Dc4eZI9KpErnmokiiV4rdXPTjPxwFq7HFwxvjkLiTdF4FrzF9gPBsNvq6Mn6GeDLY6rRSJp4Y60LmAlOH+z0hbt93XCBOcFNHZ3NaxvcKfU10HRN0rHo7wF2YAo2ACJSDDTMBql/3Q48B3xHKfWJiPQEGnZ01rQp8dw/G8t9Te4nGq9yZd8eBgU5Qo1HJY16sZKeDDEjXe6+LI+FJS5eWutiZD87q7f5UIDXF+RkbTDpNeMVCRs31MHrpebYxgy0h2rj13DNpBzW7hCOhXqUNLY5HM8Pb5HvhM8PBbh6krmHEBmL31iJBWu8sdFEj7xRzbOra6O6fS143/TFv7bBnbTbl6ZrIypyrZnoIJGBwELgYeDHSqnLQxb+e8CNQHk8i15E5gI3KqVuSnVAxcXFatOmTakerskQmYybjxSqTbu9rN3hS6lyZbza65v3+LlhmunUXvB+LU6HcFo3g817QzHmCapWWtz7t0pqIwzufCcYNkEpmDbGyeY9fn59ffcGUTeNRdPEi9Bx2iG2rA6Y92xlv6Yi8JF5BldPrl+d/PGtajbv8XPd+TnkOYXXNrg55VbkOeGG6fEn5I6WSKdpHiJSppQqjvdaqhb9Y8D9QEHEc3cDS5RSFSIJS5mOBJSIvAP0Bl5USv02zgDvBO4EGDRoUIpD0mSStqxcCQ2zZYNBKN3l5ZzBWeEInePVpvXsq1MpRejcv7CKWq/prw4EYWgfgz2HgxRkK7IcBh9/4ePQCfOajSVPRZIoDNNbR7gcQiTJRD5ZYxEgqizz4lI3dQFzgku2itHZrBqLpEIvIpcDR5RSZSIyM/Rcf+A6YGYK178AOA9wAStCs86KyIOUUvOB+WBa9OndgqY90tTWg1ZZYq9fhePMIxtr2G2pR5FY/U+DysyUvWZyHos31PLF4SB4ggjQryj1SKLYCpM2gbxswW6HymrVQOSBRnu1RsbD2w2YcEb8xiKTRjhTjiaKl31sNWvR2axdl1Qs+qnAvJAbJhvTR78N8AK7Q9Z8rojsVkoNjzn3K2C1UuoYgIgsA84FVqDpMEQmKCVzlTSVyMifuROyeWmti3c2ezjlUuRnC/6AoiBH+NYF6ZXO7V0AQQzGDzNXBldNzOKLwx665Qi1XsXcCdmNnh+5KdyzILr6V69uBkdOBsmJM5xkJRZKy72s3u7lmsk5zB6bHY4mmjQyq0FjkVSjiWJXRbGrAp3N2nVJKvRKqQeABwBCFv19kVE3oedr4og8wDvA/SKSC/iAGcCjzRyzphWJdFX0KzQY1MvGqm2m07spYp9K02hLyF9cY4rUY0trzKSnBPXRI2vIQH3IIgAi+PwqnCn71HtmwJfVzjBZBq7Vyq8oX5g00swtAMjJgsMng5zR18bnh+rFM1k0TewG7gc7PcwZlxPOXVi7w9ukxiLQcFUUb1Wgs1m7JhmPoxeReUCxUuohpVSViPwO+BAzN2SZUurNTL8n6I2nliBeU+6FJbWcN9zBmh3etIU+nabR3XMNsrOE66c1bsFbEToicO8V+eHuSC+tdZHnhKOnFL27CR6fCmfKCiRtnh0bNjlxpFk7p1eBcKxa4Q6F1keKfLKqlZbInzPEzpZ9dTgMOFip+NWLJ5g7IYfFG1wcr1ZNaiwCDVdFTVkVaDonaQm9UmolsDLO8/kRj5cASyJ+fg4zxLLF0BtPLYNlvV43JTfKz7todW24kUW618tk0+jSci8L3jfH0rPAoMajmDTCyR0XSyhCxyDLHuToqXrneaL68BaxdeKnjcli+5d1rN/pY+poJx9/4SO6b1Tj1410ewH0LIB/ntuNX754knFDHSzf7OFAZZAJw7LYss/H8Wp/kxqLAHFXRemuCjSdk06RGavLqLYMVoGzyGYVpp83iL0JHoBMNo22JvdASEAj/dEThmVRWWO6R+Yn6OwUj9gyAwXZsHaHj9ED7HxVGaB0l4/j1UF6FpjWdip9Wq26PDPGOLn7yRNUVpvPW/HwM8aYNee//5cqbEZqZYMTEW9VlO6qQNM56RRCr8uotgz9Cm2c1s0IJz3NGOMMN+WelqRZRaLrNacAW6x1PG6Igyy74PWrKH9091yDonwDpyNh2G+YBnV0RjrYWO6nX6G5gbtyq4fPDtbxnRm5LCgx7z0YVI0WSbOwyhZYkyRA30KJcns9u9JceaZStTMZuiaNJhGdQuh1GdWm09jehmV1njPYwWsb3Ly8rj6TtClWZ1MLsJmVMF3UehXdcoRrJuewoMTFJ3v9jB5o58vjQf66vIZbZpnVMJ9cXkNQwdWTU49SserolB/00yO3vp79BWOyePdjX1jk07G4rbIFFgOKjPBGbLWrjo3lXrx+K8Y+9SJsjaFr0mji0SmEXpdRbRrJ9jYswVhW5iGomm91NsXitMbo9iuunZzNkD5mnZfCPOGcIQ4+2GmWanh1nYvfv2mOXQS+lSBCJ/K6C96vJRAwx2P1VR07OItP9vrZWO5j4ogs1mwz9yms1nzpTHCRtX4AfnV9D37ydBWn3IofLzBrLCRLpNLBBZpM0CmEXi9Zm0YqexupWIjpiFJTE6keeaOaWWeblTBvmZXH396rYev+OuoCpk++e64RntyTFgqrDGCEsmX/cEcP9hwNhPuqrtnuI6jg1lk5vLbBTbXHbMt3w/TEoZ2J7j2y1s+MMU5WbfemXOtHBxdoMkmnEHrQS9amkIm9jZbuWXuwMsCi1bUYUm8dm83DFROH2zleHUxaYTO2qJkAxcMdbNrtZ9nHHq6elBue4M4Z4mDzHj8LSlxJywxEZrb2LTSiCqVNGlF/3uJS0+1ltzVe68dCBxdoMk2nEXpN+mRib6Ole9YuWl3LBaOdLP3QzSuhPYJ+PQwKsoUPd/tTroTp8cG1k7M5vbedx5bWsG2/n7GDHLzzsYfRAx0M7W02Fk+1jk5sZqs1wU0cEd2w5cbpeUmFPV6D80xMwNr1o7HQQt+FycTeRktEPEVOHpcV5/B6qZsrzsvhzTIPr25wE0yhEmbktbLsQo1SXPS1HOw2oWeBgb9OcfhkgGDQdPkdrAykVEcntuZNQY7RIMfAKpSWCokanC/d5ObqSbnh45rS4Uu7fjQWWui7MJnY22iJiKfIycMay5ub3Li8Ku0NYWvC6VtohMd55cQcnlpRyym3mYV6w7TcRv37sSUW8kLlcYryhVfXuzCM+loyh6qCKRdKg/grIqsb1+iBjiZNwNr1o4lFC30Xp7l7Gy0R8RSv41X3XCMqgzada3n9KlzU7JZZeXTLFbId4PFDRZIs1MgSCz0LDGaelcV7n3oxRDGol50vjtSxeIObCcOyeH+LB8MgpUJpkbH7VTUBoH6itBqcpzoBx7ppDmbA9aPpXGih1zSLloh4yuTkYVXC/GCHlylnml2ZjlcHMSR5PffYOvFvfeThoq/lMKSPg7+8U8OW/X4uHufk7Y+84czWxjo9Ra4MehYY3DY7jzc2unllvRvDkKhG7Kk2OG8J14+m86GFXtNsMh3xlMnJwzrnlXUu3v7YDLspzJNGM1tjy/0GglC6y0f3XGH3oTqG97VT61UoBR+Wm52ukrmUrGsahnDvFXnYDGFhSS3nDHbw0Re+8KogE26a5rp+NJ0PLfSadkmmO15lognKe594WFhSy6yznRTmCSKSci0Z65q/W1LNmf0d2CNKNlw7JZenVtQ2qxF7JFYjdp1XorHQQq/pkLRk+GC8Jihvf+ymxg2XjM/i1fVuggoK80ip5k3kNSP3Hyy/eWG+kbKrJpaEm+FNvJ6mc6KFXtPhaOnwwfhNUGpRwNsfe8ONTRK9V7yOXNY1I/cfAgFFUb7RqFsl2YSmy39oUkELvabDkanwwUQiGr8JisH105K7P2JLE6/a7uXV9W5GDbCHrznvvOzwpnBjk0YqE5ou/6FJhZSFXkRswCbgQGQrQRF5HPhuZPOROOcOArYDv1JK/V8zxqvJIB01e7K1Sjc0RTxjSxNb/y8udXPr7Lyoaybz7ac6oenyH5pkpGPR3wPswGwODoCIFAOJ2/XU8zvgrfSGpmlJOnL2ZGuUbmiqeMaWJoZQHf917rSvqfssaDJFSkIvIgOBy4CHgR+HnrMB/wvcCFzdyLlXAXuA2maOVZNBMpk92dorg/ZQuiGeH/7G6XkNShOD+XNTOnLpPguaTJGqRf8YcD9QEPHc3cASpVSFSPxOPiKSD/wMuBi4L9HFReRO4E6AQYMGpTgkTXPIlLXYFiuDti7dkMgPD/FLE7+63s2Ms1LLao1t/KI3WjWZIKnQi8jlwBGlVJmIzAw91x+4DpiZ5PRfAY8qpWoSTQYASqn5wHyA4uLihp2XNRknU9ZiW9VVacvSDY354Z/4XlH4cbLSxKk2ftEbrZrmkopFPxWYJyJzgWxMH/02wAvsDgl4rojsVkoNjzl3EvANEfkt0AMIiohHKfWHTN2ApmlkylrsqH7k5ohoY354SK00MWSu8YtGk4ykQq+UegB4ACBk0d8XGXUTer4mjsijlJoWccyvgBot8u2DTFmLzVkZtHXUT1NFNFN++I46SWo6HhmPoxeReUCxUuqhTF9bk1kyYS02p+F3pnz7rT1hpOuHT0Rz3WdtPVFqOg5pCb1SaiWwMs7z+RGPlwBL4hzzq7RHp2n3NHVlkCnffrJ2fi1BU1sExtIc91lHDo/VtD46M1bTbJqyMshU0lMq7fxaglT98I3RHPeZbi6iSQct9Jo2IVNJT4EgzB6b3ax2fm1JU91n2r+vSQejrQeg6ZpYboudB/zUBRQ7D/hZWFKbtDtTJKYVbLYItAi38+vkSUXWRBmJTqbSJEJb9Jo2IVNJT+OGOqL83Km28+vo6GQqTTpooe/itGXkRipui2SZo6+Xupk4IotFq2upqAombefXWdDJVJp00ELfhWnvkRupZo4uK/Nw6EQwaTu/zoZOptKkihb6Lkx7j9zQmaMaTWbQm7FdmPYeudHex6fRdBS00Hdh2nvkRnsfn0bTUdCumy5Me4/caO/jayq6dIGmtdFC34Vp75Eb7X18TaG9b4BrOida6Ls47X0zs72PL13a+wa4pnOiffQaTSuiN5g1bYEWeo2mFdEbzJq2QAu9RtOKZKLGj0aTLtpHr9G0Ip1xg1nT/klZ6EXEBmwCDkS2EhSRx4HvRjYfiXjtYuC/gSzAB/xUKfV+s0et0XRgOtsGs6b9k45Ffw+wA7M5OAAiUgwUNnLOMeAKpdRBETkbeAcY0JSBajQajaZppOSjF5GBwGXAkxHP2YD/Be5PdJ5S6mOl1MHQj9uAHBHRpoxGo9G0Iqluxj6GKeiRbXvuBpYopSpSvMa1wEdKKW/sCyJyp4hsEpFNR48eTfFyGo1Go0mFpEIvIpcDR5RSZRHP9QeuA36fypuIyFnA/wDfi/e6Umq+UqpYKVXcu3fvlAau0Wg0mtRIxUc/FZgnInOBbEwf/TbAC+wWEYBcEdmtlBoee3LI7bMYuFkp9XnGRq7RaDSalBClVOoHi8wE7ouMugk9X5Mg6qYHsAr4tVLqtRTf4yiwL/RjL8wN3c5MV7hH6Br32RXuEbrGfXbEexyslIrrEsl4HL2IzAOKlVIPYfrxhwMPichDoUPmKKWOJDo/cqAiskkpVZzpMbYnusI9Qte4z65wj9A17rOz3WNaQq+UWgmsjPN8fsTjJcCS0OPfAL9p1gg1Go1G0yx0CQSNRqPp5LR3oZ/f1gNoBbrCPULXuM+ucI/QNe6zU91jWpuxGo1Go+l4tHeLXqPRaDTNRAu9RqPRdHLaVOhF5F9EZJuIbBWRF0QkW0QWiMgeEdkc+jcudKyIyOMisltEPhWRc9ty7OkgIveE7nGbiNwbeq5IRJaLSHno/8LQ8x3yPhPc469E5EDE73JuxPEPhO7xMxG5pM0GngQReUpEjojI1ojn0v7dicgtoePLReSWtriXRKR5jzNF5GTE7/ShiHMuDf0+d4vIz9viXhojwX1eF/qbDYaKNEYeH/dvtL3fZ1yUUm3yD7OK5R4gJ/Tz34FbgQXAN+IcPxd4CxBgMlDaVmNP8z7PBrYCuZjhrO9h5hb8Fvh56JifA//TUe+zkXv8FWaCXezxY4BPACcwFPgcsLX1fSS4t+nAucDWiOfS+t0BRcAXof8LQ48L2/remniPM4F/xLmGLfR7HIZZlvwTYExb31sK9zkaOBMzbLw42d9oR7jPeP/a2nVjx6xoaccUiYONHHsl8Iwy2QD0EJF+rTHIZjIa8wvvUkrVYWYKX4N5PwtDxywErgo97oj3megeE3El8KJSyquU2gPsBia2wjjTRim1GqiMeTrd390lwHKlVKVSqgpYDlza4oNPkTTvMRETgd1KqS+UUj7gxdA12g3x7lMptUMp9VmcwxP9jbb7+4xHmwm9UuoA8H/AfqACOKmUejf08sOhpe+jUl/WeADwZcQlvqJj1LbfCkwTkZ4ikotp9Z0O9FH1lT8PAX1CjzvifSa6R4C7Q7/Lp6zlPx3zHiNJ93fXEe830T0CTBGRT0TkLTELFkLHvMfG6Ey/y7YT+tCX/krMZVF/IE9EbgIeAEYB52EudX/WVmPMBEqpHZiVO98F3gY2A4GYYxTQYeNcG7nHJ4AzgHGYk/kjbTPClqOj/+5SIeYeP8KsqfI1zOq1r7fVuDSp05aum4uAPUqpo0opP/AacL5SqiK09PUCT1O/pD9AvZUIMDD0XLtHKfU3pdQEpdR0oArYBRy2XDKh/636Px3yPuPdo1LqsFIqoJQKAn+lE/wuQ6T7u+uI9xv3HpVSp5RSNaHHywCHiPSiY95jY3Sm32WbCv1+YLKI5IqIABcCOyL+uATTL2jtkC8Bbg5FNkzGdPWk2vSkTRGR00L/D8L0XS/CvB8r+uIW4I3Q4w55n/HuMWZv4Wqif5fXi4hTRIYCI4CNrTneZpLu7+4dYI6IFIZWsnNCz7Vn4t6jiPQNfTcRkYmYGnIc+BAYISJDRSQLuD50jY5Kor/RjnmfbbkTDPwa2IkpAM9i7nC/D2wJPfcckB86VoA/Yu54byFih7y9/wPWANsxd+gvDD3XE1gBlGNGqRR15PtMcI/Phu7hU8wvQ7+I4x8M3eNnwNfbevyN3NcLmG4nP6Y/9vam/O6A72Ju6O0Gbmvr+2rGPd6N2Y/iE2AD5ircus5czNXq58CDbX1fKd7n1aHHXuAw8E6yv9H2fp/x/ukSCBqNRtPJaevwSo1Go9G0MFroNRqNppOjhV6j0Wg6OVroNRqNppOjhV6j0Wg6OVroNRqNppOjhV6j0Wg6Of8fZLAaXvEERYAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "## plot the simulated values, \n",
    "## as well as their numerical average (for any dataset size > 6 simulations), \n",
    "## and finally the expected value.\n",
    "\n",
    "# plot raw simulations\n",
    "plt.plot(history[:, 0],history[:, 1], 'o', fillstyle='none', color = 'cornflowerblue')\n",
    "\n",
    "min_count = int(min(history[:, 0]))\n",
    "max_count = int(max(history[:, 0]))\n",
    "\n",
    "count_range = range(int(min(history[:, 0])), int(max(history[:, 0])))\n",
    "expected_values = np.zeros(shape=(len(count_range), 2), dtype=float)\n",
    "expected_values[:, 0] = list(count_range)\n",
    "\n",
    "data = np.genfromtxt(data_path, delimiter=',', names=True)\n",
    "age = list(data[:]['age'])   \n",
    "true_mean_age = statistics.mean(age)\n",
    "\n",
    "counter = 0\n",
    "for count in count_range:\n",
    "    flag = history[:,0] == count\n",
    "    subhistory = history[flag,1]\n",
    "    if len(subhistory)>6:\n",
    "        mean_sim = np.mean(subhistory)\n",
    "        # plot average simulation by dp release of dataset size\n",
    "        plt.plot(count, mean_sim, 'o', color = 'indigo')\n",
    "    if count<=1000:\n",
    "        expected_values[counter, 1] = true_mean_age\n",
    "    else:\n",
    "        expected_values[counter, 1] = ((true_mean_age* 1000) + ((count-1000)*50))/count\n",
    "    counter += 1\n",
    "\n",
    "# plot the expected value by dp release of dataset size\n",
    "plt.plot(expected_values[:, 0], expected_values[:, 1], linestyle='--', color = 'tomato')\n",
    "#plt.xlabel('DP Release of Age')\n",
    "#plt.ylabel('n')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the banding is an artifact of the snapping mechanism. This provides protection against floating-point attacks. We observe a similar trend across a larger range of sample sizes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-15T16:41:48.688821Z",
     "iopub.status.busy": "2021-04-15T16:41:48.688134Z",
     "iopub.status.idle": "2021-04-15T16:41:55.140443Z",
     "shell.execute_reply": "2021-04-15T16:41:55.139910Z"
    },
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "\n",
    "# initialize data\n",
    "n_range = range(100, 2001, 200)\n",
    "n_sims = 50\n",
    "ns = []\n",
    "releases = []\n",
    "\n",
    "with sn.Analysis(dynamic=True) as analysis:\n",
    "    data = sn.Dataset(path = data_path, column_names = var_names)\n",
    "    age_dt = sn.to_float(data['age'])\n",
    "\n",
    "    for n in n_range:\n",
    "        for index in range(n_sims):\n",
    "            # get mean of age\n",
    "            ns.append(n)\n",
    "            releases.append(sn.dp_mean(data = age_dt,\n",
    "                              privacy_usage = {'epsilon': 1},\n",
    "                              data_lower = 0.,\n",
    "                              data_upper = 100.,\n",
    "                              data_rows = n))\n",
    "analysis.release()\n",
    "\n",
    "# get released values\n",
    "values = [release.value for release in releases]\n",
    "df = pd.DataFrame({'n': ns, 'release': values})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-15T16:41:55.155238Z",
     "iopub.status.busy": "2021-04-15T16:41:55.145060Z",
     "iopub.status.idle": "2021-04-15T16:41:55.437767Z",
     "shell.execute_reply": "2021-04-15T16:41:55.437274Z"
    },
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAeC0lEQVR4nO3df3Rc5X3n8fdXsvxD/oHBIwx4ACXYQNKWEqKyyVIoCZVTAYcUGrK0SXa2QEnaBNN6WwqbnACnyTlNu9lN7dMlJQnsNCEJtAlbl1ggbVIv7DnZEBmDYmPAA8hBBtsag0HGBo+t7/5xr/BYaOSZke69I93P6xwd3bk/5vnqauY7zzz3uc9j7o6IiKRHU9IBiIhIvJT4RURSRolfRCRllPhFRFJGiV9EJGVmJR1ANTKZjLe3tycdhojItLJx48aiu7eNXT8tEn97ezt9fX1JhyEiMq2Y2fbx1qupR0QkZZT4RURSRolfRCRllPhFRFImFYm/WCxy4403smfPnqRDERFJXCoSfz6fp7+/n3w+n3QoIiKJm/GJv1gs0t3djbvT3d2tWr+IpN6MT/z5fJ7RoadHRkZU6xeR1Jvxib+3t5dSqQRAqVSip6cn4YhERJI14xN/Z2cnLS0tALS0tLBy5cqEIxIRSdaMT/y5XA4zA6CpqYlcLpdwRCIiyZrxiT+TydDV1YWZ0dXVxZIlS5IOSUQkUdNikLbJyuVyDAwMqLYvIkJKEn8mk2Ht2rVJhyEi0hBmfFOPiIgcTYlfRCRllPhFRFJGiV9EJGWU+EVEUkaJX0QkZZT4RURSRolfRCRllPhFRFIm0jt3zWwAGAYOA4fcvcPMTgDuA9qBAeDj7v5qlHGIiMgRcdT4P+Tu57p7R/j4FuDH7r4C+HH4WEREYpLEWD0fBS4Ol/PABuAvJ/uka9asoVAojLttcHAQgGw2+45ty5cvZ9WqVZMtXkRk2oi6xu9Aj5ltNLMbwnVL3f3lcHknsHS8A83sBjPrM7O+oaGhSQVx4MABDhw4MKnnEBGZKWx0PtpIntxsmbvvMLMTgV7gRmCduy8u2+dVdz9+oufp6Ojwvr6+uuMYrdGvWbOm7ucQEZluzGxjWTP72yKt8bv7jvD3buAB4Hxgl5mdHAZ1MrA7yhhERORokSV+M5tvZgtHl4GVwGZgHTA6I0oO+JeoYhARkXeK8uLuUuCBcL7bWcB33f0hM/s5cL+ZXQdsBz4eYQwiIjJGZInf3Z8Hfn2c9XuAS6IqV0REJqY7d0VEUkaJX0QkZZT4RURSRolfRCRllPhFRFJGiV9EJGWU+EVEUkaJX0QkZZT4RURSRolfRCRllPhFRFJGiV9EJGWU+EVEUkaJX0QkZZT4RURSRolfRCRllPhFRFJGiV9EJGWU+EVEUkaJX0QkZSKbbF1ERCpbs2YNhUJh3G2Dg4MAZLPZcbcvX76cVatW1V22Er+ISIM5cOBApM+vxC8ikoCJauyj29asWRNJ2WrjFxFJGSV+EZGUUeIXEUkZJX4RkZSJPPGbWbOZbTKzB8PHl5jZ42b2hJn9XzNbHnUMIiJyRBw1/puArWWP7wQ+4e7nAt8FvhBDDCIiEoo08ZtZFrgM+GbZagcWhcvHAS9FGYOIiBwt6n78XwNuBhaWrbseWG9mB4DXgQ+Md6CZ3QDcAHDaaadFG6WISIpEVuM3s8uB3e6+ccymPwMudfcscA/w38Y73t3vcvcOd+9oa2uLKkwRkdSJssZ/AXCFmV0KzAUWmdmPgLPd/WfhPvcBD0UYg4iIjBFZ4nf3W4FbAczsYuDPgd8FdprZme7+LNDJ0Rd+RUQiV2mAtKgHR2sUsY7V4+6HzOyPgB+Y2QjwKnBtnDGISHKSHJGyGlEPjtYoYkn87r4B2BAuPwA8EEe5IjJ9xJl0K32ARD04WqPQ6JwiEpskR6SUI5T4RUQiMlHT1kS2bdsGTPxBOZ5qm8OU+GNSLBa54447uP3221myZEnS4YhIDAqFApue2srhtqU1HdcU9rTvG3ql6mOah3ZVva8Sf0zy+Tz9/f3k83lWr16ddDgiEpPDbUt54/c+FXk583/w7ar31eicMSgWi3R3d+PudHd3s2fPnqRDEpEUU+KPQT6fx90BGBkZIZ/PJxyRiKSZEn8Ment7KZVKAJRKJXp6ehKOSETSTG38Mejs7GT9+vWUSiVaWlpYuXJl0iGJzHj19KiptzcNjN+jZnBwkObXh2tqf69X89AuBt/aX9W+SvwxyOVydHd3A9DU1EQul0s4IpGZr1Ao8NRTm8i0eQ1HGQC7hx6vqazikNW0f9KU+GOQyWTo6upi3bp1dHV1qTunSEwybc5VVx2MvJwf/nD2uOuz2Sw7h16JrVdPtu2EqvZV4o9JLpdjYGBAtX0RSZwSf0wymQxr165NOgxJsbSNSDk4OMjrr1vF2vhUKg4ZB98ajLycqaLEL5JyUQyOFueF1en6wZQkJX6RlIhzRMpCocDmzZtZsGBB1ceMdnkeGBio+ph9+/ZV3JbNZtk9tDu2Nv4T28b/xtQ8tKvmXj1Ne18FYGTx8VUf0zy0C9TGLyJJWrBgAeedd16kZTz+eG29b+K2fPnyuo7btje4u39FlYkcgLYTqi5PiV9EJCL1NkFFPUS17twVEUkZ1finUKNPKyeSNsWh2nr1vLY3uBHruMW13PQVlHNiW02HJEqJPyZpmctTpFHU077+2t6gZ9GJbStqOu7Etvrb85OgxD+FNK2cSGBwcJDh4eHIL74ODw+//W16rHq+Qaflfao2fhGRlFGNX0SmXDab5dChQ7F056x03UwqU41fRCRllPhFRFKmqqYeMzsTuBNY6u6/ambnAFe4+5cijU5EatIIk49I46u2jf8bwF8A/wDg7v1m9l1AiV+kgRQKBZ7ZvJVTF55U9TEth4Iv/vu3v1pTWS8O75xw+759+2rq1bN/fzB7VGtra9XHTDRWT6Ob6EP6WB/Gk/3ArTbxt7r7Y2ZHzTJzqO5SRSQypy48if98/h9GXs5XH7un4rZ6+rSPJrv29vaajptO/eerNW/evEifv9rEXzSzMwAHMLOPAS9Xc6CZNQN9wA53v9yCT48vAVcDh4E73X1md5oVSRn1oT+2JJvIqk38nwXuAs42sx3AC8Anqzz2JmArsCh8/J+AU4Gz3X3EzE6sPlwRmcjg4CBvDA9PWBufKi8O72T+4BuRlyNTr6rE7+7PA79tZvOBJncfruY4M8sClwFfBlaHq/8Y+AN3Hwmfe3e1wdZz4QrSN8FDvWMGTde/V6RWld4jUbetN4pqe/XcBNwDDAPfMLPzgFvcvecYh34NuBlYWLbuDOA/mNmVwBCwyt23jVPmDcANAKeddhoQXLja9IunGGmtYYxqwA4GAy5tfG7ii1Hlmva/UlMZ04XGDJrZstks+w+/Glsbf2u2+olCpoOo29YbRbVNPde6+9+Z2UeAJcCngG8DFRO/mV0O7Hb3jWZ2cdmmOcCb7t5hZlcBdwMXjj3e3e8iaF6io6Pj7aHyRlpP4M33Xl5l2PWb+9SDkZcRFY0ZJDKxmVBrn4xqE/9od55LgX909y02povPOC4ArjCzS4G5wCIz+w4wCPww3OcBgm8SIiISk2rv3N1oZj0Eif9hM1sIjEx0gLvf6u5Zd28HrgF+4u6fBP4X8KFwt98Cnq0ncBERqU+1Nf7rgHOB5919v5ktAeptRPxr4F4z+zNgH3B9nc8j01SxWOSOO+7g9ttvZ8mSJUmHE7m4J+h5cXhnTb16dofXs06s8drZi8M7OYuZ1cafFtX26hkxsxeAM81sbq2FuPsGYEO4vJegp4+kVD6fp7+/n3w+z+rVq499wAw21Rfb67mZqbStCEDr6bUl8bM4fkbePJUG1fbquZ6gP34WeAL4APBT4MORRSYzUrFYZP369bg7P/rRj8jlcjO+1h/nxXbdOCXVqLaN/ybgN4Dt7v4h4H3A3qiCkpkrn89TKpUAKJVK5PP5hCMSSZ9qE/+b7v4mgJnNcfengbOiC0tmqocffvioxw899FBCkYikV7UXdwfNbDFBj5xeM3sV2B5VUDJzzZo1a8LH01mcQyLPlDtIJRnVXty9Mly83cz+DTgOUFVNajZ2GN3pPKzuWIVCgc1PPsnC2dV/mB06dBiA7Vu3VH3M8EENjCuTU/Ur1Mx+E1jh7veYWRuwjGCwttRRza5+7e3tDAwMHPU4SpX+V1F0owRYOHsW5y+NtovjY7tqGze/kSQ5Br0cUW2vntuADoJ2/XuAFuA7BHfnpk6hUODZzY9z2oLDVR8zuxRcTnlz4OdVH/PLfc01x9YoKr3BZ8+e/Y7HY9/McbzBoxizaHBwkOGDhyJPzMMHD739wVWLRh+YLC3j5DSCamv8VxL05HkcwN1fCu/eTa3TFhzmCx3RNlN8qW9BpM+fhNbWVswMd2fOnDk1zbZUj0oJS10Yj4gz4arG3hiqTfwH3d3NbHQilvkRxiQzwERv8Ouvv55CocCdd945o24AymazHB5+LZamnkpNVBNR0pVR1Sb++83sH4DFZvZHwLUE8/CK1Ky1tZVzzjlnRiX9UbU29ewPL+62zqq+WU8Xd2Wyqu3V81/NrBN4naCd/4vu3htpZCLTzGTmmT19xYrIyxIZVXWvnjDRJ5rsBwcHadr/Wixj5Tft38Pg4Pg1q2B6u+bI2+C3Dzczv8JFPM1GdkScvayg8rnQcAkyXUyY+M1smHCC9bGbAHf3ReNskxgUCgU2bdkEi2s8MBxMe9OOTdUfs7fGMmJWKBTY8outLG6tfvrmkYPBdBI7nttTU1l791c9U6hIw5ow8bt7Q/XcyWaz7HprVmwzcGWzJ1WM481DL8fSq2fuRBfxFsPIxRNOizAlmjZUO7JHMurp2rhgbv0XYOspT6SR1HsDVwZY6O6pvIFLjtDNbNNf2uZHEDD38VpyxuxUdgOXu59pZqcA/+TusdzAdcLp7/HO/3I3hUKB4QMHGZkf/Yuz6Y09LJw3e9yLaIVCgZE3hzl9YfU3cNVj+3AzTXMXVoxh31v78MXH/v9Nlu01FsxZUDGO/fv2MXuc4yophb9bajjmINC6oHIMbx0o1dTUU6+9+3czZ17LlF1cHf3QTPJi7eDgIMVikUwmU1c3UWlc93/m3290946x63UDl0zK7O0/Zc7ru2tK/PU4CPiiE2Ea9mbZsWNHxTuFR9dX+tY0b948li1bFllspVKJPa8E1zleeeUVTjppKbNm1fKRLNPRtLiB691t87nv0x9k1arvsfGFnfG08b/wMO8/4yTWfPpT79i2atX3eHP7z6Nv49+8gLntv1Exhk07NjGyIoY2/i1NvG/Z+8aN46qH/5bh/UVOjjiGl4GFrSPc9+kPvmPbqlXfY8uW2i7u7nsz6Gtfa1v/nP27+ZVfe8+456KSCadefOkYYwadupxVn/5YTTHW4qtf/SqvPb2eUqlES0sLF5xxGas/m+5Z0WaS+z8z/vpjJn4zM+BB3cDVWAYHB+G1mC687oVBH/+CZjabZW+xyHVYpCF8C2fxBAOq1WrbtmCe2WVn1NZsuIwlNZfXyNclent7j5oYp6enJ/XTYabBMRN/WNO/GliNbuCScewkSMzVGu1AWUvK3UnlnqvqP1+/zs5O1q8/UuNfuXJl0iFJDKpt6nkc2OvufxFlMFK9bDbLkA3F1p0zu2zqattDYa+exTXcrbq4zrJkYrlcju7ubgCamprI5XIJRyRxqDbx/zvgE2a2HXhjdKW7nxNJVDJtqLY9vWUyGbq6uli3bh1dXV3qzpkS1Sb+j0QahYgkJpfLMTAwoNp+ilQ7SJvm1x3jl/tqG6tn1/7gIuzS1uqbZn65r5kza45MpDaZTIa1a9cmHYbEaObMdB2jetqaD4bt2nPbq2/XPrPOshqBptgTaVxK/HVomHbtvXV05xy99aCWgUX3EsywPEXinPGp0acbFEmCEv80Ve83gdGEt2JZDeO/L6u9vImS5ujYMLfddltiFxM1v6ukWeSJ38yagT5gh7tfXrZ+DXCtu8+8iWVjUG9ttBF61OTzefr7+8nn85HfLKRau8g7xTHe7k3A1vIVZtYBRDsxqTSkYrFId3c37k53dzd79tQ2Hr6ITF6kid/MssBlwDfL1jUDfwvcHGXZ0pjy+TwjI0HPpsOHD5PP5xOOSCR9oq7xf40gwZf3YfwcsM7dX57oQDO7wcz6zKxvaGgowhAlTr29vRw6FExpeejQIXp6ehKOSCR9Ikv8ZnY5sNvdN5atOwW4Gjhmp2F3v8vdO9y9o62tLaowJWYXXnjhUY8vuuiihCIRSa8oL+5eAFxhZpcCc4FFwBbgLaAQDPpJq5kV3H16dlYXEZmGIqvxu/ut7p5193bgGuAn7n68u5/k7u3h+v1K+uny6KOPHvX4kUceSSgSkfRq7Fm0Zcbp7Oxk1qzgi+asWbM0DLBIAmK5gcvdNwAbxlk/o/rwa5iCYysfBri5uVkDg4kkQDX+mMybN093i3JkGGAz0zDAIgnRkA1TKA019qmgYYBFkqXEL7HTMMAiyVJTj4hIyijxi4ikjBK/iEjKKPGLiKTMtLu427T/FeY+9WBNx9ibrwPgcxfVVA6cVFM5jaLe+wnSci+BSNpNq8Rf/6xTwwCsOKOWRH7StJ3vdiK6l0BEzN2TjuGYOjo6vK+vr+7jG2HWKRGRuJnZRnfvGLtebfwiIimjxC8ikjJK/CIiKaPELyKSMkr8IiIpo8QvIpIySvwiIimjxC8ikjJK/CIiKaPELyKSMkr8MSkWi9x4443s2bMn6VBEJOWU+GOSz+fp7+8nn88nHYqIpJwSfwyKxSLd3d24O93d3ar1i0iilPhjkM/nGR0FdWRkRLV+EUmUEn8Ment7KZVKAJRKJXp6ehKOSETSTIk/Bp2dnbS0tADQ0tLCypUrE45IRNJMiT8GuVwOMwPAzMjlcglHJCJpFnniN7NmM9tkZg+Gj+81s2fMbLOZ3W1mLVHHkLRMJsMpp5wCwCmnnMKSJUsSjkhE0iyOGv9NwNayx/cCZwO/BswDro8hhkQVi0V27NgBwEsvvaRePSKSqEgTv5llgcuAb46uc/f1HgIeA7JRxtAIynvxuLt69YhIoqKu8X8NuBkYGbshbOL5FPDQeAea2Q1m1mdmfUNDQ5EGGTX16hGRRhJZ4jezy4Hd7r6xwi7/A3jE3R8db6O73+XuHe7e0dbWFlWYsVCvHhFpJFHW+C8ArjCzAeD7wIfN7DsAZnYb0AasjrD8hlHeq6epqUm9ekQkUZElfne/1d2z7t4OXAP8xN0/aWbXAx8Bft/d39EENBNlMhm6urowM7q6utSrR0QSlUQ//q8DS4GfmtkTZvbFBGKIXS6X45xzzlFtX0QSNyuOQtx9A7AhXI6lzEaTyWRYu3Zt0mGIiOjOXRGRtFHiFxFJGSV+EZGUUeIXEUkZJX4RkZRR4hcRSRklfhGRlFHiFxFJGSV+EZGUUeIXEUkZJX4RkZRR4hcRSRklfhGRlFHiFxFJGSX+mBSLRW688Ub27NmTdCgiknJK/DHJ5/P09/eTz+eTDkVEUk6JPwbFYpHu7m7cne7ubtX6RSRRSvwxyOfzuDsAIyMjqvWLSKKU+GPQ29tLqVQCoFQq0dPTk3BEIpJmSvwx6OzspKWlBYCWlhZWrlyZcEQikmZK/DHI5XKYGQBNTU3kcrmEIxKRNFPij0Emk6Grqwszo6uriyVLliQdkoik2KykA0iLXC7HwMCAavsikjgl/phkMhnWrl2bdBgiImrqERFJGyV+EZGUUeIXEUmZyBO/mTWb2SYzezB8/C4z+5mZFczsPjObHXUMIiJyRBw1/puArWWPvwL8d3dfDrwKXBdDDCIiEoq0V4+ZZYHLgC8Dqy24i+nDwB+Eu+SB24E7J1vWmjVrKBQK427btm0bAKtWrXrHtuXLl4+7XkRkpoq6xv814GZgJHy8BNjr7ofCx4PAsvEONLMbzKzPzPqGhoYmFcS8efOYN2/epJ5DRGSmiKzGb2aXA7vdfaOZXVzr8e5+F3AXQEdHhx9rf9XaRUSqE2VTzwXAFWZ2KTAXWAT8HbDYzGaFtf4ssCPCGEREZIzImnrc/VZ3z7p7O3AN8BN3/wTwb8DHwt1ywL9EFYOIiLxTEv34/5LgQm+BoM3/WwnEICKSWrGM1ePuG4AN4fLzwPlxlCsiIu+kO3dFRFJGiV9EJGWU+EVEUkaJX0QkZcz9mPdGJc7MhoDtk3yaDFCcgnCmewzQGHEohiMaIY5GiAEaI45GiAGmJo7T3b1t7Mppkfingpn1uXtH2mNolDgUQ2PF0QgxNEocjRBD1HGoqUdEJGWU+EVEUiZNif+upAOgMWKAxohDMRzRCHE0QgzQGHE0QgwQYRypaeMXEZFAmmr8IiKCEr+ISOrMiMRvZneb2W4z21y27gQz6zWzbeHv48P1ZmZrwsne+83svCmKYa6ZPWZmT5rZFjO7I1w/7uTyZjYnfFwIt7dPRRzhcw+Y2S/M7Akz6wvXxXY+zOyssOzRn9fN7E/j/p+Ez32TmW0O/yd/Gq6LNI4Kr8erwxhGzKxjzP63hmU+Y2YfKVv/O+G6gpndMkVx/FX4tz1hZj1mdsqx/nYzy4XnapuZ5aYghtvNbEfZ6+PShM7FfWUxDJjZE1HGUSGGXzezn4bv1381s0VRnwsA3H3a/wAXAecBm8vW/Q1wS7h8C/CVcPlSoBsw4APAz6YoBgMWhMstwM/C578fuCZc/3Xgj8PlPwG+Hi5fA9w3hedjAMiMWRfr+SgrtxnYCZyewP/kV4HNQCvBSLT/G1gedRwVXo/vAc4iGKW2o2z9e4EngTnAu4DnwnPWHC6/G5gd7vPeKYhjUdnyqrLX4Lh/O3AC8Hz4+/hw+fhJxnA78Ofj7BvruRiz/avAF6OMo8K5+DnwW+HytcBfRX0u3H1m1Pjd/RHglTGrP0owmTvh798tW/+PHvh/BDOCnTwFMbi77wsftoQ/TjC5/D9XiGM0vn8GLjEzm2wcE4j1fJS5BHjO3bcnEMN7CBLYfg9mfPs/wFVRxzHe69Hdt7r7M+Ps/lHg++7+lru/ABQIhi0/Hyi4+/PufhD4frjvZON4vezhfILX6Ggc4/3tHwF63f0Vd38V6AV+ZzIxTCDWczEqfN99HPhelHFUiOFM4JFwuRf4vShjGDUjEn8FS9395XB5J7A0XF4GvFi2X8UJ32tlZs3h18XdBP/E56g8ufzbcYTbXyOYmGYqONBjZhvN7IZwXeznI3QNR95QccewGbjQzJaYWStBrfbUBOKYSKUyo3ydftnMXgQ+AXwxoTg+FzYp3T3a1JZADKMuBHa5+7YE4tjCkcR9NcHrM/IYZnLif5sH350i77fq7ofd/VyCuYTPB86OuswKftPdzwO6gM+a2UXlG+M6HxZcz7gC+Kex2+KIwd23Al8BeoCHgCeAw3HH0Wjc/fPufipwL/C5BEK4EzgDOBd4maCZJUm/z5HKSdyuBf7EzDYCC4GDcRQ6kxP/rtGv6eHv3eH6HRz5VIUIJnx3970Ecwt/kHBy+XHKejuOcPtxwJ4pKn9H+Hs38ADBh1AS56MLeNzdd4WPY4/B3b/l7u9394uAV4Fnk4hjApXKjCOWeznStBBbHO6+K6wkjQDf4MiMfLGfi/C9dxVwX9nqOM/F0+6+0t3fT/Dh81wcMczkxL+OYDJ3OHpS93XAfwx7MXwAeK3sa3/dzKzNzBaHy/OATmArlSeXL4/vYwST0U+65mlm881s4egysJKgySPW8xEaW5OKPQYzOzH8fRrBG/y7ScQxgXXANRb08noXsAJ4jOCi3woLeoXNJmgyWzfZwsxsRdnDjwJPl8Ux3t/+MLDSzI4Pm2RWhusmE0P5dZMrCV6fozHEdi5Cvw087e6DZetii6Ps9dkEfIGgA0j0MdR6NbgRfwiSy8tAiaDN6zqC9vIfA9sIenOcEO5rwN8TfLL+grIeFpOM4RxgE9BP8EIe7SHw7vAfViBo8pgTrp8bPi6E2989RXG8m+BK/5ME7YefD9fHfT7mE3yDOa5sXawxhM/9KPBUeD4uiSOOCq/HK8Plt4BdwMNl+38+LPMZoKts/aUE31CeG/0/TkEcPwhfn/3AvwLLjvW3EzRHFMKfP5yCGL4dltFPkLROTuJchOv/J/CZcfaf8jgqnIubwud7FvhrwtEUojwX7q4hG0RE0mYmN/WIiMg4lPhFRFJGiV9EJGWU+EVEUkaJX0QkZZT4RURSRolfRCRllPhF6mBm7Wa21cy+YcE4+z3hHdsiDU+JX6R+K4C/d/dfAfZyZNwbkYamxC9Svxfc/YlweSPQnlwoItVT4hep31tly4cJZvkSaXhK/CIiKaPELyKSMhqdU0QkZVTjFxFJGSV+EZGUUeIXEUkZJX4RkZRR4hcRSRklfhGRlFHiFxFJmf8PwYThD1k2nr4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# get true mean\n",
    "data = np.genfromtxt(data_path, delimiter=',', names=True)\n",
    "age = list(data[:]['age'])   \n",
    "true_mean_age = statistics.mean(age)\n",
    "\n",
    "# plot distribution of releases by n\n",
    "plot = sns.boxplot(x = 'n', y = 'release', data = df)\n",
    "plot.axhline(true_mean_age)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SmartNoise `resize` vs. other approaches\n",
    "The standard formula for the mean of a variable is:\n",
    "$\\bar{x} = \\frac{\\sum{x}}{n}$\n",
    "\n",
    "The conventional, and simpler, approach in the differential privacy literature, is to: \n",
    "\n",
    "1. compute a DP sum of the variable for the numerator\n",
    "2. compute a DP count of the dataset rows for the denominator\n",
    "3. take their ratio\n",
    "\n",
    "This is sometimes called a 'plug-in' approach, as we are plugging-in differentially private answers for each of the terms in the original formula, without any additional modifications, and using the resulting answer as our estimate while ignoring the noise processes of differential privacy. While this 'plug-in' approach does result in a differentially private value, the utility here is generally lower than the solution in SmartNoise.  Because the number of terms summed in the numerator does not agree with the value in the denominator, the variance is increased and the resulting distribution becomes both biased and asymmetrical, which is visually noticeable in smaller samples.  Here's an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-15T16:41:55.448318Z",
     "iopub.status.busy": "2021-04-15T16:41:55.447646Z",
     "iopub.status.idle": "2021-04-15T16:42:07.267832Z",
     "shell.execute_reply": "2021-04-15T16:42:07.268277Z"
    },
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Plug-in mean privacy usage:    approximate {\n",
      "  epsilon: 1.0\n",
      "}\n",
      "\n",
      "SmartNoise mean privacy usage: approximate {\n",
      "  epsilon: 1.0\n",
      "}\n",
      "\n"
     ]
    }
   ],
   "source": [
    "n_sims = 1000\n",
    "history = np.zeros(shape=(n_sims, 2), dtype=float)\n",
    "\n",
    "for i in range(n_sims):\n",
    "\n",
    "    with sn.Analysis() as analysis_plug_in:\n",
    "        data = sn.Dataset(path = data_path, column_names = var_names)\n",
    "        age = sn.to_float(data['age'])\n",
    "\n",
    "        dp_mean = sn.dp_mean(\n",
    "            data=sn.to_float(data['age']),\n",
    "            privacy_usage={\"epsilon\": 1.0},\n",
    "            implementation=\"plug-in\",\n",
    "            data_lower=0.,\n",
    "            data_upper=100.)\n",
    "        \n",
    "    dp_plugin_mean = dp_mean.value\n",
    "\n",
    "\n",
    "    with sn.Analysis() as analysis_smartnoise:\n",
    "        # load data\n",
    "        data = sn.Dataset(path = data_path, column_names = var_names)\n",
    "\n",
    "        age_dt = sn.to_float(data['age'])\n",
    "    \n",
    "        # get dp release of the size of the dataset\n",
    "        dp_num_records = sn.dp_count(data= age_dt,\n",
    "                                 privacy_usage={'epsilon': 0.5},\n",
    "                                 lower=0,\n",
    "                                 upper=10000\n",
    "                                )\n",
    "        # get mean of age\n",
    "        age_mean = sn.dp_mean(data = age_dt,\n",
    "                          privacy_usage = {'epsilon': 0.5},\n",
    "                          data_lower = 0.,\n",
    "                          data_upper = 100.,\n",
    "                          data_rows = dp_num_records   # Here we're providing the DP release of n.\n",
    "                         )\n",
    "    analysis_smartnoise.release()\n",
    "    dp_smartnoise_mean = age_mean.value\n",
    "    \n",
    "    # store simulation results \n",
    "    history[i,0] = dp_plugin_mean\n",
    "    history[i,1] = dp_smartnoise_mean\n",
    "    \n",
    "print(\"Plug-in mean privacy usage:    {0}\".format(analysis_plug_in.privacy_usage))\n",
    "print(\"SmartNoise mean privacy usage: {0}\".format(analysis_smartnoise.privacy_usage))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-04-15T16:42:07.274753Z",
     "iopub.status.busy": "2021-04-15T16:42:07.274243Z",
     "iopub.status.idle": "2021-04-15T16:42:07.509781Z",
     "shell.execute_reply": "2021-04-15T16:42:07.509271Z"
    },
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD4CAYAAADsKpHdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA81klEQVR4nO3dd3ib1dn48e/Rsrx34izH2XsRk0AhNJQWKG0JXYyXFSgFyqaUFmhpeWnpmx8vHYyXUkYJ0JRQAoQQCBRoWCWQvXccJ7azvG15aZ3fH48syduxJUuy7891+cozzvPojuLo1hnPOUprjRBCCGGKdABCCCGigyQEIYQQgCQEIYQQPpIQhBBCAJIQhBBC+FgiHcDJysrK0nl5eZEOQwghYsqGDRvKtNbZnZWJuYSQl5fH+vXrIx2GEELEFKXUoa7KSJOREEIIQBKCEEIIH0kIQgghgBjsQxBCnDyXy0VxcTGNjY2RDkWEmd1uZ/jw4Vit1pO+VhKCEANAcXExycnJ5OXloZSKdDgiTLTWlJeXU1xczKhRo076emkyEmIAaGxsJDMzU5JBP6eUIjMzs8c1QUkIQgwQkgwGht78O0uTkRB9pLrBxeHyeqYMTcFkityHc949b4ft3oWLvhW2e4vwkxqCEH2gptHFtx77lO888RkPrtwZ6XAi4qGHHmLKlClMnz6dmTNn8uWXX4bldX7/+9+32FdKcdddd/n3H3nkER544IFO7/HUU0/x4osv9iqOwsJClFL86le/8h8rKyvDarVyyy239Ore4SIJQYg+8PxnhRRXNgCw+PNCGpyeCEfUt9asWcPKlSvZuHEjW7du5YMPPmDEiBEhfQ2tNV6vt01CiIuL4/XXX6esrKzb97rxxhu56qqreh3TqFGjePvtQI3s1VdfZcqUKb2+b7hIk5EQYVbX5Gbx5wdbHPt4bynnT82JUEQBL//4tF7f47JnvuiyzNGjR8nKyiIuLg6ArKws/7m8vDwuu+wyVq1ahcVi4emnn+bee+9l//793H333dx44404HA4WLFhAZWUlLpeL3/3udyxYsIDCwkLOO+885s6dy4YNG5gzZw4NDQ3MnDmTKVOmsGTJEiwWC9dffz1/+tOfeOihh1rEVVhYyLXXXktZWRnZ2dk8//zz5Obm8sADD5CUlMTPfvYzHnvsMZ566iksFguTJ09m6dKl1NXVceutt7J9+3ZcLhcPPPAACxYsaPP3TkhIYNKkSaxfv578/HxeeeUVLr74Yo4cOQJAaWkpN954I4cPHwbgz3/+M2eccQZr167l9ttvp7Gxkfj4eJ5//nkmTJjA4sWLWbFiBfX19Rw4cIDvfve7PPzwwz3+t2tNaghChNmq7ceorHe1OPbejmMRiiYyzj33XIqKihg/fjw33XQTH3/8cYvzubm5bN68mXnz5rFw4UKWLVvGF198wW9+8xvAGFv/xhtvsHHjRlavXs1dd91F8/K/+/bt46abbmLHjh08//zzxMfHs3nzZpYsWeK//80338ySJUuorq5u8bq33norV199NVu3buXyyy/ntttuaxP7okWL2LRpE1u3buWpp54CjOavr33ta6xdu5bVq1dz9913U1dX1+7f/dJLL2Xp0qUUFRVhNpsZOnSo/9ztt9/OnXfeybp163jttde47rrrAJg4cSKffvopmzZt4sEHH+S+++7zX7N582ZeeeUVtm3bxiuvvEJRUVG3/x26IjUEIcJs7/HaNsc+2HUcl8eL1TwwvpMlJSWxYcMGPv30U1avXs0ll1zCokWLWLhwIQAXXnghANOmTcPhcJCcnExycjJxcXFUVVWRmJjIfffdxyeffILJZKKkpITjx48DMHLkSE47rfOaTkpKCldddRWPPfYY8fHx/uNr1qzh9ddfB+DKK6/k5z//eZtrp0+fzuWXX85FF13ERRddBMC//vUvVqxYwSOPPAIYw3oPHz7MpEmT2lx//vnnc//99zN48GAuueSSFuc++OADdu4M9CnV1NTgcDiorq7m6quvZt++fSilcLkCXyjOOeccUlNTAZg8eTKHDh0KWfObJAQhwuxgWdtvjrWNboorGxiVlRiBiCLDbDYzf/585s+fz7Rp03jhhRf8CaG5KclkMvm3m/fdbjdLliyhtLSUDRs2YLVaycvL84+1T0zs3nt4xx13cMopp3DNNdecVNxvv/02n3zyCW+99RYPPfQQ27ZtQ2vNa6+9xoQJE7q83mazMXv2bP7whz+wc+dOVqxY4T/n9Xr54osvsNvtLa655ZZbOPvss3njjTcoLCxk/vz5/nPB74/ZbMbtdp/U36czkhCECLPCoIQQZzHR5PYCUFxZH/GE0J32/1DYs2cPJpOJcePGAUazx8iRI7t9fXV1NYMGDcJqtbJ69WoOHep4Jmer1YrL5WozdUNGRgYXX3wxzz33HNdeey0AX/nKV1i6dClXXnklS5YsYd68eS2u8Xq9FBUVcfbZZ3PmmWeydOlSHA4H5513Ho8//jiPP/44Sik2bdrErFmzOozprrvu4qtf/SoZGRktjp977rk8/vjj3H333f73ZebMmVRXVzNs2DAAFi9e3O33qbcGRn1ViAjxejWHKur9+9OGpfq3iyoaIhFSRDgcDq6++momT57M9OnT2blzZ5dDP4NdfvnlrF+/nmnTpvHiiy8yceLEDstef/31/mae1u66664Wo40ef/xxnn/+eaZPn85LL73Eo48+2qK8x+PhiiuuYNq0acyaNYvbbruNtLQ07r//flwuF9OnT2fKlCncf//9ncY/ZcoUrr766jbHH3vsMdavX8/06dOZPHmyv4/i5z//Offeey+zZs0KaQ2gK6q5YyZW5Ofna1kgR8SKkqoGzlj0bwCS7RbOnTyY1zaWAPCT+WP4xfkdf7CF0q5du/zt2/JgWv8X/O/dTCm1QWud39l10mQkRBgFNxflpNjJTg60FRcF1Rz6knxoi45Ik5EQYVRYHpQQUu0MSg50CBZVDpwmIxEbJCEIEUZtawiBhFBSGZkaghAdkYQgRBgdDmoWGpxiJyPBhtk3sV2Zw0m9s+86DIXoiiQEIcKotLbJv52ZaMNkUmQl2fzHiqXZSEQR6VQWIoxKHYGEkBpvjIvPTorjeI1xvKiinvGDk/s2qAdSuy7T43tXd11GRC2pIQgRJlprymqd/v3UBCMhZCYF+hFOBNUg+juz2czMmTOZOnUqP/zhD6mvN5rTkpKSQvYa1113XYupILqyePFilFJ88MEH/mPLly9HKcWyZctCFleskIQgRJjUOT00uIxprq1mRbzVDECKPVAxL3cMnITQPOnc9u3bsdls/oewQunZZ59l8uTJJ3XNtGnTWLp0qX//5ZdfZsaMGaEOLSZIk5EQYVJW27K5qHlpw9T4QB9CmcPZ5ro+dfXK3t/jhW+f9CXz5s1j69atLY599NFHPPLII6xcacR0yy23kJ+fz8KFC3nnnXf46U9/SmJiImeccQYFBQX+csHmz5/PI488Qn5+PklJSdx+++2sXLmS+Ph43nzzTQYPHtxuLJ9++ikul4umpib279/PzJkz/ec3bNjAT3/6UxwOB1lZWSxevJghQ4bwzDPP8PTTT+N0Ohk7diwvvfQSCQkJLFy4kJSUFNavX8+xY8d4+OGH+cEPfnDS71EkSA1BiDBpr/8AICU+8D2sbADVEJq53W5WrVrFtGnTulW+sbGRG264gVWrVrFhwwZKS0u7dV1dXR2nnXYaW7Zs4ayzzuKZZ55pt5xSiq9//eu89957vPnmm/6ZVwFcLhe33nory5YtY8OGDVx77bX88pe/BOB73/se69atY8uWLUyaNInnnnvOf93Ro0f57LPPWLlyJffcc0+34o0GkhCECJPgGkJaQqBWEJwcyiNdQ+hDzQvX5Ofnk5uby49+9KNuXbd7925Gjx7NqFGjALjsssu6dZ3NZuPb3zZqL7Nnz6awsLDDss1rFixdurTF/ffs2cP27dv5xje+wcyZM/nd735HcXExANu3b2fevHlMmzaNJUuWsGPHDv91F110ESaTicmTJ/un6Y4FYWsyUkqNAF4EBgMaeFpr/WirMgp4FLgAqAcWaq03hismIfpSWQc1hODtgVRDaO5D6IjFYsHr9fr3m6e37sx5553H8ePHyc/P59lnn21xzmoNNNN1NU30nDlz2LZtGwkJCYwfP95/XGvNlClTWLNmTZtrFi5cyPLly5kxYwaLFy/mo48+8p8LnqI6luaLC2cfghu4S2u9USmVDGxQSr2vtQ4eAvBNYJzvZy7wF9+fQsS84GcQ0jpICOV1Ea4h9KD9P1xGjhzJzp07aWpqoqGhgQ8//JAzzzyTCRMmUFBQQGFhIXl5ebzyyiv+a957772Qvf6iRYvarEswYcIESktLWbNmDaeffjoul4u9e/cyZcoUamtrGTJkCC6XiyVLlvinq45lYUsIWuujwFHfdq1SahcwDAhOCAuAF7WRQr9QSqUppYb4rhUippUGNQcFJ4FkuxWFUW2urHfi9nixDJCV0zozYsQILr74YqZOncqoUaP86wvEx8fz5JNPcv7555OYmMipp54altf/5je/2eaYzWZj2bJl3HbbbVRXV+N2u7njjjuYMmUKv/3tb5k7dy7Z2dnMnTuX2tq2K+PFmj6Z/loplQd8AkzVWtcEHV8JLNJaf+bb/xD4hda6w/mtZfprESuue2E9H+wy2o/vOGccc0dn+s9d/9J6ahuNJoy1vzyHQcn2du8RKi2mQ47BB9McDgdJSUlorbn55psZN24cd955Z1heqz+I2umvlVJJwGvAHcHJ4CTvcT1wPRiLcQsRC4L7B948+AveL7Hwq3nG2PsUu9WfEMpqnWFPCC3E4NPEzzzzDC+88AJOp5NZs2Zxww03RDqkfimsCUEpZcVIBku01q+3U6QECF4derjvWAta66eBp8GoIYQhVCFCLrgP4cqp97f40E+Nt1JSZcxjVF43cDqWe+rOO++UGkEfCFvDpW8E0XPALq31HzsotgK4ShlOA6ql/0D0FxVBHcYj04eRmRB4KCoSI41iabSL6Lne/DuHs4ZwBnAlsE0ptdl37D4gF0Br/RTwDsaQ0/0Yw06vCWM8QvSZRldg2gqzUmw5/m+UUpw2/BtA3z+LYLfbKS8vJzMz0z8UU/Q/WmvKy8vbjJbqrnCOMvoM6PQ3zze66OZwxSBEpFTWBz7kk+wWPiw0WkybE0JKixpC+BPC8OHDKS4u7vZTviJ22e12hg8f3qNrZS4jIcKgss7l306Ka/vfLHj6ir6Y4M5qtfqf9BWiIzL4WYgwqAqqISTb2yaE5LhADaGy3tXmvBCRIAlBiDAI/pBvr4aQFJQkgpOHEJEkCUGIMKjssoZgabesEJEkfQhChEHwt/6kOAuXzVzU4nxyixqCNBmJ6CAJQYgwCG4ySrZbSY5La3E+uBmpqsGF16sxmWQ4qIgsaTISIgxaDzv95NBKPjkUWOHLYjb5l9T0eLV/GgshIklqCEKEQXAzUHKchXcPr8TscnDV8fWUjboIR/ZMkuIs/ofXKuudpCZYO7qdEH1CaghChEHrGoLJ3UBixS6G7HmRSf++BpPL0aIfQTqWRTSQhCBEGATXEFItHhKq96EwVgOzOKvJPvhmq6Gn0rEsIk8SghBhEPyNf+LRNzC7G1qcH7zn7yTbzO2WFyJSJCEIEWIer6a6IfCNf1B52wWdEqv2MJPd/n15WllEA0kIQoRYTYOL5hmIE2xmkiq28VxDNc81VFOfMsZfbk7j5/5teVpZRANJCEKEWHDzzzBbPfa6EuKBOJOV8rxv+c+NcBX4t4PXThAiUmTYqRAhVhXUXHSKtRDc8HeLHWd8NmcmB2YczWnYD2hASaeyiApSQxAixIKbf6YpoxbwjjWOVSYzrvhsPOY4ABLcVWRjrG8sncoiGkhCECLEgr/tT/Ae8G97LXZQJpqSAsuITzAVAdKpLKKDJAQhQiw4IYxx7fNve83GsoaNSbn+YxPVYd81UkMQkScJQYgQa+5DSMFBhsdYslKj8PqaioJrCJNMRkKQJiMRDSQhCBFi1b4P91x1wn9Mm6zgW9y+MTlQQ5igjCajRpeXRt+8RkJEiowyEiLEmmsIwQnh6fgRFE25DYDGoBrCeFMxZjx4MFNV7yIn1YwQkSI1BCFCrLkPYYQq9R9zxWf7t73WJFxx6QDYcDNUlQHSbCQiTxKCECHWXg1hscfBG0c+9O874wf5t5vLSUIQkSYJQYgQa+5DGBGUEP7jrGJd5Q7/vqudhCAPp4lIk4QgRIg11xCCE4I221qUcQY1ITU3LUkNQUSaJAQhQsjrm+nUhJdhvr4BAK+p5Wpo7TUZSQ1BRJokBCFCqLbRjdaQQwU2ZQwjddtSQLX8rxbcZDTCnxCkhiAiS4adChFCVQ2+ZxBMgeYiZ/wgbJ3UEEb4O5WlhiAiSxKCECEUGHLaMiE8MOknLcq549LwmqyYvC4ylINk6qWGICJOmoyECKH2OpSDn0HwU6YWx0eoE1JDEBEnCUGIEGr+lh/coeyyZ7O0+F2WFr/boqzT3nKkkYwyEpEmCUGIEGpeS3kIFf5jTnsmW6v3srV6b4uyzoSW/QjVUkMQESYJQYgQqqzzJQRV7j/mtme2W7b1w2lVDS5082LMQkSAJAQhQsho9tEMUYEagsue0W7Z4Caj4aoUj1dT0+gOd4hCdEgSghAhVFXvJJ1a7MqoKXgs8XgtCe2WDe5Ubu5zqKyTfgQROZIQhAihynoXQ4NrB3FG7SDZkkiyJbFFWZc9y789XJUCmnJJCCKC5DkEIUKoqt7Zov/A5es/uHfCj9qU9ViT8JjtmD2NJKom0nBQIQlBRFDYaghKqb8ppU4opbZ3cH6+UqpaKbXZ9/PrcMUiRF+prHeRE1RD6KhDGQCl2tQSKuqawhmeEJ0KZ5PRYuD8Lsp8qrWe6ft5MIyxCNEnKuucDA2uIfiajF44tIIXDq1oU751P4I0GYlICluTkdb6E6VUXrjuL0S0cXm81Da5ybEGjzAyagh7HIXtX9OqhiCdyiKSIt2pfLpSaotSapVSakqEYxGiV5rnMRraTh9CR5zxwQlBaggisiKZEDYCI7XWM4DHgeUdFVRKXa+UWq+UWl9aWtpRMSEiqnnaiiG0bTLqiKvFswhl0qksIipiCUFrXaO1dvi23wGsSqmsDso+rbXO11rnZ2e3M1GYEFHAmJxOd79TGXAF1RCGSUIQERaxYadKqRzguNZaK6XmYCSn8i4uEyJqVdY7yaSGOGU8beyxJOC12AHItKW1e03rp5XLHZIQROSELSEopV4G5gNZSqli4DeAFUBr/RTwA+AnSik30ABcqmUiFxHDOnoGAeCucVe1e43HluJfFyFF1eOqqwx7nEJ0JJyjjC7r4vwTwBPhen0h+prxlHL3O5QBUAqnPRt7/REAst3HaHB6iLeZwxWmEB2K9CgjIfqNynpny0nt4gIJ4ZmDr/HMwdfavc6VMNi/PVIdo0LWRRARIlNXCBEiVXUuRnVQQzhYX9Lhdc4WCeEEFQ4nw9LiwxOkEJ2QGoIQIVJR3+op5Q6mvW7NGd+yhlAu01eICJGEIESIlDuaTr4PgbY1hDIZaSQiRBKCECFS5nB2a6W01oJrCLmm45yobQx5bEJ0h/QhCBEiFY4GBqvAsNHgp5SH2ge1d4lRLj4bLwoTmiFUUF5ZHdY4heiIJAQhQqDe6SbJWYbF7gXAbU1Bm23+87eMubTDa7XJQr01gyRXOSal8VQcAmaHO2Qh2pAmIyFCoNzRsw7lZo1BNQhrTWGowhLipHQrISilXldKfUspJQlEiHaUOppaPoPQqv/giQNLeeLA0g6vdyXm+LcTHYdDH6AQ3dDdD/gngf8C9imlFimlJoQxJiFiTlltU6tpK1rO03ik8QRHGk90eL03KZAQMp0leL0yi4voe91KCFrrD7TWlwOnAIXAB0qpz5VS1yilrOEMUIhYUNbLJiNP4hD/9hiK5WllERHdbgJSSmUCC4HrgE3AoxgJ4v2wRCZEDCl3NDFMlfn3u/sMQrPG5JH+7UmmwxyraghZbEJ0V3f7EN4APgUSgO9orS/UWr+itb4VSApngELEgjJHE7nquH8/+NmC7nDZs2jAmCo7XTmoPnEopPEJ0R3dHXb6jG8RGz+lVJzWuklrnR+GuISIKWW1TeSqQB+BK77lcwejEoZ1fgOlKLUOIdd1EAB3yVY4ZUbI4xSiM91NCL8D3ml1bA1Gk5EQA56z5jiJypiDyGWOx2NtWXH+8ajvd3mPyrjh/oRgKd0R+iCF6EKnCcG3qtkwIF4pNQtQvlMpGM1HQgjAHjRUtDEuG5TqpHT7HIkjwGFsJ1ftClVoQnRbVzWE8zA6kocDfww6XgvcF6aYhIg5SUHTWwevb9DsD/teBDpeOQ3AnTIKfN0QWXX7QhugEN3QaULQWr8AvKCU+r7Wuv3VPYQY4Oqa3GS7jvgWiAVPYtuEUO6s6vI+1qyRePYqzEqT4zkCTQ6IkzEbou901WR0hdb670CeUuqnrc9rrf/YzmVCDCglVQ3kmoI6lNupIXRHcmIiR8hiBKWY0NTt/5zEKeeGKkwhutTVsNNE359JQHI7P0IMeCWVDa2GnHY8s2lnTEqx3zLWv1+/fWWvYxPiZHTVZPRX35//3TfhCBF7iqsa+EbQkNOeJgSAwwlToXYNAIkH3wOte9RBLURPdPfBtIeVUilKKatS6kOlVKlS6opwBydELDheXkmObx0EL6Y28xgBTEjKY0JSXpf3qkudiEMbD6glNB6DY9tCGqsQnenu1BXnaq1rgG9jzGU0Frg7XEEJEUtcpQf823XWDDCZ25S5euSFXD3ywi7vlZYcz2bvmMCB3dJsJPpOdxNCc9PSt4BXtdaypJMQPvEVgWcG6rt6IrkLWUk21uugyYTXPAlVgWccnG5vr+4vRGe6mxBWKqV2Yyzj9KFSKhuQhV+FALIce/3bTcm57Zb5nz3P8T97nuvyXoOS7WxmAke1b7ZUZy28dh2NJwq49/WtTPnNu1z81zUUltWFJHYhgnV3+ut7gK8A+VprF1AHLAhnYELEAqfby0h3gX/fmz6q3XK17jpq3V1/iFvNJoZmpPCk+0K82teZXPQl1idnM33jr8nylLH2YAUXPvEZRRX1Ifk7CNHsZFZAmwhcopS6CvgBIAOkxYB3rKqBSSowM6krJa/X9xw3KJl9ejj/9HzVf8yMl8ssq1kVdw9jVAk1jW5eXFPY69cSIlh3Rxm9BDwCnAmc6vuRWU7FgHe4qIAsVQNAI7ZeDTltNm6w8XTycu+ZPOi6gh3ewFoJaaqOP1mfxIKbNzaV4PJIn4IIne7OdpoPTNZay7p+QgSp2L/Bv11mHQohWHY8PcFGToqdYzWN7NR57HTnMUMd4G7rP7HgYbrpIDeb3+RRx/f5eE8pX5/csyejhWitu7+924GcLksJMcB4jgaeE6hLbL9DGWB66nimp47v9n2/O2sYafGB1WkbB83k+NhL/PsLLe9hxc2q7cdOMmIhOtbdGkIWsFMptRZoaj6ote56YLUQ/Vh6dWDdAk9a+x3KAJcOP/+k7puVFMePzhzFxkNVJMebmTYsjSqGkV38L2yNZaQrB2eYtrGpKLXHsQvRWncTwgPhDEKIWFTtqOMU92b/KiGWwRNxhfD+CTYLZ44LfupZUTP4NLIOGQ+rXWhew09LZ1Hd4CI1qDYhRE91d9jpxxhPKFt92+uAjWGMS4ioV7L5A1JUAwDlKg1XB88gADyw6y88sOsvvX7N6pyv+LfPNa3HThNbi6t6fV8hoPujjH4MLAP+6js0DFgeppiEiAlNO972bxfYp3Y6CZ3T68Lp7X39oTF5JE0JQwBIUo2cZdrKlqKqXt9XCOh+p/LNwBlADYDWeh/Q+/F1QsQot9tDztF/+/frBs3umxdWippBgRHfp5r2sFkSggiR7vYhNGmtncr3DUgpZQFkCKoYMD7ac4KH3t5FXZOb+RMHkVOxjtsoBaCeONJHz+qzWOrTAqOVZpv28myJTC0mQqO7CeFjpdR9QLxS6hvATcBb4QtLiOjx8trD3Pt6YHjpK18e5G3bH/3164OJM7FYbH0WT0PqOP/2VHWQqppaahtdJNulY1n0TnebjO4BSoFtwA3AO8CvwhWUENGi3unm4Xd3tzh2k/lNJpqKAGjESuPkS7u8z6npUzg1fUpIYvLYUvz9CDblYao6SEGpTHYneq9bNQSttVcptRxYrrUu7c41Sqm/YayfcEJrPbWd8wp4FLgAqAcWaq1l5JKIKv9cV0RlvdEZnBRn5tnc9zmtaJn//LGRC7Cndf2k8HeHnhPSuOrTxhNXfxQwmo0OlDqYMSItpK8hBp5OawjK8IBSqgzYA+zxrZb2627cezHQ2dM43wTG+X6uB3o/Jk+IENJa89x/Dvr3nxjyHqcVPevfr0ubQP3Y70QiNOpTg/sR9nGg1BGROET/0lWT0Z0Yo4tO1VpnaK0zgLnAGUqpOzu7UGv9CVDRSZEFwIva8AWQppQachKxCxFWB0rrKKownjP4ke0D5h/7m/9cbeYMDp3yC7Spe+329+14jPt2PBay2II7lmeZ9nPghDQZid7rKiFcCVymtfZ/TdJaFwBXAFf18rWHAUVB+8W+Y20opa5XSq1XSq0vLe1Wi5UQvfb5gTIAsqjmZ+Z/+I/XZs6gaOZP0WZ7pELDmTgEj8noyB6kqig/XhyxWET/0VVCsGqty1of9PUj9NmQBq3101rrfK11fnZ2dl+9rBjgPttn/OrfYVlGvDYWCGxMHE7RjDu6XTMIG2WiIWmEfzehajdumQpb9FJXCcHZw3PdUQKMCNof7jsmRMR5vJo1BeWMVMe41Lzaf/z4+P9Cm+MiGFmAM2iqjLH6EMWVDRGMRvQHXSWEGUqpmnZ+aoFpvXztFcBVvo7r04BqrfXRXt5TiJDYeaSG2kY3PzB/gkUZ37wd6ZNxZM6IcGQBTUmBhDDZdJjDsqSm6KVOh51qrc09vbFS6mVgPpCllCoGfoOvmUlr/RTGswwXAPsxhp1e09PXEiLUtpVUA5rvmNb4j1XkntfpfEWdOTMz9E8yNwbVECaqw2yShCB6qbtPKp80rfVlXZzXGHMkCRF1dhypZroqIM90HACPJb5XtYMLcuaFKjS/xqAawlhVwsqympC/hhhYer/enxD90I4jNVxo/ty/X5N9Ktrc8+kpmjxOmjy97XZryWtNpM6aAUCcctN0bE9I7y8GHkkIQrTi8Wp2H6vmfPM6/7HqnNN7dc//3v0U/737qd6G1kZdwnD/tr1iZ8jvLwYWSQhCtFJQ6mCw+yjDlTHs1GO2U5cRmnmIQs2dMtK/nV53AKMlVoiekYQgRCs7jtRwhimwVnJ9+iQwha27rVc8KYGR23neIqrqQ7mIpxhoJCEI0cruY7V8xbTdvx+ttQMAZ3KgyWicKpahp6JXJCEI0cqB49V8JaiG4MhoM1lv1GhKGIoXYyhsrjpB8YnyCEckYll01oOFiCB9bDsZypg91GlJoSlpRBdXdO2c7Lm9vkd7tNlGtTmLdE8pJqVxlOyE2WPC8lqi/5OEIESQRpeHkY5N/v8ZdRmTe/wwWrBzBoUnIQDU2IeQXmdM+qhKdwGRmZJbxD5pMhIiyMGyOmarwHj+hozJIblvjctBjSs8axY0JAYmCY6v2h+W1xADgyQEIYLsO17LqaZAQghed6A3Fu39G4v2/q3rgj0QPNIoo74gLK8hBgZJCEIEKTu8m2xVDUCTstOUNLyLKyLPlJ7n3x7hPoTHK88iiJ6RhCBEEHPJl/7t8oTRoKL/v4hOGYYnaKTRiTIZaSR6Jvp/24XoQ4MqNvm3HWkTIhhJ92mTlVKV5d8vP7glgtGIWCYJQQgfj1cz3hl4IE1nh6ZDuS+UWwPLkTeWbItgJCKWybBTIXyOlhQxRh0BwI0JT8a4kN37m4PPDNm92lNjHwbOrQCYS2WSO9EzkhCE8Knc8zHNXcjFpuG9mu66tXlZp4TsXu1pTBoBvuUQkqr3hvW1RP8lTUZC+KhDgfUPjtnHhvTepU2VlDZVhvSewbxpef7twQ0HQGY9FT0gCUEIn/Sy9f7t6pTQPH/Q7E/7X+JP+18K6T2DxaXl0KitAKToGnCcCNtrif5LEoIQAI3VDGkwnvL1anBlTYxwQCcnNSGOIj3Iv990RDqWxcmThCAEQNFaTHgBOKRzSEtNi2w8J8mkFMfNg/37NYWbIxeMiFmSEIQAmvb+27+9h1xS460RjKZnyuMCT1W7SzZ1UlKI9klCEALw7n3Pv33INgFTCGY47Ws1iXn+7cRSeThNnDwZdipExUHiqw8A4NQWKlJD339w0ZCzQ37P1pwpI3GVmbEqDykNRVBfAQkZYX9d0X9IDUGI/R/4N7frPNKSk0L+EnMypjEnY1rI7xssNTGRQzrQj8CRjWF9PdH/SEIQYt+//JubvWPJTooL+UsUNxynuOF4yO8bLCPRxgE9NHCgRBKCODmSEMTAVl8BBR/7dzd7x5CdHPqE8GTBKzxZ8ErI7xssPcHGAR2Y08hbvCGsryf6H0kIYmDbvAQ8TQAUeHMoI4P0xNgbYQRgs5g4Zs3173uLN8gTy+KkSEIQA5fXC+ue8+++751NRqINcwysgdARd9JQHNoOgKWhFEr3dHGFEAGx+5svRG/tew8qDwLg0Hb+453KoJTQNxf1pczkeLbrUYED+9+PXDAi5khCEANTXRm8dYd/9xPvdJxYyUmxRy6mEMhKsrHFOzpwIGgElRBdkecQxMDiaoTitfCv+8FxDIA6lcgKz+kA5KSGJyFcPOy8sNy3taykON72jgkcOPQ5OOvAltgnry9imyQE0f95XLDpJdj6KhR9CdrT4vT/uS+kimQgfAlhZh8tx5mZZKOCFA57s8k1lYLHCYWfwfi+SUgitkmTkejfmhyw5Iew8k44/HnLZKDMVE6+io0e4xt1st1Coi0835EK6oopqCsOy72DJdksxFlMbNZB6zlsezXsryv6B0kIov/yuODv34OC1UEHFaQMhXHnwkVPsS0tMKVETkp82EJ5tvB1ni18PWz3b6aUYlCync+9QetB73gDao6G/bVF7JMmI9F/rX3GaCJqNvUHMHkBxKf7D+3bWuDfHpIa2yOMmg1JtbO2cgi7vSOYaCoCrxvW/w2+9stIhyainNQQRP/kKIWPFgX2Z/wXzF7YIhkA7DhS49/OzegfHa9DfP0gqzxzAgfXPWuMrBKiE5IQRP/0yf9CU7WxnTIMpv2wTZHKBidHqhsBMCvF8PTwNRn1paGpxt9jnZ5AOWnGwYYKWP4TeXJZdCqsCUEpdb5Sao9Sar9S6p52zi9USpUqpTb7fq4LZzxigHDWweZ/BPbzfwTmttNR7DxS698elhaP1dw/vh9lJNmwmRVeTPzV9c3AiX3/gn/9SpKC6FDY+hCUUmbg/4BvAMXAOqXUCq31zlZFX9Fa3xKuOMQAtP01cPo+7FOGwfBT2y2240i1fzsvKyGsIV2Z++2w3j+YSSlyUuI5XFnPVj2GY8PPJ6f4XePkmieMCf0ufBzM0oUoWgrnV6I5wH6tdYHW2gksBRaE8fWEMGxYHNgefx60s/qZ26tZX1jh3x+ZGd7+g0nJo5mUPLrrgiEyLKj561/282HEaYGTW/4B/7wSXA19Fo+IDeFMCMOAoqD9Yt+x1r6vlNqqlFqmlBrR3o2UUtcrpdYrpdaXlpaGI1bRXxzfCSW+aZ9NFhhzTrvFNh6upKbRDRjPH+RmhLeGsKu2gF21BV0XDJEx2YFFfjYW18L8e42hts32vAMrbpXmI9FCpBtN3wLytNbTgfeBF9orpLV+Wmudr7XOz87O7tMARYzZ+WZge8RpYE9tt9hHu0/4t2cMSwv7GsovHV7JS4dXhvU1go3IiMdqNv5Ox2qaOFrrgtNvhanfDxTa9qoxHFUIn3AmhBIg+Bv/cN8xP611uda6ybf7LDA7jPGIgWDXisB23pntFtlcXMXGoir//owRaeGNKQIsJhOjs4JqCYcrjaaz2de0rCm8ey9UFvZ9gCIqhTMhrAPGKaVGKaVswKXAiuACSqkhQbsXArvCGI/o78r2wQnfmAVzHAwLfL/QaMrrmvh4bymPfbjPf3xiTjIZiba+jrRPjA1qNvpoTyne5uahOTdAep6x7Wlq+byGGNDCNsxAa+1WSt0CvAeYgb9prXcopR4E1mutVwC3KaUuBNxABbAwXPGIASC4uWjYKWCNp8Hl4Y1NJXy8p5TqRleL4gk2MxdMHUJ/NWloCu/vOobToymqrGdTURWzc9PBEgdzb4R3fSPBtyyFM26HQZMiG7CIuLD2IWit39Faj9daj9FaP+Q79mtfMkBrfa/WeorWeobW+myt9e5wxiP6ud1BbfQjz6CizsnPX9vKii1H2iSDeKuZ758ynMS4/jv0Mt5q5pTcDP/+0rVF1DuNjnQGTw2qQWn4+OG+D1BEnf77v0EMLDVH4MgmY1uZaco5hf99bw+ltU3+IjazYnBKPLkZ8Zw+Oot4m7nPwrsu73t99lrB5o7OYN2hCjxeo5bw0Du7uf6s0YzMSIBZVwZGZO1aAY4TkDQoInGK6CAJQfQPe98NbOdMY/mOag6W1QGggAumDWHmiPCPJurI6MThEXndFLuVC6bm8NZWY7bTA6UOfvHaVmYMT+PK03IZPmiy0e/idcPGF+Gsn0UkThEdIj3sVIjQ2LPKv1mfM5u3twWmez53cg6n5KZHLBkAbK7aw+aqyCx4P3NEOl+f1PKb/5biKu55Yxu7ks8IHNzwAng9iIFLEoKIfU0OKPjYv/tOZS5OjxeAwclx5Oeld3Rln/lnyXv8s+S9iL3+6aOzuOGs0UwYnOw/5vZo/mdXJm6rbzRS9WEo+CgyAYqoIAlBxL4DHxrDJwFv2kjePuj1n5o/cVBEawbRZFCynYvzR3D9WaPJTjLWfnBi4UPnlEAhWV1tQJOEIGLfrsDoosOJ02hwGc0emYk2xgWNxReGwcl2rjgtl1S7MQPsx+6ghLDrLXDWRygyEWmSEERscztbdCivrA5MIDcrNw0ltYN2JcVZuXDmUAAK9BCOaN/wVKfDmOdIDEiSEERsO/gJNBmrnrkTBvFZhdFGblKK6cPSIhhY9MvLTGRSTjKg+MwzNXBCmo0GLEkIIrYFzV1UkDAdY5ApjBuUFFUPnd00+hJuGn1JpMNo46vjjcki/+MNSgj7PzTWTBADjiQEEbtcjS2mq3jfMca/PXlISiQi6tDw+MEMjx8c6TDayE62MzY7ieNksN9rNCHhdbV86lsMGJIQROzauwoaqwBwJ2TzWU0WABaTYnzQ8MposLZiG2srtkU6jHadNjoTgM+9kwMHty2LUDQikiQhiNgVtG7yvqQ5aF9z0dhBSdgs0fWrvfzoapYfXR3pMNqVl5lARqKNNd7JeJvXyyn8FGqPRzQu0fei63+NEN1Vewz2f+DffasuMFPnpChrLop2SilOyU2jkhR261zjoPbCdqklDDSSEERs+uIvxocW4MyawsZKOwBmpRg7SJ49OFnThxvzPH3qnRY4GFQDEwODJAQRe+orYN2z/t1tyfP826OyErFb+m4W0/4i0WZhwuAkvvBOokn7Rmcd3w5Ht0Y2MNGnJCGI2PPFX4wHqADScllekes/NTEnujqTY8msEek0YGedd0LgoNQSBhRJCCK2HN8J/3nUv1sz/nvsKzWmWlDAhChNCHeOvZI7x14Z6TA6NSo7kVS7hU+80wMHt7wMzrrIBSX6lCQEETtcDfD6j/0T2ZExhk9dgc7kUVmJJNii52G0YNlx6WTHRX7W1c6YlGJmbjrb9CiOaV+sjVVGUhADgiQEERsaq+Gl7xnt2gBmG8y7izUFlf4i0Ty66NOyjXxatjHSYXTplNw0lDLzrufUwMEv/gJeb8cXiX5DEoKIfgdWw1/PgsOfB47l/4gSlc3+UqMvwaRUVPcfrDr+GauOfxbpMLqUFGdl4pBkPvLOoF4bU2RTvh92vxXZwESfkIQgolfFQVh6Obx0EVQWBo7nXwsTv8Vne8v8h8YNSora5qJYc/qoTBqJ4wPvKYGD7/8a3E0dXyT6BUkIIvp4vbDmSXjy9JZz6tiSYN7dMOV7eLTm032BhDB9eGoEAu2fhqbFMzY7iTc9X8Ghjec7qCw0mo5EvyYJQUQXZz28ejW8dy+4GwLHR58NF/0FRn8VgA2FFZTVGd9Y461mxsnDaCH11fHZ1BHPMs9ZgYOrfw8l0d8PInpOEoKIHo018OKFLaa0Jj0PLvgDzLsL4gOjdN7edsy/PTs3DbNJfpVDaWhaPPkj03nfO5tD3kHGQU8T/PMqqDka2eBE2Eijq4gOTQ5Y8gMoXhc4NvE7Rn+B2dqi6M6jNew5XgsYncmz8zL6MtIeuWf8tZEO4aSdM3EwB0od/LH+B/ze+hyJqgmqiyh77Ks8nPUQezxDibOYGJGewPwJ2Zw3JSfqJhUUJ0f+9UTkedxGM1HRl4Fjc66HuTe0SQYerVn8eaF/f+rQFFLsLctEoxRrEinW2GrWsllMXDF3JE3xg3nCfRFubXxcZLmP88DRm5l25FXWHSzjtY3F3PryJr7xp49ZvedEhKMWvSEJQUSW1vDOXS1mLmXO9TDpwnaLr9xylMMVxpPJFpPi7AmD+iLKXvvwxJd8eOLLrgtGmbQEG9eckYd32Kk84r6YRm0k3wTVxO+sz/Oq7UEmqsMAHCqv55rn1/HH9/eite7stiJKSUIQkfX5Y7BhcWB/+iUdJoM1BeW8vO6wf/+MsVmkxEd/7QDgw9Iv+bA09hICGM8mXDhjKF875wLWTryH6rgh/nP5pr28HfdL7o17FStuAB77cB+/Wr4dr1eSQqyRhCAiZ/vrxvj2ZqPnw8wr2hQrdTTx3GcHefTDff5jI9ITOGNMZh8EKZolxlnIzJ1EyRn/Q+moBXiVMausGQ83qDd4N/G/GaNKAFjy5WF+uXybJIUYI53KIjIO/Btevz6wP2gKfOV2UAqPV7PneC2bi6rYVlLNwbKWk6ulJ1j54ezhMrIoQrTZxomxl1CVcyZDdz1LYtUeAMZ4DvC2/X7ubrqOt7xf4eW1RZhNit8umIpSKsJRi+6QhCD63r734ZUrjcXcAVKHw9m/5Fidh/d3HOGz/aVUN7rbvXRMViIXzRomTyVHAWfSMArz7yej6D0G71uKyevCrht53PYE+e49POS+gr9/cRivht8umIrZJEkh2sn/KhFSjiY3ZbVNuL2ajEQb6QnWwLdDrWH932DVz8Hr+8BPzKZg1j0s//QYawsr2r2nCUVeVgJzR2UwJjtJvm1GE2WiIveb1KVPYcTWPxNXbzwfcrXlfWaaDnCn6yb+8SWU1jbxyA9nkBojfT4DlYq10QD5+fl6/fr1kQ5DBNlaXMXyTUf4aM8JCtpp3pk6LJX52XV859gTDDryof9cvTWTJ+MWsr4ivs09k+KMFbxGZyeRl5mI3Rrbq6A1eZwAxJltEY4kfEzueobueJrUE2v9xxq0jSfcF/Gs5wIyU1O469wJfGfGUHleIQKUUhu01vmdlpGEIHrC7fGyavsxnv3sIFuKqtotY8LLHNNuLjGv5lumL7Apj//cIe8gFrkvpZKWU1aPyUpkzqgMRmcnYZKaQOzRmozD7zJ43z8w6cC/d5lO4TXPPD7xTqfENooxI/MYlZ3M4BQ7Wck2BqfYGZWVSE6KXWqAYSIJQYTc0eoGXt9Ywj++PExJVUOb84NUDecl7OJ0vZm5nk1kqpo2Zf7lmc1Lnq/jwmg+MJsUU4emMnd0BoOT7WH/O0TCO8c+BeCCnHldlOwf7LWFDNv+F+yOonbPe7WiERtuTHgw04CNKp1MhSmDptRRJOVOZ/yseaSPnNHm4UTRM91JCNKHIFrQWlPb5KbC4aSy3kmZw8nR6gb2n3CwvrCSnUdbfsDH4WSOZR8Xp+/nNL2FbMduaP5i2OqL3kFzHq/yDQ6YhpGZaCErKY68rAQm5aTEfJNQVz4r3wQMnITQmJzHgbm/J71kNdkHl2Ntatk/ZFKaBALTaacDQ1UFcAiqN8G212EbOLFRnTqRxJEzSciZAClDjTmtLL61Grxu8LiM/imTCayJxvnkHLBH74JJ0UoSwkDldlJbeYy9hUUcPHKCohOVHK2q41itC4dL0YSNBmw0aStuzJjxEq+amKNqyDWdYJwqJt9SwHR1AKt2Qm37L+OypVE7KJ/KYWfTmDKKb/Xt31JEkslM5YivUzn8aySVbSG5bDPx1fuw1Jdi9XRvnWYbTrKrt8LWrbD15F7eG5+ByhqPGjwZBk+FnGkwaBLERe9CSpEW1oSglDofeBQwA89qrRe1Oh8HvAjMBsqBS7TWheGMqV9yNRpLTPp/qqChChoqoaECb10ZTVXHcdYcRztOYGsoJcHrIBnjjZ8dfC8TEHcSr92qxVErM/WpY3FkTseROYPGlDxQ0oE4oCkTjuxZOLJnBY55PSjtQnm9KO3G5GnC5KyhsfIojeWHsVUfZIi7iOx2mhy7y9RQAUVfGD9BGuJzaErOxZU4BE9CNq64dFyWJFzmRNyWeLyWBLAlYk5IJT4pnaS0TNJT07FY+nctFsKYEJRSZuD/gG8AxcA6pdQKrfXOoGI/Aiq11mOVUpcC/w+4JFwxUV8BO94A7VsfVnuNH68nUPX0OI1pft1N4Ko3FnZ31oOrzth21RsfwJ4mY1I2rztwPzA+/JQJlAKU8afWgA68nu9HAy6PF682oZUJrzLjVWY0ZrzKhBfjUqU9mLwuTF4XNm8jdt227b4zJiDe9xMuVYmjqUiZTEXKJNxmXz+Aowkce8L4qrHD02TMv1RbsivCkUQ5aybkZOLOmcU2p5sTJ45jrjrIUH2CQaqKVOUgkUYseFCABxMeTHhRmPEShwuL6nz95/iGY8Q3HOu0TFdcWHApG15lQZvMKJTxX97/p7GtlQmTxYbdHo+yxIE13veTCLYE33aC0QRmjvP9aQWTBZQZTObAFyqlYNR8yBrbq9g7E84awhxgv9a6AEAptRRYAAQnhAXAA77tZcATSimlw9HTfehzeOm74G4M+a17SgH9ZRBiWl0BaXUFjD66suvCA1BKvLGi22k7fxvhSGJQFFYwrbixardRQ+48/0AT0L0Wsu751h/g1OtCeMOAcCaEYUDwEINiYG5HZbTWbqVUNZAJlAUXUkpdDzTPc+BQSp30185hySonJ0kNK63XZCfE1rC2WItZ4m2PA4BOh3h0U6y9vxB7MUdzvNWP31C5v+LHBa0OZ9Hqc7MdI7u6d0x0KmutnwaeDsW9lFLrD1V5Q/H/ss/EWswSb3jFWrwQezHHYrxdDSntjnBWxkqAEUH7w33H2i2jlLIAqRidy0IIIfpYOBPCOmCcUmqUUsoGXAqsaFVmBXC1b/sHwL/D0n8ghBCiS2FrMvL1CdwCvIcx7PRvWusdSqkHgfVa6xXAc8BLSqn9QAVG0gi3kDQ99bFYi1niDa9YixdiL+YBGW/MTV0hhBAiPKJwQJcQQohIkIQghBAC6IcJQSllVkptUkqt9O0/p5TaopTaqpRappRK6uC6e5VS+5VSe5RS50VzvEqpPKVUg1Jqs+/nqUjFG3T8MaWUo5PrIvL++l77pGOOpvdYKbVYKXUwKJaZHVx3tVJqn+/n6vbKRFm8nqAyrQec9GW8Sin1kFJqr1Jql1Lqtg6ui5b3t7vxnvT7GxPPIZyk24Fd4J9o/06tdQ2AUuqPwC1A6zmVJmN0aE8BhgIfKKXGax00oXsUxetzQGs9sw/ia611vCil8jEmrGxXhN9f6EHMPlHzHgN3a62XdXSBUioD+A3Gs28a2KCMqWIqwxqp4aTj9WmIkvd3Icbw94laa69SalDrC6Ls/e0yXp+Tfn/7VQ1BKTUc+BbwbPOxoA9XhTGdT3u96AuApVrrJq31QWA/xtQb0RpvRLQXrzLmrPpf4OedXBqR99cXX09jjoj24u2m84D3tdYVvg+p94HzQx1fa72INyI6iPcnwINaG5OSaa1PtHNpNL2/3Ym3R/pVQgD+jPGfvMXsIkqp54FjwETg8Xaua2+ajWHhCbGFP9OzeAFG+aqRHyul+mqS/T/TNt5bgBVa66OdXBep9xd6HjNEz3sM8JCvGfFPypgluLWo+h2m63gB7Eqp9UqpL5RSF4UzyCB/pm28Y4BLfLGsUkqNa+e6aHp/uxMv9OD97TcJQSn1beCE1npD63Na62swmip2Ec7ZVE9CL+M9CuRqrWcBPwX+oZQK62og7cWrlBoK/JCOk1ZE9TLmqHiPfe7F+HJwKpAB/CKccXRXCOId6Ztu4b+APyulxoQtWDqNNw5o9MXyDPC3cMbRXSGI96Tf336TEIAzgAuVUoXAUuBrSqm/N5/0tVcvBb7fzrXdmWYj1Hocr6/ppdy3vQE4AIzv63iBHcBYYL/veIIyHjJsLRLvL/Qi5mh5j5VSf9daH9WGJuB52m9ui5rf4W7Gi9a6xPdnAfARMKu9cuGOF+Pb/uu+Mm8A09u5NmreX7oXb8/eX611v/sB5gMrMWaYHus7poBHgEfaKT8F2IKReUcBBYA5iuPNbo4PGI3xi5nR1/G2c9zRQfmIvr89jDlq3mNgSNDvxJ+BRe2UzwAOYnSUp/u2oznedCDOt50F7AMmRyjeRcC1QcfXRfn72514e/T+9sdRRsEU8IKvqq8wPpR+AqCUuhDI11r/WhtTavwTY60GN3Cz7rsRMCcdL3AW8KBSyoXRtnij1rqig3tGRJS+v52K4vd4iVIqG+N3YjNwoy/efF9c12mtK5RSv8WYQwyMTseojReYBPxVKeXFaKlYpFsuntWXFmHEfCfGPOXXtY43yt7fLuOlh++vTF0hhBAC6F99CEIIIXpBEoIQQghAEoIQQggfSQhCCCEASQhCCCF8JCEIIYQAJCEIIYTw+f/3YmjP8ooXHgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import seaborn as sns\n",
    "\n",
    "data = np.genfromtxt(data_path, delimiter=',', names=True)\n",
    "age = list(data[:]['age'])   \n",
    "true_mean_age = statistics.mean(age)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "sns.kdeplot(history[:, 1], fill=True, linewidth=3,\n",
    "                 label = 'SmartNoise Mean')\n",
    "sns.kdeplot(history[:, 0], fill=True, linewidth=3,\n",
    "                 label = 'Plug-in Mean')\n",
    "\n",
    "ax.plot([true_mean_age,true_mean_age], [0,2], linestyle='--', color = 'forestgreen')\n",
    "#plt.xlabel('DP Release of Age')\n",
    "#plt.ylabel('Density')\n",
    "leg = ax.legend()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that for the same privacy loss, the distribution of answers from SmartNoise's resizing approach to the mean is tighter around the true dataset value (thus lower in error) than the conventional plug-in approach.\n",
    "\n",
    "*Note, in these simulations, we've shown equal division of the epsilon for all constituent releases, but higher utility (lower error) can be generally gained by moving more of the epsilon into the sum, and using less in the count of the dataset rows, as in earlier examples.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}