E25_kpca_sklearn_model_try.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "f6bb6243",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    " \n",
    "from sklearn.preprocessing import StandardScaler, KernelCenterer\n",
    "from sklearn.decomposition import PCA, KernelPCA\n",
    "from sklearn.utils import extmath\n",
    "from scipy.stats import f\n",
    "from sklearn.metrics.pairwise import euclidean_distances\n",
    "from sklearn.preprocessing import StandardScaler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f8c5542c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The Datas\n",
    "data = pd.read_csv('plums.csv')\n",
    "X = data.values[:,1:]\n",
    "Xstd = StandardScaler().fit_transform(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "9d9fafc0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ker_pca(X, n_components=3, gamma = 0.01):\n",
    "    \n",
    "    # Calculate euclidean distances of each pair of points in the data set\n",
    "    dist = euclidean_distances(X, X, squared=True)\n",
    "    \n",
    "    # Calculate Gaussian kernel matrix\n",
    "    K = np.exp(-gamma * dist)\n",
    "    Kc = KernelCenterer().fit_transform(K)\n",
    "    \n",
    "    # Get eigenvalues and eigenvectors of the kernel matrix\n",
    "    eig_vals, eig_vecs = np.linalg.eigh(Kc)\n",
    "    \n",
    "    # flip eigenvectors' sign to enforce deterministic output\n",
    "    eig_vecs, _ = extmath.svd_flip(eig_vecs, np.empty_like(eig_vecs).T)\n",
    "    \n",
    "    # Concatenate the eigenvectors corresponding to the highest n_components eigenvalues\n",
    "    Xkpca = np.column_stack([eig_vecs[:,-i] for i in range(1,n_components+1)])\n",
    " \n",
    "    return Xkpca"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e3532be3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0sAAAGHCAYAAABh8cmQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAABqmUlEQVR4nO3de1yUdd4//tcMAwODIOJ4zMq9ya3ooCh5oL1FxBuHwya6hbRp0eH+bWmKrtpBMzW/mneew4y93bYx3V3QbSUWOcVBXZWKErcDdZdZba1uNiKKDAzMzPX7gxghhmGAmeu6Zub1fDz2scF8mHl/5PCZ93V9Pu+3QhAEAURERERERNSJUuoAiIiIiIiI5IjJEhERERERkR1MloiIiIiIiOxgskRERERERGQHkyUiIiIiIiI7mCwRERERERHZwWSJyEUyMjIwY8aMbj/+qa+//hoKhQLHjx/v9WtNmzYNjz32WJ/iJCIizyLm33yuL8Do0aPx//7f/5M6DJIJJktEAJqamrB69WqMGTMGQUFBCA8Px1133YWXX37Z6efYuXMnDh486PT466+/HufPn8ekSZMAAN999x0UCgWOHDnS2/CJiEgE33//PRYtWoTRo0cjICAAQ4YMwa9+9SucPn3ara/717/+Fdu2bXPra3gylUoFvV7f66977LHHMG3atC6fr66uxtKlS/sfGHkFldQBEMnBE088gcrKSuzcuRNjx47FlStXUFNTg3/+859OP8fAgQN79Zp+fn4YPnx4b0MVVUtLCwICAqQOg4hIct9++y1iYmIwdOhQvPrqq7jtttvw73//G9u3b8fkyZORl5cHnU7X5+cXBAFmsxn+/v5dHgsPD+9P6NRLQ4YMkToEkhHeWSICkJeXhxUrViA1NRU/+9nPMHbsWGRkZOD555/vNC43NxcTJkxAYGAgBg8ejMTERFy6dAlAz9vuvvnmG9x6661IT09HS0tLl214119/PQAgLi4OCoUCo0eP7tUcsrKycMsttyAwMBBjxozBhg0bYDabbY//6U9/wqRJkzBw4EBotVokJyfj888/tz3eHs8f//hHJCUlITg4GKtXr8batWtx00034a233sItt9yC4OBgTJs2DV988UWv4iMi8mQLFy5Ea2srKisrkZiYiBtuuAETJ07En//8Z0yfPh0ZGRloamoCANvfzY6OHz8OhUKBr7/+GgCg1+uhUqlQWVmJqKgoqNVqlJWV2X3tn26NmzZtGh599FE899xzGDp0KMLCwrBq1SpYrVa88MILGDZsGIYMGYJVq1Z1ep7Ro0dj1apVeOyxxxAaGgqtVouVK1fCarU6nHtP68vo0aOxevVqPPHEEwgLC8PQoUOxa9cumEwmLFq0CIMGDcJ1112HXbt2dXreq1evIjMzE9dddx00Gg2ioqLw17/+1fZ4+7p04MABpKSkQKPR4D/+4z863UUaPXo0LBYLHn74YSgUCigUCgDApUuXMG/ePNxwww0ICgrCzTffjK1bt0IQBNv36LXXXsPRo0dtX9f+vD/dhtfQ0IDf/OY3GDJkCNRqNaKjo1FaWtqrOMlzMVkiAjBixAgUFxejrq6u2zGvv/465s2bh9TUVJw6dQqVlZXQ6XSwWCw9Pv8//vEPTJkyBTNnzsSf//xnu3drTp06BQB48803cf78eVRXVzsd/9q1a7Flyxa8+OKL+PTTT7Fz50787ne/w7p162xjTCYTnnvuOZw6dQpvv/02/Pz8kJycjJaWlk7P9fTTT+OBBx7Axx9/jMcffxwAcP78ebz66qv44x//iJMnT6KhoQGPPPKI0/EREXmyS5cu4fDhw3jyyScRGhra5fFnn30W33//Pd5+++1ePa/VasXTTz+Nbdu24bPPPkN0dLTTX/uXv/wFra2tOH78OLZt24aNGzciOTkZV69exd///nds2bIFGzduRFFRUaevy8rKwsiRI1FdXY3t27dj586dyMrK6vZ1nFlf2p93zJgxeP/997F48WIsWrQIs2fPxs9+9jNUV1fjySefxOLFi1FbWwug7U7aL3/5S/zjH/9Abm4uPv74YzzxxBNIT09HeXl5p+d+5pln8OCDD+LDDz9Eeno6HnvsMdvFvurqavj5+WHHjh04f/48zp8/D6Btzbv99tuRl5eH2tparF69GmvWrLElMMuXL8evf/1rTJkyxfZ1c+fOtftv8Mgjj6CkpAT79+/H6dOncffddyMlJQWfffaZ03GSBxOISDh+/Lhwww03CEqlUrjjjjuE//7v/xYOHTokWK1W25jrr79eWLhwYbfP8dBDDwnx8fFdPi4vLxcGDhwobNq0qdP4r776SgAg/P3vfxcEQRC+/fZbAYBQWVnZY7yxsbHCo48+KgiCIDQ2NgpBQUFCUVFRpzF79+4VBg4c2O1zXLx4UQAgHD9+vFM8L7zwQqdxa9asEfz8/IQLFy7YPpeTkyMoFAqhqampx1iJiDzdu+++KwAQ/vrXv9p9vP3v6UsvvSQIQtvfzYiIiE5j/v73vwsAhK+++koQBEF4/fXXBQDCsWPHenz9jn/z2z8eO3ZspzGRkZHC7bff3ulzd955p7Bs2TLbxzfeeKPwi1/8otOYZ599Vhg1apTd13J2fbnxxhuFWbNm2T62WCxCSEiIkJKS0ulzYWFhQlZWliAIglBZWSmo1Wqhvr6+03M//PDDtudqX5e2bt1qe9xsNgsDBgwQsrOzbZ/z8/MTXn/9daEnixcvFmbMmGH7+NFHHxViY2O7jLvxxhuF9evXC4IgCF988YUAQDh8+HCnMVFRUcLDDz/cqzjJM/HMEhGAu+++G19++SXee+89VFVV4dixY7j33nuRmJiI/Px8/PDDD/j222+RkJDQq+f96KOPkJiYiG3btmHhwoVuif2TTz5BU1MTfvWrX9m2HwCAxWJBc3MzfvjhBwwZMgSnT5/GunXrcPr0aRgMBttWhG+++QZ333237esmTpzY5TVGjhzZaQ/3yJEjIQgCLly4gBtuuMEt8yIi8gV33XVXn75u7NixnT4ePnx4l3Oww4cPx4ULFzp9bsqUKZ0+vvvuu/Hiiy/iypUrXe6aObu+/DQepVKJIUOG4M477+z0uaFDh9riqa6uRktLC6677rpOr9nS0oIxY8Z0+ty4ceNs/+3n54ehQ4fi+++/7/qP0oHVasVLL72EnJwcfPfdd2hubkZraytuvPFGh1/3U+13wqZOndrp81OnTkVVVVW/4yT5Y7JE9COVSoWYmBjExMRg2bJl2L9/P+bPn49jx47h1ltv7dNz3nDDDRg9ejT279+PefPm9boIhDPa95ofPHgQP//5z7s8Hh4eDqPRiISEBPziF7/A66+/jmHDhgEAbrvtti7b8IKDg7s8x0+3DbYvmj3tcyci8gY33XQTFAoFPv74Y8yePbvL45988gkA4OabbwbQlhi0X5Bq19ra2uXr/Pz8EBgY2KeYfloIQqFQ2P1cf/5OO7O+9DUeq9WKgQMH2t1y/tM1x94a1NO8tm7dihdffBHbt29HVFQUQkJCsH37dhw+fNjh1/VHX+Ik+eOZJaJutCdIFy5cwNChQzFq1KhOBzqdMXDgQLz99ttQKpWYMWOGrRiEPe1/ZJ05A9XRbbfdhsDAQJw9exY33XRTl//5+fnh008/xQ8//IANGzZg2rRpuPXWW3Hp0qUuizkREXUVHh6OpKQk7Nq1C1euXOny+Isvvohhw4bhv/7rvwDAdgel49/z9nOpUnvnnXc6fXzy5Elcd911ds9iObO+9FV0dDTq6+vR3Nzc5Xl7u2MhICCgy9p57Ngx6HQ6PPLII4iKisJNN93UpTCRva/7qdtuu832fD99/ttvv71XcZJnYrJEBCA2NhbZ2dl4//338c0336C8vBwLFixAWFgY4uLiAABr1qzB7373O6xfvx6ffvopPvnkE+zatQsGg8Hhc4eGhqKkpAQajQbx8fG4ePGi3XFarRYDBgxAaWkp/v3vfztMrDoaMGAAVq5ciZUrV+KVV17B//3f/+GTTz5BTk4Onn76aQDAjTfeCLVajaysLHz55ZcoLy9HZmZmp20VRETUvVdeeQUqlQrTp09HcXExvv32W1RXV+PXv/41KioqoNfrERQUBKCtqqnRaMTzzz+PL7/8EgcPHsQrr7wi8QzanD59GmvXrsXnn3+OP/3pT9i5cyeWLVtmd6wz60tfTZ8+HTNmzMCcOXOQl5eHs2fP4oMPPkBWVhb27NnTq+f62c9+hsrKSpw7d862Jt988804cuQIKisr8fnnn+O5557Du+++2+XrPvvsM3zyyScwGAwwmUxdnjsiIgL33XcfFixYgJKSEnz22WfIzMzExx9/jBUrVvT9H4A8BpMlIgCJiYm2ktk333wzHn74YYwZMwYnTpyAVqsF0Na8Tq/X4y9/+QvGjRuHqVOnoqioCCpVz7tZBwwYgKKiIgwaNAhxcXFd9pADbds2XnnlFRw4cACjRo1CVFSU0/GvXr0a27Ztw549ezB27Fj84he/wPbt223lx7VaLfbv34+3334bt912G5YvX44tW7ZAqeSfACIiZ9x444344IMPMGnSJPzmN79BREQEEhMTYTKZUFVV1anH0s0334w9e/bgz3/+M26//Xb84Q9/wMaNGyWM/ppFixbhm2++QXR0NBYtWoQnn3wSmZmZ3Y7vaX3pK4VCgfz8fMyZMwdLly7FLbfcguTkZBw+fBgRERG9eq6tW7figw8+wOjRo21nqFavXo3Y2FjMmjULU6ZMwaVLl7B48eJOX/foo4/irrvuQkxMDIYMGYI///nPdp//97//PWbOnIl58+Zh7NixOHHiBAoKCnDLLbf0bfLkURQC9+EQEREReb3Ro0fjsccew3PPPSd1KEQeg5eViYiIiIiI7GCyREREREREZAe34REREREREdnBO0tERERERER2MFkiIiIiIiKyo+eax0RERITTp0/j9ddfh9VqRXx8PFJTUzs9XlpaipKSEiiVSgQGBuI3v/kNRo0aJU2wRETkEl5/ZuncuXOivI5Wq+2xOakn4DzkhfOQF86j70aOHCnq67ma1WpFZmYmnnvuOQwePBjPPvssMjMzOyVDRqMRGo0GAPD++++jpKQEq1atcur5xVqrxOQtvy+9wTn7Bl+cM+D98+5uneI2PCIioh6cOXMGw4cPx7Bhw6BSqRATE4Pq6upOY9oTJQBobm6GQqEQO0wiInIxbsMjIiLqQV1dHQYPHmz7ePDgwfjiiy+6jCsuLsbhw4dhNpvx/PPPd/t8ZWVlKCsrAwBs2rQJWq3W9UFLTKVSeeW8HOGcfYMvzhnw4XlLHQAREZG30Ol00Ol0OH78ON588008+eSTdsfNmDEDM2bMsH3sjVtbvH3Ljj2cs2/wxTkD3j9vbsMjIiLqo/DwcFy8eNH28cWLFxEeHt7teHvb9IiIyPMwWSIiIupBREQEzp8/jwsXLsBsNuPkyZOIjo7uNOb8+fO2/z516hRGjBghdphERORi3IZHRETUAz8/PzzyyCPYsGEDrFYr4uLicP311yM3NxcRERGIjo5GcXExPvroI/j5+WHAgAFYuHCh1GETEVE/MVkiIiJywvjx4zF+/PhOn5s7d67tvx9++GGxQyIiIjeTTbLUU7M/ADh58iQOHjwIhUKBG2+8EZmZmeIHSm5jsVhQWVKM6sP5ULW2wOwfgIkpsxA3UwelkjtGiYhIWhaLBYWFR3HwYA2MRn9oNK1ISxuPpKRYrlNEXkoWyZLVasVrr73WqdlfdHR0p2Z/58+fR15eHtavX48BAwbg8uXLEkZMrmYwGLA1cwGSzPVYEaKGQqGAYBJQpd+BVfv1WLZzt0+WqyQiInkwGAzIyMhCbW06TKYHACgACDh+vALZ2eug1y/iOkXkhWRxGcSZZn/l5eWYOXMmBgwYAAAYOHCgFKGSG1itVmzNXICVAY2ICQ20NXJUKBSICQ3EyoBGbM1cAKvVKnGkRETki6xWKzIyslBTswMmUzzaEiUAUMBkikdNzQ5kZGRxnSLyQrK4s+RMs79z584BAFavXg2r1Yr77rsP48aNEzNMcpPKkmIkmeuh0QTafVyj8kOisR5HSksxXafr02twix8REfVVYeFR1NamAwjuZkQwamvn4q233sZ//ueEPr0Gt/gRyZMskiVnWK1WnD9/HmvWrEFdXR3WrFmDLVu2IDi48x8uqbqie0tXYynm8Y/SIiwOUTscExOixsslh5E2b55Tz9lxHhcuXMDaR+ZjpukiVgwIsG3xe2fvTqzJ2Y91r+/H0KFD+z0Pd+DPlbxwHkS+6cCBUzCZtjgcYzLFY+/elX1KlrjFj0i+ZJEsOdPsLzw8HGPGjIFKpcLQoUMxYsQInD9/HjfddFOncVJ1RfeWrsZSzMPSeNW29a47CoUCFuNVp2Nrn4fVasWq+elYGdAITYeETKFQYEqIGmPNl7Fqfjo27MuR5ZU7/lzJC+fRd911RifyBEajP65tveuOAkZj799Wddzi1/nOVfsWv8nIyFiC/Pw1slyniLydLH7rnGn2N3HiRHzyyScAgCtXruD8+fMYNmyYFOGSi5n9AyAIgsMxgiDA7B/Q6+e2bfFT+dl9XKPyQ6K5bYsfERGRPRpNKwDH6xQgQKMx9/q5nd3iV1x8rNfPTUT9J4tkqWOzv6VLl2LKlCm2Zn/vv/8+AGDs2LEICQnB0qVLsW7dOsybNw8hISESR06uMDFlFqoaTA7HnGwwYVLK7F4/93sFb2GKE1v83i041OvnJiIi35CWNh5qdYXDMWp1OR56KKbXz922xW+6wzEmUzxycj7o9XMTUf/JYhse0HOzP4VCgYceeggPPfSQ2KGRm8XN1GHVfj3GmRvt3gEymi0oUoVhQ0JCr59b1dri1BY/VWtLr5+biIh8Q1JSLLKz16GmZjLs3wFqRGRkLmbNykJdXV2vntvZLX5NTf69el4icg1Z3Fki36ZUKrFs525sbAnGiSvNti15giDgxJVmbGwJxrKdu/u0V9udW/yIiMg3KJVK6PWLEBW1BGp1Ga5tyROgVpchKmoJ9PpFfVqnnN3iFxTU2uvnJqL+k82dJfJtWq0WG/bloLK0BJsL8mzlvSfdOxsbEhL6fKh1YsosVOl3ICbUfllyADh6sQGRv47va+hEROQDtFot8vPXoKjoKHJzl6OpyR9BQa1IT58Ana7vxRfS0sbj+PGKH/s32adQlCApaUxfQyeifmCyRLKhVCoRr0tEvC7RZc8ZN1OHVfscb/ErOX8JykMH8ctf3ctKQ0RE1C2lUonk5DgkJ8e57DmTkmLx6qvrcPp091v8BCEP+/aZkZaWwnWKSGT8jSOvplQqMXHWHDz7ybc4YbjSeYuf4Qo2fvYvLPv5SCRZLrMiHhERiU6pVOKBB6KgUPwGQAU6bvFr+/gpAC/g00/TWRGPSAJMlsjrfXasEtvuuAEtVgGbPz+HrZ+fw+bPz6FVEPDCbddDq/ZnRTwiIpJMcfGXEIS9AEwAngOw5sf/bwGQBWAoK+IRSYTb8MjrqVpb4KdUIm7oQMQNHWh3DCviERGRVNoq4vkBSPzxf/awIh6RFHhnibweK+IREZGcsSIekXwxWSKv586mt0RERP3lbNPb9PQJIkVERO2YLJHXi5upQ6EqDEazxe7j7U1vp/Wh6S0REVF/JSXFIjIyB0BjNyPamt7qdFPFDIuIwGSJfIA7m94SERH1lzub3hJR/7DAA/kEdzW9JSIicgV3Nb0lov5hskQ+wx1Nb4mIiFzFHU1viah/eJmCiIiIiIjIDiZLREREREREdjBZIiIiIiIisoPJEhERERERkR1MloiIiIiIiOxgskRERERERGQHS4eTQxaLBZUlxag+nG/rTTQxZRbiZurY84GIiCRnsVhQWHgUBw/WwGj0h0bTirS08UhKiuU6RUT9xmSJumUwGLA1cwGSzPVYEaKGQqGAYBJQpd+BVfv1WLZzN7RardRhEhGRjzIYDMjIyEJtbTpMpgcAKAAIOH68AtnZ66DXL+I6RUT9wksuZJfVasXWzAVYGdCImNBAKBQKAIBCoUBMaCBWBjRia+YCWK1WiSMlIiJfZLVakZGRhZqaHTCZ4tGWKAGAAiZTPGpqdiAjI4vrFBH1C5MlsquypBhJ5npoVH52H9eo/JBorseR0lKRIyMiIgIKC4+itjYdQHA3I4JRWzsXxcXHxAyLiLwMkyWy672CtzAlRO1wTEyIGu8WHBIpIiIiomsOHDgFk2m6wzEmUzxycj4QKSIi8kZMlsguVWuLbetddxQKBVStLSJFREREdI3R6I9rW++6o0BTk78Y4RCRl2KyRHaZ/QMgCILDMYIgwOwfIFJERERE12g0rQAcr1OAgKCgVjHCISIvxWSJ7JqYMgtVDSaHY042mDApZbZIEREREV2TljYeanWFwzFqdTnS0yeIFBEReSMmS2RX3EwdClVhMJotdh83mi0oUoVhWkKCyJEREREBSUmxiIzMAdDYzYhGREbmQqebKmZYRORlmCyRXUqlEst27sbGlmCcuNJs25InCAJOXGnGxpZgLNu5mw3/iIhIEkqlEnr9IkRFLYFaXYZrW/IEqNVliIpaAr1+EdcpIuoXNqWlbmm1WmzYl4PK0hJsLsiDqrUFZv8ATLp3NjYkJHABIiIiSWm1WuTnr0FR0VHk5i5HU5M/goJakZ4+ATrdGq5TRNRvTJbIIaVSiXhdIuJ1iVKHQkRE1IVSqURychySk+OkDoWIvBAvuRAREREREdnBZImIiIiIiMgOJktERERERER2MFkiIiIiIiKyg8kSERERERGRHUyWiIiIiIiI7GDpcCIiIiecPn0ar7/+OqxWK+Lj45Gamtrp8YKCApSXl8PPzw+hoaF44oknMGTIEGmCJSIil+CdJSIioh5YrVa89tprWLlyJbZv344TJ07gu+++6zRm9OjR2LRpE7Zs2YLJkydj//79EkVLRESuwmSJiIioB2fOnMHw4cMxbNgwqFQqxMTEoLq6utOY22+/HWq1GgAwZswY1NXVSREqERG5ELfhERER9aCurg6DBw+2fTx48GB88cUX3Y6vqKjAuHHjun28rKwMZWVlAIBNmzZBq9W6LFa5UKlUXjkvRzhn3+CLcwZ8eN5SB0BEPbNYLKgsKUb14XyoWltg9g/AxJRZiJupg1LJG8REcnLs2DGcPXsWa9eu7XbMjBkzMGPGDNvHBoNBhMjEpdVqvXJejvjynC0WCwoLj+LgwRoYjf7QaFqRljYeSUmxXrdO+eL3GfD+eY8cOdLu55ksEcmcwWDA1swFSDLXY0WIGgqFAoJJQJV+B1bt12PZzt0+eaWHSEzh4eG4ePGi7eOLFy8iPDy8y7gPP/wQhw4dwtq1a+Hv7y9miESSMRgMyMjIQm1tOkymBwAoAAg4frwC2dnroNcv4jpFHks2qf7p06eRmZmJRYsWIS8vr9tx77zzDtLS0vDll1+KFxyRRKxWK7ZmLsDKgEbEhAZCoVAAABQKBWJCA7EyoBFbMxfAarVKHCmRd4uIiMD58+dx4cIFmM1mnDx5EtHR0Z3GfPXVV9izZw+eeuopDBw4UKJIicRltVqRkZGFmpodMJni0ZYoAYACJlM8amp2ICMji+sUeSxZ3FlqrzL03HPPYfDgwXj22WcRHR2NUaNGdRrX1NSEoqIijBkzRqJIqTvcJuYeJfn5SDLXQ6MJtPu4RuWHRGM9jpSWYrpOJ3J0RL7Dz88PjzzyCDZs2ACr1Yq4uDhcf/31yM3NRUREBKKjo7F//340Nzdj27ZtANq2rDz99NMSR07tfGmbmJjy8kpRW5sOILibEcGorZ2L4uJjSEqaJmJkRK4hi2SpY5UhALYqQz9NlnJzczFr1izk5+dLESZ1g9vE3Ofvb+ZicYja4ZiYEDU2FxxiskTkZuPHj8f48eM7fW7u3Lm2/169erXYIZGTuE3MffburYLJtNHhGJMpHjk5y5kskUeSxaUUe1WGflpy9ezZszAYDF0WKpIWt4m5l7LVZPs37Y5CoYCqtUWkiIiIPAu3iblXU5MK1/5Nu6NAUxPP8JFnksWdpZ5YrVa88cYbWLBgQY9jpSrH6i3lFHs7j6K8PCRbLkOjsn/3Q6PyQ5LxMk69UwXdPbNcFWaPvOX7IQQEQmgWHCZMgiDATzNA1vP1lu8H50HkeQoLj3KbmBsFBZkBCHCcMAkICmoVKSIi15JFstRTlaHm5mZ8++23WLduHQCgvr4eL730Ep566ilERER0ei6pyrF6SznF3s6j7M/7sGJAgMMxU0ICsPlPbyA65u7+huc0b/l+/GJOGqqyXkBMqP0zSwBwssGEcfcmOz1fKc6Xecv3g/Pou+5KshK524EDp2AybXE4htvE+u6hh6agsrLix7t29qnV5UhPn+D0c/J8GcmJLJKljlWGwsPDcfLkSSxevNj2uEajwWuvvWb7eO3atZg/f36XRInEp2pt4TYxN5p5zz34zasvY5y5ERqVX5fHjWYLilRh2JCQ4NTz8XwZEfkao9Ef3CbmPqmpCYiMXISamsmwf/euEZGRudDp1jj1fDxfRnIji/S8Y5WhpUuXYsqUKbYqQ++//77U4ZEDZv8ACILgcIwgCDD7O777RPYplUos27kbG1uCceJKs+3fWhAEnLjSjI0twVi2c7dTV9p4voyIfJFG04q2bWKOcJtYXymVSuj1ixAVtQRqdRmu/VsLUKvLEBW1BHr9IqfXKZ4vI7mRxZ0loOcqQx056opO4pqYMgtV+h09bhObdO9sEaPyLlqtFhv25aCytASbC/JsW+cm3TsbGxISnN6SUFlSzDLkRORz0tLG4/hx124To860Wi3y89egqOgocnOXo6nJH0FBrUhPnwCdbo3T6xTPl5EcySZZIs8UN1OHVfv1LtsmRvYplUrE6xIRr0vs83O8V/AWVrAMORH5mKSkWGRnr3PZNjGyT6lUIjk5DsnJcX1+Dp4vIzmSxTY88lyu3CZG7sXzZUTki1y5TYzci+fLSI54Z4n6zVXbxMi9zP4BEEw9lyHn+TIi8jau2iZG7nXtfBnLkJN8MFkil3DFNjFyL54vIyJf5optYuRePF9GcsRLKUQ+Im6mDoWqMBjNFruPt58vm8bzZUREJIGkpFhERuYAaOxmRPv5sqlihkU+jskSkY/g+TIiIpIzni8jOeI2PCIfwvNlREQkZzxfRnLDZInIx/B8GRERyRnPl5GcMD0nIiIiIiKyg3eWiAgAYLFYUFlSjOrD+bbteRNTZiFupo7bHoiISHIWiwWFhUdx8GANjEZ/aDStSEsbj6SkWK5T5DZMlogIBoMBWzMXIMlcjxUhaigUCggmAVX6HVi1X49lO3dDq9VKHSYREfkog8GAjIws1Namw2R6AG29mAQcP16B7Ox10OsXcZ0it2AaTuTjrFYrtmYuwMqARsSEBtqa1ioUCsSEBmJlQCO2Zi6A1WqVOFIiIvJFVqsVGRlZqKnZ8WMPpvamtQqYTPGoqdmBjIwsrlPkFkyWiHxcZUkxksz10Kj87D6uUfkh0VyPI6WlIkdGREQEFBYeRW1tOoDgbkYEo7Z2LoqLj4kZFvkIJktEPu69grcwJUTtcExMiBrvFhwSKSIiIqJrDhw4BZNpusMxJlM8cnI+ECki8iU8s+QleDif+krV2mLbetcdhUIBVWuLSBERkTfi4XzqK6PRH9e23nVHgaYmfzHCIR/DZMkL8HA+9YfZPwCCSXCYMAmCALN/gIhREZE34eF86g+NphWAAMcJk4CgoFaRIiJfwks5Ho6H86m/JqbMQlWDyeGYkw0mTEqZLVJERORNeDif+istbTzU6gqHY9TqcqSnTxApIvIlTJY8HA/nU3/FzdShUBUGo9li93Gj2YIiVRimJSSIHBkReQMezqf+SkqKRWRkDoDGbkY0IjIyFzrdVDHDIh/BZMnD8XA+9ZdSqcSynbuxsSUYJ640QxAEAG1b705cacbGlmAs27mbZwqIqE94OJ/6S6lUQq9fhKioJVCry9C2JQ8ABKjVZYiKWgK9fhHXKXILnlnycDycT66g1WqxYV8OKoqK8LT+9zAZvofJYoV6yHDMevBBhIeHSx0iEXkoHs4nV9BqtcjPX4PDhyuxa1cGvvuuFUAzrr/eH//f/5fCdYrchim4hzP7B9juBHSHh/PJGXV1dSj/016kKhuxI2Iwsm8eih2DLAh842Wsmp8Og8EgdYhE5IGuHc53hIfzqWd1dXX43e+O4YsvHkR9/R9RX/9XfPRRDpYsGYl77lnHdYrcgsmSh+PhfHIFFgohInfh4XxyBRYKIakwWfJwPJxPrsBCIUTkLjycT67AQiEkFSZLHo6H88kVWCiEiNyFh/PJFVgohKTCAg9eoP1wfmVpCTYX5EHV2gKzfwAm3TsbGxISuABRj1gohIjcqf1wflHRUeTmLkdTkz+CglqRnj4BOt0arlPUIxYKIakwWfISSqUS8bpExOsSpQ6FPJDZPwCCSXCYMLFQCBH1h1KpRHJyHJKT46QOhTzQtUIhjhImFgoh1+OlHCJioRAiIpI1FgohqTBZIiIWCiEiIlljoRCSCpMlImKhECIikjUWCiGp8MwSEQEABg0ahLj75+PA67/Hvs+/h9pPicAhw5H48BPYMFPHBYiIiCSl1Wpx6NBz2LjxVbz55itobQ2Gv38j7r33FqxcuRoqFd/Wkuvxp4qIYDAYsDVzAZLM9dgxSA1F+FAIgoCqhqso3K/H2LsmQqvVSh0mERH5MIPBgIyMLNTWpsNkehFtxR4E6PUVeO+99dDrF3GtIpfjpWIiH2e1WrE1cwFWBjQiJjTQVhFPoVAgJjQQKwMasTVzAbuiExGRZKxWKzIyslBTswMmUzyuVcVTwGSKR03NDmRkZHGtIpdjskTk4ypLipFkrodG5Wf3cY3KD4nmehwpLRU5MiIiojaFhUdRW5sOILibEcGorZ2L4uJjYoZFPoDJEpGPe6/gLUwJUTscExOixrsFh0SKiIiIqLMDB07BZJrucIzJFI+cnA9Eioh8BZMlIh+nam1x2IwWaNuSp2ptESkiIiKizoxGfzhuSAsACjQ1+YsRDvkQJktEPs7sH2ArFd4dQRBg9g8QKSIiIqLONJpWXCsX3h0BQUGtYoRDPoTJEpGPm5gyC1UNJodjTjaYMClltkgRERERdZaWNh5qdYXDMWp1OdLTJ4gUEfkKJktEPi5upg6FqjAYzRa7jxvNFhSpwjAtIUHkyIiIiNokJcUiMjIHQGM3IxoRGZkLnW6qmGGRD2CyROTjlEollu3cjY0twThxpdm2JU8QBJy40oyNLcFYtnM3m9ISEZFklEol9PpFiIpaArW6DNe25AlQq8sQFbUEev0irlXkcmxKS0TQarXYsC8HlaUl2FyQB1VrC8z+AZh072xsSEjg4kME4PTp03j99ddhtVoRHx+P1NTUTo/X1tZi7969+Oabb7BkyRJMnjxZmkCJvJRWq0V+/hoUFR1Fbu5yNDX5IyioFenpE6DTreFaRW7BZImIALRdtYvXJSJelyh1KESyY7Va8dprr+G5557D4MGD8eyzzyI6OhqjRo2yjdFqtViwYAH+9re/SRgpkXdTKpVITo5DcnKc1KGQj2CyRERE1IMzZ85g+PDhGDZsGAAgJiYG1dXVnZKloUOHAkCPpfiJiMhzyCZZ6ml7Q0FBAcrLy+Hn54fQ0FA88cQTGDJkiDTBSsBisaCypBjVh/NtW6QmpsxC3EwdbzsTEblZXV0dBg8ebPt48ODB+OKLL/r8fGVlZSgrKwMAbNq0CVqttt8xSs1isSAvrxRvvPEOmppU0GgsePDByUhN9Z2tvCqVyiu+l73BOfsOn5231AEAzm1vGD16NDZt2gS1Wo3S0lLs378fS5culTBq8RgMBmzNXIAkcz1WhKihUCggmARU6Xdg1X49lu3c7ZM/vOQdeCGAfNGMGTMwY8YM28cGg0HCaPrPYDAgIyMLtbXpMJk2oq15qICKigpERi6CXr/IJ9YprVbr8d/L3vKFOVssFhQWHsXBgzUwGv0xcKACs2ffjqSkWJ9ap7z9ez1y5Ei7n5fFd7jj9gaVSmXb3tDR7bffDrVaDQAYM2YM6urqpAhVdFarFVszF2BlQCNiQgNt2zsUCgViQgOxMqARWzMXwGq1ShwpUe8ZDAY89+D9UO/diRWmc/itcBErTOcQoN+BVfPTvfqPMnmW8PBwXLx40fbxxYsXER4eLmFE8mG1WpGRkYWamh0wmeLRligBgAImUzxqanYgIyOL6xR5JIPBgFmzXkBm5giUl29BVdUmFBdvxOLFw3HPPeu4TvkAWSRL9rY3OEqGKioqMG7cOBEik15lSTGSzPXQqPzsPq5R+SHRXI8jpaUiR0bUP7wQQJ4kIiIC58+fx4ULF2A2m3Hy5ElER0dLHZYsFBYeRW1tOoDgbkYEo7Z2LoqLj4kZFlG/8UIAATLZhtcbx44dw9mzZ7F27Vq7j0u1D9xd+zj/UVqExSFqh2NiQtR4ueQw0ubN6/frect+VM7DvSwWC0r+lo/jbx6AstUEq78a//mruZh5zz12tyTYm0dRXh6SLZehUdn/+dao/JBkvIxT71RBd88st8yjt+T6/egtb5mHmPz8/PDII49gw4YNsFqtiIuLw/XXX4/c3FxEREQgOjoaZ86cwZYtW9DY2IgPPvgABw4cwLZt26QO3e0OHDgFk2mLwzEmUzxycpYjKWmaOEGRz/vp1jmNphVpaeN7tXWuNxcC+LPtvWSRLDm7veHDDz/EoUOHsHbtWvj7+9t9Lqn2gbtrH6el8WqPlZUUCgUsxqsueX1v2Y/KeTinL+eFOp6hW9x+hq5ZQFXWC/jNqy/bPUNnbx5lf96HFQMCHMY3JSQAm//0BqJj7u7fRF2EP1d9191ecE8yfvx4jB8/vtPn5s6da/vvm266CdnZ2WKHJTmj0R/Xrrh3R4GmJvvrNpEjfUl6Op+hewDtZ+iOH69AdvY6p8/Q8UIAATLZhufM9oavvvoKe/bswVNPPYWBAwdKFKn4zP4BEATB4RhBEGD2d/ymk+in+nJeyJVb51StLU5dCFC1tvRtgkQkCo2mFYDjdQoQEBTUKkY45EXsnRcqL9/i8LyQK7fO8UIAATJJljpub1i6dCmmTJli297w/vvvAwD279+P5uZmbNu2DStWrMD//M//SBy1OCamzEJVg8nhmJMNJkxKmS1SROQN+pr0uPIMHS8EEHmHtLTxUKsrHI5Rq8uRnj5BpIjIG/Q16XHlGTpeCCBAJskS0La9YefOncjKysKcOXMAtG1vaL/DtHr1auzZswebN2/G5s2b8fTTT0sZrmjiZupQqAqD0Wyx+7jRbEGRKgzTEhJEjow8WV+TnvcK3sIUJ87QvVtwqMcYeCGAyDskJcUiMjIHQGM3IxoRGZkLnW6qmGGRh+tr0tO2dW66w+du2zr3QY8x8EIAATJKlsg+pVKJZTt3Y2NLME5cabZdiRcEASeuNGNjSzCW7dztU3X+qf/6mvS4cuscLwQQeQelUgm9fhGiopZArS7DtSvxAtTqMkRFLYFev4jrFPVKX5MeV26d44UAAmRS4IEc02q12LAvB5WlJdhckGc7iD/p3tnYkOA7ndHJdfqa9Jj9AyCYBIdf6+zWOduFgMwFSDTWI6a9WIQg4GSDCUWqMF4IIPIQWq0W+flrUFR0FLm5y9HU1Na4c86cO6DTreHvMfVaX5Oea1vnHH2tc1vn2i8EZGQsQW3t3A7bAQWo1eWIjMzlhQAfwGTJQyiVSsTrEhGvS5Q6FPICfU16JqbMQpV+B2JCA7v9upMNJky617mtc1JcCOhLBUAi6plSqURychySk+MAeE/1SJJGX5OetLTxOH684sfExr7ebJ2T4kKAK8qek+sohJ5OWHu4c+fOifI63rIocB7y4q55lBcVIqCHpOeooQFvCRoMv26ULZkAgFXz07EyoNHueSej2YKNLcHYsC+n0x90uXw/OpY9n9LhTlZVgwmFP97JclROVi7z6C+WDpcfsdYqMXnL70tvcM6uU1BQicWLhztMepTKYtx00z5cf/1wWzIBAPfcsw41NTtg/7xTI6KiliA/v++Jjju/z53Lnk/HtTtZFYiMzHG67Lk7ePvPd3frFNNTIh/kzHmh0n/XYYtW2amceF1dnceeoXNl2XMiInIvZ84LWa35+PzzfZ3KidfV1XnsGTpXlj0n15HfTwoRuZ2jwiHHDVew8bN/YemYkVAqFF2SifDwcGzYl4PWh5dic+BIbFMOxubAkTA//Fts2Jcj2RWvnriy7DkREbmXo8IhQDmApwCsRdtb2c7JRHh4OPLz1yAr63vExy9HTMwziI9fjl27LiA/f41s1ylXlj0n1+GZJSIf9dPzQg3/Pg/zue8wc2goXrjteih/cp5Jo/JDorEtmZiu03ncGbr3Ct7CCicqAG4uOITpOp1IURERUXd+el7ou++u4vPPAUGYDSALXa/5X0smkpKmdTpD5wnaKgBucTimrQLgciQlTRMnKOKdJSJf1l445Nldv0PIsBHYGDkK04YM7JIotXO2h5IcubLsORERiaO9cMgbbyzDqFEDIAjZAHTo7i2ssz2U5MiVZc/JdZgsEREA708mzP4B6KmejbNlz4mISHzenkxcqwDoiHNlz8l1uA2PyIc4Kpvtyh5KcuTqsudEROR6jspmu7KHkhy5uuw5uQbvLBH5CIPBgOcevB/qvTuxwnQOvxUudqp0d2vsdFQ1mBw+x8kGEyaleGYy4UwFwCJVGKYlJIgcGRERAW3r1KxZLyAzcwTKy7egqmpTp0p3Ot1NUKsrHD6HJycTzlQAjIzMhU43VcywfB6TJSIf4EzZ7Hfy/oJCP+9NJhxVAJR72XMiIm/nTNnsP/7xA9x6q/cmE44qAMq97Lk34zY8Ih9gK5utsb8FTaPyQ5LxMi7+6jFsfOsvSDTWI6ZD09aTDSYU/di0Vaw/0o62DPankWDHCoDtzzvp3tnYkJDABYiISCLOlM3+9NN0bNz4BRSKJaitndshqRKgVpcjMjJX1GTC0ZbB/qxTHSsANjX5IyioFenpE6DT9b2RLvUdkyUiH+B02exj5bJIJgwGA7ZmLkCSuR4r2pM2k4Aq/Q6s2q/Hsp27+9wno70CoCeVPSci8nbOls0uLCySRTJhMBiQkZGF2tp0mEwPoD1pO368AtnZ66DXL+rXOuVpZc+9GZMlIh/Qm0p3UicTHbcMdrwT1r5lcJy5ERszF2DDvhxeYSMi8hK9qXQndTLRcctg5zth7VsGJyMjYwny83knyBvwO0jkAzypbLZty6DKz+7jGpUfEs1tzXGJiMg7eFLZ7Ly80h63DLY3xyXPx2SJyAdMTJnlMZXu3it4C1Oc2DLoqc1xiYioq7S08R5T6W7v3iqYTNMdjvHk5rjUGZMlIh/gSWWzvb05LhERdeVJZbObmlTw5ua41BmTJSIf4Ellsz1pyyAREbmGJ5XNDgoyw1O2DFL/scADkY/wlLLZE1NmoUq/AzGh9sucAz9uGbxX+i2DRETkOp5SNvuhh6agsrLix9Ll9sllyyD1H5MlIh8idaU7Z8TN1GHVfj3GmRvtFnlo3zK4QQZbBomIyLWkrnTnjNTUBERGLkJNzWTYL/LQvmVwjdihkRswWSIiWbFtGcxcIIvmuK7kjka7REQkrvYtgxkZ8miO60qOGu36KoXQ0+EAD3fu3DlRXker1cJgMIjyWu7EeciLL8/DarWisrQE73XcMpgyG9Mk3DLYn+9Hx0a7UzokgFUNJhT+mAD2tYFhb0nxczVy5EhRX8/TiLVWiclb/n71BufsG9rnbLVaf9wyeOonWwanemSi1LnR7nRcSwArEBmZ4/V9o7pbp3hniYhkyRO2DDqLjXaJiLyPJ2wZdJYzjXbnzFmOv/51lc+tU741WyIiCbDRLhERyVlh4dEeG+1+9NG9Ptlol8kSEZGbsdEuERHJ2YEDp3pstNvcPN0nG+0yWSIicjM22iUiIjkzGv3BRrv2MVkiInIzNtolIiI502hawUa79jFZIiJys4kps1DVYHI45mSDCZNS2GiXiIjEl5Y2Hmp1hcMxgYEVPtlol8kSEZGbxc3UoVAVBqPZYvfx9ka709hol4iIJJCUFIvIyBwAjd2MaMQdd/wFOt1UMcOSBSZLRERuZmu02xKME1eabVvyBEHAiSvN2NgS7LGNdomIyPO1N9qNiloCtboM17bkCVCryxAVtcQny4YD7LNERCQKrVaLDftyUFlags0dG+3eOxsbJGy0S0REBLStU/n5a35stLv8J41212Do0KE+14AY6CFZslqt+Nvf/ob/+7//w6hRo3DPPfdgwIABtsdffPFFPPvss24PkojIG3hTo1254DpFROQ63tRo11UcXsr805/+hKqqKkRGRuJf//oXVqxYge+++872+Geffeb2AImIiLrDdYqIiNzJYbJ04sQJPP3000hJScGKFStw33334YUXXsDZs2cBoMdSuERERO7EdYqIiNzJ4TY8o9GI0NBQ28fTp0/HgAED8OKLL2LZsmU9NlkkIiJyJ65TRETkTg6TpREjRuCLL77ALbfcYvvcxIkTERAQgM2bN6Olhd3miYhIOlyniIjInRxuw0tMTMS3337b5fPjxo3D0qVLOy1OREREYuM6RURE7uTwzlJsbGy3j91+++24/fbbXR4QERGRs7hOERGRO7GxBxERERERkR1MloiIiIiIiOxgskRERERERGSHwzNLYjp9+jRef/11WK1WxMfHIzU1tdPjra2t2LVrF86ePYuQkBAsWbIEQ4cOlSZYIiLyOVyniIh8j8M7S++//z7+93//1+5j//u//4uamhqXBGG1WvHaa69h5cqV2L59O06cONGpAzsAVFRUIDg4GFlZWUhOTsYf//hHl7w2ERF5Lq5TRETkTg6TpcOHD+M///M/7T42depU5OfnuySIM2fOYPjw4Rg2bBhUKhViYmJQXV3dacz777+PadOmAQAmT56Mjz/+mJ3ZiYh8HNcpIiJyJ4fJ0nfffYdbb73V7mM333xzl6tqfVVXV4fBgwfbPh48eDDq6uq6HePn5weNRoOGhgaXvD4REXkmrlNERORODs8stbS0oKmpCUFBQV0ea25ulmVn9LKyMpSVlQEANm3aBK1WK8rrqlQq0V7LnTgPeeE85IXzkB9PXKcA6dYqMXnTz5mzOGff4ItzBnx43o4eHD16NN555x3ExcV1eezdd9/F6NGjXRJEeHg4Ll68aPv44sWLCA8Ptztm8ODBsFgsMBqNCAkJ6fJcM2bMwIwZM2wfGwwGl8TYE61WK9pruRPnIS+ch7xwHn03cuRItzyvJ65TgHRrlZi85felNzhn3+CLcwa8f97drVMOt+HNnj0ber0eBQUF+OGHH2A2m/HDDz+goKAAe/fuxZw5c1wSXEREBM6fP48LFy7AbDbj5MmTiI6O7jRmwoQJOHLkCADgnXfewW233QaFQuGS1yciIs/EdYqIiNxJIfRw+vSdd97BG2+80emKmlarxfz58zF58mSXBXLq1Cns3bsXVqsVcXFxmDNnDnJzcxEREYHo6Gi0tLRg165d+OqrrzBgwAAsWbIEw4YN6/F5z50757IYHfGWbJvzkBfOQ148YR4WiwWVJcWoPpwPVWsLzP4BmJgyC3EzdVAq266PedOdJcDz1ylAvLVKTJ7w++JqnLNv8MU5A66bt8ViQWHhURw8WAOj0R8aTSvS0sYjKSnWtk5Jobt1qsdkqd25c+fQ0NCAkJAQty56rsZkqXc4D3nhPORF7vMwGAzYmrkASeZ6TAlRQ6FQQBAEVDWYUKgKw7Kdu6HVar0uWWrnqesUwGTJW3DOvsEX5wy4Zt4GgwEZGVmorU2HyTQdgAKAALW6ApGROdDrF0l2Lqq7dcPhmSWTyYQ333wT3377LX72s59h9uzZ8Pf3d0uARETUd1arFVszF2BlQCM0mkDb5xUKBWJCAzHO3IiNmQuwYV+OhFG6HtcpIiLPYLVakZGRhZqaHQCCOzyigMkUj5qaycjIWIL8/DWS3mH6KYeRvPbaa/jggw9w3XXX4d1338W+ffvEiouIiHqhsqQYSeZ6aFR+dh/XqPyQaK7HkdJSkSNzL65TRESeobDwKGpr09E5UeooGLW1c1FcfEzMsHrkMFk6ffo0nnvuOcybNw/PPvssPvjgA7HiIiKiXniv4C1MCVE7HBMTosa7BYdEikgcXKeIiDzDgQOnftx61z2TKR45OfL6O+4wWTKZTBg0aBCAtn2KRqNRlKCIiKh3VK0tPVZeUygUULXKs+9QX3GdIiLyDEajP9rOKDmiQFOTvLZSOzyzZLFY8PHHH9s+tlqtnT4GgNtvv909kRERkdPM/gEQTILDhEkQBJj9A0SMyv24ThEReQaNphWAAMcJk4CgoFaRInKOw2Rp4MCBePXVV20fDxgwoNPHCoUCu3btcl90RETklIkps1Cl34GY0MBux5xsMGHSvbNFjMr9uE4REXmGtLTxOH68AiZTfLdj1OpypKdPEDGqnjlMll555RWx4iAion6Im6nDqv16jDM32i3yYDRbUKQKw4aEhE6fd6Yvk5xxnSIi8gxJSbHIzl6HmprJsF/koRGRkbnQ6dZ0+qzUfZnkvxISEVGPlEollu3cjY0twThxpRntLfQEQcCJK83Y2BKMZTt3d1pYDAYDnnvwfqj37sQK0zn8VriIFaZzCNDvwKr56T7ZR4SIiNxDqVRCr1+EqKglUKvL0LYlD2jrs1SGqKgl0OsXdVmnZs16AZmZI1BevgVVVZtQXr4FixcPxz33rBNlnXJ4Z4mIiDyHVqvFhn05qCwtweaCPNudokn3zsaGhIROC1Bv+jJ5wh0mIiKSP61Wi/z8NSgqOorc3OVoavJHUFAr0tMnQKdb02WdkkNfJiZLREReRKlUIl6XiHhdosNxtr5MGvtnnDQqPyQa2/oyTdfp3BEqERH5IKVSieTkOCQnxzkc15u+TElJ01wdpg0vFxIR+SBf7ctERESeQS59mZgsERH5IF/ty0RERJ5BLn2ZmCwREfkgs3+ArQhEd7yxLxMREXmGa32ZHHF/XyYmS0REPmhiyixUNZgcjjnZYMKkFO/qy0RERJ4hLW081OoKh2PE6MvEZImIyAfFzdShUBUGo9li9/H2vkzTftKXiYiISAxJSbGIjMwB0NjNiPa+TFPdGger4UnI05tBEpHnsvVlylyARGM9YkLUUCgUEAQBJxtMKFKFdenLRL5H6maQROS72vsyZWQsQW3tXJhM8Wg7wyRArS5HZGRul75M7sBkSSIGgwFbMxcgyVyPFe1vUkwCqvQ7sGq/Hst27oZWq5U6TCLyYr3py0S+x2AwICMjC7W16TCZHkD7m5TjxyuQnb0Oev0irlNE5Fa96cvkLkyWJMBmkEQkF872ZSLfIpdmkEREzvZlctvrS/KqPs7WDFLlZ/dxjcoPiea2ZpBERERi600zSCIib8ZkSQJsBklERHIml2aQRERSY7IkATaDJCIiOZNLM0giIqkxWZIAm0ESEZGcyaUZJBGR1JgsSYDNIImISM7k0gySiEhqTJYkwGaQREQkZ3JpBklEJDWWDpeALzeDZCNeIiL5k0szSCmwES8RdcRkSSK+2AySjXiJiDyHHJpBio2NeInop5gsSciXmkGyES8RkeeRuhmkmNiIl4js4W87iYKNeImISM7YiJeI7GGyRKJgI14iIpIzNuIlInu4DY9EwUa8RM5hERQiabARL5FzfK0ICpMlEoXZPwCCSXCYMLERL/k6FkEhks61RryOEiY24iXf5otFUJgskSgmpsxClX4HYkIDux1zssGESfeyES/5JrkVQeEdLvI1aWnjcfx4xY9l0u1jI17yZXIsgiLGXS6ueCQKNuIlckxORVAMBgOee/B+qPfuxArTOfxWuIgVpnMI0O/AqvnpMBgMbo+BSGxsxEvkWF5eqayKoBgMBsya9QIyM0egvHwLqqo2obx8CxYvHo577lnnsrWKyRKJwtaItyUYJ640QxAEAG1b705cacbGlmCvbcRL5Ay5FEHpeIcrJjTQtnW2/Q7XyoBGbM1cAKvV6tY4iMTW3og3KmoJ1OoytG3JA9oa8ZYhKmqJ1zbiJXLG3r1VsimC0vEu17Wm2cC1u1w7kJGR5ZK1itvwSDS+2IiXyFlyKYJiu8Olsb9lVqPyQ6Kx7Q7XdJ3OrbEQic0XG/ESOaupSQW5FEHpTan/pKRp/XotJkskKl9qxEvUG3IpgvJewVtY4cQdrs0Fh5gskVfypUa8RL0RFGSGXIqgtJX63+JwTNtdruX9TpZ4iYSISAYmpsxCVYPJ4ZiTDSZMSnFvERS53OEiIiJ5eeihKVCrKxyOEasIipil/pksERHJgFyKoJj9A2xnCrvDMv9ERL4nNTVBNkVQrpX6d8Q1d7m4DY9IZCzJTPbYiqBkLkCisR4x7X2WBAEnG0woUoWJUgSFZf6JCPC9xqPUs/YiKBkZS1BbO7dDYQUBanU5IiNzRSuCImapf4XQ0yVED3fu3DlRXker1XpFOV3Ow706Nh2d0uHNcFWDCYU/vhnu2MxNrvPoLc7DeVarFZWlJXivYxGUlNmY5sIiKI7mYbVasWp+elu/JztlzI1mCza2BPe639PIkSP7HK8vEGutEpO3/N73hrfMuXPj0em49oa4ApGROZ0aj3rLnHvDF+cMXJu31Wr9sQjKqZ8UQZkqWiJttVpxzz3r7PR8ateIqKje9Xzqbp1isuQi3vKLw3m4T1/ehMpxHn3BechLT/NoT+oTzd3f4epth3YmS44xWfIO3jDn3r4J9YY595YvzhmQ37yvJfXd3+XqzVrV3TrFbXhEImFJZvIULPNP5LvELMlM1B9ilfqXPFm6evUqtm/fjh9++AFDhgzB0qVLMWDAgE5jvv76a+zZswdNTU1QKpWYM2cOYmJiJIqYqG9Ykpk8Ccv8E/kmMUsyE/WXGKX+JU+W8vLycMcddyA1NRV5eXnIy8vDvHnzOo0JCAjAk08+iREjRqCurg7PPPMMxo4di+Dg7q56EMkPSzITEZHciVmSmcgTSJ4sVVdXY+3atQCA2NhYrF27tkuy1HEPYXh4OAYOHIgrV64wWfIR3lI9Ti5NR4mIyLXaK8fl5X2Cy5cFj64cd60ks/SNR4nkQPJk6fLlyxg0aBAAICwsDJcvX3Y4/syZMzCbzRg2bJgY4ZHEOlaPW9F+0NwkoEq/A6v26/t00FwqLMlM5Jmc2S4OABs2bMAXX3yBW265Bc8884wEkZIUOleOewDth8yPH69Adva6Xh8yl5qYJZmJPIEoydL69etRX1/f5fPp6emdPlYoFA6vul+6dAlZWVlYuHBht1dqysrKUFZWBgDYtGmTaH+gVCqVR/0x7I6c5mG1WrH6oV+3VY/rUBRBoVAgJjQQ48yNeGnZYrzyt+IuPw9ymke7+x6Yh4U5+zHOfLnbanhvB2rxyq9/bZuPHOfRF5yHvHjLPMTizHZxALjnnntgMplsaxB5P6vVioyMLDuV4xQwmeJRUzMZGRm9K18staSkWGRnr0NNzWR0Vw2vrfHoGrFDI5KEKMnS6tWru31s4MCBuHTpEgYNGoRLly4hNDTU7jij0YhNmzbh/vvvx89//vNun2/GjBmYMWOG7WOxShzKrZxiX8lpHuVFhUhoNkDTzZ0YjcoP/3XFgL/86U9dCiLIaR4dZW59uYemoy+jrq7ONl6u8+gtzkNepJiHJ5cOd2a7OADccccd+OSTT0SOjqTkjZXj5NR4lEgOJP9Jj46OxtGjRwEAR48exV133dVljNlsxpYtWzB16lRMnjxZ7BBJIu8VvIUpTlSPe7fgkEgR9V97SebWh5dic+BIbFMOxubAkTA//Fts2JfDq/1EMtTb7eLkO9oqx013OKatctwHIkXkGu0lmbOyvkd8/HLExDyD+Pjl2LXrAvLz13CtIp8i+Zml1NRUbN++HRUVFba94ADw5Zdf4u2338bjjz+OkydP4tNPP0VDQwOOHDkCAFi4cCFGjx4tXeDkdt5aPY4lmYnkx1XbxZ0l1ZZxMfnCdk+zWQNnKseZzRqP/Ld46KH78NBD9zkc4wvf55/yxTkDPjxvqQMICQnB888/3+XzERERiIiIAABMnToVU6dOFTs0khirxxGRWFyxXbw3pNoyLiZv2bbqiEplhDOV41Qqo9f+W/jC9/mnfHHOgPfPu7vt4pJvwyPqzsSUWahqMDkcc7LBhEkprB5HRO7jzHZx8k1paeOhVlc4HMPKcUSejckSyVbcTB0KVWEwmi12HzeaLShShWFaQoLIkRGRL0lNTcWHH36IxYsX46OPPkJqaiqAtu3i2dnZtnHPP/88tm3bho8++giPP/44Tp8+LU3AJJqkpFhERuYAaOxmRHvlOO6OIfJUkm/DI+qOUqnEsp27e6get5sVeYjIrZzZLg4AL7zwgphhkQywchyR92OyRLLWXj2usrQEmwvyoGptgdk/AJPunY0NCQlcgIiISFLtleOKio7i0KGVuHxZQFBQK9LTJ0Cn85z+SkRkH5Mlkj1WjyMiIjlTKpVITo7DQw/d59UH4Il8ES93EBERERER2cFkiYiIiIiIyA5uw7PDYrGgsqQY1YfzbWdkJqbMQtxMHfceExGR5CwWCwoLj+LgwRoYjf7QaFqRljYeSUmxXKeIiFyIydJPGAwGbM1cgCRzPVa0V18zCajS78Cq/Xos27nbJ7sXE5F34UUhz2UwGJCRkYXa2nSYTA+gvfra8eMVyM5eB71+EdcpIvJ4crkoxGSpA6vViq2ZC7AyoBEaTaDt8wqFAjGhgRhnbsTGzAXYsC+HbyaIyGNduHABzz14Py8KeSCr1YqMjCzU1OwAENzhEQVMpnjU1ExGRsYS5OezChsReS45XRTiX9IOKkuKkWSuh0blZ/dxjcoPieZ6HCktFTkyIiLXsFqtWPPwPKwMaERMaCAUCgWAaxeFVgY0YmvmAlitVokjJXsKC4+itjYdnROljoJRWzsXxcXHxAyLiMhlOl4Uuta7DLh2UWgHMjKyRFunmCx18F7BW5gSonY4JiZEjXcLDokUERGRa1WWFGOm6SIvCnmoAwdOwWSa7nCMyRSPnJwPRIqIiMi15HZRiNvwOlC1ttiusnZHoVBA1doiUkRE5M2kODf0XsFbWDEgwOGYmBA1NhccwnSdzi0xUN8Zjf64dpW1Owo0NfmLEQ4RebmO54bMZg1UKqPbzw21XRTa4nBM20Wh5UhKmuaWGDpistSB2T8AgklwmDAJggCzv+M3GkREPZGqmAwvCnk2jaYVgADHCZOAoKBWkSIiIm8l1bkhuV0U4ja8DiamzEJVg8nhmJMNJkxKmS1SRETkjToWkxH73JDZPwCCIDgcw4tC8pWWNh5qdYXDMWp1OdLTJ4gUERF5IynPDV27KOSIeBeFmCx1EDdTh0JVGIxmi93HjWYLilRhmJaQIHJkRORNpCwmMzFlFt656viuES8KyVdSUiwiI3MANHYzohGRkbnQ6aaKGRYReRkpzw3J7aIQk6UOlEollu3cjY0twThxpdl29VUQBJy40oyNLcFYtnM3y7ESUb9IWUwmbqYOJerBvCjkoZRKJfT6RYiKWgK1ugzXrr4KUKvLEBW1BHr9Iq5TRNQvUhaTkdtFIZ5Z+gmtVosN+3JQWVqCzQV5tkPXk+6djQ0JCVyAiKjfpDw3pFQqse71/Vg1Px2JxnrEtJ+XEgScbDChSBXGi0Iyp9VqkZ+/BkVFR5GbuxxNTf4ICmpFevoE6HTsr0RE/SfluaH2i0IZGUtQWzu3wzZAAWp1OSIjc0W9KMRkyQ6lUol4XSLidYlSh0JEXkjqYjJDhw7lRSEPp1QqkZwch+TkOKlDISIvJHUxGTldFGKyREQksokps1Cl34GY0MBux5xsMGHSve47N8SLQkRE1J20tPE4frzix7s69rn73JBcLgrx8iERkchYTIaIiORMbueGpMRkiYhIZCwmQ0REcsZiMtdwGx4RkQRYTIaIiOTsp+eGzGYNVCqjzxWTYbJERCQRnhsiIiI563huSKvVwmAwSB2S6HwjJSQiIiIiIuolJktERERERER2MFkiIiIiIiKyg8kSERERERGRHUyWiIiIiIiI7GCyREREREREZAdLh8ucxWJBZUkxqg/n2/qwTEyZhbiZOp+pb09ERPJlsVhQWHgUBw/WwGj0h0bTirS08XjwwV9JHRoRUb8xWZIxg8GArZkLkGSux4oQNRQKBQSTgCr9Dqzar8eynbuh1WqlDpOIiHyUwWBARkYWamvTYTI9AEABQMDx4xX4/e+X4ve/f4LrFBF5NN6akCmr1YqtmQuwMqARMaGBUCgUAACFQoGY0ECsDGjE1swFsFqtEkdKRES+yGq1IiMjCzU1O2AyxaMtUQIABUymeFRXb0FGRhbXKSLyaEyWZKqypBhJ5npoVH52H9eo/JBorseR0lKRIyMiIgIKC4+itjYdQHA3I4JRWzsXxcXHxAyLiMilmCzJ1HsFb2FKiNrhmJgQNd4tOCRSRERERNccOHAKJtN0h2NMpnjk5HwgUkRERK7HM0sS6alwg6q1xbb1rjsKhQKq1haRIiYiIl/SXeGGpKRYKJVKGI3+uLb1rjsKNDX5ixEuEZFbMFmSgDOFG8z+ARBMgsOESRAEmP0DRIyciIh8gaPCDdnZ66DXL4JG0wpAgOOESUBQUKsoMRMRuQO34YnM2cINdyXfg6oGk8PnOtlgwqSU2WKETUREPqKnwg01NTuQkZGF++4bB7W6wuFzqdXlSE+f4PaYiYjchXeW+uin2+j8ggdgbEJij/2PbIUbNIF2H9eo/JBorEcLFChUhWGcudFukQej2YIiVRg2JCS4bE5ERORdetpKZ4+zhRsUin8jMjIHNTWTuxnbiMjIXOh0a1w0GyIi8TFZ6gO72+gaDE71P3qv4C2scKJww+bCt7Bs525szFyARGM9YtpfRxBwssGEIlUYlu3czca0RGDzZiJ7nNlKZ2+taivcsMXhc5tM8ThwYDn0+kXIyFiC2tq5He5CCVCry3HnnW/i979fxN9BIvTtwgXJA5OlXuq4ja7j3aH2bXTjzI3YmLkAG/bl2P3h703hBq1Wiw37clBZWoLNBXm2N4GT7p2NDQkJ/OUiAps3E9nTcStd57s+7VvpJiMjYwny89d0WUt6U7hBq9UiP38NioqOIjd3OZqa/BEU1Ir09AmYN2876urqXDsxIg/U1wsXJA9MlnrJ2W10R0pLMV2n6/J4bws3KJVKxOsSEa9LdM0EiLxIfy9eEHmr3vRASkqa1umR3hZuUCqVSE6OQ3JyXKcR/J0j6t+FC5IHyb8rV69exfr167F48WKsX78eV69e7Xas0WjE448/jtdee03ECDvrb/+jiSmzWLiByEXYvJnIvv70QEpLG8/CDUQuwubNnk/yZCkvLw933HEHXn75Zdxxxx3Iy8vrdmxubi5uvfVW8YKzo7/9j+Jm6lCoCoPRbLH7eHvhhmks3EDUIzZvJrKvPz2QkpJiERmZA6Cxm69rL9wwtZ9REnk/Nm/2fJInS9XV1YiNjQUAxMbGorq62u64s2fP4vLlyxg7dqyY4XVh9g+AIAgOxzjqf6RUKtsKN7QE48SVZttzCYKAE1easbElmIUbiJzE5s1E9l3bSueI/R5ISqUSev0iREUtgVpd1uF5BKjVZYiKWgK9noUbiJzB5s2eT/IzS5cvX8agQYMAAGFhYbh8+XKXMVarFW+88QYWLVqEjz76SOwQO5mYMgtV+h2ICbV/Zgn4cRvdvd1vo2PhBiLXYPNmIvvS0sbj+PGKHyvU2edoK52jwg06Hc9WEDmLzZs9nyjJ0vr161FfX9/l8+np6Z0+VigUdt/0lJaWIioqCoMHD+7xtcrKylBWVgYA2LRpk8uri9z3wDwszNmPcebL3fY/ejtQi1d+/eseF5O58+Zj7rz5Lo2vv1QqlVdUZOE85MVd85hx/3y8s2u9w614VQ0t+K9HHnTJ6/P7QZ4iKSkW2dnr+tUDqbvCDUTkvP5euCDpiZIsrV69utvHBg4ciEuXLmHQoEG4dOkSQkNDu4z5/PPP8emnn6K0tBTNzc0wm80IDAzEAw880GXsjBkzMGPGDNvHBoPBNZPoIHPryz30P3rZY8ularVat/ybiY3zkBd3zWNCTAxWvToQYx00by5UDcSGyVNc8vr8fvTdyJEjRX09X9e+la67HkiRkbncSkckAldcuCBpSb4NLzo6GkePHkVqaiqOHj2Ku+66q8uYxYsX2/77yJEj+PLLL+0mSmKxt43OTzMA4+5N5jY6IhHZzgCyeTNRF9xKRyQ9XrjwfJInS6mpqdi+fTsqKiowZMgQLF26FADw5Zdf4u2338bjjz8ucYT2/bT/kbdccSbyNDwDSO529epVbN++HT/88INtnRowYECnMV9//TX27NmDpqYmKJVKzJkzBzExMRJFfA230hFJjxcuPJtC6Km0m4c7d+6cKK/jLckS5yEvnIe8cB5958nb8Pbv348BAwYgNTUVeXl5uHr1KubNm9dpzLlz56BQKDBixAjU1dXhmWeewfbt2xEc3F1vlc7EWqvE5C2/L73BOfsGX5wz4P3z7m6dYipLRETkgDMtLkaOHIkRI0YAAMLDwzFw4EBcuXJF1DiJiMj1JN+GR0REJGfOtLjo6MyZMzCbzRg2bFi3Y9xduVUOfLHqIufsG3xxzoAPz1vqAIiI5MZisaCypBjVh/NtZ6AmpszCfQ/M6/mLySP1t8VFu0uXLiErKwsLFy50eA5BjMqtUvP2LTv2cM6+QS5ztlgsKCw8ioMHa2A0+kOjaUVa2ngkJcW65RyUXObtLt1tw2Oy5GbdvemKm6njgT4iGTIYDNiauQBJ5nqsaK+uZxJQpd+BhTn7kbn1ZZ+8subt+tviAgCMRiM2bdqE+++/Hz//+c/dFarLif2Gi4j6z2AwICMjC7W16TCZHkB7hb3jxyuQnb0Oev0irlUuwr+CbmQwGPDcg/dDvXcnVpjO4bfCRawwnUOAfgdWzU/36uycyBNZrVZszVyAlQGNiAkNtN1BUCgUiAkNxFPKy9iauQBWq1XiSElM7S0uAHTb4sJsNmPLli2YOnUqJk+eLHaIfWYwGDBr1gvIzByB8vItqKrahPLyLVi8eDjuuWcd1ykiGbJarcjIyEJNzY4OpcgBQAGTKR41NTuQkZHFtcpFmCy5SU9vulYGNPJNF5HMVJYUI8lcb7fBLQBoVH5INNfjSGmpyJGRlFJTU/Hhhx9i8eLF+Oijj5CamgqgrcVFdnY2AODkyZP49NNPceTIEaxYsQIrVqzA119/LV3QTuAbLiLPVFh4FLW16bDf5BYAglFbOxfFxcfEDMtrcRuem9jedGkC7T6uUfkh0dj2pmu6TidydERkz3sFb2FFiNrhmJgQNTYXHOLvrQ8JCQnB888/3+XzERERiIiIAABMnToVU6dOFTu0funNG66kpGkiRkZEjhw4cAom0xaHY0ymeOTkLOfvrgvwzpKbvFfwFqY48abr3YJDIkVERD1RtbY4PLwPtN0dVrW2iBQRkfu0veGa7nBM2xuuD0SKiIicYTT649qd4O4o0NTkL0Y4Xo/JkpvwTReR5zH7B6CnPt2CIMDsHyBSRETuwzdcRJ5Jo2kF4HitAgQEBbWKEY7XY7LkJnzTReR5JqbMQlWDyeGYkw0mTEqZLVJERO7DN1xEniktbTzU6gqHY9TqcqSnTxApIu/GZMlN+KaLyPPEzdShUBUGo9li93Gj2YIiVRimJSSIHBmR6/ENF5FnSkqKRWRkDoDGbkY0IjIyFzqdZ52jlCsmS27CN11EnkepVGLZzt3Y2BKME1eabXeHBUHAiSvNeMk6EMt27mbvGfIKfMNF5JmUSiX0+kWIiloCtboM1+4QC1CryxAVtQR6/SKuVS7CanhuYnvTlbkAicZ6xLQ3txQEnGwwoUgVxjddRDKk1WqxYV8OKktLsLkgz9ZMetK9s/HKr3+Nuro6qUMkcon2N1wZGUtQWzu3Q/lwAWp1OSIjc/mGi0imtFot8vPXoKjoKHJzl6OpyR9BQa1IT58AnW4Nf29dSCH0dLDGw507d06U19FqtXab91mtVlSWluC9jm+6UmZjWkKCLH+Qu5uHp+E85IXzkBcp5jFy5EhRX8/TiLVW2WO1Wn98w3XqJ2+4pvZrnfKW35fe4Jx9gy/OGfD+eXe3TvHOkpsplUrE6xIRr0uUOhQiIqIulEolkpPjkJwcJ3UoRESyI79bG0RERERERDLAO0tEJCsWiwWVJcWoPpxv27o6MWUW4mbqZLl1lYiIfIvFYsHf/laBgwdrYDT6Q6NpRVraeCQlxXKd8kJMlohINgwGA7ZmLkCSuR4r2ouimARU6Xdg1X49lu3cDa1WK3WYRETkowwGA1JT1+PDD++FyfQA2ouiHD9egezsddDrF3Gd8jJMf4lIFqxWK7ZmLsDKgEbEhAZCoVAAABQKBWJCA7EyoBFbMxfAarVKHCkREfkiq9WKjIwsVFdv6VA9EgAUMJniUVOzAxkZWVynvAyTJSKShcqSYiSZ66FR+dl9XKPyQ6K5HkdKS0WOjIiICCgsPIra2nQAwd2MCEZt7VwUFx8TMyxyMyZLRCQL7xW8hSkhaodjYkLUeLfgkEgRERERXXPgwCmYTNMdjjGZ4pGT84FIEZEYmCwRkSyoWltsW++6o1AooGptESkiIiKia4xGf1zbetcdBZqa/MUIh0TCZImIZMHsH4CeemQLggCzf4BIEREREV2j0bQCcLxOAQKCglrFCIdEwmSJiGRhYsosVDWYHI452WDCpJTZIkVERER0TVraeKjVFQ7HqNXlSE+fIFJEJAYmS0QkC3EzdShUhcFotth93Gi2oEgVhmkJCSJHRkREBCQlxSIyMgdAYzcjGhEZmQudbqqYYZGbMVkiIllQKpVYtnM3NrYE48SVZtuWPEEQcOJKMza2BGPZzt1s+EdERJJQKpXQ6xfhrruWQ60uw7UteQLU6jJERS2BXr+I65SXYVNaIpINrVaLDftyUFlags0FeVC1tsDsH4BJ987GhoQELkBERCQprVaLY8e2Y9++N5GbuxxNTf4ICmpFevoE6HRruE55ISZLRCQrSqUS8bpExOsSpQ6FiIioC6VSieTkOCQnx0kdComAyRL5NIvFgsqSYlQfzrfdxZiYMgtxM3W8OkRERJKzWCwoLDyKgwdrYDT6Q6NpRVraeCQlxXKdIhIBkyXyWQaDAVszFyDJXI8VIWooFAoIJgFV+h1YtV+PZTt3Q6vVSh0mERH5KIPBgIyMLNTWpsNkegBtPX4EHD9egezsddDrF3GdInIzXpIgn2S1WrE1cwFWBjQiJjTQ1gxVoVAgJjQQKwMasTVzAaxWq8SREhGRL7JarcjIyEJNzQ6YTPG41gxVAZMpHjU1O5CRkcV1isjNmCyRT6osKUaSuR4alZ/dxzUqPySa63GktFTkyIiIiIDCwqOorU0HENzNiGDU1s5FcfExMcMi8jlMlsgnvVfwFqaEqB2OiQlR492CQyJFREREdM2BA6dgMk13OMZkikdOzgciRUTkm3hmiXrkjUUQVK0ttq133VEoFFC1togUERER9ZU3FkEwGv1xbetddxRoavIXIxwin8VkiRzy1iIIZv8ACCbBYcIkCALM/gEiRkVERL3lrUUQNJpWtDU9dZQwCQgKahUpIiLf5JmXW0gU3lwEYWLKLFQ1mByOOdlgwqSU2SJFREREveXNRRDS0sZDra5wOEatLkd6+gSRIiLyTUyWqFveXAQhbqYOhaowGM0Wu48bzRYUqcIwLSFB5MiIiMhZ3lwEISkpFpGROQAauxnRiMjIXOh0U8UMi8jnMFmibnlzEQSlUollO3djY0swTlxphiAIANq23p240oyNLcFYtnO3x+51JyLyBd5cBEGpVEKvX4SoqCVQq8vQtiUPAASo1WWIiloCvX4R1ykiN+OZJeqWtxdB0Gq12LAvB5WlJdhckGcrXjHp3tnYkJDABYiISOa8vQiCVqtFfv4aFBUdRW7ucjQ1+SMoqBXp6ROg063hOkUkAiZL1C1fKIKgVCoRr0tEvC5R6lCIiKiXfKEIglKpRHJyHJKT46QOhcgn8ZIEdYtFEIiISM5YBIGI3I3JEnWLRRCIiEjOWASBiNxN8m14V69exfbt2/HDDz9gyJAhWLp0KQYMGNBlnMFgQHZ2Ni5evAgAePbZZzF06FCxw/UptiIImQuQaKxHTHufJUHAyQYTilRhLIJARESSaS+CkJGxBLW1czuUDxegVpcjMjKXRRCIqF8kT5by8vJwxx13IDU1FXl5ecjLy8O8efO6jNu1axfmzJmDO++8E83NzT0WHiDXYBEEIiKSMxZBICJ3kjxZqq6uxtq1awEAsbGxWLt2bZdk6bvvvoPFYsGdd94JAAgMDBQ7TJ/GIghERCRnLIJARO4iebJ0+fJlDBo0CAAQFhaGy5cvdxlz7tw5BAcHY8uWLbhw4QLuuOMOPPDAA3avFpWVlaGsrAwAsGnTJmi1WvdO4EcqlUq013InzkNeOA954TyIiIh8iyjJ0vr161FfX9/l8+np6Z0+VigUdrfXWa1WfPrpp3jppZeg1Wqxfft2HDlyBNOnd21EN2PGDMyYMcP2scFg6P8EnKDVakV7LXfiPOSF85AXzqPvRo4cKerrERERuYIoydLq1au7fWzgwIG4dOkSBg0ahEuXLiE0NLTLmPDwcIwePRrDhg0DAEycOBGff/653WSJiIiIiIjIFSQ/9RgdHY2jR48CAI4ePYq77rqry5ibbroJRqMRV65cAQB8/PHHGDVqlKhxEhERERGRb5E8WUpNTcWHH36IxYsX46OPPkJqaioA4Msvv0R2djaAtoOb8+fPxwsvvIBly5ZBEIROW+2IiIiIiIhcTfICDyEhIXj++ee7fD4iIgIRERG2j++8805s2bJFzNCIiIiIiMiHSX5niYiIiIiISI6YLBEREREREdmhEARBkDoIIiIiIiIiueGdJRd55plnpA7BJTgPeeE85IXzIHKeL/6ccc6+wRfnDPjuvJksERERERER2cFkiYiIiIiIyA4mSy7iLX2fOA954TzkhfMgcp4v/pxxzr7BF+cM+O68WeCBiIiIiIjIDt5ZIiIiIiIiskMldQCe6urVq9i+fTt++OEHDBkyBEuXLsWAAQO6jDMYDMjOzsbFixcBAM8++yyGDh0qdrjdcnYeAGA0GvHb3/4Wd911Fx599FGRI3XMmXl8/fXX2LNnD5qamqBUKjFnzhzExMRIFHFnp0+fxuuvvw6r1Yr4+HikpqZ2ery1tRW7du3C2bNnERISgiVLlsjq56hdT/MoKChAeXk5/Pz8EBoaiieeeAJDhgyRJlgHeppHu3feeQfbtm3Diy++iIiICHGD7IEzczh58iQOHjwIhUKBG2+8EZmZmeIHSl7D0/8O94Wza+iGDRvwxRdf4JZbbvHYimLesk71Rk9zrq2txd69e/HNN99gyZIlmDx5sjSBupC3rOMuJVCf7Nu3Tzh06JAgCIJw6NAhYd++fXbHrVmzRvjHP/4hCIIgNDU1Cc3NzWKF6BRn5yEIgvCHP/xB2LFjh/D73/9epOic58w8/vWvfwnnzp0TBEEQLl68KPz3f/+3cPXqVTHDtMtisQhPPvmk8O9//1tobW0Vli9fLnz77bedxhQXFwu/+93vBEEQhOPHjwvbtm2TIlSHnJnHRx99ZPsdKCkp8dh5CIIgGI1G4fnnnxdWrlwpnDlzRoJIu+fMHM6dOyesWLFCaGhoEARBEOrr66UIlbyIJ/8d7itn19APP/xQqK6uFl588UURo3Mdb1mnesOZOX///ffC119/LWRlZQlVVVUSReo63rKOuxq34fVRdXU1YmNjAQCxsbGorq7uMua7776DxWLBnXfeCQAIDAyEWq0WNc6eODMPADh79iwuX76MsWPHihme05yZx8iRIzFixAgAQHh4OAYOHIgrV66IGqc9Z86cwfDhwzFs2DCoVCrExMR0if/999/HtGnTAACTJ0/Gxx9/DEFmxw2dmcftt99u+x0YM2YM6urqpAjVIWfmAQC5ubmYNWsW/P39JYjSMWfmUF5ejpkzZ9qugg8cOFCKUMmLePLf4b5ydg294447EBQUJGZoLuUt61RvODPnoUOH4sYbb4RCoZAoStfylnXc1Zgs9dHly5cxaNAgAEBYWBguX77cZcy5c+cQHByMLVu24KmnnsK+fftgtVrFDtUhZ+ZhtVrxxhtvYP78+WKH5zRn5tHRmTNnYDabMWzYMDHCc6iurg6DBw+2fTx48OAuf3w6jvHz84NGo0FDQ4OocfbEmXl0VFFRgXHjxokQWe84M4+zZ8/CYDBg/PjxYofnFGfmcO7cOZw/fx6rV6/GqlWrcPr0aZGjJG/jyX+H+6q3c/ZU3rJO9UZv1zRv4C3ruKvxzJID69evR319fZfPp6end/pYoVDYvapgtVrx6aef4qWXXoJWq8X27dtx5MgRTJ8+3V0h29XfeZSWliIqKqrTL5AU+juPdpcuXUJWVhYWLlwIpZLXC6Rw7NgxnD17FmvXrpU6lF5rv3iwYMECqUPpF6vVivPnz2PNmjWoq6vDmjVrsGXLFgQHB0sdGsmYL/4ddtWcibyJJ6/jvcVkyYHVq1d3+9jAgQNx6dIlDBo0CJcuXUJoaGiXMeHh4Rg9erTtqtnEiRPx+eefi54s9Xcen3/+OT799FOUlpaiubkZZrMZgYGBeOCBB9wZdhf9nQfQVqRi06ZNuP/++/Hzn//cXaH2Snh4uK0ACABcvHgR4eHhdscMHjwYFosFRqMRISEhYofqkDPzAIAPP/wQhw4dwtq1a2W5ha2neTQ3N+Pbb7/FunXrAAD19fV46aWX8NRTT8mmyIOzP1NjxoyBSqXC0KFDMWLECJw/fx433XST2OGSB/HWv8OOuGLOns5b1qnecHZN8ybeso67mrwv58hYdHQ0jh49CgA4evQo7rrrri5jbrrpJhiNRtt+7I8//hijRo0SNc6eODOPxYsX49VXX8Urr7yC+fPnY+rUqaInSj1xZh5msxlbtmzB1KlTZVWxJiIiAufPn8eFCxdgNptx8uRJREdHdxozYcIEHDlyBEBbBbbbbrtNdlcwnZnHV199hT179uCpp56S7RmZnuah0Wjw2muv4ZVXXsErr7yCMWPGyCpRApz7XkycOBGffPIJAODKlSs4f/68R2+HIul58t/hvnJmzt7AW9ap3nBmzt7GW9ZxV2NT2j5qaGjA9u3bYTAYOpUL/fLLL/H222/j8ccfB9CWfb/xxhsQBAH/8R//gd/85jdQqeRzQ8/ZebQ7cuQIvvzyS9mVDndmHseOHcOrr77aKWFduHAhRo8eLV3gPzp16hT27t0Lq9WKuLg4zJkzB7m5uYiIiEB0dDRaWlqwa9cufPXVVxgwYACWLFkiyze2Pc1j/fr1+Oc//4mwsDAAgFarxdNPPy1t0Hb0NI+O1q5di/nz58sqWQJ6noMgCHjjjTdw+vRpWwnnu+++W+qwyYN5+t/hvnB2DX3++efxr3/9C83NzQgJCcHjjz/ucWc9vGWd6o2e5nzmzBls2bIFjY2N8Pf3R1hYGLZt2yZ12P3iLeu4KzFZIiIiIiIisoPb8IiIiIiIiOxgskRERERERGQHkyUiIiIiIiI7mCwRERERERHZwWSJiIiIiIjIDiZLREREREREdsin4Q+RD1i4cCHq6+uhVCoRGBiIcePG4dFHH0VgYCAA4PTp0zh06BC++uor+Pv7Y9SoUfjlL39ptxHeP//5T+zbtw9nz55FQ0MDDhw4IPZ0iIjIC3GtIrqGd5aIRPb0009j3759+J//+R+cPXsWb775JoC2jufbt2/H1KlTkZ2djT179mDu3Ll4//337T6PSqXClClTujQOJiIi6i+uVURteGeJSCLh4eEYN24cvv32WwiCgL179+JXv/oV4uPjbWMiIyMRGRlp9+tHjhyJkSNH4t///rdYIRMRkY/hWkW+jneWiCRiMBhQU1OD0aNH49y5c7h48SImT54sdVhEREQ2XKvI1/HOEpHINm/eDD8/P2g0GowfPx5z5szB2bNnAQBhYWHSBkdERASuVUTtmCwRiWzFihW48847O30uJCQEAFBfX4+hQ4dKERYREZEN1yqiNtyGRyQDI0eOxODBg/HOO+9IHQoREZFdXKvIFzFZIpIBhUKBhx56CG+++SYqKythNBphtVrx2Wef4Xe/+53drxEEAS0tLTCbzQCAlpYWtLa2ihk2ERH5EK5V5Iu4DY9IJiZPnozAwED89a9/xR/+8AcEBATg+uuvxz333GN3/A8//IAnn3zS9vG8efMwZMgQvPLKK2KFTEREPoZrFfkahSAIgtRBEBERERERyQ234REREREREdnBZImIiIiIiMgOJktERERERER2MFkiIiIiIiKyg8kSERERERGRHUyWiIiIiIiI7GCyREREREREZAeTJSIiIiIiIjuYLBEREREREdnx/wOPEgRMIKmTjQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1008x432 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "kpca1 = KernelPCA(n_components=3, kernel='rbf', gamma=0.01)\n",
    "Xkpca1 = kpca1.fit_transform(Xstd)\n",
    " \n",
    "Xkpca2 = ker_pca(Xstd)\n",
    " \n",
    "with plt.style.context(('ggplot')):\n",
    "    fig, ax = plt.subplots(1, 2, figsize=(14, 6))\n",
    " \n",
    "    #plt.figure(figsize=(8,6))\n",
    "    ax[0].scatter(Xkpca1[:,0], Xkpca1[:,1], s=100, edgecolors='k')   \n",
    "    ax[0].set_xlabel('PC 1')\n",
    "    ax[0].set_ylabel('PC 2')\n",
    "    ax[0].set_title('Scikit learn')\n",
    "    \n",
    "    ax[1].scatter(Xkpca2[:,0], Xkpca2[:,1], s=100, facecolor = 'b', edgecolors='k')\n",
    "    ax[1].set_xlabel('PC 1')\n",
    "    ax[1].set_ylabel('PC 2')\n",
    "    ax[1].set_title('Our implementation')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9dce9c7f",
   "metadata": {},
   "source": [
    "# Let's Apply the TEP Dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "870b4c7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "train_normal = pd.read_csv('./data/normal_training.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b24e2793",
   "metadata": {},
   "outputs": [],
   "source": [
    "test_faulty = pd.read_csv('./data/faulty_testing.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "9dc03fff",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_train = train_normal[(train_normal.simulationRun>=1)].iloc[:500,3:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "821d3154",
   "metadata": {},
   "outputs": [],
   "source": [
    "IDV = 1\n",
    "df_test = test_faulty[(test_faulty.faultNumber==IDV) & \n",
    "                               (test_faulty.simulationRun==1)].iloc[:,3:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "84383b7b",
   "metadata": {},
   "outputs": [],
   "source": [
    "std_scl = StandardScaler().fit(df_train)\n",
    "dtrn_scl = std_scl.transform(df_train)\n",
    "dtrn_scl = pd.DataFrame(dtrn_scl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "f47b8e14",
   "metadata": {},
   "outputs": [],
   "source": [
    "dtrn_scl = dtrn_scl.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "c3786979",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(500, 52)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dtrn_scl.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "83a299a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "std_scl = StandardScaler().fit(df_test)\n",
    "dtst_scl = std_scl.transform(df_test)\n",
    "dtst_scl = pd.DataFrame(dtst_scl)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad59007c",
   "metadata": {},
   "source": [
    "### Train"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "ff9bc160",
   "metadata": {},
   "outputs": [],
   "source": [
    "gamma = 0.01\n",
    "n_comp = 16\n",
    "conf_T2 = 0.95\n",
    "nSamples = 500"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "167b28dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate euclidean distances of each pair of points in the data set\n",
    "dist = euclidean_distances(dtrn_scl, dtrn_scl, squared=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "5163edff",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(500, 500)"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dist.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "d3ca31de",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate Gaussian kernel matrix\n",
    "K = np.exp(-gamma * dist)\n",
    "Kc = KernelCenterer().fit_transform(K)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "17e8f1e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get eigenvalues and eigenvectors of the kernel matrix\n",
    "eig_vals, eig_vecs = np.linalg.eigh(Kc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "463342ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "eig_vals = pd.DataFrame(eig_vals)\n",
    "eig_vecs = pd.DataFrame(eig_vecs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "90e2fe36",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(500,)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eig_vals.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "bec5f10c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(500, 500)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eig_vecs.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "4de840a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the score matrix\n",
    "score_train = Kc@eig_vecs.iloc[:,:n_comp]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "1260e35c",
   "metadata": {},
   "outputs": [],
   "source": [
    "### T2 Process ###\n",
    "a = eig_vals.iloc[:n_comp,:]\n",
    "z = a.values.tolist()\n",
    "extract_egv=[]\n",
    "for i in range(len(z)):\n",
    "    extract_egv.append(z[i][0])\n",
    "## extracted eigen value\n",
    "c = np.diag(extract_egv)\n",
    "## eigen value with linear diagonal value\n",
    "invegvldiag = c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "e0eaabbc",
   "metadata": {},
   "outputs": [],
   "source": [
    "### T2 Score ###\n",
    "T2_train = score_train@invegvldiag@score_train.T\n",
    "T2_tfinal = np.diag(T2_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "b3e1dd4c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([9.71037800e-06, 1.78081108e-05, 1.91242739e-05, 1.07693525e-05,\n",
       "       2.26205809e-06, 4.56145652e-07, 7.09350842e-06, 3.78008491e-06,\n",
       "       1.56044916e-07, 1.19602359e-07, 4.01380024e-07, 4.46524492e-07,\n",
       "       2.40125349e-07, 3.86456764e-07, 1.21743371e-07, 1.84936978e-08,\n",
       "       9.32876375e-07, 6.45697811e-07, 2.71498107e-06, 1.02604039e-06,\n",
       "       3.97716804e-07, 4.72828137e-07, 3.30404799e-07, 5.39399803e-07,\n",
       "       8.24881799e-07, 3.82423144e-07, 1.27480228e-06, 9.32171118e-07,\n",
       "       9.50400275e-06, 1.42199333e-05, 2.63876161e-07, 1.23461078e-07,\n",
       "       4.67670475e-07, 7.19302618e-07, 1.52799701e-07, 6.32374289e-08,\n",
       "       3.42691212e-07, 3.01693659e-07, 3.33708143e-07, 3.33232491e-07,\n",
       "       2.11399969e-07, 1.90638827e-07, 2.70502883e-07, 1.16627597e-07,\n",
       "       7.81111958e-08, 9.37324754e-08, 7.73724315e-08, 2.24198894e-07,\n",
       "       2.91911562e-06, 1.89406615e-06, 4.57778426e-07, 1.92363441e-07,\n",
       "       6.33111823e-07, 4.86402925e-07, 1.58765230e-07, 1.58010025e-07,\n",
       "       3.07623927e-06, 3.52793907e-06, 5.59442726e-08, 1.43393302e-07,\n",
       "       8.24871843e-08, 8.55845668e-08, 2.06306628e-07, 1.24921226e-07,\n",
       "       5.28600668e-08, 8.23087490e-08, 1.17815702e-07, 2.14848164e-07,\n",
       "       1.82392353e-07, 1.78845067e-07, 1.15136853e-05, 1.97294341e-05,\n",
       "       2.82850617e-06, 9.84353664e-06, 3.18696348e-06, 5.11554186e-06,\n",
       "       5.59236763e-07, 1.47645931e-06, 3.15503507e-07, 6.42713274e-07,\n",
       "       6.28304802e-07, 4.93005337e-07, 5.86421007e-07, 5.81140329e-07,\n",
       "       7.94572390e-07, 4.15350164e-07, 1.29364009e-07, 5.01504750e-08,\n",
       "       2.04263853e-07, 2.52229827e-07, 6.30424327e-07, 4.17088608e-07,\n",
       "       1.17716444e-07, 2.31865172e-07, 1.91529905e-07, 2.47743389e-07,\n",
       "       4.19374189e-06, 1.88217802e-06, 1.50203114e-06, 1.84683895e-06,\n",
       "       1.06753221e-06, 8.50239832e-07, 1.11052958e-06, 9.32883966e-07,\n",
       "       2.14234441e-06, 1.50218214e-06, 1.14583375e-07, 7.26997529e-08,\n",
       "       1.68970842e-06, 7.79439424e-06, 6.07480827e-07, 1.12171815e-06,\n",
       "       1.38277104e-06, 5.79295720e-07, 1.31818443e-07, 1.16890999e-07,\n",
       "       1.64604387e-05, 1.79129030e-05, 6.77127783e-07, 5.38264953e-07,\n",
       "       3.41387991e-07, 2.07822004e-07, 1.05074067e-07, 1.41062394e-07,\n",
       "       7.48438812e-08, 1.49606877e-07, 3.13479050e-07, 1.66873540e-07,\n",
       "       1.17529334e-07, 5.72970136e-07, 2.08102945e-07, 9.21874637e-08,\n",
       "       1.12944030e-07, 1.10200724e-07, 2.09183943e-07, 1.76251121e-07,\n",
       "       8.31806248e-07, 8.44500684e-07, 9.07752802e-08, 1.63210558e-07,\n",
       "       3.29384644e-07, 8.79007763e-07, 9.45820619e-08, 1.02272519e-07,\n",
       "       5.69163439e-07, 2.20336697e-07, 8.28363993e-06, 3.42816648e-05,\n",
       "       4.74860928e-07, 7.45483991e-07, 3.57097692e-07, 2.76576987e-07,\n",
       "       1.62761352e-06, 7.81274737e-07, 5.65832842e-06, 3.58493401e-07,\n",
       "       3.44113340e-07, 1.94239044e-07, 9.39117871e-07, 1.14027406e-06,\n",
       "       6.85295512e-08, 6.87742699e-08, 1.70499625e-07, 8.88734377e-08,\n",
       "       1.25900759e-06, 5.28397450e-07, 2.28489973e-07, 1.88092485e-07,\n",
       "       2.99704598e-07, 2.00540871e-07, 3.61077515e-07, 3.31156049e-07,\n",
       "       2.29164446e-08, 7.71218337e-08, 2.80125462e-08, 6.16828633e-08,\n",
       "       4.76832506e-07, 3.82938763e-07, 4.92528552e-07, 4.43830289e-07,\n",
       "       3.06288889e-07, 3.01625513e-07, 1.38688698e-07, 1.42796383e-07,\n",
       "       1.61128511e-07, 3.50803771e-07, 2.04927849e-05, 2.02554349e-05,\n",
       "       1.93428130e-07, 2.34859918e-07, 2.29114866e-07, 1.77733649e-07,\n",
       "       6.67025235e-07, 1.58199087e-07, 9.82861898e-08, 3.05528893e-07,\n",
       "       7.79153614e-07, 1.07286746e-06, 1.20137033e-05, 1.15533663e-05,\n",
       "       2.83425511e-07, 1.55026013e-07, 3.96166449e-07, 4.43543746e-07,\n",
       "       4.22280020e-07, 1.80378935e-07, 1.05199221e-05, 1.26472061e-05,\n",
       "       1.74013301e-05, 1.46275197e-05, 2.62222906e-07, 8.34690838e-07,\n",
       "       9.33886364e-07, 1.94909460e-06, 5.06939121e-07, 2.08395148e-07,\n",
       "       1.00687805e-07, 3.47523047e-07, 1.21545457e-07, 9.40841878e-08,\n",
       "       1.10246535e-06, 1.54190889e-06, 2.90988821e-07, 2.55457550e-07,\n",
       "       8.07362330e-07, 2.87679572e-07, 1.00335489e-06, 1.45211754e-06,\n",
       "       5.85376995e-07, 2.80131360e-07, 1.34476723e-05, 1.81616279e-05,\n",
       "       2.07250099e-07, 6.77132768e-07, 3.12296599e-06, 2.82683726e-05,\n",
       "       5.60699655e-06, 1.11344714e-06, 8.17000183e-07, 7.00061116e-07,\n",
       "       3.83844479e-07, 1.85927897e-07, 3.59584701e-07, 7.31825085e-07,\n",
       "       1.20925129e-07, 2.19055948e-07, 4.43714466e-07, 5.92963709e-07,\n",
       "       4.79727964e-07, 1.08636149e-06, 3.93248004e-07, 3.81277222e-07,\n",
       "       3.51719054e-08, 1.89242172e-07, 3.92966576e-08, 4.77177442e-08,\n",
       "       1.31237634e-07, 2.22768664e-07, 2.89070270e-07, 3.70961463e-07,\n",
       "       1.73533413e-07, 2.20262566e-07, 2.12101692e-07, 2.06773281e-07,\n",
       "       1.31787756e-07, 1.03131197e-07, 6.03436011e-07, 7.73838420e-07,\n",
       "       3.21671158e-07, 3.44290353e-07, 8.69942565e-07, 1.37173803e-06,\n",
       "       1.76849887e-07, 3.72585243e-07, 4.00557405e-07, 1.57512902e-06,\n",
       "       1.48256434e-07, 3.18173761e-07, 5.01061714e-07, 6.98928228e-07,\n",
       "       7.48806520e-07, 2.78275595e-07, 6.38337108e-07, 3.72354091e-07,\n",
       "       9.79561382e-07, 1.17970142e-06, 1.41290724e-05, 1.47182203e-05,\n",
       "       3.45532679e-07, 2.14007465e-07, 2.62814663e-07, 1.61048993e-07,\n",
       "       2.89718303e-07, 4.52189716e-08, 7.04307680e-07, 3.07162152e-07,\n",
       "       1.22612374e-07, 2.19596339e-07, 1.45501951e-07, 1.27163131e-07,\n",
       "       7.12492230e-08, 1.16326801e-07, 1.03613029e-06, 6.16602891e-07,\n",
       "       6.95367865e-07, 3.03337076e-07, 2.29977799e-05, 2.23323059e-05,\n",
       "       1.53286247e-07, 1.13423481e-07, 1.40421387e-07, 8.66066421e-08,\n",
       "       1.93045813e-07, 8.18195577e-08, 5.32177346e-08, 1.52441891e-07,\n",
       "       3.66827344e-07, 8.13296721e-08, 2.54932504e-07, 1.09420318e-07,\n",
       "       9.22925917e-07, 1.20533375e-06, 1.55595859e-07, 1.75183073e-07,\n",
       "       1.19654922e-06, 5.30716939e-07, 2.75871999e-07, 4.91257689e-08,\n",
       "       6.06380018e-08, 9.44878619e-08, 1.22912529e-07, 1.27581779e-07,\n",
       "       2.73696736e-07, 1.02852797e-06, 2.29009898e-07, 5.46514276e-08,\n",
       "       3.59012875e-08, 6.74108468e-08, 1.63077094e-06, 1.21179256e-06,\n",
       "       8.34790794e-08, 4.86424267e-07, 6.99882771e-07, 7.88160690e-07,\n",
       "       1.55830884e-07, 5.98589436e-07, 1.45601034e-06, 1.13584361e-06,\n",
       "       2.12257905e-07, 4.36644548e-07, 1.16718017e-06, 7.81232431e-07,\n",
       "       1.17825064e-07, 5.51337877e-08, 5.12658385e-08, 4.12506432e-07,\n",
       "       6.22923234e-07, 1.29225571e-06, 2.32711238e-07, 3.35514638e-07,\n",
       "       3.69626162e-07, 1.69202708e-07, 1.85737159e-07, 5.41833414e-07,\n",
       "       2.41978680e-07, 2.13430779e-07, 8.41006980e-08, 4.03563426e-08,\n",
       "       2.64360337e-07, 1.55976583e-07, 4.03165346e-07, 2.12078331e-07,\n",
       "       8.46735590e-08, 1.69172204e-07, 6.05487666e-08, 6.33582793e-08,\n",
       "       1.31118914e-07, 2.68060492e-07, 8.05545845e-08, 8.93330275e-08,\n",
       "       1.59944758e-07, 5.47792067e-08, 1.15791741e-07, 2.13156526e-07,\n",
       "       3.56733740e-08, 3.97932807e-08, 3.78913845e-08, 4.27579167e-08,\n",
       "       1.14530275e-07, 9.60763359e-08, 2.89336365e-07, 4.41292444e-07,\n",
       "       2.95494439e-08, 7.43928938e-08, 4.23928027e-08, 3.65647134e-08,\n",
       "       1.08269334e-07, 4.33111808e-08, 1.58850485e-07, 1.12969968e-07,\n",
       "       2.01052887e-07, 3.05902591e-07, 1.57050667e-07, 1.94565274e-07,\n",
       "       8.00748314e-08, 1.41813020e-07, 5.69077786e-08, 1.57384000e-07,\n",
       "       1.13849948e-07, 5.41229617e-08, 1.71989592e-07, 1.92382396e-07,\n",
       "       3.39707626e-07, 2.92761225e-07, 2.45949996e-07, 2.81101708e-07,\n",
       "       2.41327090e-07, 2.77705636e-07, 3.11218911e-07, 5.85796647e-07,\n",
       "       1.88709270e-07, 1.53534668e-07, 7.50262495e-08, 1.27404462e-07,\n",
       "       7.73376268e-08, 1.12040308e-07, 3.20795083e-07, 1.07745145e-06,\n",
       "       6.42733960e-07, 4.00486585e-07, 1.09158814e-06, 7.43382412e-07,\n",
       "       7.77963869e-06, 1.00742631e-05, 4.39895559e-07, 6.82046071e-07,\n",
       "       1.94949911e-07, 1.38806422e-07, 1.95478866e-07, 8.25038153e-08,\n",
       "       3.15236973e-07, 3.86720256e-07, 2.73083122e-07, 7.12016493e-07,\n",
       "       2.52595025e-07, 5.36688522e-07, 4.93656095e-07, 2.38934163e-07,\n",
       "       3.22966213e-07, 3.46308980e-07, 2.72992090e-07, 3.57822434e-07,\n",
       "       1.37036875e-07, 4.05085240e-07, 1.15874601e-07, 4.37485034e-07,\n",
       "       4.17409669e-07, 8.09998839e-07, 2.31690370e-07, 1.70539956e-07,\n",
       "       9.28972789e-08, 1.83370064e-07, 4.73212343e-07, 8.26048414e-08,\n",
       "       1.33040355e-07, 2.53770904e-07, 2.46637574e-07, 1.01731422e-07,\n",
       "       5.55028888e-07, 2.27004343e-07, 1.62159864e-07, 3.89866862e-07,\n",
       "       3.90875662e-07, 2.31446407e-07, 4.66064703e-07, 2.32083171e-07,\n",
       "       4.45601517e-07, 7.16497428e-07, 6.07445440e-07, 2.12618727e-07,\n",
       "       3.72770303e-06, 1.41251761e-06, 5.57047393e-07, 4.83402187e-07,\n",
       "       1.04419947e-07, 9.45248953e-08, 2.40763621e-06, 1.33579412e-06,\n",
       "       4.68760375e-06, 1.47906240e-06, 3.75511153e-07, 2.69737790e-07,\n",
       "       7.85587998e-08, 3.86102165e-08, 1.95085394e-07, 5.15270795e-07,\n",
       "       2.67807729e-06, 1.22422169e-06, 3.77290444e-07, 1.70903030e-07])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T2_tfinal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "67f045ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Tlim ###\n",
    "F = f.ppf(conf_T2, n_comp, nSamples-n_comp)\n",
    "T2_trlim = ((n_comp*(nSamples**2-1))/(nSamples*(nSamples-n_comp)))*F"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "f27b6116",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "27.511202863723632"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T2_trlim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "17a1a69f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEMCAYAAADHxQ0LAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1F0lEQVR4nO2de5RU1Z3vv7+q7kYwDfSAKNI02D5IABMHiOJoIjGaq1kYJxDjI/PIjQ5mJll3XCtz52omYRwzN2bWGifOXXHGEOP1JiOIqEnU0WWiER8zgtCOhm6xFVugm1dD20AnIN1dte8f5+zT++w6p+rUu6vO97NWL6hTVefsfWqf7/79fvu39xalFAghhNQ/iWoXgBBCSGWg4BNCSEyg4BNCSEyg4BNCSEyg4BNCSEyg4BNCSEyg4BNCSEyg4JO6RETOF5FXRORFEVknIo3VLhMh1YaCT+qVXgCXKqU+CWAngKurWxxCqg8Fn9QlSql9Sqnj7sthAOlKl0FEukRkWaW+R0guKPikJhGRPxWR37p/H4hIynh9WEQmuJ+bA+AzAJ7I8/w7ReSyYr6jlFqglNpYru8Rki8UfFKTKKX+n1LqQ0qpDwH4LoAn9Wul1FSl1AkRmQzgpwC+rJQaqW6JCak+FHxSD5wH4A3zgIg0AHgIwN8ppbrDvigi/0tE9ojIkIh0i8inReSnANoAPOF6DH/tfvZWEXnX/eybIvJ54zwZ37Et96jXCvjebBF5TEQOisiAiPygFDeNxBClFP/4V9N/AHoArLCO/TGAAQAb3b9rA743D87g7unu67kAznT/vxPAZdbnrwFwOhxD6VoAvwMw03jf9x3zdT7Xsr6XhNOZfR/AyQBOAnBxte85/2rzr6EkvQYhVcIN28wF8Lp5XCn1UzjhnGykAEwAMF9EDiqldmb7sFJqg/FyvYjcBuB8AL+IUNS8rmVwPpxO5n8qpUbdYy9H/C4hPhjSIbXOxwAMAXgv3y8qpXYAuAXA7QD6ReQhETk97PMi8ici8ro7KHwYwEIA08txLYPZAHYZYk9IwVDwSa1zHoDfKKUK2slHKbVWKXUxgDkAFIB/0G+Zn3OzfX4E4OsApimlpgLoBCDm6UpxLYteAG3umAQhRUHBJ7XOebDCOVERkXkicqmbwvkBgOMYy9c/AKDd+PjJcIT5oPvd/w7Hwjexv1PotUxeBbAPwPdE5GQROUlELsqjmoR4UPBJrfMxFCj4cGLq3wNwCMB+ADMA3Oa+dyeAb7nhm79SSr0J4C4Ar8AR6HMB/Id1Pt93Cr2W+SWlVArAVQDOArAbQB+cAWNC8kYK9IQJIYTUGLTwCSEkJlDwCSEkJlDwCSEkJlDwCSEkJlDwCSEkJozryRzTp09Xc+fOrXYxSIF0dztrls2bN6/KJSEkXnR0dBxSSp1iHx/Xgj937lxs3bq12sUgBbJs2TIAwMaNG6taDkLihojsCjrOkA4hhMQECj4hhMQECj4hhMQECj4hhMQECj4hhMQECj4hhMQECj4hhMSEcZ2Hj+5uwM3l9vjiF4G/+Avg2DHgs5/N/M6Xv+z8HToEfOELme//+Z8D114L9PYCf/zHme9/4xvAVVc517755sz3v/Ut4LLLgNdfB265JfP9734X+IM/AP7zP4FvfjPz/bvvBs47D3j2WeDv/z7z/R/+EJg3D3jiCeCuuzLf/+lPgdmzgfXrgX/918z3H3kEmD4deOAB58/mqaeASZOAf/kX4OGHM9/XOfP/+I/Ak0/635s4EXj6aef/3/kO8Nxz/venTQMefdT5/223OfcIGPsNW1uBf/s35/+33DL2vuacc4A1a5z/r1oFvP22//3zznPuHwD80R8BfX3+9y+8ELjzTuf/K1cCAwP+9z/9aeDb33b+f+WVwPHj/veXLwf+6q/8ZTZh26udtvfKK/734972XGjhE0JITBjXG6AsWbJEcaZt7cKZtoRUBxHpUEotsY/TwieEkJhAwSeEkJhAwSeEkJhAwSeEkJhQMcEXkXYR+bGIPFKpaxJCCBmjKMEXkftFpF9EOq3jV4hIt4jsEJFbAUAp1aOUurGY6xFCCCmcYi38BwBcYR4QkSSAewBcCWA+gOtFZH6R1yGEEFIkRQm+UupFAO9bh88HsMO16IcBPATg6mKuQwghpHjKEcOfBaDXeN0HYJaITBORewH8vojcFvZlEVklIltFZOvBgwfLUDxCCIknFVtLRyk1AOCrET63BsAawJlpW+5yEUJIXCiHhb8HwGzjdat7jBBCSBUph+BvAXC2iJwhIk0ArgPweBmuQwghJA+KTctcB+AVAPNEpE9EblRKjQL4OoBnAGwH8LBSqqv4ohJCCCmGomL4SqnrQ44/BeCpYs5NCCGktHBpBUIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkUfEIIiQkVE3wRaReRH4vII5W6JiGEkDEiCb6I3C8i/SLSaR2/QkS6RWSHiNya7RxKqR6l1I3FFJYQQkjhNET83AMAfgDgJ/qAiCQB3APgcgB9ALaIyOMAkgDutL7/FaVUf9GlJYQQUjCRBF8p9aKIzLUOnw9gh1KqBwBE5CEAVyul7gSwvKSlJIQQUjTFxPBnAeg1Xve5xwIRkWkici+A3xeR27J8bpWIbBWRrQcPHiyieIQQQkyihnSKRik1AOCrET63BsAaAFiyZIkqd7kIISQuFGPh7wEw23jd6h4jhBAyDilG8LcAOFtEzhCRJgDXAXi8NMUihBBSaqKmZa4D8AqAeSLSJyI3KqVGAXwdwDMAtgN4WCnVVb6iEkIIKYaoWTrXhxx/CsBTJS0RIYSQssClFQghJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZQ8AkhJCZUTPBF5CMicq+IPCIif16p6xJCCHGIJPgicr+I9ItIp3X8ChHpFpEdInJrtnMopbYrpb4K4IsALiq8yIQQQgohqoX/AIArzAMikgRwD4ArAcwHcL2IzBeRc0XkSetvhvudzwH4dwBPlawGhBBCItEQ5UNKqRdFZK51+HwAO5RSPQAgIg8BuFopdSeA5SHneRzA4yLy7wDWFlxqQggheRNJ8EOYBaDXeN0H4IKwD4vIMgArAExAFgtfRFYBWAUAbW1tRRSPEEKISTGCnxdKqY0ANkb43BoAawBgyZIlqrylIoSQ+FBMls4eALON163uMUIIIeOQYgR/C4CzReQMEWkCcB2Ax0tTLEIIIaUmalrmOgCvAJgnIn0icqNSahTA1wE8A2A7gIeVUl3lKyohhJBiiJqlc33I8afAFEtCCKkJuLQCIYTEBAo+IYTEBAo+IYTEBAo+IYTEBAo+IYTEBAo+IYTEBAo+IYTEBAo+IYTEBAo+IYTEBAo+IYTEBAo+IYTEBAo+IYTEBAo+IYTEBAo+IYTEBAo+IYTEBAo+IYTEBAo+IYTEBAo+IYTEBAo+IYTEBAo+IYTEBAo+IYTEBAo+IYTEBAo+IYTEBAo+KZiOXYO45/kd6Ng1WO2iEEIi0FDtApDapGPXIL503yYMj6bR1JDAgzctxeI5LdUuFiEkC7TwSUFs6hnA8GgaaQWMjKaxqWeg2kUqCfRaSD1DC58UxNL2aWhqSGBkNI3GhgSWtk+rdpGKhl4LqXco+KQgFs9pwYM3LcWmngEsbZ9WF8IY5LWM93qt3bwbT3fuw5ULZ+KGC9qqXRwyzqHgk4JZPKelbILYsWsQj77WBwGwYlFrRYS31ryWtZt345s/2wYAeOmdQwBA0SdZqVvB79g1WBPWZ62Us5J07BrE9WtewXBKAQA2dPRh3Z+VN7yif4fVyxdg8NhwTfweT3fuy3hNwSfZqEvB79g1iOt/tMmz1MotFoVSjZhxLYQANvUMYMQVeyC/8EohHWitxu6vXDjTs+z1a0KyUZeC/9hrfRgeTQMAhkfTeOy1vnH5AFc6ZlwrIYCl7dPQmBTPwo8aXtHCfWIkjWRCcMfVCyPVbzzF7vPpsHTdxnsHTsYPdSn4KsfrqJQ73FLpmHGthAAWz2nBulUX5h3D39QzgBMjaSgAo2mF1b/oxLzTmnN+t2VSExIiAFRVY/drN+/G6l90Iq1UZE/jhgvaxuVvSMYnFRN8EVkG4DsAugA8pJTaWK5rrVzUike29mIkpdCYFKxc1Jr3OSrh5lc606WWQgC5BoSDOuOWSU0QAZTbw6eVymmtd+waxB1PdiGVVkgmBKuXL6iKdd+xaxCrf9GJ0bRT+OEayRIitUUkwReR+wEsB9CvlFpoHL8CwD8DSAK4Tyn1vSynUQB+C+AkAH0FlzgCtoVYCOPJzS8V9RIC0IO6ukNft+pCAMAdT3bB1UsIgKYI1rr+nRUApRQGjw2Xt/BZypFKj/miCZFxnyVEao+oFv4DAH4A4Cf6gIgkAdwD4HI4Ar5FRB6HI/53Wt//CoCXlFIviMipAP4JwJeKK3o4OqXvkY4+jKbSePS1vrwt9Eq4+aYX0ZAQXLNkdtlTEOshBPDoa31efH84pfDoa32YNXWiN24jAC4+ezpuueycnPdyvKRiLm2fhgmNCQyPpJFwxx9q3cAYr8Q5My6S4CulXhSRudbh8wHsUEr1AICIPATgaqXUnXC8gTAGAUwooKyRMAfutL2Ur4VeKTff9CKGUwprN+8uqHOKG7bXJnA6aG0gKzjhqij3cLxMIBsv5ah3ajUjq1QUE8OfBaDXeN0H4IKwD4vICgD/DcBUON5C2OdWAVgFAG1t+VuiposOOGKQr+VWKTdfW5e6c1IoffioXNZMNa2kFYtasaGjz7PKVyxqxaaeATj+mLNAVD6/WTknkOXDeClHPeMzsuokVJsPFRu0VUo9BuCxCJ9bA2ANACxZsiTvBBvTRU+6YZIFp0/xFveK8uNWys3XVp0OP6VSpb1euayZaltJi+e0YN2fZVrDExqrH5oh4xvTE0wr53WcKEbw9wCYbbxudY9VFds1BpC3OFXSvdZW3UrXSi3l9co18DweBrRta5ghERKFwWPDBXuC9UAxgr8FwNkicgYcob8OwA0lKVWRmGJwz/M7ChKnSrvX5bheuTyVsPNWezCsVPew2vUg5UMPjsfVE4yalrkOwDIA00WkD8DfKqV+LCJfB/AMnMyc+5VSXWUraYGMlyyMKJRaaMpl9QadNyjMU4tUO1xFykvcPcGoWTrXhxx/CsBTJS1RiamVH7hcQlMuT8U+b71siDIewlWkvIyXwfFqeJJ1ubSCzXj5gbNRqNCMl/BDkCe1oWqlKZxa8ghJ7VItTzIWgj9eRDEbhQjNeAo/1IonlYt6qQcZX9gaVC1Psu4Fv9qiGLWzKURoxlv4oRY8qSjUSz3I+CBIg6rlSda94FdTFPNZ/bAQL4Thh9JTC94gqS2CNOhrnzqrKp5k3Qt+LlEs1QNunyef1Q8LWRYXYPih1FTbGyT1SZgGBXmS5TY46l7ws4liqR7woPNEXf2w2GVxGX4oHeMtREbqg6DJoPc8v6NsepSNuhZ8s7f82qfOyni/VA940Hmirn5Y6mVxGZIoHIbISLnQhlk2Ua+EwVG3gh+ltyzVAx50nqjhllIui8uQhJ98Oz+GyEi5ySbqlTA46lbwo/SWpXrAzUXQxDpe7nV79Nr/+roMSTgU2vkxREbKSTZRr4TBUbeCH7W3LOUDrjdPz3dN+0LLoHd+0puBNCQFDcmEt+pmy6SmwFhhORhvoaTxHI8fb/dKM17LVU/kEvVyGxx1K/iVds9LJTD5PHSbegYwkhqL/6dSCl+8wFnA9NDQCdz+RBdGU+UP74zHUFKUBd4AVFzgxuO9ylYudgKlwb6P+t5WyiDT1K3gA8VZztXIic9XDJa2T0MyAbg7+0EBmDyhAQ+8stO349fwSBp3P/t2pC3/CiFXZ1cN0ci1wFtDQgCRinSIwNg92Hv4+Lj0PMLWQsrWHtkZRCPoue7eP1RQKnax1LXgF0Ixsd9iPYp8vYTFc1pw7cfb8ODm3QCc9b279h317fgFAGkAL71zCJvfex/r/qz0DSuos+vYNYg9h4+jMZGomkVrd/i++5ty9hgrxy5jNnZHY4bdTMOg3LuTtUxqwuCx4cDz699weCQNEUHLpCbf/Toxksa9L7yL82ZPLXifibhiP9ePvtaHh7f0BqZiMw+/Apg3uZjQTLHxt0K8hBWLWp1NvUfSEAFOjKQ8QUkmBG3TTsaO/t8CcBrWY6/1lbwhBeUZf+m+Teh7/xgggpkltmgLfSjMjem1hR+2y1i2a+R7fbNNpdIK154/G7OmTvR9Px9DI5/r23s8JwShIZvVyxd4VucdT3Zh9fIFaEgIhlNOx/irNw/gue0H0NSQwMpFrePSUwkiSodXTuznWoDAVGzm4ZcZneHySEef59qvXr4gw9KpFIV4CYvntGD18gX49s+3IaWAV3cOIiHAZR85FTdfciYefa3PE3wAPsu/1GXX5f2bn23zBEaUQkIEAlWSVLOoD0XQzGdzY/rbP7cQAPB0576MDc+zXaOQh9J+4Fcuas34zmOv9Xn3LJuA5nt987wAsoZsVi5qRVop7zODx4ZxzZLZWLt5d8b3FVATcxZydXiVIMgg0kaamYpd6GZN+RBbwbcbAuDc5M69R/CJs0/Br9/q9yydeac1Vy0MEYXBY8MwDAakFfDrt/px8yVnYuWiVjyytRcjKYXGpGDlota8zm2mfS44fUpOC6lj1yA2bO0d20RenAZdKssqigfWsWsQ1/9okydGev9bc2P6zr1HvKyqLTvf9/3Gm3oGvHYxPOK/RiEeoO6UgzqXoHuWTIYLaD7Xt88LOIKnBdo+V5iIa3FKG99fuai1LNtylhrzdwdQNY/Efq6DDDvm4ZcRuyEIgGRC8EhHn9f4geq7q1Hcd3vwFnBcRr1I07pVFxb0YNppn0B4SEC7y3sPH/dikwBwSvME3HBBW/4VD2Fp+zQ0JN1N6kOEUQs5MBbGWrGo1fcwHRo6EWpRt0xqGhMIAEPHR7xsinwfStuLtDuXjl2DuPvZt717JgC+sDjTAzDrH/X6m3oGfOe9bP6pXgxen9/2PLSI6xj+0vZpnjgNHR9B176jvk5rvAq9xhybSMO5D8V47qWKsQcZdszDLyPmg5NMCK5Z4qQzrnt1t68TqHQ+u0lU990evAWczitskaaojdZO+wSCQwKmu2wOSiZEMP1DE4q5BRllbpnUBCi3TEoFftY+quB/mFomNeH2xzszLGozkyYh8LymH738HpSRTRH1oQzzInXnohfNS6WVL9xgemFB6XxBk/yCsDuHr15ypq+82QTG9pCWtk/z2qLdaY1nzDoOHR/BfS+/h7RSuP3xTnTtPYIVRngt13NRSDivkNne5byvsRV8282+4YI2zxozO4EFp0/BHU92lfVHDvr+o6/1oWvPkcBUuaBzm4O3ZlwwKJYdtdEubZ+GxqRkWPh2SMB0l81Byfufn4zmk4pvYmaZEyJIK0cgtRdjlz8sjKUfpnue35FhUQPwZdIkEoK0W289wKYzLOwB1zDs+Lk2IHTnYi6aBwDnzpqC1VctiDReoD3RDR19oZlXUSzGIIGxPaQfvvAuPjZ7akamyXgP52h0+e5+9m1vjGI4pbB2825vkiSQO+so33Ce2eEn3WeylN5uIcRW8PUgnmmxAI5wivtvIQMpxYy0a6F/eGsvRi3LOpl0PI2wcwc93HY6oPZizFS7R62sHbuDWLfqwowYvnb3WyY1+dxlM767eE4LNpRA7AH/gwalkEhkHwTW5Q4TJNu7E4yJnO60Fs6agt/0HQkM+UXJ3bfj5w1JwbVLZvvalZmpAQDb9w+F1ttse0Ehq7BymIIe1RCxPaTn3urHsnkzfPfMvg+6vOOhAwgzcrxEAjh1NMN5QPbnAsgeTgu6t+ZY0Gha4ds/35bhVVSa2Ap+UG6sfpCaGhJY4VqF+cZssw34ZSPI/TdJK4WuvZkWv+2im6/NOmqLpjE5Zr0qOJaiFuiwzipbBsvq5Qu8TiDbwGwUdzmKQDca18wmLrrcQbMZzbDIIx19WPfq7oz8+Gs/3obuA10ZIT+dsZLrt7Xj519cMhv/+/Pn+urU2JDwhBsAUqloi2kFhaxykY8hsvD0KZ4oAs4g9+CxYc+g2Hv4ONa9ujv02almTn7YUuXaE00AOLd1CrbvH8pIyTVTUNdv6cXC06f4LPIwjyns3jpja+K1g5QCHty8G+u39AZa+5WYyFb3gh/U25vWaa5BvFwZFjZDx0d8A35RB4fs8IhNym2I+XQ+WjB0vXQYxLReTZGJ4rL6OpGRNJ7u3JdzBm+uNEft1aTcEMy6VRdGjjXrcwS9l+26ur6jqfD8+HmnNfvO+72ntkf+bYNSMTMwxiDMUFmueheSeRU1HKE9X9OzaXLHsey0QjOvvFTphMWKXlA9MwyGqxage/9QxjN9zZLZ3jhYKq2w+hedGQPsQVoSNHta36NFbVPRsfsw0u44DeBY+0HnrsREtroWfPsmfvnCubjv5fe8XOybLj4DzRMbMXR8BD96qScwLS4o9JPNjb/v5fe81wInZTIKZqMEHGvAJJEQLDx9Cla6sfpcA3ZApiUbZL2aIhPFm7GzHv5jxyFs2fl+1gYaJjb69/lgZMzKHU6pQHc6bDBr7ebd+PbPtyGtkNFZ5BK5XPnxdjjE/G0BJ4c/rD1kswa1SGjLLwHgorOmB3acQfXOFbIKIqqnqu+ZV66zp+PKhTMzxrFMI2jeac2+DsB8dvIpYylEL6ie9kB39/6hwHCughO206G2tFI+8TaNkts/t9A7h+kdJhOCF7r7cdcvu71B/4ak4PL5p+K5t/ozzh21rZaKuhZ82xpd81KP9yOMphXue/k93HH1Qnz/V2/7BPYjbgOwz5Hrh9jU49/MxMyUyYUpEC2TmpzGpIVQHLdaz37MZ1VOLRh2zrRtvdplyDbI9+BNS3H3s2/jP3Yc8u6tXqtH3wez3vZDqLOetGVkI4gmFh27Br0JZ0BmZ2F2TkGpeFE9h5ZJTXi6c19GzP3ldw5hc8+A18kEZdTYHofOfklaK5teuXCmJy5RHvSwDtAuu53dEyXF1/ytbrnsnKzhzy0738fq5Qsyxr4KEe9SiF62euoy64H/oPokxHlulXJmY7/Q3Y9/+tXbvt9+OKWwfsvujNnTAPDw1l68unPQV6ZUSuFjs6di2bwZvvVzzMwwM+KQTAj2Hj6Ojl2DJRf9uhZ8s/GKSMYDm1Yq8EHetucIvnTfpozd5ZMJwRu9h7HqJ1txSvOEjMGXpe3FbWZixp31AwTAFy99unNfQQ9FPgIRtDuYLSC3XHYOtux832fpb+4ZQBrwrKDmD0bRfFJDRmd2++OdGEkpJBPwZcMAjpW+4PQpkcRiU89Ahidkej46HGcuF2Bb5GH3JWygT59fh8h0JwOEZ3ms3bwbT3fuw4mRlNfBjaYUzjplEi5on5aRCRZljCIb2cZicp0vTDDtpQHMAc5v/Xyb9xk99lWIeBcy8SgszGJ7VXc/+3bowL8ZzhUA154/GweHTuC57QcyxFtz6uST0H1gyOcdbuoZQMpukIDP0zANLQAZ42Gde49440r5LrMehboW/CCr2fxhm1zLanPPgC/1MGh3ee3S/fLNA97n7JS4qFZUNoIGRc2H4MqFM7Fl5/slXZVz9fIFWVNPwwQkw9JP+a2gQ7894aVlarFZ9ZOt3udG08BZM07Gu/2/9cJp005uwsbu/lCxMB/olklNSMpY+KsxKZ7gaJxZyH5rTrcHLapAZoaJPaaiB/wuOns6JjYmfe1AZ/kEjQGt3bwb3/zZtsDfYcfB36Hv8HHnfhmeqO6gGozU4KAOIMwLKtZStjsGu10D/tm3eijCXASsEPEGMrPkstXV7JQT4pTz9b4jGdlDQUsrmMkG9pyMhadPcX+DzPIJHAG/+ZIzcfMlZ2aUyUxjNpc4MTVC/9/OABw8NoxZUyd640rlCO3UteAD/huse1c7o6Rr7xHfeiH6RzUnXDk/hL8FDI9mLjscxYrKhv2wmtkRutzFWoD2NXJ5DWECoi39zT3OBC0RBD4kmo5dg3hu+wHfsfbpJ6Nv8Jj3QO4/egL73zyAxqRA0iojJmwvb5xWQFKATxsPlikOtof2SMdYSqPAeUCDlkm2xyq0UOiw1ca3D3pituD0KYETuQAnzm9iZr8AjsB37jnihXfECDcMp9TYSqiSOcM5zKPIFsqKYhEHYbdr3dm//M4hrz7mImCbegYy2mmuBenM+qywJp8F1dXMiNPrSGnsdEuvw7bGSsw5GQDQ1jIR67fszvD6EwKs+kQ7mic2ZoRBzUwwM405V/pl0O8077Tmsi6vUPeCb2Jab+aPpictZZtwtXr5Ap81qYkyaJkPYYNOQKb7p2OP+eb12tfI5TXktNbcFSgTCUESTriisSGBZmOWrXarzecomRDPUvrGw69j58Ax7735MyfjMwtOQ8ukJq+egJENYixvDAAfmz01VAjNdMK1xmxkhfBlkm3v0O5g9fo8egJa2NIIVy6ciZfeOeRd8+ZPtuPoiVEnj33U6Uy27TmChoTguvPbvM5j2GpodmebzYoPC2UBmW1Ie735TgxaPKcFVy6ciVfeHfAlQdiLEUbpoIDsXknYe3bao4nZ6dpjEnYY1uzYdxz8XWB9GxKCyxec5pVH3wO9/IgOUV778TYvayrX8xn0O+Uzk7sQYiX42WKb9k0Ocre+84fn+gYJBcgIFeQaZDTd4qDPh4WFzEZ/YiTtK0e22ZZB6GtoS2Teac0+cdMN2iyj7dJri0anNioAKq3wxfPbvNTGb7zc4LvvZjgtaY1xrPrkmb7Qx7Ufb8O805pDt3BMhixvHDTAqMsDOPnVpkAkxBEHfU5zsMy0avVvp++fbfGOre8j2HFgCFf/4GVc+/E2T0DNGd2Ak1pphsNSaYXTp07EDRe0YWN3vy9kBGSuAZOrE7ZDWaa1a3p25sQgO1VQExZOMVcfveniM/DAKzv9GVfGYH6QaJvtK1t9wt5bPKcFd1y9EN/62TbYQ//LzjnFK2uuLQWDvBWbVFoFzjd49LU+X4jywc27sWFrL9KAFxFYv7XXN+ku1+/0tU+dVXKh18RK8HNZRUE9f5ClrafEm+KlrRoddw2KQdrhCO1R2JZVUFioZVKTZx1rF1ZTaKzPzvZZ2j4to4ympfa1T52Vc4xhpTFoN+QO2trx8FMnT8D/+PQ5vjoHCeM9z+8I3MLRFHD7QQ4K4Zh1uOPqhd76Nfre67TCoMGyoCW0gxaOg7vcw2hKeaGFN/q2eXWzLWdz4NtsYx27BrGxu9/7XDIhuOqjM/Hkb/ZlDDwHzazOJaC2Z/fKu2PeiZ0qCGQPp+jfVCmFV3oGfGIPwJe2a7eToFnj2YR5xaJWHBo6gVOa/WszzTutGee2TsEbfUd8xw8c/SCw4zYxV4INGsvT6BCvOVitO7OJjcmMz4+klK/jGLWWcQjzMCqxHHusBD+fgaQwS3vw2LAvvicCLJs3A89uPxC4RkegW5pSXjBCT7nOtRjV4LHhjPivpiEZPf1TE2QJ975/LLCM5qzSbGMMdubB/n1HMXfaydhz+LgvG2f/UWe/XbvOtjCaq1bqex00+xHwW6FhM0K19RSUkqo9lWyLwwHBa8nrNL8gnu7c55U3Sqqkvc7PtR93Ori02hs4hmIPYIYJqK5j0PiPnSqYrZ0EDcqKOGGpIMwJeua1g84bZNkGrWevny29TaAdbwf8mXZhYm96j00NCdz+uYXo3HsEh4ZOYOPbBz2vTxtwgHPtE2745+V3DqExKWhMimeY6AUETQsfQEbIUBMU1gFQto1aYiP4YQNJ2Qhy6VsmNWXEDac3T8iY0To8ksYdT3Rh4awpWLGoNcPyHE0rz2JPKQRONjLRKZ9miqn2MK5ZMjvvhmE/sOu39HqzARPi/Lljm74lgu0ZyqZ4Af7MA6UU3hv4HR56dXfGRLFcsz21YJudXFo54gTAJ6JBFrj2RoImBAVZfEHGgO2ZaEvPFiyd5qeUyuiQr1w40yunva6R9gLDPEu9zk/QPbfvl5l2aIZSvvaps7xVOYP2T73hgrbADlCfd8/h44FbMpphwfVben0dnuP1AgqCdFplTNADEHpem6ANXPSz1bn3qCf2AuCjrVMAOGKfbU0cIHMlWG28fNdd/iJsgHn18gX41s+3QSl4obDr3DBmy6QmdO494q07ZXceYfU0wzpmllY5ZtzGQvCLncFnf/+mi8/wllnVy9nqGbDmYNwbfUfwRt8RL4ZndjY/fOHdjNS+bNiDiHc82ZURRskHfb57X3gXz24/AKXGHpyLzpqOtt+b5FnHAvjqm63TNAULIoDbkBOAP4UypPGv3bwb33attoakoNE9l340zVgzEG6Bh1nQUe5vUP55MiBUZ6/v83TnPl8c+PL5p3odk9lBhHmBZjnMEJOZo223k6DZynYoxVyVczjEysyWNaMHlO0wpfb40uZkQzdjSmcxueP5gZOcws5rlsHewEXXzw7hJBOC1VctAADPclfwrxVlsrTdvxKs7SUHGXtL26dh8Niwb1Vucf/0d4PSm3NlQvmNL/+ksJpNyxSRTwD4knvN+UqpP6jUte0BTz1RJooQ2NbTyGgazRMbsf7mzKnti+e0eINx5oOvY3iNybFFuG6+5Exfap+dP26XQV9LT4oKs8qyEdTwnn+r39eAkwnxUg+1dWw3wsFjw4GTs/Q90MJ5zzMnY+fA75AUeKKoRSvoIbdnzo6mFD4zfwamN0/A+i29gVPewyxwszz6ocu1p4EtfHYnay6TEfSeme2UTAhmNE/w4sj6oTa9wLAHWgupGWJ6vrsfv3bDhhu29nqze3VqokZ3CObAbNj+qdnSB81nRg8oh4mVPdlw8NiwY0QAEOW0qXRaQURwaOhEpPPqMgRl4Ng0WAkA5raMupMJuse3f26htyxHmMEVNGal6yviTBzU4z5h+/yGjSGYZQkz5qqSliki9wNYDqBfKbXQOH4FgH8GkARwn1Lqe2HnUEq9BOAlEflDAFuKKXS+LG2f5lsJ7+GtvZGWuQ2KH9oDuGbmhv5X56YPWzE8bdlt2NqLa5bMxu1X5Q4vBXkn+rpRM4K04On0Mb3mzKYe/1IQAuDSD88A4Kw3Mu/UZpw6+SQsmzcjr0aoG/iGyRMwsSmJr3xmXqSOaVNP5szZ6c0T8N3Pn2tMhvHHmrNZ4NnuYdROUn8u2yqiUWdLBq1rFHYv7dDOc67YA047uuOJLqy+akFGaqKCI4BKOamxC2ZO9qVO3nG18/ia8eugLK+w8a4o4xAduwZ939VrWKWVwsbufjQkExgd9Q9SBhkjLZOafEJsj2EJHC/qZmtjlxWLWrFha2/WlS8BJ5SiO9+wvRVsY7Fz75HQMSKFwvf5NdtTIcZcVKJa+A8A+AGAn+gDIpIEcA+AywH0AdgiIo/DEf87re9/RSmlUw9uAHBjEWXOm8VzWny9fiqlkArIvbYxY7j2pI1cKzHqCRg6hqcbhCn8Expzi0/Q4GqupWiDymamj+nlABaePsVbLCrhWivPbj+AX791wNgu8QiWzZsRmLYZpTE2n9QQ6g3YLG2fhibLzdbhKjuLR187Ssgm6gbhYeSTI27PljTjyPovyl6wppDuPXzct5sZ4IQ09KDkTRefgXtf7PHe04sCamvRFPug7KewwcRs+ytkW7LB/u6mngHPQ0ylFS79yAw8b+wZDcBbbkMbI/qYLmYyAXz6w6f6niUAePGdg7j5kjMz7l2ulS+B6IsFmsaiDhHZY0R6vCXXGGGUiW65PIJiiCT4SqkXRWSudfh8ADuUUj0AICIPAbhaKXUnHG8gAxFpA3BEKTUU9H45sSdXBeVw29gNwpy0ERQmCouLmgOLpvBHER+7DILcS9EGCZTtsh4aOuETg0s/PJZpZHvR67fs9iZB3f7EmKUfJfd/6IPRyNtDmh2lHWoIW7U018OREQcuIPWtkBzxIJEwDYIoD7QZjtIWq4n+bZsnNvrW+jl6YhS3fvYj3gC6M66svJDU0vZpSCbFyyLJtpGMWc5sHV+u75reSv/RD3whwvVbdmcYI7OmTvR1Sum0M7nu5kvO9M1fCCvHikWtvjkXQSmnUcZ4Mo1F/xwC26OLGjGo1g5YxcTwZwHoNV73Abggx3duBPB/s31ARFYBWAUAbW2luxn2jwvkDotkaxB2z79ha2/g4JA+T/f+IXTtOYIJDQn8V+9hjKZU6CbcucodlHliEiZCGzrGvndK8wSfGOhMIz3j0OTNfUexbY8Te/dWphzNvtMS4Ij99n1HcdcvuyOHUsLEMB+xsb9nxoFT6eBF1HKVKawdhL0XJBKFWm1BHqNtrOgBSB3CmDyhAUdPjIZmwiTcf5MC3G5tqxj2XESxiMOWbtDCuGFrrzfgqkOkp04+CcDYIKweBDUHVXXu/qaegUjrSS2e0+LNuQhLOdWfAzJDsyamsWjOIdDCfemHZ0Ra/0aPt+jsHjP5oFwhHJuKZukopf42wmfWAFgDAEuWLMk9YpMHQa5nvt8xjy+bN8PLtBlJKfzwhXex5k+WZHzWXkArqZ+2kNztXGWIYpUEbdpiLgcAINAdfbpzn2c9AcBZMz6EnoO/DVwjJ1fpj34wApUl4yCf9dLzmUMR9L2wTJ6oZLPKg96zRaLYwTfT2p/ePCHDA7JDGPe+2OPlhNuZMHYnqC3/XGMduSziXIvyffLsU3xWu97DF8hcm2hTz4CXFy9A3iuKduwaxOCxYdx08Rno2nc0dPOiKOM7QSEqU7if234gUoqpPd6SVsGzd8sp+sUI/h4As43Xre6xuiRInKZbs/6ee6s/cA1rewGtVFr/GzxQlIsoYYwo4Y+w9D/TevrKRWd4A7ZJd0KJXv44VzqomSMdtAdoPgOpUdzvbN+LOlhaKgotbzayLTBmhzAAhGbChHWeUcY6srU92wszF+UbHknjN32HfZ9fMGtKhjHi7QVhtYugpU7CxoaCki3CNi+K6jna9bYHyr+wuDXn5vZBXkeUEG0pKUbwtwA4W0TOgCP018EZkK07wsRp5aJWrH91txfmUAFxQiBzAa3GpJOmVi7hyacR2+l/9szZxXMy1/COImJrN+/GXnfpXwXgyxfODY0HZ5sgY5e3kIdBfy/KYGkpKfXgW7bf1RQTPSkvaNtETdCGJeZYR5Rwo43dkejQi56Zuv/oCQBjq5SaBoO+V7aw6zraWUvZNgixJ8xlew4K8RyDhDssnGtjT3QDcodoS0nUtMx1AJYBmC4ifQD+Vin1YxH5OoBn4GTm3K+U6ipbSatI2IO2eE4LvvOH52aNE2rX8qufbEfXvqNYMHMyjp4YjbR8aqHk04iDPhsU+so3FGZ7NV37jmZcN9vAZjkotQBXmly/qxYTHesP2qQnyEvQc03CVvw0yZa/H+bV6Jmp+twfbZ2ChbOm5FVH21Nb9+pYenPQRkT20tbZYv2FeGLZZijnIt8QbSmJmqVzfcjxpwA8VdISjUOyPWi5pqZni2lmm2xlk0+8O59GXI7QA+B4NRus1/Z1SzmwWU3y3bu1UKL+VmZM2G5jdthGx5DtDUKCwnX2+jNB+fu2mNkzUxMCbN8/hG17jmSdaRw2QG56pGEzls1zBC1tHXRfi/Eci6WShkgsllYollwPWtgPli2mmU+8rpCJQ/k0onI0uBsuaMPfT/8QBn53At/9/LmB6WelHtisBsUu2xHl/PZEp2znzxb2CQrb6Bhy0FyToHPnyt+3sWfimum/5kJ0Qc9W9/6hjOP2IHzYeEOte3PlgoIfkUIaUFhMM1+BKzQlsdrMmDwBMyZPCM01Lpd3UUmKndSVjUI6k2zeqJmdo8M2dqcbJvb63HaqZJTUTDut+MV3xjJy3jkwhLt+2Q2lgAmN/k1ZTI/D3rui0oPw5aZSXiIFv4wECVohcb9CUxJrgVq2xNZu3o2HXt1d1EBnNgrp6LN1onY70mMm5mY4uc6dawu/bLNxNbp87xwYws9f3+sdPzEy5gVnG3Ct1iB8uSi3l2hCwS8zuQZAo56j1i3heqNj16CTEWMMRoYNdBZKoR19WBvL1o7szXAKjXdH6aT067t+2e07LhK8EUlYR1rLxoJJJT14Cn6NUC+N26RSbmw52NTjX3gumcg9LyFfytHRB7WjUgpO1E7KWV/Hf2zVJ9ox77RmfGTmZG82bjk60vFGJT14Cj6pCpV0Y8uBPRhpLs9bSirR0ZdScKJ2Ukvbp+GkRmfwVcQR+8sXnOZNmAKyZwzVE5X04Cn4pCqUc7CzEtRTmK3UdYnSSQVd01zsLVfGUL1RKQ+egk8qTilmdY4H6inMVo262Ne0PY24iH0lia3g13L8uNYJSg/kb0DqyWsar8RS8Gs9flzrBKUHEgLUl9c0Homl4NfqRKZ6gZYcIdUhloJfzxOZagVacoRUnlgKPi1MQkgciaXgA7QwCSHxI5H7I4QQQuoBCj4hhMQECj4hhMQECj4hhMQECj4hhMQECj4hhMQEUUrl/lSVEJGDAHYV+PXpAA6VsDi1AOscD1jneFBMnecopU6xD45rwS8GEdmqlFpS7XJUEtY5HrDO8aAcdWZIhxBCYgIFnxBCYkI9C/6aahegCrDO8YB1jgclr3PdxvAJIYT4qWcLnxBCiEFdCr6IXCEi3SKyQ0RurXZ5SoWI3C8i/SLSaRz7PRH5lYi84/7b4h4XEfk/7j34jYgsql7JC0NEZovI8yLypoh0ichfusfruc4nicirIvKGW+e/c4+fISKb3bqtF5Em9/gE9/UO9/25Va1AEYhIUkT+S0SedF/XdZ1FZKeIbBOR10Vkq3usrG277gRfRJIA7gFwJYD5AK4XkfnVLVXJeADAFdaxWwE8p5Q6G8Bz7mvAqf/Z7t8qAP9aoTKWklEA31BKzQewFMDX3N+ynut8AsClSqmPATgPwBUishTAPwD4vlLqLACDAG50P38jgEH3+Pfdz9Uqfwlgu/E6DnX+lFLqPCP9srxtWylVV38ALgTwjPH6NgC3VbtcJazfXACdxutuADPd/88E0O3+/4cArg/6XK3+AfgFgMvjUmcAkwC8BuACOBNwGtzjXhsH8AyAC93/N7ifk2qXvYC6troCdymAJ+Hsb1/vdd4JYLp1rKxtu+4sfACzAPQar/vcY/XKqUqpfe7/9wM41f1/Xd0H123/fQCbUed1dkMbrwPoB/ArAO8COKyUGnU/YtbLq7P7/hEAtbhn590A/hpA2n09DfVfZwXglyLSISKr3GNlbdux3fGqHlFKKRGpu7QrEfkQgEcB3KKUOioi3nv1WGelVArAeSIyFcDPAHy4uiUqLyKyHEC/UqpDRJZVuTiV5GKl1B4RmQHgVyLylvlmOdp2PVr4ewDMNl63usfqlQMiMhMA3H/73eN1cR9EpBGO2D+olHrMPVzXddYopQ4DeB5OOGOqiGgDzayXV2f3/SkABipb0qK5CMDnRGQngIfghHX+GfVdZyil9rj/9sPp2M9Hmdt2PQr+FgBnuyP8TQCuA/B4lctUTh4H8Kfu//8UTpxbH/8Td3R/KYAjhqtYE4hjyv8YwHal1D8Zb9VznU9xLXuIyEQ4Yxbb4Qj/F9yP2XXW9+ILAH6t3CBvraCUuk0p1aqUmgvnef21UupLqOM6i8jJItKs/w/gMwA6Ue62Xe2BizINhnwWwNtwYp9/U+3ylLBe6wDsAzACJ4Z3I5zY5XMA3gHwLIDfcz8rcLKV3gWwDcCSape/gPpeDCfO+RsAr7t/n63zOn8UwH+5de4EsNo93g7gVQA7AGwAMME9fpL7eof7fnu161Bk/ZcBeLLe6+zW7Q33r0vrVLnbNmfaEkJITKjHkA4hhJAAKPiEEBITKPiEEBITKPiEEBITKPiEEBITKPiEEBITKPiEEBITKPiEEBIT/j++NczIZv+6ywAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T2: 0.0\n"
     ]
    }
   ],
   "source": [
    "plt.figure(figsize=(6,4))\n",
    "plt.semilogy(T2_tfinal,'.')\n",
    "plt.axhline(T2_trlim, ls='--',c='r')\n",
    "plt.axvline(160, c='k')\n",
    "plt.title('$T^2$ statistic')\n",
    "plt.show()\n",
    "## print the T2 score\n",
    "print(f'T2: {(T2_tfinal>T2_trlim).sum()/T2_tfinal.shape[0]}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6643baa1",
   "metadata": {},
   "source": [
    "### Test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "ab2a556d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate euclidean distances of each pair of points in the data set\n",
    "dist = euclidean_distances(dtst_scl, dtrn_scl, squared=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "d95971cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate Gaussian kernel matrix\n",
    "K = np.exp(-gamma * dist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "96b42753",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(960, 500)"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "K.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "7cd211cd",
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "Kernel matrix must be a square matrix. Input is a 960x500 matrix.",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m~\\AppData\\Local\\Temp/ipykernel_19396/2500287805.py\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mKc\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mKernelCenterer\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfit_transform\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mK\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;32mc:\\users\\m477a\\appdata\\local\\programs\\python\\python39\\lib\\site-packages\\sklearn\\base.py\u001b[0m in \u001b[0;36mfit_transform\u001b[1;34m(self, X, y, **fit_params)\u001b[0m\n\u001b[0;32m    845\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0my\u001b[0m \u001b[1;32mis\u001b[0m \u001b[1;32mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    846\u001b[0m             \u001b[1;31m# fit method of arity 1 (unsupervised transformation)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 847\u001b[1;33m             \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfit\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mX\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mfit_params\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtransform\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mX\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    848\u001b[0m         \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    849\u001b[0m             \u001b[1;31m# fit method of arity 2 (supervised transformation)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mc:\\users\\m477a\\appdata\\local\\programs\\python\\python39\\lib\\site-packages\\sklearn\\preprocessing\\_data.py\u001b[0m in \u001b[0;36mfit\u001b[1;34m(self, K, y)\u001b[0m\n\u001b[0;32m   2231\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   2232\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mK\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m!=\u001b[0m \u001b[0mK\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 2233\u001b[1;33m             raise ValueError(\n\u001b[0m\u001b[0;32m   2234\u001b[0m                 \u001b[1;34m\"Kernel matrix must be a square matrix.\"\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   2235\u001b[0m                 \u001b[1;34m\" Input is a {}x{} matrix.\"\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mK\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mK\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mValueError\u001b[0m: Kernel matrix must be a square matrix. Input is a 960x500 matrix."
     ]
    }
   ],
   "source": [
    "Kc = KernelCenterer().fit_transform(K)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a52da7fe",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}