Linear Discriminant Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "057ba754-b67d-443a-9fc1-a243d7c8a580",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from sklearn.model_selection import RepeatedStratifiedKFold\n",
    "from sklearn.model_selection import cross_val_score\n",
    "from sklearn.discriminant_analysis import LinearDiscriminantAnalysis "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db022023-85cf-4f37-b448-17ce8afca8dc",
   "metadata": {},
   "source": [
    "## Linear Discriminant Analysis\n",
    "\n",
    "The aim of linear discriminant analysis (LDA) is the discriminate between sample groups, i.e. to classify them, and obtain a measure of the confidence with which a sample can be assigned to a particular group. The discrimination can be achieved by using one or more linear discriminant functions to divide the sample groups. A ‘linear discriminant function’ can be visualised as a surface dividing a data space into different regions, the simplest case being a ‘binary classifier’, where the data space will be divided into two regions (Haswell, 1992).\n",
    "\n",
    "A LDA has been performed on data from two sets of home brew:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "2eb7cc8d-304a-4687-bdaa-60ad20ea2cc5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      glycerol%  acetaldehyde%\n",
      "Brew                          \n",
      "1          1.38           1.64\n",
      "1          1.40           1.70\n",
      "1          1.24           1.72\n",
      "1          1.36           1.74\n",
      "1          1.38           1.82\n",
      "1          1.48           1.82\n",
      "1          1.54           1.82\n",
      "1          1.38           1.90\n",
      "1          1.56           2.08\n",
      "2          1.14           1.78\n",
      "2          1.20           1.86\n",
      "2          1.18           1.96\n",
      "2          1.30           1.96\n",
      "2          1.26           2.00\n",
      "2          1.28           2.00\n"
     ]
    }
   ],
   "source": [
    "# Read data\n",
    "\n",
    "brew_data = pd.read_csv(\"Data/LDA_data.csv\", index_col=\"Brew\")\n",
    "print(brew_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "5d574108-e7c5-4ddc-84a6-3fe45e60a94a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuIAAAIjCAYAAABGV34uAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAB5kElEQVR4nO3dd1hTZxsG8DsJW/ZwAOJAxQHuvXBrFfeqdbRVq1Wr1lVntbUOXK3V2qqlirvWvUfFVbVqxdmCuKoi4iKAbMj4/shHakjAgCQnkPt3XV6t70lOHuQY7zy873tESqVSCSIiIiIiMiqx0AUQEREREZkjBnEiIiIiIgEwiBMRERERCYBBnIiIiIhIAAziREREREQCMPsgfuvWLdy6dUvoMoiIiIjIzFgIXYDQ5HK50CXQW0RFRQEA/Pz8BK6EhMJrgABeB6TC64CKE7PviBMRERERCYFBnIiIiIhIAAziREREREQCYBAnIiIiIhIAgzgRERERkQDMfteUt1EqlXj16hXS09OhUCiELscspaenAwAePXokcCUkFFO5BsRiMWxsbODu7g6RSCRoLUREVPQxiOdBqVQiJiYGSUlJsLa2hkQiEboks+Tg4CB0CSQwU7kGsrKykJycjIyMDHh5eTGMExHRO2EQz8OrV6+QlJSEUqVKwdXVVehyzFZ2N9TGxkbgSkgopnQNSKVSPH/+HK9evYKHh4fQ5RARURHGOeJ5SE9Ph7W1NUM4Eam5urrC2tpa/eGAiIiooBjE86BQKDgdhYi0SCQSrhkhIqJ3xiBORERERCQABnEiIiIiIgEwiBMRERERCYBBnPIlPDwc/v7+ej9eqVRi+PDh+PHHHzXGMzMzsWjRIrRo0QK1atXCwIEDcf369UKtNS0tDV9++SUaNWqE+vXrY9asWUhJSdHruVlZWejRowemTZum83hmZia6deuGffv26XW+mTNnYseOHQCAwYMHY+bMmbk+dvfu3ahevbpe5zUFly5dgp+fH549e2bw15o2bRo++ugj9e/37t2L5s2bo2bNmvj999/h5+en9/dEl7i4OHTu3BlJSUmFUC0REVHeGMRJbzdu3MDo0aP1XqSWmZmJmTNn4o8//tA6tmjRIhw9ehSLFy/GgQMH4Ofnh48//hjPnz8vtHpnz56N8PBwrFmzBqtXr8bly5cxe/ZsvZ67YsUKREZG6jyWnJyMMWPGICoqSq9z/fnnn7h58yZ69+6t1+M7d+6Ms2fP6vVYczNz5kx8//336t9nf5g7cuQImjdvjnPnzqFTp04FPr+bmxu6du2KRYsWFUa5REREeWIQJ70sXboUAwcOhKenp16P/+eff9CvXz9cunQJjo6OOh/z5ZdfokmTJvDx8cGECROQmpqKmzdv5quu7G5sTs+ePcPBgwcxZ84c1K5dG/Xr18e8efNw6NCht4b98PBw7Nq1C1WqVNE6duHCBfTo0QNxcXF617h8+XIMGTIEYrF+f92y79xI2hwcHODk5KT+/evXr1G/fn14eXnB1tYWHh4esLa2fqfXGDRoEA4cOCD4XTyJiKj4YxAviKws4OFDYX5lZeWr1Pnz56NLly4aY48fP4afnx8iIyMxbdo0+Pn56fy1cuVK9XPOnj2LNWvWYPDgwXq97vnz51G/fn3s27dP510Rv/zyS7Rp0waAqsMcEhICBwcH1KxZUz0tpH///urue3h4OKpVq4ajR4/q9fpXr16FWCxG3bp11WN169aFRCJBeHh4rs9LSUnB1KlTMWvWLLi5uWkdP3nyJHr06IFff/1VrzquXbuGf/75B+3bt9cYT05Oxrhx41CzZk20atUKoaGh6mM5p6b4+flh586dGDhwIAICAtCqVSts375dfTwjIwMLFy5E69at4e/vj8aNG2P69OlIS0tTn69jx4746quvUK9ePUycOBGNGzdGSEiIRk3ff/89unfvrvPrUCqVCA0NRYcOHVCrVi10794dZ86c0fnYJ0+eYNy4cWjUqBFq1KiBNm3aaLzWgwcPMHToUNStWxf16tXD6NGj8eTJE/XxtWvXom3btvD390fHjh2xZcsW9bHsqSlPnjyBn58fZDIZZsyYob6Wck5N+e2339CxY0fUrFkTXbt2xZ49e9THLl26hICAAPz4449o2LCh+tp2cHBAkyZNsGHDBp1fHxERUWHhnTXzKysLqFYNuH9fmNf39QUiIwFLS70e3rNnT2zcuBERERHqcLd//35UrVoV1apVw8yZMzFp0iSdz7Wzs1P///79+wGoQp0+RowYodfjQkNDsXDhQohEIixcuBClSpUCoJpy0Lt3b2zevBmdO3fGrFmz0LNnT72nHTx//hyurq6wfOPPycLCAq6uroiNjc31eQsWLEBAQAA6d+6M3377Tev4rFmz9Hr9bCdPnkStWrXg7OysMX706FEMGzYMe/fuxeXLlzFv3jyULl06169v6dKlmD17NoKDg7Fx40Z89dVXaN68Oby8vLBo0SKcO3cOS5YsQenSpXHz5k31B6zs+dQPHz5EQEAA9u7di8zMTLi6umL//v0YPnw4AFXQ3r9/f64ftH7++WesXr0as2fPRt26dXHo0CGMGTNGI9hmGzVqFLy8vLBx40bY2Nhg7969WLJkCZo1a4Zq1aph8uTJqF69OubMmYO0tDR8/fXXmDFjBjZu3IiTJ0/il19+wfLly+Hj44MLFy7gyy+/RPny5VGvXj31a5QpUwbnzp1DYGAgpk6diq5du2rVsXXrVqxcuRJz5sxB9erVce3aNXzzzTcAVH8vANX0qUuXLmHHjh0aN+gJDAzEmjVr9J7KREREVBAM4sVc9erV4efnh/3792sE8Q8++ACAqvunq2NtLG3btkWjRo1w9OhRzJw5E66urggMDISfnx/GjRuH5cuX488//4SFhYU6BD99+lTd5ZfL5QCAOnXqqM957do1pKWl6ZyiYGVlhYyMDJ21hIWF4cyZMzh48GChfX03btxApUqVtMZr1aqFL774AgBQsWJF3LhxAxs2bMg1iPfu3RudO3cGAIwbNw4bN27EzZs34eXlhVq1aqFLly7qoOrt7Y2tW7fizp07GucYPXo0ypYtCwDo1asXNm3ahDt37qBKlSoIDw/H8+fP0a1bN63XViqV2LhxIz7++GP06NEDgCpsy2QypKamajw2PT0dPXv2RJcuXdQfqsaMGYPVq1cjKioK1apVw6NHj9CsWTN4eXnBwsICS5YswatXrwCoflpjaWkJT09PeHl5oW/fvvD29oaPj4/G60gkEvXt5R0cHHTe/Xb16tX47LPP1H+mPj4+ePr0KVavXq0O4gAwfPhwlCtXTuO5lStXRmxsLJ4/f67+OoiIiAobg3h+WVqqOtIxMcK8vpeX3t3wbL169UJISAi++OIL3LhxAzExMeoO4uzZs3HgwAGdzxs5ciQ+/fTTdy45L9nBsFq1avjnn3+wYcMGBAYGAlAFpN9//x0nT57Exo0b1R36kiVLYu/evQBUQXfKlCnq32ezsbFBZmam1utlZmZqdPqzSaVSfPnll1iwYIFW9/pdxMXFaUyPyfbmBwcA8Pf3x++//57recqXL6/+/+wPTln/n6bUvXt3nDt3DosXL8bDhw9x7949PH78GN7e3urniEQijd9nf0Dbt28fpkyZgv3796NFixY6A218fDxevnyJmjVraoyPHTsWgGqKRzYbGxsMGjQIhw8fxs2bN/Ho0SNERkZCoVCopxmNHz8eixYtwtatW9G4cWO0atVK/cGqW7du2LlzJzp06IAqVaqgefPmCAoK0jlNKC9SqRTPnz/HokWLsHTpUvW4TCaDXC7XuDayr8E3Zf85xMXFMYgTUdF39ixw5QowaBBQsqTQ1dAbGMQLwtISeCMYmbquXbtiyZIluHTpEo4fP46WLVuqg8348eMxbNgwnc97c1FcYcrMzMSZM2dQu3ZtdVcTAKpUqaIx7zghIQExMTGQSCT4888/0ahRIwCqKSbZHczsLfNydjRLly4NqVQKuVwOiUQCQBXCpFIpSup4Ezpz5gzi4uIwYcIE9VhGRgZEIhGOHTuGa9euFehrFYlEOneZyblwU6lUwsrKKtfz6DqmVCoBqHYSCQsLQ8+ePdGhQwdMmDABc+fO1Xq9nOfo1asXQkNDMX78eBw5cgQLFizQ+dqW+fjgl5KSgoEDB0Iul6Njx45o1KgRatWqhdatW6sfM2TIEHTu3BmnTp3ChQsXsHDhQqxbtw779u1TT5kJDw/HuXPncObMGaxbtw5z587Ndf56XjV/+eWXaNiwodZxC4v/3vpsbGy0jmd/z/RdYEtEZJKUSqB9eyAsTPX7V6+AXN7rSRgM4mbAzc0NLVu2xPHjxxEWFqYxz9nNzS3f3cZ3JZFIMHXqVIwZM0bjQ8CtW7fg6+ur/v3s2bNRqlQpTJ48GV9++SXatWun9x7m9erVg0wmw7Vr11C/fn0AqgWfCoVCY65xtvbt22t1rqdOnQoPDw9Mnjy5IF8mAMDDwwNSqVRrPOfWiFevXtU5heVt4uPjsXPnTqxcuRIdOnQAoPrAER0d/dYdbrI/oK1btw5isVj9k4icHBwc4OHhgVu3bmk8ZvDgwQgMDERAQIB67Ny5c4iMjMSlS5fUP1l48OABFAoFlEol4uPj8cMPP+CTTz5B37590bdvX9y8eRN9+/bF7du38eTJE8THx2PgwIFo0KABJkyYgE8++QTHjx/PVxB3cHBAqVKl8OTJE/Tt21c9vm3bNkRGRmp9UMkp+3um60MbEVGR8PgxkKNJhTfeD8k0MIibiV69emHy5MmwsbFBq1atCv38UqkUlpaWes03l0gkGDx4MFavXo1y5cqhYsWK2LFjB65fv65eILl3716cOnUKO3bsQMWKFXH06FFMnToVe/bs0ejsNmrUSOd+3qVKlcJ7772HmTNnYsGCBVAqlfjyyy/RvXt39VSD9PR0JCUlwdXVFfb29rC3t9c4h42NDUqUKKHVbc+PmjVr6txH/dKlS1ixYgWCgoJw9uxZHDlyBOvWrcv3+bPrDgsLQ9WqVZGcnIw1a9YgNjZW59ScN2V/QPvpp5/Qt2/fPDvyw4cPxw8//IAKFSogICAABw8exI0bNzB79myNDxplypQBABw4cABt2rTB48ePsXDhQgCqn4Q4OTnh7NmziI6OxsSJE2Fra4vdu3fD0dERFSpUwIMHD7Bo0SI4ODigXr16ePz4MSIiIjTCtL5GjRqF4OBgeHp6okmTJrhx4waCg4PVC1TzEhERAW9vb51TdYiITN7SpcCUKZpjt28DOrb7JWExiJuJVq1awcbGBkFBQXkGroLq06cPGjZsiODgYL0eP3bsWFhbW2PBggV4+fIlatSogdDQUFSrVg3Pnj3D/PnzMXToUFSvXh3p6emYOXMmevbsie+++w5Tp07V6zXmzZuHefPmYcSIEbCwsEDHjh0xY8YM9fHDhw9j+vTpCAsL05g/XZiyt+5LTEzUmOrTv39/RERE4Oeff0bp0qURHByMJk2a5Pv8lpaWWL58ORYtWoSgoCC4urqiZcuWGDp0KE6cOPHW5/fo0UO9JWNehgwZgvT0dCxZsgRSqRSVK1fG6tWrUblyZY054jVr1sQXX3yBn3/+GUuXLoWnpyd69+6Ns2fP4tatWxgwYADWrFmD4OBgDB48GJmZmQgICMAvv/wCBwcH9R7tK1euRGxsLNzc3NCrVy+9wnNOAwYMQGZmJn755Rd88803KFWqFEaPHq3Xjj6XLl1C27Zt8/2aRESCysgASpQA/r+RAQCgZk3g+nVAJBKsLMqdSJk90dRMZd9WvXbt2lrHsm/o8S4dUVMRHx+PFi1aYPv27ahRo4bQ5eRL9rZyuubyFgV9+/ZF9+7dMWjQIKFL0bJ582b89ttv6u0pTZUxrwGpVIrWrVtj//79uf7dL07vDUVJ9k+/dN3Ei8wHr4NcnDsHtGihOfbbb5yOYuK4EqmYi4+Px7FjxzBjxgz4+/sXuRBeHIwfPx5btmxRb7VoCv7++2/s27cPq1evNskPCELaunUrOnfuzJBNREVH377aITw+niG8CGAQL+ZkMhlmzpyJR48eYf78+UKXY5aaN2+O2rVrY+fOnUKXonb16lXMnj0bzZo1Q58+fYQux2TExcXh4MGDmDZtmtClEBG93fPnqiknb/77MmaMareUQtyKlwyHU1PMZGpKUVbUp6bQuzO1a4DvDcLglAQCeB2orV0LjBypOXb9OlCrliDlUMFwsSYRERFRUSGTAZ6ewMuX/435+AAPHgD/v28GFR2cmkJERERUFISHq24q+GYIX7cOePSIIbyIYkeciIiIyNR98gkQEqI59uIF8MYdqqnoYRAnIiIiMlVSKZDzDtgDBwKbNwtTDxUqTk0hIiIiMkXbtmmH8AsXGMKLEXbEiYiIiEyJQgFUrw78f4cYAIC9vao7bmkpXF1U6NgRN0OzZ8/GzJkzhS6DiIiIcoqIUC28fDOEf/89kJTEEF4MMYibEaVSie+//x7bt28XuhQiIiLKacoUIOcdsJ88AcaNE6YeMjhOTTESpVKJrKREZMZLoZBnQSyxhJWLKywdnCASiQz++tHR0ZgxYwbu3r0LT09Pg78eERER6SkpCXB01Bzr0gU4eFCYesho2BE3AoUsC6/vRSL54T1kJkohS05CZqIUyQ/v4fW9SChkWQav4erVqyhTpgwOHDgAb29vg78eERER6WH/fu0QHhbGEG4m2BE3MKVSiaR/70KelqrzuDwtFUn/3oVjpWoG7Yx3794d3bt3N9j5iYiIKB+USkCsox+algbY2Bi/HhIEO+IGlpWUmGsIzyZPS0VWUqKRKiIiIiJB/f67dgj/5htVOGcINyvsiBtYZrxUv8clSGHl6GzYYoiIiEhYVlZAVo4pqX/8ATRvLkw9JCgGcQNTyPWb/22MeeJEREQkkPh4wNVVe1yhAIywaQOZJk5NMTCxRL89P8UW3BuUiIioWJoxQzuEf/KJaioKQ7hZY0fcwKxcXJGZ+PbpKVbOOj4lExERUdGmK2g/fQqUKWP8WsjksCNuYJYOTpDY2uX5GImtHSwdnIxUERERERncpUu6Q7hSyRBOauyIG5hIJIJDhcq5bmEosbWDQ4XKRrmpT7ZNmzYZ7bWIiIjMTsWKwL//ao7t3An07i1MPWSyGMSNQGxhCcdK1VR31kyQQiHLgtjCElbOxruzJhERERlYSgpgb689LpMBEonx6yGTx6kpRiISiWDl6Ax7n4pwrOgHe5+KsHJ0ZggnIiIqDpYt0w7hXbuqpqIwhFMu2BEnIiIiehe6mmr376umqBDlgR1xIiIiooKIiMh9QSZDOOmBQZyIiIgov1q0AGrU0Bxbu1YVwon0xKkpRERERPrKzASsrbXH09N1jxPlgR1xIiIiIn2EhmqH7QYNVF1whnAqAHbEiYiIiN5G11zwmzeBgADj10LFBoM4ERERUW4ePgQqVNAe51xwKgScmkJERESkS58+2iE8OJghnAoNO+Jm4NWrV1iyZAnOnz+P9PR01KpVC1OnTkWVKlWELo2IiMj0KBS6b8KTlKT7zplEBcQgbiRyuRznT1/Gwb3HIX2VAFd3ZwT16IBmrRpCYsA7bikUCnz22WdQKpX48ccfYWdnh5UrV+Kjjz7CoUOH4OLiYrDXJiIiKnL27QN69NAc8/QEYmIEKYeKNwZxI4h7FY8xH32BiFt3IBaLoVAoIBaLcXT/SVQPqIJVoYvh5m6YQHz79m1cu3YNhw8fhq+vLwBgyZIlaNiwIc6cOYMeOd9siIiIzJWuBZnnzwNNmxq/FjILnCNuYHK5HGM++gKR/9wFoOpQv/nfyH/uYsxHX0Aulxvk9cuUKYM1a9agwhtz3ET/f6NJTEw0yGsSEREVKS9e6A7hCgVDOBkUg7iBnT99GRG37kCp0L2wQ6lQIuLWHVw485dBXt/FxQWtWrWCWPzft3rTpk1IT09H8+bNDfKaRERERcaYMUCpUppjkyerFmTqCudEhYhTUwzs4N7j6ukouRGLxTi45zhatGls8HrCwsLw7bff4uOPP1ZPVSEiIjI7SiUg1tGPfPUKcHMzfj1kltgRNzDpq4Q8QzigmqYijUsweC27d+/GuHHj8N5772HKlCkGfz0iIiKTdPq07hCuVDKEk1ExiBuYq7uzxrQQXcRiMVzdnA1ax08//YTp06fj/fffx+LFi99aExERUbHk7Ay0bq05duQI9wYnQTCNGVhQjw56dcSDenYwWA0///wzli9fjnHjxuHLL79UL9YkIiIyG4mJqjnfOTcqkMuBTp2EqYnMHoO4gTVr1RDVA6pAJNYdfkUiEaoHVEHTwAYGef3bt2/ju+++Q+/evdGvXz+8fPlS/Ss1NdUgr0lERGRSvv5a1Ql/05Ahuc8TJzISLtY0MIlEglWhi3XuI65QKFDNvzJWhS422E19Dh8+DLlcjl27dmHXrl0ax8aPH4/Ro0cb5HWJiIhMgq6fAj95Anh5Gb8WohwYxI3Azd0FW/atxoUzf+HgnuOQxiXA1c0ZQT07oGlgA4PeWXPixImYOHGiwc5PRERkkq5eBerV0x7nXHAyIQziRiKRSNCiTWOjbFFIRERk1mrUACIiNMe2bQPef1+YeohywSBORERExUNaGmBnpz2elQVYMPKQ6eEKBSIiIir6fvhBO4R36KCaisIQTiaKVyYREREVbboWZN65A1SubPxaiPKBHfE8iMViyOVyocsgIhMjl8t5UywiU3Dnju4QrlQyhFORwH9J8mBjY4OMjAxIpVKhSyEiEyGVSpGRkQEbGxuhSyEyb+3bA35+mmOrVnFXFCpSODUlD+7u7sjIyMDz58+RkJBg0G0GKXfZP5Xgn7/5MpVrQC6XIyMjAw4ODnB3dxe0FiKzJZPBz99fezwtDeAHZCpiBO2IKxQKbNu2DV27dkWdOnXQrl07LFy4EMnJyXo9f8uWLWjfvr3B6hOJRPDy8oK7uzssLS0N9jqUt6SkJCQlJQldBgnIVK4BS0tLuLu7w8vLCyJdPw4nIsPaulU7hAcEqLrgDOFUBAnaEQ8JCcHy5csxbNgwNGnSBP/++y9WrFiBe/fu4ZdffsnzucePH8fChQtRpkwZg9YoEong4eFh0NegvKWnpwMAypUrJ3AlJBReA0Skcy741atAnTrGr4WokAgWxJVKJUJCQtC/f39MmjQJANC0aVO4uLhgwoQJiIyMRLVq1bSel5iYiJUrV2Lz5s1wdHQ0dtlERERkTE+eAGXLao9zLjgVA4JNTUlJSUG3bt0QFBSkMV6xYkUAwOPHj3U+b+PGjfj999/x3XffoU2bNgavk4iIiAQyaJBWCH81bhyibt8WqCCiwiVYR9ze3h6zZs3SGj9x4gQAoFKlSjqfFxQUhJEjR8LKygpnzpwxaI1EREQkAIUC0LU4OzERcbGxxq+HyEBMateUGzduYO3atWjXrh18fX11PqZChQqF/rpyuRxRUVGFfl4qHKmpqQDA75EZ4zVAAK8Dc1HizBl4jxypMSZ3csK9S5eA2FheB0WEX86tJUknkwni4eHh+PTTT+Ht7Y158+YJXQ4REREZmV/VqlpjjzduRFrDhgJUQ2R4JhHEDx8+jGnTpqF8+fIICQmBi4uLUV9fIpHwk5sJy+568HtkvngNEMDroFh79QrQtUOZQgGfHLul8Dqg4kTwO2uuX78eEydORO3atbFlyxaULFlS6JKIiIjIWCZN0g7hY8eqdkXhfv1UzAnaEd+xYweCg4PRuXNnLFq0CFZWVkKWQ0RERMaiVAJiHf3AFy90d8eJiiHBgnhcXBzmz58PLy8vDBw4EBERERrHfXx8YGVlhXv37sHHxweurq4CVUpERESF6vx5oHlz7XHuDU5mRrAg/scffyAtLQ0xMTEYOHCg1vHFixejdOnSGDJkCBYuXIhevXoJUCUREREVKk9PIOcWhPv2Ad26CVMPkYAEC+I9evRAjx493vq4vLYnCg4OLsSKiIiIyGCSkwEHB+1xuVz3FBUiM8Arn4iIiAwrOFg7hPfrl/s8cSIzYRLbFxIREVExpWvnk4cPgXLljF4Kkanhx1AiIiIqfLdu6Q7hSiVDONH/MYgTERFR4WrQAKhZU3NswwbuikKUA6emEBERUeHIyABsbLTHMzMBS0vj10MAAIUsS32DJLEFvw+mhB1xIiIienc//6wdwlu2VAVAhnBBKGRZyEiQIunBHSRE/Y2kB3eQkShVBXMyCeyIExER0bvRNRc8IgKoVs34tRAAVQhPfRaDTOkr9Zg8PQ0pjx7A2tUDtqW9ILZgDBQaO+JERERUMPfv574gkyFcUIqsLI0Q/qYM6Uso2RU3CQziRERElH/dugGVKmmOffstF2SaiMx4aZ7HM+LjjFQJ5YU/kyAiIiL9yeWArikNKSmAnZ3x6yGdFApZnseV8ryPk3GwI05ERET62blTO4RXqqTqgjOEmxQrR5c8j1s65X2cjIMdcSIiIno7XXPBL19W7RlOJsfC1hZiK2soMjO0jomtrGFhYytAVZQTO+JERESUu9jY3BdkMoSbLLGlFRwqVoFFCXuNcYsS9nCoWAViSyuBKqM3MYgTERGRbiNGAJ6emmMzZ3JBZhEhsbKGfTlfOPn5w7FSNTj5+cO+XCVIrKyFLo3+j1NTiIiISJNSCYh19Ori4wFnZ6OXQwUntrAEeDdNk8WOOBEREf3nxAntEG5jowrnDOFEhYodcSIiIlKxsQEycizu+/13oF07YeohKuYYxImIiMxdfDzg6qo9rlDoXqhJRIWCU1OIiIjM2cyZ2iH8k09UU1EYwokMih1xIiIic6UraMfGAqVLG78WIjPEjjgREZG5uXw5973BGcKJjIZBnIiIyJxUrgw0aqQ5tnMn9wYnEgCnphAREZmD1FSgRAntcZkMkEiMXw8RsSNORERU7H33nXYI79ZN1QVnCCcSDDviRERExZmuueD37wMVKxq/FiLSwI44ERFRcRQRkfuCTIZwIpPAIE5ERFTctGwJ1KihOfbzz1yQSWRiODWFiIiouMjKAqystMfT0wFra+PXQ0R5YkeciIioONi4UTuEN2qk6oIzhBOZJHbEiYiIijpdc8Fv3gQCAoxfCxHpjUGciIioqHr0CChfXnucc8GJigROTSEiIiqK+vfXDuHBwQzhREUIO+JERERFiUKh+yY8SUmAvb3x6yGiAmNHnIiIqKjYv187hHt7q7rgDOFERQ474kREREWBrgWZ588DTZsavxYiKhQM4kRERKbs5UugZEntcYVCdzgnoiKDU1OIiIhM1bhx2iF8yhTVVBSGcKIijx1xIiIiU6NUAmIdvbJXrwA3N+PXQ0QGwY44ERGRKTlzRncIVyoZwomKGQZxIiIiU+HqCrRqpTl25Aj3Bicqpjg1hYiISGivXwNOTtrjXJBJVKyxI05ERCSkuXO1Q/iQIVyQSWQG2BEnIiISiq6g/eQJ4OVl/FqIyOjYESciIjK2a9d0h3ClkiGcyIwwiBMRERlTzZpA3bqaY9u2cUEmkRni1BQiIiJjSE8HbG21x7OyAAv+c0xkjtgRJyIiMrQff9QO4R06qLrgDOFEZot/+4mIiAxJ11zwO3eAypWNXwsRmRR2xImIiAzhzp3cF2QyhBMRGMSJiIgKX4cOgJ+f5tiqVVyQSUQaODWFiIiosMhkgKWl9nhaGmBjY/x6iMiksSNORERUGH79VTuE16yp6oIzhBORDuyIExERvStdc8GvXgXq1DF+LURUZDCIExERFVRMDODtrT3OueBEpAdOTSEiIiqIIUO0Q/jcuQzhRKQ3dsSJiIjyQ6EAJBLt8cREwNHR+PUQUZHFjjgREZG+jhzRDuGurqouOEM4EeUTO+JERET60LUg88wZoGVL49dCRMUCgzgREVFe4uIAd3ftcYVCdzgnItITp6YQERHlZsoU7RA+dqxqKgpDOBG9I3bEiYiIclIqAbGOXtWLF4CHh/HrIaJiiR1xIiKiN124oDuEK5UM4URUqBjEiYiIsnl7A82aaY7t38+9wYnIIDg1hYiIKDkZcHDQHpfLdXfHiYgKAd9diIjIvC1apB3C+/XLfZ44EVEhYUeciIjMl66dTx49Anx8jF8LEZkdftQnIiLzc+uW7hCuVDKEE5HRMIgTEZF5adQIqFlTc2zDBi7IJCKj49QUIiIyDxkZgI2N9nhmJmBpafx6iMjssSNORETFX0iIdghv2VLVBWcIJyKBCNoRVygU2L59O7Zu3YonT57Azc0Nbdu2xdixY2Fvb6/zOTKZDD/88AP27NmDhIQE1KhRA9OmTUPNnD9mJDJj0rh4yGRyWFlZwtnFyeCvY2NjDUcnHVu/EZkCXXPBIyKAatWMXwsR0RsEDeIhISFYvnw5hg0bhiZNmuDff//FihUrcO/ePfzyyy86nzN//nzs2bMHkydPhqenJ9avX4+PPvoI+/btQ9myZY38FRCZFmlcAv768xp++XELnse+QLmKZTHq849RPaAKnJwdC/F14nHu9GVsWPMrXr2MQ+WqFTF64lBU9qsIB0fdH6KJjO7BA8DXV3ucc8GJyEQINjVFqVQiJCQE/fv3x6RJk9C0aVMMHDgQc+bMwblz5xAZGan1nCdPnmD79u2YOnUqBg0ahDZt2uCXX36Bk5MTQkJCBPgqiExHvDQBS+b+gCljvsLtf+4iXpqI61f+xshBk7A1dDeSXicXyutIX8Vj9pTFmDVxAe5GPUC8NBGXL1zDR33G4uDu40hNSS2U1yF6J927a4fw775jCCcikyJYEE9JSUG3bt0QFBSkMV6xYkUAwOPHj7Wec/HiRcjlcnTs2FE9ZmVlhVatWuHs2bOGLZjIxMVEP8Ohvb/rPLZ6eSgSE14XyuvcjfoXZ8Mu6Dy25JsfkJiQVCivQ1QgcrlqKsr+/ZrjKSnA558LUhIRUW4Em5pib2+PWbNmaY2fOHECAFCpUiWtYw8ePICTkxNcXV01xsuVK4enT58iPT0dNrpWxL+FXC5HVFRUvp9HxpGaquqw8nuUO3t7e/y2aW+ux5VKJU4eP4dWHRojLS2twK9ja2uHraG7cj0uk8lx5dJ1BNT1Q0ZGRoFfJydeAwS8/TqwP3oUXjnCdma5cvj32DEgOtrQ5ZGR8P2gaPDz8xO6hCLBpLYvvHHjBtauXYt27drBV8e8vuTkZJ2LOEuUKAFA1WUvSBAnKuqUCiWSkvKeepL0OhkiXYvW8vk6KckpeT4muRBehyi//KpW1Rp7tGMH0gMCBKiGiEg/JhPEw8PD8emnn8Lb2xvz5s3T+RjlW+b2FfQff4lEwk9uJiy768HvUd5at2+OsKN/5Hq8acsG77ygWS6Xo0WbJrh84Vquj6nbsBbKlSv3Tq+TE68BAnK5Dp49A8qU0X6wUonCvQrJVPD9gIoTk9hH/PDhw/j4449RpkwZhIaGwsXFRefj7O3tkZKi3Y1LTk5WHycyV42b14Obh6vOY76Vy6NsOa93fg2JRIKOQa1z3RmlZt0a8Cjl9s6vQ6SXkSO1Q/jMmVyQSURFhuBBfP369Zg4cSJq166NLVu2oGTJkrk+tmLFikhISEBiYqLG+KNHj+Dt7Q0rKytDl0tksjxKuSN0x0r4VddcX9GwaR38tGkJ3HMJ6fl+nZJu2LDrB1SopNlvbNG6Mb5dPReubro/SBMVGqVStSBz7VrN8fh4IJefqBIRmSJBp6bs2LEDwcHB6Ny5MxYtWvTWIN20aVMAwLFjx9CvXz8AQGZmJk6fPo0WLVoYvF4iUyYSiVCugjfWbFqK16+T8TrhNVzcnOHg6ABnl8LbQ1wikaBSlQpY9+tyvH6djKTEJLh5uMLB0Z439SGDs/vzTyDnfHBbWyCV22YSUdEjWBCPi4vD/Pnz4eXlhYEDByIiIkLjuI+PD6ysrHDv3j34+PjA1dUVXl5e6NmzJ+bNm4fU1FSUK1cO69evx+vXrzF8+HCBvhIi0+Lq7gJXd8N3pd08XHOdCkNkCLoWZOLECaBtW+MXQ0RUCAQL4n/88QfS0tIQExODgQMHah1fvHgxSpcujSFDhmDhwoXo1asXAGDu3LlwdHTE2rVrkZqaiho1amD9+vWFvjiMiIhMRG4LMhUK3bevJyIqIkTKt21FUsxdv34dAFC7dm1B66DccYU88RowY82aARdy3EBq+HDg55+FqYcEx/cDKk7y1RHPzMzE9evX8eTJE8THx0MikcDNzQ1lypRB7dq1YWFhMrshEhFRUaej2/3g2DFU7NBBgGKIiAqfXsn51KlT2Lx5My5fvgyZTKa1n7dIJIKtrS0aN26Mfv36oVWrVoaolYiIzMG+fUCPHlrDUbdvG78WIiIDyjOInz17FvPnz0d0dDRq166NYcOGoUqVKihbtizs7e2hUCiQkJCA58+f4/r16+qb8lSsWBETJ05Eu3btjPV1EBFRcaBrzvc33wCzZgG8pTkRFTO5BvGxY8fiypUrGDJkCHr37p3n/t4A0LlzZwCqPb337NmDmTNnYt++fVi5cmXhVkxERMXP69eAk5P2eFYWwGmPRFRM5XpDn+rVqyMsLAyjRo16awh/U7ly5fD5558jLCyMCymIiOjtBgzQHcKVSoZwIirWcn2HGzVq1Dud2N7eHp999tk7nYOIiIo5XVNRwsOBunWNXwsRkZEJfot7IiIyQ3/8oTuEK5UM4URkNgoUxDMzM7F8+XJ06NABtWvXRseOHbFixQpkZmYWdn1ERFTciERAy5aaY6NHq0I4EZEZKdDku0WLFuHMmTPo168fnJ2dER0djU2bNiE+Ph5z5swp7BqJiKg4yMgAbGy0x1NSADs749dDRCSwXIN4cnIy7O3tdR47fvw4li9fjnr16qnH3Nzc8NNPPzGIExGRtsmTgWXLtMfZBSciM5br1JQ2bdpg9erVSElJ0Trm6uqKc+fOQaFQAADS09Px119/wc3NzXCVEhFR0SQSaYfw339nCCcis5drEF+yZAnCwsLQpk0brF27FqmpqepjX3zxBUJDQ1GvXj0EBgaiYcOGuHDhAqZPn26UoomIqAi4dSv3BZm84RsRUe5TUwIDAxEYGIiTJ0/ihx9+wLp16zBs2DAMGjQIzZo1w8mTJ3Hq1CnExcXBw8MDLVq0YEeciIhULC0BmUxzrGtXYP9+YeohIjJBb12s2aZNG7Rp0wYnTpxQB/Lhw4dj4MCB6NWrlzFqJCKiokIu130Tnrg4wNXV+PUQEZkwvbcvbNeuHfbu3YuvvvoKe/fuRdu2bREaGoqMjAxD1kdEREXFsmW6Q7hSyRBORKRDnh3x6OhoHDlyBM+ePYOLiwuaNGmCjh07omPHjjh8+DBWrVqFkJAQfPLJJxgwYACsrKyMVTcREZkSXXPBt28H+vUzfi1EREVErh3xEydO4L333sPmzZsRERGBw4cPY/DgwVi0aBEAoHPnzjh48CCmTp2KX3/9FW3btsWmTZuMVjgREZmAhATdIVyhYAgnInqLPHdN6datG86cOYNff/0VR44cwbJlyxAaGorY2FgAgEgkQteuXXHo0CFMmjQJmzdvNlrhREQksJUrARcXzbFatVRTUXSFcyIi0pDr1JQXL16gfPnyEL3xZurj4wOlUonXr1+jTJky6nGxWIwePXqgW7duhq2WiIiEp1QCYh19nOhowNvb+PUQERVRuQbxrl274ocffkBUVBS8vb2RnJyMY8eOISAgAFWqVNH5HLGuN2YiIio+btwAatfWHAsMBE6fFqIaIqIiLdcg/tVXX6FGjRo4cuQIIiIi4OLigt69e2Po0KEaXXIyT0lJyUhPTYdEIoGru8vbn0D0DqytbAEAyUkpsHcoIXA1Zqx7d+19wC9eBBo1EqYeIqIiLtcgLhaL0b9/f/Tv39+Y9ZCJS0lOwaN/n+Cn70Lx983bcHV1xqBhfdGybRO4MZBTIYt7JcXZsD+x+ZedkEoT4F+zKkZN+AjlKnijhD0DudGkpAD29trjCgXnghMRvYNc55Jcu3btnU9+5cqVdz4HmY6M9AycOfEn3g8agTNhFxD3Uoq7UQ8w54tFmPrZ14h7JRW6RCpG4l7F44sxczHni8W4G/UAcS+lOBN2Ae8HjcDZsIu8h4GxbNyoHcIXLOCCTCKiQpBrEJ84cSI+/fRT3Lx5M98nvXjxIoYOHYopU6a8U3FkWhLiX+OrqYt1Hrv85zXcvBph5IqoOLt25Rb+uqi7IfDV1MVIkL42ckVmSCQCPvxQc+zFC2D6dGHqISIqZnIN4ocPH0aFChXwwQcfoEOHDli0aBEuXryI+Ph4rcfGxcXh5MmTWLx4Mdq0aYPhw4fD19cXBw8eNGjxZFx3ox4gPT33LuTmdTuR9DrZiBVRcfU6MQlb1u3M9XhaWjru3fnXiBWZmbt3tbvd1aqpuuAeHsLURERUDOU6R9zW1hZTp07FBx98gE2bNmHnzp0IDQ1VH7O3t4dCoUBiYiJkMhmUSiUcHR3Rs2dPfPjhh/D09DTW10BGkpyUd8hOTkqBTCYzUjVUnMnlcqSmpOb5mOSkFCNVY2aGDQPWrdMcO3ECaNtWmHqIiIqxPG9xDwBly5bFjBkzMGnSJFy5cgVXr15FdHQ0EhISIBaL4ebmBk9PTzRu3Bh16tThFobFmF/1Snkeb9qyAey5gI4KQQn7EmjUrD4i/76b62Pedj1SPmVmAtbW2uMyGSCRGL8eIiIz8NYgns3a2hrNmjVDs2bNDFkPmTAXV2c0bdkAF87+pXXM1tYGfT7oCksrSwEqo+LGysoS/Qd3x68bduucDtUssCGcXRwFqKyY2rsX6NlTc2zyZGDJEkHKISIyF2xfk96cXZzwzbLpaNOxhcZe8l5ly2DdbytQqjTnjlLhKVXaA+t+WwGvsqXVYyKRCG07tcDcpdPg7OIkYHXFiJubdgiPjmYIJyIyAr074kQA4FHSDXOXTEXSl2Pw8kUc7O1LwNnFEW4errzRExUqC0sL1Kjph027f8SrV1IkJ6WgjGcpODo5wMFRx57WlD9PngBly2qOubgAUm5DSkRkLOyIU745OjnAq2wZ1K7nj0p+FeBe0o0hnAxCJBKpri+JAk6uJeBVtgxDeGGYMkU7hO/ZwxBORGRk7IgTUZGgUCiELqHok8sBCx1v+xkZgJWV8eshIjJz7IgTEZmDkye1Q/jHH6v2BmcIJyISRL474jKZDLdu3UJsbCwaNmwIGxsbyOVyODlx4RQRkUmqXh2IjNQcu3MHqFxZmHqIiAhAPjviR44cQatWrfDBBx9g0qRJuHv3LsLDwxEYGIiQkBBD1UhERAXx6pXqDpk5Q7hSyRBORGQC9A7i586dw6RJk1C+fHlMnToVSqUSAODt7Y0qVapg2bJl2Ldvn8EKJSKifFi4UPt29KGhqhBOREQmQe8gvmrVKvj7+2Pjxo3o3r27etzX1xdbt25FnTp1sGHDBoMUSUREelIqVV3wGTM0x5OTgQ8/FKYmIiLSSe8gHhkZiS5duui8hb2FhQWCgoLw77//FmpxRESUD5cvAznfo7t2VYXzEiWEqYmIiHKl92JNS0tLyGSyXI8nJCTA0pK3NyciEkSrVsCZM5pj168DtWoJUQ0REelB7454w4YNsXPnTmRkZGgde/HiBbZu3Yp69eoVanFERPQWr1+rpqLkDOEKBUM4EZGJ0zuIT5w4ES9fvkS3bt3w448/QiQSISwsDAsXLkRQUBCSk5Mxbtw4Q9ZKRERvWr0ayLl17PLl/80TJyIik6b31BRfX19s2bIF8+bNw6ZNmwAAmzdvBgD4+/tj1qxZqFatmmGqJCKi/yiV2nPBAdUt6l1cjF8PEREVSL5u6OPn54dNmzYhISEBjx8/hkKhgJeXFzxybpFFRESG8c8/gL+/5lijRsDFi8LUQ0REBZbvO2sCgLOzM5ydnQu5FCIiylPfvsDOnZpj584BzZoJUw8REb2TXIN427Zt830ykUiEEydOvFNBRESUQ1oaYGenPS6X656iQkRERUKu7+Cenp5av6RSKWJiYmBvb4+aNWvC398fTk5OiImJQVpaGurWrWvM2omIir9t27RD+Jw5uc8TJyKiIiPXjnj2gsxsO3fuxN9//42NGzeiYcOGGseuXbuGESNGwD/nvEUiIio4sVj7lvSxsUDp0sLUQ0REhUrvdsratWsxZMgQrRAOAHXq1MGHH36I9evXF2pxRERm6cED1faDb4bwcuVUv2cIJyIqNvQO4i9fvoSrq2uux+3s7JCYmFgoRRERma1RowBfX82xI0eAhw8FKYeIiAxH711T/Pz8sHPnTvTt2xd2OeYrSqVSbNmyBbV4FzciooLJygKsrHSPWxRogysiIjJxer+7jxs3DiNGjEBQUBCCgoJQtmxZpKen49GjR9i/fz+ysrLw/fffG7JWIqLi6fBhoEsXzbGxY4EVK4Sph4iIjELvIN60aVOsXbsWS5cuxdq1a9XjIpEI9evXx7Rp01CjRg2DFElEVGx5ewMxMZpj//4LlC8vSDlERGQ8+fp5Z9OmTbF79271NoYikQheXl5w4S2ViYjyJzYW8PTUHLO2BtLThamHiIiMTu/Fmn379sXGjRvx6tUruLq6IiAgAP7+/gzhRET5NWuWdgjfvp0hnIjIzOjdEReJRFiwYAEWLVqEBg0aoGvXrujQoQMcHBwMWR8RUfGhUAASifZ4WhpgY2P8eoiISFB6d8R/++03nDhxAuPGjUN8fDxmzpyJZs2aYcyYMThy5AgyMjIMWScRUdF29qx2CB8wQLU3OEM4EZFZytcccW9vb4wcORIjR47E/fv3cfjwYZw8eRITJ06Era0t2rVrh27duqFZs2YQiUSGqpmIqGipWxe4dk1zLCICqFZNmHqIiMgkFHhzWl9fX4waNQr16tXDjh07cOTIEezfvx/79+9HqVKlMHToUAwZMoSBnIjMl1QKuLlpj+e8bT0REZmlfAdxmUyG8+fP4+jRowgLC0NSUhJcXFwwcOBAdO3aFWKxGNu2bUNwcDAePnyIOXPmGKJuIiLT9u23wKRJmmM//wwMHy5MPUREZHL0DuJnz57FkSNHcPLkSbx+/Vo9FSUoKAjNmjWD5I25jzVr1sTTp0+xb98+BnEiMi9KJSDWsfzm9WuAi9uJiOgNegfxESNGwMLCAi1atEDXrl3Rpk0b2OSxwMjHxwclSpQolCKJiIqEq1eBevU0x9q1A37/XZh6iIjIpOkdxL/++mt06tQJTk5Oej3+m2++KXBRRERFTqdOwLFjmmN//QXUry9MPUREZPL0DuL9+/c3ZB1EREVTcrLuKScKBcDF6kRElId8LdbcuHEjfv/9d7x8+RJZWVlax0UiEU6cOFFoxRERmbRfftFefLl4MTBlijD1EBEVYdOmTcOePXs0xiwtLVGyZEm0bt0a48aN03tmRlGhdxD/6aef8P3338Pe3h4VK1aEtbW1IesiIjJturrdr17p3q6QiIj0Nn36dLi4uAAAMjIycO/ePWzfvh23bt3Ctm3bNDYIKer0DuK//fYb6tevjzVr1nARJhGZr9u3tW/EU6sWcP26IOUQERU37dq1g7e3t8ZY+fLl8fXXX+Ps2bNo3bq1QJUVPr1vcR8XF4euXbsyhBOR+Ro8WDuEnz7NEE5EZGCNGjUCANy9e1fgSgqX3kG8cuXKiI6ONmQtRESmKSNDNRVl82bNcZkMCAwUpiYiIjPy7NkzAKrtsQFg8ODBGDZsGL777jvUqVMHTZo0QVRUFADg3r17GDNmDOrXr49atWrh/fffxx9//KE+15gxY9CwYUMoFAr12KlTp+Dn54d58+ZpvO7o0aPRuXNng31degfxiRMnYtu2bQgLC4OSt2cmInOxcyeQ854J06erbtxTjOYpEhGZitevX0MqlUIqleL58+e4cOEC5s2bhxo1aqBNmzbqx129ehVHjhzBlClT0LNnT1SqVAlRUVHo378/7t27h5EjR2LChAmQyWQYMWIEDh8+DABo2bIlEhMTERkZqT7XpUuXAABXrlxRj2VlZeHixYsINGDDJdc54m3bttUay8rKwmeffQYbGxu4uLhAlGOx0rvsmhIZGYk+ffogLCwMpUuXzvVxt2/fxpIlSxAeHg5ra2sEBgZi0qRJKFWqVIFel4goVyVKAKmpmmMxMYCnpzD1EBGZgZ49e2qN2djYYOPGjbCyslKPpaamYsmSJahVq5Z6bN68eXB1dcWePXtgZ2cHABg0aBA+/PBDzJ8/H+3atUPLli0BABcvXkSNGjUAqIJ4qVKlEBUVhaSkJDg4OOD69etISUlBq1atDPa15hrEPXX8Q6NrrDDcv38fI0eOhEwmy/Nxjx8/xsCBA+Hu7o65c+fCyckJoaGhGDBgAPbu3QtHR0eD1EdEZubRI6B8ec2x0qWB2FhByjFVCoUckMshS08HoITE2gYiiQXE/EkBEb2DJUuWwN3dHYCqCRwTE4MtW7Zg4MCBWLt2LZo2bQpAFc4DAgLUz4uPj8fly5cxePBgpKenIz09XX2sffv2WLhwIW7duoV69eqhcuXKuHjxIoYNG4bExETcvn0bkydPxuLFi3H16lUEBgbijz/+gIODA+rlvGNyIco1iG/atMlgL5pNJpNh+/btWLZsGSwtLd/6+I0bN0Imk2H9+vXqDwWNGzdGp06dEBISgokTJxq6ZCIq7saPB1as0Bw7cAAIChKmHhOlkMmQGf8Kqc9iVNN0AAAi2LiXhE3J0hBbvP09nYhIl7p162rtmvLee++hQ4cO+Oabb3DkyBEAgLOzM8Ti/2ZZZ69l3LRpU645Nvb/DZUWLVrg119/hUwmw19//QWRSIR+/fphzZo1+OuvvxAYGIhz586hadOmsLDI12138iXfZ5bJZLh16xZiY2PRsGFD2NjYQC6XF2iD9fDwcCxduhTDhg1DqVKlMGvWrDwf/++//6JKlSoanXlra2sEBATgzJkzDOJEVHAyGaCrIZCZqXvczMlSkpAa+yTHqBLpr55DZGEBG/dSEIn1XoZERJQnFxcXNGrUCL///jsSExMBQGs/cblcDgAYOHAg2rVrp/M8lSpVAgAEBgZi3bp1uHnzJi5evIjq1auru99XrlyBVCpFREQEBg0aZMCvKh+LNQHgyJEjaNWqFT744ANMmjQJd+/eRXh4OAIDAxESEpLvF/f19cWJEyfw2Wef6bU5e5kyZfD8+XOtKSxPnjzhji5EVHDHj2uH7REjVJ1ehnAtiqwsVSc8F+kvn0HxlqmGRET5lb3LiTiXD/leXl4AVAG9adOmGr9KliyJzMxM2NraAgDq1auHEiVK4OLFi7hy5QoaNGgAAGjYsCH+/vtvhIWFAYB6Prmh6N0RP3fuHCZNmoS6deti+PDhCA4OBgB4e3ujSpUqWLZsGTw8PNC9e3e9Xzx7/o++evTogZ07d2LmzJkYP348rK2tsXHjRty9e/et88vzIpfL1VvekOlJ/f9iOX6PzJchr4EKHTrA6vFjjbEHx48jy8cH4DWnk4+XJxQZ6bkeV8rlyMrMwL1//y3U1+V7AQG8DooKPz+/Qj3fq1evcPHiRVSrVg0ODg46H1OyZEn4+/tjz549GD58uHojj6ysLMyYMQO3b9/GmTNnAACWlpZo0qQJTpw4gaioKIwfPx6AKohnZWVhzZo18Pf3z3dWzS+9O+KrVq2Cv78/Nm7cqBG2fX19sXXrVtSpUwcbNmwwSJHZ6tevjwULFuDkyZNo3bo1mjVrhgcPHmDAgAGwybm9GBFRHiSvXsGvalWtEB51+7YqhFOulEqo9lXP09uOExHpduLECezbtw/79u3D3r17sXr1avTr1w9paWmYMGFCns+dNWsWMjMz0bt3b6xatQpbtmzBhx9+iBs3buCzzz6Di4uL+rEtW7bEP//8AwDqBZnZQT86Otqgu6Vk07sjHhkZiQkTJuj8cYCFhQWCgoKwZMmSQi1Ol169eqF79+54/PgxHBwc4O7ujunTp8PZ2bnA55RIJIX+yY0KT3bXg98j81Xo18DcucCcOZpjW7YAH3wAXmVvp1DIoXBxQ6b0lc7jFiUcYGVjU+h/Z/leQACvA3OwcOFC9f9LJBI4OTkhICAA8+fPR5MmTfJ8bp06dbBt2zasXLkS69evh0wmQ4UKFRAcHKy1LWL2tBM/Pz/1zntisRj16tXD6dOnDT4tBchHELe0tMxz+kdCQoJeO5+8i/v37+Pvv/9G9+7dUaFCBfV4REQEqlevbtDXJqJiQKHQfROe1FTg//MG6e3EYglsS3pClvQaiqxMjWMiiQQlvMtBbMBdBoioeAoODlZPfX6bvHb3q1GjBlavXv3Wc5QpU0bnFKc1a9boVUNh0HtqSsOGDbFz505kZGRoHXvx4gW2bt1q0H0WAeDOnTv44osvNBZmXrx4Ebdv3851dSwREQDgwgXtEN67t2qeBUN4vkmsrOBYqSpsy3hDYmMLsbUNbEqWgWPl6hBbWQtdHhFRkaB3y2LixIno378/unXrhpYtW0IkEiEsLAynT5/Gnj17kJmZiXHjxhVqcVKpFI8fP0alSpVgb2+PwMBAeHt7Y+LEiRg7diykUikWLlyIWrVqoVu3boX62kRUjDRpAly8qDl26xbg7y9MPcWE2NIKNu6lYO3iBihV3XBuWUhEpD+93zF9fX2xZcsWlCxZEps2bYJSqcTmzZuxYcMG+Pj4IDQ0FNWqVSvU4k6fPo3+/furJ9Lb2dkhJCQETk5OmDBhApYuXYrOnTsjJCREr+0PicjMJCSoFhXmDOEKBUN4IRGJRBBbWEJsackQTkSUTyKlUn1LNL0lJCTg8ePHUCgU8PLygoeHhyFqM4rr168DAGrXri1oHZQ7LsyhAl0DK1cCOX9Kt2oVMHp0IVZGxsT3AgJ4HVDxUqDVNM7Ozu+0SwkRkcEolYCuzmxCAlCAOwATEREZSq5BvG3btvk+mUgkwokTJ96pICKiArtxA8j5062WLYH/38CBiIjIlOQaxD09PbXG/v77b6SlpcHPzw8VK1aEQqFAdHQ0IiIi4OrqiubNmxu0WCKiXHXvDuzfrzl28SLQqJEw9RAREb1FrkE85/6MO3fuxN9//42NGzeiYcOGGseuXbuGESNGwJ+Ln4jI2FJSAHt77XGFQo+7PxIREQlH7yXua9euxZAhQ7RCOKC6i9GHH36I9evXF2pxRER52rhRO4TPn6+aJ84QTkREJk7vxZovX76Eq6trrsft7OyQmJhYKEUREb2VrqD94gVQhHdxIiIi86J3R9zPzw87d+5Eamqq1jGpVIotW7agVq1ahVocEZGWu3e1Q7ifn6oLzhBORERFiN4d8XHjxmHEiBEICgpCUFAQypYti/T0dDx69Aj79+9HVlYWvv/+e0PWSkTmbtgwYN06zbETJ4AC7PJERERUUJGRkejTpw/CwsJQunTpAp9H7yDetGlTrF27FkuXLsXatWvV4yKRCPXr18e0adNQo0aNAhdCRJSrzEz41aypPS6TAbyrLhERGdH9+/cxcuRIyGSydz5Xvm7o07RpU+zevRtSqRQxMTEQiUTw8vKCi4vLOxdCRKTTrFnwmz9fc2zSJGDpUmHqISIio1EqlchKSkRmvBQKeRbEEktYubjC0sEJIiMvypfJZNi+fTuWLVsGS0vLQjlnge6s6erqmufCTSKiQqHrTfbxY6BsWePXQkRERqWQZSHp37uQp2muT8xMlEJiaweHCpUhtiicQKyP8PBwLF26FMOGDUOpUqUwa9asdz4n76xJRKZH1x0yAdWCTCIiKvaUSqXOEJ5NnpaKpH/vwrFSNaN1xn19fXHixAm4ublh9+7dhXLOfN1Zk4jI4MqUAZ490xiKGzECryZOhJ9AJRERkXFlJSXmGsKzydNSkZWUCCtHZ6PU5O7uXujn1PvOmkREBpWVBVhZaY8nJ+PVkyfGr4eIiASTGS/V73EJUqMFcUPQex9xfURERBTm6YjIXKxYoTuEK5VAiRLGr4eIiASlkGfp9ziZfo8zVXov1szMzMSKFSvwxx9/IDU1FQqFQn1MLpcjJSUFycnJiIyMNEihRFRM6Zrbd+YM0LKl8WshIiKTIJbotwjTmIs1DUHvjvj333+PkJAQJCYmwtbWFjExMShTpgwsLCzw7NkzZGVlYebMmYaslYiKkwcPdIdwpZIhnIjIzFm56Lc7n5Vz0d7FT+8gfvToUTRs2BAnT57Ezz//DACYPXs2jh07hjVr1kAmkxXanopEVMw1agT4+mqOffwxd0UhIiIAgKWDEyS2dnk+RmJrB0sHJyNVZBh6B/Hnz5+jQ4cOEIvFKFWqFNzc3HDt2jUAQGBgIHr27InffvvNYIUSUTGgUKi64Jcva47HxWnfup6IiMyWSCSCQ4XKuYbx7H3EjX1Tn8KmdxC3sbHR6Hj7+Pjgzp076t/XrFkT0dHRhVsdERUfW7fqvh29UgnwBmFERJSD2MISjpWqwb58JVg5u8LC3gFWzq6wL18JjpWqCTo/vFevXoiKikLp0qXf6Tx6L9asVq0azp49i/79+wMAKlasqO6IA6qOeVH/VEKmSSKRaCwOpiJI13vD3r1A9+5GL4WIiIoOkUgEK0fnIr1FYV707oh/8MEHCAsLwwcffIDk5GR06dIFERERmD59On7++WeEhoYiICDAkLWSmYmXJuDh/ceIefQCGakyvHweByXnEBctz5/rDuEKBUM4ERGZPb074u+99x6Sk5Oxfv162NraomnTphg4cCC2bNkCQHUnzmnTphmsUDIvz2JfYMqYr3Aj/B/1WGnPkvj+5/moXNUXFhY6pjiQaendG8h5C+AOHYBjx4Sph4iIyMSIlO/YYnz69CkSExPh6+sLK1035DBx169fBwDUrl1b0DroP3EvpRg5eDLuRN7XOmZXwhZ7ft+AMl6lBKiM9KarC/7kCeDlVaDTRUVFAQD8/HiTe3PG64AAXgdUvOg9NWXIkCH4888/tcY9PT1RrVo1nDt3Dl26dCnU4sg8xT59rjOEA0BqShoO7zsBuVxu5KpIL8eO5b43eAFDOBERUXGV69SUtLQ0xMfHq39/+fJltG/fHuXKldN6rEKhwNmzZ/HkyRPDVElm5ebViDyPXzx3Bf0H94C9A299blJ0BfBffgGGDjV+LUREREVAnkG8R48eSEpKAqBatbpgwQIsWLBA5+OVSiWaNWtmmCrJrLi4Oud53MnZEZaWei9vIEN7/Rpw0nFDBZlM93aFREREBCCPIO7q6oolS5bg1q1bUCqVWLVqFdq3b69zTpZYLIarqyunplChqN3AHxaWFpBlyXQeHzi0D6xtrI1cFenUty+wc6fmWPXqwD//6H48ERERqeXZVgwMDERgYCAA1aLM999/H7Vq1TJKYWS+nJ0dMWfhZHw5OVjrWOce7VC+QlkBqiItuqai3L4NcAEVERGRXvT++f7ChQvV///ixQvExsaiYsWKsLa2hoWFBcRivdd9EuXJ1s4WbTu1RJVqvvj5h024/c89eJR0w9DRHyCgdnW4uDkLXaJ5O3kSaNtWe5x7vBMREeVLvibahoeHY/78+YiMjAQArFu3DnK5HDNmzMC0adPQuXNngxRJ5sfeoQSq+VfBN0un48WLlxCLRShX3kfoskhXF/zzz4HvvjN6KURERMaiUCiwfft2bN26FU+ePIGbmxvatm2LsWPHwt7evsDn1buNffPmTXz88cdISUnBhx9+qB53cnKChYUFJk+ejDNnzhS4ECJdStjbITMrHekZaUKXYt7S0nSH8LQ0hnAiIir2QkJC8M0336BVq1ZYtWoVPv74Y+zduxfjx49/p/Pq3RH//vvv4e3tjd27dyM1NRWhoaEAgICAAOzfvx8DBgzAmjVr1HPKiaiYGDsW+OEH7XFORSEiIgOTy+U4f/oyDu49DumrBLi6OyOoRwc0a9UQEiPtzKVUKhESEoL+/ftj0qRJAICmTZvCxcUFEyZMQGRkJKpVq1agc+sdxK9du4bRo0fDxsYGaWma3Ul7e3v069cPK1asKFARRGSidHXBz5wBWrY0fi1ERGRW4l7FY8xHXyDi1h2IxWIoFAqIxWIc3X8S1QOqYFXoYri5uxi8jpSUFHTr1g3vvfeexnjFihUBAI8fPy5wEM/XCsu8bmGfkZEBhUJRoCKIyMRcu5b7HTIZwomIyMDkcjnGfPQFIv+5CwDqjJn938h/7mLMR18Y5U7b9vb2mDVrFurVq6cxfuLECQBApUqVCnxuvYN4rVq1cPDgQZ3HUlNTsWPHDgQEBBS4ECIyESIRULeu5ljfvpyKQkRERnP+9GVE3LoDpUL3vz1KhRIRt+7gwpm/jFyZyo0bN7B27Vq0a9cOvr6+BT6P3kF83LhxiIiIwKBBg7B3716IRCLcvHkTGzduRPfu3fHkyRN8+umnBS6EiAQmk+nugickAL/9ZvRyiIjIfB3ce/ytW2OLxWIc3HPcSBX9Jzw8HMOHD4e3tzfmzZv3TufSO4jXqVMHa9aswbNnz7Bo0SIolUp89913WLBgAdLT0/Htt9+icePG71QMEQlk/nzA0lJ7XKnUfft6IiIiA5K+SnjrlGeFQgFpXIJxCvq/w4cP4+OPP0aZMmUQGhoKF5d3m6Oer33EmzVrht9//x0RERF4/PgxFAoFvLy84O/vDwuLfJ2KiEyFri747t1Az57Gr4WIiAiAq7uzeoFmbsRiMVyNeJO/9evXY9GiRWjYsCFWrVoFBweHdz5nvhZrPn36FMuWLYO3tzfee+89dOnSBZcvX8ayZcsQFxf3zsUQkRHdu5f7gkyGcCIiElBQjw56dcSDenYwSj07duxAcHAw3nvvPYSEhBRKCAfyEcTv3LmDnj17Yv369YiNjVWPv379Glu2bEGPHj0QHR1dKEURkYGVLQtUrqw51qABF2QSEZFJaNaqIaoHVIFIrKNhBEAkEqF6QBU0DWxg8Fri4uIwf/58eHl5YeDAgYiIiMD169fVv6RSaYHPrXcQX7ZsGUqUKIFDhw6hatWq6vHJkyfj0KFDsLS0xNKlSwtcCBEZgVKp6oI/eaI5HhMDXL4sTE1EREQ5SCQSrApdjGo1VE2j7IWb2f+t5l8Zq0IXG+WmPn/88QfS0tIQExODgQMHon///hq//vjjjwKfW++J3devX8eYMWNQvnx5rWNly5bFoEGDEBISUuBCiMjAfvkFGD5ce5xdcCIiMkFu7i7Ysm81Lpz5Cwf3HIc0LgGubs4I6tkBTQMbGO3Omj169ECPHj0Mcm69g7hCoUB6enqux5VKZZ7HiUhAuuaCr14NjBxp/FqIiIj0JJFI0KJNY7RoUzx35tN7akrt2rWxfft2vH79WutYSkoKduzYgVq1ahVqcUT0jp490x3CFQqGcCIiIoHp3RH/7LPPMGjQIAQFBaFr164oV64cRCIRHj9+jEOHDuHFixdYuHChIWslovxo3hw4f15zzN0dePlSmHqIiIhIg95BvFatWur9E3/55ReNY1WrVkVwcDDq1KlT6AUSUQHo6oJHRQFVqhi/FiIiItIpX3fhqV+/Pnbs2AGpVIqYmBgoFAqUKVMGALB//37MnTsXBw8eNEihRKSHv/4CGjbUHueCTCIiIpNToNthurq6wsHBAWFhYfjxxx9x/vx5yGQyo61eJSIdqlZVdb3fNHcu8OWXwtRDREREecp3EP/777+xe/duHDp0CK9fv4ZSqYS7uzt69+6N/v37G6JGIspLaipQooT2uEwG8MMxERGRydIriMfFxWHfvn3Ys2cP7t27B6VSCdH/56COHTsWI0eOhIVFgZrrRPQuvv8e+PxzzbGgIODAAUHKISIiIv3lmp5lMhlOnjyJ3bt349y5c5DJZLCyskJgYCDat28PPz8/9OnTB1WrVmUIJxKCrgWZ9+4Bvr7Gr4WIiIjyLdcE3bx5cyQmJsLe3h7t27dH+/bt0bJlS9jb2wMAYmJijFYkEb0hMhKoXl17nAsyiYiIipRcb+iTkJAAW1tbdO3aFZ06dULjxo3VIZyIBNKqlXYIX72aIZyIiKgIyrUjvmHDBhw4cAAHDx7Etm3bIBKJULt2bXTo0AHt27c3Zo1ElJUFWFlpj6enA9bWxq+HCIBSoYBCLoMiMwNQKiG2soFIIoGYi4TNikIug1ImhzwzA2KJBGJLS4gsLNVryUg3RVYWlHIZ5FlZEFtYqH5Z6nifp2It1yDeqFEjNGrUCLNnz8aZM2dw4MABnDlzBlevXsWiRYtQvnx5iEQipKamGrNeIvOzaRMwZIjmWN26QHi4MPUQAVDIZMhMkCI19gmgVKgGRSLYuJeGjUdJiC0shS2QjEKRlYWUp4+QlZigHhNZWMC+bEVYlLCHSJzrD97NmjwzA8mP7kOe9l+GEltbw75cJUisbfghxoy8dZWllZWVeo54cnIyjh07hoMHD+LSpUtQKpWYOnUqdu/ejT59+qB9+/aw0tW1I6KC0fVmfP06UKuW0UshepMsNRmpTx9rDiqVSH8ZC5GFBDbuJSESMYQVZwq5DCkxj5D1OkFjXCmTIenfu3CqUh0SG1thijNhiqxMJD24o/pJ0pvjGRlIuh8Fx8rVIWGWMhv52u7E3t4evXv3Ru/evfHy5UscOnQIBw4cwJ9//omLFy/C0dERly5dMlStRObj8WOgXDntcc4FJxOgyMpSdcJzkf7iGaycXSHhj9mLNaVMrhXC3ziKtJfPYOfpw6lKOcgz0rVCeDalXIas5ERIXD2MXBUJpcDtCg8PD3z00UfYtWsXjh49itGjR8PZ2bkQSyMyUwMGaIfwBQsYwslkKJUKKDLScz8ulwEKhRErIiHIM3O/BgBAlpIMJa8DLVnJSXkff53IPzczUigbgJcvXx5jx47F2LFjC+N0ROZJodB9J8zXrwEHB+PXQ5QLkUikmjaV14dDznEt9kRv6XSLJBJeBjqI33LvFZGFBf/+mBFO4CMyBQcOaIfwMmVUQYchnEyNRAIrZ7dcD6sW6XE6QnEntrSCSJJ7qLRx46JdXSwdnPM8buNWkos1zQiDOJHQRCKgWzfNsT/+AJ4+FaYeorcQiyWwLVUGIkvtkCUSS1DCq9xbu35U9IktLFDCp4LOYxZ2JWDp4GTkiooGkYUFbMt46Txm7erBLQzNDN8piYTy8iVQsqT2uELBH0uSyZNYWcOpUjVkSOOQkRAHKJWwcnKGtVtJBgkzIRKJYVnCAY5VaiD9RSxkKckQSSSwdi8FKwcniHV8UCNALJHA2sUDFnb2SH8eC3lGOsSWVrApWQYWdnb8EGtm+N0mEsL48cCKFZpjEyYA334rTD1EBaAKD6Vg7eoGJUQQSyTcN9rMiMRiWNjYws6rHJQKOUQiEaej6EF1Ax8HSHxsVfvwi8QM4GaK33UiY1IqAV1B5eVLwN3d+PUQvSORSAwRO+BmTyyR6F5sTnli+Ca2LoiM5dEj3SFcqWQIJyIiMkMM4kTGsHgxUL685tihQ9wbnIiIyIzxZyJEhpSRAdjaagbuKlWA27e5IJOIiMjMsSNOZChnzwI2NpohfNcuICqKIZyIiIjYEScyiF69gD17NMcSEgAn7qtLREREKibTEY+MjESNGjXw7NmzPB8nlUoxffp0NG/eHA0bNsTIkSPx8OFD4xRJ9Daxsapu95shfNw4VVecIZyIiIjeYBJB/P79+xg5ciRkMlmej1MqlRgzZgzOnj2LyZMnY/HixXj58iWGDBmCxMREI1VLlIuffgI8PTXHbtwAvv9emHqIiIjIpAk6NUUmk2H79u1YtmwZLPW4A9fDhw9x9epVLFq0CD169AAA+Pr6ol27djh58iR69uxp4IqJdJDJgNKlgbi4/8YqVADu3uW+ukRERJQrQTvi4eHhWLp0KYYOHYrJkye/9fEZGRkAgBIlSqjHnP7/4/6EhASD1EiUp7/+AiwtNUP4hg3AgwcM4URERJQnQTvivr6+OHHiBNzc3LB79+63Pr5q1apo1KgRVq1ahYoVK8LFxQXBwcGws7NDu3btClyHXC5HVFRUgZ9PhpWamgoAJvc9KjVzJpx37dIYu/fnn5C7uKh2RqFCY6rXABkXrwMCeB0UFX5+fkKXUCQIGsTdC3A3wa+++grDhw9H586dAQBWVlZYtWoVypYtW9jlEekkjo9H5SZNNMYSe/TAs+BggSoiIiKioqhIbV94//59vP/++/Dx8cGMGTNgY2OD3377DePGjUNISAjq169foPNKJBJ+cjNh2V0Pk/gebd4MDB6sOXbpEpwaNgT3RDEck7oGSDC8DgjgdUDFS5EK4qGhoQCAdevWqeeGN2vWDB988AEWLFig1/QWogJRKFR3xLx//78xJyfg5UvVHHEiIiKifDKJ7Qv19fTpU/j6+qpDOACIRCLUq1cP9+7dE7AyKtb+/lu18PLNEP7DD6ob9DCEExERUQEVqSBeoUIF3L17F69fv9YYv3HjBry8vASqioq1SZOAgADNsZgYYMwYYeohIiKiYsOkg7hUKsX169eRnJwMAPjoo48gkUgwdOhQHD9+HGfPnsWkSZNw+fJljB07VuBqqVhJSlLdIfPbb/8b69pVdYfMnDftISIiIioAkw7ip0+fRv/+/fHPP/8AALy9vbFt2zZ4eHhg+vTpmDhxImJjY7F+/Xr1LipE72zvXsDRUXPs5Elg/35ByiEiIqLiSaRUKpVCFyGk69evAwBq164taB2UO6OtkFcqgcaNgcuXNcfT0gAbG8O+NuWJuyQQwOuAVHgdUHFi0h1xIqO5dw8QizVD+IIFqnDOEE5EREQGwCBONHcuULmy5tiDB8D06cLUQ0RERGahSO0jTlSo0tIAOzvNsRYtgDNnVAs1iYiIiAyIHXEyT7//rh3CDx4Ezp5lCCciIiKjYEfcCJKTUpD0Ohkx0bGwtbOFRyk3uLq5wMJCInRp5qljR+D4cc2xpCTA3l6YeoiIiMgsMYgbWNyreKxYvBb7dx6DXC4HALi4OmHeshmo17gW7OxsBa7QjDx5ApQtqzk2dSoQHCxMPURERGTWODXFgFJSUrHq23XYs/2wOoQDQLw0EWOHTceTR08FrM7MfPeddgiPjGQIJyIiIsEwiBvQ68Qk7N1+SOcxhUKBH79dh6SkFCNXZWYyM1XbD06c+N+Yvz8glwNVqwpXFxEREZk9BnEDinsphUwmz/V4+OWbSEtNM2JFZub8ecDaGsjI+G/s11+BW7dUe4YTERERCYhzxA3I1jbvG8HYO5SAmDt0GMaAAarQ/ab4eMDZWZByiIiIiHJiW9CAHJ0dUcarVK7H3x/SAy5uzsYryBy8eKHafvDNEP7pp6o7ZDKEExERkQlhEDcgdw9XLFw+CxaW2j94qOxXEZ17tIdEwi0MC01ICFAqxwefq1eBn34Sph4iIiKiPHBqigGJRCLUqOWH3cdD8cuPW3DpfDjsStii/+AeaPdeIDxKugldYvEglwPe3sCzZ/+NeXoCjx8D/KBDREREJopB3MCsra1RvmJZzJg7HsnJqRCLRXBxdYaYiwULx7VrQN26mmMhIcCwYcLUQ0RERKQnBnEjsbWzhS1v3lO4Pv0UWLNGc+z5c6BkSWHqISIiIsoHBnEqeuLjAVdXzbEBA4CtW4Wph4iIiKgAOD+Cipbt27VD+PnzDOFERERU5LAjTkWDQgHUqAFERPw3ZmsLJCQAVlaClUVERERUUOyIk8mzevAAftWra4bw774DUlMZwomIiKjIYkecTJtIhAo5x6KjVdsVEhERERVh7IiTaXr+XHWHzDd17Ki6QyZDOBERERUDDOJkekaNAkqX1hh68cUXwNGjAhVEREREVPg4NYVMh1IJ6LjR0d1Ll6BwcgJ3ByciIqLihB1xMg2nTmmHcEtLQKmEwslJmJqIiIiIDIhBnITn5AS0aaM5dvw4kJkpTD1ERERERsCpKSScxETA2Vl7XKHQXqhJREREVMywI07CmDNHO4QPHaqaJ84QTkRERGaAHXEyPl1BOyYG8PQ0fi1EREREAmFHnIznyhXdIVypZAgnIiIis8MgTsZRrRrQoIHm2G+/qUI4ERERkRni1BQyrLQ0wM5Oe1wmAyQS49dDREREZCLYESfDWblSO4R36aLqgjOEExERkZljR5wMQ9dc8Hv3AF9f49dCREREZILYEafCFRWV+4JMhnAiIiIiNQZxKjxt2wJVq2qOrV7NBZlEREREOnBqCr27rCzAykp7PC0NsLExfj1ERERERQA74vRuNm/WDuG1a6u64AzhRERERLliR5wKTtdc8GvXVEGciASjVCiglMugyMqCUqmE2NIKIokEYu5WRERkUhjEKf8SEwFnZ+1xzgUnEpxCLocsKREpMY+glMtVgyIRbNxLwsajNMQWlsIWSEREapyaQvmza5d2CJ83jyGcyETI09OQ/PjBfyEcAJRKpL98jvRXL6BUKIQrjoiINLAjTvpRKoF69VRTT96Ung5YWwtTExFpUMhkSIt9kuvx9FfPYe3qAYmuxdVERGR07IjT2925A4jFmiF88WJVOGcIJzIdSgVkqcm5H1cooFTIcz9ORERGxSBOefvyS8DPT3Ps4UNgyhRByiGivInEeS/IFOlaZE1ERILg1BTSLTUVKFFCc6xtW+D333XvlkJEghNJLGDt6o70V891HpfY2L41qBMRkfGwI07ajh7VDuFHjgAnTjCEE5kwkVgMG49SEFvqmAMuEqOEd3mILblrChGRqWBHnP6jVKq63qdOaY4nJ2sHcyIySWJLKzhWqooMaRwy4l9BqVTA0t4JtiVL6w7oREQkGAZxUnn0CChfXnNs1izgm28EKYeICk5saQWbkqVh5eoOEZSARAIxp6QQEZkcBnECli7VXnwZFQVUqSJMPUT0zkQiESSchkJEZNIYxM1ZRoZqysmbN/6oUwcID+dccCIiIiID42JNc/XHH4CNjWYI37kTuHqVIZyIiIjICNgRN0d9+qhuVf+mhATAyUmQcoiIiIjMETvi5uTZM1W3+80QPnasarcUhnAiIiIio2IQNxdr1gBlymiO3bgBrFghTD1EREREZo5TU4o7mQwoXRqIi/tvrFw54P59QMLtzIxNIZcDCgUgFkPMP38iIiKzxo54cXblCmBpqRnCQ0OBhw8Zwo1MIZMhKyUZKY8f4PWDKKQ8foCs1GQoZDKhSyMiIiKBsCNeXA0fDvzyi+bYy5eAu7sw9ZgxhUyG9LiXSH8e899YRjqykhJhW8oL1m4eEFvwryIREZG5YUe8uJFKVQsy3wzhgwerFmQyhAtCKZdphPA3pT2PgVLOrjgREZE5YhAvTrZsAdzcNMcuXgQ2bhSmHgIAZMTHveW41EiVEBERkSnhz8OLA4UCqFoVuHv3vzEnJ9VUFN7iWnCKrMy8j8vyPk5ERETFEzviRd0//6gWXr4ZwleuVN2ghyHcJFg65L1Hu6W9o5EqISIiIlPCjnhRNmkS8O23mmMxMYCnpzD1kE6WJewhkljonAsusrCAZQl7AaoiIiIioTGIF0VJSYBjji5qUBBw4IAw9VCeRBaWcPStiqSHd6HIzFCPi62s4VC+MkQW/MkFERGROWIQL2r27QN69NAcO3kSaN1akHLo7UQiESQ2NnD09YMiKwuKrEyILa3+/4shnIiIyFwxiBcVSiXQpAlw6ZLmeFoaYGMjTE2UL9nhGyghdClERERkArhYsyi4fx8QizVD+Pz5qnDOEE5ERERUJDGIm7q5c4FKlTTHHjwAZswQph4iIiIiKhScmmKq0tIAOzvNsRYtgDNnVHfOJCIiIqIijR1xU3TihHYIP3gQOHuWIZyIiIiomGBH3NR06gQcO6Y5lpQE2HOvaSIiIqLixGSCeGRkJPr06YOwsDCULl1a52OmTZuGPXv25HqOqKgoQ5VneE+eAGXLao5NnQoEBwtTDxGRiRHxJ4JEVMyYRBC/f/8+Ro4cCZlM+86Dbxo9ejTef/99jbFHjx5h2rRp6NevnyFLNKzly4EJEzTHIiKAatUEKYeIyJQoFQooZDL4lC4FAJBnZkJsYQGRmLMriahoEzSIy2QybN++HcuWLYOlHjc28fHxgY+Pj/r3crkc33zzDapWrYqZM2caslTDyMwEnJyA9PT/xmrUAG7eVG1XSERk5hQyGTIT4pD6LAZQKFSDYjHsSnvBytkNYguT6CcRERWIoGkvPDwcS5cuxdChQzF58uR8P//XX39FREQEvv76a1hZWRmgQgO6cAGwttYM4b/+Cvz9N0M4EREApVKJrKREpD6N/i+EA4BCgdSn0chKToRSqRSuQCKidyRoK8HX1xcnTpyAm5sbdu/ena/npqSkYMWKFejevTtq1qxpoAoNZMAAVeh+k1QKuLgIUw8RkQlSyrKQ9iwm1+Npz2JgWcIRIj1+okpEZIoEDeLu7u4Ffu6uXbvw+vVrjBw58p3rkMvlRlnoKYmLQ6VmzTTGEt5/H8+/+gp48UL1i7SkpqYCKOKLcemd8BowTz5enlBkZeZ6XJGZiYyMdDx+8MCIVZHQ+H5QNPj5+QldQpFQZCfXbdmyBW3btkWFChWELkUvouRkrRD+cPduZFSvLlBFRESmTQTukkJExVuRDOK3b9/Gw4cPCzSvXBeJRGL4T25Xrvz3/2XKANHRKC+RGPY1i4nsrgc/XZsvXgPmSSGXIcveEVnJr3Uet3RwgrWNLa8LM8P3AypOiuSqwNOnT8POzg6BgYFCl6K/evWAkyeB69eBp08BhnAiojyJJRaw8/KBSMf7pUgigZ1nWe6aQkRFWpF8B7t+/Tr8/f2L1k4pIhHQurXQVRARFSliK2s4Vq6O9FcvkJkoBQBYObnCxr0kxJZF6N8AIiIdTLojLpVKcf36dSQnJ2uM37lzB5UqVRKoKiIiMhaRSASJlTXsSnvBoowPLMr4wK60FyRW1rzTJhEVeSYdxE+fPo3+/fvjn3/+0RiPi4uDo6OjQFUREZGxicRiPH3+Ak+fv+AdNYmo2DCZqSm9evVCr1693joGADdu3DBWWUREREREBsG2AhERERGRABjEiYiIiIgEwCBORERERCQABnEiIiIiIgEwiBMRERERCYBBnIiIiIhIAAziREREREQCYBAnIiIiIhIAgzgRERERkQAYxImIiIiIBMAgTkREREQkAAZxIiIiIiIBMIgTEREREQmAQZyIiIiISAAM4kREREREAmAQJyIiIiISAIM4EREREZEAGMSJiIiIiATAIE5EREREJAAGcSIiIiIiATCIExEREREJgEGciIiIiEgADOJERERERAJgECciIiIiEgCDOBERERGRABjEiYiIiIgEwCBORERERCQABnEiIiIiIgEwiBMRERERCYBBnIiIiIhIAAziREREREQCYBAnIiIiIhIAgzgRERERkQAYxImIiIiIBMAgTkREREQkAAZxIiIiIiIBMIgTEREREQnAQugCiIio6FIqlVDK5QCUEEksIBKJhC6JiKjIYBAnIqICUWRlIjMpEZnSV1AqlbBydoG1kxvEVlZCl0ZEVCQwiBMRUb4psjLx+n4UFJkZ6rG0tFSkv3wOR9+qkFjbCFgdEVHRwDniRESUL0qFAmkvYjVCuPqYTIaUmMdQyGUCVEZEVLQwiBMRUb4o5DJkSONyPS5Lfv3/eeNERJQXBnEiIsofJQCl4i2PURqlFCKiooxBnIiI8kUkEkFsZZ37cbEEEPOfFyKit+E7JRER5YvY0hK2pTxzPW7tXgpiCfcCICJ6GwZxIiLKN0sHJ9iUKgNAc99wK1d32Lh7QMSOOBHRW7FlQURE+Sa2sICNe2lYu7hDnpYKpVIJC9sSEEkkEFvwnxYiIn3w3ZKIiApELJEAEgkkecwXJyKi3PFnh0REREREAmAQJyIiIiISAIM4EREREZEAGMSJiIiIiATAIE5EREREJAAGcSIiIiIiATCIExEREREJgEGciIiIiEgADOJERERERAJgECciIiIiEgCDOBERERGRABjEiYiIiIgEwCBORERERCQABnEiIiIiIgEwiBMRERERCUCkVCqVQhchpPDwcACARCIRuBLKjVwuB8DvkTnjNUAArwNS4XVQNEgkEgQEBAhdhsmzELoAorfhmy3xGiCA1wGp8Dqg4sTsO+JERERERELgHHEiIiIiIgEwiBMRERERCYBBnIiIiIhIAAziREREREQCYBAnIiIiIhIAgzgRERERkQAYxImIiIiIBMAgTkREREQkAAZxIiIiIiIBMIgTEREREQmAQZyIiIiISAAM4mQyIiMjUaNGDTx79izPx718+RKzZs1C69atUadOHfTq1QtHjhwxUpVkaPpeBy9evMDkyZPRpEkT1K1bF6NHj8ajR4+MVCUZkr7XwJtiY2NRr149/PjjjwasjIxJ3+tg37598PPz0/o1d+5cI1VKVHAWQhdABAD379/HyJEjIZPJ8nxcZmYmhg8fjqSkJIwbNw4lS5bEsWPH8Pnnn0MulyMoKMhIFZMh6HsdZGRkYPjw4cjIyMDs2bNhY2ODVatWYdCgQTh06BAcHR2NVDEVNn2vgTcplUrMmDEDycnJBqyMjCk/18Ht27dRrlw5LF68WGPc3d3dUOURFRoGcRKUTCbD9u3bsWzZMlhaWr718WfPnsXt27exY8cO1KxZEwDQrFkzPH36FD///DODeBGV3+vg1KlTiIqKwq5du+Dv7w8AqFy5Mtq2bYtjx46hb9++hi6ZCll+r4E3bd26FQ8ePDBQZWRMBbkOoqKiUKNGDdSuXduwxREZAKemkKDCw8OxdOlSDB06FJMnT37r40uUKIH+/fsjICBAY7xixYp4/PixocokA8vvddC8eXNs27ZNHcIBqP/RzszMNFidZDj5vQayRUdHY+nSpfjmm28MWB0ZS0Gug9u3b8PPz8/AlREZBjviJChfX1+cOHECbm5u2L1791sf36RJEzRp0kRjLCsrC2fOnEHlypUNVSYZWH6vA3t7e9StWxeA6vt///59LFq0CC4uLmjfvr2hyyUDyO81AAAKhQLTpk3De++9h5YtWxq4QjKG/F4HL168QFxcHCIiItCpUydER0fD29sbo0aNQo8ePQxfMNE7YhAnQRXGHL4lS5bg4cOHWLVqVSFUREJ4l+tg7NixOHXqFMRiMebPn4+SJUsWYmVkLAW5BjZs2IAnT55g9erVBqiIhJDf6+D27dsAgCdPnmDKlCmwtrbG3r17MXXqVMjlcvTu3dsQZRIVGgZxKrKUSiWWLFmCDRs2YNiwYWjXrp3QJZEAPvnkE3z44YfYv38/pk+fDgDo1auXwFWRod2/fx/Lly/HihUr4ODgIHQ5JBB/f3+sXr0aDRo0gL29PQDV1LW4uDh8//33DOJk8hjEqUjKzMzEtGnTcOjQIQwbNgxffPGF0CWRQOrVqwdANW0pJiYGa9asYRAv5uRyOaZPn45OnTqhWbNmGjtrKBQKyGQyWFjwnzdz4OrqitatW2uNBwYG4sKFC5BKpXB1dRWgMiL9cLEmFTnJycn4+OOPceTIEcyYMYMh3AxFRETg0KFDWuM1atTAixcvBKiIjCk2NhY3btzA3r17UaNGDfUvAFi5cqX6/6n4u3btGnbs2KE1npGRAQsLC/60hEweWwZUpMjlcowaNQo3btzAd999h06dOgldEgng4sWLWLx4MQICAuDj4wNAdW1cvHgRVapUEbg6MrSSJUti586dWuN9+vTBgAEDOB3BjFy/fh3BwcEICAhA1apVAah+KnLs2DHUrVs331thEhkbgziZNKlUisePH6NSpUqwt7fHr7/+isuXL6N///4oXbo0rl+/rn6sSCRCrVq1hCuWDCbnddCrVy9s2rQJo0aNwtixY2FjY4MtW7bgzp07WLdundDlkgHkvAZybmGarWTJkrkeo6JP13vBxo0b8dlnn+Hzzz9HiRIlsHXrVty5cwdbtmwRulyit+LUFDJpp0+fRv/+/fHPP/8AAI4dOwYA2L59O/r376/xa8CAAUKWSgaU8zpwdnbG5s2bUaVKFcydOxfjx49Heno6NmzYgEaNGglcLRlCzmuAzFPO68DJyQmbN29GzZo1sXDhQnz++edITU1FaGgoGzNUJIiUSqVS6CKIiIiIiMwNO+JERERERAJgECciIiIiEgCDOBERERGRABjEiYiIiIgEwCBORERERCQABnEiIiIiIgEwiBOR2VEoFDh8+DCGDRuGwMBA+Pv7o3nz5hg/fjyuXbum9fg2bdpg8ODBAlRaOPJT//Pnz9GoUSNER0cDAFJSUjBu3DjUqlUL3bp1w5UrV7Ses3XrVnTo0AEymUzr2MyZM7Fw4cJ3+wKIiIopBnEiMitJSUkYMWIEJkyYgMzMTAwePBhfffUV+vXrh1u3buH999/H5s2bhS5TMPPnz0eXLl1QtmxZAMCaNWtw4cIFfP755/D09MTo0aPx+vVr9eMzMzOxZs0ajB49GhYW2jdrHjNmDLZv347bt28b7WsgIioqeIt7IjIrc+bMwblz5xAcHIyePXtqHBs5ciQ+/fRTBAcHo1mzZqhQoYJAVQrjr7/+QlhYGE6cOKEeO3z4MAYMGICPP/4Yffv2RePGjXHmzBl07doVgOoutzY2Nurf5+Tp6YkuXbpg4cKF2LBhg1G+DiKiooIdcSIyG+Hh4Th06BC6d++uFcIBwNraGl999RWysrKwe/duASoUVmhoKOrVq4cyZcqox54/fw5vb28AgL29PZydnfHs2TMAQEZGBtauXYsxY8ZAIpHket6+ffvi4sWL7IoTEeXAIE5EZuPAgQMAgBEjRuT6mHLlyiE0NBSjRo3SeXzZsmXw8/PDvXv3NMYVCgVatGiBcePGqcdu3LiBTz75BPXr10ejRo0wYsQIREVFaTzv2rVr+Pjjj1GnTh3UqVMHQ4cOxc2bNzUe06ZNG8yaNQszZsxAzZo10bJlS0ilUgDAlStX8NFHH6mfP2TIEPz111/6/6H8X2xsLE6dOoV27dppjLu4uCApKUn9NSYnJ8PFxQUAsG3bNtjb2yMoKCjPc9euXRulS5fGli1b8l0XEVFxxiBORGbj8uXL8PDwgK+vb56Pa9KkCezs7HQeyw6dR44c0Tr3ixcv1FM0rly5goEDB+L+/fsYPnw4Ro0ahXv37mHIkCF48uQJAOD8+fMYPHgwkpKSMH78eIwaNQpPnz7FwIEDtRZFHjp0CFFRUZgxYwb69esHV1dXhIWFYfDgwYiNjcWoUaMwatQoxMbG4qOPPkJYWFi+/mz++OMPyOVytGrVSmO8QYMG2L17N+7du4fQ0FBkZWWhYcOGSE9PR0hICMaMGQOx+O3/lDRo0ABnz57NV01ERMUd54gTkdl49uwZKlasqDWelpaGtLQ0jTGJRAInJyetx/r5+aFKlSo4evQoxo4dqx4/fPgwHBwcEBgYCABYtGgRnJ2dsWvXLnUHOTAwEJ07d8bWrVsxefJkzJkzBwEBAdi8ebN6asegQYPQo0cPzJs3D3v37lWfPz09HT/++CNKlSoFAJDJZJg7dy5KlSqFXbt2wd7eHgDw/vvvIygoCF9//TVatmwJS0tLvf5swsPDYWdnp16kme3zzz/H0KFD0aVLF4jFYkydOhU+Pj5Yt24dHB0d0blzZ73OX6VKFRw4cADR0dFar0FEZK7YEScis6FUKqFUKrXGV6xYgSZNmmj80jWHPFvXrl1x79493LlzB4AqFB87dgzt27eHlZUV4uLicPPmTXTt2lUdwgGgQoUK2LVrFz755BNEREQgOjoa7dq1Q2JiIqRSKaRSKdLT09G6dWtERkbi+fPn6uf6+PioQzgARERE4NmzZxg4cKA6hAOAo6MjBg0ahOfPn+Pvv//W+88mOjoaXl5eEIlEGuNly5bFoUOHsGvXLpw9exYfffQRUlNTERISgs8++wxisRh79+5Fp06d0Lp1a6xcuRIKhULr/NnhO/unAURExI44EZmRkiVL4tWrV1rj/fv3R4sWLdS/Dw4ORnJycq7nCQoKwrfffoujR4+iSpUqOH/+PBISEtTTUmJiYgCo5pvnVL16dQDAn3/+CQBYvHgxFi9erPN1nj59qg7fbm5uGseyA62unV2yu/5Pnz5FnTp1cv063pSQkABHR0edx6ysrODv76/+/ebNm+Hm5ob33nsPd+/exbRp0/Dll1/Cx8cHkyZNQqlSpdCvXz+Nc2R/WIiPj9erHiIic8AgTkRmo06dOtizZw8ePXqkEZLLly+P8uXLq3/v5OSUZxD39PRE3bp1ceTIEYwbNw5HjhyBh4cHGjVqBADqjnDO7vKbsh8zfvx41K5dW+dj3pxGk3NXEl2d/ZzH9J2WAgBisVhnJzunlJQU/PLLL/j6668hEolw9OhR+Pj4YODAgQCATp064dChQ1pBPPvcee2uQkRkbjg1hYjMRrdu3QCotul7V0FBQXjw4AEePHiAU6dOoVOnTuqQmb393+PHj7Wet2TJEqxduxZeXl4AADs7OzRt2lTjl729PeRyOWxsbHJ9/eznP3jwQOvYv//+CwAoXbq03l+Pm5sbEhIS3vq4TZs2oWTJkujYsSMA4NWrV3B1dVUfd3Z2xosXL7Sel33unJ19IiJzxiBORGajadOm6Ny5M7Zt24atW7fqfMz+/fv1mlv93nvvwdLSEitXrkRCQoLGFn6lSpVC1apVcejQIY3OenR0NDZu3IhXr17B398fHh4e2LRpE1JSUtSPSU5Oxueff47p06fn2T2uUaMGPDw8sG3bNo3XSE5OxtatW+Hh4aExneRtPD098eLFC8jl8lwfk5ycjPXr12Ps2LHqbr+HhwdiY2PVXfgnT55ozGXPlj3f3dPTU++aiIiKO05NISKz8s0330Aul+Prr7/Gzp070a5dO3h4eOD58+c4fvw4oqKi4O7ujmnTpuV5HhcXFzRr1gyHDx9G2bJltaaXTJ8+HcOHD0fv3r3Rt29fiMVibN68GY6Ojvjkk09gaWmJWbNmYcKECejVqxf69OkDa2tr7NixA0+fPsXSpUt13jI+25vP7927N/r06QMA2LlzJ168eIEVK1bota1gtsaNG2P37t24e/cuqlatqvMxoaGhKF26NNq3b68ea9euHX744QfMnDkTZcuWxfHjxzFnzhyt516/fh3lypVjECciegODOBGZFXt7e6xYsQKnTp3Crl27sGPHDrx8+RL29vaoVq0aZs+ejV69esHW1vat5+ratStOnz6NLl26aB1r3LgxNmzYgBUrVmDVqlWwtrZGgwYNMGXKFHh4eABQzad2cnLCTz/9hB9//BFisRiVK1fGTz/9hNatW7/19bOf/+OPP2LVqlWwsLBArVq1MH/+fNSvXz9ffy4tWrSAWCzGlStXdAbxpKQkbNiwAQsXLtSY+161alXMnz8fq1atwqlTpzBs2DD07t1b47kKhQLXr1/Xe6tDIiJzIVLmteKHiIhydfjwYUyYMAGHDx9+602CioIxY8ZAKpVi27ZthXre8+fPY+jQodi3b1+u3XYiInPEOeJERAWgVCrx66+/olatWsUihAPA0KFDcfXqVZ2LTN/F3r170axZM4ZwIqIcGMSJiPJBJpNh3Lhx6NevHy5duoThw4cLXVKhqVevHlq3bo21a9cW2jmjo6Nx7NgxfP7554V2TiKi4oJBnIgoHywsLPDo0SM8ePAAn332GTp06CB0SYVqzpw5OH78eKF1xX/88Uf0798fNWvWLJTzEREVJ5wjTkREREQkAHbEiYiIiIgEwCBORERERCQABnEiIiIiIgEwiBMRERERCYBBnIiIiIhIAAziREREREQC+B+wdEUhlbLgkQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 759.85x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Preliminary visualisation of data (not required for LDA)\n",
    "\n",
    "# Plot data\n",
    "\n",
    "sns.set_style(\"whitegrid\")\n",
    "sns.set_context(\"talk\")\n",
    "\n",
    "g = sns.relplot(x=\"glycerol%\", \n",
    "                y=\"acetaldehyde%\", \n",
    "                data=brew_data, \n",
    "                kind=\"scatter\", \n",
    "                hue=\"Brew\", \n",
    "                height=8, \n",
    "                aspect=1.2)\n",
    "ax = g.axes[0,0]\n",
    "\n",
    "# Plot suggested binary classifier\n",
    "\n",
    "ax.plot(brew_data[\"glycerol%\"], \n",
    "        1.13*brew_data[\"glycerol%\"] + 0.41, \n",
    "        'r', \n",
    "        label='y={:.2f}x+{:.2f} (binary classifier)'.format(1.13, 0.41))\n",
    "ax.set_xlabel(\"Glycerol (%)\")\n",
    "ax.set_ylabel(\"Acetaldehyde (%)\")\n",
    "plt.legend(fontsize=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "95acdaf7-acef-4e50-a858-e339636d4052",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define predictor (X) and response (y) variables\n",
    "\n",
    "X = brew_data[brew_data.columns]\n",
    "y = brew_data.index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "aaa17464-5cfb-4a2d-b685-93e5c39bf699",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearDiscriminantAnalysis()"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Fit the LDA model\n",
    "\n",
    "lda = LinearDiscriminantAnalysis()\n",
    "lda.fit(X.values, y.values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "d2e86490-4480-4a25-8f15-b6d73806dbe8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Variance ratio shows 100% discrimination on 1 axis\n",
    "\n",
    "lda.explained_variance_ratio_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "0ee320d7-9ef4-4b62-837c-4cad776fceab",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define method to evaluate model\n",
    "\n",
    "cv = RepeatedStratifiedKFold(n_splits=6, \n",
    "                             n_repeats=3,\n",
    "                             random_state=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "0a81db12-135a-412b-aeed-9f14898bd06b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9351851851851851\n"
     ]
    }
   ],
   "source": [
    "# Evaluate accuracy of model with 3 repeats of 10-fold cross-validation (as defined above)\n",
    "\n",
    "scores = cross_val_score(lda, X, y, \n",
    "                         scoring='accuracy', \n",
    "                         cv=cv, n_jobs=-1)\n",
    "\n",
    "print(np.mean(scores))    # model achieves a mean accuracy of 93.5%"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "a438fc27-5eaa-4568-b419-f4aa8fae7f5c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The new observation is predicted to be from Brew 1\n"
     ]
    }
   ],
   "source": [
    "# Define new observation and predict Brew\n",
    "\n",
    "new = [1.3, 1.8]\n",
    "\n",
    "print('The new observation is predicted to be from Brew {}'.format\\\n",
    "      (lda.predict([new])[0]))    # model correctly predicts Brew 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "90fc6071-54f2-452a-870a-f51c68ea3065",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Predicted: Brew 1\tPrior: Brew 1\n",
      "Predicted: Brew 1\tPrior: Brew 1\n",
      "Predicted: Brew 1\tPrior: Brew 1\n",
      "Predicted: Brew 1\tPrior: Brew 1\n",
      "Predicted: Brew 1\tPrior: Brew 1\n",
      "Predicted: Brew 1\tPrior: Brew 1\n",
      "Predicted: Brew 1\tPrior: Brew 1\n",
      "Predicted: Brew 1\tPrior: Brew 1\n",
      "Predicted: Brew 1\tPrior: Brew 1\n",
      "Predicted: Brew 2\tPrior: Brew 2\n",
      "Predicted: Brew 2\tPrior: Brew 2\n",
      "Predicted: Brew 2\tPrior: Brew 2\n",
      "Predicted: Brew 2\tPrior: Brew 2\n",
      "Predicted: Brew 2\tPrior: Brew 2\n",
      "Predicted: Brew 2\tPrior: Brew 2\n"
     ]
    }
   ],
   "source": [
    "# Loop over original dataset to predict brew:\n",
    "\n",
    "for index, row in brew_data.iterrows():\n",
    "    print('Predicted: Brew {}'.format(lda.predict([[row[0], row[1]]])[0]) + \\\n",
    "          '\\tPrior: Brew {}'.format(index))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38c21ba8-330e-49c9-b866-d703bde1436e",
   "metadata": {},
   "source": [
    "##### Now we will reclassify one of the samples (observation 8) from Brew 1 to Brew 2, and re-run the LDA to see if the model reveals and corrects the missclassification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "2741a6db-9ec5-4e89-a5f1-06e1e398a277",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   glycerol%  acetaldehyde%\n",
      "1       1.38           1.64\n",
      "1       1.40           1.70\n",
      "1       1.24           1.72\n",
      "1       1.36           1.74\n",
      "1       1.38           1.82\n",
      "1       1.48           1.82\n",
      "1       1.54           1.82\n",
      "2       1.38           1.90\n",
      "1       1.56           2.08\n",
      "2       1.14           1.78\n",
      "2       1.20           1.86\n",
      "2       1.18           1.96\n",
      "2       1.30           1.96\n",
      "2       1.26           2.00\n",
      "2       1.28           2.00\n"
     ]
    }
   ],
   "source": [
    "# Reclassify one sample in dataset (observation 8 Brew 1 --> Brew 2)\n",
    "\n",
    "as_list = brew_data.index.to_list()\n",
    "as_list[7] = 2\n",
    "\n",
    "brew_data_reclassified = brew_data\n",
    "brew_data_reclassified.index = as_list\n",
    "print(brew_data_reclassified)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "bafd11ec-25f6-4769-9c29-f8875d35035b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define predictor (X) and response (y) variables for reclassified dataset\n",
    "\n",
    "X_2 = brew_data_reclassified[brew_data_reclassified.columns]\n",
    "y_2 = brew_data_reclassified.index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "966ce0fc-eb5b-4c68-803a-c9a8f520afbe",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearDiscriminantAnalysis()"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Fit the LDA model for the reclassified dataset\n",
    "\n",
    "lda_2 = LinearDiscriminantAnalysis()\n",
    "lda_2.fit(X_2.values, y_2.values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "19fe5256-439e-4898-949a-c0338117c1c0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.])"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Variance ratio shows 100% discrimination on 1 axis for reclassified dataset\n",
    "\n",
    "lda_2.explained_variance_ratio_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "400e8d6f-ed2b-423f-9717-9a11b9eff65b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define method to evaluate model for reclassified dataset\n",
    "\n",
    "cv_2 = RepeatedStratifiedKFold(n_splits=6, \n",
    "                             n_repeats=3,\n",
    "                             random_state=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "1215d13f-28fb-4f40-9d1f-4515d449ac4c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7962962962962963\n"
     ]
    }
   ],
   "source": [
    "# Evaluate accuracy of model with 3 repeats of 10-fold cross-validation (as defined above)\n",
    "\n",
    "scores_2 = cross_val_score(lda_2, X_2, y_2, \n",
    "                         scoring='accuracy', \n",
    "                         cv=cv_2, n_jobs=-1)\n",
    "\n",
    "print(np.mean(scores_2))    # model achieves a mean accuracy of 79.6% for reclassified dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "41b2e120-46a5-43b7-a86a-b39529476396",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Predicted: Brew 1\tPrior: Brew 1\n",
      "Predicted: Brew 1\tPrior: Brew 1\n",
      "Predicted: Brew 1\tPrior: Brew 1\n",
      "Predicted: Brew 1\tPrior: Brew 1\n",
      "Predicted: Brew 1\tPrior: Brew 1\n",
      "Predicted: Brew 1\tPrior: Brew 1\n",
      "Predicted: Brew 1\tPrior: Brew 1\n",
      "Predicted: Brew 1\tPrior: Brew 2\n",
      "Predicted: Brew 1\tPrior: Brew 1\n",
      "Predicted: Brew 2\tPrior: Brew 2\n",
      "Predicted: Brew 2\tPrior: Brew 2\n",
      "Predicted: Brew 2\tPrior: Brew 2\n",
      "Predicted: Brew 2\tPrior: Brew 2\n",
      "Predicted: Brew 2\tPrior: Brew 2\n",
      "Predicted: Brew 2\tPrior: Brew 2\n"
     ]
    }
   ],
   "source": [
    "# Loop over reclassified dataset to predict brew:\n",
    "\n",
    "for index, row in brew_data_reclassified.iterrows():\n",
    "    print('Predicted: Brew {}'.format(lda_2.predict([[row[0], row[1]]])[0]) + \\\n",
    "          '\\tPrior: Brew {}'.format(index))\n",
    "\n",
    "# Model correctly reclassifies 8th observation from Brew 2 to Brew 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2516eea5-4219-49f4-8546-4845703b9e2f",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "Haswell, S.J. (1992). Practical Guide to Chemometrics. Marcel Dekker Inc., New York."
   ]
  }
 ],
 "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
}