performance-curves-python3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Performance curves\n",
    "\n",
    "In this notebook we will show how to compute performance curves (ROC and PR curves to be specific) given only positive and unlabeled data.\n",
    "\n",
    "We will show several results:\n",
    "- the true curves rank, assuming all labels are known\n",
    "- the curves that would be obtained by treating the unlabeled set as negative\n",
    "- bounds on the performance curves given a correct estimate of the fraction of positives in the unlabeled set\n",
    "- bounds on the performance curves given an interval estimate of the fraction of positives in the unlabeled set that is obtained by peeking at the labels of 300 unlabeled instances.\n",
    "\n",
    "This notebook is intended to enable you to easily change the configuration of experiments, so feel free to make modifications where you see fit to see their effect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "import operator as op\n",
    "import optunity.metrics\n",
    "import semisup_metrics as ss\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import pickle\n",
    "import csv\n",
    "import util\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create or load data set\n",
    "Generate simulated data or load the results of an experiment in an existing paper. Feel free to configure the code below as you please."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fraction of positives/negatives that are known\n",
    "# known_neg_frac == 0 implies PU learning\n",
    "known_pos_frac = 0.1\n",
    "known_neg_frac = 0.0\n",
    "\n",
    "# load results from the paper?\n",
    "# http://www.sciencedirect.com/science/article/pii/S0925231215001174\n",
    "use_paper_data = True\n",
    "dataset = 'covtype'\n",
    "#dataset = 'sensit'\n",
    "\n",
    "# generate a simulated data set (requires use_paper_data = False)\n",
    "# number of instances per class\n",
    "# the number of unlabeled instances depends on these settings + known_X_frac\n",
    "\n",
    "# if we use a lot of data, computation times are a little long\n",
    "# but the resulting estimates become excellent, which corroborates consistency\n",
    "# num_pos = 20000\n",
    "# num_neg = 100000 \n",
    "\n",
    "# if you simulate a smaller data set you can see some effects better\n",
    "# but estimated bounds will be much wider\n",
    "# because the ECDF confidence interval becomes large\n",
    "num_pos = 2000  \n",
    "num_neg = 10000 \n",
    "distid = 2  # can be 1, 2 or 3, these correspond to certain curves in ROC space\n",
    "\n",
    "# generate rankings and label vectors and compute corresponding beta\n",
    "# beta is based on known_pos_frac, known_neg_frac and the number of pos and negs\n",
    "# labels is a list of {True, False, None}, where None indicates unlabeled\n",
    "# true_labels is a list of {True, False}\n",
    "if use_paper_data:\n",
    "    labels, true_labels, decision_values, beta = util.load_dataset(dataset, \n",
    "                                                                   known_pos_frac,\n",
    "                                                                   known_neg_frac)\n",
    "else:\n",
    "    labels, true_labels, decision_values, beta = util.simulate_data(distid, num_pos, num_neg,\n",
    "                                                                    known_pos_frac,\n",
    "                                                                    known_neg_frac)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data set characteristics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Known parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total number of instances: \t 40000\n",
      "number of known positives: \t 2000\n",
      "number of known negatives: \t 0\n",
      "number of unlabeled instances: \t 38000\n"
     ]
    }
   ],
   "source": [
    "print('total number of instances: \\t %d' % len(labels))\n",
    "print('number of known positives: \\t %d' % len(list(filter(lambda x: x == True, labels))))\n",
    "print('number of known negatives: \\t %d' % len(list(filter(lambda x: x == False, labels))))\n",
    "print('number of unlabeled instances: \\t %d' % len(list(filter(lambda x: x == None, labels))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unknown parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of latent positives: \t 18000\n",
      "number of latent negatives: \t 20000\n",
      "beta: \t 0.4737\n"
     ]
    }
   ],
   "source": [
    "print('number of latent positives: \\t %d' % len(list(filter(lambda x: x[0] == None and x[1] == True, zip(labels, true_labels)))))\n",
    "print('number of latent negatives: \\t %d' % len(list(filter(lambda x: x[0] == None and x[1] == False, zip(labels, true_labels)))))\n",
    "print('beta: \\t %1.4f' % beta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Estimate beta\n",
    "\n",
    "This step can be avoided if sufficient domain knowledge exist to guesstimate beta, or when using data-driven approaches as mentioned in the paper. \n",
    "\n",
    "In this notebook we use a point estimate you can provide (\\hat{beta} in the paper). Change this parameter to see its effect in the figures later on. Additionally, we will perform Bayesian inference on beta by peeking at the true label of some unlabeled instances, as described in Section 6.3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "true value of beta\t0.4737\n",
      "point estimate of beta\t0.4737\n"
     ]
    }
   ],
   "source": [
    "betahat = beta\n",
    "\n",
    "print('true value of beta\\t%1.4f' % beta)\n",
    "print('point estimate of beta\\t%1.4f' % betahat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, estimate beta by peeking at the latent labels of a few unlabeled instances. This corresponds to Section 6.3 in the manuscript."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Credible interval via PPF: 0.3947 <= beta <= 0.5066.\n",
      "Width of credible interval based on PPF: 0.11196.\n",
      "\n",
      "Credible interval via HPD: 0.3947 <= beta <= 0.5066.\n",
      "Width of credible interval based on HPD: 0.11196.\n"
     ]
    }
   ],
   "source": [
    "# determine the latent labels (requires manual labeling in practice)\n",
    "latent_labels = [true_label for label, true_label in zip(labels, true_labels) if label is None]\n",
    "\n",
    "# we will peek at 300 labels to estimate beta based on a 95% credible interval\n",
    "num_peeks = 300\n",
    "credible_interval_length = 0.95\n",
    "\n",
    "# peek at the true label of a subset of unlabeled instances\n",
    "unlabeled_ranks = [idx for idx, lab in enumerate(labels) if lab == None]\n",
    "sample = np.random.choice(unlabeled_ranks, num_peeks, replace=False)\n",
    "sample_labels = [true_labels[i] for i in sample]\n",
    "npos = sum(sample_labels)\n",
    "nneg = sample.size - npos\n",
    "\n",
    "# compute interval based on the percentile function of the posterior\n",
    "beta_interval_ppf = ss.estimate_beta_by_labeling(npos, nneg, credible_interval_length, HPD=False)\n",
    "print('Credible interval via PPF: %1.4f <= beta <= %1.4f.' % (beta_interval_ppf.lower, beta_interval_ppf.upper))\n",
    "print(\"Width of credible interval based on PPF: %1.5f.\\n\" % (beta_interval_ppf.upper - beta_interval_ppf.lower))\n",
    "\n",
    "# compute interval based on the highest posterior density region (shorter)\n",
    "beta_interval_hpd = ss.estimate_beta_by_labeling(npos, nneg, credible_interval_length, HPD=True)\n",
    "print('Credible interval via HPD: %1.4f <= beta <= %1.4f.' % (beta_interval_hpd.lower, beta_interval_hpd.upper))\n",
    "print(\"Width of credible interval based on HPD: %1.5f.\" % (beta_interval_hpd.upper - beta_interval_hpd.lower))\n",
    "\n",
    "# continue with beta based on HPD, because this is the shortest credible interval\n",
    "# note that the coverage of both intervals is equal to credible_interval_length\n",
    "beta_interval = beta_interval_hpd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Compute cumulative rank distribution of known positives"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sort the labels in descending order of corresponding decision values\n",
    "sort_labels, sort_dv, sort_true_labels = zip(*sorted(zip(labels, decision_values, true_labels),\n",
    "                                                     key=op.itemgetter(1), reverse=True))\n",
    "\n",
    "# ranks of the known positives\n",
    "known_pos_ranks = [idx for idx, lab in enumerate(sort_labels) if lab]\n",
    "\n",
    "# compute rank ECDF of known positives\n",
    "known_pos_ecdf = ss.compute_ecdf_curve(known_pos_ranks)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we determine a confidence interval on the rank CDF of known positives. We can do this in several ways, our code provides a bootstrap approach and a method based on the Dvoretzky–Kiefer–Wolfowitz (DKW) inequality. Feel free to experiment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "ci_width = 0.95      # width of the confidence band on ECDF to be used\n",
    "use_bootstrap = True # use bootstrap to compute confidence band\n",
    "nboot = 2000         # number of bootstrap iterations to use, not used if use_bootstrap = False\n",
    "\n",
    "if use_bootstrap:\n",
    "    cdf_bounds = ss.bootstrap_ecdf_bounds(labels, decision_values, nboot=nboot, ci_width=ci_width)\n",
    "else:\n",
    "    cdf_bounds = ss.dkw_bounds(labels, decision_values, ci_width=ci_width)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For reference, we will also compute the rank ECDF of latent positives. In practical applications this is impossible, as the latent positives are by definition not known."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "latent_positives = map(lambda x, y: x == True and y == None, true_labels, labels)\n",
    "sort_lps, _ = zip(*sorted(zip(latent_positives, decision_values),\n",
    "                          key=op.itemgetter(1), reverse=True))\n",
    "latent_pos_ranks = [idx for idx, lab in enumerate(sort_lps) if lab]\n",
    "latent_pos_ecdf = ss.compute_ecdf_curve(latent_pos_ranks)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the rank CDFs of known and latent positives. It may occur that the rank CDF of latent positives is not within the confidence interval of known positives, in which case the corresponding bounds on performance will not be strict."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzs3XlclOX+//HXxQAii6iAouKCiguCIOKuiPuWS+5bZaVmadlyTKtz2s75nezU6Xw9ZZpZmmWlaZnlvu8LoIgL7ooioiiK7AzD9fvj1jnuqTkOy+f5ePDImfueez630rznvu/r/lxKa40QQggB4GDvAoQQQhQeEgpCCCGsJBSEEEJYSSgIIYSwklAQQghhJaEghBDCSkJBiAeglHpXKfWdvesQ4mGTUBDFllLqpFIqWymVoZRKVkrNVkq526mWpkqppUqpy0qpVKXUTqXU01eXRSqlCq7WmaGUSlRKzVdKNblpG1oplXndepftsS+ieJNQEMVdT621OxAKNALeeNQFKKVaAGuBDUBtwAt4Huh23WpJV+v0AJoDB4FNSqkON20uRGvtfvWnrO2rFyWNhIIoEbTWycAKjHAAQCnVQym1Wyl1RSl1Win17nXLalz9Zv6UUuqUUuqCUuqt221bKeWklPpBKbVQKeV8m1U+Ar7RWn+otb6gDTFa64G3qVNrrRO11m8DM4EP/+SuC3FfJBREiaCU8sP4Zn70uqczgSeBskAP4HmlVJ+bXtoaqAt0AN5WStW/abulgUVALjBQa51303JXoAWw4AHK/hkIU0q5PcBrhXggEgqiuFuklEoHTgPngXeuLdBar9da79VaF2it44AfgLY3vf49rXW21noPsAcIuW5ZGWA5cAx4Wmttuc37l8P4/+zsA9SeBCiM0Lpm19XrEpeVUv99gG0KcVcSCqK466O19gAigXqA97UFSqlmSql1SqkUpVQaMOb65VclX/fnLOD6C9XNgYbAZH3nzpKXgAKg0gPUXgXQwPUXlMO01mWv/rz0ANsU4q4kFESJoLXeAMwGPr7u6e+BxUBVrbUnMB3jm/m9Wgl8AKxRSlW8w/tmAduAfg9Q9uPALq115gO8VogHIqEgSpL/Azoppa5dbPYAUrXWOUqppsDQ+92g1vpfGOGyRil181HGNa8DI5RSE5RSXgBKqRCl1I83r6gMVZRS7wAjgTfvtyYh/gwJBVFiaK1TgDnA364+9QLw/tVrDm8D8x9wu3/HuNi8WilV/jbLtwLtr/4cV0qlAjOApdetVlkplQFkAFFAMBCptV75IDUJ8aCUTLIjhBDiGjlSEEIIYSWhIIQQwkpCQQghhJWEghBCCCtHexdwv7y9vXWNGjXsXYYQQhQpMTExF7TWPn+0XpELhRo1ahAdHW3vMoQQokhRSiXcy3py+kgIIYSVhIIQQggrCQUhhBBWEgpCCCGsJBSEEEJY2SwUlFJfK6XOK6X23WG5Ukr9Vyl1VCkVp5QKs1UtQggh7o0tjxRmA13vsrwbEHD1ZzQwzYa1CCFEkZaTk/NI3sdmoaC13gik3mWV3sCcqxOVbwfKKqUeZHYqIYQotg4cOMDkyZMJCgpi/vwH6u5+X+x581oVjHlzr0m8+twtc9kqpUZjHE1QrVq1R1KcEELYg9aa999/n9Wr1xEXF8uVK2nWZRMnfsCAAQNQ6n4mCLw/9gyF2+3VbSd30FrPwJiUhPDwcJkAQghR7OTnw8aNF5gyZRq///4PCgryrcscHV2pXLk5n3zyi00DAewbColA1ese+wFJdqpFCCEeOa01Gzfu4W9/+ytmc2V2715Abu4lAJycPKhZszNt2vyF/v2bExYGPn/YuejPs2coLAbGXZ2nthmQprW+5dSREEIUF8nJmcybt5Ljx4+zefMWDh/eSkbGuRvWqVSpBT17/o1BgzrTqpUJZ2ew8cHBDWwWCkqpH4BIwFsplQi8AzgBaK2nY8xP2x04CmQBT9uqFiGEeNTy8yE9HdLSYNMmWLt2HYsXP09q6qHbrl++fCjDh3/K00+3IjT0EabATWwWClrrIX+wXANjbfX+QghhS6mpsGsXHDwIx49v5siROHJzL5GTc5ns7EucOhVFZuZFypSpSm5uBqmpN96y5eTkiqenP2XLNqVJkyE0adKRUaMU7u522qGrilzrbCGEsJesLAuTJ0eza9dWAgJGkZPjTrly+axZM4W4uAW3fU1m5hkATCYngoJGULXqU/j7N8DT05MaNRTVq0PLluDq+ij35M4kFIQQ4iYpKXDlCuzdm8vJkxfIzr7M4cOHWbXqX5w5sx2A3r1zSUjYwYEDK8nLy7rtdkwmJ/z8WlCnzku0adODhg1dqFIFqlQxLho7FsJP4EJYkhBCPHr5+ZCcDPv2QUwMnD4NCxaEc/Hijad9nJ1dcXZ25ddf37A+5+7uh6trRcqVq0bZsjXx9KyGu3sd/P1b4+vrTtWqMGAAOBSBbnMSCkKIEqOgwLgGkJEBR49CdjZkZZlJSrrIkSNbOHJkDbVqDaZ69Ua4usZQtmxZ0tPdMJkcyc/Pw2zOJi8vy3pk0KTJaEaOfIvmzavh5ASVburJ4OYGTk522NE/QUJBCFEixMXBmjVw7Bj88ktL0tJOkJeXhtmcfcN6qanb+O23A+Tn592yDZPJhVq1OtG162AGDnyMli3LPNLhoo+ChIIQotg6dw6iomD16r04OgaSkWEiOBg2bvQiKWnb1bUUTk4ugMZsziExMRaAChUaUrFiK3x8fPHx8cHXtyGNGjWldWsnatWy2y7ZnISCEKJY2LnTOBpIT4dLlzSZmVc4enQNu3ZNJzFxFQMG/Bdv72r88ssGzp8/fN0rtfVooXRpb/z8utG27YsMG9YEf3+oXt0++2MvEgpCiCJJa9ixw7g4fOgQnD0LO3d+wIkT33Px4lHM5v+1mjaZnFi9+iMuXfpfD04HB0dcXMrh6loBH59QGjToTv/+g6lSxYGgIChb1h57ZX8SCkKIQi8+HhITYc8e4yKxszNkZuazZ88SKlRoSI0a/tSrB7t3J5CcbIwWcnQshbOzKzk5V7BYzFy6dBqTqTStWo2lS5fOtGoVQf36pXBwgFKljPsETCY772ghIKEghCi0zGaYMQNiY1NIStrIxYvbgYtkZJxj//6lAPTu/Q/8/Qcyc+ZfOXUqBkdHZ/Lz88jPzyU/PxelHKhatRP16vWlU6cBjBzpRbly9t2vwkxCQQhRqOTmwt69kJlpXCc4ckTz++/tSE7ef8u6jo4uHD68geXLPyA3N9P6fJky1fHxaYiPT3OaNh1C8+b+tGkDfn6Pck+KJgkFIUShkZ4Os2bB7t1XcHJKxtu7DmFhCrP5dbZt+4YqVYIBuHLlLPHxq8nMTCU+fhUAtWoNoVOn16lTpybBwWWoXt3oLlqtmnG6SdwbCQUhhF1ZLHDxonGxeMmS88ye3Y9z5zbj4eHDhx8msXbtpxw+vJ6UlKMcOrT2htd6eTUiMPBx6tePoHv3tvToUThbRxQl8tcnhHjkLlyA/fuNC8hJScaponPn9vDzzx3Izr6IUg6YTM588EFTTp/ebX2dk5M7fn5N8fGpS926j9GxY1d693bA09OOO1PMSCgIIR6p3383bihLTTVG/Zw+/RV79szl0KF1gNFbyGIxc/nyGS5fPkOpUl707Pk+ffq0p1OnACpUkCFCtiShIISwKYvFmGRm3z7jNNGuXasoW9aTKlUUJ07sZNu2r244GsjLy0IpBypXbkGNGoPp1m04L79c3u7zDJQUEgpCiIdi2zbjhrK0tP9NH3n6tNGG+tSpfHJyDpGQMIv16/+NyeSExWK+ZRvu7t5Urz6CXr1eJjKyCp6e0KRJ0eguWlxIKAghHpjZbLSW2LwZDh8GDw/jA9zFBfLz87hw4TBxcdOJi/uJ9PTz1tdZLGYqVqyLv39T/PxC8PKqgatrTS5fbkidOiaeffbRTFIvbiWhIIS4b3l5sHy5cYfx2bOaAwe+4PTpb/H3b0bfvpNJTz9PcnIM06b1ueW1pUt78sQTX9Ko0QDOnjXuUC4oMFpMN2tmzDtQvrwddkoAEgpCiD+Qn2+0mzaZICfHOCo4fx5On87mxImZHD8+lxMndgBw5kwc69b9l4ICyy3bqVu3HZGRYwkI6MqhQ24cPw7lykFQEAQGQsWK4O39qPdO3ExCQQhxW4cOwdatRs+h5GSsrSHS0sycOvU5GzZ8SFra2Rtek5ubgVIOeHpWokyZipQvX52aNVsQFtYPN7fapKbCkSNGEAwdCl5eFLv5CIo6CQUhxC2+/tqYmSw5uQBYSnb2SUJCxgGQmZnBZ5+9DIBSCq01YBwJ9Or1d6pWDaVUKTfrttLSjCDw9YXKlY1TRN27y8XjwkpCQQgBGOf1jxyBpUshLu4MO3e+RELCBjIzLwIQFtaP06djOXp0i3X0kNaaatUa06LFc3h6PovJ5EBCApQubWwzN9doMdG8OfTqJaeHigIJBSFKqGsT1WdmGvcPbN9uDCE9dOh3li3raV3PwcERrS28/nrlG17v79+Mnj0/Jz8/DBcX8PeHvn2NEHBxufZaIyDkFFHRIaEgRAljscCWLcbP+fPGqB+tjeePHfs7y5a9DYBSDmhdQEFBPgCurmWpUKEOAQER1KzZgYKCzri4OBAQAJGRULu2HXdKPDQSCkKUEFrD6tVGi4nkZFCqgMqVk8jNzcTXty4AR4/q69YvwM8vhP79/0OZMs1QypW8PKNvUblyULUqtGoFISH22iNhCxIKQpQABQXwxRdw4AAolUp09GtER/+I2ZyDv38zxo79jQULXmP79m8BKFu2Mh07vk3Vqs9gNjvh6mocCZQubZweatoUypSx804Jm5BQEKKYu3gRvv0W4uPNJCV9xurV75KTcwWAUqXcuXQpkb//PYS0tLM4ODjSs+e7VK36KqVKlaZWLWjXDurVs/NOiEdGQkGIYigpCaKjjTuOc3Ph8OFTbN7cnbNnjdnL3N19yMhIITc3g9zcDACqVw/nqafmkpJSh7JloV8/49SQXCQuWSQUhChmtm2DJUuMWcy8vIwpKCtX9iU62gMPjwq4u3tz9uwBTCYnwsL6UatWK6pWDScvrymXLztQpw507SpHByWVhIIQxci5c0ab6jNnzpOZOZf0dB98fYeTmpqEi4sH6ennSU8/j7u7Nz17LqFcuaa4uBizlfn4QIcOEBZm770Q9mTTUFBKdQWmACZgptZ68k3LqwHfAGWvrjNJa73UljUJURzl58P8+caF5OXLPyQmZhIAVao05MqVZH799a/k5+fi6OhMixYTqF79ZYKDvWnTxjgicHc3hqYKYbNQUEqZgKlAJyARiFJKLdZaH7hutb8C87XW05RSgcBSoIatahKiuElPhx07YPduo3V1WtpcayCUKeNLWloSCxdOACA0dBC1an1ASIg/jRpBx47/u8lMiGtseaTQFDiqtT4OoJT6EegNXB8KGrg2sM0TSLJhPUIUC2YzHDxozGR29CgkJED58pr8/LksX/4X63pXriQDRjg8/vgnlCo1hJAQGDxYWlOLO7NlKFQBTl/3OBFodtM67wIrlVIvAm5Ax9ttSCk1GhgNUK1atYdeqBBFhdkM06YZYWA2GxeRmzaFhQsnsnLlR9b1HBwcadduHMHBj+Pq2oorV0wEBMCYMdKITtydLX89bjeQTd/0eAgwW2vtB3QHvlVK3VKT1nqG1jpcax3uI9MxiRIqP9+4AW3vXggIKKBcuWh8fCAnJx1X13I4OBgXBdzdfRg5chkBAf/BxSUCNzcT7drBSy9JIIg/ZssjhUSg6nWP/bj19NCzQFcArfU2pZQL4A2cRwgBGNcMjh2D+HhISkrjwIG/8P33C8jOTmP8+FV8//0Yzp8/ChidTNu2nYqTU0Xq1zdGE1Wt+gdvIMR1bBkKUUCAUsofOAMMBobetM4poAMwWylVH3ABUmxYkxBFyoYNxj0HFgt4e5vZurUfBw+uAcBkcubTT7thsZipVCmQxx+fTHBwT2JjoX17GDjQzsWLIslmoaC1zldKjQNWYAw3/VprvV8p9T4QrbVeDLwGfKmUegXj1NIIfW3GDiFKMK1hzRrjRrSMDPDx2cv8+e9w8OAaTCZHLJZ8LJY8AJo0GcLw4TNITnbnwAGjR9Fjj9l5B0SRpYraZ3B4eLiOjo62dxlC2ERuLvz6K5w4YcxtUKYM+PpmMH68xw3rOTqWokGDLjRv/iQBAX3Zt09RuzaEh0OXLnYqXhRqSqkYrXX4H60ndzQLUUjk5cGXX0JcHFSpYsxj7OAAeXkOuLt7k5l5Ea01tWu35vnnF3H5shcpKUbDu5YtjQluvLzsvReiqJNQEKIQ2LwZdu0yRhYFBhawadO/8PCoQJUqwXzxRX8yMi4AEBHxHCEhn3HypCO+vtCtG7RpY8xvIMTDIKEghB2dPm3MiRwfD6VKwcWLU3jvvQ9JSzsLGKeJ8vNz8fMLYejQz1GqJU5Oxt3IzZsbrxHiYZJQEMIOsrONXkWHDsGlS1CzJqxaNYFVqz4GjJvPCgryyc/PpWXLpxk8+AtOnXJCa+OaQdu2dt4BUWxJKAjxiF2+DN99B7GxUKMG5OVt4Isv3uTYsa3WdQoK8nFz86Jnzw8pU+ZZzpwxrjO0aAHNbu4LIMRDJKEgxCOitREG8fHGHMmNGhmnf44fP8Xx49ut69Wu3Zru3f9KvXrt2bHDiaAgGDVKmteJR0NCQYhHID0dvvnGaGJXsyZYLFspVaol2dlpHD68Aa0LAGjV6hmGDJlKVpYLMTEQEAC9e0sgiEdHQkEIG1u82Dg6OHgQata8wOLFE9my5WsGD/6MVas+4uLFBEwmJ/r3/5jGjccRF+eAr69xmmjIEGOuAyEeFQkFIWzkwAH4/XdITAQPD9D6O95++wnr8t9+e4fMzItUrRrK009/i6dnEMePQ8OGxukiZ2c7Fi9KLOmZKMRDpjXs3Am//Qb790NgoCY+fjLffGMEgpubFyaTE5mZF6lWLYxXXllDcnIQp05B/fowdqwEgrAfOVIQ4iHat88Ig7NnQSnjTuPZs0ewffscAJRyIDPzIgCtWo2kRYupxMc7U6eOcUdy9er2rF4ICQUhHoo1ayA62hhump5unAK6pmLFOiil0FqjdQH16nWgffu3ycuLwMvLuBGtVSu5diAKBwkFIf4Es9mY+ObYMeNxQACkpMSzfv06IiNfYO/eJaxb9ylaa1xdy/HUU19Ts2YfDh+G0FAYOVLuShaFi4SCEH/C3LlGv6IGDSA7O5EPP+zDqVMxAOzfv5y4uN8A8PdvzrPPziUhoSaXLxtTaI4YAY7yf6AoZORXUogHNHu20cTOzS2BH398l23bZluXVagQQFzcbzg6OtOr1z+IiHiZ2FgnAgNh8GCZDU0UXhIKQjyAmBijb1F8/Cds2vQOubkZANSs2ZL09HOcP38ER8dSjB+/glKl2rJnDzRubNx3IO2tRWEmoSDEfTp40JgmMzlZ4+qahdmcg69vPXx96xMX9xsFBflUqlSfPn2+Ij29BWXLGqeKGjSQQBCFn4SCEPdhzx7jDuWTJ6FZM4WDw1/x9vZnwYLXiI39BaUUERFjaNz4X5jNHjRsaASCtKkQRYWEghD3aPly4wjh8uU8qlQ5iMVSlx9+eJmNG6cDUKNGU3r2/A4HhwAqVDBGF3XtaueihbhPEgpC/IH8fFi1CjZtgpycy2zb9hRxcb9TunRZsrJScXJyISLiH9SsOZ7y5R0JD4fHHrN31UI8GAkFIe4iMdE4XbRvH5Qpk87q1cPYt28pAFlZqfj61icsbDpNmkTQpAl07gxOTnYuWog/QUJBiNvIyYH16yEqCk6cAF/fY0yZ0oSsrEvWdXr1eh9v7zeoV8+R4cOhUiX71SvEwyKhIMRNUlPhq6/g+HFjtJC39xY++qgD+fm5ALi4lOGpp2YBffH1NWZDk0AQxYWEghBXZWTA2rXGDWmnTkHz5rB37xKmT+9nDYRmzYYzYMAn5OT4kJlpXDsID7dz4UI8RBIKQmDMfbBkCRw+DD4+UKnSLubMmcqWLV8DULasH8899xM1azYnJcU4pdS4sQSCKH4kFESJt327MRlOerrxQT9tWi/27l2K1hZrm4pOnV7FwcHEuXOQkgJNmsATT/zxtoUoaiQURIm2caNx/0FWFlSseJJ33unM+fNHAPDwqMC4cUuoUcM4HDh1CiwWaN0aBgww5ksQoriRUBAlUlaWcbro4EG4fLmAxMR/MGXKO9blwcE9GDv2t6vzIBinlQB69ID27e1UtBCPgISCKHFycozRRXv2QKVK2axe3Y0jRzYAxsxoAwf+h/btXwKMG9diYox5EiIijJnUhCjOJBREiXLgAKxbZ8yBEBoK+/Yt48iRjYBxMXnMmAX4+zdDa7hwwbig3KAB9O4NtWrZuXghHgGbhoJSqiswBTABM7XWk2+zzkDgXUADe7TWQ21ZkyiZzpwxhppu3w5paVC+/CEWL57J6tX/BjSNGw/giSe+pHRpT7Q21q1a1bh+MHQoODvbew+EeDRsFgpKKRMwFegEJAJRSqnFWusD160TALwBtNJaX1JKVbBVPaLkOnwY5s2DhATjgz4x8Z/8+utfMb6HQI8ef6Nnz/dQSnH6NFy8CHXqwPjxMjOaKHls+SvfFDiqtT4OoJT6EegNHLhunVHAVK31JQCt9Xkb1iNKoLg4WLMGjh6FkJALzJgxiEOH1gJQoUJtnnlmLv7+TcnJgdhY8Pc3rh307y+BIEomW/7aVwFOX/c4EWh20zp1AJRSWzBOMb2rtV5+84aUUqOB0QDVqlWzSbGi+Dl4EBYsgHPnLKSm/j9ef/0fWCxmAMqXr8bf/haHs3NpUlLg3Dlo1AieeQbc3OxcuBB2ZMtQuN0obn2b9w8AIgE/YJNSKkhrffmGF2k9A5gBEB4efvM2hLhBVpZx7WDrVjh3TrNmTSTHjm0GjNFFXbu+Qa9e7+Pg4MDJk8YIo1atYNAg+9YtRGFgy1BIBK6fntwPSLrNOtu11mbghFLqEEZIRNmwLlFMWSwQHW10Nz1xAtzccjl6dKI1EGrWbMHo0fMpV86PggJjJFJBgXExeeBA+9YuRGFhy1CIAgKUUv7AGWAwcPPIokXAEGC2Usob43TScRvWJIqphATjVNHJk8bjcuV28s03Izh7Nv7qFJnPM2TIZyilKCiAHTsgONhoeteihV1LF6JQsVkoaK3zlVLjgBUY1wu+1lrvV0q9D0RrrRdfXdZZKXUAsAATtNYXbVWTKJ4WLDDuOzh5EurVO8vMmf05fnwrABUr1uXJJ2dSu3ZrwGiLHR9v3HswfDhUkPFuQtxAaV20TtGHh4fr6Ohoe5chCoGsLPj5Z+OO47JlNdu2vcmqVR9TUJAPQEBAW8aPX46TkwtXrhizqJUqBSEh0KyZ3IwmShalVIzW+g/7+sqgO1HkaG3cd7B/P1y6BGXLZjJ1akNSUowzj0opOnZ8la5dJ+Hk5EJamtHMrl4944JyaKidd0CIQkxCQRQ5ixbBtm3GrGiZmd8yffqz1qGmfn6hjBr1I76+dQFITjbuZg4Nheeft2fVQhQNEgqiSFm50pg32cUFUlNXMnfuc1gsZpRyoE+ff9K160TrumlpxpFE167QpYsdixaiCJFQEEXG9u2waRPs2fMt+/f/l4QE49pSWFh/hg//Aje38oAxzPTUKaNdRWgoPP64PasWomiRUBCFXm4uLFxotLres+cbVqwYAUCpUu506/YmnTtPwGQyfpXT0oyWFrVqGe0qZO4DIe6PhIIotNLTYc4cOHvW+NYfH/8eGza8C4C7uzdvvBGFt3cNAPLyjFFIfn5Gu4rRo8HBwX61C1FU3XcoXO1+OlhrPdcG9QgBwL59xrzJR45A7dqXWbZsFLt2LQDAxaUMb74Zg5eX0Qfr6FG4csUIg4EDwdfXnpULUbTdMRSUUmWAsRiN7RYDq4BxwF+AWEBCQTx027cbdxsfO2Z80/fx2c5///s0yckHAXBz8+Ldd/dTpkxFCgqMQHB3h8ceM+5OlnkPhPhz7nak8C1wCdgGjAQmAM5Ab6117COoTZQg33wDp0/D+fPg5AS+vsn89ttEtm+fY12nfPlqvPLKGsqUqUhSkjEzWsWKxnWDVq3sWLwQxcjdQqGm1joYQCk1E7gAVNNapz+SykSJkJkJ331nzHtQsSIEBsLBg6v44INhZGSk4OhYiqCg7mRkXOCFF37Fza0cCQlGZ9PWrWHwYHvvgRDFy91CwXztD1pri1LqhASCeJgSE2HuXOMUkL8/aH2C//73OeLjVwFQq1Yrhg//gkqVAlHK6MSelGRcgG7SRAJBCFu4WyiEKKWu8L95EUpf91hrrcvYvDpR7OTmGhePT5wwmtgdOQIBASeYN288cXG/WddzcyvP2LG/4eZWzvrcpUtGIDRtajSzE0I8fHcMBa216VEWIoq3ggKYP98IgeRkMJmgXDmoWfMY770XTF5e9tU1FZ06vUa/fv+yHh2AceHZYoGwMBg2zD77IERJcLfRRy7AGKA2EIfR+jr/URUmio/MTPjiCzh8GCpXNrqUAuTmZvLPf/a0BoKfXwhPPz0HP7+GgBEkB41BR7i6QsuW0Lu3PfZAiJLjbqePvsG4rrAJ6A40AMY/iqJE8ZGRYQTC/v1Gu2qlICvrMjk5GUyb1ofk5HgcHEwMGvRf2rZ93np0kJdnnF6qU8d4Xb164O1t550RogS4WygEXjf66Ctg56MpSRQXmZnGUNP4eGjSpIBduxayatUnnDwZhbu7N+np5/D2rsmLLy7B17ceANnZxhGFq6sxEc7YscapJiHEo3Gvo4/yrz+/K8Qfyc+HWbNg1y4oW3Y7H374orWBHUB6+jlq1mzB2LG/4e7uRWoqpKQYp4wCAozhpuF/OB2IEOJhu1sohF4dbQTGiCMZfSTuSWIi/PCDMdTUyyuGqVO7kp2dhsnkjMWSB0BoaB+efPJr3NzKcfas0aYiMBCCgozTRUII+7hbKOzRWjd6ZJWIYmH9eli3zgiGoKBU3nmnB9nZaTg4mLBY8ihTxpfhw2cQEtITMK4bODoafYuefNK+tQsh7h4KRWvyZmFXeXnw9dd1k0suAAAgAElEQVTG/QcFBca3/dxcF8xmY2RRQYGFsLD+DBs2DXd344rxnj1QqRL07AkNG9qzeiHENXcLhQpKqVfvtFBr/YkN6hFFUGqqcWfyvn1QrRp4e2uiouazcOEEcnKu4OpajoED/0PTpk9y+rTi1CljFJKfnzEr2rUhqkII+7tbKJgAd/53R7MQt0hKgpkzISEBQkI0Gzf+m/Xrp3Lx4kkAKlduwMiRP2AyBVtDo0MH478BAfatXQhxq7uFwlmt9fuPrBJR5GzZAjt3GqeMgoMvMGVKL44f3waAUg4MGjSFtm1fIDnZgStXjG6m/frZuWghxF3dLRTkCEHc0Y8/GoHg7Ay+vvt45512ZGRcAMDZ2ZUxY36mQYMuJCcb9x60aCGBIERRcLdQ6PDIqhBFyi+/GIFQrhzs2fMxCxdOsC6rUKE2Y8f+jq9vXU6cALMZ2rY1LiYLIQq/uzXES32UhYjC78oVmD3bOF3k5gYm0wl+/nmidXlQUA9GjJiFh4cPSUlGR9QOHaBHD/vVLIS4P/c9R7MomSwW+Oor476CsDDIzU3hyy9HonUBDg6OPPnkTFq0eAqzGaKijMZ37dpJIAhR1EgoiD905gwsXGh0LG3YMIOlSz8gOnoeKSnH8PDwYezY3/H3b8q5c8ZopOBg6N/fCAYhRNEioSDu6qefYPdu48O+fPktvPvuQNLSkgCoVi2MF15YTLlyVTh1CnJyjKOIESOMeZaFEEWPhIK4LYsFfv0VoqON6wfe3hv45JMOaG0BIDT0cZ59di7OzqWJjzdCoHlz4whBCFF0SSiIW5w5AwsWGC2vfX3h/PnfmDbtcWsgdOv2Fr17/x2tFXFxRquK3r2NZnZCiKLNwZYbV0p1VUodUkodVUpNust6/ZVSWiklzZLt7NIl4x6EqCioVy+P+fM7M3VqLwoKLIBi6NBp9OnzD/LyFLGxRiC0bi2BIERxYbMjBaWUCZgKdAISgSil1GKt9YGb1vMAXgJ22KoWcW+2bTM6nB45YnzQp6YmceHCCQAcHUvxwguLadCgM3l5xlFEgwYwbpzRx0gIUTzY8vRRU+Co1vo4gFLqR6A3cOCm9f4O/Av4iw1rEX9g2zZYvhwuXoTatY+TnJzNlCmduXw5idKlyzJ+/HL8/ZuRlATnz0P9+hIIQhRHtgyFKsDp6x4nAjdMn6KUagRU1Vr/rpS6YygopUYDowGqVatmg1JLroICmDEDDh0yLhafO/cJM2e+g9YF5OVlERAQwejRP1GmTAUOHzbmPujYER5/XAJBiOLIlqFwu48M6xwNSikH4D/AiD/akNZ6BjADIDw8XOZ5eIg2bjRuSKta1cKyZePYuHG6dVlYWH+efnoOzs6lOXXKCI127aBzZzsWLISwKVteaE4Eql732A9Iuu6xBxAErFdKnQSaA4vlYvOjs3On0enUYsng++973xAI4eGDGDXqR7Quzc6dxlFBhw4SCEIUd7Y8UogCApRS/sAZYDAw9NpCrXUa4H3tsVJqPfAXrXU0wubWrYMNG+D48XhWrepHcnK8ddm1IwSlTOzbZ0yC07+/MdJICFG82SwUtNb5SqlxwAqMCXu+1lrvV0q9D0RrrRfb6r3FnWkNn38OR4+C1uksWtSW9PQU6/KePd+jR4+/kZCgSE+HunXh+eeNawlCiOLPpv+ra62XAktveu7tO6wbactaBKSnw5w5cOCAMetZ6dLuNGv2BKtXGzOrDh06jbZtxxAVBVWrQtOm0KePnYsWQjxS8v2vhEhOhu++MwKhXr0cLlyIZ9OmGWzcOB2lHHjiiS9p2vQZtm2D2rVh+HCoXt3eVQshHjUJhRLgwAGjj9GhQ1C27Eo+/HA0V64kk5+fi8nkxOjRP9GgQW/i4iA0FAYPhooV7V21EMIeJBSKubg4o+31mTOg9WymTHmGayODXV3L8dJLy9G6KQcPGncojx0LJpN9axZC2I+EQjG2ebMxyujCBcjJ+YK5c8dYl/n41GL8+BXk5dWioAC6dTN+5IY0IUo2CYViymw2AuH8eTNRUU8RFfWDdZmXVw3+8peNKFWZo0eNltfdu9uxWCFEoWHTLqnCPgoKYOZMSEwET8+9HDq01rqsbNkqvPbaepSqzMmTxggjmQNBCHGNHCkUM1rDtGlG64patSA3151rHUeqVQvj+ed/4cIFo39UeDg8/bQdixVCFDoSCsWI1jB9ujHaqHr1dDIyTvCf/3QgI+MCdepEMnbsEvbvd8XXF3r0MI4ShBDiehIKxcj33xujjUqV2smbbzZDKYXWmgYNujJq1AJ27XIlMNDoXxQaau9qhRCFkYRCMbF0KezfDzk5W5k5sz0AWmvq1GnLE0/8woEDLjRsCE8+Cd7ef7AxIUSJJaFQDMTFGcNPExI28dNPXcjPzwXA378Zw4b9xpkzRiCMHm20vxZCiDuRUCjCcnJgwQLjCOH8+f0sWNAVszkbgJCQPvTq9QPJyUYgjBnzBxsTQghkSGqRtmYNbN0KLi6aRYu6k5eXBUBgYGdGjfqBS5dcaN3aOEIQQoh7IaFQROXmQny80dI6I2Mzly8b8xfVrdueF15YTHS0C35+MGQIOMi/shDiHsnpoyLIYoGvv4a1a9dSo0YmM2YMpaAgn/DwQYwYMYedO52pXx969pS2FUKI+yOhUMRYLDB1KmzYEMOyZR0wmZyxWPJo2nQoTz89h61bTQQFwdCh4Odn72qFEEWNhEIR8/XXcPAgHDz4NwAsljwiIsYwZMhU4uIcCAyUQBBCPDg521yEzJ0L+/bB/v1vcODAMgAaNOjK4MGfsnevA15e0LGjBIIQ4sHJkUIRMXUqnDwJGzY8w549swCoUaMp48b9zsGDJsqXh759oVEj+9YphCjaJBSKgGXLjFNG+/ZNtAZCtWphjBu3lR07TAQEGEcIEghCiD9LQqGQS0w07lhOTT3A5s1TAPDy8ueNN3YSE2MiNBSeeALKl7dzoUKIYkFCoRBLTIRvvoG4uCRWrOhEfn4ulSs34K23dpORYcLNzWhsJ4EghHhY5EJzIZWXB/Pnw9atO9m0qTNpaUnUrt2GSZN2cPGiE4cPQ9260LatvSsVQhQncqRQCKWmwuzZsHXrQVatisRszsbHpxZjxiwkN9eNS5egWTPjtJEQQjxMcqRQyJjN8N13sGPHWRYvboTZnI2Dg4mxY3/D0dGH48ehZk149llwdrZ3tUKI4kZCoRDJzjam0ty5M5nFi+tjNufg4ODIq6+uo3z5+uzZAw0awAsv2LtSIURxJaePCon8fPjyS9i27QwLF9bBbM5CKcXYsb/i59eGQ4egYUPjbmVpcCeEsBUJhUIgI8O4hnDggIWVK1tdDQQHXnxxCYGBXdmxA8LCYNw4e1cqhCjuJBTsLDvbOEKIi9McO/YaqakJgOL553+lfv2uREdD/frQv7+9KxVClAQSCnZ0/jz88APs3atJT/8PmzZNwdHRmaefnkNIyGNs2WLcpTxunDFvghBC2JpNP2qUUl2BKYAJmKm1nnzT8leBkUA+kAI8o7VOsGVNhcXZs7BwIaxff5QtW7qRknIMgKeemkVIyCC2b4egIGOSHAkEYUtms5nExERycnLsXYp4CFxcXPDz88PpASdkt9nHjVLKBEwFOgGJQJRSarHW+sB1q+0GwrXWWUqp54F/AYNsVVNhkZ9vHCGsXbubVavakJeXCUCvXn+nadOhbNsGwcEwahR4eNi5WFHsJSYm4uHhQY0aNVAyK1ORprXm4sWLJCYm4u/v/0DbsOU4lqbAUa31ca11HvAj0Pv6FbTW67TWWVcfbgeKfdPnnByj4+mWLbEsX97cGgiPP/4h3br9lS1bIDAQhg+XQBCPRk5ODl5eXhIIxYBSCi8vrz911GfLUKgCnL7uceLV5+7kWWDZ7RYopUYrpaKVUtEpKSkPscRHq6AAZs6E9ev3smpVFyyWPAC6dJlI166vs22bccpo0CCoUMHOxYoSRQKh+Piz/5a2DIXbVaZvu6JSw4Fw4KPbLddaz9Bah2utw318fB5iiY9OZiZ8+ilERSWzenVH0tPPA9Cs2XB69/6A6GgICDCuIVSrZudihRAlli1DIRGoet1jPyDp5pWUUh2Bt4BeWutcG9ZjN+vWwb/+ZbTALlv2EDk5aYAxa9qgQbPZvVtRv77Ry0hmTRMlzcmTJwkKCrJ3GQ/dyJEjOXDAuIT6z3/+84ZlLVu2tEdJ98SWoRAFBCil/JVSzsBgYPH1KyilGgFfYATCeRvWYjfffw8rVsDly2lUrryHGTN6kZ+fS7VqYfTrN58TJ0yEhMBLL0H16vauVgjxsMycOZPAwEDg1lDYunWrPUq6JzYLBa11PjAOWAHEA/O11vuVUu8rpXpdXe0jwB34SSkVq5RafIfNFUm//w47duSxfPmTfPVVLT7+OIKcnCuEhfVn5MjtnDvnQWAgjB0LJpO9qxXCoJS648+MGTOs682YMeOu6z6I48eP06hRI6Kiopg9ezZ9+/ala9euBAQE8Prrr1vX++GHHwgODiYoKIiJEycCMH/+fF599VUApkyZQs2aNQE4duwYrVu3BqBGjRq88847hIWFERwczMGDB2+pYfbs2fTu3ZuuXbtSt25d3nvvPeuyTz75hKCgIIKCgvi///s/ADIzM+nRowchISEEBQUxb948ACIjI4mOjmbSpElkZ2cTGhrKsGHDAHB3dwdg0KBBLF261Lr9ESNGsHDhQiwWCxMmTKBJkyY0bNiQL774AoCzZ88SERFBaGgoQUFBbNq06YH+nu/GpiPgtdZLgaU3Pff2dX/uaMv3tyeLBdavT2DevL6cO7fL+nxoaB+GD/+W48edaNwYRo60Y5FCFCKHDh1i8ODBzJo1i9DQUPbv309sbCy7d++mVKlS1K1blxdffBGTycTEiROJiYmhXLlydO7cmUWLFhEREcFHHxmXJTdt2oSXlxdnzpxh8+bNtGnTxvo+3t7e7Nq1i88//5yPP/6YmTNn3lLLzp072bdvH66urjRp0oQePXqglGLWrFns2LEDrTXNmjWjbdu2HD9+nMqVK7NkyRIA0tLSbtjW5MmT+eyzz4iNjb3lfQYPHsy8efPo3r07eXl5rFmzhmnTpvHVV1/h6elJVFQUubm5tGrVis6dO/Pzzz/TpUsX3nrrLSwWC1lZWbds88+S1mo2UFAA//znKWbMaM+5c7tQyvhrDg5+jP7953PkiAt16hj3IcigD1HYaK3v+DN69GjreqNHj77ruvcjJSWF3r1789133xEaGmp9vkOHDnh6euLi4kJgYCAJCQlERUURGRmJj48Pjo6ODBs2jI0bN+Lr60tGRgbp6emcPn2aoUOHsnHjRjZt2nRDKPTt2xeAxo0bc/LkydvW06lTJ7y8vChdujR9+/Zl8+bNbN68mccffxw3Nzfc3d3p27cvmzZtIjg4mNWrVzNx4kQ2bdqEp6fnPe93t27dWLt2Lbm5uSxbtoyIiAhKly7NypUrmTNnDqGhoTRr1oyLFy9y5MgRmjRpwqxZs3j33XfZu3cvHjYYty6hYAOzZsGsWeNJTz+OyeSM1gVUrhxEv34/kJTkRMOG0txOiOt5enpStWpVtmzZcsPzpUqVsv7ZZDKRn59/18Bp0aIFs2bNom7durRp04ZNmzaxbds2WrVqdcs2r23vdm4+/aWUuuP71qlTh5iYGIKDg3njjTd4//33776z13FxcSEyMpIVK1Ywb948Bg8eDBjB/OmnnxIbG0tsbCwnTpygc+fOREREsHHjRqpUqcITTzzBnDlz7vm97pWEwkN07U7lVavWcuLEIgAsljyqVWtCx44rycpyp3FjGD1ajhCEuJ6zszOLFi1izpw5fP/993ddt1mzZmzYsIELFy5gsVj44YcfaHt1XtqIiAg+/vhjIiIiaNSoEevWraNUqVL39e0dYNWqVaSmppKdnc2iRYto1aoVERERLFq0iKysLDIzM/nll19o06YNSUlJuLq6Mnz4cP7yl7+wa9euW7bn5OSE2Wy+7XtdO2W2adMmunTpAkCXLl2YNm2a9TWHDx8mMzOThIQEKlSowKhRo3j22Wdv+15/lnTVeYhmzIDdu+Hy5WicnV3Jy8uiQoUg2rVbTXh4Gdq3N+5WFkLcys3Njd9//51OnTrh5uZ2x/UqVarEBx98QLt27dBa0717d3r3NpoltGnThtOnTxMREYHJZKJq1arUq1fvvmtp3bo1TzzxBEePHmXo0KGEh4cDxoXgpk2bAsaQ00aNGrFixQomTJiAg4MDTk5OTJs27ZbtjR49moYNGxIWFsbcuXNvWNa5c2eefPJJevXqhfPV6RRHjhzJyZMnCQsLQ2uNj48PixYtYv369Xz00Uc4OTnh7u5ukyMFdb/n/uwtPDxcR0dH27uMW0yfDgcPgo9PJl980ZrTp2MpV64W/fqto2XLqgwdCtcdCQtRaMTHx1O/fn17l1FozJ49m+joaD777DN7l/LAbvdvqpSK0VqH/9Fr5UjhIVizBpYvX8nhw9OBDE6fjqVChdp0776VNm18ePJJe1cohBD3RkLhTzp7Fn7+OYply/qTl5cOgKtrOSIjf6VuXR+uDksWQhQRI0aMYMSIEfYuw24kFP6ExET49NOjfPllZ8xmIxCqVw+ndevvqVs3gD595KY0IUTRIqHwgBISYMGCPL799nnM5ssAVKgQQGTkMvz8vOnZEx7g+pYQQtiVhMIDuHQJFiyw8O9/9+Ts2dUAuLl50abNckJCvOnRA+rUsXORQgjxACQU7tOBA0aDu40b93H27EoATCZnIiOXER5ekyFDwNfXzkUKIcQDkpvX7sP+/fDTTxAXl0lMzFMAeHpWolev1bRt24RnnpFAEOJBXGsQdyeXL1/m888//1PvMXv2bJKSbune/9AkJSXRv39/AGJjY29odLd48WImT558p5cWKhIK9+jiRZg+PYaff55EdPQAzpzZQ4UKAYwYsZ+aNdvQtq3MliaErRSFUKhcuTILFiwAbg2FXr16MWnSJJu998MkoXAPEhPhww+PMGNGJ+Li/sW+fctwd/emU6ff8fAoR+fOcPUmRyGKtLu1wv4zP/cqIyODDh06WFtb//rrrwBMmjSJY8eOERoayoQJEwD46KOPrK2l33nnHcCYsKd+/fqMGjWKBg0a0LlzZ7Kzs1mwYAHR0dEMGzaM0NBQsrOzb3jfyMhIXn75ZVq2bElQUBA7d+4EIDU1lT59+tCwYUOaN29OXFwcABs2bCA0NJTQ0FAaNWpEenq6dbKgvLw83n77bebNm0doaCjz5s1j9uzZjBs3jrS0NGrUqEFBQQEAWVlZVK1aFbPZzLFjx+jatSuNGzemTZs21rbeP/30E0FBQYSEhBAREfEn/nXv0d26HBbGn8aNG+tH6ccftX7uuYPaw6OGxphOVLu6ltW9em3T772n9apVj7QcIR66AwcOWP987Xf8Yf/8ETc3N6211mazWaelpWmttU5JSdG1atXSBQUF+sSJE7pBgwbW9VesWKFHjRqlCwoKtMVi0T169NAbNmzQJ06c0CaTSe/evVtrrfWAAQP0t99+q7XWum3btjoqKuq279+2bVs9cuRIrbXWGzZssL7XuHHj9Lvvvqu11nrNmjU6JCREa631Y489pjdv3qy11jo9PV2bzeYbapw1a5YeO3asdfvXP+7Vq5deu3at1lrrH3/8UT/77LNaa63bt2+vDx8+rLXWevv27bpdu3Zaa62DgoJ0YmKi1lrrS5cu/eHfpdY3/pteA0Tre/iMlQvNdxEXV8CMGf9l7dpXrM95e9ekS5cV1KpVm4EDZdipKF60ndveaK1588032bhxIw4ODpw5c4Zz587dst7KlStZuXIljRo1AowjjCNHjlCtWjX8/f2t7bfv1h77ZkOGDAGMpnpXrlzh8uXLbN68mYULFwLQvn17Ll68SFpaGq1ateLVV19l2LBh9O3bF7/7mEd30KBBzJs3j3bt2vHjjz/ywgsvkJGRwdatWxkwYIB1vdxcY3biVq1aMWLECAYOHGht+21LEgp3sHWr5tlnH+fgwf9NBle6dFkiI5cTFFSbrl0lEIR42ObOnUtKSgoxMTE4OTlRo0YNcnJybllPa80bb7zBc889d8PzJ0+evKXd9s2niu7kXttlK6WYNGkSPXr0YOnSpTRv3pzVq1fj4uJyT+/Tq1cv3njjDVJTU4mJiaF9+/ZkZmZStmzZ207EM336dHbs2MGSJUsIDQ0lNjYWLy+ve3qvByHXFG4jPh5+/70AFxc/6y+Kg4OJyMiFtGwZwMsvw3XzgAghHpK0tDQqVKiAk5MT69atIyEhAQAPDw/S09Ot63Xp0oWvv/6ajIwMAM6cOcP583ef5v3mbdzs2jSamzdvxtPTE09PTyIiIqxdTdevX4+3tzdlypTh2LFjBAcHM3HiRMLDw2+Z1vNu7+Xu7k7Tpk0ZP348jz32GCaTiTJlyuDv789PP/0EGKG3Z88ewJhOtFmzZrz//vt4e3tz+vTpu+7nnyVHCjfJytLMn6/JzHTA2TkNrTUuLp60a/cDHTu257nn4B6/EAgh7tOwYcPo2bMn4eHhhIaGWttee3l50apVK4KCgujWrRsfffQR8fHxtGjRAjA+aL/77jtMd+krM2LECMaMGUPp0qXZtm0bpUuXvmF5uXLlaNmyJVeuXOHrr78G4N133+Xpp5+mYcOGuLq68s033wDwf//3f6xbtw6TyURgYCDdunXj7Nmz1m21a9eOyZMnExoayhtvvHFLLYMGDWLAgAGsX7/e+tzcuXN5/vnn+cc//oHZbGbw4MGEhIQwYcIEjhw5gtaaDh06EBIS8mB/ufdIWmffpFWrkTg7d8TD4zC//fYOjo6lGT58E61bN+aZZ2RyHFH8SOtsY/TRxx9/bJ03oaiT1tkPyQcfbGHr1q/w8FhGenoSoBg8+Efee68x1arZuzohhLA9CYWrvv/+EB9+aEx8YAQC9Ogxhddf7yWBIEQxd/1pnJJOQgFITk7jxRe7k5Z2HFCAplOn9/n55xe5OjueEEKUCCV+9JHWmh49RpGaehylHABNo0ZPsHDhWxIIQogSp8SHwkcfTWPXrp8AhdYF1K7di08++QoPjxL/VyOEKIFK9Cef1pCW1hmTqTSg8fNrzsyZ84mMdLJ3aUIIYRcl+prC9Ol5zJ49CoslG09Pfz7/fBFt25b64xcKUYxNnQoPs5lo5cowduzD297DcvLkSbZu3crQoUPv63UjRozgscces7bJfphatmzJ1q1bH/p270eJDAWtNc899zc2b04gKWk9pUv78uqrK+nZs6K9SxPC7pKSoHr1h7e9qzclFzonT57k+++/v+9QuBf5+fk4Ot7/x6u9AwFK6OmjJUuW8OWX/4/4+O8wmVwYOfI3Xn+9tr3LEqLE+u6772jatCmhoaE899xzWCwWEhISCAgI4MKFCxQUFNCmTRtWrlzJyZMnqVevHk899RQNGzakf//+ZGVlARATE0Pbtm1p3LgxXbp0sd5lfPToUTp27EhISAhhYWEcO3aMSZMmsWnTJkJDQ/nPf/6DxWJhwoQJ1nbcX3zxBWB8iRw3bhyBgYH06NHjju00IiMjefPNN2nbti1TpkwhJSWFfv360aRJE5o0acKWLVsASElJoVOnToSFhfHcc89RvXp1Lly4APxvsiGtNRMmTCAoKIjg4GBrC47169cTGRlJ//79qVevHsOGDXvoTQxLXCikpqYyZMiTVx8pBg6czaRJ4dK6Qgg7iY+PZ968eWzZsoXY2FhMJhNz586levXqTJw4kTFjxvDvf/+bwMBAOnfuDMChQ4cYPXo0cXFxlClThs8//xyz2cyLL77IggULiImJ4ZlnnuGtt94CjPYZY8eOZc+ePWzdupVKlSoxefJk2rRpQ2xsLK+88gpfffUVnp6eREVFERUVxZdffsmJEyf45ZdfOHToEHv37uXLL7+867f5y5cvs2HDBl577TXGjx/PK6+8QlRUFAsXLmTkyJEAvPfee7Rv355du3bx+OOPc+rUqVu28/PPPxMbG8uePXv+f3vnHhtVlcfxz8/SUoIUURFbiwvuIkK1IArZDbVViFDRBasGJGrQVWQRwsOti8QElbgBMRtZBGlERWR3obQshIhkeW5YHisFClhqhFq7bJEIAhbkIX2c/eOeuUxLZ5iZzqPQ3ye5mXPvPfec7/3dx2/OOTO/w/r163nllVdcB1dcXMzs2bMpLS2lvLzcdTbhokV1H124cIH77hvATz+dBGDo0PeZOXMEKSkxFqYoLZgNGzawa9cu+vbtC8C5c+e4yU5j+MILL1BQUEBeXl69CKKdO3emf//+ADz99NPMmTOH7OxsSkpKePDBBwGora0lOTmZ06dPc/jwYXJycgB8RjNdu3Yt+/btc2dPq6qq4uDBg2zevJmRI0cSFxdHSkoKAwYM8HkuI0aMcNPr16+ntLTUXT916hSnT59my5YtrFixAoDs7Gw6dOhwSTlbtmxx6+zUqRNZWVkUFRWRlJREv3793FDdvXv3pqKigoyMDH8mDoqIOgURyQb+AsQBHxpjZjbY3xr4FLgHOA6MMMZUREKLMYbHHnuc0lIn8uCddz7P+PG/138rK0qMMcYwatQoZsyYccm+s2fPUllZCThzJrRr1w7wHeY6LS2N7du319t36tSpgHW89957DB48uN72zz//PODZ49q2beum6+rqGg28F0h3j788DUOD19TUBKQtUCLWfSQiccA84CGgJzBSRHo2yPY8cNIY8yvgXeDtSOlZt24dq1d/BkCnTv2YMCEP+4VCUZQYMnDgQAoLC92++hMnTrghs6dMmcJTTz3F9OnTGT16tHvMoUOH3Jf/kiVLyMjIoHv37hw7dszdXl1dzf79+0lKSiI1NZWVK1cCzuQ1Z8+ebTQc9/z586murgbgwIEDnDlzhszMTJYuXUptbbSXfvMAAAm0SURBVC1Hjhxh06ZNAZ3XoEGDmDt3rrvuaelkZGSwbNkywGmdnDx58pJjMzMzyc/Pp7a2lmPHjrF582b6RWnO30i2FPoBZcaYcgARWQoMA0q98gwD3rDpQmCuiIgJ88hJXV0dubm5ALRu3Z4ZM9bw3HMtqudMUQImJSW8vxi6XPdsz549eeuttxg0aBB1dXXEx8czb948KioqKCoqYuvWrcTFxbF8+XIWLlzIAw88QI8ePVi0aBFjxoyhW7dujB07loSEBAoLC5kwYQJVVVXU1NQwadIk0tLSWLx4MWPGjGHatGnEx8dTUFBAeno6rVq1olevXjz77LNMnDiRiooK+vTpgzGGjh07snLlSnJycti4cSN33XUXt99+O1lZWQGd95w5cxg3bhzp6enU1NSQmZlJXl4er7/+OiNHjiQ/P5+srCySk5PdFpCHnJwctm/fTq9evRARZs2axc0333zJvA0RIZA5O0NZgCdwuow8688AcxvkKQFSvda/AW5spKwXgZ3AzltvvTWgOUq9qa2tNbm5uaZNm/Zm6tSdQR+vKFczjc3n25xpOF/zlcb58+dNdXW1McaYbdu2ufM+h5PmOkdzY51wDVsAgeTBGPMB8AE48ykEK+Saa67hnXfeYfLkyaToqLKiKDHk0KFDDB8+nLq6OhISEliwYEGsJdUjkk6hEujstZ4KNPyfpCdPpYi0AtoDJyIlSB2Colz5dOnShZKSkljLCJlu3bpRXFwcaxk+ieT/FIqAbiLSVUQSgCeBVQ3yrAJG2fQTwEbbzFEUJYroY3f10NRrGTGnYIypAcYD/wS+ApYZY/aLyHQRGWqzfQTcICJlwMvAq5HSoyhK4yQmJnL8+HF1DFcBxhiOHz/u878YgaBzNCtKC6e6uprKykrOnz8faylKGEhMTCQ1NZX4+PrRnnWOZkVRAiI+Pp6uXbvGWobSTGhxsY8URVEU36hTUBRFUVzUKSiKoiguV9xAs4gcA0L9E/6NwA9hlBMuVFdwqK7gaa7aVFdwNEXXL4wxHS+X6YpzCk1BRHYGMvoebVRXcKiu4Gmu2lRXcERDl3YfKYqiKC7qFBRFURSXluYUPoi1AB+oruBQXcHTXLWpruCIuK4WNaagKIqi+KeltRQURVEUP6hTUBRFUVxajFMQkWwR+VpEykQkKtFYRaRCRL4UkT0istNuu15E1onIQfvZwW4XEZlj9e0TkT5e5Yyy+Q+KyChf9fnR8bGIHBWREq9tYdMhIvfY8yyzxwY0y7kPXW+IyGFrsz0iMsRr31Rbx9ciMthre6PX1oZt/8Lqzbch3APR1VlENonIVyKyX0QmNgeb+dEVU5uJSKKI7BCRvVbXm/7KEpHWdr3M7u8Sqt4QdX0iIt962au33R7Nez9ORIpF5LPmYKt6BDI925W+AHE4U33eBiQAe4GeUai3ggbTiwKzgFdt+lXgbZseAqzBmY3u18AXdvv1QLn97GDTHYLUkQn0AUoioQPYAfzGHrMGeKgJut4AchvJ29Net9ZAV3s94/xdW2AZ8KRN5wFjA9SVDPSx6XbAAVt/TG3mR1dMbWbP4Vqbjge+sHZotCzgJSDPpp8E8kPVG6KuT4AnGskfzXv/ZeDvwGf+7B4tW3kvLaWl0A8oM8aUG2MuAEuBYTHSMgxYZNOLgEe9tn9qHP4DXCciycBgYJ0x5oQx5iSwDsgOpkJjzGYundEuLDrsviRjzHbj3K2fepUVii5fDAOWGmN+NsZ8C5ThXNdGr639xjYAKGzkHC+n64gxZrdNn8aZD+QWYmwzP7p8ERWb2fP+ya7G28X4KcvbjoXAQFt3UHqboMsXUbmOIpIKPAx8aNf92T0qtvKmpTiFW4D/ea1X4v9hChcGWCsiu0TkRbutkzHmCDgPOXDTZTRGSnu4dNxi0+HUN9423z8W20UTgq4bgB+NM9lTyLpsc/1unG+ZzcZmDXRBjG1mu0P2AEdxXprf+CnLrd/ur7J1h/0ZaKjLGOOx15+svd4VkdYNdQVYf6jXcTbwR6DOrvuze9Rs5aGlOIXG+vmi8Vvc/saYPsBDwDgRyfST15fGaGsPVke49c0Hfgn0Bo4Af46VLhG5FlgOTDLGnPKXNZraGtEVc5sZY2qNMb1x5mLvB/TwU1bMdInIncBU4A6gL06X0JRo6RKRR4Cjxphd3pv9lBP1+76lOIVKoLPXeirwXaQrNcZ8Zz+PAitwHpbvbbMT+3n0MhojpT1cOiptOiz6jDHf2we5DliAY7NQdP2A0/xv1WB7QIhIPM6L92/GmH/YzTG3WWO6movNrJYfgX/h9Mn7Ksut3+5vj9ONGLFnwEtXtu2GM8aYn4GFhG6vUK5jf2CoiFTgdO0MwGk5NBtbRXSgtbksODPMleMMyHgGX9IiXGdboJ1XehvOWMA71B+snGXTD1N/kGuHuTjI9S3OAFcHm74+BD1dqD+gGzYdQJHN6xlsG9IEXcle6ck4/aYAadQfWCvHGVTzeW2BAuoP3r0UoCbB6R+e3WB7TG3mR1dMbQZ0BK6z6TbAv4FHfJUFjKP+4OmyUPWGqCvZy56zgZkxuvfv5+JAc0xtVU9XsC+XK3XB+WXBAZy+zteiUN9t9oLsBfZ76sTpD9wAHLSfnptLgHlW35fAvV5l/Q5nIKkMeC4ELUtwuhWqcb5JPB9OHcC9QIk9Zi72n/Ih6lps690HrKL+C+81W8fXeP3Kw9e1tddgh9VbALQOUFcGTpN7H7DHLkNibTM/umJqMyAdKLb1lwDT/JUFJNr1Mrv/tlD1hqhro7VXCfBXLv5CKWr3vj32fi46hZjaynvRMBeKoiiKS0sZU1AURVECQJ2CoiiK4qJOQVEURXFRp6AoiqK4qFNQFEVRXNQpKEqEESeKaW6sdShKIKhTUJQgsOGV9blRrlr05laUyyAiXcSZw+B9YDfwkYjs9I7Rb/NViMibIrLbxti/o5GyRovIGhFpE81zUJRAUaegKIHRHSes8t3AH4wx9+L8YzZLRNK98v1gnCCI84F6XUYiMh74LfCoMeZclHQrSlCoU1CUwPivcWLsAwwXkd04IRTScCY88eAJnrcLJ66Th2dwouU+bpxAbIrSLFGnoCiBcQacaRNxWgADjTHpwGqc+DQePC/8WpzgZB5KcJyEd1RNRWl2qFNQlOBIwnEQVSLSCefbfyAUA2OAVSKSEilxitJU1CkoShAYY/bivOD3Ax8DW4M4dgtOK2O1iNwYGYWK0jQ0SqqiKIrioi0FRVEUxUWdgqIoiuKiTkFRFEVxUaegKIqiuKhTUBRFUVzUKSiKoigu6hQURVEUl/8DDIRXC6HCUgQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# convenience plot functions\n",
    "def plot_proxy():\n",
    "    p = plt.Rectangle((0, 0), 0, 0, color='blue', alpha=0.4)\n",
    "    ax = plt.gca()\n",
    "    ax.add_patch(p)\n",
    "    return p\n",
    "\n",
    "def fix_plot_shape(fig):\n",
    "    ax = fig.add_subplot(111, aspect='equal')\n",
    "    axes = fig.gca()\n",
    "    axes.set_xlim([0,1])\n",
    "    axes.set_ylim([0,1])\n",
    "\n",
    "xs = list(range(len(labels)))\n",
    "plt.figure(1)\n",
    "plt.fill_between(xs, list(map(cdf_bounds.lower, xs)), list(map(cdf_bounds.upper, xs)), \n",
    "                 color='blue', alpha=0.4)\n",
    "plt.plot(*zip(*known_pos_ecdf), color='black', linestyle='dashed', linewidth=2)\n",
    "plt.plot(*zip(*latent_pos_ecdf), color='black', linewidth=2)\n",
    "plot_proxy()\n",
    "plt.xlabel('rank')\n",
    "plt.ylabel('TPR')\n",
    "plt.legend(['known positives', 'latent positives', 'expected region'], loc=4)\n",
    "plt.title('Rank CDF')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Compute contingency tables for each rank\n",
    "These tables are computed based on estimate of beta and the rank CDF of known positives. We will show the results for betahat and a the Bayesian credible interval on beta."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first, compute contingency tables based on the point estimate betahat\n",
    "# presorted = True is a computational shortcut\n",
    "# we can use this because we already sorted by decision values earlier\n",
    "tables_point = ss.compute_contingency_tables(labels=sort_labels, decision_values=sort_dv,\n",
    "                                             reference_lb=cdf_bounds.lower,\n",
    "                                             reference_ub=cdf_bounds.upper,\n",
    "                                             beta=betahat, presorted=True)\n",
    "\n",
    "# second, compute contingency tables based on the credible interval \n",
    "# we start with a convenience function for readability\n",
    "compute_tables = lambda beta: ss.compute_contingency_tables(labels=sort_labels, decision_values=sort_dv,\n",
    "                                                   reference_lb=cdf_bounds.lower,\n",
    "                                                   reference_ub=cdf_bounds.upper,\n",
    "                                                   beta=beta, \n",
    "                                                   presorted=True)\n",
    "\n",
    "# tables based on lower bound on beta, only lower bound tables are to be used\n",
    "tables_interval_lb = compute_tables(beta_interval.lower)\n",
    "# tables based on upper bound on beta, only upper bound tables are to be used\n",
    "tables_interval_ub = compute_tables(beta_interval.upper)\n",
    "# merge into a single variable\n",
    "# ss._lb_ub is a collections.namedtuple with fields 'lower' and 'upper'\n",
    "tables_interval = ss._lb_ub(lower=tables_interval_lb.lower, upper=tables_interval_ub.upper)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Compute and plot performance estimates based on the contingency tables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.1 ROC curves\n",
    "\n",
    "Start off with the easy stuff: the true ROC curve (unknown in practice), and the curve we obtain by assuming beta=0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute the true ROC curve (we use Optunity's ROC function)\n",
    "_, roc_true = optunity.metrics.roc_auc(true_labels, decision_values, return_curve=True)\n",
    "# compute the ROC curve that would be obtained by treating all unlabeled as negative\n",
    "_, roc_neg = optunity.metrics.roc_auc(labels, decision_values, return_curve=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute bounds based on the proposed method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can directly use the contingency tables we already computed anyways\n",
    "roc_bounds = lambda tables: ss._lb_ub(lower=ss.roc_from_cts(tables.lower),\n",
    "                                      upper=ss.roc_from_cts(tables.upper))\n",
    "roc_bounds_point = roc_bounds(tables_point)\n",
    "roc_bounds_interval = roc_bounds(tables_interval)\n",
    "\n",
    "# an alternative without all these intermediate steps would be:\n",
    "# roc_bounds_point = ss.roc_bounds(labels, decision_values, beta=betahat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, what if we would compute bounds on the labeled instances and the ones we peeked at to estimate beta? The indices of the random sample of unlabeled data we peeked at is in sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total number of known positives after peeking:\t2135\n",
      "total number of known negatives after peeking:\t165\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\hdemeule\\AppData\\Local\\Continuum\\anaconda3\\lib\\site-packages\\optunity\\metrics.py:211: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return float(FP) / (FP + TN)\n"
     ]
    }
   ],
   "source": [
    "# labels and decision values of peeked instances in sample_labels and sample_decision_values\n",
    "sample_decision_values = [decision_values[i] for i in sample]\n",
    "\n",
    "# labels and decision values of known instances before peeking\n",
    "filtered = filter(lambda x: x[0] is not None, zip(labels, decision_values))\n",
    "apriori_known_labels = list(map(op.itemgetter(0), filtered))\n",
    "apriori_known_dvs = list(map(op.itemgetter(1), filtered))\n",
    "\n",
    "# all known labels and decision values after peeking\n",
    "all_known_labels = apriori_known_labels + sample_labels\n",
    "all_known_dvs = apriori_known_dvs + sample_decision_values\n",
    "\n",
    "# data statistics\n",
    "print('total number of known positives after peeking:\\t%d' % sum(all_known_labels))\n",
    "print('total number of known negatives after peeking:\\t%d' % (len(all_known_labels) - sum(all_known_labels)))\n",
    "\n",
    "# compute the resulting ROC curve\n",
    "_, roc_peek = optunity.metrics.roc_auc(all_known_labels, all_known_dvs, return_curve=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(nan, 0.0),\n",
       " (nan, 0.0033333333333333335),\n",
       " (nan, 0.01),\n",
       " (nan, 0.016666666666666666),\n",
       " (nan, 0.023333333333333334),\n",
       " (nan, 0.03666666666666667),\n",
       " (nan, 0.043333333333333335),\n",
       " (nan, 0.06),\n",
       " (nan, 0.07666666666666666),\n",
       " (nan, 0.09333333333333334),\n",
       " (nan, 0.10333333333333333),\n",
       " (nan, 0.11666666666666667),\n",
       " (nan, 0.13666666666666666),\n",
       " (nan, 0.15),\n",
       " (nan, 0.17),\n",
       " (nan, 0.18666666666666668),\n",
       " (nan, 0.20333333333333334),\n",
       " (nan, 0.22333333333333333),\n",
       " (nan, 0.23666666666666666),\n",
       " (nan, 0.26),\n",
       " (nan, 0.27),\n",
       " (nan, 0.27666666666666667),\n",
       " (nan, 0.2866666666666667),\n",
       " (nan, 0.30666666666666664),\n",
       " (nan, 0.3333333333333333),\n",
       " (nan, 0.35333333333333333),\n",
       " (nan, 0.36666666666666664),\n",
       " (nan, 0.38333333333333336),\n",
       " (nan, 0.4),\n",
       " (nan, 0.4166666666666667),\n",
       " (nan, 0.43666666666666665),\n",
       " (nan, 0.47),\n",
       " (nan, 0.5133333333333333),\n",
       " (nan, 0.5333333333333333),\n",
       " (nan, 0.56),\n",
       " (nan, 0.5933333333333334),\n",
       " (nan, 0.6333333333333333),\n",
       " (nan, 0.6733333333333333),\n",
       " (nan, 0.6966666666666667),\n",
       " (nan, 0.7433333333333333),\n",
       " (nan, 0.7833333333333333),\n",
       " (nan, 0.8333333333333334),\n",
       " (nan, 0.8533333333333334),\n",
       " (nan, 0.8633333333333333),\n",
       " (nan, 0.8833333333333333),\n",
       " (nan, 0.9266666666666666),\n",
       " (nan, 0.96),\n",
       " (nan, 0.98),\n",
       " (nan, 0.9833333333333333),\n",
       " (nan, 0.9866666666666667),\n",
       " (nan, 0.99),\n",
       " (nan, 0.9933333333333333),\n",
       " (nan, 0.9966666666666667),\n",
       " (nan, 1.0)]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "roc_peek"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the resulting curves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\hdemeule\\AppData\\Local\\Continuum\\anaconda3\\lib\\site-packages\\matplotlib\\patches.py:75: UserWarning: Setting the 'color' property will overridethe edgecolor or facecolor properties.\n",
      "  warnings.warn(\"Setting the 'color' property will override\"\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAAEWCAYAAAD1t5d8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzsnXlcVcX7x9/DLosI4obIosiqkmmY5paWablC7qm5Lz+zzPpm6zerb4tpi7ZpbuWWaZqlplmuqam4gCDglogCCoIssnPn98e54GUVzIuo83697ot7z8w585xzD+dzZ+aZ5xFSShQKhUKhuF8xudMGKBQKhUJxJ1FCqFAoFIr7GiWECoVCobivUUKoUCgUivsaJYQKhUKhuK9RQqhQKBSK+5r7VgiFEMOFEL/faTtqEkKIDCFE0zvQrrsQQgohzKq7bWMghIgQQnS9hf3UPalQ3AFqhBAKIc4LIbL0D+IEIcQyIYStMduUUq6UUvYwZhuGCCE6CCF2CCHShRCpQohfhRB+1dV+GfbsEkKMM9wmpbSVUp4zUnteQoi1Qogk/fmHCSFeFEKYGqO9W0UvyJ7/5hhSSn8p5a6btFNK/Kv7nlQoFBo1Qgj19JFS2gIPAK2BV++wPbdEWb0aIUR74HdgI+AMeAChwD5j9MBqWs9KCNEMOAjEAi2llPbAQKAtYHeb27pj517TrntVudvtVyhuGSnlHX8B54HHDD7PBjYbfLYE5gAXgMvAN0Atg/J+wHEgDTgL9NRvtwcWA/HAJeA9wFRf9izwl/79N8CcEjZtBF7Uv3cGfgISgX+AaQb13gbWASv07Y8r4/z2Al+Vsf034Hv9+67AReA1IEl/TYZX5hoY7PsKkAAsBxyATXqbU/TvXfT1/wcUANlABvCFfrsEPPXvlwFfApuBdDQha2ZgTw8gGkgFvgJ2l3Xu+rorDL/PMsrd9W2P0p9fEvC6QXkgcAC4pv8uvwAsDMol8H/AaeAf/bbP0YQ3DTgCdDKob6q/zmf153YEaALs0R/ruv66DNbX7412f10D9gOtSty7rwBhQA5ghsH9rLc9RG/HZeAT/fYL+rYy9K/2GNyT+jr+wHYgWb/va+Vcv1rAXCBG/338pd/WFbhY3v8ape/dt4AswNGgfmv992Gu/zwGiES7p7YBbnf6+aFe6vVvX3fcAClL/XO6ACeAzw3KPwN+ARzRehC/Ah/oywL1//yPo/VwGwM++rKfgQWADVAfOARM1JcVPXSAzvqHptB/dtA/EJz1xzyif0hYAE2Bc8AT+rpvA3lAf33dWiXOzRpNdB4t47xHA/H6912BfOATNNHrgvZA9q7ENSjc9yP9vrWAukCwvn07YC3ws0HbuyghXJQWwmT99TUDVgI/6Muc9A/OIH3Z8/prUJ4QJgCjK/j+3fVtf6u3PQBNVHz15W2Ah/VtuaM9iF8oYfd2/bUp/HHwjP4amAEz9DZY6cteRrvHvAGhb69uyWug//wgcAVohyago9DuV0uDe/c4mpDWMthWeD8fAEbo39sCD5c4ZzODtp7lxj1phyb6MwAr/ed25Vy/L/XfZ2O9jR3090FXbi6Exe5dYAcw3qD+x8A3+vf9gTOAr/66vgHsv9PPD/VSr3/7uuMGSFn0z5mB9utcAn8CdfRlAk0QDHsj7bnxy38B8GkZx2ygf5ga9hyHAjv17w0fOgLtF3pn/efxwA79+3bAhRLHfhVYqn//NrCngnNz0Z+TTxllPYE8/fuuaGJmY1D+I/BmJa5BVyAX/YO+HDseAFIMPu/i5kK4yKDsSSBK/34kcMCgTKD9kChPCPPQ99LLKXfXt+1isO0QMKSc+i8AG0rY3e0m91gKEKB/Hw30K6deSSH8Gni3RJ1ooIvBvTumjPu5UGz2ALMAp3LOuTwhHAocq8T/jgnaj7aAMsq6cnMh3FOifBw37v3C77Xw/+I3YGyJtjNRvUL1ustfNWmOsL+U0g7tn9cHrdcBUA+tV3NECHFNCHEN2KrfDtov8bNlHM8NMAfiDfZbgNYzLIaUUgI/oD18AIah9YAKj+NceAz9cV5DE9pCYis4rxRABzQqo6wR2rBTUV0p5XWDzzFovdKbXQOARCllduEHIYS1EGKBECJGCJGG9kCuU0XnlASD95loPRr0NhWds/76XazgOFcp+/wr1Z7e0WaT3pEqDXifG/dHIcW+AyHEDCFEpN4x5xraMHnhPuXdM2XhBswo8f03QbsGZbZdgrGAFxAlhDgshOhdyXYra6MTWo+xsudTkpK2rwPaCyGc0UZKJNrQPmjX4nOD65CMJpaNb7FthaJGUJOEEAAp5W603sgc/aYktF+8/lLKOvqXvdQca0D7R25WxqFi0XqETgb71ZZS+pfT9GrgaSGEG1ov8CeD4/xjcIw6Uko7KeWThmZXcD7X0YbHBpZRPAit91uIgxDCxuCzKxBXiWtQlg0z0Ib+2kkpa6M91EB7cFVocyWIR+vpagcUQhh+LoM/0IZpb5WvgSiguf5cXuPGeRRSdD5CiE5o83aDAAcpZR204fPCfcq7Z8oiFvhfie/fWkq5uqy2SyKlPC2lHIr2A+wjYJ3+O77Z9a+sjUloc71l1b2O9gMKAP2PoHol6hSzQ0p5Dc2xaxDaD8LV+h86hTZNLHEtakkp91fCToWixlLjhFDPZ8DjQogHpJQ6tLmjT4UQ9QGEEI2FEE/o6y4GRgshugshTPRlPlLKeLR/6LlCiNr6smZCiC5lNSilPIbmWLII2KZ/IIA2RJcmhHhFCFFLCGEqhGghhHioCuczExglhJgmhLATQjgIId5DG96cVaLuLCGEhf5h3htYW4lrUBZ2aOJ5TQjhCPy3RPlltPnOW2Ez0FII0V/vafh/QMMK6v8X6CCE+FgI0VBvv6cQYoUQok4l2rNDm5PMEEL4AJMrUT8f7fs0E0K8BdQ2KF8EvCuEaC40Wgkh6urLSl6Xb4FJQoh2+ro2QoinhBCV8nYVQjwjhKin/w4L76kCvW06yv8ONgENhRAvCCEs9fdNu5KV9MddAnwihHDW35/thRCWwCnASm+vOdqcnmUlzF6FNvwdrH9fyDfAq0IIf/252QshyvqBp1DcVdRIIZRSJgLfo82Pgfbr/gzwt35o7A+03g5SykNoTiefov3q3402hAPaP7MFcBJtiHIdFQ/RrQYew+CfX0pZAPRBm2P7B+0X+CK0obbKns9fwBNoziXxaEOerYGOUsrTBlUT9HbGoQ3NTpJSRt3sGpTDZ2jOD0nA32hDqYZ8jtYDThFCzKvsuejPJwmthzsbbdjTD80zMqec+mfRRN8diBBCpKL1uEPQ5oVvxktovZN0NGFac5P629Dms06hXetsig8BfoI2//o7msAuRrtWoM2bfacf/hskpQxBmzP+Au27OYM2l1dZeqKdcwbaNR8ipcyWUmaiee/u07f1sOFOUsp0NAewPmj3xWng0XLaeAnN+ecw2nDlR4CJlDIVmIJ2v15C6yFWNIRdyC9Ac+CylDLUwKYN+mP/oL8Hw4FelTieQlGjETdGPRR3EqFFIlkhpaxoiLFGIoQwQXvADpdS7rzT9igUCkVVqJE9QkXNRwjxhBCijn4IrnDO7u87bJZCoVBUGaMJoRBiiRDiihAivJxyIYSYJ4Q4I7RwWw8ayxaFUWiP5qmYhDZ8119KmXVnTVIoFIqqY7ShUSFEZ7S1gd9LKVuUUf4k8Bza+rR2aAvoSzkDKBQKhUJhTIzWI5RS7kGbuC+PfmgiKaWUf6OtcavMWjOFQqFQKG4bdzLIbmOKe/Jd1G+LL1lRCDEBmABgY2PTxsfHp1oMVCjuB/LzC4iLiyMtLZ2cnMLRbZEspa5uhTsqFPcId1IISy6IhnIWGUspFwILAdq2bStDQkKMaZdCcU9z/TrExMDixetYuvQlUlJiipV7ez9BdPQfMeXsrlDcc9xJIbyIFkaqEBe09XMKhcIIpKdLli0L5ccf13P48Hxycq4VlQlhyiOPvM7gwS8wZowDNjZCdwdNVSiqlTsphL8AU4UQP6A5y6Tqo8EoFIrbQFYWJCVJduw4xPr1P7Fnz09cu3Yj77IQpri7d2LkyCkMHToQa2to1AjMVFZCxX2G0W55IcRqtADaTkKIi2hhtswBpJTfAFvQPEbPoAVYHm0sWxSKe5mcHG2oMyMDUlLg4sUCDhzYx5EjPxEdvZ709BvBZMzMLGnXbhTOzsFMndqVRx6xwLQqYdgVinsQowmhPtBwReUSLUalQqEwoKAApNRe+fmQl6e9NyQnB65dg3Pn4MABHfHxkVy6dIC4uP2cObOZjIwrRXWFMEXKAhwd3ejRYzZuboMICoLAwKrZdeTIkfpmZmaLgBaoYByKuwsdEJ6fnz+uTZs2V0oWqkEQhaIGoNNpvbpDh+DECW2blCD0LmXCwLUsO/saly4d5OLFA3rx+5vs7LRix3NyaoqtbQMSEqJwc3sUX9+n8fUdhI+PKe3bQ/PmVbfRzMxsUcOGDX3r1auXYmJiomIzKu4adDqdSExM9EtISFgE9C1ZroRQobiD6HRw9ixs3QpxcWBjAy4uYGJSWK4jISGKc+cOFL0SEiIpLxCGlVUdnnxyEd7eQVha5tC0qSVNmwoaNYIGDbTj/wtaKBFU3I2YmJjIevXqpSYkJJQK7gJKCBUKo6PTwaVLEB8P//yjDWtKCampkJSkDX06OICHB+h0BcTGHicq6k9OndrFuXMHyMy8Vux4JibmuLm1oW7dZhw//iO+vk/i7R1E/fqdcXd345FHBN7e4ORkdbtPxUSJoOJuRX/vljmkr4RQoTAChUJ34QLs2AGJidrwprX1Da9Mc3No1EiSlHSa0NA/iYz8g+jonWRmphQ7lp1dY5ycWqDTZZGQcAw/v74MGbICgP79v8DRsQ4ODvDQQ9C06Y3epEKhqBxKCBWKWyQ5Ga5e1Tw109O1nl56uvY5MfFGz8/BAdzcbux37VocUVF/Fr1SUoqnCLS3d8Pd/TE6deqGtfVVdu9ey8GD2zAzMyMoKIgpU8bTpSi9dGXyGt8bJCUlmS5atMhx5syZiXfaFsW9hRJChaIKZGVBeDgcPKgNdRY6sZiZgamp9tfSEhwdtR4fQHp6IkeO7ObUqd1ER/9JfHxksWPa2DjRtGk3XFy607hxAI6Oobz22ljs7U2ZMWMGGRnJvP/++4wZM4YGDRpU8xnXHK5evWq6ePHi+mUJYX5+PmZ3cAFkXl4e5oVfuOKuQwmhQlEJkpPhyBH46y9tSUPJXp4haWlXCAvbw6lTuzh1ahdxcRHFyi0sbGjSpDPu7t3x938MHx9Pjh1bTUjIKrZtm4xOp6NvX38eeeQRPvjgA+bOnVsNZ1jzmTFjhktsbKylj4+PX5cuXdL69OmT+u677zaqX79+3smTJ623bNlyunfv3s1Pnz4dAfDWW281yMjIMP3kk0/iIiIiLCdNmuSanJxsZmVlpVu0aFFM69atsw2Pn5qaajJ27FjXsLAwa4DXXnst7tlnn71mbW3dOjMz8xjA0qVLHTZt2mT/008/nQ8ODnZ3cHDIP3HihHWLFi2ytm7dWic0NPSkk5NTAYCrq2uLffv2RZmamjJ69Gi3S5cuWQB88sknF3r06HG9eq+eoiKUECoUBly7pone9euQna31AM+dg1OntN5egwY3enqFpKVd4fTp3URHa8IXH3+yWLm5uRWurh1o3Lgrbm6P0r59IF26WODmBrGxUXTu7EFiYiLe3t68/vrrDBo0CH9/fwAsLCyq69SrhBCijTGOK6U8Ul7Z3LlzL/bu3btWVFTUSYBNmzbZhYWF2Rw7dizCx8cnNzo6utyLNW7cOLeFCxfGtGzZMmfHjh02kydPdv37779PGdaZOXNmo9q1axecOnXqJEBiYuJNQw2cPXvWat++fafMzMwYPXo0K1eurPP8889f3bFjh42Li0tukyZN8vv06ePx4osvXn7iiScyTp8+bfHEE080P3fuXMTNjq2oPpQQKu57MjIgNhb+/hvOnNGcTXQ6bdjTxERzcHF1veGEUlCQx9mz+zlxYgvh4VuIiyuee9rMrBYuLh1wc+uKm1tXXFwewtvbklattF5kevoFoqKi8PPrgaenJ3369GHkyJF07twZIcqKRa8oj1atWl338fHJrahOamqqybFjx2wHDhzYrHBbbm5uqQu9Z8+e2j/88ENRDLp69eoV3Kz9oKCglMIh2WHDhiW/8847zs8///zVlStXOgYHBycD7Nu3r/bp06drFe6TkZFhmpKSYuLg4KDiudYQlBAq7htyciAzU3t//bq2fu/4cbiijzNhZwdNmpTtdZmWdoWIiN84cWILJ09uIysrtajMzKwWTZo8QosWXejYsSudOj1E3bqWWFhArVrg5KTNHx4/fpxp0+ayevVqGjRowIULFzAzM2Px4sXVcPa3l4p6btWJtbV1kZiYmZlJne6GtmRnZ5sAFBQUYGdnl1/YkywPKWWZP0QMt2VlZRWrYGtrW9Rg9+7dr48dO9YyLi7ObOvWrXX+97//xRUeNyQkJNLW1lYtPamhKCFU3NMUFGg9vZCQG0sYQPPmNDEpf65Pp9Nx4cJRTpzYTHj4FmJiDhdbxN6ggTfNmz9FkyZPMmlSRzp0sCw3WPXOnTt55ZVXOHz4MDY2NkybNo0XXngBUxXks0rY29sXXL9+vdzFIS4uLvnJyclmCQkJpvb29rpt27bZd+/ePc3R0VHn4uKSu2TJEocxY8ak6HQ6Dh48WKt9+/ZZhvt37do17ZNPPqm/ZMmSWNCGRuvVq1dQt27dvKNHj1oFBARkb9y40cHW1rbMnqKJiQm9evW6NmXKlCaenp5ZDRs2LADo2LFj2kcffVT/3XffvQywf//+Wh06dMgq6xiKO4MSQsU9S1ISbNigzfHVr6/19ioaeSwoyCM6eifHjq3n+PGNpKUlFJWZmVni6fko3t5P4uHxJFZWzXBxgQEDwNm5+HGklBw5cgQnJyfc3d0RQpCZmcnnn3/OiBEjcHBwMNIZ39s0bNiwoE2bNhnNmzf379atW2qfPn1SDcstLS3ljBkz4gMDA31dXFxyPD09i5xhVq9efW78+PFuH330UaP8/HwxYMCA5JJC+MEHH8SPHj3atXnz5v4mJibytddeixs1atS1WbNmXerXr59no0aN8nx8fLIqEuPhw4cnd+nSxXfevHnnC7ctXLgwdty4ca5eXl5+BQUFol27dukdOnS4cBsvjeJfIsoL1VRTUYl5FRUhJUREwIEDcP68NjRZv3759XNzs4iM3M7Roz8RFvZLsSguDg5NaNr0KZo1e4qmTR+lcWMbatUCe3to21brSRoKq06nY9myZcyfP5/jx4/zwgsv8Omnnxb1JO+m+T8hxBEpZVvDbaGhoecDAgKS7pRNCsW/JTQ01CkgIMC95HbVI1TcM+h0WszO3buhbl3NwaUs7cnOTufEiS0cO7ae8PDN5OTc8GRv2NAPb+9gPD2DaNAggA4dBA89pB2vopHMf/75hzFjxrBr1y4CAgL46quvGDZsGHB3CaBCcT+ihFBx1yOlNv+3ezccPQru7qVFKysrlePHN3L06DpOnvyd/PycojJn5zb4+gbTrNkAmjb1wc9P6+25ukKdSgZu+fDDDzly5AiLFi1izJgxSvwUirsIJYSKu5LkZG3oMyZGW+OXlqat83Nzu+H1mZNznbCwTYSE/EB4+G9F4ieEwMWlIw88EMRTTwXh7OyGmRn4+GgiWtlYnRcuXCAzMxMfHx8++ugjXn/9dVxdXY1yvgqFwngoIVTcVVy6pEV3CQ3Vhj1r1dKWPRT6n+TlZRMa+hshIWsIC/uV3FxtvYQQAi+vLvj5DcLTcwDTpzcq5eRSWaSULFmyhOnTpxMQEMDevXupU6cOdSrbfVQoFDUKJYSKu4K4OC2LQ0RE6QXu+fm5nDjxByEhP3D8+M9kZ6cX7efh8TCtWg3B3X0g1tbOWFjAhAmlPT0ry8WLFxk/fjxbt26la9euLFmy5DacnUKhuJMoIVTUSHQ6LcRZcrI29xcRofX+Coc+8/JyiIzcxdGj6zh69KdiqYtcXR+kbdvBtGkzCBsbdxIToWdPTUAbN9bCpN0KoaGhdO7cmfz8fObPn8+UKVMwUTmPFIq7HiWEijtGUpIW3SUyEi5f1oQvL08rK1zwrtOBlZXWA7x+PZGDB7cQFvYrERHbyMnJKDqWs7M/bdsOoW3bwTRo0BzQFtOfPw/Dh0OrVlW3LyUlhVWrVmFtbc3o0aPx9fVl4MCBzJw5E09Pz9twBRRVITo62sIwqPbNmDdvXt2+ffumubu75/3btl999dWGK1eudDIxMWHu3LkXgoOD0/7tMRU1ByWEimonLQ127tRSGQkBtWuDjY32t7CDJYQ2FxcfH8mxY78SFvYr584dQMobIbRcXFrRqlVf2rYdTOPGLQBN/JKTtVBq2dnQuXPVRFCn07Fz504WL17M+vXrycnJoX///owePRoLCwsWLVp0Oy+FwoisWLHC6YEHHsj6t0J45MgRq/Xr1ztGR0dHxMTEmD/++ONe/fr1C7+TaZ8Utxf1TSqqjawsLdzZzp2a0DVpUnqZQ0FBHqdP7yUsTBO/xMSzRWWmpuZ4ez9Gq1Z9aNWqD3XruhU79uXLmueop6c2hFq/PjRvXjUbJ06cyKJFi3BwcGD8+PGMGTOG1q1b/5vTVtxG8vPzCQoKcg8PD7du2rRp9tq1a88fP37c6sUXX2ySmZlp4uDgkL9y5crzO3bssA0PD7ceOXJkUysrK11ISEjk22+/3XDr1q11cnJyTNq2bZuxcuXKmMoMba9bt65OUFBQcq1ataSPj0+um5tbzq5du2wee+wxlUrpHkEJoeK2U7iuLzlZ+5yVBdHREBWl5fJr1Kh4KiOdTsepU7vYv38pYWG/FgtobWNTl5YtnyIgoC9+fj2wsrIr1V5qqvYaNAj8/KAqmYuklKxfv5727dvj7OzM6NGj6datGwMGDMDKyupWL8F9QWBgoHfJbUFBQckzZ85MTE9PN+nevXupnyHPPPNM0rRp067Gx8eb9evXr5lh2aFDh6Jv1ub58+etFixYcL5Hjx7XBw4c6D579ux6mzZtcti8efMZZ2fn/G+//dbhpZdearx27drzX3/9df05c+bEdu7cORPg5ZdfvjJnzpx4gP79+3v88MMP9sOGDUt98803G6xdu7ZuybYefvjh9GXLlsVeunTJ4uGHHy4ah3d2ds6NjY21AJQQ3iMoIVTcVo4ehe3bteFPIbSXTqcNfdarV1wAU1IuceDAMvbtW0JSUlH2Gxo18qVVq760atWHpk0fxsSkdEgXwyFQGxuYNAlcXKpm6+7du/nPf/7DoUOHePPNN3nnnXfo0KEDHTp0uNXTVxiZhg0b5hYmtR0xYsTVDz/8sNHp06drdevWzQu0H1X16tUrcyj0t99+s/vkk08aZmdnm1y7ds3Mz88vC0h99913LxcGxC6LssJQCiHurtiUigpRQqi4bRw6BOvXQ8OG5WdvLyjIIzT0V/btW0xExNaiOT8HhyZ06DCahx8eQb16nly/rjnOpKRoPUzQBLWgQEunZG6uLYBv1Qo8PDSP0MoSFhbGq6++ypYtW3BxcWHx4sWMHDnyX579/UdFPTg7OztdReWNGjXKr0wPsCQlI/bY2NgUeHp6Zh0/fjyqov0yMzPFjBkz3A4ePHjS09Mz78UXX3QuTNN0sx6hi4tLYQ8QgLi4OAsXF5d/7YCjqDkoIVT8a3Q6be7vl1+0XllZQ5MJCVH89ddi/v77e9LTtQSApqbmBAQE0bHjOHx9H8PExJSsLM3Ts1EjrQdpYXGjZ2lhccOD1MNDe38rzJkzh/379zN79mymTp1KrVq1br6TokYQHx9v8ccff9g89thj11etWuUYGBh4ffny5U6F23JycsSJEycs27Ztm21ra1uQmppqCpCZmWkC0LBhw/zU1FSTX3/91aFPnz4pADfrEQYHB18bPnx407feeutyTEyM+fnz5626du2qhkXvIZQQKm6Z7GzNQWXLFrhwobQI6nQ6Dh1axZ4933D27L6i7Y0a+dGx4zjatXsGO7t6RdsTErTe3+DBEBBQccqkqiCl5Mcff6Rly5b4+fnx8ccf8/nnn6t0SHchTZs2zV6yZEndKVOmuHl4eOTMnDkz9qmnnkqdNm2aa3p6umlBQYGYPHny5bZt22aPHDky6bnnnnN7+eWXdSEhIZHDhw9P9PPz83dxcckNCAiotJC1bds2u3///sleXl7+pqamfPLJJzHKY/TeQqVhUlSaixe18GbJyXDtmpblHcDWVsvOYMiFC8dYtWoy//xzEABLS1seemgIjzwyFg+PdsWGuPLytNBpTZvC009XPtB1Zbh8+TJTpkxh/fr1TJgwgQULFty+g9/DqDRMinsRlYZJccsUFGgCuG2bNhdnba2t+XNyKl03KyuNX355i5075yOlDnv7RvTpM4uHHhqKlZUteXmaR2lurnZc0Ob7evSATp0qTnVUFQp7gf/3f/9Heno6H374ITNmzLg9B1coFPcUSggVN+W33zQhbNKkuNenIVJKQkJ+ZO3a6aSmxiOECd27v0CfPrOoVas2eXlapghzc2jZUpv/q1VL++vioq3/u50sX76cUaNG8dBDD7Fs2TL8/PxubwMKheKeQQmhokIiImDv3rJz/BVy+fIpVq+eSmTkdkALdD18+Nc0afIAoK0dvHABnnwS2rev2jq/qpKfn4+ZmRkDBw4kIyODCRMmoOZzFApFRagnhKJcEhNh7VrNg7MsEczLy+a33z5g27YPyc/PxdragQEDPqRjx3FFwagLRfCJJ6BLF+Pau3LlSj744AN2795N3bp1mTJlinEbVCgU9wRKCBWliIjQhkIvXLgxJ1iShIQovv12CBcvhgLQvv2zBAfPLuYFmpqqOdb07GlcEZRSMmvWLGbNmkWXLl1UdniFQlEllBAqihETAytWgKOjtl6vpKZIKdm/fxk//DCV3NxM6tVrxqhRS2nSpBMZGZo3qZQ3gmlPnqwdx1hkZ2czduxYVq1axahRo1i4cCEWxhx7VSgU9xwqmZqiiOvXYfVqbSmEvX1pEczKSmPJkmeSWQKdAAAgAElEQVT4/vsx5OZmEhg4jNdfP0qdOp1ITtYCXD/xhJb2aPJkeP5544ogwMsvv8yqVat4//33Wbp0qRLB+5QlS5Y4NG3a1L9du3Ze+/fvr7VmzRr7f3tMa2vrCqOtR0dHWzRv3ty/KscMDg52X7p0qQPA4MGD3Y4cOaIC2tYAjNojFEL0BD4HTIFFUsoPS5S7At8BdfR1ZkoptxjTJkXZ5ObCzz9rAbIbNy5dfv78YRYtGkpi4lksLKwZOvRL2rcfxcWLgvr1Ydiw0msJq4M333yT7t27079//+pvXFFjWLp0qdPnn39+oU+fPunz5s2rGxISYjN48ODUm++pkZeXh3l5LtFGYs2aNTHV2qCiXIzWIxRCmAJfAr0AP2CoEKKkD/sbwI9SytbAEOArY9mjKJ+kJFiwQEuQ6+xcvEyn0/H773P46KMOJCaepUmTB3j99aN06PAsiYmCxo1hwoTqFcGoqCjGjh1Lbm4u9evXVyJ4H/HYY4818/f39/X09PSfM2eOE8BLL73U6MiRI7bPPfec29ixY5t88MEHzr/++quDj4+P37fffuuQlpZmMnDgQPcWLVr4+vr6+q1YsaIOaIl7e/Xq1bRbt26enTp18iqvzdTUVJP27dt7+fn5+Xp5eRXtDzfSQnl5efn17NmzaXp6ugnA3r17rR966CFvf39/344dOzaPiYkppbKBgYHee/bssQat9/ncc8819vb29gsICPCJjY01A4iIiLAMCAjwadGihe8LL7zgfLNequLWMGaPMBA4I6U8ByCE+AHoB5w0qCOB2vr39kCcEe1RlCA+XssW8fffN2J4GpKXl82SJSM4enQdAN26TSMoaDbm5pbk5GjBr59+Giwtq8/m8PBwunfvjhCCixcv0rRp0+prXHEDIdoY5bhSHqmoeOXKlecbNGhQkJGRIVq3bu33zDPPpMyZMyd+z549tQtTLs2bNy8zJCTE5vvvv78AMHXq1MaPPvpo2tq1a88nJSWZtm3b1rdv375pAEePHrUNCwuLaNCgQUF5bVpbW+s2b958xtHRURcfH2/Wrl07n2HDhl2D0mmhPv7443qvv/76lWnTprmWlRqqvDaysrJM2rdvnzF//vxLkyZNcpk/f3692bNnx0+dOrXJlClTrkycODF59uzZ9crbX/HvMKYQNgZiDT5fBNqVqPM28LsQ4jnABnisrAMJISYAEwBcjT3pdJ9QGCTb3FzLFlFyVCgrK42vvurHqVO7qFXLntGjvycgoC+gLYm4dAmCgsqOLmMsjh07xuOPP46lpSU7duxQIngf8tFHHzXYvHlzHYCEhATziIgIq4YNG1YYN3TXrl21t23bVmfevHkNAXJycsSZM2csADp16pRWkQgC6HQ68cILL7j8/ffftiYmJly5csXi4sWLZlA6LdS8efPqh4WFpVY2NVQh5ubmcsiQIakAbdq0uf7HH3/UBjh27Jjt77//fgZg3LhxV99+++0qJhtTVAZjCmFZPuwlA5sOBZZJKecKIdoDy4UQLWRhbp7CnaRcCCwELdaoUay9jzh/XhNBZ+eyF7enpV1m3rxexMYew96+EdOmbcXFpRV5eRAXp60p7NwZ2rYtva+xOHz4MD169KB27drs2LGDZs2a3XwnhfG4Sc/NGGzatMlu9+7ddiEhIVF2dna6wMBA76ysrJtO70gpWbdu3ZmAgIAcw+1//fWXjbW1ta68/QpZsGCB49WrV81OnDgRaWlpKRs3btyysN2SS3WEEEgpRWVSQxliZmYmC9fempmZkZ+fr9YAVSPG9Bq9CDQx+OxC6aHPscCPAFLKA4AVUI19jPsLnU7ryRUujyhLBBMTzzF79iPExh6jfv3m/Oc/+3BxaQVo+3bpAq++qkWJMalGn2NTU1OaNWvGnj17lAjep1y7ds3U3t6+wM7OTnfs2DGr0NBQm7Lq1a5duyAjI6Po7nz00UfT5s6d20Cn0zRv3759Vcq7lZqaaurk5JRnaWkpf/31V7u4uLii/5zCtFAAq1atcuzQoUNGq1atspOTk80Kt+fk5IiQkJBb8g594IEHMpYtW+YAsGTJEsdbOYbi5hjzUXYYaC6E8BBCWKA5w/xSos4FoDuAEMIXTQgTjWjTfcvu3fDuu/DVV5qA1a5duk5sbCizZz9CYuJZXF0f5OWX/8LJyQOAjAxtSUXXrlqM0OoiJkZzrHvwwQc5fPgwbuVl/FXc8wQHB6fm5+cLLy8vv9dee825vFRKvXr1Sj916lStQmeZDz/8MC4/P1/4+Pj4NW/e3P+NN94owy+6fMaNG5ccGhpq06JFC98VK1Y4enh4ZBeWFaaF8vLy8ktJSTF76aWXEq2srOQPP/xwdubMmS7e3t5+/v7+frt377a9lXOeP39+7Pz58xu0bNnSNz4+3tzW1rbCYVzFrWHUNExCiCeBz9CWRiyRUv5PCPEOECKl/EXvRfotYIs2bPofKeXvFR1TpWGqOufPwzffVBw0+9Sp3Xz5ZV+ys9Pw9u7G5MkbqFXrhlr+8w+MGAHVGbt63759PPHEE8ydO5eJEydWX8MKlYaphpCenm5iY2OjMzExYeHChQ5r1qxx/PPPP8/eabvuVu5IGib9msAtJba9ZfD+JPCIMW2438nKgh9/1JY3lCeC4eG/8fXXA8jPz6FNm4GMHr0cc3PNFVSn0+YFPTzAx6f67A4JCeHJJ5+kcePG9OvXr/oaVihqEPv27bN+/vnnXaWU1K5du2DZsmXn77RN9yIqxNo9jE6npVBKT9d6g2WhiWB/8vNz6dx5EkOHfoGJiRZhOyVFewUEVO+cYFhYGD169KBu3br8+eefNGzYsHoaVihqGD179syIjo4+efOain+DEsJ7FClh+3Y4dEhLoVQWhiLYtetUhgyZhxCC5GQtYLaLi7ZOsDp9U9LT0+nRowc2Njb8+eefuLgob3GFQmFclBDeg0gJu3bBjh2aCJbVkytPBDMytHWCkyaVHXTb2NjZ2fHpp5/Spk0bPDw8qrdxhUJxX6KE8B4jKwt+/VWLGOPmVnYewfJEUKeDK1dg7Fht3+okJiaGs2fP0q1bN4YOHVq9jSsUivsaJYT3ECkpsHSplgqpqj1B0NYJPvKIlkWiOomLi6N79+5cv36ds2fPYl1WAkSFQqEwEkoI7xEyMmDZMsjM1Ob2yuLUqT3limByMjg4wOOPV5/NAFeuXKF79+5cvnyZ7du3KxG8i/jgA+rFxXHbUjY4O5P36qs1bx1xdHS0xc6dO20nTZqUXJX9goOD3Xv37p06evTolNttU+vWrX2OHTtW6cg1VWXw4MFu//nPfy63adMm++a1tQDihbFeK1N///79tWJjYy1uliGkMJNIYdzYyvDOO+/Unz59epKdnd1NowYVooTwHiA7W4sWk5paOntEIVevxrBgQTD5+bl06TK5mAhmZ2sCOnq0Fny7ukhOTubxxx8nJiaGrVu38vDDD1df44p/TVwc5m5u5N6u48XEUCOTSZ4+fdpyzZo1jlUVwspwq+mfjCmCYPwUUSEhIdZVTZVVWRYsWNBg/PjxyVURQpWY9y6noAB++kkb1ixPBHNyrvPVV/3IyEjCz68HgwffEMHCANrBwdCgQTUaDixYsIDo6Gg2btxI586dq7dxxV3JV1995diyZUtfHx8fv2HDhrnl5+dz6tQpCzc3txbx8fFmBQUFtGnTxnv9+vW1o6OjLTw8PPyrkiYpPDzcskOHDl7e3t5+fn5+vhEREZavv/5645CQEFsfHx+/WbNm1c/Pz2fixIkuLVq08PXy8vL7+OOPnUALrj1y5EjXZs2a+Xft2tUzKSmpzI5GYGCg99SpUxs/9NBD3u+9916DuLg4syeeeKJZixYtfFu0aOH7+++/2wDExcWZdejQobmfn5/vsGHD3JydnVvGx8ebwY2kwTqdjokTJ7o0b97c38vLy+/bb791AC0ua2BgoHfPnj2benh4+Pft29ejMMRcIUePHrVq2bKlb+Hn6OhoCy8vL79CGwtTRA0fPty1RYsWvp6env7Tp08v5ykDy5Ytq9u6dWuf5s2b++/cudMaoKwUWNnZ2aJkqqydO3dat27d2sfX19evdevWPqGhoUU5bRISEsw7derU3M3NrcWkSZOKxrvKsuu9996rf+XKFfMuXbp4tWvXrtzUWiVRQngXI6W2TjA8vPzhUCkl3303mosXQ6lf35Nx437A1FT7/8zIgAsXoGdPba1gdfPKK69w6NAhHq/u8VjFXcnRo0et1q1b5xgSEhIVFRV10sTERH7zzTd1vby8cp9//vmE0aNHu7799tsNvL29s4OCgtJAS5M0adKkxFOnTp20s7PTffzxx/VycnLEtGnTXDdu3Hg2IiIictSoUUkvvfRSY4Bhw4Z5TJo06Up0dPTJkJCQKFdX17z//e9/l9q2bZsRFRV18r///e+Vzz77zMne3r4gPDw8MjQ0NPK7776rFxUVZbF8+fI6Z86csYyOjo5YtmxZzNGjR8sNq3bt2jXTw4cPR8+aNevyxIkTm7z44ouXw8PDIzds2HB20qRJ7gAzZ8507tKlS/rJkycjg4KCUuLj40v1mL///vs6J06cqBUZGRnx559/nnrrrbdcCkU9MjKy1pdffhl75syZiAsXLlhu3769mD0PPvhgdl5enjh58qSF/liO/fv3LzWM+8knn1wKDw+PjIqKiti3b5/dwYMHywyymJmZaXLs2LGoefPmxUyYMMED4LXXXmv06KOPpoWHh0fu3bs3+o033nDJzc0Vr776alyfPn1SoqKiTo4fPz4lICAg+9ChQ1GRkZEn//vf/176z3/+U/REO3nypPXPP/98LjIyMuKXX35xOHPmjHl5dr3xxhtX6tevn7d79+5TBw8ePHWTW6oINTR6F3PgAOzdq0V9KW+Zw2+/fcCRI2uxsrJjypRfsLFxALRkvAUFMH589a4TzMrKYsqUKbz11lt4eHjQqlWr6mtccVezdetWu/DwcOuAgABfgOzsbJP69evnA7z44otJ69evd1i2bFm9sLCwogXoVUmTlJKSYnL58mWLkSNHXgOwtraWlM6Ywx9//FE7KirK+pdffnEASE9PNz158qTV7t277QYNGpRsZmaGu7t7Xvv27dPLO5ehQ4cWDbPu27ev9unTp4vEJSMjwzQlJcXk0KFDtj///PMZgKeffjqtdu3apeKM7t27t6jNJk2a5Ldr1y7jr7/+sra3t9e1bNnyerNmzfIA/P39M8+ePVtKSPv375+8YsUKx/fffz9hw4YNDmvWrDlXss53333nuGzZMqf8/HyRmJhoHhoaatWuXbuskvWGDRuWDNCrV6+MjIwMk6SkJNOKUmAZkpycbDp48GCP8+fPWwkhZF5eXtETrWPHjml169YtAPD09Mw+e/aspaenZ15l7aoMSgjvUk6d0lIpubqWH/ElNPRXfvnlDYQQjBmzkkaNtFGQtDRtSHTyZKhXjak+c3JyCA4OZuvWrfTq1UutE1RUCSmlGDhw4NUvv/zyUsmy9PR0k4SEBAuAtLQ0UwcHBx1ULU1ScnJypUbIpJRi7ty5F4KDg9MMt2/atMm+ZHvlYTh/JaUkJCQk0tbWtpjoViYOdEV1LC0tiwpNTU3LTO00YsSIlIEDBzYdMmRIihCCli1bFktVFRUVZfHFF180OHLkSGS9evUKgoOD3bOzs8u8TuVc63JTYBl+fuWVVxp36dIlffv27Wejo6MtunXr5l1YZmFhYXgeMi8vT1TFrsqghkbvQi5fhpUrtTm98ubZ4+JOsmTJcKSU9O37HgEBfQBtnWFKCjz7bPWLYFBQEL/99hsLFy5k0KBB1de44p6gZ8+eaZs2bXK4dOmSGcDly5dNT506ZQFaFvqnn3766muvvRb37LPPFq2CrUqaJEdHR13Dhg1zly9fXgcgKytLpKenm9jb2xdkZGQUrch9/PHHU7/++ut6OTk5AiAsLMwyLS3NpEuXLulr1651zM/PJyYmxvzvv/+2q8x5dezYMe2jjz6qX/h5//79tQACAwMzli9f7giwfv362mlpaaVWBXfp0iV93bp1jvn5+cTFxZkdOnTItlOnThUmKjbE398/x8TEhLfeest5wIABpZyBUlJSTGvVqqVzdHQsiI2NNdu1a5d9ecdavXq1A8C2bdts7ezsCurWrVtQXgqskqmy0tLSTF1cXHIBFixYcNNUfBXZZWNjU5CamlolbVM9wruM8HDNOaZWLbApMxsbXL+ewtdf9yM7O502bQbRq9erAOTlQXw8PPOM1pOsLgpFcMuWLSxYsIBx48ZVX+MKo+HsTN7t9PR0dqbCLO5t2rTJfuONNy51797dS6fTYW5uLufNm3fh9OnTFsePH7dZvHhxlJmZGRs2bHD4/PPP6/bs2TO9ME3SlClT3Dw8PHIM0yRNmzbNNT093bSgoEBMnjz5ctu2bbNXrFjxz/jx493effddZ3Nzc7l27dqzgYGBWWZmZtLb29tv2LBhSW+88caV8+fPW7Zs2dJXSikcHR3ztmzZcnbEiBHX/vzzz9re3t7+Hh4e2YGBgeUOjRqycOHC2HHjxrl6eXn5FRQUiHbt2qV36NDhwocffhj39NNPN/Xz83No3759Rr169fLq1KlTbHh0xIgR1/bv32/r6+vrL4SQs2bNuujq6pofFhZW6eseFBSU/O6777p89NFHpXra7du3z2rRokVm8+bN/V1dXXPatGmTUd5xHBwcClq3bu2TkZFhunDhwn8APvzww7gJEya4+vj4+EkphYuLS87OnTvP9OrVK33OnDmNfHx8/GbMmBH/yiuvJIwbN85j3rx5DTt16pRWXhuVsWvUqFFJvXr1al6/fv28ys4TGjUNkzG4X9MwZWbC1q1a7NCGDaG85XYFBfl88cVTnDz5O02aPMDLL/+FpaUNOp2WSql3b+jUqXptT01NpUePHowdO5YJEyZUb+OKW+JeSMMUHR1t0bt37+anT5+OuNO23ApZWVnCzMxMmpub88cff9hMnTrVLSoqSgXg/hfckTRMitvD+fOwZg1cv15+xJhCNmyYycmTv2Nr68TkyT9jaal1G2NjNQHs2LFaTAYgOzsbKSX29vbs27cPMzN1uykUleXMmTMWgwYNalbY+12wYMH5O23TvYp6MtVwrl+HJUu0jPI3S8Tw99/L2b59LiYmZkyc+BN162pTJZcva2mYevasviDa2dnZBAcHU1BQwJYtW5QIKqodb2/v3Lu1NwjQsmXLnMjISNUDrAaUs0wN5+JFbZlD7doV1zt//jDLl48HYMiQ+Xh5aQvUs7I0D9GBA6G6tKhQBLds2UJwcDAm1ZXIUKFQKG4B9YSq4Zw+DZaWFddJTY3XxxDNoXPniXTpMgnQFszHx0NQkJahvrp48cUX2bJlCwsXLmT8+PHV17BCoVDcAmq8qgYjpeYlWqdO+XXy8nL45ptgrl2Lw9OzE4MHz0On08Km2djAqFHg61v+/rebPXv28PXXXzN9+nQlggqF4q5ACWENJjlZ69U5OpZdLqVk1arJnDt3AAeHJkycuA4zMwvOn9dCpvXrV71BtAFsbW3p27cv7777bvU2rFAoFLeIGhqtwcTGVly+bdts9u9firl5LaZM+Znatetz8SJ4eWnDodUtggAPPvggGzduxKa8RY4KxV1GdHS0xTfffFPOz9HyCQ4Odl+6dKmDMWxq3bq1jzGOW8jgwYPdjhw5UuknyAsvvOD8888/VxhAYNOmTXbbt283+oNh06ZNdo8++qhnVfZRQliDOXmy/PWCR46sY8OGmQghGDt2JY0aPcj585pn6ZAh5UecMRZHjx5l8uTJpKXddC2sQnFXUZiGyRjHzsurMIZAuVRHGqbK5iIE+Oyzz+L69+9fYQCBHTt22O3du7fcQORlcavXp6ooIayh5Odr8UTLmh/8559DLF06AoABAz7Cy2sAcXHQqxeMGVP9PcG8vDzGjBnDxo0bKZnqRaG4nag0TDUzDZNh77dx48Ytp0+f7uzn5+fr5eXld+zYMavo6GiL77//vt4333zTwMfHx2/r1q225Z33iy++6Dx06FC3Rx55pHlQUJBHq1atfEJCQoqeaoGBgd579+61rih1U1VRQlhDOXQIcnJKL3m4ejWGr77qS15eNh07jqNHj5e4elWbD+zUqfqWSBjy8ccfExoayldffUWdijx7FIp/gUrDVLPTMBni5OSUf/LkycgxY8Ykfvjhhw28vb1zR44cmThp0qTLUVFRJ3v27JlR3nkDhIWFWW/btu3Mr7/++k9wcHDyypUrHQFiYmLMr1y5Yt6pU6fMilI3VRXlLFMDOXxYyyzRpEnx7VlZaXzxRW/S0i7j49OdYcO+Kor47u5e/XYCREZGMmvWLAYOHEj//v3vjBGK+wKVhqlmp2EyZNiwYSkAgYGBmYXXqSTlnTdAz549rxVm4xg5cmTKY4895vXpp5/Gff/99w59+vRJgYpTN1UVJYQ1jOPHtaDaLi7F5/kKCvL59tvBxMWF07ChDxMnrsPU1JycHG2ZRHVmkjDk+eefx8bGhvnz598ZAxT3DSoNU5m2lFtW3WmYDLGyspIAZmZmsqx2C20v67wBbGxsiq6Ph4dHXp06dfIPHjxYa/369Y4LFiyIgYpTN1UVNTRagzh0CH74AZydwaLEb7effvoPERFbsbV1YurUzVhba0OQKSng7199odNKsnDhQtasWUODBg3ujAGK+waVhqlmp2G6GXZ2dgXp6elF51DeeZfF008/nfz+++83TE9PNw0MDMyCqqduqgglhDWEo0dh/XqtJ1jS2eXvv5fz55+fYmpqzuTJG6hXr2lRWU4OeN/y76BbQ59sE51Oh7u7O48//nj1GqC4LzFMw+Tl5eXXrVs3r9jYWPPNmzfbHj9+3Oa9995LmDx5crK5ubn8/PPP6wIUpmHy8vLyS0lJMTNMwzRz5kwXb29vP39/f7/du3fbAqxYseKfL7/8sr6Xl5df27ZtfWJjY80M0zDNmjWr/vTp05N8fHyyW7Zs6du8eXP/8ePHu+Xl5YkRI0Zca9q0aY63t7f/2LFjXauShuno0aM2Xl5efs2aNfP/4osv6oGWwmjHjh21/fz8fDdv3mxfXhomf3//LF9fX/+uXbt6FaZhqsp1DQoKSt64caPjiBEjSs0PGqY7GjFihHtFaZhuRnBw8LXNmzfXKXSWKe+8y+KZZ55J2bx5s2O/fv2KxPqVV15JePvtt10efPBBn4KCUqPGVUKlYaoBZGbCnDmah2hJEbxw4SizZz9CXl42w4d/Q+fOE4vKpNTWGr76avm5CW83UkpmzJjBp59+yurVqxkyZEj1NKyoVlQapjuPSsN0+1FpmGowBw5Abm5pEUxPT+TrrwfoPUTHFxNBrRwaN65eEZw+fTqff/4506ZNY/DgwdXTsEJxH6LSMFUfSgjvMKmpsHs3NGpUfHtBQR4LFw4iOfkCTZu2Z8iQ0s4oV6/Cww9Xj51SSp577jm+/PJLpk+fzty5c0s5IigUNQmVhklRWdQc4R1ESvjzT83RpWQkmHXrXubUqV3Urt2QiRPXYW5efK1oYqImnm3bUi2cPHmSRYsW8dJLLykRvH/R6XQ69cUr7kr0926ZET9Uj/AOcvCg9vLwKL79wIHv2bHjc0xNzZk06Sfq1CkezCEjQ4s8M3x49UWR8ff3JzQ0FC8vLyWC9y/hiYmJfvXq1Us1MTG5u5wLFPc1Op1OJCYm2gPhZZUrIbxDnD0LGzdqi+YN89aePx/CihUTABgy5AuaNetQbL/sbK03OG5c9eQY3Lp1KwkJCTz77LN4V7d7qqJGkZ+fPy4hIWFRQkJCC9RokuLuQgeE5+fnjyur0KhCKIToCXwOmAKLpJQfllFnEPA2WgSHUCnlMGPaVBNIS4OVK7VF8IbrBdPSrvDNNwOKEux27jyh2H55eRAXB4MHQ7NmxrfzzJkzDB06FDc3N4YNG4ZFycWNivuKNm3aXAH63mk7FIrbjdGEUAhhCnwJPA5cBA4LIX6RUp40qNMceBV4REqZIoSoX/bR7h2k1MKn6XRgaxD5T3OOGUhKykWaNevA4MHziu1XUKAtlXjqKWjd2vh2pqen079/f0xMTNiwYYMSQYVCcc9izOGNQOCMlPKclDIX+AHoV6LOeOBLKWUKgJTyihHtqRGEhWlZ50t6ia5dO4PTp/dQp45zUYLdQnQ6iImB7t2hY0fj26jT6Rg1ahSRkZH8+OOPeJScxFQoFIp7CGMKYWPAMLXsRf02Q7wALyHEPiHE3/qh1FIIISYIIUKEECGJiYlGMtf4pKbCzz9rImjob7J//zJ27pyPmZkFEyf+hL39DZWUEs6f1wTwsceqJ5Tazp072bBhA3PmzKF79+7Gb1ChUCjuIMacIyzrkV3S08wMaA50BVyAvUKIFlLKa8V2knIhsBC0yDK331TjIyX8+qv2vpZBRL3o6J2sXKktlB869EuaNi2+MDAxEXx84Mknqy+eaPfu3dm3bx/t27evngYVCoXiDmLMHuFFwDCRkAsQV0adjVLKPCnlP0A0mjDecxQOiTZseGPbhQvH+OqrfuTn59Kt2zQ6dizu0KTTaeHXnniiuGepsThy5Aj79+8HoEOHDmqZhEKhuC8w5uP1MNBcCOEhhLAAhgC/lKjzM/AogBDCCW2otFQ+rLudK1dKD4kmJp5l/vxeZGen07btYAYO/LTUfpcvQ5s2pecTjUFsbCy9e/dmzJgx5OdXKWavQqFQ3NUYbWhUSpkvhJgKbENbPrFEShkhhHgHCJFS/qIv6yGEOAkUAC9LKa8ay6Y7QWgorFsH1tY3hkRTUxP47LMepKVdxtf3MZ599jtMSnT58vK016OPGt/G9PR0evfuTWZmJn/88QdmdyLNvUKhUNwhjPrEk1JuAbaU2PaWwXsJvKh/3VNICb//Djt2aIGxCyPAZGWlMn9+Lx90pdsAACAASURBVJKSzuHm1pZJk9aXCp+Wnw8XLkCfPuDoaFw78/PzGTx4MBEREfz222/4+/sbt0GFQqGoYaif/kbAUATd3cFUn4qyoCCPBQueJjb2OPXrN+e557ZgZVU8d2d+vrZU4sknq2epxNKlS/ntt99YsGCByiuoUCjuS5QQGoH9+0uLoJSS1aunEhn5B7VrN+D553/Hzq50HspLl7RlEl26VI+tY8aMoVGjRvTu3bt6GlQoFIoahooXeJvJy4OdO7VM84UiCPDHH5+yd+9CzM2tmDLlF5yc3Evtm5amJec1tghKKfnkk0+4cOECpqamSgQVCsV9jRLC28zp09qSB0uDab/Q0F/46aeXABg9+ns8PAJL7SclJCVB//7F44/ebqSUvPnmm8yYMYNFixYZryGFQqG4S1BCeJvZswfs7W98jo09zuLFw5BS0q/fe7RpM7DM/S5dghYtoLmRV1G+/fbb/O9//2PcuHG8/fbbxm1MoVAo7gKUEN5GEhK0wNgODtrnzMxrfP11f3JyrvPwwyPo1eu1UvsUFGjOMY0bQ18jx/WfPXs277zzDmPGjGHBggWllmwoFArF/YhylrlNXL8Oe/feyDQvpWTVqilcvRqDm1tbnnnm21KRWnQ6LY7oo49qDjLGXL6Xm5vLmjVrGDx4MN9++60SQYVCodCjhPBfkpurOcfs368tfXBx0bYfPLiCw4dXY2lpw9ixq0qtFQRtrWDHjtCzzFDjtxcLCwv27NmDiYmJEkGFQqEwoMpPRCGEqRBiuDGMudvIyIDvvoPdu6F+fXBz0zxFExPPsmrVFAAGD55HgwalJ/6uXNGy0z/xhHFtTEhIYNKkSWRkZGBjY0Mtw4jfCoVCoShfCIUQtYUQrwohvhBC9BAaz6HFAh1UfSbWPNLT4cgR+PJLLWO8u/uNIdGCgjwWLx5OTk4GDz74NB06jC61f0EBZGdrmeaN6SGanZ1N//79Wb58OefO3XMhXBUKheK2UNHQ6HIgBTgAjANeBiyAflLK49VgW43k1ClYvlxb7uDoeMMxppCNG9/kn38O4uDgwjPPLCwzg0NSEjzwgHHDp0kpGTduHAcPHuSnn36iVatWxmtMoVAo7mIqEsKmUsqWAEKIRUAS4CqlTK8Wy2ogOTmwYYMmfra2pcv37FnAtm0fIYQJo0cvx8bGoXQltN5gmzbGszM3N5cxY8awcuVK3n33XYKCgozXmEKhUNzlVDRHmFf4RkpZAPxzP4sgwF9/adFfyhLB48d/LpoXHD78G7y9u5Z5jKwsbX83N+PZmZCQwPbt23n//fd5/fXXjdeQQqFQ3ANU1CMMEEKkcSPTfC2Dz1JKWdvo1tUQsrK0Be87doCzc+nyM2f2sWjRUKTU0bv323TqNL7cYyUlaUsljOG4mZiYiJOTE66urkRFReFQctxWoVAoFKUoVwillKblld0vFBTATz9pOQWFADu7G04xhcTHR/Lll33Iy8umY8fx9O79VtkHQxPUggIwRqajEydO0LNnTyZMmMB///tfJYIKhUJRScoVQiGEFTAJ8ATC0BLr3jepy3U62LQJjh3TvELL8HkhOzuDr7/uT2ZmCgEBfRk27KsynWOSk28E1H78cXByur22HjhwgCeffBJra2sGDBhwew+uUCgU9zgVDY1+hzZPuBd4EvAHnq8Oo2oCu3bBgQPliyDAjz8+z+XLp3B29mfcuNWYmpa+nNeuaRkppk4te1j13/L7778zYMAAnJ2d2b59O+7u7re/EYVCobiHqUgI/Qy8RhcDh6rHpDvPtWvafKCra/lzeYcP/8C+fUswN7di/Pg1WFhYFyvX6bT5wPx8mDDBOCKYlJREUFAQnp6ebNu2jYYNG97+RhQKheIepyIhNPQazS9ryO9e5fBhrRdYXuzPpKR/WLFiIgCDBn2Gs3PxSb8rV7QlEn5+WhzRRo2MY6eTkxPr168nMDCQOnXqGKcRhUKhuMepSAgf0HuJguYpel94jWZmasskyutc5efn/n979x6dVXXue/z7hBCEcItELiECUQiKoEJRYQ+5iFbBzWUfZAjequdYUbS0XvbWqrVibbW7dW+PjnIOFVtbRPFCjwoUKx6VgpaLSCBguIiAkA0hXEICCYGEzP3Hyu1NQgjkXXnf9eb3GSNjZM213rUeMgIPc645n8mrr95KcXEBAwdOZNiwqSHni4u9XuDDD/u3YH7Hjh1s3LiRcePGcf311/vzEBGRZqK+RLjeOTewySKJEhkZXiKrOTu0wrvvPsKOHStJSjqfO+6ovaNETo63nZJfSbCkpIRbb72VzZs3s337ds0OFRFppPoSoWuyKKJExU4SXbrUff4f//gTS5f+jvj4BO69910SE0OzXUGBV3XGz6oxM2bMYOXKlbz99ttKgiIiYVBfIuxsZg+f6qRz7j99iCeiMjO9tX6dO9c+t3PnGt544z4Abrnl/5CWdhXg1Rw9dMjbiQLgzjtP3ZtsrM8++4znn3+eu+++m5tvbtZ1z0VEwqa+RNgCaEtVZZmYVloKn3wC551X+1xBQS6zZk2ktPQ4w4ffx9VX3w14yS83F/r0gTFjvF3m6/p8OBQWFnL77beTnp7OSy+95M9DRESaofoS4V7n3C+aLJII27QJ8vO9dYPVOeeYM+du8vJ2c8EFQ5k82UtChYXeUOgPfgAXX3zqtYbhkpiYyNNPP80VV1xBYmKivw8TEWlG6kuEzaInCF7ZsyVL6q74snLl62zYsIjWrTswdeo7xMd7GwgePAjjxnlLJPzknOO7776jV69eTJ069fQfEBGRM1Jf6edrmyyKCFu/3lv8XnNXicOH9/DOO14xnZtv/t8kJaUCXuIEf2qG1vTCCy9wySWXkJWV5f/DRESaoVMmQufcoaYMJFKOHoW//rX2ukHnHHPnTqWo6DADBvwzQ4feWXnuwAG47LK6t2MKp4ULF/LYY48xduxYLrroIn8fJiLSTPmwGVCwfPKJVwu0devQ9rVr57Nhw19p3boDt932+5D1gsXFcMUV/sa1fft27rjjDgYOHMhrr71GnB/7NomISPNOhNnZsHJl7TqgJ04cY/78fwVg4sR/Jympe+W5o0e9tYI9evgXV3FxMZMmTcLMmD9/Pm3atDn9h0RE5Kw020RYVuYNibZrV7uw9scfv8ChQ7tITb2Mq6/+YWW7c7B/P9xwgz8b61aIi4tj2LBhzJkzh7S0NP8eJCIi9c4ajWlbtsDOnVAzzxw6tJsPP3we8CbIxMVV7U+ck+NNkOnf37+4ysrKSEhI0FpBEZEm0ix7hCUl3qa7dS2XeO+9xykpOcagQZPo23dkZXtxsdcjHDvWvzWDWVlZXHrppWRmZvrzABERqaVZJsKvvoK8PG9YtLrs7ExWr36D+PgEbrrptyHn9u2DG2/0dpn3Q2FhIZMmTSI3N5dOnTr58xAREaml2Q2NFhTARx/VvUfgwoVPAzB8+H0kJ/eqbC8q8pZKDPRpLw7nHPfddx+bN29myZIldO/e/fQfEhGRsPC1R2hmo81si5ltM7Of1nPdJDNzZjbYz3jAWy5RVgatWoW279y5hnXr3qdly9aMHv14yLmKCTJ+FdOePXs2c+fOZcaMGVx33XX+PEREROrkWyI0sxbATGAM0A+4xcxqFSQzs3bAj4FVfsVSITsbVq8+VW/w5wBcc82P6NChanV9RW9wwAB/YnLO8f7773P99dfzs5/9zJ+HiIjIKfk5NHolsM05tx3AzN4CJgA1a4U9C/wG+FcfYwHg44/rXi6xY8cqNm78kFat2nLDDY+GnMvNhZtu8q83aGYsWLCAoqIiLZoXEYkAP//l7Q7srnacXd5WycwGAuc75xbVdyMzm2pma8xszf79+88qmMOHYdu2uneO//TTlwGvN9i2bdVU0rw8b5Peyy47q0fWq7i4mOnTp5Obm0t8fDzt27cP/0NEROS0/EyEdS0yqNz13szigBeBR053I+fcK865wc65weed5YZ/mzZ5yx5qdrqOHNnP2rXzMTOGD7+vsr2szEue//IvEB/mfrNzjnvuuYff/e53rFixIrw3FxGRM+JnIswGzq92nArsqXbcDugPLDWzncAQYIEfE2acgxUr6u4NfvHFHyktPUH//v9Mp049K9v37vXqifbsWfszjfXrX/+auXPn8uyzzzJhwoTwP0BERBrMz0T4JdDHzNLMLAGYAiyoOOmcy3fOJTvnejnnegErgfHOuTXhDiQnx9s/sOZuEWVlZSxf/nsARoyYVtleXOz1HP2YwLlhwwaeeuopJk+ezJNPPhn+B4iIyBnxLRE650qBHwEfAZuAd5xzX5vZL8xsvF/PrcvGjdCiRe32rKwlHDiwg06denHJJTdUtufkeMsl/HhtN2PGDJKSkpg5c2bIjhYiIhIZvi6od84tBhbXaPv5Ka4d6UcMRUXesGjnzrXPLVv2fwEYPvzeypqihw971w72aUXja6+9RlZWlqrHiIhEiZifr/+Pf3i1RRMSQtsPHdpFZuYiWrRoyT/90/8CvHeJeXkwfnz4J8jk5+dTUlJC+/btGTJkSHhvLiIiZy2mE2FeHixdWnv3eYDly2fjXBmDBt1E+/ZedzE/35scc8EF4Y9l2rRpDBkyhNLS0vDfXEREzlpMJ8JPPvEWwtdcDH/yZAmff/4qAMOHV02Syc+HoUPDH8drr73GvHnzGD9+PPHh7mqKiEijxGwiPHoUMjK8BfE1rVv3PgUFOaSkXEKfPsMAOHnSmynap09448jIyOD+++9n1KhRmiUqIhKFYjYRHjhQ9wJ6gL//fRbg7TJRMXPz4EGvnmibNuGLIS8vj5tuuolOnToxb9489QZFRKJQzCbCnJy6N9DNydnCli2fkpDQhiFD7qhsP3YMvve98MZQUFBAcnIy8+fPp3Nd01ZFRCTiYraLsn07JCbWbl+2zOsNXnnlbbRu3QHwFtC3bh3+KjI9e/Zk1apVWi8oIhLFYrZHuHNn7UoyJSXHWbHiTwCMGFFVVzQnx6siE66Ry/3793P//feTm5urJCgiEuViMhEeOQKFhbXXDm7e/AlFRYdJTb2MHj0GAf4soH/88ceZPXs2Bw8eDN9NRUTEFzGZCPfvr/v9YEbG/wNg0KCbAG+Hibw8mDAhfL3BVatW8Yc//IGHHnqIiy++ODw3FRER38RkIszJqd128mQp69d/AMDAgRMBb2Zp//6Qlhae5548eZIHHniAlJQUnnrqqfDcVEREfBWTk2W2b6/9fnDbts85evQAXbqk061bP8CbJHP55eF77quvvspXX33Fm2++Sbt27cJ3YxER8U1MJsKdO6Fjx9C2imHRgQMnhkxg6dYtfM8dN24cubm5TJkyJXw3FRERX8VcIiwo8Hp61SfKlJWVhSRCgBMnvCUTSUnhea5zTkOiIiIBFHPvCA8dqj1R5rvv1nD48H+RlJRKz57e9ND8fK+cWjhWN3z44YeMHDmSnLpeToqISFSLuR5hfr43G7S6uoZFi4ogPb3xzyssLGTatGm0adOGc889t/E3FBGRJhVziXDfvtDdJpxzZGT8BYDLL/8fle1mkJLS+Oc988wzfPfddyxbtoyEmgsXRUQk6sXc0Ojevd67v6rjLHJzt9G2bTK9e18NQGmplyyTkxv3rN27d/PSSy9x1113MWzYsMbdTEREIiLmEuG+faGJMCvrYwD697+RFi28DnBBAVx4Yd07U5yJF198EfB6hSIiEkwxlQhLSrwk16pVVdvWrZ8B0LfvNZVtR49C376Nf96zzz7LokWL6NGjR+NvJiIiERFTiTA/33v3VzETtKzsJFu3/h0ITYTQ+PeDzjkSExP5/ve/37gbiYhIRMVcInSu6nj37nUcO5ZPcnIanTp5eyxVzCg977yzf86mTZvo168fa9eubUS0IiISDWIqER4+HJoIt2ype1g0NbX2zhRnYsaMGWRnZ2tIVEQkBsRUIszJgXPOqTretu1zANLTR1a2HTniLaQ/W6tXr+add97hwQcfJLmx005FRCTiYioR7t0bmgizs9cDVFaTATh5Es62I1dWVsb06dPp2rUrjz76aGNCFRGRKBFTC+pzcqpqhx47ls/BgzuJj29F585VXUDnoEuXs7v/ggULWL16NX/+85+1u4SISIyImURYXAzHjnm7zQNkZ2cCkJLSv3L9YHExtG/vfZ2N8ePH8/777zNu3LhwhCwiIlEgZoZGK5ZOVKgYFj3//Msq244cgd69z+7+xcXFxMXFMWHCBOIauxJfRESiRsz8i56XF3pckQhTU6sSYVHR2U2U2bp1Kz169GDJkiWNCVFERKJQzCTCrKzQijK7dnlr/KonQrMzfz/onGP69OkcP36cSy+9NByhiohIFImJd4QnTkBmZtX7wRMnjpGdnYlZHD16fA/wZou2aHHmhbbfffddlixZwssvv0zXrl3DHLmIiERaTPQId+70kmF8eVrfvTuDsrJSUlIu4Zxz2gJeDdK0NC8ZNlRBQQEPPvggAwcOZNq0aeEPXEREIi4meoTr10ObNlXHO3asBiAt7arKtqNH4Ux3Slq4cCE5OTm89957xMfHxI9KRERqCPy/7sePw8aNUH3UcufOVQD06nVlZZtzXmm1M3HbbbdxxRVXkB6OrexFRCQq+To0amajzWyLmW0zs5/Wcf5hM8sys0wz+8TMep7pM3bs8LZfqj7kuX37CqCqR1hRaLuhE2VKSkrIzPTWISoJiojENt8SoZm1AGYCY4B+wC1m1q/GZRnAYOfcpcB84Ddn+pwVK6Bt26rjAwd2cPDgd7Rpk0RKSn/AWzbRtWvorNL6PPfccwwaNIgtW7acaTgiIhIwfvYIrwS2Oee2O+dOAG8BE6pf4Jz7zDlXVH64EjijwcvcXPjmG+jUqaqtYseJ9PQRlQvfCwoavhFvRkYGv/zlL5kyZQp9w7F7r4iIRDU/E2F3YHe14+zytlO5G/iwrhNmNtXM1pjZmv3791e2f/kltGwZWlFm69alQOjWS6Wl0LMBg64nTpzgzjvvJDk5mZdffvn0HxARkcDzc7KM1dHm6mjDzG4HBgMj6jrvnHsFeAVg8ODBDqCwEFatCp0k45yrtQdhxfvBhiwBfOaZZ9iwYQMLFy7k3HPPPf0HREQk8PxMhNnA+dWOU4E9NS8ys+uAJ4ERzrnjDb15ZqaX5Kqvati//1vy8rJJTOxEt26XAHDoEKSnN6zQdps2bbjnnnsYO3ZsQ8MQEZGA8zMRfgn0MbM04L+AKcCt1S8ws4HA74HRzrncM7n5unXQsWNoW1VvcGTl+8HCQhg6tGH3fPLJJ3Guzk6riIjEKN/eETrnSoEfAR8Bm4B3nHNfm9kvzGx8+WW/BdoC75rZOjNb0JB7l5V5ew8mJoa2V7wfTE/3hkWPHfNmlF5wQb1xMnXqVP72t78BYFbXiK6IiMQqXxfUO+cWA4trtP282vfXnc19jxypqh1a7V4hPUKA/fth9Oj6y6rNmjWL2bNn069fP0aPHn024YiISIAFstbo4cO12/bt20p+/l7atetMt279KC312gcMOPV9du3axSOPPMKYMWP4yU9+4k+wIiIS1QKZCPPyvJJp1VUNi44EjF27YNQo6NDh1Pd54okncM4xa9YsDYmKiDRTgUyEe/fWrhJTfVh0zx7o3x9Gjjz1PTIzM3njjTd4+OGH6dGjh3/BiohIVAtk0e3s7NDdJpxzlT3C1NRrSEyEiRMhrp40P2DAAN566y1uvPFGf4MVEZGoFshEuHdvaFm1nJzNFBTso337rrRs2ZchQ0ITZU1lZWXExcUxefJk/4MVEZGoFrih0bIybxPeli2r2kJni1q95dSKi4sZNGgQr7/+uq9xiohIMAQuEVbMBq2uIhH26XMNZtCt26k//9JLL7F+/XpSUlJ8ilBERIIkcInw5MnQ47Kyssr3g927X0OPHpCQUPdn9+3bx69+9SvGjRvHtdde62+gIiISCIFLhCUlofVF9+7N4ujRA3TsmEJCQm/61dzxsJqnn36aY8eO8cILL/gfqIiIBEIgE2H1iTBV+w9eAxinWgmRnZ3Nq6++ygMPPKBd50VEpFLgZo2WlEDr1lXH1d8PxsWderul1NRUli9frs12RUQkROASoXNV6wPLyk7yzTd/ByAlZSRpaaGzSSvk5+fToUMHhjZ0GwoREWk2Ajc0Wt327SspLDzEeeddSMuWF9T5fjA3N5fevXszc+bMpg9QRESiXqAT4YYNiwC49NJxgHHhhbWveeihh8jPz2fUqFFNG5yIiARCoBNhZuZCANLTx3LuuZCcHHp+8eLFvPnmmzzxxBNcfPHFEYhQRESiXWAT4YEDO9mz52vOOac9SUnDGDgw9PzBgwf54Q9/SP/+/Xn88ccjE6SIiES9wE2WqbB58/8H4KKLrgUSqLkiYvXq1Rw9epTFixfTquZWFSIiIuUCnAg/BaBPn2tp3bp2WbUxY8awa9cuOnbsGIHoREQkKAI5NOqcY8sWLxF26TKKAQOgRQvv3I4dO5g3bx6AkqCIiJxWIBPhnj0by7dd6kZi4kWV7wdPnDjB5MmTmTZtGgcOHIhskCIiEgiBHBrduNFbNtGz52jGjKkqq/bEE0/w5Zdf8pe//IXkmlNIRURE6hDIHuHGjd6yiZEjxzFsmNeWlZXFiy++yL333svEiRMjGJ2IiASJOeciHcMZ6dbtcrdvXyZxcS355puDpKW1BWDChAksXbqUb7/9Vr1BkUYys6+cc4MjHYdIUwjc0GhJSRHOOVJSriQ1tW1l+6233srYsWOVBEVE5IwELhGWlh4HoEuX9JAC25MnT45QRCIiEmSBe0dYkQh79eoNwKJFi3juuec4fvx4JMMSEZGAClwiPHmyGIC+fXtTUlLCI488wty5c2lRsZBQRETkDAR2aHTAgN7MmTOHrVu38sEHHxAfH7g/ioiIRIHAzRo1i3PgWL/+MFOmDKVVq1asXbsWM4t0aCIxQ7NGpTkJYDfKkZjYhc2bv2DTpk28/vrrSoIiInLWAveOECApqTc9e3bmjjvu0GxRERFplAD2CKFbt95cddVgrrpqTqRDERGRgAtkjzAhoYxt27ZFOgwREYkBgUyEGRkf8Oijj0Y6DBERiQGBTIRFRQX8+Mc/jnQYIiISA3xNhGY22sy2mNk2M/tpHedbmdnb5edXmVmvhtz3ssuuYMSIEeEOV0REmiHfEqGZtQBmAmOAfsAtZtavxmV3A3nOud7Ai8C/N+Tejz32pJZMiIhIWPjZI7wS2Oac2+6cOwG8BUyocc0E4M/l388HrrUGZLhRo4aENVAREWm+/Fw+0R3YXe04G7jqVNc450rNLB/oBByofpGZTQWmlh8e79q160ZfIvZHMjX+PFEsSLFCsOINUqwAfSMdgEhT8TMR1tWzq1nPrSHX4Jx7BXgFwMzWBKn0U5DiDVKsEKx4gxQrePFGOgaRpuLn0Gg2cH6141Rgz6muMbN4oANwyMeYREREQviZCL8E+phZmpklAFOABTWuWQDcWf79JOBTF7Qq4CIiEmi+DY2Wv/P7EfAR0AL4o3PuazP7BbDGObcA+APwupltw+sJTmnArV/xK2afBCneIMUKwYo3SLFC8OIVOWuB24ZJREQknAJZWUZERCRclAhFRKRZi9pE6Fd5Nj80INaHzSzLzDLN7BMz6xmJOKvFU2+81a6bZGbOzCI27b8hsZrZzeU/36/N7M2mjrFGLKf7XehhZp+ZWUb578ONkYizPJY/mlmumdW5Ltc8L5f/WTLNbFBTxyjSJJxzUfeFN7nmW+ACIAFYD/Srcc39wKzy76cAb0dxrNcAbcq/nxapWBsab/l17YBlwEpgcLTGCvQBMoCk8uPO0fyzxZuEMq38+37AzgjGOxwYBGw8xfkbgQ/x1vsOAVZFKlZ96cvPr2jtEfpWns0Hp43VOfeZc66o/HAl3prKSGnIzxbgWeA3QHFTBldDQ2K9B5jpnMsDcM7lNnGM1TUkXge0L/++A7XX1jYZ59wy6l+3OwGY4zwrgY5m1q1pohNpOtGaCOsqz9b9VNc450qBivJsTa0hsVZ3N97/siPltPGa2UDgfOfcoqYMrA4N+dmmA+lm9oWZrTSz0U0WXW0NiXcGcLuZZQOLgelNE9pZOdPfbZFA8rPEWmOErTxbE2hwHGZ2OzAYiOQeUvXGa2ZxeDuB3NVUAdWjIT/beLzh0ZF4Pe3lZtbfOXfY59jq0pB4bwH+5Jz7DzMbireOtr9zrsz/8M5YtPwdE/FVtPYIg1SerSGxYmbXAU8C451zx5sotrqcLt52QH9gqZntxHs3tCBCE2Ya+nvwgXOuxDm3A9iClxgjoSHx3g28A+CcWwGcg1eQOxo16HdbJOiiNREGqTzbaWMtH2r8PV4SjOQ7LDhNvM65fOdcsnOul3OuF947zfHOuUgUYW7I78H7eJORMLNkvKHS7U0aZZWGxLsLuBbAzC7GS4T7mzTKhlsA/KB89ugQIN85tzfSQYmEW1QOjTr/yrNFKtbfAm2Bd8vn8+xyzo2P4nijQgNj/Qi43syygJPAvznnDkZxvI8As83sIbxhxrsi9B84zGwe3pBycvk7y6eBlgDOuVl47zBvBLYBRcD/jEScIn5TiTUREWnWonVoVEREpEkoEYqISLOmRCgiIs2aEqGIiDRrSoQiItKsKRFKo5jZSTNbV+2rl5mNNLP88h0WNpnZ0+XXVm/fbGYvRDp+EZGoXEcogXLMOXd59YbyLbGWO+fGmlkisM7MKuqWVrS3BjLM7D3n3BdNG7KISBX1CMVXzrlC4Cvgwhrtx4B1qIiziESYEqE0Vutqw6Lv1TxpZp3w6pV+XaM9Ca8m6LKmCVNEpG4aGpXGqjU0Wm6YmWUAZcCvy0uNjSxvzwT6lrfnNGGsIiK1KBGKX5Y758aeqt3M0oHPy98Rrmvq4EREKmhoVCLCObcVeB54LNKxRcUk/QAAAExJREFUiEjzpkQokTQLGG5maZEORESaL+0+ISIizZp6hCIi0qwpEYqISLOmRCgiIs2aEqGIiDRrSoQiItKsKRGKiEizpkQoIiLN2n8DyOxHiBH1eNUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# convenience function for plot legend\n",
    "def plot_proxy2():\n",
    "    p = plt.Rectangle((0, 0), 0, 0, color='none', \n",
    "                      edgecolor='red', hatch='xxx', alpha=0.8)\n",
    "    ax = plt.gca()\n",
    "    ax.add_patch(p)\n",
    "    return p\n",
    "\n",
    "xs = [float(x) / 100 for x in range(101)]\n",
    "roc_point_up = ss.zoh(*zip(*roc_bounds_point.upper))\n",
    "roc_point_lo = ss.zoh(*zip(*roc_bounds_point.lower))\n",
    "roc_interval_up = ss.zoh(*zip(*roc_bounds_interval.upper))\n",
    "roc_interval_lo = ss.zoh(*zip(*roc_bounds_interval.lower))\n",
    "fig = plt.figure(2)\n",
    "fix_plot_shape(fig)\n",
    "plt.plot(*zip(*roc_true), color='black', linewidth=2)\n",
    "plt.fill_between(xs, list(map(roc_point_lo, xs)), list(map(roc_point_up, xs)), color='blue', alpha=0.4)\n",
    "plt.fill_between(xs, list(map(roc_interval_lo, xs)), list(map(roc_interval_up, xs)), \n",
    "                 color='none', edgecolor='red', alpha=0.8, hatch='xxx')\n",
    "plt.plot(*zip(*roc_neg), color='black', linestyle='dashed')\n",
    "plt.plot(*zip(*roc_peek), color='red', linestyle='solid', linewidth=2)\n",
    "plot_proxy()\n",
    "plot_proxy2()\n",
    "plt.xlabel('FPR')\n",
    "plt.ylabel('TPR')\n",
    "plt.legend(['true curve', 'beta=0', 'after labeling',\n",
    "            'expected region via betahat', 'expected region via interval'],\n",
    "           loc=\"upper left\", bbox_to_anchor=(1,1))\n",
    "plt.title('Receiver Operating Characteristic curve')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Including the resulting ROC curve based on the peeked data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.2 PR curves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can directly use the contingency tables we already computed anyways\n",
    "pr_bounds = lambda tables: ss._lb_ub(lower=ss.pr_from_cts(tables.lower),\n",
    "                                      upper=ss.pr_from_cts(tables.upper))\n",
    "pr_bounds_point = pr_bounds(tables_point)\n",
    "pr_bounds_interval = pr_bounds(tables_interval)\n",
    "\n",
    "# an alternative without all these intermediate steps would be:\n",
    "# roc_bounds_point = ss.roc_bounds(labels, decision_values, beta=betahat)\n",
    "\n",
    "# compute the true ROC curve (we use Optunity's ROC function)\n",
    "_, pr_true = optunity.metrics.pr_auc(true_labels, decision_values, return_curve=True)\n",
    "# compute the ROC curve that would be obtained by treating all unlabeled as negative\n",
    "_, pr_neg = optunity.metrics.pr_auc(labels, decision_values, return_curve=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Normalizing weight for negatives: 15.815\n"
     ]
    }
   ],
   "source": [
    "# PR curve based on all known labels after labeling 300 unlabeled instances\n",
    "# we compute this by constructing the PR curve with weighted instances\n",
    "# the weighting is used to fix the empirical class contributions to our estimate of beta\n",
    "\n",
    "# compute the weight\n",
    "npos_prior = len(list(filter(lambda x: x == True, labels)))\n",
    "neg_weight = float(npos + npos_prior) / npos\n",
    "print(\"Normalizing weight for negatives: %1.3f\" % neg_weight)\n",
    "\n",
    "all_known_labels, all_known_dvs\n",
    "all_known_labels_sorted, _ = zip(*sorted(zip(all_known_labels, all_known_dvs),\n",
    "                                         key=op.itemgetter(1), reverse=True))\n",
    "known_after_labeling_cts = ss.compute_labeled_cts(all_known_labels_sorted)\n",
    "\n",
    "# weigh contingency tables\n",
    "weighted_cts = map(lambda ct: ss.ContingencyTable(TP=ct.TP, FN=ct.FN,\n",
    "                                               FP=ct.FP * neg_weight,\n",
    "                                               TN=ct.TN * neg_weight),\n",
    "                   known_after_labeling_cts)\n",
    "pr_after_labeling = ss.pr_from_cts(weighted_cts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAAEWCAYAAAD1t5d8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzsnXdYk9f7/98nCWGFDYKAslcAF4ijWrdg66porVZt/Tnrx2qrttrlt1bbaltt1bbWVamjzmq1aN1bq4giyh4KooCyh8wk5/fHSRCZQVnR87quXJJnnft5hLxzn3MPQikFh8PhcDgvK4KWNoDD4XA4nJaECyGHw+FwXmq4EHI4HA7npYYLIYfD4XBeargQcjgcDuelhgshh8PhcF5quBC+IBBCIgkhfes5pj0hpJAQImwms5ocQkgSIWSg8ucvCSHbW9omDoejWXAhbGKUH9TFSgF6SAjZQgiRNPY4lFJPSunZeo65RymVUErljT2+UoTKlfeZSwi5TAjp0djjcDgcTmPDhbB5GEYplQDoAqArgM+rHkAYmv7/sVt5n+YAzgDY28L2NDqEEFFL28DhcBoXTf/g1SgopQ8A/AvACwAIIWcJIV8TQi4BKALgSAgxIoRsJoSkEUIeEEKWVZ7KJIRMI4REE0IKCCFRhJAuyu2Vpwj9CCGhhJB8pRe6SrndnhBCVR/mhBBrQsghQkg2ISSBEDKt0jhfEkL2EEK2KseKJIT4qnmfMgA7ANgQQiwqXXMoIeRmJY+xQ6V97Qgh+wkhGYSQLELIz8rtToSQ08ptmYSQHYQQ42d5/oSQEcrx8wkhiYSQgKrPrtK9b6/yzKYQQu4BOE0IOUoImV3l2uGEkFHKn90JISeUzzWWEPLms9jL4XCaBy6EzQghpB2A1wCEVdo8EcB0AAYAkgH8AUAGwBlAZwCDAUxVnj8GwJcAJgEwBDAcQFYNQ60GsJpSagjACcCeWkzaCeA+AGsAowF8QwgZUGn/cAC7ABgDOATgZzXvU6y0MQtAjnJbFwC/A5gBwAzAegCHCCHaSqEPVt6/PQAb5bgAQAB8q7TRA0A75TNoEIQQPwBbAXykvJ9XASQ14BJ9lOP7A/gTwLhK15YCsANwmBCiD+CE8pg2yuN+JYR4NtRmDofTPHAhbB7+JoTkArgI4ByAbyrtC6KURiq9KFMAQwB8QCl9TCl9BOBHAG8pj50K4DtK6TXKSKCUJtcwXjkAZ0KIOaW0kFJ6peoBSlHuBWAhpbSEUnoTwCYwYVZxkVJ6RLmmuA1Ax3ru803lfRYDmAZgtPK+oHy/nlJ6lVIqp5T+AaAUQHcAfmBC95HyvksopRcBQHmPJyilpZTSDACrwESpoUwB8LvyWgpK6QNKaUwDzv9SaVsxgAMAOhFC7JT73gawn1JaCmAogCRK6RZKqYxSegPAX2BfNDgcTiuEC2HzMJJSakwptaOUzlJ+mKpIqfSzHQAtAGnK6cNcMM+pjXJ/OwCJaow3BYArgBhCyDVCyNAajrEGkE0pLai0LRnMG1ORXunnIgA6hBARIeRtZVBMISHk30rH7KGUGgOwBBABwKfKvc1X3Zfy3top7WgHILmSaFZACGlDCNmlnCbOB7AdbA2yoaj77Gqj4v9J+cwO48kXlLfApoIBdp/dqtzn2wCsnmNsDofThPCF/5ancvuPFDAvybwmUVDud6r3gpTGAxinDL4ZBWAfIcSsymGpAEwJIQaVxLA9gAdqXH8Hnnzw17Q/kxAyA8A1QsiflNI0pe1fU0q/rnq8Mrq0PSFEVMN9fwv2jDpQSrMIISOh5hRtFep6do8B6FV6X5NoVW3TshPA/xFCzgPQBQsOUo1zjlI66Bls5HA4LQD3CFsRSsE4DmAlIcSQECJQBouopgI3AVhACPFhQabEudL0XAWEkAmEEAtKqQJArnLzUykTlNIUAJcBfEsI0VEGrkxBHQLXwHuJAXAMwMfKTRsBzCSEdFPark8IeZ0QYgAgBEAagOXK7TqEkFeU5xkAKASQSwixAVvjexY2A5hMCBmgfK42hBB35b6bAN4ihGgpA4LUmcY8Aub9fQUWLatQbg8G4EoImai8nhYhpCshxOMZ7eZwOE0MF8LWxyQAYgBRYIEm+wC0BQBK6V4AX4MFYhQA+BtsXbEqAQAiCSGFYIEzb1FKS2o4bhxYcEoq2LrX/1FKTzTivXwPYDohpA2lNBRsnfBn5X0lAHgXAJRrkMPAAoTugQXwjFVeYwlY2kke2HTk/mcxhFIaAmAy2JprHtharepLxBdg3mKOcrw/1bheqdKWgZWPV3rXg8GmS1PBppdXANB+Frs5HE7TQ3hjXg6Hw+G8zHCPkMPhcDgvNU0mhISQ3wkhjwghEbXsJ4SQNYQlct9S5plxOBwOh9OsNKVHGAS2VlUbQwC4KF/TAaxrQls4HA6Hw6mRJhNCSul5ANl1HDICwFZlYvgVAMaEkLZNZQ+Hw+FwODXRknmENng6mfy+clta1QMJIdPBvEaYAT72zWEdh/MScx3IpJRa1H8kh6P5tKQQkhq21RjCSindAGADABBCaLfXXsOYMYcRFwfo6jauUQoFUFYGiETAzJlAW+6jcl5CCCE1le7jcF5IWlII74OVvVJhC5Z3VS+nTp1C//6pMDGxhpFRk9iG1FTgyhXgjTea5vocDofDaR20ZPrEIQCTlNGj3QHkKSur1ImxsTFKS0uxZ89HIDX5lI2EpSVw4wZQUFD/sRwOh8PRXJoyfWIngP8AuBFC7iv7uc0khMxUHnIEwB2wCiMbAcxS57q2trbQ0dFBSMifSE291SS2A4BQCFAK3LzZZENwOBwOpxXQZFOjlNJx9eynAP7X0Otqa2vj3XffxW+//Ybw8F1wd+9Q/0nPSJs2wIULQPfugJZWkw3D4WgE169fbyMSiTaBNZbmxTg4moQCQIRMJpvq4+PzqOpOjew+MXbsWKUQ7sGbb34N0kRzpDo6wMOHwI4dwMCBgKEhEBfHAmk6dWqSITmcVotIJNpkZWXlYWFhkSMQCHhtRo7GoFAoSEZGhjQ9PX0TWMPxp9BIIezduzcMDCyRlZWIlJQwtG/fdEVp2rcHHjwAfv0VIIRNl+rqAl5eTBA5nJcILy6CHE1EIBBQCwuLvPT0dK8a9ze3QY2BUChEp06BAICwsANNOhYhgIUFYG/PRNHeHiguBpKSmnRYDqc1IuAiyNFUlL+7NWqeRgohAEilrwEAoqMbs2tQ3ahmYPX0gNDQp/fl5wObNgHJPPuKw+FwNAqNFUJHx74QCrWQlHQNjx/XVcmt8TEzAyIigMeP2fv8fOD334H4eGDvXqC0tFnN4XBeCjIzM4XLly/n1W44jY7GCqFYrA9Hx16gVIGYmFPNOrYqtSIqCrh9G1i3jomhkxOQmwucPdus5nA4LwVZWVnCzZs3t6lpn0wma25znqK8vLxFx+c8HxorhAoF4OY2GAAQFXW82cc3Nmbe386dLGjGyoptt7EBzp1jlWk4HE7jMX/+fNuUlBRtd3d36YwZM2yDg4MNunXr5jps2DAHNzc3z9jYWLGLi4un6vjFixdbzps3zxoAIiMjtXv37u3i6enp4ePj4xYWFqZT9fp5eXmC0aNH27u6ukpdXV2lQUFBxgCgp6fXWXXMli1bTAIDA+0BIDAw0H7q1Km23bp1c505c2Y7Gxsb78zMTKHq2Pbt23ulpKSIUlNTRf7+/k5eXl4eXl5eHsePH9dvwsfEeQY0Nu5RoQDc3QcjOPgTREUdB6W0ydIoasLIiEWPisVPbxcKAW1tICQEGDmy2czhcJoVQohPU1yXUnq9tn0rV668P3ToUN2YmJgoAAgODja4deuWflhYWKS7u3tZbGysuLZzp06dardhw4Zkb2/v0tOnT+u/99577a9cuRJX+ZhFixa1NTQ0lMfFxUUBQEZGhrDmqz0hMTFR59KlS3EikQiTJ0/Gjh07jOfOnZt1+vRpfVtb27J27drJhg0b5jBv3ryH/v7+hfHx8WJ/f3+XO3fuRKr/VDhNjcYKIaWAjU0nSCTmyM6+h4cP42Bl5dasNlQVQRUWFkBYGDBkCBNFDofTNHTo0OGxu7t7WV3H5OXlCcLCwiRjxoxxUm0rKyur9q35/Pnzhrt27bqjem9hYSGvb/xRo0bliJR5VOPHj8/+6quvrOfOnZu1Y8cO08DAwGwAuHTpkmF8fHxFe4DCwkJhTk6OwMTERKHWTXKaHI0VQoUCEAoF8PAYhGvXduLmzQMICFjU0mYBYFOl5eVAYiIglba0NRxO41OX59ac6OnpVYiJSCSiCsUTbSkpKREAgFwuh4GBgUzlSdZGbbNKlbcVFxc/dYBEIqkYcMCAAY+nTJminZqaKjp69Kjx119/naq6bmhoaLREIuGpJ60UjV0jpJSlM3ToMBQAcODAJ/jll+FITW0dMw6Ghqx7BYfDaRyMjIzkjx8/rvUzy9bWVpadnS1KT08XFhcXk2PHjhkBgKmpqcLW1rbs999/NwEAhUKB//77r1oDt759++avWrWqIhhHNTVqZmZWfuPGDR25XI6DBw+a1Da+QCDAkCFDcmfNmtXO2dm52MrKSg4AvXr1yl+xYkXFdS9fvtzIzeM4z4vGCqFCwYTQ1/ctDBu2BNra+rh16x989VUHnD79a0ubB2Nj5hHm5ra0JRzOi4GVlZXcx8en0MXFxXPGjBm2Vfdra2vT+fPnp/n5+XkMGDDA2dnZuUS1b+fOnXe2bNli7ubmJnVxcfH866+/jKue/+2336bl5uYKXVxcPN3c3KRHjhwxAIAlS5Y8GDFihHOPHj3cLC0t6wwPffvtt7MPHjxoOnr06BzVtg0bNqTcuHFD39XVVerk5OT5888/8xSQVgZhta81B19fXxoaGoqvv2ZioyqGnZ//EMHBS3Du3DoAQEDAOvj7z4SeXsvZeu8emxrt2RNo144F0nA4mgAh5Dql1LfytvDw8KSOHTtmtpRNHM7zEh4ebt6xY0f7qts1do1QNTWqwtDQEuPH/wqx2AMnTszB0aPvQU9PBDe3qTA0BExNm99GGxsgIQGIjARsbYF33kGLCjOHw+FwqqPRU6M1beva9X0sX74KAHDgwHQUFS2FXF6OwsJmNhDMA7SyAuzsgPR0YPt2oKSk/vM4HA6H03xorBBW9QgBth7n6gosXPghvv/+e1BKsXr1YgQF+SEmpmU77NrYAPfvA7t21SziHA6Hw2kZNFYIaxKTggLAR5nmu2DBApw8eRL29vZISLiJzZu74vTptc1rZBVsbYHYWCA6ukXN4HA4HE4lNFYIq3qEqihSB4cn2wYMGIDbt29j1qz3oVDIsHfvh0hNrTOVqMmxsACCg4GyOlOAORwOh9NcaKwQqoSPUuYJJiezaVH9KlX8JBIJfvllDfr0mQmFQo49ez5AS0bKSiSsQHdoKLO9qAho4XrBHA6H81Kj8VGj9+8D5ubA4MFAx461H79w4VJcvboL0dEnEB5+CJ06jWg+Y6vQti1w8CDw77/sPgCWgC8SsWAauZz9LJEALi6sGbCREWBgwISeECAzEzh9mr0fOJCXcuO8+MTGxoqHDh3qEh8fr1bVjDVr1pgNHz48397e/rlbQ3zyySdWO3bsMBcIBFi5cuW9wMDA/Oe9Jqf1oJFCSOkTAaEUGD6ciUVddOxojt69v8KJE3Owd+88eHr6Q0urWgH6ZkEsZvYSAggE7B7KypiXa2TEtisUrEzblSvAxYtPvF8tLSb8qansOuXlbN1x2DCgTRsmqKpj4+OB8+eBzp2BLl2qBxdxOC8y27dvN+/UqVPx8wrh9evXdfbv328aGxsbmZycrDVo0CDXESNGRKhqjHI0H438n6y6PiiR1H+OpSXQvft7iIhYj7S0SJw8+SOGDPmk6Yysh8rJ9YTU7NGJxdWnemUyoLiYBd6orpGTA/zxx5PrGhqynzMzmRe5dy/LZRwyhK1Rcjiaikwmw6hRo+wjIiL0HB0dS/bu3Zt08+ZNnXnz5rUrKioSmJiYyHbs2JF0+vRpSUREhN6kSZMcdXR0FKGhodFffvml1dGjR41LS0sFvr6+hTt27EgWCOpfHdq3b5/xqFGjsnV1dam7u3uZnZ1d6dmzZ/UHDhz4uBlumdMMaKwQVv65qljUhFAIuLmJ8Nprq7F580AcObIMtrYd4O39etMZ2gSIREzcKmNiwl4Am1YtL2f/qgKHzMyApCTgp58AT09g0CAuiJznx8/Pr1q7l1GjRmUvWrQoo6CgQDBgwACXqvsnTJiQOWfOnKy0tDTRiBEjnCrvCwkJia1vzKSkJJ3169cnDR48+PGYMWPsv/vuO4vg4GCTw4cPJ1hbW8s2btxosmDBApu9e/cmrVu3rs0PP/yQ8uqrrxYBwEcfffTohx9+SAOAkSNHOuzatcto/PjxeV988YXl3r17zaqO1b1794KgoKCUBw8eiLt3716RiWxtbV2WkpIiBsCF8AVBo4VQoWBTizpqznA6OQHR0QPwyitTcOnSZvzyyzAMH74MQ4Z80qy9DJsSobB6KTdCWGK/QsHqn965A0yfzqZSORxNwsrKqmzw4MGPAWDixIlZy5cvbxsfH6/bv39/V4AV1LawsKhxKvTff/81WLVqlVVJSYkgNzdXJJVKiwHkLV269OHSpUsf1jZmTcF1hBDNqk3JqRONFsLSUlZvVF0Ns7Jix06YsAFmZvb455/FOHjwM6SkhOGdd7ZAR0eNOVYNRiBgU8TZ2cDmzcC0aWy9kcN5Fury4AwMDBR17W/btq1MHQ+wKlW/sOrr68udnZ2Lb968GVPXeUVFRWT+/Pl2V69ejXJ2di6fN2+etapNU30eoa2trcoDBACkpqaKbW1tnzsAh9N60Mj0CZUQlpU1rIaoubnKixTg9dc/x3vvHYSOjiFu3NiHFSt6ID29wX+XGompKXuGq1cDx48DeXktbRGHox5paWnikydP6gPAn3/+aern5/c4OztbpNpWWlpKQkNDdQBAIpHI8/LyhABQVFQkAAArKytZXl6e4J9//qlop7R06dKHMTExUVVfQUFBKQAQGBiYu3//ftPi4mISExMjTkpK0unbty+fFn2B0GghLC1tmBBKJCyQpLSUve/YcRg++eQqLC3dkJoagW++8cW1a7sb3+BWiLk585AvXgRWrAB+/BE4dozVROVwWiuOjo4lv//+u5mrq6s0JydHtGjRoke7du1KXLRoka2bm5vU09NTeu7cOQkATJo0KfP999+3c3d3l+ro6CjefvvtDKlU6jlkyBDnjh07qi1kvr6+JSNHjsx2dXX1DAgIcF21alUyjxh9sdDINkwXL4Zi2TI2zRkQALzyivrn797N1skqTwmWlBRg27ZpCA1lItinzyyMGbMKWlovR3KeKrE/N5cF2tjZsXQMa+uWtozTUvA2TJwXkdraMGm0R0gpy7trCA4O7EO/Mjo6Bpg6dSfGjfsFIpEY5879im++8cWpU6uRnZ3SOEa3Yghhkbc2Niy/MSsL+PVXICTk6QhdDofDeRHRSP++ch6hOqkTlbG0rHk7IQR9+86Cvb0fNmwYg9TUCOzZ8wH27PkA9vZ+6NIlEB06DIOVlfsLE2FaG+bmbPp4/34gJYV5h2Jx/edxOByOJqKxQqgSw4YKYX35c/b2vvi//4tAePghhIX9hdu3jyApKQRJSSHYv38hDA2t4O7eH25u/eDm1h/m5g4vpDBqazPvMCwMePgQGDOG5x5yOJwXE40VQlUZMXWqylRGT48FzJSU1J5/qK2tDz+/cfDzG4eysiJERBxFWNh+xMScRH5+OkJC/kRIyJ8AACOjtnB07AEHh+5wdOwOOzsfiMUNa0OvUCjw+HEW8vMfIj//IUpLCyCRmMPQ0AqGhpbQ1pa0iNgKBED79sCjRywZ384O6N6d1T/V1W12czgcDqdJ0FghVCXTP8sHsoMDEBenXiK+WKyHLl1GoUuXUaCUIi0tGrGxZxAbexpxcWeRl5eGsLD9CAvbDwAQCEQwMbGFSCSGUCiu+FcoFIFSBRQKOeRyGSiVQy4vR2FhJgoKMkBp7d16xWI9GBpawcrKHc7OveDs3Av29l2brVZqmzbsmefmsmAjoRDw9QX69ate5YbD4XA0DY0VwpISFtX4LI6Stzdw/TrzDBvStYEQAmtrKaytpejX73+glOLhwzjcvXsFd+5cwd27V3D//i1kZSU12CZ9fVMYGlrCwKANtLUNlAL5EHl5aSgrK0Jm5h1kZt5BRMQRAIBIJIadnS8cHLrDxsYb1tZesLaWNtgbVRdCWPECY2NWvu36deDmTbZ+6OLScM+cw+FwWgsaK4SUPqmv2VDc3YGxY5l307bts0/zEUJgZeUGKys39OjxDgCgpKQQBQWPIJeXQy4vg0zGXnJ5OQQCYZWXCPr6ZjAwsIBIVHM0CqUUJSUFyMtLxb17YUhIuIjExIt48OA2EhMvIzHx8lP2mJs7wsGhO7p0CYSnZwDE4safwxQKWYRpcTGwbx/bZmzMRNHdvdGH43Dq5ffffzdZtmyZtYWFRfmPP/6YkpKSIh47duxzlYrQ09PrXFRUFFbb/oa2hQKAwMBA+6FDh+ZNnjw5Z+zYsXYff/zxQx8fn5LnsZPz/DSpEBJCAgCsBiAEsIlSurzK/vYA/gBgrDxmEaX0SH3XVQlhQ5Lpq9K5M/Nydu0CHB3VP4/SJ/0Ca0JHR9KopdoIIdDVNYSuriGsrNzh5zcOAFBUlIvExMtISQlDamoEHjyIQHp6DDIyEpGRkYiQkB3Q1taHt/dQdOkyGp6eAY1eQk5Xl60bAkBhIRAUxDpc9O7Npq05nOZiy5Yt5qtXr743bNiwgjVr1piFhobqN0QIy8vLoaWl1ZQmVmP37t3JzTogp1aaTAgJIUIAvwAYBOA+gGuEkEOU0qhKh30OYA+ldB0hRArgCAD7+q6tym173jqZHTuyyiq5ucyjUYcHD1hpN2Pj5xPi50VPzxje3q/B2/u1im0yWRnS02MQGXkUN27sQ1LSNYSG7kZo6G4IBCI4OvaAh8dAuLsPgIODH4TCxvvDl0jYmuvRo8DZsyzIpl07tr5oZcULfHMah4EDBzqlpaWJS0tLBTNnzny4YMGCzAULFrS9fv265P3337f7+++/844cOWJcUlIicHd3l8yfPz9t7NixeVOmTGkfHR2tK5fLyWeffZY6YcKE3DVr1pj9+++/RqWlpYKioiLBlStX4moaMy8vTxAQEOCcl5cnlMlkZPHixakTJkzIBWpuC2VgYKC4cOGCXtXWUHZ2dk/VJ/Xz83NTdcfQ09PrPGXKlEfHjx830tHRUQQHBye0a9dOFhkZqT1+/HgHuVxOBg4cmLdhwwbLurxUzrPRlB6hH4AESukdACCE7AIwAkBlIaQAlN3zYAQgVZ0LU8qm5xqaTF8VQlh39z/+qC6Ecjmb+isqYvvEYiaAhABTp7IO86mprav6ikgkhq1tB9jadoC//8fIzExCWNhfuHHjL9y9exUJCReQkHAB//zzf9DWlsDVtW+FMFpbez53ZKpIxAKRZDIgIwNITmaVagQCYMIEwK1a0x6OxkKIT5Ncl9Lrde3esWNHkqWlpbywsJB07txZOmHChJwffvgh7fz584YqUVmzZk1RaGio/tatW+8BwOzZs2369euXv3fv3qTMzEyhr6+vx/Dhw/MB4MaNG5Jbt25FWlpaymsbU09PT3H48OEEU1NTRVpamqhbt27u48ePzwWqt4X6/vvvLT777LNHc+bMaV9Ta6jaxiguLhb06NGjcO3atQ9mzpxpu3btWovvvvsubfbs2e1mzZr1aMaMGdnfffcdT2BqIppSCG0AVC7Lch9AtyrHfAngOCHkfQD6AAbWdCFCyHQA0wGgffv2UChqblr7LLi4MM+ysPBJwEdxMRM5Ozvm1Vy/ztYSHz0CXn8dcHYG3nuP1eesKw2jpTE3t8egQfMxaNB8FBXlIi7uLKKjTyEm5iTS02Nw+3Ywbt8OBgAYGlqhQ4eh6NHjXTg59XwuURSJngTWAOzLxB9/AOPHA15ejXFnnJeVFStWWB4+fNgYANLT07UiIyN1rKys6qwbevbsWcNjx44Zr1mzxgpghbkTEhLEANC7d+/8ukQQABQKBfnggw9sr1y5IhEIBHj06JH4/v37IqB6W6g1a9a0uXXrVp66raFUaGlp0bfeeisPAHx8fB6fPHnSEADCwsIkx48fTwCAqVOnZn355Ze29T8lTkNpSiGs6ZO0asGucQCCKKUrCSE9AGwjhHjRKrkElNINADYArNYowD5sG0OAhEJgwABgxw4megIBkJYGTJrEmtgCgIcH229qytIGAJaP2Lcv8O+/bBqwtaOnZ4xOnUaiU6eRAICcnAeIiWGiGB19Enl5abh4cRMuXtyENm1c0KPHO+jefRJMTds1wtisos/OncD//te6vGjOM1KP59YUBAcHG5w7d84gNDQ0xsDAQOHn5+dWXFxc72o0pRT79u1L6NixY2nl7RcvXtTX09OrPW9Jyfr1602zsrJEt2/fjtbW1qY2NjbeqnGrfmEkhIBSStRpDVUZkUhEBcqFdZFIBJlM9uJV6WjFNGVIw30AlT9FbVF96nMKgD0AQCn9D4AOgHpX/hQKNkXZWAEZXl4s4jE/H7h/Hxg58okIqvZPmwa8+ebTpcY6dmR2yGSNY0dDUU0/pqdXr59aHyYmNujRYxImT96KFSseYPHiW/D3Xwgjo7Z49CgeBw9+jk8/tcPata/j9u0jUCjq/byoE11dlnO4Zw+bYuZwGkpubq7QyMhIbmBgoAgLC9MJDw+vcU7I0NBQXlhYWPHp0K9fv/yVK1daqn6HL1261KBQ6ry8PKG5uXm5trY2/eeffwxSU1MrPgWqtoXq2bNnYYcOHUpqaw3VUDp16lQYFBRkAgC///57C0YlvNg0pRBeA+BCCHEghIgBvAXgUJVj7gEYAACEEA8wIcxQd4DGKrYiEgGvvgp8/DEwezarnlIVR8cnEZIqJBKgWzc2ZdrYFBYCd+8+3StQLmdTsfn5bF+bNky0u3Zl2zLUfnJPQwiBjY03Ro1ajm+/vYf33/8Xvr5jIRSKERFxBD///DoWL3bFiRMr8fhx9jPfk6kQkcY2AAAgAElEQVQps/H06We+BOclJjAwME8mkxFXV1fpp59+al1bK6UhQ4YUxMXF6bq7u0s3btxosnz58lSZTEbc3d2lLi4unp9//rlNQ8adOnVqdnh4uL6Xl5fH9u3bTR0cHCrSHaq2hVqwYEGGjo4Ora01VENZu3Ztytq1ay29vb090tLStCQSSZ3TuJxno0nbMBFCXgPwE1hqxO+U0q8JIV8BCKWUHlJGim4EIAGbNv2YUnq8rmv6+vrS4OBQrFvHRKul61+qyo+pRNnEhCXqq0uJ8k+q8jRvWRmLTn39deDGDebxCQSAlhZbF9XTYwLcufMTrzgrC9i2jU3rikRsytfa+vm85sLCTFy8uBnnz69DVhaL9BaL9dCnz3sYNGgBjIysGnxNuZx5sZMn8+CZ1gxvw9Q6KCgoEOjr6ysEAgE2bNhgsnv3btNTp04ltrRdmkptbZg0sh+hSgjff//5Uygag4QE5sGVlAChoUyMDA3rT/jPywMKCpi45eezbaqfR45knimlQGbmEwGsi6IiIEoZk3v/PvDff2z9sryceWKEMC+yoWurCoUct28fwZkzaxEdfQIAoKWlg969p2Pw4I9hYtKgL9goKgKys4FZs1hqBaf1wYWwdXD06FHJ3Llz21NKYWhoKA8KCkry8vIqrf9MTk28UEJ46FAo1q8H5swBzMxa2qLqpKSw4BqZrPb8ufx8Jp7TprEqLYWFwL17zAM0Nwf8/Z9v6pdSliMZHMzSTF57DXj8GDh1inmcNjbMa2wo9+7dwOHDy3Dz5gEAgEikDX//jxEQsKhB5d1yctj9zZzZMA+a0zxwIeS8iLxQQnjwYCg2bADmzm3ZpPa6yMlhKQMqT6yyqFHKvLtJkwDbJg6GTk1lz0jlBZaUsIT3c+dY8IqZ2bMJ7v37t3DkyDJcv74XAGBmZocxY35Ep04j1U69ePiQieCUKbxWaWuDCyHnRaQ2IdTYWqNA6y7jZWICzJjBpklVbaPEYvaSSFgUZXN0VqqaqqCjAwQEsEjYY8eAxERmk4lJw2qu2tp2wPTpe5CQcAm7ds1GSspN/PbbKEilg/HGG8vRvn3neq9hacmeT1AQ8O67XAw5HE7LoLFCWNXLao3o6jasjmlzYmvLPLG0NNZFIjLySfSrsbH6VXucnV/Bp5+G4vz59Th48DNERR1HVNRxuLn1x6BB8+HpGQBBHd9Y2rZlXuvWrUwM61sH5XA4nMZGY4UQaP1CqAm0bcteQ4awmquJiUBICEvPUHfqVCAQom/fWfDxGYOjR5fj4sWNiI09jdjY02jbVopBg+aje/dJEApr/nWztmZRstu3s+niysE8qi89HA6H01RotBC25qlRTcTYGPDxYa/kZJbvl5jIhEguZ5GrxsYsjaMmDAwsMGbMSgwduhgXLmzEqVM/IS0tClu3TsHJkz/irbfWwM2tX43n2tiwSNevv2ZCqKXFSt0JBMDo0ay6D6d18e23sEhNRaNVbre2Rvknn6ifR9xcxMbGis+cOSOZOXNmg5JoK7dcamybOnfu7B4WFqZ25ZqG0tAWUZULiKtz/OXLl3XVaZWl6iSiqhurDl999VWbDz/8MNPAwEDtKiAaLYTcU2g67OxYrl95OQv8efQIiIkBoqNZwI1QyMRKS4uJY+UvJbq6Rhg8eAH695+D0NDdOHRoMVJTI7BqVX906RKIwMAfYG5uX21MW1tWNUguZ/8aGAClpWwNcdAgVvRAXHPbRk4LkJoKLTs7NFqdoORktMr/3fj4eO3du3ebNlQI1eFZ2z81pQgCTd8iKjQ0VK+hrbLUZf369ZbTpk3LbogQaqRPxYWw+dDSYikgXl7MM/v0UxYENGgQ0KULC3hJTn66Ao4KkUiM7t0nYsmSaIwYsQxisR5u3PgLX37pge3bZyA29gwUiqcLZagKB2hrM7HV0wPs7YEzZ4AffgCuXWPiyHk5+fXXX029vb093N3dpePHj7eTyWSIi4sT29nZeaWlpYnkcjl8fHzc9u/fbxgbGyt2cHDwHDVqlL2rq6s0ICDAsaCgQAAAFy5c0Ovataubp6enR69evVySk5O1ACAiIkK7Z8+erm5ublKpVOoRGRmp/dlnn9mEhoZK3N3dpUuWLGkjk8kwY8YMWy8vLw9XV1fp999/bw6w4tqTJk1q7+Tk5Nm3b1/nzMzMGh0NPz8/t9mzZ9t07drVbdmyZZapqakif39/Jy8vLw8vLy+P48eP6wNAamqqqGfPni5SqdRj/PjxdtbW1t5paWkigDUNVo05Y8YMWxcXF09XV1fpxo0bTQBWl9XPz88tICDA0cHBwXP48OEOVcsk3rhxQ8fb27tiriU2Nlbs6uoqVdl4/vx5PQB4++2323t5eXk4Ozt7fvjhh7VWCg4KCjLr3Lmzu4uLi+eZM2f0ACA/P18wZswYey8vLw8PDw/p9u3bjUtKSsi3335r/c8//5ioqv+cOXNGr3Pnzu4eHh7Szp07u4eHh2urrpuenq7Vu3dvFzs7O6+ZM2dWxNnXZNeyZcvaPHr0SKtPnz6u3bp1c1Xz10qzPUI+Ndr8CIUsSb9yofG7d4EDB2pvS6WlpYPXXvsMPXq8g/37FyIk5E9cuLABFy5sgKGhFXx8xsDXdywcHXvUGFgjFDIPtbgY+Ptv4PBhVmzAyYn9Dhgbt858Uk7jcuPGDZ19+/aZhoaGxmhra9MJEya0/+2338xmz56dNXfu3PTJkye379q162M3N7eSUaNG5cfGxoob2iZp/PjxDgsWLEifNGlSblFREZHL5eTrr79+sHLlSsszZ84kAMAPP/xgbmRkJI+IiIguLi4mXbt2dR82bFj+1atX9RISErRjY2Mj79+/r+Xt7e357rvvZtV0L7m5ucJr167FAsCwYcMc5s2b99Df378wPj5e7O/v73Lnzp3IRYsWWffp06fg22+/Td+3b5/hzp07q5UP2bp1q/Ht27d1o6OjI9PS0kR+fn4egwcPLgSA6Oho3Zs3b96xt7cv9/HxcT9x4oTE39+/UHVuly5dSsrLy0lUVJRYKpWWbd261XTkyJHVpnFXrVr1wNLSUi6TydCzZ0+3q1ev6nbr1q246nFFRUWCsLCwmH///Vcyffp0h/j4+MhPP/20bU0tsD755JPUylOe2dnZgpCQkBgtLS38/fffBh9//LHtsWPHEgEgKipKLzw8PEpXV1fh7OzstWDBgofOzs7lNdn1+eefP1q3bp3luXPn4tq2bat2FWiNFkLuEbYOHByA6dOBTZvYFGptRQRMTGwxZcoODBnyKUJCdiI0dBcyMhJx5sxanDmzFmZm9ujTZxZ69ZoKff3qZXl0dZkglpezqjkXL7LthABjxgAdOjThTXJanKNHjxpERETodezY0QMASkpKBG3atJEBwLx58zL3799vEhQUZHHr1q2KnqcNaZOUk5MjePjwoXjSpEm5AKCnp0dRvWMOTp48aRgTE6N36NAhEwAoKCgQRkVF6Zw7d87gzTffzBaJRLC3ty/v0aNHQW33Mm7cuIpp1kuXLhnGx8dXJC8VFhYKc3JyBCEhIZK///47AQBGjx6db2hoWK3O6IULFyrGbNeunaxbt26FFy9e1DMyMlJ4e3s/dnJyKgcAT0/PosTExGpTzyNHjszevn276TfffJN+4MABk927d9+peswff/xhGhQUZC6TyUhGRoZWeHi4Tk1COH78+GwAGDJkSGFhYaEgMzNTWFcLrMpkZ2cLx44d65CUlKRDCKHl5eUVn+69evXKNzMzkwOAs7NzSWJiorazs3O5unapAxdCTqMgkQDvvAOsX88q6xgZsW01ee3W1p4YOXIZRoxYiuTk6wgN3Y3Q0N3IykrC/v0fIzj4S3TvPhH9+s2BtbW02vlaWiy4RkVJCavkc+8e4O3NhLi4mEXBtmnD8xNfFCilZMyYMVm//PLLg6r7CgoKBOnp6WIAyM/PF5qYmCiAhrVJys7OVmuOiVJKVq5ceS8wMDC/8vbg4GAjdYtJVF6/opQiNDQ0WiKRPCW66hQ7qesYbW3tip1CobDG1k4TJ07MGTNmjONbb72VQwiBt7f3UwsPMTEx4p9//tny+vXr0RYWFvLAwED7kpKSGp9TLc+61hZYld8vXLjQpk+fPgUnTpxIjI2NFffv37+iErFYLK58H7S8vJw0xC510MjJRU3JI3zZMDFhnuGgQSzCNCnpyZeWmiCEwN7eF6NHf49vvknC7NnBkEoHo6ysCOfPr8eSJZ746adBuHUruM42UDo6bB0xJATYsAFYtgxYtQrYvBlYsYJNpzZFhxBO8xIQEJAfHBxs8uDBAxEAPHz4UBgXFycGWBf60aNHZ3366aep7777bkWfmIa0STI1NVVYWVmVbdu2zRgAiouLSUFBgcDIyEheWFhYUZBw0KBBeevWrbMoLS0lAHDr1i3t/Px8QZ8+fQr27t1rKpPJkJycrHXlyhUDde6rV69e+StWrKiYR7l8+bIuAPj5+RVu27bNFAD2799vmJ+fX60oYp8+fQr27dtnKpPJkJqaKgoJCZH07t27zkbFlfH09CwVCARYvHix9RtvvFEtGCgnJ0eoq6urMDU1laekpIjOnj1ba4bxzp07TQDg2LFjEgMDA7mZmZm8thZYVVtl5efnC21tbcsAYP369fVWkK7LLn19fXleXl6DtE2jPUK+Rtj6MDVl0Z29e7Mk+aQklqdYHwKBAN7er8Pb+3Wkpkbh7Nmf8d9/fyA6mjUOtrBwQr9+76Nnz8nQ1a1enFQorLlcnVwOhIczkfT1BXr2ZAE+Vb9EUcoKoOfmskAdExMeoVof1tYob8xIT2tr1NnF3cfHp+Tzzz9/MGDAAFeFQgEtLS26Zs2ae/Hx8eKbN2/qb968OUYkEuHAgQMmq1evNgsICChQtUmaNWuWnYODQ2nlNklz5sxpX1BQIJTL5eS999576OvrW7J9+/a706ZNs1u6dKm1lpYW3bt3b6Kfn1+xSCSibm5u0vHjx2d+/vnnj5KSkrS9vb09KKXE1NS0/MiRI4kTJ07MPXXqlKGbm5ung4NDiZ+fX61To5XZsGFDytSpU9u7urpK5XI56datW0HPnj3vLV++PHX06NGOUqnUpEePHoUWFhblxsbGT02PTpw4Mffy5csSDw8PT0IIXbJkyf327dvLbt26pfZzHzVqVPbSpUttV6xYUc3T7tGjR7GXl1eRi4uLZ/v27Ut9fHwKa7oGAJiYmMg7d+7sXlhYKNywYcNdAFi+fHnq9OnT27u7u0sppcTW1rb0zJkzCUOGDCn44Ycf2rq7u0vnz5+ftnDhwvSpU6c6rFmzxqp37975tY2hjl3vvPNO5pAhQ1zatGlTfvXq1Th1noFG1hrdsSMUu3cDn3/OxbA1k5vLWlRVrnXaEB4/zsGlS5tx9uzPFW2gtLUl6N59Ijp1egOurn0gEqn3OaxQsHZW5eVsqtTVlQXZiMUsVzI5me1TCaRCwQJyAgJqz5t8kXkRao3GxsaKhw4d6hIfHx/Z0rY8C8XFxUQkElEtLS2cPHlSf/bs2XYxMTFR9Z/JqY0XstYonxpt3RgbA8OHA7t3syjThgqKvr4JBg9egIEDP0R4+CGcPr0GcXFnce7cOpw7tw46Oobw9PRHhw7D4eU1BBJJ7aGjAsGTiNaSEpYTWV7OPEaJhPW1FFX6a1AoWEBOaipLG+FRqZzmJiEhQfzmm286qbzf9evXJ7W0TS8qGukRbt0ain37gMWLW9oaTn2o2kEdPswCV0Qi1gaK0iddOBoyBXn//i1cu7YLt279g9TUiIrthAjg6NgDnp7+kEr9YWfnA4HgGfpMVeHhQyacZmaAVMpqx1pbv/gBOC+CR8jhVOWF8gjriJvgtDIIYeuFVlbAvn3s/87Y+In3df8+ExotLXZMfV6+rW0H2Np2wBtvfIPMzLu4desfhIcfQlzcOSQmXkJi4iUcOrQY+vqm8PAYBKnUH1Lp4AY3D1Zhacn+LSoCrlwBLlxgNjo6sqlTZ2e+lsjhaDpcCDnNgosL8Mkn1bdTyno2Hj0KREWxgBd1hcXc3AH9+89B//5zUFych5iY08ruF8eQmXm3Ii0DYCkbUqk/PD394ezcG2JxA3pOgXmuqs4YCgWLQt2+ndnavTurhdq2LRdFDkcT4ULIaVEIYVOmEyawyM6jR9nanYEBIJOxlzrBNrq6Rujc+Q107vwGKKV49CgBUVHHEBl5DHFxZ5CaGonU1EicPLkKWlo6cHHpA09Pf7i59YONjXeDplEFAmaTqenTyf0CAeDuzqJTHRxeziAbDkcT0UghlFerr8DRdAQC5ll16sQ8w8hI1r1eLAZu3GBrdappU5GI5Snq6jKxqTqdSgiBpaULLC1d0K/fbJSXl+LOncuIjDyGqKjjSEkJQ1TUMURFHQPAIlEdHbvDyekVODn1hIND9xpTNGqicnK/XM7KzUVEsDXEwYOBjh25IHI4rR2NFELuEb646OiwYt5dujzZNmgQq1ZTVMQ8sIwMlu7w6BHbRikTG0tLlk9YFS0tbbi59YOLSz+89tpy5OU9RELCCcTFHUdCwkVkZt6tyFcEWOCNjY037O390L59Z7Rv3wU2Nh3qnU4VCln0qYUFq2xz4ADwzz9sTbRtWyaMpqaN+bQ4zQFvw1Q/H3zwgXXfvn0LRo4cWWvuZHBwsIG2trZi0KBBaif8PwvBwcEGlWvDqoNGCqFczlMnXiZEIjbVWBOq/MCrV5nnqAqC1tVlL5Uw5uSwaFVTU8DQ0BJOThPQrt0ETJgAFBWlITHxMhISLuHOncu4d+8G7t8Px/374RXjECJA27YeaNeOCaOTU0/Y2fnWOqWqqosql7NxY2OBuDhg/HgWaMPzXzUH3oapfn766afU+o45ffq0gUQikTdECJ/1+TQUjfxz5B4hR4UqP/CNN1iBhenTWe6ihwerEKPqbejrC3z4IbBgATB3LvDxx8DAgcyrzMtrC3f3QLz55iosWnQFP/2UhwULzmPMmB/RvfskWFt7AQBSUyNx9ep27N07D8uXd8dHH1lh8+YJuHp1BwoLa84qEAqZKLZty6ZLN20CvviC2Xr0KFsD5agPb8PUOtswBQYG2m/ZssUEAGxsbLw//PBDa6lU6uHq6ioNCwvTiY2NFW/dutXit99+s3R3d5cePXpUUtt9z5s3z3rcuHF2r7zyisuoUaMcOnTo4B4aGloRJeDn5+d24cIFvbpaNzUUjfQIuRByakJb+0mLKD+/uo/V0QH69gW6dgVu3wbOnwcyM9lMg0KhCxOT3nBw6F2R5lFWVowHD27j3r0buHfvOqKjTyArKxkhITsQErJDWTfVD97eQ9Gp00hYW3tWK0IskTzJP5TL2ZgpKYC/PxNLQ0N2D5ya4W2YWncbpsqYm5vLoqKiopcvX26xfPlyy927dydPmjQpQyKRyL/66quHdd03ANy6dUvv6tWrMRKJhC5ZsqTNjh07TH19fVOTk5O1Hj16pNW7d++iulo3NRQuhJyXGn19FqTTvTtbb8zOZlOtcXFsOlMmY+Kkra0La2s/ODgwhaWUIj09BhER/yIy8l/Ex5/H3btXcffuVRw69AUsLJzQseMIdOo0Ek5OPatNoQqFrFB4ejqwcSPbZmQETJ7Mq9jUBm/D1LrbMFVm/PjxOQDg5+dXpHpOVantvgEgICAgV9WNY9KkSTkDBw50/fHHH1O3bt1qMmzYsByg7tZNDUUjhZBHjXKaAlWuoK0tm0otLQXu3GGv7GzW5unxY+bVlZQQKBQe8PLyQL9+8yCTFSIm5hRu3TqE8PBDyMhIxMmTq3Dy5CoYGFigS5fR6NNnFmxsvJ4aU5WwD7AgoI0bgSlTWMAN52l4G6Yabal1X3O3YaqMjo4OBQCRSERrGldle033DQD6+voVz8fBwaHc2NhYdvXqVd39+/ebrl+/Phmou3VTQ9HINUIeLMNpDrS12Vrj668DEycCCxeyYJd27YD+/YFJk4DOnYG0NCAjQwKpdAQmTdqM779Px0cfXcCgQfNhYeGEgoIMnDu3Dl995Y2VK/vi+vV9kMurN1qwsGDBPuvWMW+U8zS8DVPrbsNUHwYGBvKCgoKKe6jtvmti9OjR2d98841VQUGB0M/PrxhoeOumutBIIZTJuBBymh+RCPDyYqLYty/rYDF8OBPIAQPYGmNyMlBcLISzcy8MG/YDFi2KxxdfhKNv3/9BW1uCuLhz2LBhDD791B6HDy9DcXHeU2OYmbFiAlu2AEeOAIW1Nr15+ajchsnV1VXav39/15SUFK3Dhw9Lbt68qb9s2bL09957L1tLS4uuXr3aDABUbZhcXV2lOTk5osptmBYtWmTr5uYm9fT0lJ47d04CANu3b7/7yy+/tHF1dZX6+vq6p6SkiCq3YVqyZEmbDz/8MNPd3b3E29vbw8XFxXPatGl25eXlZOLEibmOjo6lbm5unlOmTGnfkDZMN27c0Hd1dZU6OTl5/vzzzxYAa2F0+vRpQ6lU6nH48GGj2toweXp6Fnt4eHj27dvXVdWGqSHPddSoUdkHDx40nThxYrX1wcrtjiZOnGhfVxum+ggMDMw9fPiwsSpYprb7rokJEybkHD582HTEiBEVYr1w4cL0L7/80rZLly7u8uecJtTIottLloTi5k3gs89a2hoO5wlFRU8Cb3JzmaAJhSx1g02B5uPKlW04e/ZnpKezWTk9PWMMGrQA/fvPgY7OEwdCLmedLwAW0DNoEAuoaS5ehKLbvA0TpyovXNFt7hFyWht6ekC3bky4SkrYe7kcuHkTOHYMKCw0hKPj/+DgMAvp6adw6dLXiIs7i4MHP8fJkz/C3/9jpeeoD6GQTcHK5azYd1ERMHYs/71/meBtmJoPjRRCvkbIac0IBE8KdAuFgI8Pe1HKkuuTkgguXx6INm0GoqTkNIKDv0Bi4mXs378QJ0+uwvTpe+Hi0rvifDs7JqYODkxoOerh5uZWpqneIAB4e3uXRkdHcw+wGeBrhBxOM0EIC8BxcwPeeYd5jrq6/bFgwUXMmXMU9vZdkZ//ED/9NBDXr+996jxbW+DQIRaY04IoFAoF/8vjaCTK390ak+80Ugj51ChH0xEIgBEjWFHu5GQCIyN/zJ37H/r0mQWZrAwbNryJkyd/rDheLGZe5pEjT8rItQARGRkZRlwMOZqGQqEgGRkZRgAiatqvsVOjHI6mIxIBb74J9OwJxMcDp08LMXLkzzAzs8P+/Quxd+88ZGUlY8yYlRAIhLCwABISWIcLR8fmt1cmk01NT0/flJ6e7gUN/RLNeWlRAIiQyWRTa9rZpEJICAkAsBqAEMAmSunyGo55E8CXYBUcwiml4+u7Ll8j5LwoCARPysJZWgLbthEMGPAxTExsERT0Lk6fXo3i4jy8887vIITA2Jh5hbNmNX/hbh8fn0cAhjfvqBxO09Nkf0qEECGAXwAMASAFMI4QIq1yjAuATwC8Qin1BPCBOtfmQsh5EfHyYnVHk5IAR8fxmD37GMRiPfz3XxD2718EADAxAR48YJ02eKlBDqdxaMrvlH4AEiildyilZQB2ARhR5ZhpAH6hlOYAAKX0kToX5muEnBeVvn1ZxRoHB0BHpx8mT/4LAoEIx49/hxMnVgIArKyAv/4Cfv2VTZNyOJznoymF0AZASqX395XbKuMKwJUQcokQckU5lVoNQsh0QkgoISQ0IyODrxFyXlgEAuYZjhsHBAYCZmYBePfdIADAvn0L8N9/W6Gry4SyqAjYsAEIC2tZmzkcTacphbAmn61qvJsIgAuAvgDGAdhECDGudhKlGyilvpRSXwsLC54+wXkp6NiRdcfw9n4bY8awCNKtW/8fLl3aAuBJ5/vdu4HLl1mOIofDaThqCyEhxIYQ0pMQ8qrqVc8p9wG0q/TeFkDVLsb3ARyklJZTSu8CiAUTxjrhU6OclwGxGAgIYM2DBw78AK+/vhgKhRxbt/4/HDjwKRQKBXR0WI5hcDDw7bfA4cOsawaHw1EftYSQELICwCUAnwP4SPlaUM9p1wC4EEIcCCFiAG8BOFTlmL8B9FOOYQ42VVqtH1ZVeLAM52XB25t5fg8fAsOHL8Hbb/8GgUCIo0e/xaZN41BWVgyxmPU2NDcHLl1i64d8+YDDUR91PcKRANwopa9RSocpX3WGUVNKZQBmAzgGIBrAHkppJCHkK0KI6txjALIIIVEAzgD4iFJaY0fnyvA/cs7LgpYWa9Zra8v6IkqlMzBz5mHo6Bjg+vU9WLWqP/LzH1Uca2fHCn8fO9bChnM4GoRa3ScIIf8CGEMpbfGmML6+vvTtt0MhEABz57a0NRxO80Ap61F48yZLqk9OjsCePa8jP/8ezM0dMHv2YbRt6wGAfVFMSgLeeAPo3v3Zxqup+wSH86KibkJ9EYCbhJBTACpWICilc5rEKjUQVmtRyeG8uBACuLuzl0IBZGR4YdSoq5g2bTgePLiGFSt6YMaMv+DhMaCic8XBg6zZr5NTS1vP4bRu1J0aPQRgKYDLAK5XerUIlLLyVBzOy4hAwKrQDBlihf/+Owsvr1EoLs7DmjUBuHhxEwAWaGNhAWzbBqSk1HNBDuclRy0hpJT+AWAnngjgn8ptLUZzl5ficFoj7drp4ezZvejX7yMoFDJs2zYN+/cvgkKhgETCmvn++iuwYweQqTEtdTmc5kXdqNG+AOLBSqb9CiBOjfSJJoNSPjXK4agwMxPg8OHvMH36BhAixLFjK7Bx41iUlRXDxIQF0CQmAlu2sIbBHA7nadT1q1YCGEwp7UMpfRWAP4Af6zmnyeBCyOE8ja4u8Ntv07B27b/Q1jbEjRv7sHr1YJSUFEAgYGXZ8vOBkydb2lIOp/WhrhBqUUpjVW8opXEAtJrGJPXgQsjhPA0hwP/+NwjHj1+GsbEtEhIuYu3a11BSwoK9ra1ZnuGdejN1OZyXC3WFMJQQspkQ0lf52ggeLMPhtEpefVzTlb4AABfrSURBVNUTZ86cgaGhzVNiKBQCZmbAzp18vZDDqYy6QvgegEgAcwDMBRAFYGZTGaUO3CPkcGqnUydnBAefhYGBDRISLmDVqn7Iz38IQ0PmOW7eDGRnt7SVHE7rQN2o0VJK6SpK6ShK6RuU0h8ppS1W0ZCvEXI49dO7tzN27DgLIyMHJCeHYvny7khPj4G5OVBezoJneF1SDqceISSE7FH+e5sQcqvqq3lMrA4XQg5HPYYNc8bBg1fQtq0fsrKSsGJFTyQkXEKbNkBWFhAS0tIWcjgtT30eoaqI2VAAw2p4tRg8j5DDUY8+fdrg+PEz8PQcgaKiHKxe7Y+EhIuwtmZRpAUFLW0hh9Oy1CknlNI05Y+ZAFIopckAtAF0RPWWSs0GpazAMIfDUQ8vLz2Eh/+F0aMnoKzsMVavDkBy8gUAwPnzLWwch9PCqOtXnQegQwixAXAKwGQAQU1llDpwj5DDaRhCoRC7dgVh5MiJKCt7jLVrhyAp6Q+cOydDfHxLW8fhtBzqygmhlBYBGAVgLaX0DQDSpjOrfrgQcjgNRygUYt++LejVaxJKSx9j69Z3ERTkiQUL/kREBO9vxnk5UVsICSE9ALwN4LByW4tm8vFgGQ7n2RAKhThwYAtGjgyCiYkjMjLicOjQ2+jb1xtLl25DeXl5S5vI4TQr6grhBwA+AXBA2VzXEayRbrOTlsaWLblHyOE8O+bmAmzc+A5+/TUGr722CUZGdsjKisbixZMgFotb2jwOp1lRqzFva4IQQpcsyUffvgZ4tcXKfnM4Lw7Z2UB2dhlWrvwKv/32tXKrXiSlj71a1DAOp5moL4/wJ+W//xBCDlV9NY+JNdkl5B4hh9NImJoCzs5irFu3DHl5eUhKSgJQxPtUcF4a6lvn26b894emNqQhKBQyvkbI4TQBhoaGMDQ0bGkzOJxmpU4hpJSqCmuHAiimlCoAgBAiBMsnbBFkslLuEXI4HA6nUVBXTk4B0Kv0XhdAi3Y24x4hh8PhcBoDdYVQh1JaqHqj/FmvjuObFIFAxD1CDofD4TQK6srJY0JIF9UbQogPgOKmMal+BAIeLMPhcDicxkHdpPgPAOwlhKjqi7YFMLZpTFIHyoWQw+FwOI2CWkJIKb1GCHEH4AaAAIihlLZY+QmhUIuvEXI4HA6nUVDLryKE6AFYCGAupfQ2AHtCyNAmtawO+Bohh8PhcBoLdeVkC4AyAD2U7+8DWNYkFqmBXF7GPUIOh8PhNArqCqETpfQ7AOUAQCktBpsibRHKy4u5R8jhcDicRkFdOSkjhOgCoABACHECUNpkVtUDIbzoNofD4XAaB3WjRv8PwFEA7QghOwC8AuDdpjKqPggRcCHkcDgcTqNQrxASQgiAGLCmvN3BpkTnUkozm9i2umzia4QcDofDaRTqFUJKKSWE/E0p9cGTprwtCqU8j5DD4XA4jYO6cnKFENK1SS1pADo6hlwIORwOh9MoqLtG2A/ATEJIEoDHYNOjlFLaoakMqwvej5DD4XA4jYW6QjikSa1oIDJZCViOP4fD4XA4z0d9Hep1CCEfAPgIQACAB5TSZNWrvosTQgIIIbGEkARCyKI6jhtNCKGEEF91jC4rewzSYlmMHA6Hw3mRqG+C8Q8AvgBug3mFK9W9sLJ57y/K86QAxhFCpDUcZwBgDoCr6l5bICBcCDkcDofTKNQnhFJK6QRK6XoAowH0bsC1/QAkUErvUErLAOwCMKKG45YC+A5AiboX5lGjHA6Hw2ks6pOTig4TlFJZA69tAyCl0vv7ym0VEEI6A2hHKQ2u60KEkOmEkFBCSOiTbQ20hsPhcDicGqgvWKYjISRf+TMBoKt8r4oaNazj3JqkilbsJEQA4EeoUaGGUroBwAblecoyb/WdxeFw/n97dx9cVX3ncfz9zQ3kJiQSSCCEBAnKw4JokOFBBLQ7ZWhtB5hWa61lalalO7iuLa5uq3a0q7OlPuzYkVotlo6Lnd1idUbQoVJ07YZSUKShiLRVFAYIaB6AEMgDIfnuH/dyN4kBLnBvbk7yec2c8Z6HHD4JcT6c3zn3d0Xk7M5YhO5+IfO37AdGtFsvBg60W88BJgK/j0xewzBgjZnNd/d3OYOsrMEaGhURkYRIZp1sAcaY2Sgz6w/cBKw5tdPd69w9391L3L0E2AyctQQB0tJCuiIUEZGESFoRRu8p3gmsA/4CvOju75vZw2Y2/0LO3dx8TEUoIiIJEe8b6s+Lu68F1nba9uBpjv1cvOdtaTlOWlr2hYUTEREhuUOjSePuuiIUEZGECGQRgp4aFRGRxAhkEZrpilBERBIjkEUI6O0TIiKSEIGsk+zsAl0RiohIQgSuCM2MtLQ0FaGIiCRE4IrQ3WlsrFMRiohIQgSuCAFOnDime4QiIpIQgawTdz/7QSIiInEIZBGKiIgkSkCLUFeEIiKSGAEtQhERkcQIXBGaGbm5xamOISIivUTgihAiZSgiIpIIgStCd+f48dpUxxARkV4icEUIkQ/mFRERSYRAFqGeGhURkUQJZBHqFqGIiCRKIItw9+6NqY4gIiK9RCCL8IMP3kx1BBER6SUCWYTTppWlOoKIiPQSgSxCs0DGFhGRHiiQjbJr11upjiAiIr1EIIswJ6cg1RFERKSXCGQRioiIJEogi/DQoT2pjiAiIr1EIIuwoeFwqiOIiEgvEcgiDIdzUh1BRER6iUAW4dSp30p1BBER6SUCWYRmoVRHEBGRXiKQRbh//9ZURxARkV4ikEV40UWFqY4gIiK9RCCLUEREJFFUhCIi0qcltQjN7Itm9jcz22Vm3+9i/91mttPMtpvZm2Y2Mpl5REREOktaEVrk0c6ngeuACcA3zGxCp8MqgCnufgXwEvBYPOfOyMhOZFQREenDknlFOA3Y5e4fu/sJ4NfAgvYHuPtb7t4QXd0MFJ/tpKFQOiUl0xMeVkRE+qZkFmERsK/d+v7ottO5DfhtVzvM7Ntm9q6ZvdvW1pbAiCIi0telJ/Hc1sU27/JAs4XAFODarva7+3JgefRY37t3CzA1QTFFRKQvS2YR7gdGtFsvBg50PsjM5gAPANe6e3M8J25uPpqQgCIiIskcGt0CjDGzUWbWH7gJWNP+ADO7Evg5MN/dq+I9sXV1rSkiInIeklaE7n4SuBNYB/wFeNHd3zezh81sfvSwx4Fs4Ddmts3M1pzmdCIiIkmRzKFR3H0tsLbTtgfbvZ5zPudtbj52gclEREQiAjmzTG7u8FRHEBGRXiJwRRgK9eOSS/TEqIiIJEbgihBg37532bFjR6pjiIhILxDIIly+/OssXbo01TFERKQXCGQR5uYO4+DBg6mOISIivUBAi7BQRSgiIgkRuCIcPvwKBg0q5JNPPkl1FBER6QUCV4QAgwYN48iRIzQ2NqY6ioiIBFzgivDIkUpmzvw669evJz09qfMBiIhIHxC4JmloOERR0WjmzBmd6igiItILBO6KEKCl5Tgvv/wyH330UaqjiIhIwAWyCE+cOM4NN9zA888/n+ooIiIScIEbGgXIyxvK9ddfz6OPPsro0aOprq7m9ttvJzc3l8rKSmpra5kwYYLuIYqIyFkF7orQLI20NHj22WfJz8+nrKyMe++9Nzbl2sqVKyktLaWwsJAf//jH1NfXpzixiIj0ZIErwuHDLyMtDfLz8/njH/9IeXk5tbW1zJo1C4CFCxfyq1/9iqlTp3LfffcxcuRInnjiiRSnFhGRnipwRQiQFk1dUlLC7NmzGTx4cGzfiBEj+OY3v8natWt55513mD17Nm+//XZs//Tp07n//vupq6vr7tgiItIDBa4IDx/eRygU37FTp05l9erVrFq1CoCTJ09SXFzM0qVLufTSS/nJT35CVVVVEtOKiEhPF7gibGw8ErsijFda9AvS09N5+eWX2bp1K5MmTWLJkiUUFBSwadMmAHbv3s2WLVv49NNPcfdERxcRkR4okI9VnmsRdjZ58mTWr1/P1q1b+d3vfseVV14JwC9+8Qt+9KMfAZCRkcGQIUMIhULs2rWL9PR0Hn74YVatWkUoFIotWVlZlJeXA/Daa69RWVlJSUkJBQUF5ObmxpZEqK+vp62tjfT09A4Z0i70ByIi0ocFsAgt7qHRM57FjClTpjBlypTYtkWLFjF9+nT27t3L3r17qampoa2tLVY0hYWFjB8/ntbW1tjSr1+/2NevWLGCV155pcOfU1xczL59+wC4/vrrKS8vx8xIS0vDzLjssst44403APjKV77Ctm3bcHeam5tpampiypQprF+/HoBp06bx17/+tcP5v/CFL/D6668DMGbMGA4cOEAoFCIjI4Pc3FwWLFgQe1ho8eLFNDc3k5aWRlpaGqFQiNmzZ3PzzTcDkSLPyMggKyuLzMxMsrKyGDJkCHl5ebg77q7SFZFeJ4BFeOFXhKdTUlJCSUnJafcvWrSIRYsWnXb/Sy+9xIEDB9izZw/V1dXU1dURatfa1157LQUFBbFSaWtro6ioKLZ/0qRJ5OTkABAOhwmHw4wZMya2//7776e6uprW1lZOnjxJa2srl1xySWx/WVkZhw8fprW1laamJo4cOUJhYWFsf3l5eeyqsq2tjZMnT5KRkcHNN99MS0sL8+bN+8z3dM899/D4449TX1/PwIEDCYfDZGVlxcpyyZIlLF68mNraWsrKymL7Bg0aRGFhIXPnzqW0tJSWlhaOHj3K4MGDMbMz/0WIiHQjC9q9sIyMbN+w4RjTpqU6Se/S2trK9u3baWho6LCMGzeOyZMnc/z4cZ544onP7P/a177GV7/6VSorK5k3b15se21tLQ0NDfzsZz9j8eLFbN++ndLSUvr378/QoUMZMGAA4XCYpUuXct1117Fjxw4efPBBMjMzyczMJBwOk5mZya233sr48ePZvXs3r7/+eod9mZmZTJ06lUGDBlFXV0d1dXWHrw+Hw3Ffwbq7CrodM9vq7lPOfqRI8AXuinDYsL9L2hVhXxYKhWL3SrsyYMAAHnroodPuLyoq4k9/+lOHbfX19bEiGjp0KE8++SSffPIJVVVVNDQ00NjYSHZ2NgDHjh3jgw8+oKmpicbGRhobG2lqamLOnDmMHz+eiooK7rjjjs/8uRs2bGDWrFmsXr2aW2655TP7KyoqmDRpEs888wx333137Gr41PLxxx8zatQoHnvsMR544IEOJRoOh9m8eTN5eXmsWLGCF198scO+cDjMsmXL6N+/P+vWraOioiJW0qf+e+ONNwLw4Ycfcvjw4Q5fm5WVRX5+/tn/ckQkqQJXhEBC7hFK8p0a5gUYNmwY3/3ud0977FVXXRWbHagrX/7ylzl48GCsIE+V5eWXXw7ArFmzWLlyZYf9TU1NDB8+HIDS0lLuuuuu2P3RU8vAgQMBmDFjBt/73vdoamrq8PXhcBiAlpYW6uvrqaqq6vDnP/300wCsXr2aZ555pkPmcDgcK8JHHnmEF154ocP+vLw8ampqgMiw+4YNG8jPz48tl156Kffddx8AO3fuZMCAARQXF3cYbheRCxe4odGcnCH+hz9UU1qa6iQi/6+trY0TJ07ECrSpqYkTJ04wbtw4AN577z327dvXYX8oFKKsrAyAp556io0bN1JTU0NNTQ3V1dUUFhaydetWAK6++mo2bdpE//79KSwsJDs7mxkzZvDcc88BkfvHLS0tFBUVxZZRo0Z1uEd8LjQ0Kn1J4IqwX7+wV1Q0MXFiqpOIdJ+NGzeyc+dOPvzwQz799FOOHTvG2LFjWbp0KQDXXHMN77zzDs3NzbGvWbBgQewp5rlz58buzxYUFHDxxRdz9dVXUxr9F2X7p6NBRSh9SyCHRnWPUPqamTNnMnPmzNPuLy8vx905dOgQlZWVVFZWctFFFwGRB4EyMjI4cOAA27Zto6qqipaWFu68806WLVtGU1MTWVlZ5OTkMHDgQL7zne9017cl0iOoCEV6CTMjLy+PvLw8rrjiig7bX3311di6u7N///4O6z/4wQ84evQodXV1jBw5sltzi6RaIIdG33+/ibFjU51EpPfS0Kj0JYG7tkpPz9BToyIikjCBK8KhQ8doaFRERBImkJWiIhQRkUQJXKXU1u7R0KiIiCRM4IrwxInjuiIUEZGECWSlqAhFRCRRklopZvZFM/ubme0ys+93sT/DzFZF979tZiXxnFdDoyIikihJK0IzCwFPA9cBE4BvmNmETofdBhx299HAk8Cj8ZxbV4QiIpIoyayUacAud//Y3U8AvwYWdDpmAfCf0dcvAZ+3s3woXL9+mSpCERFJmGROsVYE7Gu3vh+Yfrpj3P2kmdUBeUBN+4PM7NvAtyNraS2DB6e9B4GZEiefTt9PDxakrBCsvEHKCjAu1QFEuksyi7CrK7vO5RXPMbj7cmA5gJm9G6Spn4KUN0hZIVh5g5QVInlTnUGkuyRzkHE/MKLdejFw4HTHmFk6MBA4lMRMIiIiHSSzCLcAY8xslJn1B24C1nQ6Zg1wS/T1DcD/eNBmARcRkUBL2tBo9J7fncA6IAT80t3fN7OHgXfdfQ2wAnjBzHYRuRK8KY5TL09W5iQJUt4gZYVg5Q1SVgheXpHzFriPYRIREUkkvRFBRET6NBWhiIj0aT22CJM1PVsyxJH1bjPbaWbbzexNMxuZipzt8pwxb7vjbjAzN7OUPfYfT1YzuzH6833fzP6ruzN2ynK234WLzewtM6uI/j58KRU5o1l+aWZVZrbjNPvNzJ6Kfi/bzWxyd2cU6Rbu3uMWIg/XfARcAvQH/gxM6HTMHcCz0dc3Aat6cNa/B7KirxenKmu8eaPH5QDlwGZgSk/NCowBKoBB0fWhPflnS+QhlMXR1xOAPSnMew0wGdhxmv1fAn5L5P2+VwFvpyqrFi3JXHrqFWFSpmdLkrNmdfe33L0hurqZyHsqUyWeny3AI8BjQFN3husknqyLgKfd/TCAu1d1c8b24snrwEXR1wP57Htru427l3Pm9+0uAFZ6xGYg18wKuyedSPfpqUXY1fRsRac7xt1PAqemZ+tu8WRt7zYi/8pOlbPmNbMrgRHu/lp3ButCPD/bscBYM9toZpvN7Ivdlu6z4sn7Q2Chme0H1gL/3D3Rzsu5/m6LBFIyp1i7EAmbnq0bxJ3DzBYCU4Brk5rozM6Y18zSiHwSSFl3BTqDeH626USGRz9H5Ep7g5lNdPcjSc7WlXjyfgN43t3/w8xmEHkf7UR3b0t+vHPWU/4fE0mqnnpFGKTp2eLJipnNAR4A5rt7czdl68rZ8uYAE4Hfm9keIveG1qTogZl4fw9Wu3uLu+8G/kakGFMhnry3AS8CuPsmIExkQu6eKK7fbZGg66lFGKTp2c6aNTrU+HMiJZjKe1hwlrzuXufu+e5e4u4lRO5pznf3VEzCHM/vwStEHkbCzPKJDJV+3K0p/188efcCnwcws/FEirC6W1PGbw3wrejTo1cBde5+MNWhRBKtRw6NevKmZ0tV1seBbOA30ed59rr7/B6ct0eIM+s6YK6Z7QRagXvdvbYH5/0X4DkzW0JkmLEsRf+Aw8z+m8iQcn70nuVDQD8Ad3+WyD3MLwG7gAbgH1KRUyTZNMWaiIj0aT11aFRERKRbqAhFRKRPUxGKiEifpiIUEZE+TUUoIiJ9mopQ4mZmrWa2zcx2mNmrZpab4POXmdlPo69/aGb3JPL8IiJdURHKuWh090nuPpHIezf/KdWBREQulIpQztcm2k3AbGb3mtmW6OfW/Vu77d+Kbvuzmb0Q3TYv+hmSFWb2hpkVpCC/iAjQQ2eWkZ7NzEJEpglbEV2fS2R+z2lEJmpeY2bXALVE5led6e41ZjY4eoo/AFe5u5vZ7cC/EplxRUSk26kI5Vxkmtk2oATYCqyPbp8bXSqi69lEirEUeMndawDc/dSk6MXAquhn2/UHdndLehGRLmhoVM5Fo7tPAkYSKbBT9wgNWBq9fzjJ3Ue7+4ro9q7m8FsG/NTdLwf+kcjE0yIiKaEilHPm7nXAXcA9ZtaPyCTTt5pZNoCZFZnZUOBN4EYzy4tuPzU0OhCojL6+BRGRFNLQqJwXd68wsz8DN7n7C9GPFNoU/XSNY8DC6Ccv/Dvwv2bWSmTotIzIp7T/xswqiXzM06hUfA8iIqBPnxARkT5OQ6MiItKnqQhFRKRPUxGKiEifpiIUEZE+TUUoIiJ9mopQRET6NBWhiIj0af8HT0Qs13NiaAAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pr_point_up = ss.zoh(*zip(*pr_bounds_point.upper))\n",
    "pr_point_lo = ss.zoh(*zip(*pr_bounds_point.lower))\n",
    "pr_interval_up = ss.zoh(*zip(*pr_bounds_interval.upper))\n",
    "pr_interval_lo = ss.zoh(*zip(*pr_bounds_interval.lower))\n",
    "fig = plt.figure(3)\n",
    "fix_plot_shape(fig)\n",
    "plt.plot(*zip(*pr_true), color='black', linewidth=2)\n",
    "plt.plot(*zip(*pr_neg), color='black', linestyle='dashed')\n",
    "plt.fill_between(xs, list(map(pr_point_lo, xs)), list(map(pr_point_up, xs)), color='blue', alpha=0.4)\n",
    "plt.fill_between(xs, list(map(pr_interval_lo, xs)), list(map(pr_interval_up, xs)), \n",
    "                 color='none', edgecolor='red', alpha=0.8, hatch='xxx')\n",
    "plt.plot(*zip(*pr_after_labeling), color='red', linestyle='solid', linewidth=2)\n",
    "plot_proxy()\n",
    "plot_proxy2()\n",
    "plt.xlabel('Recall')\n",
    "plt.ylabel('Precision')\n",
    "plt.legend(['true curve', 'beta=0', 'after labeling',\n",
    "            'expected region via betahat', 'expected region via interval'], \n",
    "           loc=\"upper left\", bbox_to_anchor=(1,1))\n",
    "plt.title('Precision-Recall curve')\n",
    "plt.show()"
   ]
  }
 ],
 "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}