{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Working with Unknown Dataset Sizes\n",
    "\n",
    "This notebook demonstrates the features built into OpenDP to handle unknown or private dataset sizes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 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": [],
   "source": [
    "# Define parameters up-front\n",
    "# Each parameter is either a guess, a DP release, or public information\n",
    "var_names = [\"age\", \"sex\", \"educ\", \"race\", \"income\", \"married\"]  # public information\n",
    "age_bounds = (0., 120.) # an educated guess\n",
    "age_prior = 38. # average age for entire US population (public information)\n",
    "size = 1000 # records in dataset, public information\n",
    "\n",
    "# Load data\n",
    "import os\n",
    "import numpy as np\n",
    "data_path = os.path.join('..', 'data', 'PUMS_california_demographics_1000', 'data.csv')\n",
    "age = np.genfromtxt(data_path, delimiter=',', names=var_names)[:]['age'].tolist()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By looking at the private data, we see this dataset has 1000 observations (rows).\n",
    "Sometimes the number of observations is public information.\n",
    "For example, a researcher might run a random poll of 1000 respondents and publicly announce the sample size.\n",
    "\n",
    "However, there are cases where simply the number of observations itself can leak private information.\n",
    "For example, if a dataset contained all the individuals with a rare disease in a community,\n",
    "then knowing the size of the dataset would reveal how many people in the community had that condition.\n",
    "In general, any given dataset may be some well-defined subset of a population.\n",
    "The given dataset's size is equivalent to a count query on that subset,\n",
    "so we should protect the dataset size just as we would protect any other query we want to provide privacy guarantees for.\n",
    "\n",
    "OpenDP assumes the sample size is private information.\n",
    "If you know the dataset size (or any other parameter) is publicly available,\n",
    "then you are free to make use of such information while building your measurement.\n",
    "\n",
    "OpenDP will not assume you truthfully or correctly know the size of the dataset.\n",
    "Moreover, OpenDP cannot respond with an error message if you get the size incorrect;\n",
    "doing so would permit an attack whereby an analyst could repeatedly guess different dataset sizes until the error message went away,\n",
    "thereby leaking the exact dataset size.\n",
    "\n",
    "If we know the dataset size, we can incorporate it into the analysis as below,\n",
    "where we provide `size` as an argument to a DP mean measurement 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"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DP mean: 44.65291940000049\n"
     ]
    }
   ],
   "source": [
    "from opendp.transformations import *\n",
    "from opendp.measurements import make_base_laplace, make_base_discrete_laplace\n",
    "from opendp.mod import enable_features, binary_search_chain\n",
    "\n",
    "enable_features(\"contrib\", \"floating-point\")\n",
    "\n",
    "dp_mean = (\n",
    "    # Clamp age values\n",
    "    make_clamp(bounds=age_bounds) >>\n",
    "    # Resize with the known `size`\n",
    "    make_bounded_resize(size=size, bounds=age_bounds, constant=age_prior) >>\n",
    "    # Aggregate\n",
    "    make_sized_bounded_mean(size=size, bounds=age_bounds) >>\n",
    "    # Noise\n",
    "    make_base_laplace(scale=1.)\n",
    ")\n",
    "\n",
    "print(\"DP mean:\", dp_mean(age))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Providing incorrect dataset size values\n",
    "\n",
    "However, if we provide an incorrect value of `size` we still receive an answer.\n",
    "\n",
    "`make_mean_measurement` is just a convenience constructor for building a mean measurement from a `size` argument."
   ]
  },
  {
   "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 (n=200):  44.59649541816684\n",
      "DP mean (n=1000): 45.723444891570175\n",
      "DP mean (n=2000): 46.28379378485838\n"
     ]
    }
   ],
   "source": [
    "def make_mean_measurement(size):\n",
    "    return make_clamp(age_bounds) >> \\\n",
    "           make_bounded_resize(size=size, bounds=age_bounds, constant=age_prior) >> \\\n",
    "           make_sized_bounded_mean(size=size, bounds=age_bounds) >> \\\n",
    "           make_base_laplace(scale=1.0)\n",
    "\n",
    "lower_n = make_mean_measurement(size=200)(age)\n",
    "real_n = make_mean_measurement(size=1000)(age)\n",
    "higher_n = make_mean_measurement(size=2000)(age)\n",
    "\n",
    "print(\"DP mean (n=200):  {0}\".format(lower_n))\n",
    "print(\"DP mean (n=1000): {0}\".format(real_n))\n",
    "print(\"DP mean (n=2000): {0}\".format(higher_n))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analysis with no provided dataset size\n",
    "If we do not believe we have an accurate estimate for `size` we can instead pay some of our privacy budget\n",
    "to estimate the dataset size.\n",
    "Then we can use that estimate in the rest of the analysis.\n",
    "Here is an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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 count: 999\n",
      "DP mean: 44.80287323701519\n"
     ]
    }
   ],
   "source": [
    "# First, estimate the number of records in the dataset.\n",
    "dp_count = make_count(TIA=float) >> make_base_discrete_laplace(scale=1.)\n",
    "dp_count_release = dp_count(age)\n",
    "print(\"DP count: {0}\".format(dp_count_release))\n",
    "\n",
    "# Then reuse the count to create a dp_mean measurement that resizes the dataset.\n",
    "dp_mean = make_mean_measurement(dp_count_release)\n",
    "dp_mean_release = dp_mean(age)\n",
    "print(\"DP mean: {0}\".format(dp_mean_release))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "There is an interesting trade-off to this approach, that can be demonstrated visually via simulations.\n",
    "However, before we move on to the simulation, let's make a few helper functions for building measurements that consume a specified privacy budget."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "from functools import lru_cache\n",
    "\n",
    "@lru_cache(maxsize=None)\n",
    "def make_count_with(*, epsilon):\n",
    "    counter = make_count(TIA=float)\n",
    "    return binary_search_chain(\n",
    "        lambda s: counter >> make_base_discrete_laplace(scale=s),\n",
    "        d_in=1, d_out=epsilon, \n",
    "        bounds=(0., 10000.))\n",
    "\n",
    "@lru_cache(maxsize=None)\n",
    "def make_mean_with(*, data_size, epsilon):\n",
    "    mean_chain = (\n",
    "        # Clamp age values\n",
    "        make_clamp(bounds=age_bounds) >>\n",
    "        # Resize the dataset to length `data_size`.\n",
    "        #     If there are fewer than `data_size` rows in the data, fill with a constant.\n",
    "        #     If there are more than `data_size` rows in the data, only keep `data_size` rows\n",
    "        make_bounded_resize(size=data_size, bounds=age_bounds, constant=age_prior) >>\n",
    "        # Compute the mean\n",
    "        make_sized_bounded_mean(size=data_size, bounds=age_bounds)\n",
    "    )\n",
    "    return binary_search_chain(\n",
    "        lambda s: mean_chain >> make_base_laplace(scale=s),\n",
    "        d_in=1, d_out=epsilon, \n",
    "        bounds=(0., 10.))\n",
    "\n",
    "@lru_cache(maxsize=None)\n",
    "def make_sum_with(*, epsilon):\n",
    "    bounded_age_sum = (\n",
    "        # Clamp income values\n",
    "        make_clamp(bounds=age_bounds) >>\n",
    "        # These bounds must be identical to the clamp bounds, otherwise chaining will fail\n",
    "        make_bounded_sum(bounds=age_bounds)\n",
    "    )\n",
    "    return binary_search_chain(\n",
    "        lambda s: bounded_age_sum >> make_base_laplace(scale=s),\n",
    "        d_in=1, d_out=epsilon,\n",
    "        bounds=(0., 1000.))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "\n",
    "In this simulation, we are running the same procedure `n_simulations` times. In each iteration, we collect the estimated count and mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Status:\n",
      "0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%\n"
     ]
    }
   ],
   "source": [
    "import random\n",
    "\n",
    "n_simulations = 1000\n",
    "\n",
    "history_count = []\n",
    "history_mean = []\n",
    "\n",
    "print(\"Status:\")\n",
    "for i in range(n_simulations):\n",
    "    if i % 100 == 0:\n",
    "        print(f\"{i / n_simulations:.0%} \", end=\"\")\n",
    "        \n",
    "    # See https://github.com/opendp/opendp/issues/357\n",
    "    random.shuffle(age)\n",
    "    \n",
    "    count_chain = make_count_with(epsilon=0.05)\n",
    "    history_count.append(count_chain(age))\n",
    "    \n",
    "    mean_chain = make_mean_with(data_size=history_count[-1], epsilon=1.)\n",
    "    history_mean.append(mean_chain(age))\n",
    "\n",
    "print(\"100%\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Now we plot our simulation data, with counts on the X axis and means on the Y axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABn30lEQVR4nO2deXjU1bn4P+8smcnMhCyIEECUVUBRNgmKoKCCRVGpWpfaiytWr716u916b1urbW/b25+tdtFKEaWudRcRFWQRFAgSoCCyBEFRCIpkIZkts5zfH2dmMpNMkiFkI5zP8+TJzHc95zsz5z3nXUUphcFgMBgM9bF0dAMMBoPB0DkxAsJgMBgMaTECwmAwGAxpMQLCYDAYDGkxAsJgMBgMabF1dANaixNOOEGdcsopHd0Mg8FgOKYoKSn5WinVI92+LiMgTjnlFNavX9/RzTAYDIZjChH5rLF9RsVkMBgMhrS0uYAQEauIbBSRhbH3IiK/FpGdIrJNRP6jkfNmiUhp7G9WW7fTYDAYDKm0h4rpbmAb0C32/kbgJGCoUioqIifWP0FECoD7gLGAAkpEZIFSqqId2mswGAwG2ngFISJ9gUuAuUmb7wAeUEpFAZRSX6U5dRqwRClVHhMKS4CL27KtBoPBYEilrVVMDwE/BqJJ2wYC14jIehF5S0QGpzmvD/B50vsvYttSEJHZseusP3jwYCs222AwGAxtpmISkUuBr5RSJSJyftIuBxBQSo0VkW8C84CJLbmHUmoOMAdg7NixJuugwdAKFJcGWVQSoKwiQmG+leljnBQNdnR0swwdQFvaICYAl4nIdMAJdBORp9GrgVdix7wKPJHm3H3A+Unv+wIr2qylBoMB0MLhtWI/sya7GdTLxq4DYeYv9wIYIXEc0mYqJqXUvUqpvkqpU4BrgWVKqRuA14DJscPOA3amOf0dYKqI5ItIPjA1ts1gMLQhi0oCzJrsZmgfOzarMLSPnVmT3SwqCXR00wwdQEfEQfwWuFJEtgC/AW4FEJGxIjIXQClVDvwS+DD290Bsm8FgaEPKKiIM6pWqWBjUy0ZZRaSDWmToSNolkloptYKYikgpVYn2bKp/zHpiwiL2fh7aPmEwGNqJwnwruw6EGdrHnti260CYwnxrB7bK0FGYSGqDwZBg+hgn85d72b4vRDii2L4vxPzlXqaPcXZ00wwdQJfJxWQwGI6euCH6uVW+hBfTFUXZxkB9nGIEhMFgSKFosMMIBANgVEwGg8FgaAQjIAwGg8GQFiMgDAaDwZAWIyAMBoPBkBYjIAwGg8GQFiMgDAaDwZAWIyAMBoPBkBYjIAwGg8GQFiMgDAaDwZAWIyAMBoPBkBYjIAwGg8GQFiMgDAaDwZAWIyAMBoPBkBYjIAwGg8GQFiMgDAaDwZAWIyAMBoPBkBYjIAwGg8GQFiMgDAaDwZAWIyAMBoPBkBYjIAwGg8GQFiMgDAaDwZAWIyAMBoPBkBYjIAwGg8GQFiMgDAaDwZAWIyAMBoPBkBYjIAwGg8GQFiMgDAaDwZAWIyAMBoPBkBYjIAwGg8GQFiMgDAaDwZAWIyAMBoPBkBYjIAwGg8GQFiMgDAaDwZAWIyAMBoPBkBZbRzfAYDA0T3FpkEUlAcoqIhTmW5k+xknRYEdHN8vQxWlzASEiVmA9sE8pdamIPAmcB1TFDrlRKbUpzXkRYEvs7V6l1GVt3VaDoTNSXBrktWI/sya7GdTLxq4DYeYv9wIYIWFoU9pjBXE3sA3olrTtR0qpl5o5z6+UGtlmrTIYjhEWlQSYNdnN0D52AIb2sTNrspvnVvmMgDC0KW1qgxCRvsAlwNy2vI/B0JUpq4gwqFfqXG5QLxtlFZEOapHheKGtjdQPAT8GovW2/1pENovIH0WksSmQU0TWi8haEbmiLRtpMHRmCvOt7DoQTtm260CYwnxrB7XIcLzQZgJCRC4FvlJKldTbdS8wFDgLKAD+q5FLnKyUGgtcDzwkIgPT3GN2TIisP3jwYCu23mDoPEwf42T+ci/b94UIRxTb94WYv9zL9DHOjm6aoYvTljaICcBlIjIdcALdRORppdQNsf1BEXkC+GG6k5VS+2L/d4vICmAU8Em9Y+YAcwDGjh2r2qQXBkMHE7czPLfKl/BiuqIo29gfDG1OmwkIpdS96NUCInI+8EOl1A0iUqiUKhMRAa4APqp/rojkAz6lVFBETkALm/9rq7YaDJ2dosEOIxAM7U5HxEE8IyI9AAE2Ad8FEJGxwHeVUrcCw4DHRCSKVoP9Vin1cQe01dAExjffYOjatIuAUEqtAFbEXk9p5Jj1wK2x16uBEe3RNkPLML75BkPXx6TaMLSIZN98m1USvvmLSgId3TSDwdBKGAFhaBHGN99g6PoYAWFoEcY332Do+hgBYWgRxjffYOj6mGyuhhZhfPMNhq6PERCGFtOVfPONy67B0JBmBUQsoO3bwACl1AMi0g/opZRa1+atMxjagbZ22W2p8DFCy9DRZLKCeASdbG8K8ABQDbyMzqVkMBzztGU67ZYKHxNnYugMZCIgipRSo0VkI4BSqkJEstq4XQZDu9Goy255hPuerzqqGXxLhY+pAWHoDGQiIEKxqnAKIJYmo376boPhmCXushsfjAHeWO/HYoHrJrqOagbf0ngRE2di6Axk4ub6J+BV4EQR+TXwPvC/bdoqg6EdSeey+87GANNGOY86Uryl8SImzsTQGWh2BaGUekZESoAL0An2rlBKbWvzlhkM7UQ6l91IFGaMzU45riUz+LjwqW9LuKIou03OywRj/DZkSiZeTAXAV8BzSdvsSqlQWzbMYGhP6rvs3vd8VQO1U0tm8C2NF2nqvKMZ4I3x23AkZGKD2ACcBFSgVxB5wAER+RK4LU3FOIPhmKc1Z/AtjRdJd97RDvDG+G04EjIREEuAl5RS7wCIyFTgSuAJtAtsUds1z2DoGDprpHhTA3x8f1MrC2P8NhwJmQiI8Uqp2+JvlFKLReT/KaVuFxEz5TB0WTpjpHhjA/z+8khGK4t0HlvG+G1ojEy8mMpE5L9E5OTY34+Br2Kur8bd1WBoRxrzbrJZyag+h0myaDgSMllBXA/cB7wWe/8BcC1gBb7VNs0yGAz1KS4N4g9GefD1arrnWLh8XDb5Hgvzl3uJRMhIddRZVWeGzkkmbq5fA9+LvxcRJzBDKfUisKsN22YwHPOk8ziC9LaCpryT4sbpmy/0UFET4fV1AeYt9ZLnFq46x8WikkDGqqPOqDozdE4yyuYaUydNA64DLkIHy73Yhu0yGI550nkczVlcA8DsqZ4UW8EnZWG27A01akNINU7bOftUJ9v3hVK8j9oqbsJw/NKkgBCR89AqpunAOmACOqurrx3aZjAc06TzOMqySeJ1/P+syW4eXljN3ZfmNOp+2pz3kVEdGdqCRgWEiHwB7AUeBX6olKoWkT1GOBgMmZFuUD9UHUUk9bhBvWyEm7EhZOJ91JTqyERPG1pCUyuIl4ArgGuAiIi8Tixhn8FgaJ50g3r3nIaOg3EvpKYEwJEG7iULhFyXEIk2VGuBiZ42NE2jbq5KqXuA/sCDwPnADqCHiHxLRDzt0jqD4RgmnUtpbVhRG1YN3EwnDnM06X5aNNjBFUXZPLfKx51zKnhulS+tCqm4NMiP51cwd4mXYEhx0xQ3VosgAlW+6FElHjQcfzRpg1BKKWA5sFxE7NQZqh8BTmj75hkMnY9M1TXp7ALXnOtqsC05x1JTuZf2l0ewJRLvp2/Xa8V+LBbhnhlaMMxf7uVQdZR7Znh44QN/ok0metqQCRnXpI4l51sILBQR4xph6BA6Wpd+pLmQGrMLZHps/H7jBmcRDCmmjHCwbEuAkf3tvFbsT7lW3Cj+hwXVnNrbTsnuWlRMmDy1wseh6rq4VhM9bciEjAVEMkopf2s3xGBojs6QifRIk901J9Ca2//Sah9Wi7CoJED3HAu5bgs3TvHw3Cpfg/vuL4/w7EovSsGP5lciArdc4OahhTUEaqNYBNbsCCaC64wLrKE5WiQgDIaWcjQrgM6QifRIkt01J9Ay2V/pVRTELH5KKf75vo+rz3El2hG/b3FpEKsFzh3mICdbeHKZj+wsYfu+EN1zLNSGFS6HYt5SL70LjAusITMaNVKLyFOx/3e3X3MMXZn4gHjdRBePzM7nuokuXiv2U1wazOj8zpCJ9Egqvb282odS8IcF1fzyxcNU+aIpxuFkgZfOePzyah8WgQvOcFKYb+GCM5yIwAsf+BLtiN93UYmugLf8oyD5HitKgUUUb28MEokqrjnXxe9n5SMC91+ba4SDISOaWkGMEZHewM0i8g90LYgESqnyNm2ZoctxtCuAjspEmrzqyXMJcxbXMHtqXcqLQ9VR8txCcWkwJTVGhVdxzww3p/a2J1YHl53lTAi0uEroQGU0sZoaMyArsb/Cq/A4hRdX+/E4hbc2BOjT3cKOfRFqAhEeeqOaScP1/coqIvzs6m70LrDy3CofCrDbBEHx+1n5AGzfFzJ2B8MR0ZSA+BuwFBgAlJAqIFRsu8GQMUe7AmjLMpyN0Vi6jL+9XY03qOMabr7ATb7HwpzFNby82kelT2G1gMcpWC2Ssjp47J0arBa47RE9v+rXw8rPv5Wb6Mv+cq16i6+qpo10kOu2JATRjn3xwDkLowZksa60luLSYNqVTaVX4XYK4YgyqTcMLaJRAaGU+hPwJxF5VCl1Rzu2ydBFOdoVQEekk4gbif+woDoxy5891cPDC6v5weV1qTGKS4NYBCwW4ZHZedzxWAXYSaw2BvWysXVvLTUBxfQxTjburqVfDysflobonuNnxthsJp/u4JW1fm66QKuZ3E54d3OQWy/y8Kvrc/mvf1RS7VfkuoUHrssDYFhfO8+t8jF9jJM5i2sQgZsvcIOCxxZ78QYVdz5WQaGxOxhaQCbZXO8QkTOBibFNK5VSm9u2WYauSGusADLJRNparrBxI/E9l7o57I/y+roAc5foDKr1U2MsKglw0xQ3D79Zg80q9C6wMmFoFu/+K5AQaBYBtxNmFrl4a0OAn38rl+45fhZvCrCoJIDNApEoiZiHmy9w8eJqP08u81JeE024rF51jitx3/gKrGiwg5dW+whH4KE3ahL7PU7BZtV2B4PhSGm2YJCI/AfwDHBi7O8ZEfle02cZDA3JNBr4aDhaQ3gycdfSnWVhFnwY4MYpbu651EM0Fk7wxvo6b++yiggIKakxlm0JUOlV/Ozqbnz/shwiUbhmgh7c46upGWOzCUfghG4WZo7PpjDfwnUTXVgtsGF3CJtF529SCkTA40x16U1egVV5FaDIdQn3XOrhnhkeLKKo9KoW9d9gyMTN9VagSCnlBRCR3wFrgD+3ZcMMXZO2rkXQmq6wZRURbpriYv5yH98cn82gXjZ27A9x2K8Y2d/OOxsDDOtrZ1AvGwUeC/OWevlWTAAUDXawvzzC4k0B7pxTQWG+lXy3kO9Jza00+XQHNgtMPt3B8o+CjOhnZ967NUSisGlPCKcdbpzswhtUvLTGTziijc3pVmBWK1hEuPUiT6L/00Zl8/JaP4tKAka9ZDhiMhEQAiRbESPU82gyHFt0dDRyW9KcITzTvsfjCuYt1cmLl/wrwEtr/BR4LOS7hdunerjjMb0K2l+u1UdRBa+u9RONKvI9VtaV1nLjFHfi+s+u9PLQG9VEotrI3Le7lZfX+IkqeHG1H5cD3t8WJJJUyDcQgieX6za4nYCSRm0w4Qgc9isiUZUwTC/bEiQaxaTVMLSITATEE0CxiLwae38F8HibtcjQpnSGaOS2pClDeKZ9jx83daST4p21hMKKaFTxjdFO1pXWcuU5LnYdCNO7QAuY+DUraqK8vs7PvKU+3A7BboPH3/WyqCRAnkvYsT9MJAo2C5RVRPmyMkqWDdxOC985z8Xcd72EIgqnXQsGEUCB0w55HgvXnOvioTdquG5SeqHWu8DKYV+Up1b4KK/RrrNFQ7Io3lmLw27mdIYjp1kbhFLqD8BNQHns7yal1ENt3C5DG9FccNaxTroMqvGsqJn2PX7czCIXM8dnY7XAYT+8szHAZWdlk+uypFxz3OAsnlvl44llXhx2YWR/O4GQ4pYLPTwyO5+R/e18/EWYsQPtPHp7PnfPyOGEbhZcDsHt1LaHp1f6qAmkZuGLG6X9IThQGcVqEbrnWBr9rKaP0YF04Yji7ks8fGtCNu9vC1IbVomssAbDkZBRqg2l1AZgQxu3xdAOdIZo5DhHq+pq6vx0apjH3/Um+p58rlKkBLklP6O4zWTNjgDzlmohEL8m6GC3/eURuudYGNrbyo4y/R5gwye1DO1jZ9OeEAAln4RYt0vbI8YNzmJRSYCagErcd+4SvZoJhGBkf32ezQKhKKC0y+zV57h4Ypk37bN4ebWPar+WKg8t1J5M+W7hygmuLrE6NLQ/bZ6LKVbPej2wTyl1qYg8CZwHVMUOuVEptSnNebOAn8be/kopNb+t23o80B7RyJkM/Eer6mru/HTXiPe9yhdNnBuJKp5a4UvJjJruGeV7rPQusCbcReP3t1rgiiInn5RF2PRpiKLBdopLQ1gE3tsa5MvKOoERjsLNF7jI91h5MjbId8+xJJ5XMr3yLbAnJhxi1Ph1LiXQgXZ5bkm4vP7zfR8icM8MDzu+CPPOpgBRU97LcJSIUm37LRKR7wNjgW5JAmKhUuqlJs4pQAuVseio7RJgjFKqorFzxo4dq9avX9+qbe+KNDawtpa7aabXv+/5Kq6b6EoZhLfvC/HcKl9an/36QscXjHLLhZ6Mz09uW6BW4bAL5dVRLBaYNNzB6IFZiXPr9+GN9X7e2RggEtV6/lN721i1LUgkor+cHif4amHCqXY+/iJCRU2UqIKiIXaKd4awWmBQoYXS/VFysrUK6sPSWny12tvDmQW3T/Mwd4m3gZopGWssTuKswTZ2fBEhFAGbFaJRhT8EKoqO0YjCxaOcvLUhgK9WEY1tv+ocs5IwNERESpRSY9Pty2gFISInA4OVUu/GakHYlFLVGZzXF7gE+DXw/SNo8zRgSTzfk4gsAS4GnjuCaxwxXdm7J05bRyNn6mZ6tFlRH3y9moqaaEbnJ+MPRvEGoSagyHXBhGHa8Ny/Z925Kc+oPILFAiNOtidWA/EVwV9uzeO+fx7GF4wSjcKqbaGUe23ard9HorBjX5STe1j5ujrCe1trE26AbqfgCyoefqOGXLcgQLZDb7OIFgrhqLZHxL2bZl/Uje37Qjy5TFeN88ZCHOLtsSjFV1WRhLC551IPT73na1A/wmBojkwC5W5D16d+LLapL/Bahtd/CPgxEK23/dcisllE/igi6b6tfYDPk95/EdtWv22zRWS9iKw/ePBghk1KT2sGWHV2igY7uP/aXObcUdDqmT0zGfjjLqR3PlbBfc9XJZ5xY6qudMbl7jkWXl+XWpakKVVZ/PN1Zlko8AhXn5ON3Wahd4GVWZPdvL7On3Ju/BkVFliZNsrJF4ciXD/JTWG+hfNOy9Lt2hhg5vhssmx6uHfYoFu2NiSffIKVYCw1kkXAboXPDkbwxjRJ3xjjpHeBFW9Q8R+XeigssPL7WfkM62vDF6xbRZzWz079Rf5tj5Tzh9erOVQdTQiB7jkW9hyMcKhalxV9b2stFtErh1P72CmviXYpZwRD+9CsgAD+HZgAHAZQSpWiI6qbREQuBb5SSpXU23UvMBQ4CygA/utIGpyMUmqOUmqsUmpsjx49WnoZoOt797QXzaXDTnYhLcgRJgzN4tW1fl4t9qXUYE4mndC5fJyTQ9XRRms41yf++ZbXRBlxsp1X1vr5+nCUuUu8PBgbbE/tnXqP4tIg+8sjLCoJoJSu6XygMsrVZ2u9/zsbA+S6LMyMGayDYQhFFDPGOqnw1s2JoqpONQRaYKz6OMj0MU79XJTu47MrvezYH045b/On2lBdn/qKqPKaKA++rhf18cpxcZXUY4t1gsA/LKhmf3mkS056DG1DJiqmoFKqVkTPkkTERqNVcVOYAFwmItMBJ9BNRJ5WSt0Qv66IPAH8MM25+4Dzk973BVZkcM8W05m8e45lmsu3lCyIexdYWVSis5Qu3hRICSpLVvdZLTqtxcyiuhxE+R4dmZypqiz++ea6hE17QozoZ2fzZyGiKjbDt8GWvaGEN1NxaZB/vu9LDLKgDcF5LuG9j4MU5ls4UBFNBMoBDOtrpbQskghss1shFIErz3biDSgWb9IDc5ZNx0IUDXbwSVmYhxfWoIDlH9UN3EKdeqk5Y3NBjlBenf6gE3MtbNoTYvoYJ0N62xoY5A2GpshEQLwnIv8NZIvIRcCdwBvNnaSUuhe9WkBEzgd+qJS6QUQKlVJloiXOFcBHaU5/B/hfEcmPvZ8av1ZrEx+IlIKfPlvFzPF1g4yp23vkNGfjSOdCGo4o7pxTkSIc6huJ396gV3IzxmYnhM6VR2B0ja9sBK3v/PxQhJlFTlZsrSUUUViFFFvJS7FiPdNGOXl7Q4DT+9nYuCeEr1anvHDY9Czpy0otHCwC/Xva+f5l2kB+x2PlnDM0i/e21rJ6ey0HKqPkOIXqgE7BXV6tWLMjyPvbgg1mW0479Myz8tnBzCYnycIhWaCFIvDJgVgE+c5aFpUEyHcLZw91mNQbhozIRED8BLgF2ALcDiwC5h7FPZ8RkR7oSdIm4LsAIjIW+K5S6lalVLmI/BL4MHbOA21RoCh5IKqoifLyGh8vfOBLpEow+fNbRlP5ljJxs61v6J5Z5EqoehaVBLBZYeKwzHM6FZcG8QfrVDCThmex8uNaXl6rB8yrzs7miWW+lBWjzuLq4bR+WfiDipUf16XAsFshGLNHux1CVCkG9bKnCLFwBFZurcXjFK6Z4MJqFZ5cVgMKPA7hUHWdy2rPXAtfVkXJtkMgrOMgPjsYSURUHwlWgajAWQPtrNtVd3IwpJJcbGsaXXEYDMlkku47Cvwd+HvM/bSvOkLfWKXUCmIqIqXUlEaOWY9ODBh/Pw+YdyT3OVLqD0QWi86lM2+pz9TtbSMySfldX91XXBrk86/1wP3o7fmJc5KD2xojPgm4+UIP63cFeW9rLSs/rgVg5Cl2/n16TqLSWn1B9cyfN7F13gdEymuw5nvwXHY2rnFDCEV0Gozpo53MLHIl3GsvHq1XG3G7VZZd15F+aGENbieEw3rV8dnX0YT6CeDLKi15/PWEwZEKB4DaCAzto/viytLut0P7WDnsh7NP1faZKSOcvFrsb+ZKBkMGAkJEVgCXxY4tAb4SkdVKqf9s47a1OfUHoqLBDsYMyOLOORUmf36M1nb9zcTNtv4qY1FJgCkjHHywvTbFgSCTDK3xVBjz3q1JlPAMRxQi8K9PQ8x9t5pdZRGKhmQlBNWzK7341u1k87PLIaSNxpGKGqqeXQ6Ae9wQlIKlmwP0LrAmyoRePMpBVMEPLs9JrEgVMKyvhe37IglvpGThANr7KRim1dixL4KiLqPm519H8Aa191P3HAvBkCJsTGuGDMhExZSrlDosIrcC/1BK3SciXaJgUEfVOD5WaKvEfs2l/K6/yigrj7BsS4CZ4xsWymmq7fHCOwcqIrgcOsrYahEee6eamphzWvFOPU3ftCfEFUXZfFIW5r2tQaoXrEkIhwShMNUL1vDfvxzNvHe9VAcUL632keuyUJhv5YUP/AlvocJ8Kx6nsK88ymGfbmcs916KcAA9628N4hll48v7b4zRK5pQRNsm4sbuYEjHWMx+tPyIhX5LJgzHQ3xRVyUTAWETkULgW8D/tHF72pWOqHF8LNGatRWOhPqrDKsVTupu46XVvkS+IqddD7jpBrlkwfbwG9VMGKaNxRt311JaFkkIh5svcPHEMh9z7ihI3PvJZV7OGmzntYoa0hGtqOGvi2pw2AWltK3iL4uqiUbrBv48F9gtis++jjJukLYF+NbtpHrBGqIVNVjyPeTE1FVAgziHllLf2ymu6gqFtbCYMTabN9b7eWtDALdD+P2svMR3/pOyMDv2h1s9PUpXzx7c1clEQDyA9ip6Xyn1oYgMAErbtlntQ0fUOD6W6EjX3+RVxrMrvaz4KEiWDQo8QnmNIhDSs+IbJ7vJ91gSgw7oQT4S0Z9rOAof7Q3jtGuj8VXnZLNsS5BD1VGeXOZDKZ32Iz4ghiOw/YsQlnwP0TRCwprvIRSBSNJoHDdYO+1wyRgnSzYH+ezrKL3zLZTs1sKhKkldFU1SV8WFRGsRX6WktNkKq7cHeWtDAKsFpo108M6mYEJdN25wFu9sDHDPjJwmB/GWTBg6apJhaB0yMVK/CLyY9H43cGVbNqo9aesKZ8cynUUFt2N/mCy7Hl/LaxRWCwwutLJzf4RX1vr4/ax8Zk1283isEltcvx6PT/A4hArtQMSyLUH6FlipqImSnaVzGV030ZUYEG1WqPJBzmVnpwzqANhtdLvibHzrduJ9Yw3h8hqsBR5yZpxN9rghXD/JxYIPA0wb6eTF1X72V2jjc1PqqqYERFOrjsaICwe7Vcdb+EPas+mwT+8Jx8p9JX+GG3fXEo3qQMBfvniYsooIBR4LL61ueXqUoznH0HnIxEjtRLu5noYOeANAKXVzG7bL0AYcqS64s6jg4im5426nsx8t565v5HDX3EoqvYrZj5bTK89CpVeRnaV18TOLnJzUw8Yjb9XwWcwDKjdmwvjXZyE8TuFbE3Tq7ORZbc9cC/vKo4mBuP4AHY2SIjgi5TVUzl9C1UureP9HU7ji+tN5abUv0XaLkHYlAo1vB45q1eGwgcOua1gA1CqtxvrLrXn89zNVvLMxyPmnJQ/8UfLckqIK2rEvxEMLaxLR3Y0FLO46ECbPJdz3fFXa71VnmWQYWkYmKqangO3oBHoPAN8GtrVlowytT0t0wZ1FBVeYb9WrAal7/9RKb0Kd8sjsfJZtCfDiaj9Wi3DxaAfvfVzLrMl2bp/m5s9v6n5W+cBm1bP6uP49PlDFjeGWpLQWrnFDGgzGX/50fsPVAKC8Ad779WK2fh7CPrrunKiiUXWVJd/T6CqhpasO0B5RklRALm7juGtuJS6HfmjLPwqyfV+IUQOydMR2BG65sE4VZLUKHqew8uNgQvUUt1+s3BqkJqDonmPBG4hitwnXTXSl/V61xiTDGLk7jkwExCCl1NUicrlSar6IPAusauuGGVqXluqCO4MKbvoYJ/Pe9fL3xTXcdpGHM06x8faGIHYLuLOFXQfCiTQVNQHFjLHZ9C6wpqTBAOieI0wZ4eSVtX6eXFbDJ19GEgPVrgNhrFaIRHSCuyqfIscpHPanavSbmvUTClP+yhp6xgRE8uCfjqg3QNXTSxOhz8mrhKZWHWX//tdmVU7JMRSCjsmIJxWcNjqLtzcEKauI8tXGAKf2tvHxF+GUWtbzl3uxWXXT4t+Z3gVWnHbw16qEAAqGYcoZjka/V0c7yTBG7o4lEwER/6pVisjpwAEySNZn6Fwcy7rgeM6i97YGE5XSACIxL6LnVvm4oiibl1b7qPQqlm0JMGWEk1yXhT8uqNbG7RxrokRoJArFpTo/0ZgBWYlEf+EIFOZb6JlrZdOnoQbCwSLg7JFD4GDjme7jA3vOjl0ceG45qraJAId0+0Jhqp56t9lnciQqJ0V8MaLwOC2cdlIWmz8Nc/+1ubxa7GPxJu3t9Kc3a1BRKIwFiT6+xEthviXFZbjAIwRqFXPu1J5fsx8pZ+Pu2hS1U/3v1dFMMoyRu2PJJJvrnFhOpJ8BC4CPgf9r01YZWo3i0iD3PV+VyDWVnMnzWNIFXz/Jzc0XuuldYEVEG5MvHu3k1ovcADz+rjdhnH55jZ87HqvgiaU1RJU+Nq6WuP/aXB69PV/nedkT4s45FQkB07vASs88K5s/CzGwV8Pnku0Az2XjdWa/RrDkewD45OkPmhYOTZFpKbiYyqk57FbIzdbeVoeqdcqR/eURbnuknLc3BAhH4OYL3OQ4hRyXcPEoB7kuCxYL9My1JtLgA0w+3YHFQuJ71CvfQllFajb/1vxeHcsTm65AJl5M8bxL7wED2rY5htaks+Saai0dcnwmGq+/vKgkgEVg6kgH35qQzbylXpx2nTLFF9QeTxaBQb3sLCoJ8Pi7up70yP52CmPlQ59d6eW9rcFEfMX+8ggn5gqfftVwAAqHwT56CLlhqHpxpb5JMnYb+TPPxiLaeN0eRCtq8K3b2eQqIhSBijq7eUpCv1P7aG8wiwWunuBKSTUzabiDlR8H+eb4bAb1stE9x8K7m4NMG+VMxFjEa2D/5KlKLh+XnXA5bq3vlTFydyzNlhwVkZ7A/wK9lVLfEJHhwNlKqcfbo4GZYkqONqR+Wc/i0iCvrvVzqDpK74L2Mfa1donT5Os9vqSaqNJ2gnhJzVyXJaXsaFwAXHl2NucNd/Dex0FeXuPnvJgXz4qtQVxZega9vyLCy2v0wJdsP7AVeHDPaKjvT2dgdo8bwqTTHLz4zTntJiSw28i9fnLGMRXJSQAFXdWuJqC0zSFWRvXWi3Tq9dseKdepzSuj5Ll0OdNbLnTz8Bs1dHMJChjV385He8Mcqo42Wtq0pZOEti6Razj6kqNPAk9QF0W9E/gn0KkExFHh80JFmop0J/QChxO81VB5qOH+E3uDPQtqDkNVmmSzPfuAzQ7VlXC4suH+wpPAYoWqCqipari/98naHaXykG5DMiJ6P0D5QfB7U/dbrJRV6B8Vh76CgI8iF4w5H3754mHuP78AesXyTX19AIL1CiPZsqBnb/36YBnU1pstZzmgR6F+/eV+CNem7ndmQ/eeLCoJcOuISgZSCQd0pajbTgvz6gdC0eCB+tj9e3VB5WRcHsg/Qb/e92li84fv1XDbGCcD3BGqfEJhnsJT/gXWWsg+mE1AKdgX4Ad/ysFzYgGBQIirT/ma7evDfLAsSo9uFs4vsLBhs4tKSzfsKsylvcs5w+7kjBNhszrM1x9+Ss3Lm4jWCHZrlMLoV/DmQjyMw3XGKSgFFfZ8nBOGkDf2ZE4I6u+Of8uneH+3hNWHfWRFc/AD2fYwPdz1ng3wVY2DQNiKyx7mhDT7v6x2EIxYcWeF6e5K3X/Qm4U/lPTTDYWpenFlRgLMNW5IiulDQaIq3ZiBdrZ9HiKiSNSM6F1gbTDJeGqFL3FeNAqlZRFmjs9OCOfWiMCO01k86Y5XMhEQJyilXhCRewGUUmER6VoKwI9L4G//23D7vX+EgcNg42p48o8N99//N+hzCqxdBs//reH+3/0Dup8I770Fr81vuP/hF8GdA0tehbdfaLj/bwvBZoM3n4PlC1P32ezwt1hZjlefhDVLU/d7cikcP0cvz1/7m+4D+gO/H+ChXvDbJ/Wx8x+GbRtTz+/bn+JvP6SNhGt+zQDvrtT9A4fDvX/Qr//6Cz3IJ3PaGPjPX1NWEWHA0z9LEcADgPPzxqG9poHffh989WbbE6bCTbEy5g/8u57aAncBrIUDYy7FYrmeiUOsXLjsx/q4WGWRkUD0G9ewc/wNPPbSPi58/x4uTL06Wb2v5e1el5Fbe4gLX/5PeFlv/3HsIT160gDe3FZI31w/f565KXbWem2BA+ad/F22uCfxvZFfMuDJn5B4uBfpl79+dyhrDndn2InVPHDxx9Tnp2+dxqb9eYzqU8m9F+xosP8HC85gx8Eczjn5EHdPSn32NUErP154Bnsr6wzD+IIpqqam4ihyzx5CNM0veMMnIWxWcNolYQiu76aa67JQG9YCJc9tobwmSjCk+Of7Pq4+x5XWNnC0hubO4El3vJKJgPCKSHdiQZoiMh5IM909hhkwDL773w2394yVwT71zPT782NlTkecBXkFDfd7uun/oydArwYltfXqBKDofDhlUMP9caf8c6fBqWek7pMk/4LJl8GZRan7rXame/SP+7ujL6fvWedzoDLCio+CjBucxeCTPXXHXnodnPeNlNO3l2clZn0nn34z+/ZX1Z1baAN3t7qDvzUbAr6U8+mmaz0V5lv54hvf5aRudQPHvvIIG3Z7GBPfcNP3IVLPoNu9Z0ItUXjyXeS7hdEDstiwu5ZzhzlY8Gk+0850smxbgB4X/5DtX4QJhBSBkGLKCCeDR/VnaB87rjw3zwy/h3GDs1hXqmfiw0+y8Wl5IUSgNjuXf138Q7bvCzN6QBbLtgQo//s77D6kjd9f1Tj4zdJT6x67zUr2uadR5hlGTRCe3lnAsDPuYffflxOtqVuF7Tyon+/ucnfK+XE+LdeD+46DOWn37z+svxuby3JT9lstiuE9D/NFldbx9/AEOFjjACTh/dRUHEXV/CVUL1iTWE2c0A2+PhzbHdEV7Py1imdXeimriLKoJMDXh6M8vLCacEQb/OPOAIMKrdw4OZddB8LMXVLDCx/40toGjKH52CUTG8Ro4M/A6eg5Wg/gKqVUp8roamwQ6Wmp7re+/QJI1D1oLhV68j1zXUIkCrOnejLSISe7VFotuqJbchW5Ef3sbNkb4uvDUf5yax6LNgZ4Z2MgYXQF+PuddcJ6zY4A85b66J5j4TvnuXh4YQ0FOcLM8S7W76pl055Qom5CPBvqlz+d33S8g8vByX+8NUVVU/bvf23y+AYG7VbAnRXm71eXsKfczdzi/uwpd2d+chN2i8J8C/16WCneqV2B40n+3t4Q4OLR2kB91mA7JbtCXDxa7393s5+X1wQStov6JWOnjnSmuMJm+l0ytD1HZYNQSm0QkfOAU9E2rR1KqRaUMjF0BC1dnrd01pdO3zxncQ3z3q2h0qea1CEnn/vsSi/nDnOw/KMgvQusFA12MG5wFos3BRIz2LvmVtK7wMpNF+hB6SdPVTa4Zrx29aHqKA8trMFmhaIhjkTtjzseq8AfU/HHvUv7XnsOe+cuSxsxDYAvSOXqnRRMGILbaeFQdbTRaGlcjob5vVuJQMjCsxv6cf3ovTx8xSaW7DyRp0pOptKf1fzJSS6y9e0U5ROGUFYRxWmHNdv181+8KUBUkaia92FpCLcT3koqkBQn05KxJmty56dRASEi32xk1xARQSn1Shu1ydAJaMq9sKlVSTp98+ypHp5b5ePKc5wJd9NFJYEGq5nkcw9URHl/W5CvD0d5cplOR735s1oiEV1VLj7ojOxvTwS7xXXj2/eFUlYrI/tnsWpbkHAEPE5h1cdBhvW1E4nqdBGV3ijTRjnZtCfEyP52FjGYEXcLW/6wuNGYhOoFa5h1z5lMGeHkp89W4b9cJ/dLjn2QLBuIoBoTNEdJRFlYuK2Q5Z/04LpRn3PJsDIm9v+aO14ezSFf85OChF2inp1CBLLPGsJtF7n58yIv/3zfRzgCE4fZWbVNzw0FXSp195dRiobY2fZ5mGjMuB2OKG650JNSMhZg8aYAb20ItJqh2aTgaHuaWkHMaGKfAoyA6MI0lkNnRD97kx4pja48yiPNerLEzy0uDWKxwLnDHJw33MH35lay8uMgYwfZ+dwexWaVtIPONefqbckeL3GV1NSRTop31jJlhIO3Nvh57J0aagKKfLdWgcXTc8S9d0LDB5H7nShV85ekfT7Rihre3xbkpTV+lNIDanYWfPli6my8sfMTxPVaTWG3kT1+KMGPPk1c23H6Kfg/2ApRhbfWxtzi/iza1ouzTzmUEA79C7zsKXdRV1suDWnsFIdfX8PAacM4UBnFatG5nCwCq7aFsAicfrKNLZ+F2f2l1ut9uCuE2yFcO1G7GT/4enWD78CMsdm8tSGQUnvjaDApONqHRgWEUuqm9myIoXPRmHthcx4pja08rFaa9WSJn7uoJMC0Uc5EfqV4tO76XSFuuqBOz97YoJM8QNz3fBXjBmexaU+I8uoorxb7Y8ZWldCX3/d8FbsOhBPnzV3ipcqvjb2HX1qF8tZzAQasBR4OVEQpLLDiC0YJR6DmzCH0GzOEYLjOxbQxLPkeev5qVgNvI71TkGwHyhuol3PpvJRrZA0oTBFA+w9n8/LmvgD0zAnwh8v+xa6vPfy9uD87D+Y02pb6RCtqGNnfzstr/ESVdmfNzgJ/LfQ/0crWveGUIkfRKFxzritWT0PvaOvgNpOCo33IJN33MREoZ2h90tkvHn/X26RtorGVRzhCk+cVlwbxB3UaCICpZ2qbwytr/USjYLNoD5vk9mQy6OwvjxAMKW6cUteeJ5d5OVQdTVwruc1jBmQx3+olFNEpKgq/PZED85cTDaaqjrp/82wKcoQJQ7N4a0MgEUsQDEN4w04ON5WHyW4j57KzCW3YiXfhWi0cYiuJuEDIPWdIoyaQOK5xQxpdoRyscfDI6oH825jP+MNlm1mx6wSeXH8KX3ubHzwt+Z4GdgV/LWRZYc9XkcSCx+0Ab1AnN4w/y10HwuS7pc3TxBvPqPbBBMoZUmhOr9tc6oOmVh5N2TReK/Zz84UeKmoiPLncx5PLfeS5JWGAfrXYxzsbAw3sC80NOjYrTBmRmm10yggHr8ZUSfH+xl0545HEdqtOBqjOHELe9eBbuJbAwWqsBR6uuu98VtlPZuZ4F4tK6oRDnEOvrGlUOMQFAED500mCJ6oSgsM1rnnhkHy9dMbxqBKW7OzJ+3u6c9UZ+5h5+j7O6lfBTc+PxRdq4mcfa0ODzVbIsgs2pbBZhcM+xekn2yneGSIcISUL7JXnNFT1tXZwm0nB0T6YQDlDgkz0upnk92/Mc6qx81LVBXYsFuGFD3xYLZIwQK8rrWXScMcRDzqRCCzbEqBfD1vivsu2BIhEGu/v14e1yuimC1z8Y4UPGTuE3pOHcqha69xHz/CweZk3SZDWJaybODyLF5pwke35q1lAzJU22FD/H6/3IJJZreq0le+S8IdsPFVyMm/v6MlpPQ8nhMOYvhVs3JdHVNXZJ8TtpNtVE8kZP4RIlIT7b/cc4VC1IhRRuB0QqNUNK94ZwipaBXXHYxXYrDBxWN1n35aqns5SzKqrYwLlDAky0eu2NPVBU+fVV1sVDXYQjSrmLfVx55yKlGOLS4M8u9LL/vIIc5d4mfeul/NOc3D9pPQxAIUFOjlf8n2LhjjYtCfUaH//uKAaBfzzfR8XneHg/e1BvIEogs5w8tg7NdisdUVylm0JJPT1y57+qFHDczzTKzRfZS4T4QA0qHyHywH+YIPC1AdrnKyo0cF3A7vXcP+0j9lbkc0TH57Ch5/nA4Jk2XCNG4LHCROGOVm9PYivVnGouu5itWHttetxCmMG2vlwVy2BWrhxSl1t8OLSYJvbAUwKjvYhEwHxfXSa74Ei8gGxQLk2bZWhQ8hUr9vS2IrGzkunLsj3WOkdy7gap7g0yNMrvATDcOV4J727W3n8XS8rtmpjdjohMX2Ms9Fkb43ZU+KDs68W3toY1LPmeOSP0mmzLx3rTFxv+UdBLh7t5KVHN+vZfDqvpCTVjVW0kTtdMr9kIZIp9SvfpTV8J/HJITe/encoN531KfdN3cbm/d2Y9+Ep7Ppa76/ywdLNgQZCBupCOnKyheKdtdSGdXf/scLLjVPc7WooNik42p5m60EopTag3SfOAW5H16bO3CXCcMwQH6iTaQ+9blxdsH1fiHBEJQr4TB/jTDluUUmAUASuOjubi0e7OONkB3dcnEOOU1i1LX2kctFgBwO/3MN/j/wrV9h+xX+P/CueHbtYVBJovEZGgTVlbPQGwe0Qbr7ADbHFwfvbgim1JGaMzU6f3gLAIilRy1l28Mw4W8dJJCFZ6fX/R4pr3BByr5+cEDaWfA/iTn6WwtrPunPny6N4dPUATs738cC0j3Hk1LUnGNIG9zg5SadbRNexDoRgRD8bf7k1j3BEx0BU1ESNobgL0VSgnBX4FtAHeEsptVVELgXmANnAqPZpoqG96Ci9bqbqgrKKCErBecPrtg/qZWtQ+S2Z5c9s4c2fvk0ooEe72q+rWft/S7j2QRsXXzSUl9b4efEDH9EoCRXJiH52vqyM8M3x2UwZ4UxEg7/wgQ8VK0A0akBWg9QRjabniKqEXl/QtRjqq4aaKyF6pNRfVaRLBRJRFt7cVsiyXT3oX+Aj6IviX7edy0/exfvdz8dnq1vNVCc5NSUvkP71WZi75lYiQKA2yrylulRpOjVTZw9s6+zt6wiaUjE9DpwErAP+LCL7gTHAvUqp19qhbYZ2piP1upmoCwrzrXxZGeGplV4+/1rPVAs8Fpx27QKbjrk/XpYQDnFUbZhFv3mPZ+4aicUisSI5XnoX1HlcxeMw+vWwUVETJRpLYnfzBS4OVEbTpo5w9czB92XDcqSWfA/ZWbrmgqKu0mj9QbwtaTQVCNqQ/fGX3YAofT54myttHzL9wOu8Unsur74RpfaQL60As1sgFNUFiKJRveoAnXcpHnAIeuVXVq6LEk0b5eRnV3frdIFtJvAuPY0m6xORj4AzlFJREXGia1EPVEqlKYzQ8ZhkfV2f4tIg85fpGIVT+1j5siJCZSyJbJ8CC7+4Nq/BOZfKL9NfTGBh9GeAdtG8c05FIuBu9qPl3DTFxStr/VR69e/D7QRvoC4RYLyWcyRKYrbpW7eTP926MFUg2W2c8G+T6VY0BH8jGcxGnGzjo3rBZ61Nc3aJZM76/aXM/HguZ1g+4cBhB/PXn8yqPScAQvbE08m7ti5gL9keb4l5XhUWWCkrjyQCuJXSQmTsIDufHIgkJh3b94V4/N0aXA5Lh8/ajyY55bFOU8n6mrJB1Cqlq7gopQLA7s4qHAzHB0WDHbidgs0KO/Zp4WARGNjLSllFlB/Nr2D2o+Xc93ydXcFakN7oa00yBte3s+S5hJfW+LnlQg+P3q7TlseDqe94rJxnV3qZMTabSKy6GugAwhXWflzyq4vJ7a1Todu653DirMl4xqUKh/qJL/oWWNtUOECdXcLWPWY+tKRPv2Et8PCF62R++togfvb2cPwhKzeM2YtVdAP9qz7Ct25n4vjsWF7A7013E1Xarn3dRBezJrtQClxZ2nYTicKusgjjBmclgvAqaqJUehXXTXTxyOx8rpvo4rVif4pNqL0wgXfpaUrFNFRE4im9Be3FtDn2Wimlzmj8VIOhbajyKXrlWbh+Up176podAT45oOMmHpmdl6IeOOm6Ceybt5SQPzUSuvBb56QEdyXbWRR1g/jz76fWuXBlCSu2Bqmo0SU4G6ol+vMfK4YD8PJqHxVelWLs7ZYt1IZVnVcU2lOqJTjsdWqdTHCNG8IzT47ntkfK068o7DY8M3RN7WhFDRsr8rl7fx7dXbVElAWnLcLtZ+/m1eUQGjeEfidYqfRGAcXfl3ixiLbjDO1j5yfL9PO/7SI3L6z207vAyoShWby/LciByijFpUGejB3z4OvV5LuFK89xdVi6DBN4l56mBMSwdmuFwZAhhflafZE823t9XYACj1BeoxP5JcdvfOfuM3kioqh5Yy1V+w+T27sbzulFdBs/pEGMRZwqn+KmKfr8/eV6BpmdJdoGcaGbv71dw6ZPQ3qmhOLB16sTtZjj973/2txE3MZLq30JVVU4qrDbhGBIpfMiTZBlhdqkyWs8BiPZQGy16HiE+pHcTbF4k59sO9CEkTyqdNCc8gaIKuFgLD3HwO41TDjlEFMGfcW7z5Xyz91DCV40hdzxQxICb8ZYrTo6VB3FIvDwwhoUMPl0B8u2BDhUreuHP/Oel6iCaSOzGHZSFvOWemNV6bI7ZNZuAu/S02zBoGMFY4M4PiguDfLEUi9jB9nZ/kWIqtgE32EDl0P4v1laJbRmR5B5S72IQK5LEEjUo6iv567vveILRhPpqm97pJzvTXfzyZeRhM0hy1Y3cy/wCCd1t7Jlr053neuCwz6Yc2dBiuGzokanLVfAxSMdfLirlq+rm//taSEEV5+TTa7bwrx3vTp7bCzKOdOIa4Bhfa3s2BdJETITh9v5155wiieYb91OXZ0uTTxHN2eIb535OZcMO4BS8MaOPrx+6k04xg1LtCXPrVNx9O9pZfeBCAotzLp7hK8Op16ze44lUc86vqJw2KVD9P7HqxdTUzYIIyAMxxx/XHCYj7/QqpFkI6ndCrOmaJvAix/4sFmFX12fm5KqfMf+cMoAAKQtcAS6Ct6Dr1czcZidbfsiidnk3CV6ICvwCIN72/iwNMRZg+zsOhAhFFZ4g4qbLnCzqCSQYvh8dqU3kaE2HXFh0BhWCymV8wDOGmzjw9KjrzeRfO9mK+oBJ3oCfHv0Xgpctfx6w3jy75uFqChKLCnXsln1Kie+gupTYGFfuY5KnzXZxavFfhTwzSInTy7Xnk/xLLuG9qGlRmqDoVNS6VM47HVeM3luwWHTtZJf+sDLq2v1oDNzfHZC5TRucBYrPw42MIi+vNqXSLcRP3b2VA82i3b3BV0HYWAvKyNPtvNCkk3iiiIX63eFmDg8i71fRzhUHeWwXzFtlDMxE01WhQ0stKXYhi2iVRs/uDwHmxWmjkytBGerZ0eOxuIoBvas04tnKhwczeRMiA/ovnU7mxUOAF/VOPnjyiH84p3hBL6qYXLvw/xq6/eZUr4UUVpFJKI/k0qvIjcWMvJlVZTuORYKcix072bl6gku7FZJCAfQbrEdYag2NKRJASEiI0XkKhEx9gjDUVNcGuS+56saeBodKWUVEYIh+I9LPcy5s4Dfz8rn2+e5UUClDw5VR7nq7FS7wsbdtUSjpAiCWZPdVHhVWu+VSp/i/mtzufUiN047rN8V4q65lRz2KwQ9K873WIhEYce+EGUVurhOnluYMTY7sUqJR6bH1U1nDbJjtWjhYLPA8i1B/r6khglDs3hvq6596ozZScMx91Br7FfqsOuB/LODegC2W7WKJh3JgsjtSI2KtlnqtieTMFw3RhrPp4iykH1iDic6gxy253Hdp49z/8c/ZnRFMXdPd3PxaAcWi07fYYkJjMvHZVMbVsxdUoPHIYwdqJ9/dszjaWR/O08s9TL7kaP7nhiOnkYFhIj8HHgBuBJ4U0Rua7dWGboc8QHyaF0ai0uDiQHzqRW+xPkHKiOJ8ctmhQOVqbqYsooovfJTv+5xwdBUepGiwQ5uON9NzzwrIvraZ/a3IwJ/X1KTUDNZRNtArjrHlTg/OYXIm+v9TD7dwSdfRpg03EF2llAb0ZlRLQLvba1NGHoDIX0fgO9N9yTUSnGjdTgKZ55iIxTRSfKSx21BD/wThtZ548RrbkvsuPj1vPUefaOpQkBXtZtwGthtDbaPvWMCL+zuwe+G3MdfB/6AqFi4Y8/D5P35h5Rs9xKNQjeXJFYxGz6pxWrRzgAPLazhnU21WAS+fZ4LiwXWldbyzfHZ9Mq3NPk9aa0Jh6FxmlpBXAOMVEpdB5wFzG6fJhm6IsmZU5Nn8PUL0zRFXMhMHaltB95AlBc+8PHnRYd5qySA3aZn8FNHOnl7Q4BXi32J3E5Wi06PkcyuA2HyYsVtmsoDVTTYwf3X5jLnjgJunOLmi0MRJg53YEFRXqNY9XEIS2ymP3eJlwdfr2Z/eYRFJQFG9ItnktU1tq8oyub6SW4euiWf807L0qser8Jm1Z4+Nqs2SCdsDaJdY60W+Ott+Yk2xYXb1r21KbbkU3pa8AZJ1I4GEsbzAo/+uVsa+dU3pVrKvX4yedee1yDHU+71k9l14sB4qmdKTxzDL4b9jidOvp1NjuEc9NmwWeDCPpX4Q3DyCVY2fRriojOd/OXWPK4+R9t1Jg7PomiwI/E9mTLCyYHKaKPfk9aacBiapinNZFAp5QNQSh0SEWOvMLSYowlEinuX7C+P0D3HQu8CK5NPd7Bia5BASLH50zDZWWCz6hl8XLWUXK960nAH60prGdbXnuLGeNURFreJb19UEqDKD70LrOS5hB37w1R6VSJieFdZmJH97awrrU2k76gfqTt2kIPSskgDj52X1/jJcQrhiOLRt2uoDekEf8+/700YgF8p1gPm2xuDKZ5Me75smHNEBHxBsFr0vmjSIVk2fW4oooMHI2mEhK27h/HfGs7mT8O4xg0h/5whWK16sRFK+vj6dhe+OKRALJT2P5/VMY+lPtW7mfbPn9Gz+3jec1wDnMgH22t5aY2fwnwrHqfw0V69col/T5JXcem+J3FBUuWL8ssXDyfSrry82pQcbU2aEhADRGRB7HU8UC7+HqXUZZncIJb0bz2wTyl1adL2PwE3K6UahLqKyCnANmBHbNNapdR3M7mfoXPS0kCkZFfRPyyo5jvnu3j6PZ1BdWChjTfX+ymriJLntnDJ2LqBPV296uLSIPPeraEi5lGT59Y6lyNNG13/+Puer+KeGTk8t8qXEALxNA19C6wJr6cHX69mZH87t0/1NOpnH09Z/t7WYGJl4Ii51cZtFA6bVjcVDdEV3Yb01u6r9b2gTj7BwhflUbKs2gYRV80lHxOO1HmBeS47u2GpVLuNXlefw5bP6raFonDBmXqVluxF9sUhlWjDwZhwOGuwnYmnDGDZM5czaf8izlizjpI+F1B0/SzI6w7UuSRv3xeiV56FZVsCLP8oyIh+du57vor95ZEGCQDLKiJU1ERY8GEg4YG2Y3+Ih96oaZd6FMcLTQmIy+u9/38tvMfd6MG+W3yDiIwF8hs9Q/OJUmpkC+9p6GS0NBApWTVVmG/FapGUYLRcl4WHF1anRFZD48LHahV+cLmnVROyxWe9yaukQb1s7C+P6JxEwF9uzeOplV6Kd4a447GKRGLAdPe9fpKb6ye5KS4NJqKxQdcCctiECq/C4xR27NODdmmZvkfywO92grcWzh1WZ/wOpIm6djmEc4dn8c7GIO6iIbgdcODF1OC5vHNO5VB1FKddXyMahdXbg2nLXvTME76srAsC/LA0xO4DFi6/7iZWfz0D3niWSfvfhQc+ZN3sv/Pmpgj7y7X9KC68X1nrp2+BJSEkrRY4vZ+d14r9fFIWZsd+nbdq/nIf00Y561YR5RGsFnjJrCJajYziIESkB4BS6uARXVykLzAf+DXwfaXUpbEVxbvA9UBpEyuIhUqp0zO9l4mD6Py0JBBp9qPlPDI7H5tVEquJGya5eGhhDT+4PCcR37Blb6hBLIPNkhocl07Nc7QJ2eIpI8JJGpDesSp2i0oCnHdaVooaafEmP68W+3n09oK012rs+SSvpB58vZp7LvXw9Epfou/xZ5K8ishzC6ecaONfe0KJbfFBPsum1QK14brjzzvNwQ3nubnv+SpOzLWwaY+WKH+5NY8f/6MSX9zYHbtHvCRp73wL+yuiie3xVcXNF7iwWCTRp155ujTrqbavOLfbfv7pG4kFGPX5u3xUOIFau4trznXxSVk4ES9SmG9h1IAs1pXW0re7lS2fhbhnRg4VNRHmLfWRnQVWi/CN0U6WbQkwbpCDtzYGTCzFEdBUHERT9SAE+DnwPbQxW0QkDPxZKfVAhvd+CPgxqQWG7gIWKKXKRBq6zSXRX0Q2AoeBnyqlVqVp42xixvN+/fpl2CRDR9GSCmDJqqn4uU+9p2MR4sV64ikt4naEXJf+Xp091MHG3bXsL4/wxFIvkShHlZCt/gB+am8bW/aGOL2fnc2f1ql6urlIpAP/aG+YmePrVknnDXfw4mp/g2s/u9LLyo+DRKPQK9/CyP72RMrsZOPt0D52ehdYsVrrVlJXFGUnnolCu7FmZ+liPqu2hbBZtOfTN0Y5+ObZbl4t9vH2hkBiBZDvFhQwdpA24pdVRPjZ1d347SuH+exghLvmViY8pezWOruDr1bP7vfHXHyjioS0cTuEfI815XOLC+OJY/rz+JITUcCAmp3c8Pk8gl+9yOLCy3ht1VSyXNpB4NHb9cQAYFhfOw+/UU0k5qoMdl5Z6+ewTxFVig+21zJzvItcl4V1u2pZVBIwAqIVaErF9J/AucBZSqk9ACIyAHhURP5TKfXHpi4cKy70lVKqRETOj23rDVwNnN9Mu8qAfjHj+BjgNRE5TSl1OPkgpdQcdAEjxo4d2zVCwg0p1FdN5bosiDSMtk0WPvc9X5UwEMfPW7YlwIur/byx3p9S6CfThGzp6gU89EY100Y52bQnxMWj9X8FbP+iTuAM6mVNtKmsIkKOUxqEExSXBln5cTClQNH85d5E5tOiwY4U9VX8mdxwnismEPUzyXcLZw91sHp7kEqvSngyxQXBqX31yilun7lnhoeH3qjh/2bl61TqsecctwN4g4rJpzvYsT+cUANl2YRQpO6nFolqI3jc68rtgHGDHWzcU8uDr1fTPcfC6f10OvND1VHy3MInZWEUumzshWcW8WnJg9S+OJ8Znz7DxH1v8mavmRzsMTmhips+xsmYAVmEo3pFEeeqc1zMXaIN98k1Ji4f5+SJZalJFg0toynPpO8A18WFA4BSajdwA/BvGVx7AnCZiHwKPA9MAbYCg4Bdse0uEdlV/0SlVDCeWlwpVQJ8ArRPZRVDp6JosIMrirJ5bpUvpcRnU7PDsooIG3fXprjVThmhZ6XvbAw0W9o0HencdKNRHYBXVhFhxths7r82V5clTaK4NMRTy70EaqMoBYf92tMp2R1zUYnO8TRlhDPFBTh+bUgtBxt/Jk+t0BXu4s+kwqsSwiGOoAVEvlsShuAd+0MUeCzMW+olP8lQH3/OZRVRXlnrZ9zgLL41wcV1E12c0M3Ceac5yHYI3XMs3HOph+ljnNrwrXTAXr5buG6Smy17Q9xyoYebL3ATDCne21pLRY1WQVktkqghflIPGzarcMq404je82seO/M+Dmb15OIv38Bh1X2o8Ud5fpWXJ5dr76pkV+WiwQ48Tu22m/zdyPdYj/ssrK1FUysIu1Lq6/oblVIHRcSe7oR6x90L3AsQW0H8MNmLKba9Rik1qP65MZtHuVIqElu1DAZ2N3dPQ9fkSFVThflW9tfL+KpXCloH3pKKeencdHvFrte7QA/eVb4oL6/x0T3HwnfOc/HUez4OVUcJhiFYowXDiH529nwV5qXVeoYbd9+1WkhZ3QzqZeNARZTCAj3QNbeSKi4N6oSEMeHQLVvP9GvDenZ/8ok2dpWF+fviGg77FR4niAhXTqhbTSU/57g6Le4mHHfVjbuWPr1SpygZ0tvGUyt8iMCV57hSBCnAK2t9OOxaLRU3d8b/z1vq5baLhEG9bESiivXWUyk59eecSCVWRxaXDBSGv/5rluVPZp1/HMNPymrgqiwiuBw6b5bJwtr6NCUgalu4r0WIyGXAWKXUz4FJwAMiEgKiwHeVUuWtfU9D29GRmTGnj3HyxFIvy7YEmDLCyRvr/byzUc/SbVZa1JZ0brqjBmTx1cYAI/vbeXKZl1BYoVAUDXHw9EofM8dnM3eJlwKP8Lt/q3Pa27q3locW1iRUVs+u9NKvhzWljOmyLQEsFhKrm+bKwb682qfrWAh8c7yTk06wMWdxDVGlU3R8tDdEOKL7L0A3V9OfSTqh/Pi7Xgb1siXsAvG2KFUnqOLHxKn0KnKcUBvSNoVdB8L8cUG1Li6k4Mll3oSHFIBCCHkKGH1KFju3fsa4QDnf3fMnPss+hcDA71A5bkzKM7jm3COLY2mM4zWTa3M0JSDOFJHDabYL0PyaPAml1ApgRZrtnqTXC4AFsdcvAy8fyT0MnYeOru9bNNjBJ2VhXlnr58XVfh28NlBnWy0akpVi/M2UdG6660prmTTcwaY9uv4B6B/Hpj2hlBlseY3ivuer6gadmA1i3OAsnlvl40BFlK+qovTvaeWdjQEWlQSwWmDScEejdpb6xF1hlYLV22s5UBkg3y34gipRAzteLrWl1HcYiJcNTS7wk06Qjhxg55MDdbU6xgy08eGuMBOHaycC0J5VI/vb2bQnxOQRDhZvCoKjkMP/9VciG1eQ/faznPzSL6nuM5SiH/0SPKleZ0cjEHJdQiTacBXSkut2NRoVEEopo8QztIj6aobkAj7t9YO7fpKbgYU2nlzmJRKBzw9FmTlezyyH9bUfcVuam8FDal3juJDMztJBbhOGZvHqWl2A6P1tWgdfvLOWG6dogZO8yokbZo/0WRV4BBFJxITEVyoFHgsOe5MegwnSeWrFU6TnuYQ5i2uaVOfUF6QA738c4htjnKzZEeT1df6EMF2zPUilT8Uy2TqZWeTiR/MrePdfWl1mswIWKw+XF2E9awz/mb+G8uL1DHd7KC4Nsmr1fnb6chud8Te2Kqg/gfnps1VElaLK17DglBEQBkMr01nq+8ZVHo8kuUseTVuas4UkD47x5HxvbwyglGLZliCHqhVvbwzgdmhPpikjHAlhsmlPKOEJ1BLhkOeWWFK+KHMWVzN1pJO3NgRjsQ4qxdW2MeoXOHrhAx/LPwricQo3TXGT77EwZ3ENj79bQ1UsvqS+kKwvSC2xBIfvbQ3gDWjbSHaWNpwf9usi1uEI9MrT89GrznHxz/d9VMcKGD25TFeeu/rcbpwwYAY/LzuXW3bVsmTVl/z3+v+Ewaezt/91PFbcJ+X+Ta1i609gymui3H2phxc+8CfON/WoNUZAGFqdzlTftzXakql+OnVw1Mn54jryRSUBXTI0Ctec6+LxJd6EK2nxzlqmjHCwdLOf8hpapAKLD6xKCdV+xctrtD1DgBq/SiS7a+qayUboBR/6cdiFaSMdrNhay4IP/VxRlM3sqR6eW+Xj97MaDyxMzldVVh4hFNFCQAC7TdfkzrIJl53l5P1tQWrD8PIaHxZL3bnxuBWHXRIrv+37QhTmW1lUEuD6yflYCr8N77zIKVt/wP8MGcvTy6+kaPColL6kW8XWn8AU5ltBkSIQWvp97Wq2DCMgDK1OZ6rve7RtOVJ7SnyVkaxuim9P1tcvKtHG7cWbAoQj8MH2WsafqmMprpvoarEKbFFJAG9Qq4PCR6hXjw+cv3zxcCL31eTTnbxSHOD7l+lcUz+7ultiIM1UhfPGen9CQDnsQqBWx1O8tMaPUjqLbcnuWl5d62fMgCxyXRZcDr3iu26ii0G9bGzfF2LO4hqsFm34fmatYmif6ewZO5HhOxcxbfebzA7/N8x8Anr0anIVW3/SMH2MkyeWeSnwWAhHVIu/rx1te2sLjIAwtDqZ6OuPlba01J7SnGCaPsbJa8V+whGdymLPwUhif2uowOrbQxaVBDh0OJqo+5yu7fGBMz64FuZbee/jIIX51kSb4jPrVHVUhNfXBZi7xMtLq30IcHOspjfAzCIXxTtrE/2eu8TLbRe5QHRNjy17Q4zqb+e9rbXcOacirXdSXiw6/pYLPdrr6wQr720NcvHobsy47iaeWXwJgQ0lrH8xi94FVUyveIfPN42h/5i68Kn6dTqSXYajCqJRlbh/S76vncH21toYAWFoE1qSVqOtOJq2HIk9pf6Muq4WREPBFP//5DIv35tbSWFS8r64KuVoiLc7eSDv38PKXXMrG1VhxQdOt0P4n2cqKa9RvLjaz8j+dnbsC+F2SMJFde4SLx6nsOGTWrbsDXHjFDeRiErEflTURIA6td7l45zMW+rj1bV+CjzCvvIIyz8KMnN8NrkuC08u89K7wNogJ1Y6gXfJ2GyeWOpl4vCshBfU6k/tnFU0kcKDEb59lqLvipdwPTqfw8Mn4LrqBnZZ+yYEcLpJwzXnujL+jjS2cuostrfWxAgIg6EJMrVhxAficYOzUEpRVh7hy0pdPa6xRIDxASle+CauSjladVw8aO6OxyoAndm1yhdlz0HtIdWYCivuHrxia5CagPaKOukEK5s/DbFpT0j7t2fB7dN0ig6LBVZsDfKN0U6G9rETjijKa3TN6dfXBTj71Dpv+HyPlXy3JDyYXlnrZ9oonUJjx37tJnzrRalR6MkkD75Fgx08vsTLjn0hDlQqDm4KcPFoJzPGZnPnnAqGDCxg5w/nsveZFzl3x1tYHlhNqOd4vvXNWxg1OC9xjZZMGppSI3Um21trYQSE4bjiSI2ImdowFpUEGDc4q0H+p1fW+hlYaGsyIO2TsjAPL6xOBLJNHNbyFU9xaZB/vq+jl+1WOOzX2VWfec+LzSpcc66ryVntjv1hvn9ZDlW+KItKAmz+LEyBx0JFTZR8j4Ubp9QlDJwwNIsXV/vZuLuWmUV1pVYvHuVMpPVIfmbxSOvrJroS139rQ4ACj07T0VSf6w++hbH7f7C9LtVJ8mA8oH8e/+/Eq7jwB9fC4lcYsfItGBj7zCIRsLZs0G5KjdSZbG+thREQhuOGlhgRM7Vh6IhilTJ4TBnh5MXV/iYzixaXBtmyN8Tdl+aktKmlRW8WlQTIsgmzp3qo8mmbQ7VfZ1112GlWhRWfqdusdQN2OKK447EKDlVHUxIGvrrWFzsnmrLyyXVZyHNLo88sPogmJ9i78hxX2vbEqT/4juxv5+U1fi4e7UQpxZPLa1i/S7sKx5M15rqE+xYqyiou46SzLmHqQQ9FuQr+cC8U9IBLr4eefY7o+TalRupMtrfWwggIw3FDS42ImagjCvOtlDWa/6lxHXRrGzbjqS+SU2K8ulYHp1V6VbMqrMbUJDYr5LosKZHU+8sjibThTy7zcvk4LRzipVybdwXObBBNNrLHV1q9C6ycd5qOYi+riFJWEaVosJ3vnOfmvY+DvLTGj8MGt1zoSnhSPbHUyxOLq/hO4GTGf/gO1rXLYfxkuORa6HVSRs+3seeT55JExt7CfCu3XNg16lEYAWE4bmhLI2L9/E/JKbvjhXfqU1waZH95hAdfr06Jnj6aNhXmWwmGVIOBPB6l3Vw23MbUJBOHOVj/SS1/X1KjM9Yq+GBbELdTGDswix37wzyxzJvRgH8k+v/GVn3JqsH4imHTnhDfe7ySwnwrbgc47JaEF9e60lq+OT6b97cF6X7FbH779iXcHH6bwpK3YO0yuOsXcGZRs+1J93zmLNaZZuN2pK7g3hrHCAjDcUNbGhHr538qzLckbBLpZuvxgS+e+dVqlcSgkuuytLhN08c4+ef7vgYDucshGXnqNDXDH1ho4+XVPh56Qw+IeW7hmgmZe/+0hExWWPECR8l1Pm57pBxvMJpyjUG9bLy0xs/QPnauvLgPf1t1Hff/9lpYtgCGnqlP3P4vcOfASQPStifd87FatPttV3JvjWMEhOG4oa2NiPH8T3EjeDxpX7pBIl3q7BvOcyVSZ7e0TcnZXVs6kDc2w+8I1+VMVn3pBH/3HEuDayRPBhLX6FYAV+jyNsWlQfrOeYw+h3ezrcdYwtOvY8TE0xq0qf5zmP1oeZdzb41jBIThuKE9jIiZDqLJxuDkNiWnzm7rNnQmGvMuy2TVl07w14YVobDiJ09VohT81z8qUZAIwKt/jfiKbsIV97Nj1QKK9r6Je/4P2PfuGPrcchP0a1C2JkFXdG+NYwSE4biiswyemaTObgnHYi6gprzLMln1pRP8Ywdmsf4THUQnQE1A4bDrXFjpDPVxN+UPSmHWjf/G1kMz+fq115m0702efmEngy8/iaKBdl2+rh5d0b01jhEQBkMH0JqDSkIolEewWGDaKGeKCyl0bmNpU3aGeJBhc6u++oL/vuermD21zi5QXBrk1bV+5i3VEdv1r5Hsplzli/LaZgs33PRt7n1tKp7cbLYW++m9+gVOOrgVZlwPp56Rcu9M2ngsYgSEwdABtNagkjz7fnall3OHOVj+UZDeBVaKBjtaZCxt71VIc3aGlqz6klONxPvSK8+CQNrI9uQytfFkhQD5J7g5UBnl+5e52fi8m5PKPoPf/xiGjNCCYuhIEOk0K9PWxggIg6GDaI1BJXn2faAyypQRTvr1sCWEwpEaS9s7I2lxaRCrRacFSXb1PVodfmG+lTfW+9NGtqcLQkx2U95fHiESUTy90pdwUx7Uy8YfXOdz2W8ugVVvw1svwIP3wkXfhGtmH+1j6LQ0VKgZDIZjhuTZd9yukSwUjnSgTRY4ydXV4um6W5O4MJo60kn3HEui6t6rxb5ErENLmT7GyTsbA0w+3ZEQdMs/CjJtlDNtX4oGO5g03MEra3Uiw6fe8yXclKePcdY9xywHXHA5/OYJ+PZdcNYkfYFDX0HJ+9rI0YUwKwiDoZ1pTRVOsrE7bteYfLqDXnmWFiX+a+1gwqb6miyMehfoQkDl1YrFmwLcOOXoPbnmLvHy/jYdVR1X4Y0ZkMVbG9ILu7ib8surdUbalVuD2Kw6c228RngCexZMvlT37/kqztr6KpeWvYq/oC/Zl18DRZPBduwPr8d+DwyGY4jWVuEkG7vHDMhif3mEV9b6M4qaTkdrumw219f6GVqLBjsIR3RNhtZQZ8Uz1yb3pblU6vF2PLvSy8qPg0SjUJhvYdQAvZpIVk+l9O/KW9i3bDCWRc+T/cSDsOApuOQ6mPSNo+5HR2IEhMHQjhxJ7qVMVhrpjN03XdDy2Xdrelc119e2jh84mr7s2B/mnhk5KW0b1tee8jnV71+faVPYftq5LH3tA26oXgCfbKsTEMEAOFquMusojIAwGNqRTFU4R7LSaE0PmtZ02Wyur20dP3A0qdTTtb2iJsL+8gizHy1Pm5wRYFChnT9Yz+SGn5wP4VgOrj074I//A1NmwAVXQE7j9bw7G0ZAGAztSKaz5o4sX9laAqe5vrZ1/MDRpFKv3/biUm3L6J5j4VfX57LrQJiH3qjmjfX+lBxQif6JaDsFgCNbx00sfA4WvwLnTYepV0L+Ca3Sz7bECAiDoR3JdNbcFcpXZhoF3VYC72iEbP22v7rWjwAzx2cnvLumjdKeUsP62pteAfXuB//+c9j/mXaPXfo6rF0Ov38KbPa09+8sGAFhMLQSLbUZpJs1d4X8Ph0dYXw0QrZ+25WCm+vZdmaMzWZRSSDz/vU+GW75EVz2Hdi3RwuHaBRengfnXAh9TmlxX9sKIyAMhlagtW0GXSW/T0dGGB+tkE1u+33PV5HvSQ0b23UgTO8Ca6M1xxulRy/9B3DgC1jxJrzzEowcD9/4FgwcfmTXa0NMoJzB0Aq0doBZ0WAHVxRl89wqH3fOqWiRy+rxTlzIbt8XIhypq6bXkgC81rxWCr37we/+AZd/B0q3wm++D7/7IVR8fXTXbSXMCsJgaAXawmbQEbPvYzEbbGO0porrSK91RM/RkwMzvq3Tdrz/jo7IzsnT+/Z/Bif2aTTorq0/LyMgDIYmyPQH2BVsBu2dh6k9aG0X4Eyu1eLn6MyGC6/QfwChWp3vyWqDqd+EiRenxFK0x+dlVEwGQyPEf4DXTXTxyOx8rpvo4rVineytPm2mgmhH2jMPU1em1Z6j1Qb/djd07wHP/w1+/B14/Smormrd+zSBWUEYDI1wJG6SHe2x0xp0BdfazkCrPUeLBc4s0n+7Poa3X4A3noFBp8FpoykrD7f552UEhMHQCEf6Qz/WawJ0BTVZZ6BNnuOg4XDXL+DLfXBibwD+7dA/8f6litzb/gNc7ta5Tz2MislgaIT4Dz2ZrjxgZqImKy4N6mptj5Zz3/NVadVtxzttqm7s2UdHaQOn9Mnm8KdfsP2Qrc3UmmYFYTA0QleJRciU5tRkXdGI3Ra0l7qx7803UbzzWhZ9EKCswtsm9xGlVKtdrCMZO3asWr9+fUc3w9DF6Epun0fLfc9XpU2fnVw72nDsISIlSqmx6faZFYTB0ATHul2hNTFG7OMPY4MwGAwZcbzZZAztICBExCoiG0VkYb3tfxKRmibOu1dEdonIDhGZ1tbtNBgMTdPRsR7GQN7+tIeK6W5gG9AtvkFExgL5jZ0gIsOBa4HTgN7AuyIyRCll1rIGQwfRkbEexkDeMbSpgBCRvsAlwK+B78e2WYHfA9cDMxs59XLgeaVUENgjIruAccCatmyvwWBomo6yyXRkAaXjmbZeQTwE/BjISdp2F7BAKVUmMX/eNPQB1ia9/yK2LQURmQ3MBujXr18rNNdgaB2M91PrYgzkHUObCQgRuRT4SilVIiLnx7b1Bq4Gzm+Neyil5gBzQLu5tsY1DYajxahDWp+uHOXd2GSiM0wy2nIFMQG4TESmA060DWIrEAR2xVYPLhHZpZQaVO/cfcBJSe/7xrYZDJ0eow5pfbpq0GJjk4lPysJs2Rvq8ElGmwkIpdS9wL0AsRXED5VSlyYfIyI1aYQDwALgWRH5A9pIPRhY11ZtNRhaE6MOaX26QjLEdDQ2mXh4YTV3X5rT4ZOMThMoJyKXAWOVUj9XSm0VkReAj4Ew8O/Gg8lwrNCV1SEdSVcMWmxsMhGO0CkmGe0SKKeUWlF/9RDb7kl6vUAp9fOk979WSg1USp2qlHqrPdppMLQGHR0vYDh2aCz40GalUwQldpoVhMHQVeiq6hBD69OYbWXiMEensLkYAWEwtAFdUR1iaH2amkwUlwY7fJJhBITBYDB0II1NJjrDJMMICEOXpzP4kxsMxyJGQBi6NCZozWBoOSbdt6FLk+xnbrNKwp98UUmgo5tmMHR6jIAwdGlM0JrB0HKMgDB0aUyRG4Oh5RgBYejSmKA1g6HlGCO1oUtjgtYMhpZjBIShy9MZ/MkNhmMRo2IyGAwGQ1qMgDAYDAZDWoyAMBgMBkNajIAwGAwGQ1qMgDAYDAZDWkQp1dFtaBVE5CDwWStc6gTg61a4Tmelq/cPun4fTf+OfTpTH09WSvVIt6PLCIjWQkTWK6XGdnQ72oqu3j/o+n00/Tv2OVb6aFRMBoPBYEiLERAGg8FgSIsREA2Z09ENaGO6ev+g6/fR9O/Y55joo7FBGAwGgyEtZgVhMBgMhrQYAWEwGAyGtBx3AkJE/lNEtorIRyLynIg4ReRJEdkjIptifyNjx4qI/ElEdonIZhEZ3cHNbxYRuTvWt60ick9sW4GILBGR0tj//Nj2Y65/0GgffyEi+5I+w+lJx98b6+MOEZnWYQ1vAhGZJyJfichHSduO+HMTkVmx40tFZFZH9CUdR9i/80WkKumz/HnSORfHPsddIvKTjuhLOhrp39Wx72hURMbWOz7td7LT9U8pddz8AX2APUB27P0LwI3Ak8BVaY6fDrwFCDAeKO7oPjTTv9OBjwAXOpX7u8Ag4P+An8SO+Qnwu2Oxf8308RfAD9McPxz4F+AA+gOfANaO7keadk4CRgMfJW07os8NKAB2x/7nx17nd3TfWtC/84GFaa5hjX1+A4Cs2Oc6vKP71kT/hgGnAiuAsc19Jztj/467FQR6UMkWERt6kNnfxLGXA/9QmrVAnogUtkcjW8gw9GDhU0qFgfeAb6L7MT92zHzgitjrY61/0HgfG+Ny4HmlVFAptQfYBYxrh3YeEUqplUB5vc1H+rlNA5YopcqVUhXAEuDiNm98Bhxh/xpjHLBLKbVbKVULPB+7RoeTrn9KqW1KqR1pDm/sO9np+ndcCQil1D7g/wF7gTKgSim1OLb717Hl+h9FJF5dpg/wedIlvoht66x8BEwUke4i4kLPNE8CeiqlymLHHAB6xl4fa/2DxvsIcFfsM5wXV1dwbPYxzpF+bsdaXxvrH8DZIvIvEXlLRE6LbTvW+tcYx8znd1wJiNigcTl6WdcbcIvIDcC9wFDgLPTy/L86rJFHgVJqG/A7YDHwNrAJiNQ7RgHHrG9zE318FBgIjEQL/wc7poVtw7H+uTVHvf5tQOcHOhP4M/BaR7XreOe4EhDAhcAepdRBpVQIeAU4RylVFluuB4EnqFNB7KNudgrQN7at06KUelwpNUYpNQmoAHYCX8ZVR7H/X8UOP+b6B+n7qJT6UikVUUpFgb9zDH+GSRzp53as9TVt/5RSh5VSNbHXiwC7iJzAsde/xjhmPr/jTUDsBcaLiEtEBLgA2Jb0JRW0HjTuibAA+LeY18h4tEqqLM11Ow0icmLsfz+0bv5ZdD/iHi2zgNdjr4+5/kH6Ptazncwk9TO8VkQcItIfGAysa8/2HgVH+rm9A0wVkfzYanlqbFtnJW3/RKRX7LeIiIxDj1OHgA+BwSLSX0SygGtj1zjWaOw72fn615EW8o74A+4HtqMHkKfQngTLgC2xbU8DntixAvwV7VmwhSRPhM76B6wCPkZ7QFwQ29YdWAqUor1+Co7V/jXRx6difdiM/lEVJh3/P7E+7gC+0dHtb6RPz6FVYyG07vmWlnxuwM1oo+cu4KaO7lcL+3cXsDX2+a5Fr/Lj15mOXhV/AvxPR/ermf7NjL0OAl8C7zT3nexs/TOpNgwGg8GQluNNxWQwGAyGDDECwmAwGAxpMQLCYDAYDGkxAsJgMBgMaTECwmAwGAxpMQLC0K6ISCSWoXNrLJXCD0TEEtuXnMVzm4jc1wr3uyeWkiP+fpGI5LXCdUdKUsbYVrje72PP5Petdc2OQkTyROTOjm6H4egxbq6GdkVEapRSntjrE9GBfB8ope4TkfPRGVkvFRE3Oo3GNUqpDUdxv0/RcQJfH23b6133xth172ql61Wh4wAiTRxjUzpBYadGRE5BZ2M9vaPbYjg6zArC0GEopb4CZqOT7Em9fV6gBJ3KOwUR+ZGIfBhLzHd/bJtbRN6MrUo+EpFrROQ/0Dm3lovI8thxn4rICSJyiohsF10LZKeIPCMiF4rIB6LrE4yLHT9ORNaIyEYRWS0ip8aiXB8Aromtdq6J3X+eiKyLHdsgC2cs8vn3sfZtEZFrYtsXAB6gJL4t6ZxfiMhTIvIB8FSs3ctifV8qIv1ExCq6nonEZu8REZkUO3+liAyud02riPy/WDs2i8j3YtsviLV9S6wvjuRnFns9VkRWJLVtnoisEJHdsecN8FtgYOzZHPMrouOajo7UM3/H1x9Qk2ZbJTqT5/nE6gCgo2w/BU6rd+xUdMF3QU9wFqJz8V8J/D3puNzY/0+BE5K2fwqcAJwChIERseuUAPNi170ceC12fDfAFnt9IfBy7PWNwF+Srvu/wA2x13noaFh3vbZfiU7BbY31dy+xiO90zyW2/RextsVrmLwBzIq9vjmpnW8DpwGXolM2/A86S8CeNNe8A3gpqV8FgBOdSXRIbNs/gHvqP0NgLLAiqW2rY/c5AZ0Owx57th+l64/5O7b+zArC0NmYKCIb0dlaf6uU2lpv/9TY30Z01s+h6Fw2W4CLROR3IjJRKVWVwb32KKW2KJ3gbyuwVOmRbwt6kAPIBV4UXSnsj+hBOB1TgZ+IyCZ0gRgn0K/eMecCzymdVPBLdC2LszJo5wKllD/2+my0Wg50epFzY69XoQXlJOA3se1noYVFfS4EHlMxdZVSqhxd2GaPUmpn7Jj5sWs1x5tK1zX4Gp1sr2dzJxiOHWwd3QDD8Y2IDECn6/4KXQxolVLq0qZOAX6jlHoszbVGo3PZ/EpEliqlHmjm9sGk19Gk91Hqfhu/BJYrpWbGdOsrmmjXlSp9gZijxZvBMSvRK4PewM+BH6FXZKta4f5h6tTRznr7kp9hBDOmdCnMCsLQYYhID+BvaFVNpt4S7wA3i0jc0N1HRE4Ukd6ATyn1NPB7dPlHgGog5yiamUtdyuUbk7bXv+47wPfithQRGZXmWqvQdgtrrO+TOPLMsqvRWT4Bvk2dAFgHnANElVIBtIH/drTgqM8S4HbRVRURkQJ00rhTRCRu8/kOeoUDWsU0Jvb6ygzaeLTP3NBJMALC0N5kx4yXW9EZPBejM+xmhNIVAJ8F1ojIFrQuPQdtS1gXU/HcB/wqdsoc4O24kboF/B/wm5jaK3l2vBwYHjdSo1cadmBzrG+/THOtV9HZZv+FziD8Y6XUgSNsz/eAm0RkM3oQvxtA6Vomn6Ozn4IWHDlodVl95qLtH5tF5F/A9TGhchNanbYFvYr6W+z4+4GHRWQ99QpQpUMpdQj4IGYEN0bqYxjj5mowGAyGtJgVhMFgMBjSYgSEwWAwGNJiBITBYDAY0mIEhMFgMBjSYgSEwWAwGNJiBITBYDAY0mIEhMFgMBjS8v8B+A1fHhVpR6AAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import statistics\n",
    "\n",
    "true_mean_age = statistics.mean(age)\n",
    "\n",
    "# The light blue circles are DP means\n",
    "plt.plot(history_count, history_mean, 'o', fillstyle='none', color = 'cornflowerblue')\n",
    "\n",
    "def compute_expected_mean(count):\n",
    "    count = max(count, size)\n",
    "    return ((true_mean_age * size) + (count - size) * age_prior) / count\n",
    "\n",
    "expected_count = list(range(min(history_count), max(history_count)))\n",
    "expected_mean = list(map(compute_expected_mean, expected_count))\n",
    "\n",
    "# The dark blue dots are the average DP mean per dataset size\n",
    "for count in expected_count:\n",
    "    sims = [m for c, m in zip(history_count, history_mean) if c == count]\n",
    "    if len(sims) > 6:\n",
    "        plt.plot(count, statistics.mean(sims), 'o', color = 'indigo')\n",
    "\n",
    "# The red line is the expected value by dp release of dataset size\n",
    "plt.plot(expected_count, expected_mean, linestyle='--', color = 'tomato')\n",
    "plt.ylabel('DP Release of Age')\n",
    "plt.xlabel('DP estimate of row count')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this plot, the red dashed line is the expected outcome,\n",
    "and each of the points represents a `(count, mean)` tuple from one iteration of the simulation.\n",
    "Due to the behavior of the resize preprocess transformation,\n",
    "underestimated counts lead to higher variance means,\n",
    "and overestimated counts bias the mean closer to the imputation constant.\n",
    "On the other hand, underestimated counts are unbiased, and overestimated counts have reduced variance.\n",
    "Keep in mind that it is valid to postprocess the count to be smaller,\n",
    "reduce the likelihood of introducing bias by imputing.\n",
    "If the count is overestimated, the amount of bias introduced to the statistic\n",
    "by imputation when resizing depends on how much the count estimate differs from the true dataset count,\n",
    "and how much the imputation constant differs from the true dataset mean.\n",
    "Since both of these quantities are private (and unknowable), they are not accounted for in accuracy estimates."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "In the next plot, we see the range of DP means calculated as a function of the resized row count.\n",
    "Note that the range of possible DP mean values decreases as the resized count increases, and that the DP mean gets\n",
    "closer to the prior for the true value: 38."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAgyklEQVR4nO3de3RddZ338fc3adqktNBLWoKNpUKjwjjII9ElIoiDLbSwKCgqM+KT8fJUGbWMLAdFfYA6OoqOjxiWCoiXOF7AGcWpQKWxGsClo7T2Qi/UHKBoim2aQKU1veTyff7YO+1Jek5yTnL2Pudkf15rdWXf9zcnp9/zO7+99/dn7o6IiCRHRbEDEBGReCnxi4gkjBK/iEjCKPGLiCSMEr+ISMJMKnYAuaitrfUFCxYUOwwRkbKyfv36LnefM3x5WST+BQsWsG7dumKHISJSVszsmUzL1dUjIpIwkbb4zWwnsB/oB/rcvdHMZgH3AguAncDb3P35KOMQEZFj4mjxv9Hdz3b3xnD+Y8Bad28A1obzIiISk2J09SwDWsLpFuCKIsQgIpJYUSd+B9aY2XozWx4uO9nd/xxO7wZOzrSjmS03s3Vmtm7v3r0RhykikhxR39XzenffZWZzgVYzeyJ9pbu7mWWsEufudwF3ATQ2NqqSnIhIgUTa4nf3XeHPTuA+4DXAHjM7BSD82RllDCIiMlRkLX4zOwGocPf94fRi4FPAKqAJ+Fz4878Lcb7m5mZSqVTGdR0dHQDU19cft27hwoWsWLGiECGIiJSFKLt6TgbuM7PB83zf3X9mZo8BPzSz9wDPAG+LMAYADh48GPUpRETKhpXDQCyNjY0+nid3B1v0zc3NhQpJRKTkmdn6tFvpj9KTuyIiCaPELyKSMEr8IiIJo8QvIpIwSvwiIgmjxC8ikjBK/CIiCaPELyKSMEr8IiIJo8QvIpIwSvwiIgmjxC8ikjBK/CIiCaPELyKSMEr8IiIJo8QvIpIwSvwiIgmjxC8ikjBK/CIiCaPELyKSMEr8IiIJo8QvIpIwSvwiIgmjxC8ikjBK/CIiCaPELyKSMEr8IiIJo8QvIpIwSvwiIgmjxC8ikjBK/CIiCaPELyKSMJEnfjOrNLMNZnZ/OH+Rmf3ezDaa2a/MbGHUMYiIyDFxtPivA7anzX8NeIe7nw18H/hkDDGIiEgo0sRvZvXApcDdaYsdODGcPgl4NsoYRERkqEkRH/824AZgetqy9wIPmtlB4AXgtZl2NLPlwHKA+fPnRxuliEiCRNbiN7PLgE53Xz9s1YeBpe5eD3wL+H+Z9nf3u9y90d0b58yZE1WYIiKJE2WL/zzgcjNbClQDJ5rZA8DL3f234Tb3Aj+LMAYRERkmsha/u9/o7vXuvgC4GvgFsAw4ycxeGm62iKEXfkVEJGJR9/EP4e59ZvZ/gB+Z2QDwPPDuOGMQEUm6WBK/u7cBbeH0fcB9cZxXRESOpyd3RUQSRolfRCRhlPhFRBJGiV9EJGGU+EVEEkaJX0QkYZT4RUQSRolfRCRhlPhFRBIm1pINE11zczOpVCrjuo6ODgDq6+szrl+4cCErVqyILDYRkUFK/DE5ePBgsUMQEQGU+AtqpBb74Lrm5ua4whERyUh9/CIiCaPELyKSMEr8IiIJo8QvIpIwSvwiIgmjxC8ikjBK/CIiCaPELyKSMEr8IiIJo8QvIpIwSvwiIgmjxC8ikjBlVaRtpLLHI2lvbwdGLqKWiUoli8hEVFaJP5VKseHxbQxMnZXXfnbEAVj/5O6c96noeS6vc4iIlIuySvwAA1NncejMyyI/T/W2+yM/R1J1dXWxcuVKbrnlFmbPnl3scEQSR338EruWlhY2b95MS0tLsUMRSSQlfolVV1cXq1evxt1ZvXo13d3dxQ5JJHGU+CVWLS0tuAfXXAYGBtTqFykCJX6JVWtrK729vQD09vayZs2aIkckkjyRJ34zqzSzDWZ2fzhvZvYZM/uDmW03M90vmSCLFi2iqqoKgKqqKhYvXlzkiESSJ44W/3XA9rT5fwReDLzc3c8A7okhBikRTU1NmBkAFRUVNDU1FTkikeSJNPGbWT1wKXB32uJrgU+5+wCAu3dGGYOUltraWpYsWYKZsWTJEt3OKVIEUbf4bwNuAAbSlp0OvN3M1pnZajNryLSjmS0Pt1m3d+/eiMOUODU1NXHWWWeptS9SJDk9wGVm84BT07d390dG2ecyoNPd15vZhWmrpgCH3L3RzN4MfBM4f/j+7n4XcBdAY2Oj5xKnBEYqbdHR0QFAfX39ceviKlFRW1vL7bffHvl5RCSzURO/md0KvB3YBvSHix0YMfED5wGXm9lSoBo40cy+C3QAPw63uQ/41hjiljE6ePBgsUMQkSLLpcV/BfAydz+cz4Hd/UbgRoCwxf8Rd7/GzD4HvBF4GngD8Id8jiujG6nVPriuubk5rnBKkspGSJLl0sf/FFBVwHN+DniLmT0OfBZ4bwGPLZITlY2QJMulxd8DbDSztcDRVr+759wZ7O5tQFs4vY/gTh+ZwMZ6nQGiv9YwvGxEU1OTWv2SKLkk/lXhPwmNZVwAjQlwTLGvM2QqG3H99dcXNSaROI2a+N1d34WHSaVS/GHL75k/rX/0jUOTe4NetUM7H8t5nz8eqMw7tlJRytcZMpWNUOKXJMnlrp4Ggr74MwnuzgHA3U+LMK6SN39aP59sPBDpOT69blqkx0+qRYsW8eCDD9Lb26uyEZJIuXT1fAu4GfgSwd0470LF3aRMZOqW6+3tPdri7+vro729/bhvKBOxi01kUC4JvMbd1wLm7s+4+y3o4qyUsaqqKiZNCto8s2bNOlo0TiQpcmnxHzazCqDdzD4I7ALUByFlIVur/dprr2Xnzp3cfffduqNHEieXFv91wFRgBXAOcA2gIitS1qqqqmhoaFDSl0TK5a6exwDMbMDd3xV9SCIiEqVRW/xmdq6ZbQOeCOdfaWZfjTwyERGJRC5dPbcBFwPdAO6+CbggwphERCRCOd2W6e5/GrYo9yeXRESkpORyV8+fzOx1gJtZFccPpShFMJayEaDSESKSW+J/P/BlYB7BrZxrgA9EGZSMLpVKsWHrBpiR547hWGgbdm3IfZ992VepbpFI+cnlrp4u4B0xxCL5mgEDFw6Mutl4VbRl7xFMpVI8sXEjdfkcL/y5b+PGnPfZPcK6OD98QB9AUv5yqdXzEuBDwAKGDr14eXRhSTmpA96DRXqOb5B99M1UKsXWx7czY+rcnI83cCSId9eT3XnFsa+nM6/tR6MBYaQYcunq+QnwDeCnDB00XaRkzJg6lze+/OrIz/PLJ+4p6PHSB4RRhVCJSy6J/5C7J3ucPilpHR0d/KVnf8GTcib7ejrxjsKMJ6ABYaRYckn8Xzazmwku6qaPwPX7yKIqcR0dHfx1f2XkZZOf2V/JCeFoVTLxaEAYKRYbfONl3cDss8A7gSc51tXj7v53Ecd21KxTz/BFH/8mqVSK/QePMHBC9K2iir92M71mMgsXLjxu3R/X/ge2fw/VEY+TcqgffPrJzL/onRlj8P2dMGnkv19B9Bk2fW7WOOyFTiZHHMIRwE/MHEMqleLwwd68+vjHal9PJ1NqqjK+L/K1+fHNDPQf6z2tqKzgrL89a9zHFRn0w/e/br27Nw5fnkuL/63Aae5+pPBhlacjp57LwKH9zJ0e7XNsz++vpKJ6eqTnkOKZNXMW3d3duDtmxqyZs4odkiRELol/C8Hd4oW9nSEPp805gXvfdy4rVvyA9U/v5tCZl0V+zuqnH+Kc0+toft/xLcwVK37AoWcei34Eri3TqF7w6swxbP0BG/r3xHY75/+qexHN7zs3Yxz7dnfGclfPjJdmiWHFD9j1dHc8F3efXsu802dn/Jvkq6uri6uuuoqBgQEqKir40Y9+pD5+Kagfvj/z8lwS/wzgCTN7jKF9/LqdU0rGvp7OvC7uHjj0PADTqmfmfZ55KDlLecsl8d8ceRQi4zCW/vb29ucAmHd6fkl8HrML0r8PwcXdioqKoy1+XdyVuOTy5O7DcQQi5Ws3Iz9gNdzgI1P5pNzdZK9OMZanaAf3aW4u3p3Kra2t9PX1AcHYv2vWrFHil1jk0uIvGR0dHVT0/IXqbfdHfq6Knm46OvoiP0+5G0vrd29YLmFGQ0PO+8wY47lKQbaSEjU1NfT09AyZ16DvEoeySvxSesq1tV0K6urq6O4Ovv+YGXV1+VQ8Ehm7ERO/mZ0NLAS2unvRSzHX19ez5/CkeO7q2XY/9fX6jyjjN9KH45VXXkl3dzfLli1TN4/EJmviN7ObCAZWXw983sw+6+5fjy2yEvfHA/k9ubunJ6hJefLU3G+//OOBSl6ad2RSTurq6jh06BBNTU3FDkUSZKQW/9uBs929x8xmAz8DlPgZW1/zkbBfu3pB7v3aLx3tXPtGLpmc0eCjB/lUm9hHMBqDFFxVVRUNDQ26f19iNVLiP+zuPQDu3m1meWaYiasU+rXHeqFzsA59w7zcP4CYV74XVuOkQWmkXIyU+E8zs1XhtAGnp83rAa4iG+t/el1YjU5bWxtdXV1MstyfYu4Pa2Vt2bQp53363Ono6FDilzEbKfEvGzb/71EGIjIRTDJj+uRob5bbf0S3Gcv4ZH2Hpj+4ZWZzwmV78z2BmVUC64Bd7n5Z2vJm4N3uHm1tY0m0bN0vo3WxjKUrpb6+nn3d+Y3o1dMXFPqbOim/Uq/19fV5bS+SbqS7egy4iWDYxYpwUR9wu7t/Ko9zXAdsB05MO3YjkF+RFCkrI/V3R5F081VTU1PwY46tdETwWpyax8NsYz2XyKCRvpN+GHg98Gp3fxrAzE4DvmZmH3b3L412cDOrBy4FPgNcHy6rBL4A/ANw5fjCl3IURdLNJs5+8FK46C+Si5ES/zuBRe7eNbjA3Z8ys2sIRuMaNfEDtwE3AOlF5T8IrHL3P9sIF8HMbDmwHGD+/Pk5nEpKiS48HlPq334keUa6RbMqPekPCvv5q0Y7sJldBnS6+/q0ZS8iGNjl9tH2d/e73L3R3RvnzJkz2uYiZammpibWb0AiMHKLf6QRt3IZjes84HIzWwpUE/TxbyWo6Z8KW/tTzSzl7uqwlAlLLXYpNSMl/lea2QsZlhtBIh+Ru98I3AhgZhcCH0m/qydcfkBJX0QkXiPdzhnxUOIiIlIMsZRldvc2oC3Dct3DLyISM9XfEUm4rq4uPvShDx0dG0AmPiV+kYRraWlh8+bNtLS0FDsUiYkSv0iCdXV1sXr1atyd1atXq9WfEEr8IgnW0tKChxVCBwYG1OpPCCV+kQRrbW2lt7cXgN7eXtasWVPkiCQOGmxdJCEylY6oqamhp6dnyPzwB85UNmLiUYtfJMHq6uqOTpvZkHmZuNTiF0mIbK32K6+8ku7ubpYtW8b1118fc1RSDEr8IhPIWMb97enpoaKigvb2do39mxBK/BPQWMsA6z9x+Wtra+O5rm6mTJqc8z5H+vswjKd25PeBcbjviMb+LVNK/AmjEsAT35RJk3nx9Oj76v+0f3fk55BoKPFPQGqBJVd9fT079m3Pa5/OnucAmDt11pjOJ+VHiV9kAhnLWLy97cF4S1NPzW8Y7JcxU2P/liklfpEJROP+Si6U+EUSIttFf437mzxll/grep6jetv9ee1jh4KBxLz6xLzOA3qYRSa+KVOm8MILL9Db20tV1ajDacsEUFaJf6z9ie3t+wFoOD2fRF6n/kuZULK12r/4xS+yatUqGhoa9ABXQpRV4h/r1031YYpkNrwsc1NTE7Nnzy52WBKxskr8IlJYmcoyR9nqH+nhwo6ODiD7LaK61lA4SvwFNNYnZkFvaimOTGWZi9Xdc/DgwaKcN4mU+GOiJ2alFC1atIgHH3zw6IXdxYsXR3q+kRo36pKNjxJ/AanFLuWmqamJ1atXA1BRUUFTU1ORIyqurq4uVq5cyS233DKhr3Uo8YskWG1tLUuWLGHVqlUsWbKkYMluLFVCR+sOzaaQ3aTpA89P5DuclPhFEq6pqYmdO3cWtLXf1tZGV1cXlZWVOe8zMDAAwOOPP57zPv39/QWrEJqkO5w0ApdIwtXW1nL77bdP2CSXq5aWlqMfPv39/RN64Hm1+EWk4C688MIxd/U0NDTktV+hHrRsbW2lr68PgL6+vqLe4RQ1JX4RKbhyLBZ3/vnn89BDDx2dv+CCC4oSRxyU+EUkNiNd9N2xYweHDx/m2muvzVgzSM+6FI76+EWkJAwMDDAwMMDu3cUZ2evRRx8dMv/II48UJY44qMUvIrHJ1mLv6uri6quvBuDAgQPcfPPNsV9sXrRoEQ888AB9fX1MmjQp8ofZikktfhEpukw1g+LW1NRERUWQEisrKyf0w2xq8YtI0cVdMyjbtQYzA2DatGmsXLnyuPUT5TpD5InfzCqBdcAud7/MzL4HNAK9wO+A97l7b9RxiEjpiqJm0GiVQDMVhTty5AgAPT09R28vHb5ftmMW8kMh6tIRcbT4rwO2A4PDX30PuCac/j7wXuBrMcQhIiUqippBqVSKbds2UDvHj1s3eUrw73hBi/+kGX/NctT9dO7tPG5p114be6AZ3HnnnWzatIk777yTj3/84wU9NkTcx29m9cClwN2Dy9z9QQ8RtPgzF98WkcQYrBlkZgWrGTRY3z8fJ81wTppx/AdFVOfLpKuri9bWVgDWrFlDd3d3QY6bLuoW/23ADcD04SvMrAp4J8E3guOY2XJgOcD8+fOji1BESkIUNYN6e/NrjYcP7jIpz8zYW8DO6jvvvPNo6YiBgYFIWv2RJX4zuwzodPf1ZnZhhk2+Cjzi7o9mWIe73wXcBdDY2Di2j2ARKRuDNYMKJc6yEVC40hFr164dMv/zn/+8fBI/cB5wuZktBaqBE83su+5+jZndDMwB3hfh+UUkwcqxbARw9LbWbPOFEFnid/cbgRsBwhb/R8Kk/17gYuAidx+I6vwiItlku+OnFIZIPeWUU4ZcLzjllFMKfo5i3Md/B/AM8Jvwntkfu/unihCHiMgQcQ6Rmu3D59lnnz1ufviHzXg/gGJJ/O7eBrSF03poTESKqpQfwpo5c+aQO3lmzpxZ8HMoCYuIRGQsQ1DW1dUdTfxmRl1d3XHbpFKpjB9euX4TUOIXEYlIKpViw7bt9M85Oa/9KidVQV8vA1NPYNO+/bnts3dPzsdXkTYRkYiM9aGu/tq5MHkKA9NPiuR8avGLiESp90herXEA+nphYIDK7r1gOT6A1nsk58Mr8YuIRGSkh8iyFYoD+OuRwwBUDvRTXV193Pqamhrq64+vdpPrQ2RK/CIiERnpQmu2C7+9vb1s3boVgP7+fhYsWHDcUJRlcTuniIgMlS1x33TTTUPm586dm3FsgPHQxV0RkRLy8MMPD5lva2sr+DmU+EVESkgctXqU+EVESsjUqVNHnC8EJX4RkRIyWIs/23whKPGLiJSQiy++eMj8JZdcUvBzKPGLiJSQpqamo7dvVlVVFXREskFK/CIiJaS2tpalS5diZlx66aUFGX94ON3HLyJSYqIYfzidEr+ISIkp9PjDw6mrR0QkYZT4RUQSRolfRCRhlPhFRBJGiV9EJGGU+EVEEkaJX0QkYZT4RUQSRolfRCRhlPhFRBJGiV9EJGGU+EVEEkaJX0QkYZT4RUQSRolfRCRhIk/8ZlZpZhvM7P5w/iVm9lszS5nZvWY2OeoYRETkmDgGYrkO2A6cGM7fCnzJ3e8xszuA9wBfG+9JmpubSaVSGde1t7cDsGLFiuPWLVy4MONyEZGJKtIWv5nVA5cCd4fzBvwd8F/hJi3AFVHGAFBTU0NNTU3UpxERKQtRt/hvA24Apofzs4F97t4XzncA8zLtaGbLgeUA8+fPH/VEarWLiOQmsha/mV0GdLr7+rHs7+53uXujuzfOmTOnwNGJiCRXlC3+84DLzWwpUE3Qx/9lYIaZTQpb/fXArghjEBGRYSJr8bv7je5e7+4LgKuBX7j7O4BfAleFmzUB/x1VDCIicrxi3Mf/UeB6M0sR9Pl/owgxiIgkVhy3c+LubUBbOP0U8Jo4zisiIsfTk7siIgmjxC8ikjBK/CIiCWPuXuwYRmVme4FnxnmYWqCrAOGUewxQGnEohmNKIY5SiAFKI45SiAEKE8ep7n7cg1BlkfgLwczWuXtj0mMolTgUQ2nFUQoxlEocpRBD1HGoq0dEJGGU+EVEEiZJif+uYgdAacQApRGHYjimFOIohRigNOIohRggwjgS08cvIiKBJLX4RUQEJX4RkcSZEInfzL5pZp1mtiVt2SwzazWz9vDnzHC5mVlzOObvZjN7VYFiqDaz35nZJjPbamYrw+UZxxg2synhfCpcv6AQcYTH3mlmj5vZRjNbFy6L7fUws5eF5x7894KZ/XPcf5Pw2NeZ2Zbwb/LP4bJI48jyfnxrGMOAmTUO2/7G8Jw7zOzitOWXhMtSZvaxAsXxr+HvttHM1pjZi0b73c2sKXyt2s2sqQAx3GJmu9LeH0uL9FrcmxbDTjPbGGUcWWJ4pZn9Jvz/+lMzOzFtXSSvBQDuXvb/gAuAVwFb0pZ9HvhYOP0x4NZweimwGjDgtcBvCxSDAdPC6Srgt+HxfwhcHS6/A7g2nP4n4I5w+mrg3gK+HjuB2mHLYn090s5bCewGTi3C3+QVwBZgKkFBwp8DC6OOI8v78QzgZQTFChvTlp8JbAKmAC8Bngxfs8pw+jRgcrjNmQWI48S06RVp78GMvzswC3gq/DkznJ45zhhuAT6SYdtYX4th678I3BRlHFlei8eAN4TT7wb+NerXwt0nRovf3R8Bnhu2eBnBmL4wdGzfZcB3PPA/BAPDnFKAGNzdD4SzVeE/J/sYw+nx/RdwkZnZeOMYQayvR5qLgCfd/ZkixHAGQQLr8WDgn4eBN0cdR6b3o7tvd/cdGTZfBtzj7ofd/WkgRVC99jVAyt2fcvcjwD3htuON44W02RMI3qODcWT63S8GWt39OXd/HmgFLhlPDCOI9bUYFP6/exvwgyjjyBLDS4FHwulW4C1RxjBoQiT+LE529z+H07uBk8PpecCf0rbLOu5vvsysMvy62EnwR3yS7GMMH40jXP8XgvEJCsGBNWa23oKxi6EIr0foao79h4o7hi3A+WY228ymErRqX1yEOEaS7ZxRvk8/Y2Z/At4B3FSkOD4Ydil9c7CrrQgxDDof2OPu7UWIYyvHEvdbCd6fkccwkRP/UR58d4r8vlV373f3swmGlHwN8PKoz5nF6939VcAS4ANmdkH6yrheDwuuZ1wO/OfwdXHE4O7bgVuBNcDPgI1Af9xxlBp3/4S7vxj4HvDBIoTwNeB04GzgzwTdLMX09xxrnMTt3cA/mdl6YDpwJI6TTuTEv2fwa3r4szNcvotjn6oQwbi/7r6PYIjJcwnHGM5wrqNxhOtPAroLdP5d4c9O4D6CD6FivB5LgN+7+55wPvYY3P0b7n6Ou18APA/8oRhxjCDbOeOI5Xsc61qILQ533xM2kgaAr3NsYKbYX4vw/96bgXvTFsf5Wjzh7ovd/RyCD58n44hhIif+VQRj+sLQsX1XAf87vIvhtcBf0r72j5mZzTGzGeF0DbAI2E72MYbT47uKYEzicbc8zewEM5s+OA0sJujyiPX1CA1vScUeg5nNDX/OJ/gP/v1ixDGCVcDVFtzl9RKgAfgdwUW/BgvuCptM0GW2arwnM7OGtNllwBNpcWT63R8CFpvZzLBLZnG4bDwxpF83uZLg/TkYQ2yvRehNwBPu3pG2LLY40t6fFcAnCW4AiT6GfK8Gl+I/guTyZ6CXoM/rPQT95WuBdoK7OWaF2xrwFYJP1sdJu8NinDGcBWwANhO8kQfvEDgt/IOlCLo8poTLq8P5VLj+tALFcRrBlf5NBP2HnwiXx/16nEDwDeaktGWxxhAe+1FgW/h6XBRHHFnej1eG04eBPcBDadt/IjznDmBJ2vKlBN9Qnhz8OxYgjh+F78/NwE+BeaP97gTdEanw37sKEMN/hOfYTJC0TinGaxEu/zbw/gzbFzyOLK/FdeHx/gB8jrCaQpSvhburZIOISNJM5K4eERHJQIlfRCRhlPhFRBJGiV9EJGGU+EVEEkaJX8qemd1tZmeO8xhXjPcY5cDMFpjZPxQ7DikuJX4pKeEDRHm9L939ve6+bZynvoKgIuKYpD2dXeoWAEr8CafEL0UXtkJ3mNl3CB4uerGZ/YuZPRYW8hoc2+AEM3vAgjEPtpjZ28PlbWbWaGaX27H66jvM7Olw/Tlm9nBYtO6h4RU3zex1BDWFvhDue7qZnW1m/xOe/760QmLp+33bzO4ws98Cn8+0j5nNDeuwDNZe9/ApYszsSQuKx6Ufc5qZfcuC+uybzewt4fK/D5dtMbNb07Y/kDZ9lZl9Oy22ZjP7tZk9ZWaDT49/jqBw3UYz+/DY/2pSzpT4pVQ0AF91978hqFvfQFDD5WzgHAsKzV0CPOvur3T3VxAUXjvK3Ve5+9keFMrbBPy7mVUBtwNXeVAP5ZvAZ4bt92uCJ0j/Jdz/SeA7wEfd/SyCp0xvzhJ3PfA6d78+0z4e1EuqtmCAjfOBdQSJ91Sg0917hh3v/xKUS/jb8Di/sGCwlFsJSnyfDbzazK7I4TU9BXg9cBlBwodg/IFHw9/zSzkcQyagcvl6KhPfMx7UgYegHsxighIYANMIPggeBb4Ytnjvd/dHMx3IzG4ADrr7V8zsFQQDsrRaMNxBJcFj81mZ2UnADHd/OFzUQoYKo6H/dPf+Ufb5NXAewUAc/0bwAWbh7zPcmwjqrwDg7s+HH3pt7r43jO974bF+MtLvAfzEg0Jo28zs5FG2lQRR4pdS8de0aQM+6+53Dt/IgiEBlwKfNrO17v6pYevfRFDXfLAUtQFb3f3caMIeEnc2jxC09k8lKAj3UYJS0A8U4PzpNVeqh607nDYd5SA/UmbU1SOl6CHg3WY2DcDM5oV95S8Cetz9u8AXCIaxOyrsPvkK8FZ3Pxgu3gHMMbNzw22qzOxvMpxzP0E9dNz9L8DzZnZ+uO6dBKN3ZTXKPo8C1wDtYQv8OYIPr19lOFQr8IG032kmQRG/N5hZrZlVElQ9HTz2HjM7I7wgfuVIMQ7/PSW51OKXkuPua8zsDOA3YffMAYLEuZDgAuwAQYXDa4ft+o8ElTd/Eu73rLsvDS9sNofdMZOA2wgql6a7B/i6ma0gKJPdBNwRXnx9CnhXDqFn3Mfdd1oQ0OAQe78C6j0YynC4TwNfsWBA7n5gpbv/2IJBtX9J0HJ/wN0HS0l/DLgf2Etw/WDaKDFuBvrNbBPwbfXzJ5Oqc4qIJIy6ekREEkaJX0QkYZT4RUQSRolfRCRhlPhFRBJGiV9EJGGU+EVEEub/A1DM6Ex5z0SWAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "\n",
    "releases = []\n",
    "# X axis ticks\n",
    "n_range = range(100, 2001, 200)\n",
    "# Number of samples per boxplot\n",
    "n_simulations = 50\n",
    "\n",
    "for n in n_range:\n",
    "    mean_chain = make_mean_with(data_size=n, epsilon=1.)\n",
    "    for index in range(n_simulations):\n",
    "        # See https://github.com/opendp/opendp/issues/357\n",
    "        random.shuffle(age)\n",
    "        \n",
    "        # get mean of age at the given n\n",
    "        releases.append((n, mean_chain(age)))\n",
    "\n",
    "# get released values\n",
    "df = pd.DataFrame.from_records(releases, columns=['resize to row count', 'DP mean'])\n",
    "\n",
    "# The boxplots show the distribution of releases per n\n",
    "plot = sns.boxplot(x = 'resize to row count', y = 'DP mean', data = df)\n",
    "# The blue line is the true mean\n",
    "plot.axhline(true_mean_age)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "The results from this approach have a similar interpretation as in the prior plot.\n",
    "Underestimated counts lead to higher variance means,\n",
    "and overestimated counts lead to greater bias in means.\n",
    "Thankfully, the count is a low-sensitivity query, so count estimates are usually very close to the true count."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### OpenDP `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\n",
    "terms in the original formula, without any additional modifications, and using the resulting answer as our\n",
    "estimate while ignoring the noise processes of differential privacy. While this 'plug-in' approach does result in a\n",
    "differentially private value, the utility here is generally lower than the solution in OpenDP.  Because the number of\n",
    "terms summed in the numerator does not agree with the value in the denominator, the variance is increased and the\n",
    "resulting distribution becomes both biased and asymmetrical, which is visually noticeable in smaller samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Status:\n",
      "0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%\n"
     ]
    }
   ],
   "source": [
    "n_simulations = 1_000\n",
    "history_plugin = []\n",
    "history_resize = []\n",
    "\n",
    "# sized estimators are more robust to noisy counts, so epsilon is small\n",
    "# the less epsilon provided to this count, the more the result will be biased towards the prior\n",
    "resize_count = make_count_with(epsilon=0.2)\n",
    "\n",
    "# plugin estimators want a much more accurate count\n",
    "plugin_count = make_count_with(epsilon=0.5)\n",
    "plugin_sum = make_sum_with(epsilon=0.5)\n",
    "\n",
    "print(\"Status:\")\n",
    "for i in range(n_simulations):\n",
    "    if i % 100 == 0:\n",
    "        print(f\"{i / n_simulations:.0%} \", end=\"\")\n",
    "    \n",
    "    # See https://github.com/opendp/opendp/issues/357\n",
    "    random.shuffle(age)\n",
    "\n",
    "    history_plugin.append(plugin_sum(age) / plugin_count(age))\n",
    "\n",
    "    resize_mean = make_mean_with(data_size=resize_count(age), epsilon=.8)\n",
    "    history_resize.append(resize_mean(age))\n",
    "    \n",
    "print('100%')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABBVUlEQVR4nO3deXzU9bX4/9fJZF/IHrYAAQRkDxrBDeuOW9W2LnivttZabdXa6m3vz/be9lpv7fXb2tbr0loVL27FBdeq1B3FjQICYd+3BEhCIHsyycy8f398PpmZJJNkksxMMsl5Ph558NnnMIQ5897FGINSSqmhK6a/A1BKKdW/NBEopdQQp4lAKaWGOE0ESik1xGkiUEqpIS62vwPoqZycHFNQUNDfYSilVFRZs2bNEWNMbqBzUZcICgoKWL16dX+HoZRSUUVE9nV2TquGlFJqiNNEoJRSQ5wmAqWUGuKiro1AKTWwtLS0UFJSQlNTU3+HooDExETy8/OJi4sL+p6wJQIRGQM8DQwHDPCYMeZ/211zJvA6sMc+9Iox5p5wxaSUCr2SkhLS0tIoKChARPo7nCHNGENlZSUlJSWMHz8+6PvCWSJwAf9mjPlKRNKANSLynjFmc7vrVhhjLgljHEqpMGpqatIkMECICNnZ2VRUVPTovrC1ERhjDhljvrK3a4EtwOhwvZ5Sqv9oEhg4evNvEZE2AhEpAOYAKwOcPkVE1gMHgZ8aYzYFuP8m4CaAsWPHhjFSpdraX9mAxxgKclL6O5SoUHDXW2F79t77Lg7bs4e6sPcaEpFU4GXgJ8aYmnanvwLGGWNmAw8BrwV6hjHmMWNMkTGmKDc34MA4pULu9XWlnP2H5Zz9h+W8uPpAf4ejuuBwOCgsLGTGjBl8/etfp6qqqsfPWL16NbfffnufYxERrr32Wu++y+UiNzeXSy4ZuDXgYU0EIhKHlQSeM8a80v68MabGGFNnb78NxIlITjhjUioYmw5W87OXinF5DB4Dv3hlA+sPVPV3WKoTSUlJrFu3jo0bN5KVlcUjjzzS42cUFRXx4IMP9jmWlJQUNm7cSGNjIwDvvfceo0cP7FrxcPYaEmARsMUY88dOrhkBlBljjIjMxUpMleGKSalgPfPFPprdHu++y2NY/Ple/nR1Yf8FFWWWfP/kPj/jmse/7PE9p5xyCsXFxQDs2rWLW2+9lYqKCpKTk3n88cc5/vjjeemll/j1r3+Nw+EgPT2dTz75hOXLl3P//ffz5ptvctFFF3Hw4EEA9uzZw4MPPsi1117LXXfdxfLly3E6ndx6663cfPPNAWO46KKLeOutt7jiiitYsmQJ11xzDStWrACgvr6eH/3oR2zcuJGWlhbuvvtuLrvsMvbu3ct1111HfX09AA8//DCnnnoqy5cv5+677yYnJ4eNGzdy4okn8uyzz4a0XSacbQSnAdcBG0RknX3sF8BYAGPMo8AVwA9FxAU0AguNrp2p+pkxho+2lXc4/vH2CjweQ0yMNowOVG63mw8++IDvfe97ANx00008+uijTJo0iZUrV3LLLbfw4Ycfcs899/DOO+8wevTogNVIb7/9NgBr1qzhu9/9LpdffjmLFi0iPT2dVatW4XQ6Oe200zj//PMDdtNcuHAh99xzD5dccgnFxcXccMMN3kRw7733cvbZZ/Pkk09SVVXF3LlzOffcc8nLy+O9994jMTGRHTt2cM0113jnVVu7di2bNm1i1KhRnHbaaXz22WecfvrpIXvfwpYIjDGfAl3+jzHGPAw8HK4YlOqNTQdrKKtxApCS4EAQ6pwujtY3U1xaTeGYjP4NUHXQ2NhIYWEhpaWlTJ06lfPOO4+6ujo+//xzrrzySu91Tqf173raaadx/fXXc9VVV/HNb34z4DOPHDnCddddx4svvkh6ejrvvvsuxcXFLF26FIDq6mp27NgRMBHMmjWLvXv3smTJEi666KI25959913eeOMN7r//fsDqfrt//35GjRrFbbfdxrp163A4HGzfvt17z9y5c8nPzwegsLCQvXv3RkciUCpafbzd1we7cEwmMQIrdhwBYPm2ck0EA1BrG0FDQwMLFizgkUce4frrrycjI4N169Z1uP7RRx9l5cqVvPXWW5x44omsWbOmzXm3283ChQv51a9+xYwZMwCrpPjQQw+xYMGCoGK69NJL+elPf8ry5cuprPTVeBtjePnll5kyZUqb6++++26GDx/O+vXr8Xg8JCYmes8lJCR4tx0OBy6XK6gYgqWJQKl2NpRUe7dnjBqGI0a8iWBjaXVnt6l2elO/31fJyck8+OCDXH755dxyyy2MHz+el156iSuvvBJjDMXFxcyePZtdu3Yxb9485s2bx7JlyzhwoG2vsLvuuotZs2axcOFC77EFCxbwl7/8hbPPPpu4uDi2b9/O6NGjSUkJ3LX4hhtuICMjg5kzZ7J8+fI2z3nooYd46KGHEBHWrl3LnDlzqK6uJj8/n5iYGJ566incbndY3qNAdNI5pdrZctjXy3lcdgrjsn3/0bccqu2PkFQPzJkzh1mzZrFkyRKee+45Fi1axOzZs5k+fTqvv/46AD/72c+YOXMmM2bM4NRTT2X27NltnnH//ffz7rvvUlhYSGFhIW+88QY33ngj06ZN44QTTmDGjBncfPPNXX4zz8/PD9gd9Ze//CUtLS3MmjWL6dOn88tf/hKAW265haeeeorZs2ezdevWThNMOEi0tc0WFRUZXZhGhUud08WM/3oHgBiBxd+diwDXL16F22P9Xym++3yGJQY/oddgt2XLFqZOnQrogLKBwv/fpJWIrDHGFAW6XquGlPKz7bDvG/+ojCTiHDHe7QNHGwDYfriWooKsfolvoNMP6+ikVUNK+dlyyK9aKCvZuz3Wb3vLYa0eUoOLJgKl/PiXCMZ2kgi2HW4/U4pS0U0TgVJ+9lbWe7dHZST5bfu68u2rbIhoTEqFmyYCpfzsP+r7kB8+zPfhPzwtMeA1Sg0G2lislK3F7aH0WKN3P29YQsDt0mONuNweYh36PaqDu9PD+GwdwxEu+puslO1gVSMuu4toRnIcCbEO77mEWAcZSVaXUZfHcKha1+cdSPynob7yyitpaLBKbampqSF7jRtvvJHNm9svsNi5xYsXIyK8//773mOvvfYaIuKdpmKg0ESglM2/7n+EX7VQK/9SgVYPDSz+01DHx8fz6KOPhvw1nnjiCaZNm9aje2bOnMnzzz/v3V+yZEmHwWsDgVYNKWXb10n7QKu8tES2l9UBViI4LWKRRanvvNn3ZzzV88Vc5s+f752GupX/FNMAt912G0VFRVx//fW8/fbb3HnnnaSkpHDaaaexe/du73X+zjzzTO6//36KiopITU3lxz/+MW+++SZJSUm8/vrrDB8+PGAsK1asoKWlBafTyc6dOyksLPSeX7NmDXfeeSd1dXXk5OSwePFiRo4cyeOPP85jjz1Gc3Mzxx13HM888wzJyclcf/31DBs2jNWrV3P48GF+97vfccUVV/T4PWpPSwRK2fb79RjKS0vocF5LBAOfy+Vi2bJlzJw5M6jrm5qauPnmm1m2bBlr1qwJetH3+vp6Tj75ZNavX88ZZ5zB448/HvA6EeHcc8/lnXfe4fXXX+fSSy/1nmtpaeFHP/oRS5cuZc2aNdxwww38x3/8BwDf/OY3WbVqFevXr2fq1KksWrTIe9+hQ4f49NNPefPNN7nrrruCirc7mgiUspX4NRR3ViJopYlgYGmdhrqoqIixY8d61yPoztatW5kwYYJ3KulrrrkmqPvi4+O9S0+eeOKJ7N27t9NrFy5cyPPPP8/zzz/f5vnbtm1j48aNnHfeeRQWFvKb3/yGkpISADZu3Mj8+fOZOXMmzz33HJs2+ZZyv/zyy4mJiWHatGmUlZUFFW93tGpIKdvBKl8iyE6N73A+x+/YYW0sHlBa2wg6Exsbi8fjW3Guqan7f78FCxZQVlZGUVERTzzxRJtzcXFx3hXCupsWeu7cuWzYsIHk5GQmT57sPW6MYfr06XzxxRcd7rn++ut57bXXmD17NosXL24ze6n/lNShmitOE4FSttIq34dDTmrHqqHsFN8xTQRB6EX9friMGzeOzZs343Q6aWxs5IMPPuD0009nypQp7N69m71791JQUMALL7zgveedd94J2evfd999bdYXAJgyZQoVFRV88cUXnHLKKbS0tLB9+3amT59ObW0tI0eOpKWlheeeey7sax5rIlAKcLrcHKmzVq8SgczkjiWCrBS/EkFNE26PwaHLVkaFMWPGcNVVVzFjxgzGjx/PnDlzAKsk8ec//5kLLriAlJQUTjrppLC8/oUXXtjhWHx8PEuXLuX222+nuroal8vFT37yE6ZPn85///d/M2/ePHJzc5k3bx61teGd30qnoVYK2FdZz9d+vxyA7JR4Hv6XEwJe9/2nV1PntKoBVv7inIBtCUNNmymPo3BAWV1dHampqRhjuPXWW5k0aRJ33HFHWF4rUnQaaqV6odSvfSBQtVCr7NR4byI4VN2kiaC9KBz9+/jjj/PUU0/R3NzMnDlzuPnmm/s7pIjTRKAUcMivfSArQENxq+yUeO/As0NVjbp+8SBwxx13RH0JoK+0+6hStO0xlJPSeSLI8msw1mkmfKKtinkw682/hSYCpYCD1UFWDbVrMFaQmJhIZWWlJoMBwBhDZWVlhx5K3dGqIaVo++0+q8sSge+cfyliKMvPz6ekpCToUbkqvBITE8nPz+/RPZoIlKLtuIDMIBNBeY0zrDFFi7i4OO/IXBWdtGpIKaCsJrgSQUZynHe7vFarhtTgoIlADXlOl5tjDS0AxAikJ8Z1eq3/QLPyWqfWi6tBQROBGvL8q3gykuOJ6WK0cHK8gziHdb6h2e0dU6BUNNNEoIY8/94/mcmdlwbAmla4falAqWiniUANeYeD7DHUqk0i0AZjNQhoIlBDXlmbEkH3iUAbjNVgo4lADXnB9hhqlaElAjXIaCJQQ95hvw9zLRGooShsiUBExojIRyKyWUQ2iciPA1wjIvKgiOwUkWIRCTz3r1Jh1KZqqIdtBGVaIlCDQDhHFruAfzPGfCUiacAaEXnPGLPZ75oLgUn2zzzgL/afSkVMeQ96DbW/RksEajAIW4nAGHPIGPOVvV0LbAHar7d2GfC0sXwJZIjIyHDFpFR7xpg2XUCDqRpKT/IlgiN1zWGJS6lIikgbgYgUAHOAle1OjQYO+O2X0DFZICI3ichqEVmtE1upUKpzumhodgMQ5xCS4x0drnG6mnC6fN/8/RuLK3QcgRoEwp4IRCQVeBn4iTGmpjfPMMY8ZowpMsYU5ebmhjZANaS1Lw2IdBxVnBCbSEKsb1rftIRYWi+rbmzB6XKHPU6lwimsiUBE4rCSwHPGmFcCXFIKjPHbz7ePKRURwYwheG/3Ut7bvdS7HxMjbeYjqtTqIRXlwtlrSIBFwBZjzB87uewN4Nt276GTgWpjzKFwxaRUe/5VOxmdNBSvLH2flaXvtzmW7netVg+paBfOXkOnAdcBG0RknX3sF8BYAGPMo8DbwEXATqAB+G4Y41Gqg552HW3VtsFYE4GKbmFLBMaYT4HOp3G0rjHAreGKQanu+I8Mzkzqvutoq4wkLRGowUNHFqshrcy/sbgHJQLtOaQGE00EakjzrxrKCGIMQSutGlKDia5ZrIa0ijbdRwNXDf3n/Ec7HPNPBBWaCFSU0xKBGtLKe1kiyNBeQ2oQ0USghqw6p4t6v1HFKQFGFQO8teNZ3trxbJtjOs2EGkw0Eaghq7zdYLJAo4oB1h7+lLWHP21zrE0i0BKBinKaCNSQVVbT/WCyzqQkxOKwE0et00VTi04zoaKXJgI1ZPlPIR3MrKP+YkQYluTra6E9h1Q000Sghqy200v0LBGAthOowUO7j6ohqyzIBWniHQkBj2s7gRosNBGoIas8yBLBv5/6vwGP66AyNVho1ZAasoItEXTGP3loIlDRTBOBGrKCXaLy1a2LeHXrog7HtY1ADRaaCNSQ1Wbm0S4mnNtUsYpNFas6HNdpJtRgoYlADUl1Thd1ThfQ9ajirmhjsRosNBGoIan9EpWdjSruijYWq8FCE4EakoJZq7g7ulylGiy0+6gaktq2D3TdYygtPj3g8dSEWGIEPAZqmlw4XW4SYntexaRUf9NEoIaknpQIfjzv/wU8HiNCelIcxxpaAKvn0OiMpNAFqVSEaNWQGpL8J5zrbdUQaIOxGhw0Eaghqcx/wrlu1ip+YdMjvLDpkYDntMFYDQZaNaSGpPIejCrecXRDp+d0EXs1GGiJQA1Jh0PQawi0RKAGB00EasgxxlBW7fvQzuqmaqgrbUYXa4lARSlNBGrIOVrfTLPbA0ByvIPEuN53+dT5htRgoG0EasjpabVQVlJep+cydFCZGgQ0Eagh53C1LxFkB1EtdEvRPZ2e04nn1GCgVUNqyGlTIuhD+wBARpL2GlLRTxOBGnLK/EoEwTQUP1P8R54p/mPAcykJDmJjrAnr6pwuGppdoQlSqQjSRKCGHP8SQTCJYF/1dvZVbw94TkS0nUBFPU0Easg55F8i6MMYglY6qExFO00EasgpC2EbAUCGX4NxuSYCFYXClghE5EkRKReRjZ2cP1NEqkVknf3zq3DFopS/Q1V+vYZSQ1Ei0KohFd3C2X10MfAw8HQX16wwxlwSxhiUaqOmqYVavyUq0xK6/y8wMnVsl+fT/XoOlftNZqdUtAhbIjDGfCIiBeF6vlK90aY0kJIQ1BKV35vziy7Pa4lARbv+biM4RUTWi8gyEZnez7GoIeBgdaN3OxTVQtA2EWgbgYpGQSUCEXlFRC4WkVAmjq+AccaY2cBDwGtdvP5NIrJaRFZXVFSEMAQ11LQtEQSXCBat/S2L1v620/M6qExFu2A/2P8M/AuwQ0TuE5EpfX1hY0yNMabO3n4biBORnE6ufcwYU2SMKcrNze3rS6sh7FCbEkFCcPfU7edQ3f5Oz/uvZ+C/8plS0SKoRGCMed8Y86/ACcBe4H0R+VxEvisiXa/q0QkRGSF2Ba2IzLVjqezNs5QK1sFelAi6k+6XCCrrnbjsmU2VihZBNxaLSDZwLXAdsBZ4Djgd+A5wZoDrl9jHc0SkBPgvIA7AGPMocAXwQxFxAY3AQmOM6cPfRaluHazqeYmgO7ExMQxLiqOmsQVjrOmoR6QnhuTZSkVCUIlARF4FpgDPAF83xhyyT70gIqsD3WOMuaarZxpjHsbqXqpUxLSpGgpRiQAgK9lKBGANWNNEoKJJsCWCx+16fC8RSTDGOI0xRWGIS6mQ83hMm6qhnCBLBOPSJ3d7TWZyPHsrG4C2I5eVigbBJoLfAG+3O/YFVpuBUlHhSJ3TuzJZSoKDpPjgVia7btad3V7jP1WFJgIVbbpMBCIyAhgNJInIHKB19M0wIDnMsSkVUiV+7QO5IWofaKU9h1Q0665EsAC4HsgH/CdkrwW6Hm6p1ABTesyXCIKtFgL482prGqyuVirzX/JSSwQq2nSZCIwxTwFPici3jDEvRygmpcKixD8RpAWfCI42lnd7TZtEoIPKVJTprmroWmPMs0CBiHSoKDXGBF62SakBqLSqwbsdTNWQuJ1k71tGctVOxLQwbvW9VI67iLrcOR2u9W8jKNcSgYoy3VUNpdh/poY7EKXCrSdVQ/H1pRz/0c2kHNtMXFI6AKO2LGLUlkUcnnQN+078BZ64FO/1/m0E/gvfKBUNuqsa+qv9568jE45S4VPq31jcRdVQQu0+ZvzjSuKbjgQ8P2LHEpKrd7Dl7MV44qw+E8OS4nDECG6PobqxhcZmd9C9kpTqb8FOOvc7ERkmInEi8oGIVIjIteEOTqlQMca0bSPoZOZRcTczecWPvUnAiIPJqeOZmF1ITe6J3uuGla9mysc/QDzWILIYkTbLXvoPXFNqoAt20rnzjTE1wCVYcw0dB/wsXEEpFWqV9c00NLsBSIpzkNrJgjRj1v+J1MpiADziYO8Jd3Fl4b+zcOrNHJh9J4cn/6v32oxDn5Jf7Bsc7z+ttVYPqWgSbCJo/V9zMfCSMaY6TPEoFRYHjvoaivOGBV6QJqF2HyM3L/Lul09aSEOW3zIZIlSOu5jy8d/wHhq98RHSylcBkOXXYOw/p5FSA12wieBNEdkKnAh8ICK5gH7lUVFjv38i6KR9YMz6B4gx1jKW9RlTqBx7IQD/s20R/7PNlyAqJn6L+sypAIjxMPGLnyPu5jZzFx3WEoGKIsFOQ30XcCpQZIxpAeqBy8IZmFKh1KZEkNZxQrikY9vI2fOGd79s0jVgr8NU66qn1lXvu1hiKJlxC+7YJOvemt2M2PY0WSm+BHNQE4GKIj1Zcex44GoR+TbWFNLnhyckpULvwFFfVU2gEsGoLYsQrFnQa3Pm0JjR9URzrsRsKiZ807ufX/wg+fG13v3D2lisokiwvYaeAe7HWn/gJPtHZx1VUWN/uzYCf7HOY+Ts/bt3v2J8cIXdo2MW4EwZZT2jpY7TDz/rPaeNxSqaBDv7aBEwTReOUdHqwLHOq4Zydy4lxm1NC9GYVkBj+qSgnmliYik77hrGrv8DAJMPvEguJ1NBpjYWq6gSbNXQRmBEOANRKlyaXR7vB7PQblSxMYzY8Tfv7tEx50O7HkWz0iczq5M1CWpzT6AxbTwADo+TW+OskkVNk4uappYQ/i2UCp9gE0EOsFlE3hGRN1p/whmYUqFScqwBj12WzU6NJz7W92ufVrGGxNp9ALhjk6kecUqH+xfmX8DC/AsCP1yE8olXeHevdnxEBlZbgf+UFkoNZMFWDd0dziCUCqd9lb5qoeHD2lYL5ex53btdPfxkjKPn6xTU5RTSmDaOpNp9JOHkWsf7POz+BqXHGpk6cljvA1cqQoLtPvox1ojiOHt7FfBVGONSKmT2Vvq6fo7wSwTibiZ771ve/eqRpwW8/+4tf+HuLX/p/AXsgWatvhP7Dgk0U+LXLqHUQBZsr6HvA0uBv9qHRgOvhSkmpULKv0SQ55cIMg6uIK65CoDmxGwaMqYEvL/Z00Kzp+v6/urhJ9OSkAVArtRwqePzNnMbKTWQBdtGcCtwGlADYIzZAeSFKyilQqmzEkH2ft8y3NUjTvUOIOuVmFgqxy7w7n7b8S6lWiJQUSLY33ynMaa5dUdEYgHtSqqiQts2AqsNQNzNZJZ84D1eM/zkPr9O1agzcYvV7DYzZi/JFev6/EylIiHYRPCxiPwCaxH784CXgL93c49S/a7F7WkzvURrY/GwspXENtcA0JyYQ1NaQZ9fyx2fRmWuL6GcVfN6F1crNXAEmwjuAiqADcDNwNvAf4YrKKVCpeRYIy6772hWSjyJcdZiMVn7/+G9pibvpA5jB/ydlDmdkzKnd3reX80438wr55ovqK2q6E3YSkVUUN1HjTEeEXkNeM0Yo7/ZKmrsOVLn3R6ZbrcPeNxkHXjPe7w276Qun/GNUecE/XrOjOPYzyjGcpBEaeHoP5eQdv7tPQtaqQjrskQglrtF5AiwDdhmr072q8iEp1Tf7K7wNRS3JoLUyvXeFchc8cNo6GaCuZ5alzTXu52y6TnQmVnUANdd1dAdWL2FTjLGZBljsoB5wGkickfYo1Oqj3Yf8U8E1rTR/o3EtTlzuu0t9ItND/KLTQ8G/ZoH0k/CaazCdnr1Vji4tichKxVx3SWC64BrjDF7Wg8YY3YD1wLfDmdgSoXC7oqOVUOZJR96j9XmnhDy10xJS+dLz1TfgXV/6/xipQaA7hJBnDHmSPuDdjtBXHhCUip09viVCEZlJBFfV0pK1TYAPBJLfdbMkL9mZnIcH3tm+w5seAlczpC/jlKh0l0iaO7lOaX6XW1TC2U11gewI0bISU0gy69aqD5rOp7YjquV9VVmcjxbzDjKTbp1oKkKti0L+esoFSrdJYLZIlIT4KcWCP1XKaVCaFe7hmJHjJBR+pH3WDiqhQAykuMwCJ94ZvkOavWQGsC67D5qjHFEKhClQm1nua99YHRGEjEtDaQf/sJ7rC6nMKjnnJ49p0ev64iJISs5nk8aZnGFY4UdzPtQVwGpuT16llKR0IfJVbomIk+KSLmIbOzkvIjIgyKyU0SKRSQ8X8/UkLWj3LeG8OiMJNIPf0GMx6rRbErJpyUpuA/li0bM56IR83v02lkp8ZSTyVbPGOuAccPGl3v0DKUiJWyJAFgMdLKaBwAXApPsn5uALub5VarndvmVCEZlJJFR6ustVJcb/Ld8p7sZp7tnTWKtq6Ct8PjVoBY/36NnKBUpYUsExphPgKNdXHIZ8LSxfAlkiMjIcMWjhp62VUOJZJYu9+7XBlktBPDrrY/y662P9ui1s1PjAfjSMxU3dg3rwbVQsb1Hz1EqEsJZIujOaOCA336JfawDEblJRFaLyOqKCp3hQnWvqcXNfnuyOQEmcYCEhkOAtSRlQydrEIdKjp0I6kliS+zxvhPFL4T1dZXqjf5MBEEzxjxmjCkyxhTl5mpjm+reroo67zrFuWkJ5B5e7j1Xlz0LYsLbDyI7xbfk5QdOvwnrNrykU06oAac/E0EpMMZvP98+plSfbS/zNRSPzUom07/baA+qhXorOd5Bkj3T6Sr3BDxxKdaJqn1w4J9hf32leqI/E8EbwLft3kMnA9XGmEP9GI8aRLYd9rUPTEl3kVZhLbFtkKC7jfaFiJCbZpUKXMRSmeM3w6lWD6kBJpzdR5cAXwBTRKRERL4nIj8QkR/Yl7wN7AZ2Ao8Dt4QrFjX0+JcITjXrEOMBoDF9Iu74YT161jm58zgnd16PY8hL81UPbU326x296VVw6cB8NXAEtR5BbxhjrunmvMFaC1mpkNt22JcIZtR/6d3uTWngnLyeJwGAvDTf9BVrnaOZn5wDDUeg8Sjs/ggmL+jibqUiJyoai5XqidqmFkqrGgGIj/Ew6sinvnM5PR+3WNNSR01LXfcXtpM3zFciKDnmhPFn+E4Wv9jj5ykVLpoI1KCz1a80sCBtL3HN1QC0JGTRlDaux8+7b/uT3Lf9yR7fl+tXNVRa3UhLwZm+k9veBmfPk4tS4aCJQA06Ww7VeLcviFvn3a7NPaHLtYlDLTHWQXqiVfvq9hhKGQ4ZY62TLQ2w9a2IxaJUVzQRqEFn80FfIpjX7GsfCNdso10ZYa+KBrCnsh4mnOk7uUGrh9TAoIlADTqtJYKJUkpOcwkAbkcC9ZnTIh7LiHRfg/G+ygYYf6bv5K6PrBlJlepnmgjUoOJye7xtBAtiVnuP12fPwjjiIx6PfyLYc6QeUvMgz05Ixg2bXol4TEq1p4lADSq7j9TjdFljBi6K8yWCmryTOrulWxcOP50Lh5/eq3tH+iWCvZX1eIyBCWf5LtDeQ2oA0ESgBpUNJVYPoZFUMoNdABhx9Glaifk5JzC/F91OAdIS4kiNtxqMnS4PB6saoeB0iLGH8JSuhspdvY5NqVDQRKAGlQ2lViI43+FXLZQ5DU9caq+fWeE8RoXzWK/vH5nhazDeWV4HCWkwush3gZYKVD/TRKAGlY12IrjIsdJ7rC/VQgB/2vkMf9r5TK/vH+2fCCrssQP+vYeKX9AZSVW/0kSgBg23x7DpYA15HOMk2QZYk8zV5hV1c2d4jc5sVyIAGDMX4pKt7WN7oGRVP0SmlEUTgRo0dlfU0dji5mLHl8SI9Q27PmsaroSMfo1rVIavwXj/0UarMdsRb7UVtFqvy1iq/qOJQA0aaw9UAXCJwzeIrGb4yf0UjU9irIMce6EajzHs8lYP+fUe2vSKzkiq+o0mAjVorDtQxWgqODFmB2D1Fupr+0CojMnyVQ9tO2yPfB4+HVLsFfcaj8GOd/shMqU0EahBZN3+Ki5zfObdr8ua0eO1BwK5fORZXD7yrO4v7MKYrGTv9tbWRXMkBiac7bto/ZI+vYZSvaWJQA0Kjc1utpXV8C3HCu+xqpHzQ/LsuVkzmZs1s0/PGOuXCLaX1eJu7SU00S/BbH8HGo726XWU6g1NBGpQKC6pYqbZycQYa7VTtyOJ2rwTQ/LsksYyShrL+vSMjKQ4UhOsQWSNLW72VzZYJ9LzIWeyte1pgY0v9+l1lOoNTQRqUFi97xjfcnzi3a8ZPg/jSOjijuD9efcL/Hl339YZFhHGZvtKBZv9pspm4jm+7XV/69PrKNUbmgjUoLBhz8E27QOhqhYKpfFZKd7t1oFv1okzfFNOHPwKyrdGODI11GkiUFHP4zEM3/82w8RanrIhcTgNmcf3c1QdFeT4EsHWw7W4PXY7QUIajPFbF3ndcxGOTA11mghU1NtVUcc3PL6ul9VjzonoSmTBykyOIy3R106w64jfUpX+1UPFL4DbFeHo1FCmiUBFvR3rP6MwxprB04WD6lFndHNH/xARJuT4Jr8rPuBXPTT6REjKtLbrymDn+xGOTg1lmghU1Mva+H/e7X2phSEZO+DvqtELuGr0gpA8a4Jf9VCxfztBjAMm+o0pWNv7Se6U6ilNBCqqmdoyTqjxfXs+MuaCkL9GYcYUCjOmhORZ43N9iWBHWS11Tr8qoOPO9W1v/4cuY6kiRhOBimrVKx4jHuvDdKcZTfLo0K9LvLu+hN31JSF5Vkp8LKPsVcsMUFziVypIHwO5diO3xwXFOhGdigxNBCp6NdeTuG6Rd3d18hnEhKGR+Im9r/DE3tCtLXxcXpp3+6t97UYSTzrft73mKV2nQEWEJgIVvb56msRma+WwCpNO7Yi5/RxQcCYP9yWCtQeqfd1IAQrmQ6w9QV3lDtj/JUqFmyYCFZ1czZjPHvTuvuE+hYKc0DYSh8uIYQnebqT1zS62ldX6TsYlWQPMWn31VISjU0ORJgIVndY9i9QeBKDKpPC5zGFEemI3Nw0MIsJkv+qhlbsr214w2a+H0qZXdSI6FXaaCFT0aWmCj3/v3X3TfTKjstPD0j4QLlNH+kovK/ccxePfFpA9CbImWtuuJl29TIWdJgIVfdb8H/iVBt71FFGQndLNTb133dhLuG7sJSF95rjsZJLjHQBUNbawrcxvlLEITLnQt7/m/7TRWIWVJgIVXZpq4JP7vbuvuk+nmTjG54QvEUxNm8DUtAkhfWaMCMf7NRp/urPdmIHxZ1jtBQBHtsPeFSgVLpoIVHT57AFoOAJAhRnGh545JMU7yE0LzZTTgWyp3c2W2t0hf+7M0ene7S93HaXZ7fGdjEtuu3rZyr+G/PWVahXWRCAiF4jINhHZKSJ3BTh/vYhUiMg6++fGcMajolx1CXzxiHf3efdZtBBLQVZyWNsHntn/Js/sfzPkz83PSiY9MQ6weg+t3X+s7QXHX+zb3vY2VO0PeQxKQRgTgYg4gEeAC4FpwDUiEmjY5wvGmEL754lwxaMGgXf/02o8BQ7H5vO5ZwYABX4TuUWTGBFm5ftKBR9ubVc9lDEWRs62to0HVi1CqXAIZ4lgLrDTGLPbGNMMPA9cFsbXU4PZ7o+trpS2J1vOxWCVAsbnJHd214BXOCbDu72+pIryWmfbC47/um97zWJw1qFUqIUzEYwGDvjtl9jH2vuWiBSLyFIRGRPoQSJyk4isFpHVFRU6EdeQ42qGZf/u3a0ZeRrFLfkADEuMJSs5vr8i67OM5Pg2M5K+v6Xd2sj5J0HaSGu7qUoXrVFh0d+NxX8HCowxs4D3gIDDKI0xjxljiowxRbm5uRENUA0An/0vVNjLN8Ym8Wm671vy+JxUJIrGDwRSVJDp3f5gSzlNLrfvZIwDpl3u2//iEV20RoVcOBNBKeD/DT/fPuZljKk0xrSWhZ8ATgxjPCoaHdkJn/gGjzHnWlZX+D74C7LDXy10Y8E3ubHgm2F7/qS8NDKSfI3Gn2xvV+o97hxIsAegVe1rU0WmVCiEMxGsAiaJyHgRiQcWAm/4XyAiI/12LwW2hDEeFW08bnj9VnDb3xWyJ+E87iK2+w2+Kgjj+IFWE1LymZCSH7bnx4gwb3yWd//vxYfaTkQXm9i2B9EnvwePX1dTpfoobInAGOMCbgPewfqAf9EYs0lE7hGRS+3LbheRTSKyHrgduD5c8agotPJROGDPvikOOPVHbKuox2V/SOakJDDM7n4ZTuuqtrGualtYX6NwTCZJcdZI44paJ5/vPtL2gqmX+g0w2wZbXg9rPGpoCWsbgTHmbWPMZGPMRGPMvfaxXxlj3rC3f26MmW6MmW2MOcsYszWc8agoUr4VPrjHtz/rKsiawAa/5R0LItRb6MXSd3ix9J2wvkZ8bAxz/UoFr6wpxe0/rURCWtseRB/9j7YVqJDp78ZipTpyNcMr3/eOGSBzPMy8CoCNpTXeyyZEoFookuYWZJEQa/2XPFTTxOe72pUKpl3mW6vgyDbtQaRCRhOBGng++g0cLra2Y+Jg/r+BI446p4s9R+q9l40L40Rz/SExztGmrWDpmhJvNZh1QTrM+JZv/6PfQnM9SvWVJgI1sOz8wOou2uqE70BmAQCbDvqqhUZnJJFo16kPJvPGZ5NolwrKapx83L4H0bTLIclOFnWH20zAp1RvaSJQA0dtGbz6A9/+qBNg2qXe3WK/9oHxEeg22h8S4xycOjHHu//yVyXtJqNLhDnX+fY/fxDKtbOd6htNBGpg8Ljh5e9Bfbm1n5gBp98B4vsVLS7xSwS5kZtf6JYJV3PLhKsj9nonFWSREm8tZXm0vpn3NrcbbXzcOZBnT9vlccHff6wNx6pPNBGogeGj3/rNuS9wxk8hyTfi9nBNExX2PDxxDmFMZlLEQstPGk5+0vCIvV58bAynH+crFby+rpTGFr/RxhIDJ99idakFOLASPv1jxOJTg48mAtX/ti2DFX513bMXwsjCNpf4lwbGZaXgiIncr+4/j27gn0c3ROz1AE4Yl0G6vcB9TZOLf2w63PaCzAKYfY1vf/l9sPezyAWoBhVNBKp/Ve6CV2727Y+aA7MWdrhs/QHfXP0TcyPbW+i1Qx/x2qGPIvqasTExzJ/km1frzfWHqG9uV/0z80pfFZFxw4vXwbG9kQtSDRqaCFT/cdbC8/8CTvvbfkouzP+pNdGan2a3hw1+4weOy0tjKJg1Jp3MZN8cRMs2tCsVxDjgjJ9Z7SkADZXw7BVQVx7ZQFXU00Sg+ofHA6/90DeraEwcnPUfVl/5drYcqvH2nMlKjicrJXqnne4Jh8Rwhl+p4K2Nh6hrXypIybXetxirGonKHfD0ZVCn07Wr4GkiUP3jk9/Dlr/79k+5DbKPC3jpV/v8qoXyonM1st6aMTrdm/gam928XXyo40V5U62SVGsPq/LN8OT5cHRPBCNV0UwTgYq8za/D8t/69qdeanWJDMBgWOWXCKYMHxrVQq1iRPiaX6lg2cbD1DkDdBUtOB1O8+tue3Q3PHEO7FnR8Vql2tFEoCLr4Lq2jcMjZ0PR9zq9fFdFPUfrmwFIjIthbHbkuo22uuO467jjuOu6vzBMpo0aRk5KAgCNLW7+vv5g4AsnngVfu8uqZgOrzeDpy6yR2jptteqCJgIVOTUHYclCcDVa+2kj7Q+uzqeKWLmn0rs9KS8Nh0T+VzY3IZPchMzuLwyTGBHOmOwrFfxj02GqGlsCXzzuVFhwr68B2bjhvV/Bkqu13UB1ShOBigxnLfztaqi167jjUuCcX1nTK3fCbQyf7fAlgmkjh4U7yoBWHPmKFUe+6pfXbjV1ZBp5aVapwOny8Nra0s4vzpsGlzwAuVN8x3a8C385FXa+H95AVVTSRKDCz90CL13vm1FUHHDWzyF9TJe3bT5Yw9EGq1ooOd7BxLz+mW10WdmnLCv7tF9eu1WMCGdNyfPuv7eljLIaZ+c3pOTAgvtg2jd8x+rL4dlvwT9+Di1NYYxWRRtNBCq8jIG//6TtN9GTb+kwcjiQ97f45tiZMSq9X6qFBpJJeamMybQm23N7DEv+ub/rGxxxcNL34Nxf+6qKAL78s9WQrJPVKdvQ/p+lwu+9X8G6Z337sxbC5AXd3lZe6+Sfe4569wvHZoQhuOgiIpwz1Vcq+HJPJZsP1XRxh230iXDpw5B/ku9Y2UZ47Ez45+NWslZDmiYCFT4f/86aJrnVcedC4b8Gdetr60pp/Xgan5PC8LTE0McXhcZkJjN9lK+t5MnP9tDiDuKDPCkDzv4VzPshOOwBea4mePun8LerdDTyEKeJQIXHij/AR/f69secDKf8CES6vXVPZT0fbvV9MJ0yITscEUatc48fTpzDeh9LjjXy5oZOupO2JwLHXwwX/8m72A9gNST/+WTY/Ebog1VRQROBCi1j4MN72y48P2oOfO3fu+wm2qq+2cWDH+zw7k/MSWFiBNceCOSuyTdw1+Qb+jUGf8OS4jhzsq+K6OU1Jew/2hD8AzLHwcV/tFY7a9VQaU1at/R7UF/Z6a1qcNJEoELHbS+S8snvfMdGzLLmwnF0Pz/QvqMN/PqNzRyqtnq0xMYI508fEa5ogzYsLpVhcQNraou547MYlW5Vl7k8hoc/3Nl2JbPuOOLhpBvhvN9Asm/tAzYuhUdOgvXPa9vBECImyv6xi4qKzOrVq/s7DNVe/RFYegPs+dh3bPSJcOYvIDYh4C3Nbg/bDtey+WANG0qr2VlR1+b8N+aMZsaojpPQRdoH5SsBOCdvXj9H0taRWiePfbobt73A/XlT8/je6RN6/iBnHax6DHZ92Pb4uNPgwv8HI2aGIFrV30RkjTGmKNC52EgHowah3cvhtVugxm+Q04Sz4NTbrS6M7RyuaeKt4oN8urOy7cpbthgRLp45YkAkAYAPKgZmIshJS+C8acP5x0Zreur3tpRzXF4aX/MbhRyUhFQ4/U4omG91La23RyDv+wwenQ+F/wJn/hwyuh73oaKXJgLVe3UV8OE98NXTfgfF+uCYtbBDw3CTy83SNSUs23AYd4CSqACTh6dx5pRc8rSXUFCKxmay70g9Ww7XAvD4p7sZkZ7Yu8n58k+Cyx6BdUtgyxvW9BQYWPccbHgJTvgOnP4TSM8P6d9B9T9NBKrnag7Cyr/CqkXQXOs7njDMWnDev7+6bXdFPQ9+uIPDNW1HtGYkxTExL5WC7BQKspNJjtdfyZ4QEb4+ezQVtXs4Uu/E5Tb84d1t3HPpDEak9yKZxiVbg9AmnQern4RSuxrW3QyrHoc1i2HWVVZpL+/4kP5dVP/R/3UqOFX7rTrkLX+3/jTtGibHnAyn3NpmwXmwppF+f3M5T32xF5fHVwoYk5nM1ybnUpCdjATRpVR1LiE2hoVzx7Do0z00tripaXJx77It/Prr03u/iE/GWDj3bjhUDGufgopt1nFPi1VCWPecNS5k3g9h4tkQwTWkVehpIlBtNR6z5rKv3A1HtkHZZji0rm39v7/0sVB0A+R3bINqcrlZtGIPK3Ye8R6Ld1g9gQrzMzQBhFBmcjxXnzSGZ77ch9tjqKh1cs+bm/mPi6eSmxq4sT4oI2fBiPvh4Fooft5a9KbVzvetn6yJVg+kwms6fBFQ0UF7DQ01Hg/UHrRWrzq2x/rQP7rHWvT82F5oqgruOSNmwfGXwNiTfYuh+NlTWc/DH+6ktKrRe2z4sASuOGFM1C016XRbE98lBNEFtr9tL6/lpdUleOz/11nJ8fz0/ClMyA3RhH3lm2HTq7D/S6DdZ0dsIkz/htWWMPbkoAYPqsjpqteQJoLBqLkBqvbBsX2+D/hje3wf+O4uZq3sTGwi5E61uoSOPRnSAvfvd7o8vLaulDfWHWzTIDw7P4MLZ4wgzqFVCOG29XANL39V6k0GsTHCt04YzcWzRhEfqve/5hBsfdMqEbTUdzyfNRFmL4SZV0BWL7q0qpDTRDAYNdVA5U6o3GV/q99tfdgf2wt1Zd3e3ilHgvUhnzbS6h2SMdb6j5w+psuRwc1uD8u3VfDa2lLv1NEAcQ7hwhkjmJ0fvVUGbx+2lnu8aMT8fo4keLuP1LN0zQGcLl9bTlpiLHPHZzEhJ5W0hFgcMYIBHCIkJzjIS0sgIzkOoQff5FuaYM9y2Pa29TsYyMjZMPXrMPkCGD5DSwr9RBNBNGs8BuVboWKL1WBXsdX6szbAIubBShhmfdCnjYBho6w/U+0P/6TMHv1HLa918uHWMj7cWk5NU9u1dMdkJvP1WSPJ7ksd9QDwi03WxHm/nX57P0fSM5V1Tl5dW8qhdj21upKRFMfM/HTmFmQxe0xG8CUIY6ByB2x/B/augJZOprxIyYPx863OBaPmwPBpEN8/60wMNZoIBjpnLVSXWj1zju62vukf2W799OYDXxyQkmt/sx8BqcN9H/xpIyC+b9Ml1DW7WLXnKCt2VLD5UG2H8ynxsZw5JZfCMRnEDIJvf9GaCADcxsPa/VV8tvNIh0TdndSEWOZNyOLk8dkcP2KYd6K7brmccOBLa6DhwbXg6ep1xRqoljneKnWm5llfRhJSrepIcVg91NzN0NIIzXXW/xdnjf1nrZV0XE77dcQaxBifYq3BkJJj/c6nj7FKt5njrWOD4Peyp/ptZLGIXAD8L+AAnjDG3NfufALwNHAiUAlcbYzZG86YQs7jtr61Nxy1/myqsqptnDXQXG/9kjbXW7/ELQ32L3KddV1DpTU1Q3Ndd6/SUUwspI2yvtG3/rR+s0/JDWqCt544UuekuKSa1XuPUlxa3aYraKthibHMG5/NCWMziY/VtoCBwCExFI3L4oSxmew/2sD+ow0crW/G6fLg8RhErEVuGprdHK130uw3pXWd08UHW8r5YEs58Q5hbHYKI9ISSEmIRURodnmoc7qoaWyhuqmFOqeLxhY3xgOOmCTSEr9OftrFzI3bxTTPNvJqN+Noaf+7bqwvQFXdLLITSvGpkDUeMsZZCSJ9tFVSSc6ylk6NS7amRRH7/5Dx2D9u+0/7PYpxWP8PHfEQl2Qln7iUqOxKG7ZEICIO4BHgPKAEWCUibxhj/Pqf8T3gmDHmOBFZCPw/4OpwxUTDUavHQ2sfeGN8/8Aet9VH2u2yGlNdTt8HubPW+nBvqoLGKusD3NXY1SuFXEVMLsfiR1CTMIrGpFE4k/OIdcTicAgON0g1UO1E2AeyD6DTul7T2tvDtPb7MHg84DHQ4nbjdHmoc7qpaWyhos5JdbuF0sdYD/dKinMwZXgaY7OTEMpxlpXTi+boAcvttKo5akuje0WvbCA7GUgOfN4AR+pa2FdZ36a3FwAecFdAaUXgexPtH//raYCaBnifeN5nJsIMxkkZx8t+jo89xHBzJPDDwq25Dg5vsH7CLS7ZKpkkDrOqZBOHWYkoPsVKHrGJVtJxJIAj1kos4rCSjMT4euSJAAITzoTsiSEPM5wlgrnATmPMbgAReR64DPBPBJcBd9vbS4GHRURMOOqr9n0Oz3zDWowjCuV6KshtqoCmDVAd4RfvOF1QR2X2zyA0LMma8+jkzf/dz5FEUDD/5n0RXTXSvdfSYP3UBrlmRDAu/oM1biOEwpkIRgMH/PZLgPazdnmvMca4RKQa64tLm68KInITcJO9Wyci27p43Zz29wOMTpMRI1JldI/+BmFU0WDITR749ZTREGf4Y7SqMwJWrvZANLyXoHGGWqjjrH7o5mM7j36/ky5aXRrX2YmoGFlsjHkMeCyYa0VkdWcNIgOJiKzeV+XROEMgGmIEjTPUNM7QCWerRil2dbIt3z4W8BoRiQXSsRqNlVJKRUg4E8EqYJKIjBeReGAh0H5R1DeA79jbVwAfhqV9QCmlVKfCVjVk1/nfBryD1X30SWPMJhG5B1htjHkDWAQ8IyI7gaNYyaKvgqpCGgA0ztCJhhhB4ww1jTNEom5AmVJKqdCKvpEPSimlQkoTgVJKDXFRlQhExCEia0XkTXt/kYisF5FiEVkqIgEn0RGRn4vIThHZJiILBmKcIlIgIo0iss7+eTTScfodf1BEOp33or/fz2DiHAjvp4gsFpE9fjEUdnLfd0Rkh/3znUDXDIAY3X7XtO/0EYk4RUTuFZHtIrJFRAJO/BTJ97KPcUb0/exOVIwj8PNjYAswzN6/wxhTAyAifwRuA9rPZzQNqxF6OjAKeF9EJhtj3AMpTtsuY0xhGONqr32ciEgR0Omc0QPk/ew2Tlu/v5/Az4wxSzu7QUSygP/CGq9mgDViTcVybKDEaGvs5/fyeqyu5scbYzwiktf+hn54L3sVpy3S72eXoqZEICL5wMXAE63H/D5cBUgi8MD1y4DnjTFOY8weYCfW9BcDLc6IChSnWPND/R749y5u7ff3M8g4IypQnEFaALxnjDlqf2C9B1wQ6vigTzFGVCdx/hC4xxhrojBjTHmAWyP2XvYxzgEnahIB8ADWf/w2q6aLyP8Bh4HjgYcC3BdoqotwTjXxAL2LE2C8Xcz8WETCvQrKA3SM8zbgDWNMV3NfD4T3M5g4of/fT4B77SrBP4k12257kXw/H6B3MQIkishqEflSRC4PU3ytHqBjnBOBq+0YlonIpAD3DYTfzWDihMi+n92KikQgIpcA5caYNe3PGWO+i1VFsYVwzlwahD7GeQgYa4yZA9wJ/E1EhgW4Lixxisgo4Eo6T1IR18c4+/X9tP0cK/GfBGQB/184Xj8YIYhxnD11y78AD4hI6KfA7DrOBKDJjuFx4MlwvH6wQhBnRN7PYEVFIgBOAy4Vkb3A88DZIvJs60m7fvp54FsB7g1mqot+j9Ouaqm0t9cAu4DJkYoT2AQcB+y0jyeLNdCvvX59P4ONs7/fTxF51hhzyFicwP8RuAotUu9nX2LEGFNq/7kbWA7MCUOMncaJ9e3+FfuaV4FZAe4dCP/Xg4kzku9ncIwxUfUDnAm8iTUj/nH2MQHuB+4PcP10YD1Wph4P7AYcAzDO3Na4gAlYv8BZkYozwPG6Tq7v1/ezB3H2+/sJjPT7d38AuC/A9VnAHqyG70x7O6xx9iLGTCDB3s4BdgDTIvxe3gfc4Hd81UB4L3sZZ7+8n139RFuvIX8CPGUX9wXrw+mHACJyKVBkjPmVsaa1eBFrHQQXcKsJbw+XXsUJnAHcIyItWHWOPzDGHI1gnJ0aYO9npwbg+/mciORi/buvA35gx1lkx3OjMeaoiPw31txcYDU0RjLObmMEpgJ/FREPVi3CfabtAlORcJ8d6x1Y84Lf2D7OAfBeBhUnA+P9bEOnmFBKqSEuWtoIlFJKhYkmAqWUGuI0ESil1BCniUAppYY4TQRKKTXEaSJQA5bfDI2bxJq99d9EJMY+d6aIVNvnt4jIfwW4338G0s0i8rSIxHXzmotF5Ipw/Z16SkRut/9+z3Vy/gERKW19X5TqDf3lUQNZozGm0BgzHTgPuBBrdslWK4w1g2MRcK2InBDgGa0zkM7EGml6VXhDDrlbgPOMMf/a/oT94f8NrPl1vhbpwNTgoYlARQVjzeJ4E3CbPYur/7l6YA3W1BOd3e8G/ok9CZmInGhPRrdGRN4RkZHt7+nsGhH5voisskspL4tIsn38ShHZaB//xD7mEJHf29cXi8jNgeITkTvtezeKyE/sY49ijYpeZg9Qau9MrCk3/gJc4/esXBF5zy5JPSEi+0Qkxz53rYj80y4l/VWsmVzVUNefw5r1R3+6+iHA9BFAFTCctsP6s4G9wPR21xYAG+3tROAjrLlf4oDPgVz73NXAk/b2YuCKbq7J9nuN3wA/src3AKPt7Qz7z5uA/7S3E4DVwPh2cZ5o35sCpGJ9uM+xz+0Fcjp5fx4HrsOaC78UiLOPPwz83N6+AGva8xysEa1/97vuz8C3+/vfWX/6/yeap5hQar6IrMWaQuI+Y8ymANdMFJF1WPMivWWMKRaRGcAM4D27cOHAmq3U35QurpkhIr8BMrA+uN+xj38GLLan4GideOx8YJZfu0M6MAlrHpxWpwOvGqtkg4i8AswH1nb2FxeReOAi4E5jTK2IrMSaj/9N+3nfADDG/ENEWhdmOQcr6ayy/05JQFTMl6/CSxOBihoiMgFwY314TcVqI7ikm9t2GWMK7aqRz+z5iPYAm4wxp3T1cl1csxi43BizXkSuxyqdYIz5gYjMw1qsZI2InGg/50fGmHcCPKcvFmAlog32h3oy0IiVCDojwFPGmJ+HOBYV5bSNQEUFe2K0R4GHjTE9niDLGHMEuAtr/v1tQK6InGI/O05Epre7patr0oBDdg8kbyOuiEw0xqw01qR3FVhTIr8D/LC1t5KITBaRlHavtQK4XESS7XPfsI915RrgRmNMgTGmAKvEc57dXvEZdqO4iJyPb0nPD4ArxF4+UUSyRGRcN6+jhgAtEaiBLMmu1onDmun0GeCPfXjea8DdwDysdoAHRSQd6//BA1h18wAYY5rt6pxA1/wSWIn1Yb8SKzEA/F6sFakE60N3PVCM1Vbxld3IXQFc7h+UMeYrEVmM1ZgN8IQxpqtqoWSsuv8f+D2jXkQ+Bb4O/BpYIiLXAV9grYxXa4w5IiL/Cbxr9zhqAW4F9nX3xqnBTWcfVWqQEWu5SbcxxmWXaP5iBtBC6Wrg0RKBUoPPWOBF+1t/M/D9fo5HDXBaIlBKqSFOG4uVUmqI00SglFJDnCYCpZQa4jQRKKXUEKeJQCmlhrj/H+u2lpjBXygUAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "sns.kdeplot(history_resize, fill=True, linewidth=3,\n",
    "                 label = 'Resize Mean')\n",
    "sns.kdeplot(history_plugin, 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",
    "leg = ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "We have noticed that for the same privacy loss,\n",
    "the distribution of answers from OpenDP'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,\n",
    "but higher utility (lower error) can be generally gained by moving more of the epsilon into the sum,\n",
    "and using less in the count of the dataset rows, as in earlier examples.*"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "3220da548452ac41acb293d0d6efded0f046fab635503eb911c05f743e930f34"
  },
  "kernelspec": {
   "display_name": "Python 3.8.13 ('psi')",
   "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}