1_sklearn_at_scale.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "elegant-basketball",
   "metadata": {},
   "source": [
    "# Notebook 1: scikit-learn at scale\n",
    "\n",
    "The goal of this notebook is to give you a sense of some of the computational limitations one may encounter when training or applying ML models on large scale datasets, as well as to explore some of scikit-learn's tools to improve runtime.\n",
    "\n",
    "__If you are taking this class for credit, you will need to turn in your completed notebook on Campus.__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a8349e8-c2be-41fd-9891-7bd32be57cf1",
   "metadata": {},
   "source": [
    "This notebook was prepared on Jupyter Lab 3.2.9 by [Chloé Azencott](http://cazencott.info), using Python 3.9.7 with `joblib 1.1.0`, `matplotlib 3.5.1`, `numpy 1.21.2`, `scipy 1.7.3`, and `sklearn 1.0.2`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "willing-london",
   "metadata": {},
   "outputs": [],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46b1c4c1-b207-42ee-a157-c4c4c0d960a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rc('font', **{'size': 12}) # règle la taille de police globalement pour les plots (en pt)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "exceptional-county",
   "metadata": {},
   "source": [
    "## 1. Data setup"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "hearing-system",
   "metadata": {},
   "source": [
    "We will work with the [Forest Covertypes](https://scikit-learn.org/stable/datasets/real_world.html#forest-covertypes) data set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "retained-sentence",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "compliant-invitation",
   "metadata": {},
   "outputs": [],
   "source": [
    "forest_dataset = datasets.fetch_covtype()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adolescent-ribbon",
   "metadata": {},
   "outputs": [],
   "source": [
    "X = forest_dataset['data']\n",
    "y = forest_dataset['target']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "listed-zoning",
   "metadata": {},
   "outputs": [],
   "source": [
    "X.shape, y.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "relevant-paragraph",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_samples = X.shape[0]\n",
    "n_features = X.shape[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "together-hometown",
   "metadata": {},
   "source": [
    "#### Creating subsets of varying sample sizes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "accessible-naples",
   "metadata": {},
   "source": [
    "The data set contains about 500k samples. Let us create subsets of the data, to explore how sample sizes affects complexity. We'll use the `train_test_split` function of sklearn to that effect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "behavioral-agreement",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import model_selection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "close-monitoring",
   "metadata": {},
   "outputs": [],
   "source": [
    "subsample_sizes = [500, 1000, 5000, 10000, 50000, 100000, 200000]\n",
    "subsample_data = {} # key: number of subsamples, value:(X_subsample, y_subsample)\n",
    "for ssize in subsample_sizes:\n",
    "    X_tr, X_te, y_tr, y_te = model_selection.train_test_split(X, y, shuffle=True, random_state=57, stratify=y, train_size=ssize)\n",
    "    subsample_data[ssize] = (X_tr, y_tr)\n",
    "del X_te, y_te\n",
    "\n",
    "# Adding the full data set to the dictionary\n",
    "subsample_data[n_samples] = (X, y)\n",
    "subsample_sizes.append(n_samples)\n",
    "\n",
    "subsample_sizes = np.array(subsample_sizes)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "353b2a23-f977-4508-a37b-15a6d4d73cb2",
   "metadata": {},
   "source": [
    "__IMPORTANT NOTE:__ If later in the notebook cells are taking too long to run, reduce the values in `subsample_sizes`. For example, consider \n",
    "```python\n",
    "subsample_sizes = [500, 1000, 5000, 10000, 50000, 75000, 100000]\n",
    "```\n",
    "or\n",
    "```python\n",
    "subsample_sizes = [500, 1000, 5000, 10000, 25000, 50000]\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aggregate-print",
   "metadata": {},
   "source": [
    "#### Size of each data set in memory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5049bfad-f476-4f0f-99c8-46967f5d5d35",
   "metadata": {},
   "outputs": [],
   "source": [
    "for ssize in subsample_sizes:\n",
    "    X_current, y_current = subsample_data[ssize]\n",
    "    print(f\"Dataset with {ssize} samples takes {(X_current.size * X_current.itemsize):.2e} bytes of memory\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "colonial-cigarette",
   "metadata": {},
   "source": [
    "Our datasets take from 216 KB to 251 MB of space in memory."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "speaking-newport",
   "metadata": {},
   "source": [
    "## 2. Decision trees"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "leading-costume",
   "metadata": {},
   "source": [
    "### 1.2 Training a decision tree: complexity w.r.t. the number of samples"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "suitable-words",
   "metadata": {},
   "source": [
    "The computational complexity for building a decision tree with CART is about $\\mathcal{O}(p n \\log n)$ if the features are binary and $\\mathcal{O}(p n^2)$. This is in the case where the tree that is built is balanced (about the same number of training samples go to each branch at each node). \n",
    "\n",
    "In the worse case scenario (a single training sample goes to one of the two branches at each node), this complexity becomes $\\mathcal{O}(p n^2)$ if the features are binary and $\\mathcal{O}(p n^3)$ if they are continuous.\n",
    "\n",
    "For details, you can consult [this document](http://cazencott.info/dotclear/public/lectures/Decision_Tree_Complexity.pdf)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "divided-referral",
   "metadata": {},
   "source": [
    "##### __Question 1__\n",
    "How many of the features are binary? Do you expect the runtime to be loglinear or quadratic in the number of samples?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f229b9d-851d-4e5d-b060-44e462fc77f4",
   "metadata": {},
   "source": [
    "__Answer:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "increasing-guidance",
   "metadata": {},
   "source": [
    "#### Empirical runtime as a function of the number of samples\n",
    "\n",
    "We will now evaluate this runtime empirically, on our subsets of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "binary-savings",
   "metadata": {},
   "outputs": [],
   "source": [
    "import time # to compute runtimes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "olympic-yahoo",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import tree "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "moral-marking",
   "metadata": {},
   "outputs": [],
   "source": [
    "runtimes = {} # key: number of subsamples, value: time to train a decision tree with default parameter\n",
    "for ssize in subsample_sizes:\n",
    "    X_current, y_current = subsample_data[ssize]\n",
    "    dt_classifier = tree.DecisionTreeClassifier() # default parameters\n",
    "    \n",
    "    start_time = time.time()\n",
    "    dt_classifier.fit(X_current, y_current)\n",
    "    runtimes[ssize] = time.time() - start_time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "deadly-cookbook",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(subsample_sizes, [runtimes[ssize] for ssize in subsample_sizes], 'o-')\n",
    "plt.xlabel(\"Number of samples\")\n",
    "plt.ylabel(\"Time (s)\")\n",
    "plt.title(\"Time to train a decision tree\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "conventional-peoples",
   "metadata": {},
   "source": [
    "To better visualize the first data points, we can switch to a log scale on the x-axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "nonprofit-germany",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.loglog(subsample_sizes, [runtimes[ssize] for ssize in subsample_sizes], 'o-')\n",
    "plt.xlabel(\"Number of samples (log scale)\")\n",
    "plt.ylabel(\"Time (s, log scale)\")\n",
    "plt.title(\"Time to train a decision tree\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adjustable-kansas",
   "metadata": {},
   "source": [
    "#### Fitting a complexity model\n",
    "\n",
    "It is hard to conclude on the computational complexity of the algorithm in practice from the above plot. In this section, we will check what type of computational complexity function fits our observations best. For this purpose, we'll use the [curve fitting algorithm of `scipy.optimize`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "plain-relation",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import optimize"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "embedded-tragedy",
   "metadata": {},
   "source": [
    "Let us define now 3 types of computational complexities we are willing to consider:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dramatic-treat",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f_linear(n, a):\n",
    "    \"\"\"Linear complexity\"\"\"\n",
    "    return (a * n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "formed-parliament",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f_loglinear(n, a):\n",
    "    \"\"\"Loglinear complexity\"\"\"\n",
    "    return (a * n * np.log2(n))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "overhead-builder",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f_quadratic(n, a):\n",
    "    \"\"\"Quadratic complexity\"\"\"\n",
    "    return (a * n**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "practical-millennium",
   "metadata": {},
   "source": [
    "Now we use `curve_fit` to fit the data to the models. What is happening under the hood is an ordinary least squares approach (minimization of the sum of squares)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "worldwide-orleans",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_opt_linear, p_cov = optimize.curve_fit(f_linear, subsample_sizes, [runtimes[ssize] for ssize in subsample_sizes])\n",
    "p_opt_loglinear, p_cov = optimize.curve_fit(f_loglinear, subsample_sizes, [runtimes[ssize] for ssize in subsample_sizes])\n",
    "p_opt_quadratic, p_cov = optimize.curve_fit(f_quadratic, subsample_sizes, [runtimes[ssize] for ssize in subsample_sizes])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "median-category",
   "metadata": {},
   "source": [
    "And now we plot the learned models as well as the observed runtimes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "inappropriate-palestinian",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10, 4))\n",
    "\n",
    "# Linear axes\n",
    "plt.subplot(1, 2, 1) # grid 1x2 subplots, subplot 1\n",
    "plt.plot(subsample_sizes, [runtimes[ssize] for ssize in subsample_sizes], 'o-', label=\"empirical runtimes\")\n",
    "plt.plot(subsample_sizes, f_loglinear(subsample_sizes, p_opt_loglinear[0]), 'v-', label=\"loglinear fit\")\n",
    "plt.plot(subsample_sizes, f_linear(subsample_sizes, p_opt_linear[0]), 'd-', label=\"linear fit\")\n",
    "plt.plot(subsample_sizes, f_quadratic(subsample_sizes, p_opt_quadratic[0]), 'x-', label=\"quadratic fit\")\n",
    "\n",
    "plt.xlabel(\"Number of samples\")\n",
    "plt.ylabel(\"Time (s)\")\n",
    "plt.title(\"Empirical fits of runtimes\")\n",
    "plt.legend()\n",
    "\n",
    "# Log axes\n",
    "plt.subplot(1, 2, 2) # grid 1x2 subplots, subplot 2\n",
    "plt.loglog(subsample_sizes, [runtimes[ssize] for ssize in subsample_sizes], 'o-', label=\"empirical runtimes\")\n",
    "plt.loglog(subsample_sizes, f_loglinear(subsample_sizes, p_opt_loglinear[0]), 'v-', label=\"loglinear fit\")\n",
    "plt.loglog(subsample_sizes, f_linear(subsample_sizes, p_opt_linear[0]), 'd-', label=\"linear fit\")\n",
    "plt.loglog(subsample_sizes, f_quadratic(subsample_sizes, p_opt_quadratic[0]), 'x-', label=\"quadratic fit\")\n",
    "\n",
    "plt.xlabel(\"Number of samples (log scale)\")\n",
    "plt.ylabel(\"Time (s, log scale)\")\n",
    "plt.title(\"Empirical fits of runtimes\")\n",
    "plt.legend()\n",
    "\n",
    "# Adjust spacing between plots\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adopted-macintosh",
   "metadata": {},
   "source": [
    "##### __Question 2__ \n",
    "Is the relationship between the time required to train a decision tree and the number of samples linear? loglinear? quadratic? Does this match your expectations?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a8db873-c27f-492d-8bdb-c46cc7bc0a1e",
   "metadata": {},
   "source": [
    "__Answer:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "caroline-beginning",
   "metadata": {},
   "source": [
    "#### Size of a tree in memory"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "superb-world",
   "metadata": {},
   "source": [
    "We will now use `joblib` to save estimators to disk and see how much space the models take to store.\n",
    "\n",
    "You can read more about `joblib` [here](https://joblib.readthedocs.io/en/latest/), and about *model persistence* (ie. how to save a trained model) in scikit-learn [here](https://scikit-learn.org/stable/modules/model_persistence.html).\n",
    "\n",
    "If `joblib` isn't installed on your system, you should be able to install it with `conda install joblib`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "known-calibration",
   "metadata": {},
   "outputs": [],
   "source": [
    "import joblib\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "turned-background",
   "metadata": {},
   "outputs": [],
   "source": [
    "tree_joblib_size = {} # key: number of subsamples, value: size of file on disk\n",
    "for ssize in subsample_sizes:\n",
    "    X_current, y_current = subsample_data[ssize]\n",
    "    dt_classifier = tree.DecisionTreeClassifier() # default parameters\n",
    "    dt_classifier.fit(X_current, y_current)\n",
    "    file_name = \"dt_classifier_%d.joblib\" % ssize\n",
    "    joblib.dump(dt_classifier, file_name)\n",
    "    print(\"Tree size: %.2e bytes\" % os.path.getsize(file_name))\n",
    "    tree_joblib_size[ssize] = os.path.getsize(file_name)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "computational-territory",
   "metadata": {},
   "source": [
    "##### __Question 3__ \n",
    "Complete the code below to plot, as a function of the number of samples, the size in bytes of 1) the data set and 2) the model file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "brave-briefs",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.loglog(subsample_sizes, # TO COMPLETE, \n",
    "           'o-', label='Dataset size')\n",
    "plt.loglog(subsample_sizes, # TO COMPLETE, \n",
    "           'v-', label='Model size')\n",
    "\n",
    "plt.xlabel(\"Number of samples (log scale)\")\n",
    "plt.ylabel(\"Size in bytes (log scale)\")\n",
    "plt.title(\"Sizes of models and datasets\")\n",
    "plt.legend()"
   ]
  },
  {
   "attachments": {
    "567a50c8-7e3b-4e45-8910-bde5e19de09d.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAACKCAIAAAD5baOGAAAAA3NCSVQICAjb4U/gAAAgAElEQVR4Xu2dB1gTyRfAdzcEpITeOwgKSpcOijQFxQIKKhawYMHesXLqWQ/E81AQRO70OP8qxUOUs57YQVHU8ywgICpdWuhJdv+bhJKQAhcQUWc+P79k5s2bN78J86btDvz69WtI0IBhWF1dXWNjI41GE1QHyAcIAAKAACAACHxJAkK9KRyGYSlG6I0SkBcQAAQAAUAAEPiCBJAvWDYoGhAABAABQAAQ+OIEgCP84k0ADAAEAAFAABD4kgSAI/yS9EHZgAAgAAgAAl+cAHCEX7wJgAGAACAACAACX5IAcIRfkj4oGxAABAABQOCLExD68OHDFzcCGAAIAAKAACAACHwpAnBNTQ3/srOzs83MzPjLgFRAABAABACBb4BATk7Od9XhM+sLlka/gZ8uqAIgAAgAAoCA4ASAIxScHcgJCAACgAAg8A0QAI7wG2hEUAVAABAABAABwQkARyg4O5ATEAAEAAFA4BsgABzhN9CIoAqAACAACAACghPo1Uu3BS8W5AQEAAFAABD4BgigZXdPRJ9/QzDxWzLLUg7+OmskqCPE6u9HbU34lyYxzHfVwlHKbWqoOfEb4rJbEc0Jm9aNVQazza/zNwGsBgQAAUCgKwG0NvfWpSv3XxSWkVuRQdJKWkPMRo5xMZavf/XoeUUTDcp+lDfFUk6sa7av43svnRVW/zI5LjW3Cfs6agusBAQAAUAAEPjPBLDGl+cORSbdffmxpolKo7U2fHr/b+b93DoEgqWG21tpysrpONgbiP5nvQMlg6Azwg77MWrpzd+T9NfPNJb4SifFA6UpgB2AACAACAxIAljt42sPKmkYLGnoMdXNQApu+lT08q2QpTYBN1featYGqwFpds+N6rUjhCAM/ZT5v7OGmnNHdL2gF/1469f/3covqyI30wji8lpGIydMdtIVxx0mVv08LSXjVXFpBT1NiKRiaO9kLJR378HzoioKUUbTzN3Xx159UJtrpdW8unEh/e7zomqKsIyGoZ3nJDdDGXoLgAAIAAKAACDwuQlgtVU1NHzdj6Bs7GA+VBrvmAcPMbJmlkorPP/jwWsVkNq4jRvHqWFvzob+cquGbY2QMHjq9lWj8e3DgduT93JpFCKoDB0ijUC1T5ISs2s5FkjhhvKi0prGVhpKayWX5d5Piv49k4EIrX37OOfVu7a0ltoPOekJpy5kvq1opNAojZVv7505lvKymcEZIz87fSgqNfNtJT2toTL/UVr0L8mvwWrs5/7tA/2AACAACDAIwLJKisK4+6Plpv4Sm3o/91MLR2fPmxRMEBEhDvCevNczQpLp1CnKcVG3ynKSk57omLPRQJRsZ60yHSRDGkRo/njrVOxf+U2vnrxqtLEVbxNDFEcuWDha4tX52KTnZEhIZ8ySGVbIi8TY1FeNNc9yCvyGGRKohddSsipRRMk+YP4EA7GqrP9FJ/1Tcf96jscQOxJYjOX96wMpgAAg8P0SuHr16p9//sm1/uMYgWsSr0iYZDHB62FB8uv65rJnVxKeX0tRNhk9YbK7iTzu4dgCQc9764EJuJvEap/8fvj0szpIwtB7irUkTC0YyD15rx0hBIkNneg/+s3hGyU5yX8SFNkGCsggofpnmdnFn+pbYBqGI2tG6+vqMajdEUJEkpyikpqszdC0fx41Y0KSiirKMjI2gy+9fk5pqqnBZ31iJf/8+wmFIGFaaWZqQhY+JKnGsWPUkg+lKEQCy6O8frggHhAABL5nAu7u7nj1OX2hAF6QgZGoOnpJiNrda9dvZf1b1og2luSkH3/73n/dQltpdswIcZAoEUIrbiedf14HERTsZ89xVCJA6PsB3ZP3gSOEIJHBntOcXvzyd+njhyxXWWBNr5PCozNKKWy+EcM459QEMdFBuJOEqFQqnbioGO4yKRiVSsNdXkM97jghrLUy/0UlC3AqhS4KAiAACAACgABXApy+UFAvyFRPkNYfNVV/5KTqvHtp51KzilvIL67dfW89nrNw6sfrp/581QgJa44JnDKccYxygPfkfeIIIWjQYM8pdk+j7tAnb+2h+cWte2UUiKjjsWzBGG2xsr/CwtI/cjJjxHSscTK8ZNs3usuExcTxwzVkiGSzMHS2Me4uQQAEAAFAABDoEQFWX9g7L9hRHEyU0XeaMen9i+gH9Vh1VTUKdT0kSXl/NSG9oBkSHTJhjoeWMDPnAO/J+8gRQrCowfhJZk/jH5PbJ3wYtaUVxb8gwhKSYsIEAlFIkIM5iMpQfalrpdX1jxJPKqHuRkoilLrK4gqigb2hDNgh7NHfAhACBACB75cA7gspFApe//+6L8iGDH1/Ne5SjY7pUHV5CSKtJvf260Y8HSZJSXbthqkfrp2+8r4VnxyNmWwt1drURJcTGjRoYPfkfeUI6UzMJowd/E9SXivTFcLiegaawq/yW14n7VyXKiIMUZvx6eJ/dobEIWMmmj4/lVPz6Vlq7LNURuPAopZB1oYyXbdpv9/fOqg5IAAIAAK8CPTKBTKUou+fPHj+vOzp8wyWMmDRwaPsNAnQe5Y4rP7JpRvv6dthzW//PLCReVoHUfFYv8lrQPfk/9kx8WKNxyMK9hMcFZH2EQKiNHrufE8LbXkxIkZpaUGFJGTV9IwM1cS6DiH4qMSTYBmrwLVLJ9sNVZEaJIQQiGIyKrqG6uL0h1pAAAQAAUAAEPj8BGSNXFwth6rLiQsTYERokKSirrnbrJWLnJW7HlhsaWhomwp1NWpA9+RwTQ3L+ZauptO/Z2dnm5mZcUsBcYAAIAAIAALfFIGcnJzvqsNn1rcvZ4Tf1M8BVAYQAAQAAUDg+yAAHOH30c6gloAAIAAIAAI8CABHyAMMiAYEAAFAABD4PggAR/h9tDOoJSAACAACgAAPAsAR8gADogEBQAAQAAS+DwLAEX4f7QxqCQgAAoAAIMCDAHCEPMCAaEAAEAAEAIHvgwBwhN9HO4NaAgKAACAACPAg0KMH6nnkBdGAACAACAACgMBXT6BH7xp1cXH56isKKgAIAAKAACAACHAjAJZGuVEBcYAAIAAIAALfDQHgCL+bpgYVBQQAAUAAEOBGADhCblRAHCAACAACgMB3Q6BHe4RdaNy9e5dKpfJBJCQk5ODgwEeALYnW2goLC/faI6PNVcVViLKqtCBV6qmt37YcRm2lIsJEzqbgzZZnlm+bFKgdIAAIfEsEODu97mvX2tpqzzfgAt1rYUi03trhNmp+QhF+ZW/vAlaRtNJjyo/Xa/vkmkL0Q+IyN9dFpwuZdmF1z09vC/QabTvCwsIx5Hoz9d/YWc6eIekVvTa7d5XuUe4udeGVByNf2ezmFJxYzlknVras2vhk4VUKiAcEAAFA4PMR6GF3x2GAINMnFKV3lo8fP8b/t7Cw6PIZ/8oU4CiLSwTaRK5nXuOI1eYkxqVRXDf6mwpy8zyG9YkHbDeR2tzU3ExpptC/Y/W3w9eGXSe5Ba5brD6IqGoqjBU3NTQ1NbbQuNRIwCjq458mBf1R3OGFhIxWnPttrjZjnIJWZp3YH3H2XkG9mJad78qN8+0V22/D5JPUYQhrXfhYR2usq6dfLM0Z2NiyauOdhVMJiAEEAIFvmQDW8PbKb8cS0rNyy1tElYeNmrZ89QwzKfwWdqw2dcWY0DtdJkcirvv+DhsrykqkuSgj4cQfl+7886EOImmYugesXjZBX5zjHvfW26FjVt0fE3VxszXTVdCK00Lm/nBfeWF01MJhrB0U9dmxBZv+l1tR10SFCINI8ur6ZiMnz507bogEh1JBHGG3TqdbAS6/B6wqO+V0mprRen8uif0dhWjPis2YAREIdEdEyblxq0rVe9+uhRbCbYbIr0i6uxQmEDhwCmwohntegr7PpnnWJIZSWGqoAlM75c2va1ZHfdD3nrdcrexGQsy6NXDciSBDuil8kjoNYauLwPa1Z2TVxtVr9roEoAAQAAS+PgKt2ZErt6YJO4yfsUQVLspIPB++qkb8zMFJSggsZj5r2w73jkE+Wnw9KvaBlK5GlxkPRi3MuPyGZD9jufqg6peXE8/uWNEqn7jNXpwvDKw26+c1ezKEx+8NW2AkBkGsXTdaW1xYihkFbJs4mNBSV1749MbF37bfzi6Pj52r36VwQRwhc8LHORfsmCP2fEbIt4pfNBEhtM26sKaqqkZEUV2FDVVHch8ZiTU2NCKa1l4eYwexayTfij/1r/i4sJ+3jMZHV5OH0qauSPj1pt++MVIwnyR2HX1rbN9q6yN+QA0gAAh8SQLCFoF7j/kPtdBgzBamDKf5rEi/er924mQZmKhh46XRbhtaen5VeLXW9P3zhnVxPrDEqK3/G4kgjAkANlGveUJIxvVnIfZ2fJYIm9/8vmHDH+Xm62K3uCq07fOxd1CwqK6D1/gRjLJmzPY8Omv2idT013P0jdpX1ZiGCbJH2K2f4yvQlJ9+cLmvm4ONrZNX4LY/81lP3bRc22hrTg9W02NyOTer2lnSVwsf/rZ93sTRdjZ2o8f5hyQVsAljdXcPBkx0dbS2tHZw91v58/X3zGk5WpX966bZniNtrGxHeUxdfDyHHs0tEqs5v8zGev7pUrpaDF8Jpjw+MG4E3TCLkVszWtH84/5WtusuNzENohTfit4w09PRxnaUV+D208/b9ylbXiXvXjnH293RxtLa0X121FPeB4wwch0ZE5cksTcOrr712Z3MOmkHT3v6GgMEiVuNHSVfn3mbbjmfpE5QuPmsdcETaGUP4rcGeDnZ2jiOnb7m8JXC5k5pWumNsMXeLrbWds6TF+5OfFHHMenrqg3qPgurMeAzIAAIfJsEECXjNi+I109UTVMBQZsam7r2II1ZcTH3iW5L55l2GfEzoLR5QfwTLCIuJgTDCD8HRX2funXF4Rea8w7unza4bbWOo4Nig01UVJQjQNy20QSZETJXPvnsEfJeGkXLLm1buO2WkMWkwCl6olWv/r6QjkJ6HbYSTWbvnG2C2wSTBqvyXHfE6h78FLTybI2ep2+wsRxU/YGqIssmDItKaph4BE5Qk4TKMs+ePLlls6TOb3N1aq/sXRv5ZLDvku0W8tinghJZWSEIq+YSyflDJehODlnsIA1DMFFluBBU1ymBVWfsWbjuIs1q+rLN2g0Pz8SHLa8T+9/BScoI1vTqRuqdcgv/oBk6EmgtRUuVd6Ni5Bq6+8y++BdFc4iJsa5M2yAIrSgorCeo62q2txNBU0cDaXhXWIaOJPBMctLkVRJWe2//wlXJDcMnzVqpT/x4O+mPTQsKmk+GT2TSxsreFhtPmrNKCft473zi3iXvWk9G+WvzbAg6Be5ZOBw6J1IQAwgAAt8oAbTkwb18TG26qSJ7T4SWpZ9Kr9QLmEtf3uIeMGpzfW1Z7u2Th69QTOZPMec1HaSVXtu1e/cdCZ+wiMVmnHt+ncoxSiOZTCZQm2pL/r0SE/eIaLhwzBCODkoQR8h3wkc3gKcA7XVi3E3ysCUJR+frMSo4dTB59NbSDqMRRSNnNzcR7oTaY9HCxMikD1oBJ+JWGHeMKtBy1kxEY9+1xswILzO4YEbsnczyAM2Kl6/rJezmrPJ36iyB+oFLJNR1GAMh0kMcXd2U21oV7XSEaF5idFrZ4LknDwXTt+0mWIuVTo1ITC+aMJfpPwhqo2bO9mvPyLNiWDNBWlM8P2nflvgWjKhgOfuHPcH28ghEnymisDKpo6ERkqQEglXWkjGMwDOJZzFoQWLk+WLtwF+jlw2nM/B101w2/UD0qcceG0fQMwmZBoXvmaFCr6afp84Kv92/nsz02W7LUx/PLPbchnv81IA0QAAQ+EYItL79387oJ6Qx+2Z2Wf1E36UlP8SsQiYN5nBE7VVH351a4Hv4BRVCZO3XRU83aD+W0RVNTfq+0GbUYHnCekc5XsN+Rh60LGmFS1J7doRkGrSXm1a+OroW3vadOeHD9wjxgM8L8cD6GRfiNSNEK3Iev4eHuLjq8nLzPEpkjcZXRTPfYPpjvYZ339niq5oEzaF6Yhi5tg4jaJmbypKvHdoen/G2tv3AJ9fIHljRJoKWZ2e9hXSdnLUh/JmR1laqnIm5Bpb/Krel5zrokgSDgCPnUtNv3rmZErXORfx5fEjoeca6LI8A8xpP4ZNp3kloxcPMPBzdOMO2kQBBzWO8JbE868HbjgOw7bkRJVcPc2LV0ycF3Z6NFSALj2qBaEAAEPi6CdBKLu9cdegf1Tl7N7nKsfdFtDfp6W9EbMa7tO/mcakpojJuy5HDB7Ytcpd7FrZg2clcHk/iiepZmspTX8TvOPLgE5+OEu8PZV3XR8fGxsYciwwLDR6vURi7ZP6hR/Vd5zqCOEKeE772evESwOpqySgiJYMvMQoesJqqahSRlWNfDWXTRyu7FxsyZ4KLvZWlle3YH/5uwH0zTkTSKeTIzqmqL+PW+HpMWHQgLa+RR2TPjcPqampR6psYf3sbZrCfdTyPSiWTG7uC7plOIZKWrf+OUF/1hkd/3cIfU4QlJEkIVk/uaDcU/4zC+LQQ5pPEqyystrqGHR0sKicnjtXg0Rx5YDEZmUH4ki2ZM4lDti1CgCy8VIF4QAAQ+OoIYFX3floeehX23BmxzIJ5/r2zDrS8GzffiVi5OvBcFqXLiigZWo9091m8N3aHp3DOibjb9VwpiBjMOnwy3Fcl7+Ta5Ycf8nl8HBZWGDrC0tLSytrBdXLQjqhtY0XfnkvI6Hr8QfClUT57hLwcISwlI4WgnyoqMYjTFzL2RXnNJVlYwCQpEoLifTc3JXQ59GPi9vXHCoxnL99jpy0t1JB5LOR4NUMBLK4/fu0Rz6C8vxMOhceHLqfIJ2+xFeUSacOVPZdIhjFCQ2cc2DxGvtO7wyQ1vm3NRRFrFFFLV12IVlGFt66Soo62BO12fhEVUmcsEtCK8t+j4lbaSvi2Ms8kXvqZ/Ks+VXegw5o+fWrAo6U5B0RYfUVFMyMHL3Uc8QJk4dABIgABQODrJEB7n7R9a2K1VUjcVlcljm4DLcnKKiQY+VlL9mgeBJOGGWkTLn8oqkQhCQ5ldEBCyk4bjoYhwetOrV8rEx0ZYND9EiE+s9DWUUau1VTjUwu2PpprCd00Q7feipcAomBpowfnXkrO5piZ4haSJBC0ohyvNv+AKJpZaGBvLl98yXLcEaK7URijUPBHwmkFL1+3KI9ZtNR7tLW5mYWNsSqRjTwiqee6ZG+wPbHy+fOP7aVxjeRvCD0VN8ZcA3r/pkx2mHFnMNLGj+EIHDDy85y3VBE1DfpviWjiYE2qufvXfeYIpuHh5duV4taOZrhX5JMEYQ3Fr96UsQJiWMvkfyX9VdvKLVpy5VI2RcHKmnPNvjUvOeUxRcnKRofAwpZvpViy8JUDiYAAIPDNEcA+pUdEZopP2vHjFG1uW1+N/+S8wTRMjPks5bEywV/nlZNPE1FW4bcFiMjZr43Y7CzyJHLDT7fx8X23obXwzv0CVFxdQ7aL5xOkw2ZO+AR5jhDR9V3mfXH1mTVzy6dMsNEQo1VlvaVCjAcmYRlTc20s5uSuMNRFGy5/J+q62rfrgybMahKGTg/2SNt4Yum8jz4eZqqDmsvfi45aOVFJRZx2JSU60WCpncEQ4QuXft6v6G2pKk5ofVHEfHcNWnQhIqFE3VhLlthckpWYQ5Gw01XCii6Ec0TCULf7Yu3EcWMWu1/cFBG06PUEp2FKorT6siIhuxUzLXowOuloNbQ0ZXvoIzlrIw384cB3WWkpGdVa/jOZB6tIToEzDW5F71otWjhGrezaqfPVQ4ICnRmDGZ5JWN21nf4h14jekRe3sj2Eg/NfOunC6viVwdXT3PWIH+8knssUcdweYIn7VXydGKIV3zlzmqgu0lL29K/ESy8lXXfNMRfGj/B2sF02RbvDauYH7lm6CIGvgAAg8G0TwD79ffFevbKzbuPD69faqgpL6NrY6DJP+tEKXuW1ihgO1WI7J4NVJC722vcp4PczwQZo9tFVv1cPN9NTFMdqC7Munc9oGDxn5igSf24E1fE79ubODf59xx7LP/Z7KnSVxpry715I/QjjZ0cr3z25np7xlmi6fLYd2ytt8DyCOEJeE74OE3gLwFL2G2OPahw5lpQWfauagohLKxva2+rTV5MRvVmha9/9GJ8a+RNVXGW4jxX+CErXVWZmCYica2hchNYvsanJR6/WwyRVw8lmTcK2c9dOyo+8nnjXb8qi0D2V4TEXY/ckN2IiJBl5dRNrbQkYq28seZ6S+md5PVVYUtVw5Krw1c4kKI9LJAzVdMXJ8zsi7/7D8UHakXFp56LSyLRBcprDxpvipg/q0fy/TS0mrSD+/sbpmyU1VGFZDSP3VSHB/pZtL1QgGsyPCIf2Hzp77GCDmJbNgrCQBcPajlLxSoJFFDRUSFKi6vgzM/RA3yBlmgNLOWw8dkjxUFTKbxF/opLaltP2rFo8lvFch4i25UjT0hcpkVlkCoGkYmC36ODSwFGM864inWx9tDtJwHyydEqBT4AAIPDtE0BLPhRTaB/+Ctv4V0dlCXpBCf8LHsrohSilxRWQnLMiu8fBaBQqhkhK0Xt6ioisfFNGWvyFigZITF7HxCtk+6IppvjLYroJsLhFcOiczPlx4eFXrfeybmshksqaitizUz8+pGGwkJi0kqaB06KgOTPd2p867NQM19R00+lnZ2d3uaH+7Nmz3t7efKxLSUnx8/PjIwCS+o8Amhvj739C5YerEeO7GVv1n02gJEAAEAAEIIj6PMJn3hWHqNSNltyWU/sPkSAzQvyWpeTkZD42Eom9rhTW2lDXyPkOaBgWFpcSE8RoPub2T1L/1wgtfXz1Uf7rS4m5RJPJZvzf2Nc/DEApgAAgAAh0EEBLnuSUy49yNu61w+gtVEF8io+PT2+L7S4/Vn1xg+fOexzPkCByPpHp2/i9fK47zV8q/QvUqOlF8v6dV2lKJn6h26aoCXIs6kvBAuUCAoDA90BgkNHkFRtGWHbzDpV+ICHI0mg/mAVRK3OfFdZynCCFiQr6Jlo9O3/bH2b+hzK+vRr9h8oDUUAAEAAEBi4BQWaE/VEbIXl9C/n+KKjfyvj2atRv6EBBgAAgAAh8TgJgwexz0gW6AQFAABAABAY8AeAIB3wTAQMBAUAAEAAEPicB4Ag/J12gGxAABAABQGDAEwCOcMA3ETAQEAAEAAFA4HMSAI7wc9IFugEBQAAQAAQGPAHgCAd8EwEDAQFAABAABD4nAUEen7h79y6VSuVjFf7qGQcHBz4CbEm01lZYWLjXHhltriquQpRVpQWpUk9t/W7lMGorFREmcrYSb+w8s3y3EEHFAQFAYGAS4OzZurcTv4ndnm/ABbrXwpBovbXDbdT8hCKOJ+d7mL9DDKtIWukx5cfrfG5o/A8q0Q+Jy9xcF50uZNqF3whyelug12jbERYWjiHXm6n/xs5y9gxJx2/O/apDl2ryqgtGvrLZzSk4sZyzuqzYWbXxycKrFBAPCAACgEAvCfSwT+MoRRBHyLyGCb+Yl3k3L66T9TP+ldfFvBylQ2gTuZ55SRJWm3Pu4I9/PKVwCvUkhveFFz3JzSFDbW5qxv8xjMHqb4evDbtOHu6/bsePO7dNNxXGqE0NTU2NLT2+rIlDP48ItCJjl7ed5/5stvk2Wpl1fP2MMQ629u7T1sbcK2cpVrCkjtJZq8nDJHo0rbGunvO9rww2GMsdYKzaeGfhUw5IAgQAgW+AQGveqaBRI2yWp3bMS7CG3LQDS6e42tnYjp4wd/vvT6o5R9X0/qQnYq23Q0ePGLMnq8NV0IrT1o+1dAyI+Qe/S461F6I+OxY43tnBGp/AjLCyH+3pF7Tpl4tvOC/DxUsWZB2xW6fTrQCXtsaqslNOp6kZrffnktjfUYj2rNiMGRCBQB8nUHJu3KpS9d63a6FF2/VHkPyKpLtLYQLhv1y11F0dWkuzEg7sjL75kapkyypLefPrmtVRH/S95y1XK7uRELNuDRx3IsiQbopgSZ3K2arZnX3dprNq68EVmd3qAwKAACDw9RHAyA9/CYnMJqOdL9JGyy+FLv4hW90nKDRYrubh6WMRy0sICcdmaLHPw3ooxo4Eq836ec2eDOHxe8MWGOF3NrF23WhtcWEpZhSwbeJgQktdeeHTGxd/2347uzw+dq5+l9d8C+IIBb+Y9ytqVYTQdoEk1lRV1YgoqquwoepI7pMqYXV39gdsOFci5+Bk3nCjmFUn+Vb8qX/Fx4X9vIV+Ue/kobSpKxJ+vem3bwx+h69ASewG9209+lZbn6AFSgABQKD/CGC1dw7uSILdPY2utF/PC6Fl15JvNVhuCNs0VQGfO7hYEIsm7Lt4o2jaXG1WT9hDMba6NL/5fcOGP8rN18VucVVoU8beC8Giug5e40cweu8Zsz2Pzpp9IjX99Rx9I7YLgiHBl0b5sOW7NNqUn35wua+bg42tk1fgtj/zWVcBW65ttDWnB6vpMblc585tpaKVD3/bPm/iaDsbu9Hj/EOSCtiEsbq7BwMmujpaW1o7uPut/Pn6e+aeJVqV/eum2Z4jbaxsR3lMXXw8hx7NLRKrOb/Mxnr+6VK6WgxFIcrjA+NG0A2zGLk1oxXNP+5vZbvuchPTGkrxregNMz0dbWxHeQVuP/28fT2g5VXy7pVzvN0dbSytHd1nRz3lecAIltCxnhQccfZMuN8Q9ht9W5/dyayTdvC0Z9xJD4lbjR0lX595m265YEmszcZWTTyBVvYgfmuAl5OtjePY6WsOXyls7pSmld4IW+ztYmtt5zx54e7EF3Uck76u2qDus7AaAz4DAoDAV04AIz+IPHBFOnDnEnNxluWyhoYGDBYVbbtjApaQIgnh+2ccPUgPxToYUd+nbl1x+IXmvIP7p7VftcvRC7ERJSoq4peVs+7ntCcLMiNkrnwyNwgtLCxwVayf8a+8l0bRskvbFm67JWQxKXCKnmjVq78vpKOQXoetRJPZO2eb4MJPnSAAACAASURBVDbBpMGqPNcdsboHPwWtPFuj5+kbbCwHVX+gqsiyCcOikhomHoET1CShssyzJ09u2Syp89tcndore9dGPhnsu2S7hTz2qaBEVlYIwqq5RHL+Ggm6k0MWO0jDEExUGS4E1XVKYNUZexauu0izmr5ss3bDwzPxYcvrxP53cJIygjW9upF6p9zCP2iGjgRaS9FiXATPPSBqLoGz8CRqOXs6WlFQWE9Q19VsbyeCpo4G0vCusAwdSRAkyUmTlxFY7b39C1clNwyfNGulPvHj7aQ/Ni0oaD4ZPpHZEFjZ22LjSXNWKWEf751P3LvkXevJKH9tnm1Erwb3LOwDMe44QCwgAAh8fQRa/ok/eEncPyZg6KALLNYjmiMslY8lnjz2wHatnRy18M+kTGz4PJcuC6NQD8XaFNNKr+3avfuOhE9YxGIzCT4dEUZpJJPJBGpTbcm/V2LiHhENF44ZwtELCeII+U746GbyFKC9Toy7SR62JOHofD3GGu3UweTRW0s7mCGKRs5ubt1cToUWJkYmfdAKOBG3wnhQe1aUzYMQjX3XGjOTvMzgghmxdzLLAzQrXr6ul7Cbs8rfqbME6gcukRDHWAWRHuLo6qbc5kTQTkeI5iVGp5UNnnvyUDB9226CtVjp1IjE9KIJc5lOgqA2auZsv/aMHRXt6QeMXEdGYWVSR0MjJEkJBKusJWMYQZAkngWjBYmR54u1A3+NXjacjsfXTXPZ9APRpx57bBxBzyRkGhS+Z4YKnYCfp84Kv92/nsz02c62m9lVN/cs9h1t1lUefAcEAIGvlwD6IeXw2Yax+wOGicBVbNUQNl+0Ozh39ZFl3vesjEVzs8kOO6Jn6nAMyXsoxlBdk74vtBk1WJ6w3lGOQxFr4WhZ0gqXpPYYhGQatHe6Qfthj05Bvjp4tAlzwofPBfHAPC/K+hnPxGtGiFbkPH4PD3Fx1e3FhcT4qmjmG0x/rNfw7ntUfFWToDlUTwwj19ZhBC1zU1nytUPb4zPe1rafvOQayaPeXKLR8uyst5Cuk7M2hD8z0tpKlTMx18DyX+W2cBHuwyiY9whIsCS04mFmHk51nGHbIIGg5jHeklie9eBtxyHV9jIRJVcPc2LV0ycF3R6bFSBLH1ICqgABQKC/CDQ8iP/t5fCABXZcpme02vd5JVQNZ99x+qIYRCt5kHjpHy6HN3soRq+RqJ6lqTz1RfyOIw8+8dtEg2BZ1/XRsbGxMcciw0KDx2sUxi6Zf+gRR+GCOEKeE7525LwEsLpaMopIyeBLjIIHrKaqGkVk5dhXQ9n00cruxYbMmeBib2VpZTv2h78bcN+ML7dKOoUc2TlV9WXcGl+PCYsOpOXhh225RvbcOKyuphalvonxt7dhBvtZx/OoVDK5kWNS2XOlLJKwhCQJwerJHe2G4p9RGJ8WwoIl8bICq62uYacKi8rJiWM1eDRHHlhMRmYQRq4hcyZxyLZFCJCFlyoQDwgAAgONAFp66eRlitMUF+nmRjw0UfAtQLS1qakV7yPQD+d+2HNLZckv+5cHrQ1LSDrir/A8ZseJf7ucmeihGLPmIgazDp8M91XJO7l2+eGHfB4fh4UVho6wtLS0snZwnRy0I2rbWNG35xIyup5xEHxplM8eIS9HCEvJSCHop4pKDOL0hTCCO2Vec0mWZodJUiQExTtobkrocujHxO3rjxUYz16+x05bWqgh81jI8WqGAlhcf/zaI55BeX8nHAqPD11OkU/eYivKJdKmpz8zhjFCQ2cc2DxGvtO7wyQ15tmWnqrhKYco6mhL0G7nF1EhdcZ0nlaU/x4Vt9JWQhBYkCReJTGbpupTdQdVrOnTpwY8WppzrITVV1Q0M3LwUscRL0AWDh0gAhAABAYogaYn93MaP7VsHHOJxcDdni53d18/aHX/5nOK4erRaoydOUTaavZU0zMHcnLK0eEs5yawTz0S61QvpOy04WgYErzu1Pq1MtGRAQbdLxFCsIS2jjJyraYan1qw9dGdXRnWWJr7tqInS3rdeiteAoiCpY0enHspOZtjZopbSJJA0Iryyu5mGYiimYUG9ubyxZcsZxohuhuFMQoFf+6bVvDydYvymEVLvUdbm5tZ2BirEtlmoIiknuuSvcH2xMrnzz+2l8Y1svvfHG6MuQb0/k2Z7DDjzmCkjR/D6ZtANHGwJtXc/es+cwTT8PDy7Upxa0cz3CsKloQ/s1r86k0ZKzvGr5PZNFfSX7X9AtCSK5eyKQpW1oM59pVb85JTHlOUrGx0CCzY+daXJQtfOZAICAACXyMBUdvlx+JOtIe4QwHDhYgm83+JXmEnCgmLiuLzn/LK9p0UWlUFvqgnLiHKvjLYQzFWOoic/dqIzc4iTyI3/HQbH8R3G1oL79wvQMXVNWS7DOKFHrZCVnivitXdPBC085bIhJ9OrxnRzQ6e4M8RIrq+y7wvrj6zZm75lAk2GmK0qqy3VEicbj4sY2qujcWc3BWGumjD5e9EXVf7DuPqTghDpwd7pG08sXTeRx8PM9VBzeXvRUetnKikIk67khKdaLDUzmCI8IVLP+9X9LZUFSe0vihivrsGLboQkVCibqwlS2wuyUrMoUjY6SphRRfCOSJhqNvNr3biuDGL3S9uigha9HqC0zAlUVp9WZGQ3YqZFj0YnXTbanQsJKfAmQa3onetFi0co1Z27dT56iFBgc6MwYwgSVjdtZ3+IdeI3pEXt9qxtjPeNEsnXVgdvzK4epq7HvHjncRzmSKO2wMs8R8HvoQM0YrvnDlNVBdpKXv6V+Kll5Kuu+aYC+PmdWBfNkW7S4W4Z+lJrYEMIAAIfFUEEBkdU5kOi7FPbyUQmKhpbKqDd1WY4yRXuXV/bN4gFDjOQKLuZfrJU2+lRi9wxJcFsYrExV77PgX8fibYQJK3GB8UBNXxO/bmzg3+fcceyz/2eyp0FcWa8u9eSP0I42dHK989uZ6e8ZZounw27p7Zg5Asc8APC8urKUlIiqlxWQnjUM3tOQxWIV4zQrxfl7LfGHtU48ixpLToW9UURFxa2dDeVp+E9+uI3qzQte9+jE+N/IkqrjLcx6oJg+jxnAGRcw2Ni9D6JTY1+ejVepikajjZrEnYdu7aSfmR1xPv+k1ZFLqnMjzmYuye5EZMhCQjr25irS0BY/WNJc9TUv8sr6cKS6oajlwVvtqZBOVxiYShGs5SecQg8u4/HB+kHRmXdi4qjUwbJKc5bLwpbjr744A8MvckmmgwPyIc2n/o7LGDDWJaNgvCQhYMazv0JEASLKKgoUKSElXHH6ehB/reKRMyLOWw8dghxUNRKb9F/IlKaltO27Nq8VjG0oWItuVI09IXKZFZZAqBpGJgt+jg0sBRjKOwIp3YfbQ7qwPzydKTSgMZQAAQ+FYIwNLOW6P3yEbEJYdtqqYOktex8N2xfLG7HN7xoDQKFUMkpeg9PR8xviRgcYvg0DmZ8+PCw69a72Xd1kIklTUVsWenfnxIw2AhMWklTQOnRUFzZrq1P3XI0mPV1HTT6WdnZ7u4uLBacvbsWW9vbz62paSk+Pn58REASQOCAJob4+9/QuWHqxHjSQPCIGAEIAAIfE8EqM8jfOZdcYhK3WjZzTrkZ6bCdfGxmzLxW5aSk5P5CBGJva4U1tpQ18j5omcYFhaXEhPEaD7m9k/SgKoRWvr46qP815cSc4kmk80Ya9MgAAKAACDQrwTQkic55fKjnI177TB6a7YgPsXHx6e3xXaXH6u+uMFz5z2O25wQOZ/I9G1sm1vdqRog6QOrRk0vkvfvvEpTMvEL3TZFreeHPwcITGAGIAAIfAsEBhlNXrFhhGU371Dph5rCAiyN9oNZELUy91lhLccJUpiooG+iJcl167A/zOpFGd9ejXoBA2QFBAABQGDgEBBkRtgf1gvJ61vI90dB/VbGt1ejfkMHCgIEAAFA4HMSAKtin5Mu0A0IAAKAACAw4AkARzjgmwgYCAgAAoAAIPA5CQBH+DnpAt2AACAACAACA54AcIQDvomAgYAAIAAIAAKfkwBwhJ+TLtANCAACgAAgMOAJAEc44JsIGAgIAAKAACDwOQkI8vjE3bt3qdQud0mx2Yi/esbBwaGnZtNaW2Fh4V57ZLS5qrgKUVaVFqRKPbX1u5XDqK1URJjI2Uq8sfPM8t1CBBUHBACBgUmAs2fr3k78JnZ7vgEX6F4LQ6L11g63UfMTijienO9h/g4xrCJppceUH6/zuaHxP6hEPyQuc3NddLqQaRdW9/z0tkCv0bYjLCwcQ643U/+NneXsGZJe0Wuz/4NNn0G0SzV5lYCRr2x2cwpOLOesLit2Vm18svAqBcQDAoAAINBLAj3s0zhKEcQRMq9hwi/mZd7Ni+tk/Yx/5XUxL0fpENpErmdekoTV5pw7+OMfTymcQj2J4X3hRU9yc8hQm5ua8X8MY7D62+Frw66Th/uv2/Hjzm3TTYUxalNDU1NjS48va+LQzxFBffzT+BHmncFqdnybF8ZxVmYdXz9jjIOtvfu0tTH3ylmKFSypo3TWanKY1BlBa6yr53zvK4MN60UkrNp4Z+FTDkgCBACBr4IAWpGxy9vOc38268og3hvFhfh7ONrajho/Z0v8Q467ZVtLH56L2LTQ19N59NLEii7XB2K1qcttWLpA5kfbdZebuhBpvR06esSYPVkdroJWnLZ+rKVjQMw/+IVxrL0Q9dmxwPHODtb4BGaElf1oT7+gTb9cfMN5GS5egCDriN06nW4FuLQ1VpWdcjpNzWi9P5fE/o5CtGfFZsyACAT6OIGSc+NWlar3vl0LLdquP4LkVyTdXQoTCH33qjcMd7sEfZ9N86yZV0/BUkMVmNopb35dszrqg773vOVqZTcSYtatgeNOBBnSTREsqRMmWzV7zZhVWw+uyOx1eUABIAAIfAECraVZCQd2Rt/8SFWyZSm+9WXcqlXHyk3mrNk5FH15/vjRlauoJ04EGbR1mmjFrbCVW87kiRi6jJu21MDAsuur/mEx81nbdrh3LDuhxdejYh9I6WrwfyE3Vpv185o9GcLj94YtMBKDINauG60tLizFjAK2TRxMaKkrL3x64+Jv229nl8fHztXvolUQRyj4xbxfoNUELRIhtN3MjjVVVTUiiuoqbKg6kgXV3yUf1tjQiGhae3mM7XKjL/lW/Kl/xceF/bxlNH7F5eShtKkrEn696bdvjBQsWBJ7wX1bj77V1kdogRpAABDoIwJY3Z39ARvOlcg5OJk33Chm0dp4P+GPV6QJEeErRuIuboyjat3UFad/v+//oxPd41HzE0I2nykeuvDYwSBzKe7rkEQNGy+Ndo1o6flV4dVa0/fP4349e5tc85vfN2z4o9x8XewWV4U2tey9ECyq6+A1fgSj954x2/PorNknUtNfz9E3auvf2xRxN4k/tW5XPvkKNOWnH1zu6+ZgY+vkFbjtz3zWuXXLtY22jBmx1fSYXM4dqU6z0MqHv22fN3G0nY3d6HH+IUkFbMJY3d2DARNdHa0trR3c/Vb+fP09c88Srcr+ddNsz5E2VrajPKYuPp5Dj+YWidWcX2ZjPf90KV0thqIQ5fGBcYyVS4uRWzNa0fzj/ladc3ZK8a3oDTM9HW1sR3kFbj/9vH2fsuVV8u6Vc7zdHW0srR3dZ0c95X3ACCPXkTFxSRJ74+CFtz67k1kn7eBpz7iTHhK3GjtKvj7zNt1ywZJYm5atmngCrexB/NYALydbG8ex09ccvlLY3ClNK70Rttjbxdbaznnywt2JL+o4Jn1dtUHdZ2E1BnwGBACBgU4AltCxnhQccfZMuN8QtrvH0Y+v8+oJQ8xM8EkZPYhbOFqI1T3NzqN3elhdRmx8Dmy7fDdPL9il4o1ZcTH3iW5L55l2mRiwylHfp25dcfiF5ryD+6e1X7XL0QuxKSYqKuI3knO7WF6QGSFz5ZO5QWhhYYGXxPqZXm9uJTEMQssubVu47ZaQxaTAKXqiVa/+vpCOQnodthJNZu+cbYLbBJMGq/Jcd8TqHvwUtPJsjZ6nb7CxHFT9gaoiyyYMi0pqmHgETlCThMoyz548uWWzpM5vc3Vqr+xdG/lksO+S7Rby2KeCEllZIQir5hLJxo7xhaA7OWSxgzQMwUSV4UJQXacEVp2xZ+G6izSr6cs2azc8PBMftrxO7H8HJykjWNOrG6l3yi38g2boSKC1FC3Gbe/cA0auobvP7It/UTSHmBjryrRN3NGKgsJ6grquZns7ETR1NJCGd4Vl6EiCIElOmryMwGrv7V+4Krlh+KRZK/WJH28n/bFpQUHzyfCJzIbAyt4WG0+as0oJ+3jvfOLeJe9aT0b5a/NsI3o1uWfh8PXciYBYQAAQGHgEEDWXwFm4WdTyLraJi4vBaE11DQYxx+xC8vLS0Kdyxj5hw4PLd2olnXw9lGGUhiKMHSd+AS1LP5VeqRcwl74KxiPQSq/t2r37joRPWMRiMwmeYng3RGkkk8kEalNtyb9XYuIeEQ0XjhnC0QsJ4gj5TvjoZvMUoL1OjLtJHrYk4eh8PUZXP3UwefTW0o6qIopGzm5u3VxOhRYmRiZ90Ao4EbfCuGO4gLK1C9HYd60xU6uXGVwwI/ZOZnmAZsXL1/USdnNW+Tt1lkD9wCUS4pjtINJDHF3dlNsaEO10hGheYnRa2eC5Jw8F07ftJliLlU6NSEwvmjCX6SQIaqNmzvZrz8ijTfHGaiZIa4rnJ+3bEt+CERUsZ/+wJ9heHoHoM0UUViZ1NDRCkpRAsMpaMoYRBEniaQFakBh5vlg78NfoZcPpeHzdNJdNPxB96rHHxhH0TEKmQeF7ZqjQCfh56qzw2/3ryUyf7ax7BByquWex5zPE41ABIgABQOBrIIAo240c+suRMz//z2TtOD3RxtJXWXkN+IyI3pXSCl6+wU+8vIiZ5bqpqLpVWN7Qedaa9bMsZHg4RPRdWvJDzCpk0mAOf9WBoiZ9X2gzarA8Yb2jHA81TFm0LGmFS1J7PoRkGrR3evu+JQtYvjp4NACzdvhcEA/M86Ksn/FMvGaEaEXO4/fwEBdXXf77nzzKbatY5cPMN5j+WK/h3feo+KomQXOonhhGrq3DCFrmprLka4e2x2e8rW0/eck1km/5bIloeXbWW0jXyVkbwp8ZaW2lypmYa2D5r3Jbeq6DLkkwCDhyLjX95p2bKVHrXMSfx4eEnmesy/IIMO8RkGBJaMXDzDyc6jjDtkECQc1jvCWxPOvB245Dqu1lIkquHubEqqdPCro9NitAFh41BtGAACAwgAkgujN3bh4n8TB8gZezk6uX/4b4x9WQlAx9P5BWXVmNiSgMcwvcHB4VFbbaCbrz86qtKcU8Ojjam/T0NyI2413aN/241VpUz9JUnvoifseRB5946GFmg2Vd10fHxsbGHIsMCw0er1EYu2T+oUccJ0f7dUaI1dWSUURRBl9iFDxgNVXVuBI59tVQNn20snsnIqJTs/JKalqEpKSEGiBVfGACSzqFHNkpfehY3BrfGCUr7+D1y7z0xLhG9tg6rK6mFqXmxvjbx7DkEVYgN2JQ936asxghkpat/47Q0n/nn/nrVsVkPwUJSRKC1ZPxdmNCQ/HPKIxPC2EYESSJs0hmDFZbXcNOFRaVkxPHnuPRkEqXXLCYjMwg7GUNGYXkeCnknQXiPczroTIgBggAAgOOgLDOxJ1nx658V1DaRJRVrD49d9FZ9SG6+F87ig/OYSnLKXMmWdD9jbWVFvnFtGOXMsq8mUtM7DWh5d24+U7Eaq4D72VRPIOIwazDB6b8vHLTybXLsUNRK614ScPCCkNHWDIPy0AOLi661MkbziVkzBsxni2H4I6Qzx4hr6VRmD5AQD9VVLZ366wEYAQfO/CaS7IIwiQpEsJcjObuUNGPidvXHyswnr18j522tFBD5rGQ49UMBbC4/vi1RzyD8v5OOBQeH7qcIp+8xVaUS6QNe9Pw/sYwRmjojAObx8h3eneYpMarWXir6kwhaumqC9EqqvBdQyVFHW0J2u38IiqkzjiGTCvKf4+KW2krIQgsSBKv8plNU/WpuqNpsKZPnxrwaGnORQOsvqKimZGDlzqOeAGycOgAEYAAIDCwCcAictoG+Oi44d6P1z4Km860k8G7RYKcgixUX1ODLyAx/A2ioq5KQAur8Q6u6xgb95olWVmFBCM/a8nupktCyk4bjoYhwetOrV8rEx0ZYNCDqQcsoa2jjFyrqcanFmx9dM+7ss4G6NZb8RJAFCxt9ODcS8nZHDNTCJYgSSBoBXNzlV9jI4pmFhrYm8sXX7KcaYTobhTGKBT8uW98Rfp1i/KYRUu9R1ubm1nYGKsS2ZAiknquS/YG2xMrnz//2D6r5hrJzwxmGm6MuQb0/k2Z7DDjzmCkjR/DEThg5Oc5b6kiahpKeOMQTRysSTV3/7rPPKXZ8PDy7Upxa0cz3CsKlgRhDcWv3pSxsmP8NplNcyX9VduiLlpy5VI2RcHKmnOdvjUvOeUxRcnKRofAgp1vfVmy8JUDiYAAIPC1E2h9f2nP3j8rNLznezEOFRB0LMykm3PuP8Gfd6eHlpfPXlMHaWkzEruGxn9y3mAaJsZ8Vvw6syBy9msjNjuLPInc8NNtfBDfbWgtvHO/ABVX15DtUrYgHbbgzxEiur7LvC+uPrNmbvmUCTYaYrSqrLdUiPFkJSxjaq6NxZzcFYa6aMPl70RdV/tyf4KEMHR6sEfaxhNL53308TBTHdRc/l501MqJSiritCsp0YkGS+0MhghfuPTzfkVvS1VxQuuLIua7a9CiCxEJJerGWrLE5pKsxByKhJ2uElZ0IZwjEoa63fxqJ44bs9j94qaIoEWvJzgNUxKl1ZcVCdmtmGnRg9FJR6uhpSnbQx/JWRtp4A8HvstKS8mo1vKfyTwxRXIKnGlwK3rXatHCMWpl106drx4SFOjMGMwIkoTVXdvpH3KN6B15casd604t3jRLJ11YHb8yuHqaux7x453Ec5kijtsDLHGXS//90orvnDlNVBdpKXv6V+Kll5Kuu+aYC+OnezuwL5ui3eVnyD1LFyHwFRAABL4FAlhd3r17/7zLf3b7r8tZpVJOGw6utG57aF7Uxn/msCtRO9ZLBk0yRN6mHT9XrjHtB3oHh1UkLvba9yng9zPBBowtE1rBq7xWEcOhWj3dQCGojt+xN3du8O879lj+sd9ToStKrCn/7oXUjzB+drTy3ZPr6RlviabLZ9uJdpETxBHymvB1qOYtAEvZb4w9qnHkWFJa9K1qCiIurWxob6tPf5sKojcrdO27H+NTI3+iiqsM97FqwiDmW1a6Vg2Rcw2Ni9D6JTY1+ejVepikajjZrEnYdu7aSfmR1xPv+k1ZFLqnMjzmYuye5EZMhCQjr25irS0BY/WNJc9TUv8sr6cKS6oajlwVvtqZBOVxiYShmq5l8vyOyLv/cHyQdmRc2rmoNDJtkJzmsPGmuOlsD9nwzN2WgEkriL+/cfpmSQ1VWFbDyH1VSLB/+5sXiAbzI8Kh/YfOHjvYIKZlsyAsZMGwtpc1CJAEiyhoqJCkRNXxx2nogb53ypwvw1IOG48dUjwUlfJbxJ+opLbltD2rFo9lPPIhom050rT0RUpkFplCIKkY2C06uDRwFOMorEgndh/tzmrCfLJ0SoFPgAAg8K0QQD9e+XnnH1UKusOt5++f5uuiT+qcdAnpBxw8BIUfPhex7RQqpWMT+NP6hSPoTxxiNAoVQySlOnp6SmlxBSTnrPgfHBMsbhEcOidzflx4+FXrvazbWoiksqYi9uzUjw9pGCwkJq2kaeC0KGjOTLf2pw5Zeqyamm46/ezsbBcXF9bWOnv2rLe3N5/2S0lJ8fPz4yMAkgYEARQ/5eN/QuWHqxHjSQPCIGAEIAAIfE8EqM8jfOZdcYhK3WjZiycJ+oDYf3C8HaXhtywlJyfzKZxI7HWlsNaGukbOFz3DsLC4lJggRvMxt3+SBlSN0NLHVx/lv76UmEs0mWzW9a1//QMElAIIAALfNwG05ElOufwoZ+NeO4zechTEp/j4+PS22O7yY9UXN3juvMdxmxMi5xOZvo1tc6s7VQMkfWDVqOlF8v6dV2lKJn6h26aocduzHiDYgBmAACDw7RIYZDR5xYYRlt28Q6Uf6g8LsDTaD2ZB1MrcZ4W1HE9KwkQFfROtbg/W9oeF/7WMb69G/5UAkAcEAAFAYEASEGRG2B8VEZLXt5Dvj4L6rYxvr0b9hg4UBAgAAoDA5yQAVsU+J12gGxAABAABQGDAEwCOcMA3ETAQEAAEAAFA4HMSAI7wc9IFugEBQAAQAAQGPAHgCAd8EwEDAQFAABAABD4nAeAIPyddoBsQAAQAAUBgwBMAjnDANxEwEBAABAABQOBzEhDk8Ym7d+9SqVQ+VuGvnnFwcOAjwJZEa22FhYV77ZHR5qriKkRZVVqQKvXUViDHSQCjtlIRYSJnA/JuEZ5ZOLWDGEAAEAAEPjcBzu6r+xLxm9jt+QZcoHstDInWWzvcRs1PKOJ4cr6H+TvEsIqklR5TfryOX3HVBwH9kLjMzXXR6UKmXVjd89PbAr1G246wsHAMud5M/Td2lrNnSHpFr83uA1s/j4ouBHgVgpGvbHZzCk4s5yTB2iKs2vhk4VUKiAcEAAFAoCcEethxcagSxBEyr2HCL+Zl3s2L62T9jH/ldTEvR+kQ2kSuZ16ShNXmnDv44x9PKZxCPYnhfeFFT3JzyFCbm5rxfwxjsPrb4WvDrpOH+6/b8ePObdNNhTFqU0NTU2NLjy9r4tDPIwKtyNjlbee5P5ttvo1WZh1fP2OMg629+7S1MffKWYrt86QOw1gJ8LCWHk1rrKvnfCUsAxvGMiZh1cY7C59yQBIgAAh8eQLNRRlxPyya4uZgY+3gNiV4f2puQ9epR2vpw3MRmxb6ejqPXppY0Z7Kp6diqRX18U/jR5h3BqvZ8W2TERah1tuho0eMdv5sVAAAD+dJREFU2ZPV4SpoxWnrx1o6BsT8g18Yx9rVUJ8dCxzv7GCNT2BGWNmP9vQL2vTLxTecl+HiygVZR+zW6XQrwKVBsarslNNpakbr/bkk9ncUoj0rNmMGRCDQxwmUnBu3qlS99+1aaNF2/REkvyLp7lKYQOjuDuX/YndraVbCgZ3RNz9SlWxZ81He/LpmddQHfe95y9XKbiTErFsDx50IMqSb0udJneWyEfgv1eAqy6qt698N1wwgEhAABAYeAYxamHH5Dcl+xnL1QdUvLyee3bGiVT5xm337a/vRilthK7ecyRMxdBk3bamBQftdcnx6KrZKYvjsg6Dvs2meNfNeJlhqqEI3nSxWm/Xzmj0ZwuP3hi0wwm92Yu260driwlLMKGDbxMGElrrywqc3Lv62/XZ2eXzsXP0ur/kWxBEKfjHvwGtZnhYhhLabIbGmqqpGRFFdhQ1VRzJPBf8lAau7sz9gw7kSOQcn84YbxaxZybfiT/0rPi7s5y30eywnD6VNXZHw602/fWPwO3z7Oond5L6tYt9q+y9wgSwgAAj0DQFYYtTW/41EEIZzwibqNU8Iybj+LMSeeQ0CNT8hZPOZ4qELjx0MMpdiWWzk01Ox24U1NjQimtZeHmN7eLF585vfN2z4o9x8XewWV4W2Etm7GlhU18Fr/AhG7z1jtufRWbNPpKa/nqNvxH7zr+BLo3zI8l0abcpPP7jcF59b2zp5BW77M591FbDl2kZbxrzYanpMLue2U2eRaOXD37bPmzjazsZu9Dj/kKQCNmGs7u7BgImujtaW1g7ufit/vv6euWeJVmX/umm250gbK9tRHlMXH8+hR3OLxGrOL7Oxnn+6lK4WQ1GI8vjAOMaU3WLk1oxWNP+4v5XtustNTIMoxbeiN8z0dLSxHeUVuP308/Z9ypZXybtXzvF2d7SxtHZ0nx31lOcBI1hCx3pScMTZM+F+Q9hv9G19diezTtrB055xJz0kbjV2lHx95m265X2exNqibATwBFrZg/itAV5OtjaOY6evOXylsLlTmlZ6I2yxt4uttZ3z5IW7E1/UcUz6umqDus/Cagz4DAgAAgOEQJsXxK2BRcTFhGAYafMhWF1GbHwObLt8N7sX5NtTsVcKI9eRMXFJUg9vp6e+T9264vALzXkH909rv2qXo6thK4GoqIjfSM66adOeLMiMkLnyydwgtLCwwFWxfsa/8l4aRcsubVu47ZaQxaTAKXqiVa/+vpCOQnodthJNZu+cbYLbBJMGq/KcEmN1D34KWnm2Rs/TN9hYDqr+QFWRZROGRSU1TDwCJ6hJQmWZZ0+e3LJZUue3uTq1V/aujXwy2HfJdgt57FNBiaysEIRVc4nk/MkRdCeHLHaQhiGYqDJcCKrrlMCqM/YsXHeRZjV92Wbthodn4sOW14n97+AkZQRrenUj9U65hX/QDB0JtJaixbjtnXtA1FwCZ+FJ1HL2dLSioLCeoK6r2d5OBE0dDaThXWEZOpLQx0lOmrzsw2rv7V+4Krlh+KRZK/WJH28n/bFpQUHzyfCJzDbCyt4WG0+as0oJ+3jvfOLeJe9aT0b5a/NsPnoNuWfp4e+fO0IQCwgAAv1CAKM219eW5d4+efgKxWT+FHPmKmPDg8t3aiWdfD2UYZSGIoxtJUbg04l16XMwcg19FpF98S+K5hATY10ZPtcU0kqv7dq9+46ET1jEYjMJPr0NRmkkk8kEalNtyb9XYuIeEQ0XjhnC0dUI4gj5TvgYNcfnUFwD7XVi3E3ysCUJR+frMeo4dTB59NbSDllE0cjZza2by6nQwsTIpA9aASfiVhh3TKBRNg9CNPZda8zU6mUGF8yIvZNZHqBZ8fJ1vYTdnFX+Tp0lUD9wiYQ4pjSI9BBHVzfltrZFOx0hmpcYnVY2eO7JQ8H0bbsJ1mKlUyMS04smzGV6AoLaqJmz/dozcoXCL5I+REJhZVJHQyMkSQkEq6wlYxihj5N42oEWJEaeL9YO/DV62XA6OV83zWXTD0SfeuyxcQQ9k5BpUPieGSp0OH6eOiv8dv96MtNnO9tGZ1fd3LPY93A9pKs28B0QAAT6jQD67tQC38MvqBAia78ueroB8+QEreDlG3yJ7EXMLNdNRdWtwvKGzrPWrJ9lIYPw6cS62Iw1E6Q1xfOT9m2Jb8GICpazf9gTbC/PdXxek74vtBk1WJ6w3lGOq0C7arQsaYVLUvs3hGQatLfdZNbS+erggZY54cPngnhgnhdl/Yxn4jUjRCtyHr+Hh7i46vLx9DwK7YjGV0Uz32D6Y72Gd99t4quaBM2hemIYubYOI2iZm8qSrx3aHp/xtrb95CXXyO5M6ExHy7Oz3kK6Ts7aEP7MSGsrVc7EXAPLf5Xb0nMdgkjCvEdAfZ6EVjzMzMOBjzNsGz8Q1DzGWxLLsx687Ti/2m4OouTqYU6sevqkoNsTtQJkEYQUyAMIAAJ9SgBRGbflyOED2xa5yz0LW7DsZC5j44lWXVmNiSgMcwvcHB4VFbbaCbrz86qtKcU85kR4Ds6eimAQcORcavrNOzdTota5iD+PDwk9z9ie4gyiepam8tQX8TuOPPjEXaItDyzruj46NjY25lhkWGjweI3C2CXzDz3iODkqiCPkOeFrN5eXAFZXS0YRKRl8iVHwgNVUVaOIrBz7aiibPlrZvdiQORNc7K0srWzH/vB3A+6bce6STiFHdk5VfRm3xtdjwqIDaXn4YVuukT03DqurqUWpb2L87W2YwX7W8TwqlUxu5JhU9lwpiyQsIUlCsHpyR7uh+GcUxqeFcJ8n8TIQq62uYQcOi8rJiWM1eDRHHlhMRmYQvsJB5kzikG2LECALL1UgHhAABD47ARElQ+uR7j6L98bu8BTOORF3u55eJO7XYCnLKXMmudpb2zj7bdoZMLTp4aWMMpRPT8XDVCGSlq3/jlBf9YZHf93i/rS2iMGswyfDfVXyTq5dfvghn8fHYWGFoSMsLS2trB1cJwftiNo2VvTtuYSMrgcZBF8a5bNHyMsRwlIyUgj6qaISgzh9IWPTlddckgUYTJIiISjeC3NTQpdDPyZuX3+swHj28j122tJCDZnHQo5XMxTA4vrj1x7xDMr7O+FQeHzocop88hZbUS6RNjwaiCOaYYzQ0BkHNo+R7/TuMEmNebaFQ/6/RiCKOtoStNv5RVRInbECQSvKf4+KW2kr4XvWfZzEyzZmq1V9qu4AjjV9+tSAR0tzDqOw+oqKZkYOXuo44gXIwqEDRAACgEB/E4BJw4y0CZc/FFWikARBTkEWqq+pwZeCGE4FUVFXJaCF1bUYosuzp+JnMVFLV12IVlGFOzklrnJCyk4bjoYhwetOrV8rEx0ZYND9EiEES2jrKCPXaqrxqQVbH93z/qrTlm69FS8BRMHSRg/OvZSczTEzxS0kSSBoRTnOlH9AFM0sNLA3ly++ZDm4CNHdKIxRKPjD3fhi9esW5TGLlnqPtjY3s7AxViWyzUARST3XJXuD7YmVz59/bC+NayR/QxitrWhmrgG9f1MmO8y4Mxhp48dw+iYQTRysSTV3/7rPHME0PLx8u1Lc2tEM94p9ngRhDcWv3pSxYmXUkdlqV9Jfta33oiVXLmVTFKysB3NsObfmJac8pihZ2egQWFqELwqWLHzlQCIgAAgMKAL4G7dy8mkiyir0XTqCjoWZdHPO/Sf4Ohs9tLx89po6SEsbPzvAp6fiXR+M/DznLVVETUOJj49C5OzXRmx2FnkSueGn2/hIvdvQWnjnfgEqrq4h20WrIB224M8RIrq+y7wvrj6zZm75lAk2GmK0qqy3VIjxNCYsY2qujcWc3BWGumjD5e9EXVf7DuNqHWHo9GCPtI0nls776ONhpjqoufy96KiVE5VUxGlXUqITDZbaGQwRvnDp5/2K3paq4oTWF0XMd9egRRciEkrUjbVkic0lWYk5FAk7XSWs6EI4RyQMdbvD1U4cN2ax+8VNEUGLXk9wGqYkSqsvKxKyWzHTogejk25bjY6F5BQ40+BW9K7VooVj1MqunTpfPSQo0JkxmOnjJKzu2k7/kGtE78iLW5nPBbUZiLfa0kkXVsevDK6e5q5H/Hgn8VymiOP2AEvcG9N/9bTiO2dOE9VFWsqe/pV46aWk66455sK45R0tsmyKdpe6cs/SEyBABhAABL4UAUr20VW/Vw8301MUx2oLsy6dz2gYPGfmKBLdHlEb/5nDrkTtWC8ZNMkQeZt2/Fy5xrQf6E8/8+mpsIrExV77PgX8fibYAC5N2R76SM7aSAN/RvpdVlpKRrWW/0yGAj6BoDp+x97cucG/79hj+cd+T4WuolhT/t0LqR9h/Oxo5bsn19Mz3hJNl8+2E+0ix9XVdNXV5TuvCV+HGG8BWMp+Y+xRjSPHktKib1VTEHFpZUN7W336awQQvVmha9/9GJ8a+RNVXGW4j1UTBjFfL9DVHETONTQuQuuX2NTko1frYZKq4WSzJmHbuWsn5UdeT7zrN2VR6J7K8JiLsXuSGzERkoy8uom1tgSM1TeWPE9J/bO8niosqWo4clX4amcSlMclEoZqupbJ8zsi7/7D8UHakXFp56LSyLRBcprDxpviprM/Dsgze/cJRIP5EeHQ/kNnjx1sENOyWRAWsmBY2xtu+jYJFlHQUCFJiarjT9rQA31blfkbhKUcNh47pHgoKuW3iD9RSW3LaXtWLR7LeBpERNtypGnpi5TILDKFQFIxsFt0cGngKMYpWZHOFvHR7qwnzCdL9zSABCAACHwxAhQRWfmmjLT4CxUNkJi8jolXyPZFU0zx97nQg5B+wMFDUPjhcxHbTqFSOjaBP61fOIKZxrOnwmgUKoZISjF6ekxaQfz9jdM3S2qowrIaRu6rQoL9299Nw6fKsLhFcOiczPlx4eFXrfeybmshksqaitizUz8+pGGwkJi0kqaB06KgOTPd2p867NQK19R00+lnZ2e7uLiw2nH27Flvb28+lqWkpPj5+fERAEkDnQCaG+Pvf0Llh6sR4xnDPRAAAUAAEOhrAtTnET7zrjhEpW607MWTBH1glSAzQvyWpeTkZD6FE4m9rhTW2lDXyPk2ZxgWFpcSE8RoPub2T9LXUiO09PHVR/mvLyXmEk0mm7W/RLB/GIFSAAFA4DsigJY8ySmXH+Vs3GuH0VtogvgUHx+f3hbbXX6s+uIGz533OG5zQuR8ItO3se1gdadqgKR/NTVqepG8f+dVmpKJX+i2KWp8NqoHCFhgBiAACHytBAYZTV6xYYRlN+9Q6YfaCbI02g9mQdTK3GeFtRwnSGGigr6JliT/7dP+sO+/l/Ht1ei/MwA5AAFAABAYgAQEmRH2RzWE5PUt5PujoH4r49urUb+hAwUBAoAAIPA5CYClr89JF+gGBAABQAAQGPAEgCMc8E0EDAQEAAFAABD4nASAI/ycdIFuQAAQAAQAgQFP4P9XJXFcRDtDwwAAAABJRU5ErkJggg=="
    }
   },
   "cell_type": "markdown",
   "id": "talented-orchestra",
   "metadata": {},
   "source": [
    "You can check that the sizes you have computed match those given by your operating system.\n",
    "\n",
    "For example, here's a screenshot of my file manager (for the folder in which I am running the notebook):\n",
    "\n",
    "![image.png](attachment:567a50c8-7e3b-4e45-8910-bde5e19de09d.png)\n",
    "\n",
    "Note that my file system gives file sizes in KiB and MiB (1 KiB = $2^{10}$ B = 1024 B) and not KB and MB (1 KB = $10^3$ B = 1000 B). "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eea72112-7d8d-4ae4-a36b-1854b0b93b34",
   "metadata": {},
   "source": [
    "__Comment:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "polish-leeds",
   "metadata": {},
   "source": [
    "You can now remove the joblib files with the following command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "smart-charger",
   "metadata": {},
   "outputs": [],
   "source": [
    "for ssize in subsample_sizes:\n",
    "    file_name = \"dt_classifier_%d.joblib\" % ssize\n",
    "    os.remove(file_name)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "great-organizer",
   "metadata": {},
   "source": [
    "#### Limiting the depth of the tree"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "yellow-sacrifice",
   "metadata": {},
   "source": [
    "##### __Question 4__ \n",
    "Using the `max_depth` parameter of `DecisionTreeClassifier` to limit the depth of the learned decision tree, evaluate again:\n",
    "* the runtime, as a function of the number of samples, of training a decision tree of maximum depth 4 \n",
    "* the space taken in memory by such a model, as a function of the number of samples.\n",
    "\n",
    "Compare with the previous section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e2a218e-f490-45b4-b373-8f2ae92edbff",
   "metadata": {},
   "outputs": [],
   "source": [
    "### TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b91266ca-893c-456f-8bdf-28f0ef6c4a6d",
   "metadata": {},
   "source": [
    "__Comment:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "exotic-cover",
   "metadata": {},
   "source": [
    "#### Impact of limiting the depth of the tree on tree performance."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "accurate-wheel",
   "metadata": {},
   "source": [
    "##### __Question 5__ \n",
    "Using [`sklearn.model_selection.cross_val_score`](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.cross_val_score.html?highlight=cross_val_score#sklearn.model_selection.cross_val_score), complete the code below to compute and compare the 5-fold cross-validated [balanced accuracy](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.balanced_accuracy_score.html#sklearn.metrics.balanced_accuracy_score) of all these decision trees (with/without limiting depth, for all sample sizes). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "stupid-emperor",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import model_selection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "solid-parliament",
   "metadata": {},
   "outputs": [],
   "source": [
    "cv_scores_default = {} # key: number of subsamples, value: 5-fold cv balanced accuracy\n",
    "cv_scores_depth4 = {} # key: number of subsamples, value: 5-fold cv balanced accuracy\n",
    "\n",
    "for ssize in subsample_sizes:\n",
    "    X_current, y_current = subsample_data[ssize]\n",
    "    dt_classifier = tree.DecisionTreeClassifier() # default parameters\n",
    "    cv_scores_default[ssize] = # TO COMPLETE    \n",
    "    \n",
    "    dt_classifier = tree.DecisionTreeClassifier(max_depth=4) # default parameters\n",
    "    cv_scores_depth4[ssize] = # TO COMPLETE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "forward-spain",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.errorbar(subsample_sizes, [np.mean(cv_scores_default[ssize]) for ssize in subsample_sizes], \n",
    "             yerr=[np.std(cv_scores_default[ssize]) for ssize in subsample_sizes], label='Default')\n",
    "plt.errorbar(subsample_sizes, [np.mean(cv_scores_depth4[ssize]) for ssize in subsample_sizes],\n",
    "             yerr=[np.std(cv_scores_depth4[ssize]) for ssize in subsample_sizes], label='Max depth = 4')\n",
    "\n",
    "plt.xlabel(\"Number of samples (log scale)\")\n",
    "plt.ylabel(\"Balanced accuracy\")\n",
    "plt.title(\"Cross-validated performance of various decision trees\")\n",
    "\n",
    "plt.xscale(\"log\")\n",
    "\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f843a126-a722-4474-b164-e5bd3361b280",
   "metadata": {},
   "source": [
    "__Comment:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "incomplete-virginia",
   "metadata": {
    "tags": []
   },
   "source": [
    "##### __Question 6__ \n",
    "* How does limiting the depth of the tree affect performance? \n",
    "* Was this expected? \n",
    "* Comment this result knowing that our ultimate goal is to train random forests."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89aedb33-f87a-4765-b8d4-05a23d894c72",
   "metadata": {},
   "outputs": [],
   "source": [
    "### TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5eea9e7-e141-4bc2-bfe7-c0ea9e8b281f",
   "metadata": {},
   "source": [
    "__Comment:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "vanilla-representation",
   "metadata": {},
   "source": [
    "### 1.3 Computational complexity of a decision tree (with respect to the number of features)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "starting-produce",
   "metadata": {},
   "source": [
    "Here we'll work with the data set that has 10,000 samples (to avoid longer runtimes on the full data set) and change the number of features. If 10,000 samples takes too long to process on your data set, use a smaller dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adult-northern",
   "metadata": {},
   "source": [
    "#### Creating subsets of varying numbers of features"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "false-amber",
   "metadata": {},
   "source": [
    "Here we will merely add random binary features to the data. They will not be useful, but the algorithm will still have to explore them to know that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "collected-vanilla",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_subsamples = int(1e4)\n",
    "X_current, y_current = subsample_data[n_subsamples]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "under-teens",
   "metadata": {},
   "source": [
    "Let us now define a list of feature numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "honest-regard",
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_numbers = [n_features, 100, 500, 1000, 1500] # You can change the values here \n",
    "feature_numbers = np.array(feature_numbers)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "general-france",
   "metadata": {},
   "source": [
    "Now we create the copies of the dataset \"augmented\" with bogus binary features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "linear-essence",
   "metadata": {},
   "outputs": [],
   "source": [
    "augmented_data = {} #key: number of features, value: X_augmented\n",
    "augmented_data[n_features] = X_current\n",
    "for fnum in feature_numbers[1:]:\n",
    "    # Create an array of random binary features\n",
    "    X_aug = np.random.randint(2, size=(n_subsamples, fnum))\n",
    "    # Fill the first 54 columns with the true data\n",
    "    X_aug[:, :n_features] = X_current        \n",
    "    augmented_data[fnum] = X_aug"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "mobile-artist",
   "metadata": {},
   "source": [
    "#### Size of each data set in memory"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "starting-aside",
   "metadata": {},
   "source": [
    "##### __Question 7__ \n",
    "Compute the sizes of each data set in bytes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2409d5ba-6304-48ee-8b5c-49493fabaf5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "### TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ed5f5b7-c241-4c4c-ab87-461df9b5e866",
   "metadata": {},
   "source": [
    "__Comment:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "negative-jewel",
   "metadata": {},
   "source": [
    "#### Tree complexity as a function of the number of features"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "israeli-witness",
   "metadata": {},
   "source": [
    "##### __Question 8__ \n",
    "Compute and plot\n",
    "* the runtime, as a function of the number of features, of training a decision tree with no limit on the depth\n",
    "* the space taken in memory by the model, as a function of the number of features\n",
    "* the runtime, as a function of the number of features, of training a decision tree of maximum depth 4 \n",
    "* the space taken in memory by the model, as a function of the number of features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9966025-abb7-43f4-a626-3bda70609167",
   "metadata": {},
   "outputs": [],
   "source": [
    "### TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f197006-b6b8-4956-87c6-39c8b3b180b3",
   "metadata": {},
   "source": [
    "__Comment:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "instant-boating",
   "metadata": {},
   "source": [
    "##### __Question 9__ \n",
    "* What is the relationship (logistic, linear, loglinear, quadratic, etc.) between the number of features and the runtime? \n",
    "* Does it match your expectations?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21431856-4661-429c-b81e-2927d17a9f11",
   "metadata": {},
   "outputs": [],
   "source": [
    "### TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1840aab3-7a80-4170-baa1-b4b0b0a6cf1e",
   "metadata": {},
   "source": [
    "__Comment:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "alpha-teddy",
   "metadata": {},
   "source": [
    "## 3. Random forests"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "unavailable-integration",
   "metadata": {},
   "source": [
    "### 3.1 Complexity of training a random forest\n",
    "\n",
    "Training a random forest is done by constructing `n_trees` decision trees, each trained on a data set constructing by\n",
    "* subsampling the observations; typically, this is done with bootstrap sampling, where $n$ observations are sampled *with replacement* (i.e. some will be sampled several times and others never, resulting in approximately selecting 2/3 of the samples);\n",
    "* subsampling the features; typically, this is done by randomly sampling $p$ features."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "heavy-holder",
   "metadata": {
    "tags": []
   },
   "source": [
    "##### __Question 10__ \n",
    "* What do you expect the time complexity of training __one__ of the decision trees of a random forest?\n",
    "* How much time do you expect training a random forest of 100 trees on the full data set will take if you're not limiting tree depth? If you're limiting tree depth?\n",
    "* How much space do you expect that random forest model to take in memory if you're not limiting tree depth? If you're limiting tree depth?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf263163-1269-4251-a0e4-b93adffc8d98",
   "metadata": {},
   "source": [
    "__Answer:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "balanced-worth",
   "metadata": {
    "tags": []
   },
   "source": [
    "##### __Question 11__ \n",
    "* Set up an experiment to verify this complexity using the [`sklearn.ensemble.RandomForestClassifier`](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestClassifier.html#sklearn.ensemble.RandomForestClassifier) class. Use at least 100 trees in your random forest. \n",
    "Make sure to use the training times you have observed before to decide on the numbers of samples to use. Do not limit tree depth. Do the numbers match your expectations?\n",
    "* How much space do your trained random forests take in memory? Does this match your expectations?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "illegal-polish",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import ensemble"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61825540-801e-41ef-9a9c-3897de110161",
   "metadata": {},
   "outputs": [],
   "source": [
    "### TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2951862-2f75-44e4-9332-86e49f0042b0",
   "metadata": {},
   "source": [
    "__Comment:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "constant-liberal",
   "metadata": {},
   "source": [
    "### 3.2 Prediction latency\n",
    "\n",
    "Doing predictions in bulk (many samples at the same time) rather than one sample at a time is more efficient. This is due to several reasons, that include CPU caching and [branch prediction](https://en.wikipedia.org/wiki/Branch_predictor). For many models, this can also be due to the efficiency of linear algebra libraries."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "infectious-norway",
   "metadata": {},
   "source": [
    "Here is an example for one random forest, trained on 5,000 samples, with 100 trees:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "minor-footage",
   "metadata": {},
   "outputs": [],
   "source": [
    "ssize = 5000\n",
    "X_current, y_current = subsample_data[ssize]\n",
    "rf_classifier = ensemble.RandomForestClassifier(max_depth=4)\n",
    "\n",
    "# Training\n",
    "rf_classifier.fit(X_current, y_current)\n",
    "\n",
    "# Predicting on 1 sample\n",
    "start_time = time.time()\n",
    "rf_classifier.predict(X_current[:1, :])\n",
    "mytime = time.time() - start_time\n",
    "print(\"Time to predict on one sample: %.3f s\" % (mytime))\n",
    "print(\"Time to predict on one sample (x %d): %.2f s\" % (ssize, (mytime * ssize)))\n",
    "\n",
    "# Predicting on 5000 samples\n",
    "start_time = time.time()\n",
    "rf_classifier.predict(X_current)\n",
    "print(\"Time to predict on %d samples: %.2f s\" % (ssize, (time.time() - start_time)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acquired-pattern",
   "metadata": {
    "tags": []
   },
   "source": [
    "##### __Question 12__ \n",
    "For several values of n, compare the time it takes to predict on n samples to n times the time it takes to predict on 1 sample. Produce a plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f85c8e3e-c260-444a-9502-98672c730da0",
   "metadata": {},
   "outputs": [],
   "source": [
    "### TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2645f46-195c-4d37-914a-20369214f92b",
   "metadata": {},
   "source": [
    "__Comment:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "protected-january",
   "metadata": {},
   "source": [
    "### 3.3 Parallelization"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "unnecessary-immunology",
   "metadata": {},
   "source": [
    "Random forests are embarassingly parallel, meaning that very little effort is required to parallelize them: all of the decision trees that compose it can be trained and applied independently from each other. If your machine has several computing cores, each of them can be used to compute some trees in parallel. `scikit-learn` facilitates this with the `n_jobs` parameter of the `RandomForestClassifier` class."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "quality-semiconductor",
   "metadata": {
    "tags": []
   },
   "source": [
    "##### __Question 13__ \n",
    "* Compare the runtime of training a random forest using a single core (`n_jobs=None`, default), using two cores (`n_jobs=2`), and using as many cores as possible (`n_jobs=-1`). You are free to choose wether you're limiting depth or not, how many samples you use, and so on.\n",
    "* Do the results match your expectations?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06b1a527-67c7-4736-bc88-258865d2cd88",
   "metadata": {},
   "outputs": [],
   "source": [
    "### TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69a94c06-55d6-440d-b088-8061c85a8b35",
   "metadata": {},
   "source": [
    "__Comment:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "entitled-marijuana",
   "metadata": {},
   "source": [
    "## 4. (Optional) Incremental PCA\n",
    "\n",
    "Another technique to reduce computation overheads is `incremental learning`, or `mini-batch learning`, or `online learning`. The idea here is to process only a small number of samples at once (one such set of samples is called a batch). This is typically how neural networks are trained, with gradient updates computed sample per sample or mini-batch per mini-batch.\n",
    " \n",
    "The technique can also be applied to other examples, such as Principal Component Analysis. The principle is described in [Weng et al. (2003), *Candid Covariance-Free IncrementalPrincipal Component Analysis*, TPAMI](http://www.cse.msu.edu/~weng/research/CCIPCApami.pdf). Here we will focus on observing its effects on efficiency.\n",
    "\n",
    "Note that there's a variant of random forests, called Mondrian forests, that allows for online learning, but it is not yet mainstream and it is not implemented in scikit-learn. See [Lakshminarayanan et al. (20140, *Mondrian Forests: Efficient Online Random Forests*, NeurIPS](https://papers.nips.cc/paper/2014/hash/d1dc3a8270a6f9394f88847d7f0050cf-Abstract.html) for classification and [Lakshminarayanan et al. (2016) *Mondrian Forests for Large-Scale Regression when Uncertainty Matters*, AISTATS](http://proceedings.mlr.press/v51/lakshminarayanan16.html) for regression."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "curious-investment",
   "metadata": {},
   "source": [
    "#### Memory profiling\n",
    "\n",
    "Profiling memory usage in a Jupyter notebook can be done thanks to the magic `%memit`.\n",
    "\n",
    "In order to use it, you may need to install [`memory_profiler`](https://pypi.org/project/memory-profiler/) with `conda install memory_profiler`.\n",
    "\n",
    "Here we'll just test this is working:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "solid-words",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext memory_profiler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "incomplete-potter",
   "metadata": {},
   "outputs": [],
   "source": [
    "%memit print(\"Hello World\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "contemporary-murder",
   "metadata": {},
   "source": [
    "Notice that the code was run several times to obtain the estimate.\n",
    "\n",
    "`peak memory` is the total memory usage of your system, including that consumed by running the profiled command. Although this is informative, and can tell you how close you are of running out of memory to run your code, what we're interested in here is `increment`, which tells you by how much global memory usage increased with respect to the memory usage just before running the profiled command."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "breeding-diana",
   "metadata": {},
   "source": [
    "#### PCA vs IPCA\n",
    "\n",
    "The following code can be used to compare classical PCA and incremental PCA with batch size 10 on one of our data sets (with 1500 features). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "blond-location",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import decomposition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "valid-liverpool",
   "metadata": {},
   "outputs": [],
   "source": [
    "fnum = 1000\n",
    "n_subsamples = int(1e4)\n",
    "X_current, y_current = subsample_data[n_subsamples]\n",
    "X_current = augmented_data[fnum]\n",
    "\n",
    "n_components = 2 # number of principal components to compute"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cordless-nicaragua",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Classical PCA\n",
    "pca = decomposition.PCA(n_components=n_components)\n",
    "\n",
    "# Timing the run\n",
    "start_time = time.time()\n",
    "X_pca = pca.fit_transform(X_current)\n",
    "mytime = time.time() - start_time\n",
    "print(\"Time to run classical PCA: %.3f s\" % (mytime))\n",
    "\n",
    "# Evaluating memory usage\n",
    "%memit X_pca = pca.fit_transform(X_current)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "sexual-laptop",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Incremental PCA\n",
    "ipca = decomposition.IncrementalPCA(n_components=n_components, batch_size=10)\n",
    "\n",
    "# Timing the run\n",
    "start_time = time.time()\n",
    "X_ipca = ipca.fit_transform(X_current)\n",
    "mytime = time.time() - start_time\n",
    "print(\"Time to run incremental PCA: %.3f s\" % (mytime))\n",
    "\n",
    "# Evaluating memory usage\n",
    "%memit X_ipca = ipca.fit_transform(X_current)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "biblical-sunday",
   "metadata": {},
   "source": [
    "__Questions:__\n",
    "* Compare the two approaches in terms of runtime and memory consumption.\n",
    "* Do you expect the number of samples or the number of features to have the largest impact on memory consumption for PCA? On runtime?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "funky-marathon",
   "metadata": {},
   "source": [
    "__Instructions:__ \n",
    "* Evaluate the impact of the batch size on both the runtime and the memory consumption of iterative PCA. \n",
    "* Optionally, evaluate how runtime and memory consumption vary as a function of number of features and/or number of samples."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "rapid-committee",
   "metadata": {},
   "source": [
    "The following code allows you to plot the data projected onto the 2 PCs computed by PCA, on the one hand, and iterative PCA, on the other hand. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "subject-turtle",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1, 2, figsize=(11, 4))\n",
    "\n",
    "# Custom color map with 7 discrete values\n",
    "base_cmap = plt.cm.get_cmap('tab10')\n",
    "color_list = base_cmap(np.linspace(0, 1, 7))\n",
    "custom_cmap = matplotlib.colors.LinearSegmentedColormap.from_list('tab10_7', color_list, N=7)\n",
    "plt.cm.register_cmap(name='tab10_7', cmap=custom_cmap)\n",
    "plt.set_cmap('tab10_7')\n",
    "\n",
    "# Projection on the PCs computed by classical PCA\n",
    "axs[0].scatter(X_pca[:, 0], X_pca[:, 1], c=y_current, marker='o')\n",
    "axs[0].set_aspect(1.)\n",
    "axs[0].set_xlabel(\"PC1\")\n",
    "axs[0].set_ylabel(\"PC2\")\n",
    "axs[0].set_title(\"Classical PCA\")\n",
    "\n",
    "# Projection on the PCs computed by incremental PCA\n",
    "ppp = axs[1].scatter(X_ipca[:, 0], X_ipca[:, 1], c=y_current, marker='v')\n",
    "axs[1].set_aspect(1.)\n",
    "axs[1].set_xlabel(\"PC1\")\n",
    "axs[1].set_ylabel(\"PC2\")\n",
    "axs[1].set_title(\"Incremental PCA\")\n",
    "\n",
    "# Adjust space between plots\n",
    "fig.subplots_adjust(wspace=.4)\n",
    "\n",
    "# Colorbar\n",
    "fig.colorbar(ppp, ax=axs, label=\"class\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "graduate-strategy",
   "metadata": {},
   "source": [
    "__Question__ How much did using incremental PCA affect the output of the algorithm?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "exotic-inside",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 5. (Optional) Support Vector Machines\n",
    "\n",
    "__(Optional) Instructions:__ \n",
    "Analyze the runtime and memory usage of training an SVM, as well as the space of the model in memory, for varying numbers of samples and features, in the three following variations:\n",
    "* Linear SVM, optimized via the primal\n",
    "* Linear SVM, optimized via the dual\n",
    "* Kernel SVM.\n",
    "\n",
    "Compare to theoretical complexities."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f37d59fd-502d-4436-9443-29ace1c11915",
   "metadata": {},
   "source": [
    "## 6. (Optional) Decision tree depth and number of nodes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ea2ecf9-7943-4626-8a85-dbee6a3707f9",
   "metadata": {},
   "source": [
    "This section complements Section 1.2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e46157a6-30dd-4eb9-b815-45c2475d734c",
   "metadata": {},
   "source": [
    "Once a decision tree has been trained, its depth and number of leaves can be accessed with its `get_depth()` and `get_n_leaves()` methods. In addition, the full tree structure is stored in its `_tree` argument, which attributes you can find in the help as below, and contain its number of nodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "patient-salem",
   "metadata": {},
   "outputs": [],
   "source": [
    "#help(tree._tree.Tree)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1393e6ef-7e82-4f97-bd53-6dbaaff955a0",
   "metadata": {},
   "source": [
    "__(Optional) Instruction__ Complete the code below to compute, as a function of the number of samples in the data, the number of nodes of the tree as well as its depth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "typical-portrait",
   "metadata": {},
   "outputs": [],
   "source": [
    "tree_attributes = {} # key: number of subsamples, value: {nodes: number_of_nodes, depth: depth}\n",
    "for ssize in subsample_sizes:\n",
    "    X_current, y_current = subsample_data[ssize]\n",
    "    dt_classifier = tree.DecisionTreeClassifier() # default parameters\n",
    "    dt_classifier.fit(X_current, y_current)\n",
    "    tree_attributes[ssize] = # TO COMPLETE"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67aa0555-2704-40f1-b77b-15280ae76803",
   "metadata": {},
   "source": [
    "__(Optional) Questions:__\n",
    "Given the number of nodes in a tree, what depth do you expect the tree to be if it is \n",
    "* fully balanced?\n",
    "* fully unbalanced?\n",
    "\n",
    "Refer to the [document on decision tree complexity](http://cazencott.info/dotclear/public/lectures/Decision_Tree_Complexity.pdf) for help."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f9e7faa-c696-48be-b363-0491e3b4f469",
   "metadata": {},
   "source": [
    "__(Optional) Instructions:__ Plot, as a function of the number of samples: \n",
    "* the observed tree depth \n",
    "* the expected tree depth for a fully balanced tree with the same number of nodes\n",
    "* the expected tree depth for a fully unbalanced tree with the same number of nodes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc39f4ba-5f95-424d-a172-cd9c12d86f8c",
   "metadata": {},
   "source": [
    "__(Optional) Question:__ Do you think the trees you trained are closer to balanced or to fully unbalanced?"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}