Interactive_PCA_Rotation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<p align=\"center\">\n",
    "    <img src=\"https://github.com/GeostatsGuy/GeostatsPy/blob/master/TCG_color_logo.png?raw=true\" width=\"220\" height=\"240\" />\n",
    "\n",
    "</p>\n",
    "\n",
    "### Interactive Workflow of Principal Component Analysis as a Rotation\n",
    "\n",
    "#### Michael Pyrcz, Associate Professor, University of Texas at Austin,\n",
    "    \n",
    "##### [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1) | [GeostatsPy](https://github.com/GeostatsGuy/GeostatsPy),\n",
    "\n",
    "#### Principal Component Analysis (PCA)\n",
    "\n",
    "This was a basic demonstration of PCA as a orthogonal transformation (a rotation as we preserve the inner product space, the pairwise distances and angles between all points).\n",
    "\n",
    "* by rotating the orthogonal components in 2D we can maximize the variance explained by a single component while reducing the absolute correlation between the components to 0.0. \n",
    "\n",
    "Appreciation to Professor Martin H. Trauth for the suggestion to also include the correlation and to comment on the decorrelation of components at the same time as maximizing the variance explained by the first component.\n",
    "\n",
    "First, some more on Principal Component Analysis one of a variety of methods for dimensional reduction:\n",
    "\n",
    "Dimensional reduction transforms the data to a lower dimension\n",
    "\n",
    "* Given features, $𝑋_1,\\dots,𝑋_𝑚$  we would require ${m \\choose 2}=\\frac{𝑚 \\cdot (𝑚−1)}{2}$ scatter plots to visualize just the two-dimensional scatter plots.\n",
    "\n",
    "* Once we have 4 or more variables understanding our data gets very hard.\n",
    "\n",
    "* Recall the curse of dimensionality, impact inference, modeling and visualization. \n",
    "\n",
    "One solution, is to find a good lower dimensional, $𝑝$,  representation of the original dimensions $𝑚$\n",
    "\n",
    "Benefits of Working in a Reduced Dimensional Representation:\n",
    "\n",
    "1. Data storage / Computational Time\n",
    "\n",
    "2. Easier visualization\n",
    "\n",
    "3. Also takes care of multicollinearity \n",
    "\n",
    "#### Orthogonal Transformation \n",
    "\n",
    "Convert a set of observations into a set of linearly uncorrelated variables known as principal components\n",
    "\n",
    "* The number of principal components ($k$) available are min⁡($𝑛−1,𝑚$) \n",
    "\n",
    "* Limited by the variables/features, $𝑚$, and the number of data\n",
    "\n",
    "Components are ordered\n",
    "\n",
    "* First component describes the larges possible variance / accounts for as much variability as possible\n",
    "* Next component describes the largest possible remaining variance \n",
    "* Up to the maximum number of principal components\n",
    "\n",
    "Eigen Values / Eigen Vectors\n",
    "\n",
    "* The Eigen values are the variance explained for each component. \n",
    "* The Eigen vectors of the data covariance matrix are the principal components and the Eigen  \n",
    "* Out of scope – just making the linkage\n",
    "\n",
    "#### Getting Started\n",
    "Here are the steps to get setup in Python with the GeostatsPy package:\n",
    "1.\tInstall Anaconda 3 on your machine (https://www.anaconda.com/download/).\n",
    "2.\tFrom Anaconda Navigator (within Anaconda3 group), go to the environment tab, click on base (root) green arrow and open a terminal.\n",
    "3.\tIn the terminal type: pip install geostatspy.\n",
    "4.\tOpen Jupyter and in the top block get started by copy and pasting the code block below from this Jupyter Notebook to start using the geostatspy functionality.\n",
    "\n",
    "You will need to copy the data file to your working directory. They are available here:\n",
    "   - Tabular data - unconv_MV_v2.csv at https://git.io/fjmBH.\n",
    "\n",
    "#### Install Packages\n",
    "\n",
    "For this interactive workflow to work, we need to install several packages relating to display features, widgets and data analysis interpretation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os                                                 # to set current working directory \n",
    "from sklearn.decomposition import PCA                     # PCA program from scikit learn (package for machine learning)\n",
    "from sklearn.preprocessing import StandardScaler          # standardize variables to mean of 0.0 and variance of 1.0\n",
    "import pandas as pd                                       # DataFrames and plotting\n",
    "import pandas.plotting as pd_plot                         # matrix scatter plots\n",
    "import numpy as np                                        # arrays and matrix math\n",
    "import matplotlib.pyplot as plt                           # plotting\n",
    "from matplotlib.ticker import AutoMinorLocator            # gridlines\n",
    "from matplotlib.gridspec import GridSpec\n",
    "import seaborn as sns\n",
    "import ipywidgets as widgets\n",
    "from ipywidgets import interactive, interact              # widgets and interactivity\n",
    "from ipywidgets import widgets                            \n",
    "from ipywidgets import Layout\n",
    "from ipywidgets import Label\n",
    "import matplotlib.transforms as transforms\n",
    "import math\n",
    "from ipywidgets import VBox, HBox\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "cmap = plt.cm.inferno\n",
    "plt.rc('axes', axisbelow=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you get a package import error, you may have to first install some of these packages. This can usually be accomplished by opening up a command window on Windows and then typing 'python -m pip install [package-name]'. More assistance is available with the respective package docs.  \n",
    "\n",
    "#### Declare Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_grid():\n",
    "    plt.gca().grid(True, which='major',linewidth = 1.0); plt.gca().grid(True, which='minor',linewidth = 0.2) # add y grids\n",
    "    plt.gca().tick_params(which='major',length=7); plt.gca().tick_params(which='minor', length=4)\n",
    "    plt.gca().xaxis.set_minor_locator(AutoMinorLocator()); plt.gca().yaxis.set_minor_locator(AutoMinorLocator()) # turn on minor ticks\n",
    "    \n",
    "def add_grid2(sub_plot):\n",
    "    sub_plot.grid(True, which='major',linewidth = 1.0); sub_plot.grid(True, which='minor',linewidth = 0.2) # add y grids\n",
    "    sub_plot.tick_params(which='major',length=7); sub_plot.tick_params(which='minor', length=4)\n",
    "    sub_plot.xaxis.set_minor_locator(AutoMinorLocator()); sub_plot.yaxis.set_minor_locator(AutoMinorLocator()) # turn on minor ticks   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Mutlivariate Dataset\n",
    "\n",
    "Let’s load multivariate dataset, from subsurface energy with 1,000 unconventional wells including:\n",
    "\n",
    "* porosity\n",
    "* log transform of permeability (to linearize the relationships with other variables)\n",
    "* accoustic impedance (kg/m^3 x m/s x 10^6)\n",
    "* brittness ratio (%)\n",
    "* total organic carbon (%)\n",
    "* vitrinite reflectance (%)\n",
    "* initial production 90 day average (MCFPD).\n",
    "* scaled production\n",
    " \n",
    "all samples have the support volume of a well (one measure per well).\n",
    " \n",
    "Note, this dataset is available on my GitHub in my [GeoDataSets](https://github.com/GeostatsGuy/GeoDataSets) repository."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAGoCAYAAAC+M2lpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABr00lEQVR4nO3deVwU9f8H8NeK3CCKihcIKElKoimEaB54XyRqZloeeStqYt+ALEFRU8lETTwqbyu10sqjlDxIMxAQxbxRFBXJROVSQOHz+4N2fix7MLPs7Owu7+fjwQOZ/ezn2B33vTPz+bxHxhhjIIQQQgSoJXUHCCGEGB8KHoQQQgSj4EEIIUQwCh6EEEIEo+BBCCFEMAoehBBCBKPgQQghRDAKHoQQQgSj4EEIIUQwCh4SkMlkvH5OnDih974VFxdj7dq1eP3111GvXj1YWFigWbNmeOuttxAfH6/3/lSlR48e6NGjh1bPXbduHbZu3arT/sjJZDIsWLCgynJ37tzBjBkz0KpVK1hbW8PR0RFt27bF5MmTcefOHVH6Jvftt99i1apVKh/j2//qPFfdft+gQQOt2uXj008/xU8//SRa/TVJbak7UBP99ddfCn8vWrQIx48fx7FjxxS2t2nTRp/dwsOHD9G/f3+kpaVhwoQJ+PDDD+Ho6Ih79+7h559/Rq9evZCSkoJ27drptV9iWbduHRo0aIDx48dL0v7du3fRoUMH1K1bFx988AE8PT2Rm5uLS5cuYc+ePbh58yZcXFxEa//bb7/F33//jTlz5ig99tdff8HZ2Vm0tuXefPNNfPDBBwrbzM3NRWvv008/xZtvvomgoCDR2qgpKHhIoFOnTgp/N2zYELVq1VLarm9jx47F+fPncfjwYfTs2VPhsbfffhtz585FvXr1qt3O8+fPIZPJULu28u739OlT2NjYVLsNY/DVV1/h4cOHOHPmDNzd3bntQUFBmDdvHsrKyiTrm772xUaNGkm+3+vCs2fPYG1tLXU39IpOWxmgESNGwMvLS2FbYGAgZDIZvv/+e27b2bNnIZPJsH//fm7b33//jSFDhqBevXqwsrJC+/btsW3btirbTElJwa+//oqJEycqBQ45X19fNG/eXFBbJ06cgEwmw44dO/DBBx+gWbNmsLS0RHp6OsaPHw87OztcuHABffv2hb29PXr16gUAKCkpweLFi/Hyyy/D0tISDRs2xHvvvYd///23yrEsXLgQfn5+cHR0RJ06ddChQwds2rQJFXOAurm54eLFi4iPj+dOl7i5uXGP5+Xl4X//+x/c3d25U3dz5sxBYWGhQlt5eXmYPHky6tevDzs7O/Tv3x/Xrl2rso8AkJOTg1q1asHJyUnl47Vq/f9/T/lrdfHiRfTq1Qu2trZo2LAhZs6ciadPnyo8LzY2Ft26dYOTkxNsbW3Rtm1bREdH4/nz51yZHj164ODBg7h9+7bCKSO5yqee/v33X8yYMQNt2rSBnZ0dnJyc0LNnT5w8eZLXWLV1/fp1jB49Gk5OTrC0tETr1q0RGxurUKaoqAgffPAB2rdvDwcHBzg6OsLf3x8///yzQjmZTIbCwkJs27aNG6/8lOeCBQsUxi+3detWyGQy3Lp1i9vm5uaGwYMHY+/evXj11VdhZWWFhQsXAgCys7MxdepUODs7w8LCAu7u7li4cCFevHihUO/69evRrl072NnZwd7eHi+//DLmzZung1dMf+jIwwD17t0bP/zwA+7fv48mTZrgxYsXiI+Ph7W1NeLi4jBixAgAwO+//47atWtz/wGuXr2Kzp07w8nJCWvWrEH9+vWxc+dOjB8/Hv/88w9CQ0PVtnnkyBEA4H04L7Stjz76CP7+/tiwYYPCB2ZJSQneeOMNTJ06FeHh4Xjx4gXKysowZMgQnDx5EqGhoejcuTNu376NyMhI9OjRA8nJyRq/5d26dQtTp07lAl1CQgJmzZqFe/fuISIiAgCwb98+vPnmm3BwcMC6desAAJaWlgDKj366d++Ou3fvYt68efD29sbFixcRERGBCxcu4Pfff4dMJgNjDEFBQTh9+jQiIiLg6+uLP//8EwMGDOD1Gvr7+yM2NhbDhg3D3Llz4e/vjzp16qgt//z5cwwcOJB7rU6fPo3Fixfj9u3bCl8gbty4gdGjR3OB7/z581iyZAmuXLmCzZs3Ayg/ZTdlyhTcuHED+/btq7Kvjx49AgBERkaicePGKCgowL59+9CjRw8cPXpU6+tOjDGlD1YzMzPIZDJcunQJnTt3RvPmzfH555+jcePGOHz4MGbPno2HDx8iMjISQPl1ukePHuF///sfmjVrhpKSEvz+++8YNmwYtmzZgrFjxwIoPxXXs2dPBAQEYP78+QCg8fXW5OzZs7h8+TI++eQTuLu7w9bWFtnZ2XjttddQq1YtREREoGXLlvjrr7+wePFi3Lp1C1u2bAEA7Nq1CzNmzMCsWbOwYsUK1KpVC+np6bh06ZJWfZEMI5IbN24cs7W15f5OT09nANj27dsZY4ydOnWKAWChoaHM3d2dK9enTx/WuXNn7u+3336bWVpasszMTIX6BwwYwGxsbNiTJ0/U9mHatGkMALty5QqvPvNt6/jx4wwA69atm8pxA2CbN29W2P7dd98xAOzHH39U2J6UlMQAsHXr1nHbunfvzrp37662n6Wlpez58+csKiqK1a9fn5WVlXGPeXl5qXzu0qVLWa1atVhSUpLC9h9++IEBYIcOHWKMMfbrr78yAGz16tUK5ZYsWcIAsMjISLX9YoyxsrIyNnXqVFarVi0GgMlkMta6dWsWEhLCMjIyFMrKXyt1bZ06dUrj+Ldv387MzMzYo0ePuMcGDRrEXF1dVT6vqv6/ePGCPX/+nPXq1YsNHTpU0HMrllP189VXXzHGGOvXrx9zdnZmubm5Cs+bOXMms7KyUhiLqr5NnDiRvfrqqwqP2drasnHjxik9JzIykqn6ONyyZQsDoPB+uLq6MjMzM3b16lWFslOnTmV2dnbs9u3bCttXrFjBALCLFy9y/a9bt67qF8WI0GkrA9SyZUu4ubnh999/BwDExcWhbdu2ePfdd5GRkYEbN26guLgYp06dQu/evbnnHTt2DL169VK6yDp+/Hg8ffpU6UJ9dQhta/jw4WrrqvzYgQMHULduXQQGBuLFixfcT/v27dG4ceMqZ6EdO3YMvXv3hoODA8zMzGBubo6IiAjk5OTgwYMHVY7twIEDeOWVV9C+fXuF9vv166cwC+748eMAgHfeeUfh+aNHj66yDaD8NMqGDRtw8+ZNrFu3Du+99x6eP3+OmJgYeHl5qZzdpq4teV8AIDU1FW+88Qbq16/PjX/s2LEoLS3lfUpNlQ0bNqBDhw6wsrJC7dq1YW5ujqNHj+Ly5cta1/nWW28hKSlJ4ScoKAhFRUU4evQohg4dChsbG4X3YeDAgSgqKkJCQgJXz/fff48uXbrAzs6O69umTZuq1TdNvL290apVK4VtBw4cQEBAAJo2barQX/mRqPz9fO211/DkyROMGjUKP//8Mx4+fChKH8VGwcNA9erVC0ePHgVQfnqqT58+aNu2LRo1aoTff/8df/75J549e6YQPHJyctCkSROlupo2bco9ro78FE9GRgav/gltS1VZALCxsVE6dfDPP//gyZMnsLCwgLm5ucJPdna2xv9sZ86cQd++fQGUX5D+888/kZSUhI8//hhA+YXNqvzzzz9IS0tTatve3h6MMa79nJwc1K5dG/Xr11d4fuPGjatsoyJXV1dMnz4dmzZtwvXr17F7924UFRXhww8/VCinqS35652ZmYmuXbvi3r17WL16NU6ePImkpCTuOgGf8auycuVKTJ8+HX5+fvjxxx+RkJCApKQk9O/fX+s6gfLJIj4+Pgo/DRo0QE5ODl68eIEvvvhC6X0YOHAgAHDvw969e/HWW2+hWbNm2LlzJ/766y8kJSVhwoQJKCoq0rpvmqjan//55x/s379fqb/y65fy/o4ZMwabN2/G7du3MXz4cDg5OcHPzw9xcXGi9FUsdM3DQPXq1QubNm3CmTNnkJiYiE8++QQA0LNnT8TFxeH27duws7NTmKlSv3593L9/X6murKwsANA4f75fv36YN28efvrpJ/Tv37/K/gltS9XFSHXbGzRogPr16+O3335T+Rx7e3u1/dq1axfMzc1x4MABWFlZcduFzO1v0KABrK2tuesDqh4Hyl+DFy9eICcnR+FDPTs7m3dbqrz11ltYunQp/v77b4XtmtqSb/vpp59QWFiIvXv3wtXVlSt37ty5avVp586d6NGjB9avX6+wPT8/v1r1qlOvXj2YmZlhzJgxCA4OVllGPkNt586dcHd3x+7duxX2p+LiYt7tyfeV4uJi7toXALVfVNTtt97e3liyZInK58i/WAHAe++9h/feew+FhYX4448/EBkZicGDB+PatWsK75sho+BhoHr16gWZTIb58+ejVq1a6NatG4Dyi+kffvghbt++jW7duinMie/Vqxf27duHrKwshR11+/btsLGx0TglskOHDhgwYAA2bdqEt956S+WMq+TkZDg5OaF58+bVaqsqgwcPxq5du1BaWgo/Pz9Bz5VPATYzM+O2PXv2DDt27FAqa2lpqfJb8+DBg/Hpp5+ifv36ClNoKwsICEB0dDS++eYbzJ49m9v+7bff8uqrfEJEZQUFBbhz547C6yqnri35BWv5h1rFD0DGGL766iulutSNXxWZTKZQJwCkpaXhr7/+EmUtio2NDQICApCamgpvb29YWFho7JuFhYXCB3p2drbSbCtA/ZjlM+3S0tLg6+vLba84EaEqgwcPxqFDh9CyZUveU9ptbW0xYMAAlJSUICgoCBcvXqTgQarHyckJr7zyCo4cOYKAgABu7UPv3r3x6NEjPHr0CCtXrlR4TmRkJHfeNSIiAo6Ojvjmm29w8OBBREdHw8HBQWOb27dvR//+/TFgwABMmDABAwYMQL169XD//n3s378f3333HVJSUtC8efNqt6XJ22+/jW+++QYDBw7E+++/j9deew3m5ua4e/cujh8/jiFDhmDo0KEqnzto0CCsXLkSo0ePxpQpU5CTk4MVK1YoffABQNu2bbFr1y7s3r0bLVq0gJWVFdq2bYs5c+bgxx9/RLdu3RASEgJvb2+UlZUhMzMTR44cwQcffAA/Pz/07dsX3bp1Q2hoKAoLC+Hj44M///xTZaBSZcmSJfjzzz8xcuRItG/fHtbW1sjIyMDatWuRk5ODzz77TKG8hYUFPv/8cxQUFMDX15ebbTVgwAC8/vrrAIA+ffrAwsICo0aNQmhoKIqKirB+/Xo8fvxY5fj37t2L9evXo2PHjqhVqxZ8fHxU9nXw4MFYtGgRIiMj0b17d1y9ehVRUVFwd3dXmi2lK6tXr8brr7+Orl27Yvr06XBzc0N+fj7S09Oxf/9+blGtfNrsjBkz8Oabb+LOnTtYtGgRmjRpguvXryuN+cSJE9i/fz+aNGkCe3t7eHp6YuDAgXB0dMTEiRMRFRWF2rVrY+vWrYJW+UdFRSEuLg6dO3fG7Nmz4enpiaKiIty6dQuHDh3Chg0b4OzsjMmTJ8Pa2hpdunRBkyZNkJ2djaVLl8LBwUEhcBk8qa/YE+XZVnIhISEMAFuyZInC9pdeeokBYGlpaUrPuXDhAgsMDGQODg7MwsKCtWvXjm3ZsoV3X549e8bWrFnD/P39WZ06dVjt2rVZ06ZN2bBhw9jBgwcFtyWfbfX999/zHjdjjD1//pytWLGCtWvXjllZWTE7Ozv28ssvs6lTp7Lr169z5VTNttq8eTPz9PRklpaWrEWLFmzp0qVs06ZNSrNmbt26xfr27cvs7e0ZAIWZRwUFBeyTTz5hnp6ezMLCgjk4OLC2bduykJAQlp2dzZV78uQJmzBhAqtbty6zsbFhffr0YVeuXOE14yghIYEFBwezdu3aMUdHR2ZmZsYaNmzI+vfvz83oqvxapaWlsR49ejBra2vm6OjIpk+fzgoKChTK7t+/n3vdmjVrxj788ENuZtjx48e5co8ePWJvvvkmq1u3LpPJZAqzjSr3v7i4mP3vf/9jzZo1Y1ZWVqxDhw7sp59+YuPGjVOascVn7PJywcHBGstkZGSwCRMmsGbNmjFzc3PWsGFD1rlzZ7Z48WKFcsuWLWNubm7M0tKStW7dmn311VcqZ1CdO3eOdenShdnY2DAACvvOmTNnWOfOnZmtrS1r1qwZi4yMZF9//bXK2VaDBg1S2d9///2XzZ49m7m7uzNzc3Pm6OjIOnbsyD7++GPufdq2bRsLCAhgjRo1YhYWFqxp06bsrbfeUvn/2ZDJGKuwcooQYpDGjx+PH374AQUFBVJ3hRAANNuKEEKIFih4EEIIEYxOWxFCCBHMaI481q9fD29vb9SpUwd16tSBv78/fv31V6m7RQghNZLRHHns378fZmZm8PDwAABs27YNn332GVJTU5Uy0BJCCBGX0QQPVRwdHfHZZ59h4sSJUneFEEJqFKNcJFhaWorvv/8ehYWF8Pf3V1uuuLi4yhQFZWVlePToEerXr682hQYhhJgixhjy8/PRtGlThfvH8H2y0UhLS2O2trbMzMyMOTg4KC1aq0y+SIh+6Id+6Id+1P/cuXNH8OexUZ22KikpQWZmJp48eYIff/wRX3/9NeLj49Xe65vPkUdubi6aN2+OxMREjdlQ5flw+NxqUoyyQuosKCjA+fPnuTuV6bt9McqKMSYhZcUaP99xSf360/5nmvtfdnY2/Pz88OTJE8EphYzqtJWFhQV3wdzHxwdJSUlYvXo1Nm7cqLK8paWlypxGqjRu3FjhFquVyVf2VrXjiFVWSJ1PnjzB9evX0bRpU9StW1fv7YtRVowxCSkr1vj5jkvq15/2P9Pc/+S0OWVvNFN1VWGMCUq7TAghRDeM5shj3rx5GDBgAFxcXJCfn49du3bhxIkTau/5QAghRDxGEzz++ecfjBkzBvfv34eDgwO8vb3x22+/oU+fPlJ3jRBCahyjCR6bNm2SuguEEEL+Y9TXPAghhEiDggchhBDBjOa0FSGESC05ORn3799Hq1at4OfnJ3V3JEXB4z9FRUUa79L29OlT3nWJUVZInfJFQs+ePUPt2prfYjHaF6OsGGMSUlas8fMdl9SvP+1/z7B92zbs3beP2xYSEoKoqCi9tC/W+IuKiniXraxGBY/Y2FjExsYqbCstLZWoN4QQY3H27Fns3bcP0YOA4C5A7J9AaEwMAgMD4evrK3X3JFGjgkdwcDCCg4MVtuXl5cHBwQFWVla8VmQKWbUpRlk+5V68eAGgPD2BLusVUk7XZcUck5Cyuq5T6Lho/5Nm/Hfu3AFQHjhsLMp/hx4E7t69i4CAAL31Vdd1WllZ8a6vMrpgTgghVZCnRYr9E3haUv4bAFq1aiVhr6RVo448CCFEGz4+Phg2dChC9+1D6MHybeFhYTX6ojkFD0II4WHsuHGYPmMGzbb6DwUPQgjhycfHp8qsujUFXfMghBAiGAUPQgghglHwIIQQIhgFD0IIIYLRBfP/UHoS3bVP6UkoPYn8N+1/ummf0pNIjNKTEEKIbtSo4EHpSfTTvhhlKT2JOO0LLUv7n2ntf5SehBBCiF5R8CCEECIYBQ9CCCGCUfAghBAiGAUPQgghglHwIIQQIhgFD0IIIYJR8CCEECJYjVokqAmlJ9Fd+5QegtKTyH/T/qeb9ik9icQoPQkhhOhGjQoelJ5EP+2LUdZU00NQehLdty9GWVPd/yg9CSGEEL2i4EEIIUQwCh6EEEIEo+BBCCFEMAoehBBCBKtRs60IIeJKTExEamoq8vPz0b17d6m7Q0REwYMQohPhYWFYHh3N/X3v3j2sWrVKug7pCAVE1Sh4EEK0kpiYiGvXrqFVq1YAgOXR0YgeBAR3AWL/BEJXr8aoUaPg5+cncU+1Z6oBURfomgchRLDwsDB06tQJY8eORadOnbBo0SIA5YHDxqL8NwBcu3ZNwl5WT2JiIhcQCz8FogcBq1evRmJiotRdMwh05PEfym2lu/Ypt5Bp57ZKSkpSPso4eBBA+b/l2wDA2dlZJ/+vpBh/WloaAMWAGHqwfLuXl5de+0q5rSRGua0Iqb709HQAyh+q/fv1Q+jBwwgtjyOYOXMmfH19Jexp9Xh4eABQDojy7TVdjQoelNtKP+2LUdZUcwsZY24rb29vAMofqgsWLsSChQu5i8uTJ0826v0vICAAYaGhCI2O5gLinDlzEBAQoJf29VFndXJb1ajgQQipPj8/P6UP1fCwMO7CuKenJ+Lj4yXsoe4sW74cQ4cNUwiIpBwFD0KIYPIPVflsK2OeUVUVPz8/kwqIukLBgxCiFT8/P5MOGkQzmqpLCCFEMAoehBBCBDOa4LF06VL4+vrC3t4eTk5OCAoKwtWrV6XuFiGE1EhGEzzi4+MRHByMhIQExMXF4cWLF+jbty8KCwul7hohhNQ4RnPB/LffflP4e8uWLXByckJKSgq6desmUa8IIaRmMprgUVlubi4AwNHRUW2Z4uJiFBcXa6wnLy8PAFBQUIAnT56oLSdf8i9f1KWJGGWF1Jmfn6/wW9/ti1FWjDEJKSvW+PmOS+rXn/Y/09z/NKWOqYqMMca0frZEGGMYMmQIHj9+jJMnT6ott2DBAixcuJBXnd9++y1sbGx01UVCCDF4T58+xejRo5Gbm4s6deoIeq5RBo/g4GAcPHgQp06dgrOzs9pyfI88XFxccPHiRTRt2lRtOXk05xNgxCgrpM78/HycPXsWHTp0gL29vd7bF6OsGGMSUlas8fMdl9SvP+1/prn/ZWVlwcvLS6vgYXSnrWbNmoVffvkFf/zxh8bAAQCWlpawtLTkVa+dnR3q1q2r9nF5xlM++WLEKCukTjl7e3uNYxKrfbHKArodk5CyYo4JqHpcUr/+tP+VM7X9T37aXhtGEzwYY5g1axb27duHEydOwN3dXeouEUJIjWU0wSM4OBjffvstfv75Z9jb2yM7OxsA4ODgAGtra4l7RwghNYvRrPNYv349cnNz0aNHDzRp0oT72b17t9RdI4SQGsdojjyM8Lo+IYSYLKM58iCEEGI4KHgQQggRjIIHIYQQwYzmmofYioqKNC7Vly+84UOMskLqfPbsGfdbPudbn+2LUVaMMQkpK9b4+Y5L6tef9j/T3P+Kiop4l62sRgWP2NhYxMbGKmwrLS2VqDeEEGK8eAWPX375RXDFffr0Mbj1F8HBwQgODlbYlpeXBwcHB1hZWfFakSlkha0YZfmUkydEs7a21mm9QsrpuqyYYxJSVtd1Ch2XFK9/YmIi0tLS4OHhgYCAgCrrov1P9+2LVaeVlRXv+irjFTyCgoIEVSqTyXD9+nW0aNFCmz4RQgxEeFgYlkdHc3+HhYZi2fLlEvaIGAreF8yzs7NRVlbG64ey0xJi/BITE7E8OhrRg4DCT4HoQcDy6GgkJiZK3TViAHgFj3Hjxgk6BfXuu+8KztBICDEs165dAwAEdwFsLMp/V9xOajZewWPLli1VpiGuaP369WjQoIHWnSKESK9Vq1YAgNg/gacl5b8rbic1G63zIISo5Ofnh7DQUIQeBGznAaEHy6+B+Pn5Sd01YgAET9UtLCzEsmXLcPToUTx48ABlZWUKj9+8eVNnnSOESGvZ8uUYOmyYoNlWpGYQHDwmTZqE+Ph4jBkzBk2aNIFMJhOjX4QQA+Hn5wcvLy+pu0EMjODg8euvv+LgwYPo0qWLGP0hhBBiBAQHj3r16sHR0VGMvkiK0pPorn1KD0HpSeS/jWn/S0pKQnp6Ojw8PODr66tQzlT3P72mJ1m0aBEiIiKwbds2o1vPQelJCCGqREREICYmhvs7JCQEUVFREvbI8AkOHp9//jlu3LiBRo0awc3NDebm5gqPnz17Vmed0zVKT6Kf9sUoa6rpIYwhPYnQcsa2/yUlJSEmJgbRg8rXssT+CYTGxMDLywsWFhZo1aoVPD09AZje/id6epKKhKYqIYQQQ5aeng5AcTFk6MHyyUFyffv2Re/evdG9e3epumlwBAePyMhIMfpBCCGS8PDwAFB+xCE/8gCAmZ0Bi9rAyj+AI0eO4MiRI7h37x5WrVolXWcNiNaLBFNSUrBz50588803SE1N1WWfCCFEb3x9fZUWQwLAm+3KA0fF3F6rV6+m3F7/EXzk8eDBA7z99ts4ceIE6tatC8YYcnNzERAQgF27dqFhw4Zi9JMQokFiYiKuXbuGVq1a0ZoMLcgXQ167dg0lJSWYNGkS1p8uf6zy6axr167RKntoceQxa9Ys5OXl4eLFi3j06BEeP36Mv//+G3l5eZg9e7YYfSSEaBAeFoZOnTph7Nix6NSpEyIiIpTKJCYmYseOHfStWQM/Pz+MGTMGEydORFhoKHafL99Oub1UExw8fvvtN6xfvx6tW7fmtrVp0waxsbH49ddfddo5QohmqtKmx8TEICkpiStTObiEh4VJ2GPjsGz5ciQkJGDQoEEKp7PmzJlDRx3/EXzaqqysTGl6LgCYm5sr5bkiRCpJSUm4e/cuWrVqZdL/2VWlTQ89WD6DKCAgQCG4cNNQo6MxdNgwo3hd5KfjnJ2dlRbuic3Pzw8HDhxAYmIiUlNTkZ+fj8mTJ+u1D4ZM8JFHz5498f777yMrK4vbdu/ePYSEhKBXr1467Rwh2oiIiEDPnj1rxDdtdWnT5TOIjPmeHBWPmHr27KnydJw++Pn54e2336bTVZUIPvJYu3YthgwZAjc3N7i4uEAmkyEzMxNt27bFzp07xeijXlB6Et21L2V6EpULvqKj0a9/f5XfXKUef3XTk3h5eSEkJAShMTHcLKFZs2bBy8sLBQUFcHZ2BqA8DdXZ2RkFBQUGu/8lJSUpHzHFxCAwMLDKIxBKj2Og6UlcXFxw9uxZxMXF4cqVK2CMoU2bNujdu7fWndAXSk9i+tQt+EpPT9f7aQ99iYqKQmBgIJeXqeJsK19fX6XgMnfuXIN/LWri+2hsBAcPuT59+qBPnz667IvoKD2JftoXoyzfMXl7ewNQ/qbt7e2t8XnGnp4kICCAu9eG/AhaXnblypUYOXIkN5VX1bUOQ9v/tH0fhfaVb1lTTY8jenqSNWvWYMqUKbCyssKaNWs0lqXpukRKfn5+St+06e535a+LMb0G3F0Mo6MVjpiMaQymjlfwiImJwTvvvAMrKyuFzJOVyWQyCh5EcvLTODVhtpUpq7hwT4rZVkQzXsEjIyND5b8JMVS+vr50y1QTID9i0jSZhUhD8FTdqKgolVfznz17RvnvCSGkhhAcPBYuXKjyW8DTp0+xcOFCnXSKEGI6KDWKaRIcPBhjkMlkStvPnz9vkrenJYRoj1KjmC7eU3Xr1asHmUwGmUyGVq1aKQSQ0tJSFBQUYNq0aaJ0khBifJKTk6udGkXK9CREM97BY9WqVWCMYcKECVi4cCEcHBy4xywsLODm5gZ/f39ROkkIMT7qFvrxTWkeHhaG5dHR3N8hISFYuXKlaP0lwvAOHuPGjQMAuLu7o0uXLlUu0SeE1Gzq7tDHJ0eUyoSOMTEYOXIkTb02EIIjQGFhIY4ePYp+/fopbD98+DDKysowYMAAnXVOnyi3le7ap9xC+sttpa+y2oxJVd6tuXPncnm3NNWblpYGQPmoJS0trcqbXdH+Z6C5rcLDw7Fs2TKl7YwxhIeHG3TwoNxWhOhX5bxbfK9bqDtqkW8n0hMcPK5fv442bdoobX/55Ze5c5yGinJb6ad9Mcqaam4hXeW20ldZbcZUMe8W33oDAgJUpicRsvCT9j8DyG1VkYODA27evAk3NzeF7enp6bC1tdW6I4QQUhGlJzFsgtd5vPHGG5gzZw5u3LjBbUtPT8cHH3yAN954Q6edI4TUbPL7ilPgMDyCg8dnn30GW1tbvPzyy3B3d4e7uztat26N+vXrY8WKFWL0kRBCiIHR6rTV6dOnERcXh/Pnz8Pa2hre3t7o1q2bGP0jhBBigLRarCGTydC3b1/07dtX1/0hhBBiBLQKHoWFhYiPj0dmZiZKSkoUHhPzfh5//PEHPvvsM6SkpOD+/fvYt28fgoKCRGuPEGL6kpKS6N4vWhAcPFJTUzFw4EA8ffoUhYWFcHR0xMOHD2FjYwMnJydRg0dhYSHatWuH9957D8OHDxetHUJIzRAREaFwg7uw0FAsW75cwh4ZD8EXzENCQhAYGIhHjx7B2toaCQkJuH37Njp27Cj6BfMBAwZg8eLFGDZsmKjtEEJMX2JiImJiYhA9CCj8FIgeBCyPjqbU8TwJPvI4d+4cNm7cCDMzM5iZmaG4uBgtWrRAdHQ0xo0bZ1Af7MXFxSguLtZYJi8vDwBQUFCAJ0+eqC0nX/IvXyykiRhlhdSZn5+v8Fvf7YtRVowxCSkr1vj5jkvq198U97/U1FQAyilQUlNT4enpqVDWVPe/6tyhUXDwMDc359KxN2rUCJmZmWjdujUcHByQmZmpdUfEsHTpUt43qDp//jyuX78uco/06+zZs1J3QedMcUyA8HFdu3YN9+7dQ7NmzXglGpSCob9X8kBQOQVKfn4+4uPjVT7H0McklJA8WJUJDh6vvvoqkpOT0apVKwQEBCAiIgIPHz7Ejh070LZtW607IoaPPvoIc+fO1VgmLy8PLi4uaNeuHZo2baq2nPxFtrGxqbJdMcoKqTM/Px9nz55Fhw4dYG9vr/f2xSgrxpiElBVr/HzHVbHOBQsWYPXq1dxj77//PhYsWCBqX01x/+vevTtu37qF0HXruBQoc+bMweTJk5XKmur+l5WVVWUZdQQHj08//ZSL2IsWLcK4ceMwffp0eHh4YMuWLVp3RAyWlpawtLTkVdbOzg5169ZV+7g8kyaffDFilBVSp5y9vb3GMYnVvlhlAd2OSUhZMccEVD0ueZ0XL17E6tWrFVOVr16NUaNGcTOFpB6/nDHsf8uWL8fwN9/kPdvK1PY/+Wl7bQgKHowxNGzYkEuJ3LBhQxw6dEjrxgkhwly7dg2A9jdYIsp8fX0FJVwk5QTNtmKM4aWXXsLdu3fF6o9GBQUFOHfuHM6dOwcAyMjIwLlz5wzuWgupuRITE/Hdd98hKSlJlPrl1zdi/wSelgi7wRIhuiQoeNSqVQsvvfQScnJyxOqPRsnJyXj11Vfx6quvAihP0fzqq68iIiJCkv4QUlF4WBg6deqEKVOmoGfPnggPC9N5G35+fuWpyg8CtvPKjzrCw8I0HnUkJiZix44dNAWV6JTgax7R0dH48MMPsX79erzyyiti9EmtHj16gDGm1zYJ4UPlbVOjozF02DCdn06qmKq8qvP0le8DLl8El5iYWK1U5/Ln06rsmktw8Hj33Xfx9OlTtGvXDhYWFrC2tlZ4/NGjRzrrHCH6VJ0PRH1fi/Dz86uyXnUBLSsrCzt27uTKhYSEYOXKlbzbVheQDEVSUhLS09Ph7e1NgU1EgoPHqlWrROgGIdKq7gdixWsRFdcMSHktQl1A27Fzp2JAiYnByJEjeX3Q6vMISxuGHthMCa/gMXfuXCxatAi2trZwd3dH586dq7wJPCHGIikpqdofiNy1iAq3Ta3qWoTY1AU0QPsjJEOe7WXogc3U8IoAX3zxBcLCwmBra4uAgADcv38fTk5OYvdNr4qKijQu1ReyElOMskLqfPbsGfe7qiAvRvtilBVjTPKyly5dAqD8gZiWlsZNS+dT5yfz56Nf//64dOkSWrRoga5du1aZ/oHvuLR5Tb28vBASEoLQmBguoI0aNQrfffedUkBxdnbmtf87OzsDUA5IFZ8v1f6XlpYGQPP7qG37Yu5/uiwntGxRURHvspXxCh5ubm5Ys2YN+vbtC8YY/vrrL9SrV09lWUO+KVRsbCxiY2MVtpWWlkrUG2IoWrRoAUD5A9HDw0NwXb6+vmo/qKQQFRWFwMBApKenw8PDA76+vmjcuLFCQJk1axbvi+a+vr5KAWnu3LkGcZtY+ftV+X08duwYRo0aJWHPTBTjYd++faxRo0ZMJpOxWrVqMZlMpvKnVq1afKozKLm5uQwAu337tsZy+fn5LD8/n1edYpQVUufjx4/ZTz/9xB4/fixa+wkJCWz79u0sISGhWn3lW1aMMVUsGxYaygBwP+FhYdWukw++49J1+/L379ixYyw/P1/j+6mqTk3l9bH/qTNmzBiF9zGgZflvvuNSR+z9T4o6GWPs9u3bDADLzc3lVb4iXkceQUFBCAoKQkFBAerUqYOrV6+a3Gkrwp8pXpQUMv1VH8SeCiufrVVQUKDVPS34zPaSQp8+fbBjxw58NQJo2xho26R8PYwhXJMxNYIWCdrZ2eH48eNwd3eHg4ODyh9i2ipelDS1eyD4+flhzJgxkn/IyBcbjh07Fp06dapysWF1FgEmJSWZ1D0t5JMEHj8tDxyGMOvNVPEKHhWTZ3Xv3r3KC0Z8ct4T46Rqtk3F7aR61AVndelOhAaaytLT0wHo5/3Ux0p3Pz+/8msyAlbgE+3wCh716tXDgwcPeFfarFkz3Lx5U+tOEcNFuZXEpS44yz/kKxIaaFSpeJFZzPezukFOiKioKBw7dgzbt29HQkICli5bJlpbNRmvax6MMXz99de80zE/f/68Wp0ihssQ1zPog75WLatbm6Fq5pe6NRfp6enVmj2l6/dT0/oLsWamUaZc8fEKHs2bN8dXX33Fu9LGjRvD3Nxc604Rw2ZoF5fFps8JAuqCs6pgICTQaBIVFYWRI0eK9n5qWlhoSNOaiTC8gsetW7dE7gYxNoY620bXpFi1rCo4q1rAJyTQVEXM99MQU7eQ6qt2jpEXL16gqKhI0B3GCDEWUqXj4PthzjfQSEnTqU5D6yvhj3fwOHToEHJycjBmzBhu25IlS7Bo0SK8ePECPXv2xO7du9WuPDd0lJ5Ed+0bW3oSTfik46hO+7pIT+Ll5cWd/ikoKDDI/U+euqXiSnchfa2p+5+YdQLVS0/Ce53HihUrFKbsnj59GhEREZg/fz727NmDO3fuYNGiRVp3RB9iY2PRpk0bhR9DSKtADBd3QbnC1E9DScdhbHx9fTFq1Ch67UwE7yOPv//+G59//jn39w8//IA+ffrg448/BgBYWVnh/fffF3RfAH0LDg5GcHCwwra8vDw4ODjAysqK16k3IafnxCjLp9yLFy8AANbW1jqtV0g5XZcVc0xVlV25ciWXH0rIbCsx3ittx1TVinXa/wx3/xOzTisrK971VcY7eOTn56N+/frc36dOncKbb77J/e3l5YWsrCytO0KIIfP19YWvr6/W1/akvPOeKaaTIdLjfdqqadOmuHz5MoDy86rnz59Hly5duMdzcnJgY2Oj+x4SYuT0uUCuMm3TydB9z0lVeAePN998E3PmzMGOHTswefJkNG7cGJ06deIeT05OhqenpyidJMRYSZ0LTJt0MlIGO2I8eAePyMhI+Pj4YPbs2Th37hx27twJMzMz7vHvvvsOgYGBonSSEH0Q49u21LnAhKaT0UXKE1Iz8L7mYWNjgx07dqh9/Pjx4zrpECFSEOu6gNQL5ISmk9FFyhNSMwhKyU6IKap4D3Ndn1riPrxFzvKq6ahp2fLlSEhI4JUoUN2RijZ3VSSmrdorzAkxdqrSkutyFbnYucAqHzWFhIQgKipKoQzfFeu6THlCTBsFD1Ljqbv3tS5PLYmVO0pl7q2YGAQGBvLKKqsqW7AxpDwh0qPg8R9KT6K79o0tPYSXl5dSWvK5c+fCy8uL2yfEGn9105OkpaUBUD5qunTpUpVHC5VvP1vxiKU6KU9o/6sZ6UlqVPCIjY1FbGyswrbS0lKJekMMSVRUFLeKXJ57yRioO2pq0aKFxudt374dMTExmNkZWD5Y8YjFWMZOpMU7eBw7dgwzZ85EQkIC6tSpo/BYbm4uOnfujA0bNqBr164676SuUHoS/bQvRll9pIcICAio8lSPrsdf3fQkAQEBStco5s6di65du6qtr+I1krWnAVtLIKJP+RHL3bt3Nb4GtP9RehI53sFj1apVmDx5slLgAAAHBwdMnToVK1euNOjgQYgpqnyNQtMNllReIzkIFBaXP17xOk/FlCq6ummTlGlaiG7xnqp7/vx59O/fX+3jffv2RUpKik46RYgxMYRUHn5+fhgzZkyVH8jqFi2uPa04hbjyKvOIiIhq95FWrpsW3sHjn3/+0Xhr2dq1a+Pff//VSacIMRbG9oGobh3H119/za3/ULXKPCYmplqrzGnluunhHTyaNWuGCxcuqH08LS0NTZo00UmnCDEGxviBqGrR4ty5czFx4kSujLqjE/l6GG2IUachMYSjT33jHTwGDhyIiIgIlVO7nj17hsjISAwePFinnSPEkBnrB2LFFefHjh3DwoULFR4XY5V5VXUa84evGKf4jAHv4PHJJ5/g0aNHaNWqFaKjo/Hzzz/jl19+wfLly+Hp6YlHjx5xN4YixJTJP+hKSkoAGF8qj4oXrVVNy1V3dFKdKbzq0rT4+voiIiLCqE79VSTGKT5jwXu2VaNGjXD69GlMnz4dH330ERhjAACZTIZ+/fph3bp1aNSokWgdJcQQVE4F4uf3GkIPnjGaVB58UpkAwmZw8aVq5frx48cRExOjOPsrOhpDhw0zitlYmhJJ8lnhb8wELRJ0dXXFoUOH8PjxY6Snp4Mxhpdeegn16tUTq3+EGIyKCRT/f5rrGXz99dewsLAQNZWHqjQiQglNZVIxpYquxlQ5TYvYecXEpi5rsqEffeqCVll1ZTIZZDIZzMzMIJPJdN0nQgySqg86ALCwsOA1TVZb4WFh6NmzJ6ZMmVKt0zqGeI2m4gp5PvcbMTRinOIzFoKOPG7duoXg4GAcPnxY4bRV//79sXbtWri5uYnRR72g3Fa6a99Ucws5OzsDUP6W6ezsrLDv6PK9Unm0Ex2Nfv37q/2AUtc+3/4LqVMVIe+Vl5cXZs6cidC1a9XmFRPavr73v0/mz0e//v251DZeXl46/39l1Lmt7ty5g06dOsHc3ByLFi1C69atwRjD5cuXsX79evj7+yMpKYnbQQ0R5bYi1dGxY0eVCRTF/Jap7rSONjdn8vX1Ver/rFmz0LFjR113W5D58+dj2LBhRpdXrCJfX1+u30I+vI0Z7+ARGRkJT09PHD58WCEfytChQxESEoL+/fsjMjISmzZtEqWjukC5rfTTvhhlDSW30MqVKzFy5EheKTZ08V55e3sDUD5a8Pb2rrJ+VY9X7r/8QrjU+x+fvGJ82xejrKHsf0aZ2+q3337Dnj17VDZmbW2NRYsW4e2339a6I4QYC7HuzaGuLSG3keVbp64vhJOah3fwyMnJ0XhNo0WLFsjJydFFnwghFSxbvpw7p16d2VaE6BLv4NG0aVNcvHhR7TWNv//+m9KTECIS+Tl1IactiHCU9Zc/3lN1hwwZgg8//FBl8sMHDx4gLCwMQUFBuuwbIYTojbEluZSaoAvmhw4dQsuWLfHuu+/i5ZdfBlB+u8tvv/0WjRs3rjE5XQghpkXlAkojWukuBd5HHvXq1UNiYiLeeecd7Nq1C3PmzMGcOXOwZ88ejB49Gn/99RccHR3F7CsAYN26dXB3d4eVlRU6duyIkydPit4mIcbMmJMO6ou6BZTy7UQZ7+CRmZmJunXrYv369cjJyUF2djays7ORk5ODDRs2oH79+mL2EwCwe/duzJkzBx9//DFSU1PRtWtXDBgwAJmZmaK3TYguyT/Qk5OTRW2HTsXwoy7rr7GsdJcC7+Dh7u7OXe+QyWRwcnKCk5OTXtOTrFy5EhMnTsSkSZPQunVrrFq1Ci4uLli/fr3e+kBIdVX8QO/Tpw+2b9smSjvq7jdCRyDK1GX9pVNW6vG+5iFPRyKVkpISpKSkIDw8XGF73759cfr0aZXPKS4uRnFxscZ68/LyAJTPd3/y5InacvJVo/LFQpqIUVZInfn5+Qq/9d2+GGXFGJOQsrqqMzk5Wfnc+r59GDd+PLp161ZlnQkJCdxKbB8fH41lU1NTASivTk9NTYWnp6egcdWE/S/8o4/Qu08fhddX/plgKvtfZdVZ5yMot5WUHj58iNLSUqW0740aNUJ2drbK5yxdulTpRjfqnD9/HtevX692Pw3J2bNnpe6Czhn7mI4fPw5A+QP922+/xbFjx9CsWTO1p0q2b9uGvfv2cX8PGzoUY8eNU9uW/IOu8ur0/Px8xMfH62hE6hnre9WkSRMUFhaqfI2MdUzqVCeViqDg8fXXX1c5z3z27Nlad4aPyqfJGGNqT5199NFHmDt3rsb68vLy4OLignbt2qFp06Zqy8lfZBsbmyr7KEZZIXXm5+fj7Nmz6NChA+zt7fXevhhlxRiTkLK6qtPW1harV69W+kDfvXs3V+b999/HggULFJ536tQp7N23T+mIZfqMGUpHIPL2u3fvjnv37iF09WpudfqcOXMwefJkweOi/c809r/KsrKyqiyjjqDgsWHDBpiZmal9XCaTiRY8GjRoADMzM6WjjAcPHqi9CZWlpSUsLS151W9nZ4e6deuqfVyeSZPPIi0xygqpU87e3l7jmMRqX6yygG7HJKSsvNzFixerXESmqc7evXsrpRsBoBgUVq/GqFGjFOq/f/8+AOUjlvv37yu9HhXbX7VqFUaNGqW2z7T/Gdf+p+sxyU/ba0NQ8EhOToaTk5PWjVWHhYUFOnbsiLi4OAwdOpTbHhcXhyFDhkjSJ1KzREREICYmhvs7LDQUy5YvF1xPxTvq/f3334iOjq7yZkgV73tR8YiFz2wgfebiIjUH7+BhCDd9mjt3LsaMGQMfHx/4+/vjyy+/RGZmJqZNmyZ114iJS0pKUnm71JdatVK4iyBf8g/033//HdHR0VUGBVXp1Gk2EJGS0cy2AoCRI0ciJycHUVFRuH//Pl555RUcOnQIrq6uUneNmDh199WYNGkSV0abIxEfHx8MGzoUofv2VRkUoqKiFNKpA8COHTsoDxORhKD0JIaQlG3GjBmYMWOG1N0gNUzl00Zd/7unmKp0FvJ7ZPA1dtw4TJ8xA/fv368yEMiPWMLDwrA8Oprbru0pNEK0xXuRYGRkJK+r94SYIu600X+LyM7eK9+uq3QWPj4+vO+DTov/iCHgHTwIqemioqKQkJCAqVOnctukSGehLg/Tli1bKIAQvaHgQch/+CQQ9PPzw3vvvQcA6NESCuksxvI8cqgudXmYNm7ciE6dOlF2a6IXRrPCXGxFRUUal+oLWYkpRlkhdT579oz7LZ/zrc/2xSgrxpgqlq08DTckJARRUVEq6/Ty8kJISIhC+VGjRiF23ToUFBSI8l6par/izKuAlsCBif9de4mJQZ8+fdC1a1eVdSUlJXEpOPhen6H9T9z9T1flhJYtKiriXbayGnXkERsbizZt2ij8+Pr6St0tIrGUlBRuGq78GkJMTAySkpLUPicqKgrHjh3Dl19+yf3WJ3n78qOgAxMVT2HdvHlT5fMiIiLQs2dPTJkyBT179sSiRYv01WViYngdebz66qu813kYcu6X4OBgBAcHK2zLy8uDg4MDrKyseM0mEzLjTIyyfMrJE6JZW1vrtF4h5XRdVswx3b17F4DyNNzdu3fDxsZG4VRUxToDAgIQEBBQrfaFjqty+zY2NtiyZYvSOpE2bdoo1ZeYmKi8VmXtWgwbNkzjOMQeE996hZTTdVkxxySkrK7rtLKy4l1fZbyOPIKCgjBkyBAMGTIE/fr1w40bN2BpaYkePXqgR48esLKywo0bN9CvXz+tO0KIVCpOw1V1DcGQ74GhKpX43LlzVR5Rq7vQLl/DQogQvI48IiMjuX9PmjQJs2fPVjrcjYyMxJ07d3TbO0L0wNfXVynflMI1hOho9OvfX++nOBMTE7kFgZquTVRMd6KpbMUL7RWPUuTBkxAhBF/z+P777zF27Fil7e+++y5+/PFHnXSKEH1btny5wjTcytcQ9P3tvPIdAKuaQeXn51flOhEhRymEVEVw8LC2tsapU6eUtp86dapa588IkVrFabiVT2Hp89u5qkWAVV3A50seJLdv346EhATe97shpDLBU3XnzJmD6dOnIyUlBZ06dQJQfnezzZs30/xyYvS4b+cVTmGFh4Xp9du5qmsToQfLj374XtjWpGKW3ercSY7UbIKDR3h4OFq0aIHVq1fj22+/BQC0bt0aW7duxVtvvaXzDhJSUXJyMq8cUNVR+RqCn5+fXj9k6doEMQZaLRJ86623KFAQvat8G1axkgFWvFBdMUAlJSXh7t27ggKXuro0UXX0Q9cmiKHRapHgkydP8PXXX2PevHl49OgRgPL1Hffu3dNp5wiRS05O5m7DKmYywMoXquXTdOWL6ypv15TSRF1dfBjCtQk+6VpIzSX4yCMtLQ29e/eGg4MDbt26hUmTJsHR0RH79u3D7du3sX37djH6KTpKT6K79sUoe/nyZQDK1wHS0tKUpqZq235SUhJ3obpimvXmrq4qbwR1OzMTu3bt4p4vT2ny9OlTpKSkqKyr8pRfTe+Vl5cXNzZ9v/6V07XMnDkT8+fP51WnKe5/lJ5EmeAjj7lz52L8+PG4fv26wuyqAQMG4I8//tC6I/pA6UmMV4sWLQCIOwtK1Q2fgPKjHlXbd+3apXZGlDw9iDEuyKt410T52NauXYuUlBSpu0YMiOAjj6SkJGzcuFFpe7NmzZCdna2TTomF0pPop30xynbp0kXlHfeqmx6kYllvb28Ayhequ3TpojIFCKB8JHT37l34+vqiTZs2Kuvy9vZW6Fd10pOIVVZdupa7d+/W2P2P0pMoE3zkYWVlhby8PKXtV69eRcOGDbXuCCFVGTtuHOLi4rjrAEuXLdNp/aoW0YWHhWHixIkKN4KSp18H1N/Pg1u1XqkuY7hdrLqU7zTbi1Qk+MhjyJAhiIqKwp49ewAAMpkMmZmZCA8Px/Dhw3XeQUIq8vHxQd26dUWrX9U0XaA8i21gYKDCbKsmTZoozIiS389Dfu1MXV2GjmZ7ET4EB48VK1Zg4MCBcHJywrNnz9C9e3dkZ2fD398fS5YsEaOPRALaTDE1FRUX0VXk6+urcJps2fLlyLp/Hzt27AAAbN+xA02aNMEnFS4sq6vL0PHNl0VqLsHBo06dOjh16hSOHTuGs2fPoqysDB06dEDv3r3F6B+RQHhYGJZHR3N/h4WGKnwgknLyqaxVzagyVrQSnWgiOHhkZmaiUaNG6NmzJ3r27MltZ4zhzp07aN68uU47SPSrYl4lfX4ginWkI+YRlKY0IqYQPAjRRPAFczc3N3To0AE3btxQ2P7gwQO4u7vrrGNEGlLc8yEiIkLrxXSaVGeRHh90YZnUZFqtMG/dujVee+01HD16VGE7Y0wnnSLS0fcHoqo1BbpYOV5xwZ9YK9LVzc6iow5SEwgOHjKZDOvWrcMnn3yCQYMGYc2aNQqPEWnoKpWEvj8Q1S3Mkx8BGVq9lVVOI6Lr6cOEGCrB1zzkRxchISF4+eWXMWrUKKSlpRl9OnZjTk9SOZVEcHAw+vTpo3UqhU/mz0e//v2Rnp4ODw8P+Pr6ijZ+Z2dnAMqL6ZydnRXeD6HpIfjWq4vXv2IakYKCAlFSyWiqMykpSfT3qqanx6H0JMq0yqorN2DAAJw+fRpvvPEGzpw5U52q9CI2NhaxsbEK20pLSyXqjW5UPO3DXeCOjYWrqytee+01rev19fUV9fSL/APP2dkZM2fOROjatTpdU9CxY8fyhX0xMQr1AsB3333HfdAau8pfHEJCQhAeHi5hj0hNITh4dO/eHRYWFtzfbdq0wZkzZzB06FCDv+ZhiulJ1KWSuHfvnsGmh6g8FTgkJAQJCQkaZ0Vpkx5i5cqVGDlyJFfvvr17FWYIVpyCLFV6iOqkJ0lMTFT+4hATg8DAQPj6+vLuq5BU85SehNKTyAm+5nH8+HGlFb6Ojo6Ij49HWVmZ1h0h2lF3gbtZs2YS9ko9dbdYBVDlPbi1Ib+3NwCVF9B1cWtXqehiZpy6VPOEVIVX8KiYyyovL0/jD9EvVRe458yZwwUVQ6PuA2/Lli06v29ExUkEUkxBFlt1Z8ZVPHIR8x4pxDTxCh716tXDgwcPAAB169ZFvXr1lH7k24n+VZ7xExkZKXWX1Kr8gTdoU/n2jRs36vSbb+W1I3FxcQrtmsKajOrOjFMXUHU9I42YJl7XPI4dOwZHR0cA5aetiOGpmEriyZMn0nZGA1VJ9yqvZh86bFi1Tl+pnESwYwfGvPsuQnfuVEjpLuSiuTa3oRVbde63ru5e6YZ61EoMC6/g0b17d5X/JkQb8g+8LVu2YOPGjUoX+69du1atD2dVazxCDwJ9+vZF8MyZaj9oNaUyqTyrSaz7p2tD2+SLfn5+SjPSjCVtPJEer+CRlpbGu0L5DXVIzVLxg5dPBlb5B9TGjRt1/s1XfipKVb3qPmgrzwAb8+676NO3L9cXVbehre4RkiFQlWqeED54BY/27dtDJpOBMVblKnJjXzdBhFM19TYqKqrK54n1zdfX11dQvZXvXT5oE7Bj507s2LkTADBo0CAAykcy1T1CMhSVU80Twgev4JGRkcH9OzU1Ff/73//w4Ycfwt/fHwDw119/4fPPP0d0hQ8QUjOozML731oDPh9IYn3zjYqKUljjoaneiqe5LtwHTtyodB3mYHkEomsDhPw/XsHD1dWV+/eIESOwZs0aDBw4kNvm7e0NFxcXzJ8/H0FBQTrvJDFcmtKS8/02K9Y3X77XAiqe5qpnU76t8nj69euH0IOH6doAIf8RvML8woULKlOvu7u749KlSzrplBSMObdVZfrMLcQ3h1R12xczt5CXlxd3mkuu8njef/99vP/++7h79y6X2kTd+PSd20ofZQ11/9NXWcptpUxw8GjdujUWL16MTZs2cUvbi4uLsXjxYrRu3VrrjuiDKea2kpqq6wuzZs1Cx44dpe2YQPLTZ+np6Th27BhCd+1SyIklH0/Xrl0l7CV/KSkpCoGOEF0THDw2bNiAwMBAuLi4oF27dgCA8+fPQyaT4cCBAzrvoC6ZYm4rVfSdW6hyDin5bCtjyy0UEBCAgIAATJ48GXPmzFE5pdeQclupo820YmPe//RRlnJbKRMcPF577TVkZGRg586duHLlChhjGDlyJEaPHg1bW1utO0KMm6nd71rbtRNSU5ks0USmFRPDIih4PH/+HJ6enjhw4ACmTJkiVp8IIVpSN4HBVKYVE8MhKKuuubk5iouL6Y6BhBgodckSaVox0TXBKdlnzZqF5cuXc+cACSGGg1t4WSlZIh11EF0TfM0jMTERR48exZEjR9C2bVul6xx79+7VWecIIcJRyhGiD4KDR926dTF8+HAx+qLRkiVLcPDgQZw7dw4WFhYGnTmWGB9NSRGNEaUcIWITHDy2bNkiRj+qVFJSghEjRsDf3x+bNm2SpA/ENBlyxlxCDJXgax5SWbhwIUJCQtC2bVupu0JMSMV7f9Dd9AjhT/CRBwD88MMP2LNnDzIzM1FSUqLw2NmzZ3XSMV0oLi5GcXGxxjLyW+cWFBRoPBUmX/LPZ6KAGGWF1Jmfn6/wW9/ti1FWjDEB5el2AOWprampqfD09NSqTjHeK6lff9r/xNn/pB5/ddZkCQ4ea9aswccff4xx48bh559/xnvvvYcbN24gKSlJafW21JYuXYqFCxfyKnv+/Hlcv35d5B7plyEFcl3R9ZjkXy7CDgDtmgHn75Vvz8/PR3x8vE7b0oTeK+NgamMSkgerMhljjAl5wssvv4zIyEiMGjUK9vb2OH/+PFq0aIGIiAg8evQIa9eu5V3XggULqvxwT0pKgo+PD/f31q1bMWfOHF4XzPkeebi4uODixYto2rSp2nLyF9nGxqbKdsUoK6TO/Px8nD17Fh06dIC9vb3e2xejLN8xJScn49KlS2jRogVef/11Xu0HDh6Ms6mp3DYfn46Ii/tdq34KLct3XFK//rT/6X5MQsqKNf6srCx4eXkhNzcXderUqbJ8RYKPPDIzM9G5c2cA5Xle5IdxY8aMQadOnQQFj5kzZ+Ltt9/WWMbNzU1oFzmWlpawtLTkVdbOzg5169ZV+7g8kyaffDFilBVSp5y9vb3CmFTNKBKjfbHKAspjqqjyTan4XPhOTU3F2dTUSvfvSMHVq1cVZl2JOSZA87jEbF+f+5++2pdq/xOzfbHGJD9trw3BwaNx48bIycmBq6srXF1dkZCQgHbt2iEjIwMCD2LQoEEDNGjQQGgXiBaEfrAa29RVlTel4pHTSd39zimdByGaCZ5t1bNnT+zfvx8AMHHiRISEhKBPnz4YOXIkhg4dqvMOymVmZuLcuXPIzMxEaWkpzp07h3PnzplEEj6xVfxg5TOjKDwsDJ06dcLYsWPRqVMnhIeF6bnHwqnK6VRxuzoVbwRF6TwI4U/wkceXX36JsrIyAMC0adPg6OiIU6dOITAwENOmTdN5B+UiIiKwbds27u9XX30VAHD8+HH06NFDtHZNgaZkefL06XLqvsH3699f7X0hDOEopWJOp+Au5RfAASjNBqxM6P3OCSHlBB951KpVS+FOWm+99RbWrFmD2bNnw8LCQqedq2jr1q1gjCn9UOCompBkeeq+wctP71Qm1lFKYmIiduzYwXu9hZ+fH8JCQ7mcTmtPl2+fNGlSlX2KiopCQkICtm/fjoSEBCxdtqy63SfE5PE68khLS+Ndobe3t9adIeLgPlijo5W+XVc+7Vf5G7w80MhP71SkzVEKH6quz4R/9FGVz1u2fDleatUKkyZNUuqTlbU1BgwYoPaIwljv30GIVHgFj/bt20Mmk4ExVmU6drqtq2Fatnw5hg4bVuXpJXWBRlUwUHc6LD09XevgoS4g9e7Th9fz5Ue/8j5l/7ema+HChVi4cCGlHiFER3idtsrIyMDNmzeRkZGBH3/8Ee7u7li3bh1SU1ORmpqKdevWoWXLlvjxxx/F7i+pBj8/P4wZM6bKb9jLli/ndRpH3ekwVUcpfFV12iw5OVnj6ayKfYq/Aaz8A5R6hBAR8DrycHV15f49YsQIrFmzBgMHDuS2eXt7w8XFBfPnz0dQUJDOO6kPRUVFGmduCVmJKUZZIXU+e/aM+13x+pSQer28vLiL6QUFBSrLeXl5KV1snjt3Lry8vLQel7OzMwDl02YuLi7YvGkT9u7bx5UNCQlBVFSUUp9mzpyJ0ArrjSofGaWlpXFjE+P1F+O9qmn7X3XaF6OsGGMSUlas8RcVFfEuW5ng2VYXLlyAu7u70nZ3d3dcunRJ647oQ2xsLGJjYxW20Wm26pHfOyI9PR0eHh7w9fWtVsoDVbOf5s6dCwDYu2+f4umsmBgEBgYqnSKbP38+Bg4ciD/++APLli3jdf2GECKM4ODRunVrLF68GJs2bYKVlRWA8jQgixcvRuvWrXXeQV0KDg5Wyr+Vl5cHBwcHWFlZ8VqRKWSFrRhl+ZSTJ0SztrbWab3qygUEBKi8d4S241+5ciVGjhypcH1mw4YNAJSPIu7evauy7a5du2LAgAFgZWVK12+q01ddv6dC3yva/6QZv5hjElJW13XKP8O1ITh4bNiwAYGBgXBxcUG7du0AlCcVlMlkOHDggNYdIaSiyrOfKi7mq3gUUdViPr4TBfTBENbDEKIrgoPHa6+9hoyMDOzcuRNXrlwBYwwjR47E6NGjlW5JS4iu+Pj4YNjQoQjdt0/wYj5DmIarTd4tQgyZVvfzsLGxwZQpU3TdF0I0GjtuHKbPmIH79+8b1bd3TXm3Kt4zhBBjolXwuHbtGk6cOIEHDx5wqUrkIiIidNIxQlTx8fGpMqupodGUHoaCBzFWgoPHV199henTp6NBgwZo3LixwqJBmUxGwYMYjcTERKSlpcHDw0PlRXRdUbdqP+7IEQQGBorWLiFiEhw8Fi9ejCVLliDMCDKtEqKOPq9ByBdnhu7YwV2vCWgJ7Ni5E2PHjROlTULEJjgx4uPHjzFixAgx+kKIXghNUa8Lff5Lr/LVCCBhFnBgYvl2dQknCTF0goPHiBEjcOTIETH6QoheaHvvj+qQn7p6/BRo24QWLBLjJ/i0lYeHB+bPn4+EhAS0bdsW5ubmCo/Pnj1bZ53TJ0pPorv2DT09hLoUKM7Ozmr3geqOSVMqlzNnzlB6Eh22b+j7nzZlTSI9yZdffgk7OzvEx8cjPj5e4TGZTGbQwYPSkxi2pKQkhTQnYlGXAkXMNgHVqVzy8/NFbZMQsQgOHhkZGWL0Qy8oPYl+2udbtuJsp8O//abxArauxyRPgSJ0tlV1x185lQulJ9F9+2KUpfQkyrRa50FIdVWe7QRA5SI6MRcC+vn5Kd2GlxDCj1bB4+7du/jll1+QmZmpdI/olStX6qRjxHRVXnE9YTew+7zqRXTGsoqckJpGcPA4evQo3njjDbi7u+Pq1at45ZVXcOvWLTDG0KFDBzH6SExM5dlO0zuXBw+hSQ8JIdIRPFX3o48+wgcffIC///4bVlZW+PHHH3Hnzh10796d1n8QXirfgfBMZvn20IOA7bzy33yTHhJCpCH4yOPy5cv47rvvyp9cuzaePXsGOzs7REVFYciQIZg+fbrOO0lMi7r7pAcNHarTlOWUAp0Q8QgOHra2tiguLgYANG3aFDdu3OAuOj58+FC3vSMmS36fjcqznfh+yFcVGCIiIhATE8P9TSnQCdEtwcGjU6dO+PPPP9GmTRsMGjQIH3zwAS5cuIC9e/eiU6dOYvSRmChtZztVlZcqKSkJMTExep+9RUhNIviax8qVK7n/gAsWLECfPn2we/duuLq6YtOmTTrvICEVJScnV5mXSp4vSp/pRwipaQQfebRo0YL7t42NDdatW6fTDkmF0pPorn0x00NcvnwZgPK03rS0NO4oRkj6EanHz/e9kvr1p/2P0pNUJvjIo0WLFsjJyVHa/uTJE4XAYohiY2PRpk0bhR+xU1IQ3ZLvY/KZWqoSDHbs2BEzZ85UmL2lj/QjhNQkgo88bt26pTIfVHFxMe7du6eTTomF0pPop30xysrH1KVLF5UztSqnF1m6dCneffdd3rOtpBo/pSfRfftilKX0JMp4B49ffvmF+/fhw4fh4ODA/V1aWoqjR4/Czc1N644Qw5eUlIS7d+9KPvVVPlOrqsDg5+dHF8gJEQnv4BEUFASgPHPuuEp3PzM3N4ebmxs+//xznXaOGA5Dm/pKgYEQafG+5lFWVoaysjI0b94cDx484P4uKytDcXExrl69isGDB4vZVyKRxMREbuqrvu68RwgxbIIvmGdkZKBBgwZi9IUYKCnuvEcIMWy8g0diYiJ+/fVXhW3bt2+Hu7s7nJycMGXKFG7lOTEtlXNRUeJCQgjv4LFgwQKkpaVxf1+4cAETJ05E7969ER4ejv3792Pp0qWidJJIy8/Pr/zOexIlLkxMTMSuXbvoSIcQA8L7gvm5c+ewaNEi7u9du3bBz88PX331FQDAxcUFkZGRWLBggc47SfRDU74o+S1U9T3bqnIqknv37mHVqlV6aZsQoh7vI4/Hjx+jUaNG3N/x8fHo378/97evry/u3Lmj294RvQkPC0OnTp0wduxYdOrUCeFhYUplfH19MWbMGL0ecVRORbJ69Wq6UE+IAeAdPBo1asTdv7ykpARnz56Fv78/93h+fj7Mzc1130MiuqSkpCrzRUmBLtQTYrh4n7bq378/wsPDsXz5cvz000+wsbFB165ducfT0tLQsmVLUTqpDzU5t9WlS5cAaM4XJcX4heSoEqN9sesEKLeV1K8/5bbSPrcV7+CxePFiDBs2DN27d4ednR22bdsGCwsL7vHNmzejb9++WndEH2JjYxEbG6uwTVWqlZqmYr6oih/SFfNFScHX17f8Qn1MDJeKZObMmZSjihADwDt4NGzYECdPnkRubi7s7OxgZmam8Pj3338vKO+KFCi3lWpdu3bllS9KSD91VXblypUYOXIkUlNTkZ+fj8mTJ5tUbiHKbaX79sUoS7mtlAlOjFgxp1VFjo6OWneCSI9vvigp+Pn5wdPTE/Hx8VJ3hRDyH8HBg5guyhdFCOFLcHoSQgghhIIHIYQQwSh4EEIIEcwogsetW7cwceJEuLu7w9raGi1btkRkZCRKSkqk7hohhNRIRnHB/MqVKygrK8PGjRvh4eGBv//+G5MnT0ZhYSFWrFghdfcIIaTGMYrg0b9/f4U8Wi1atMDVq1exfv16Ch6EECIBowgequTm5la5tqS4uLjKe4zk5eUBAAoKCvDkyRO15eRL/uWLhTQRo6yQOvPz8xV+67t9McqKMSYhZcUaP99xSf360/5nmvufpjQ/VZExxpjWz5bIjRs30KFDB3z++eeYNGmS2nILFizAwoULedX57bffwsbGRlddJIQQg/f06VOMHj0aubm5qFOnjqDnSho8+Hy4JyUlwcfHh/s7KysL3bt3R/fu3fH1119rfC7fIw8XFxdcvHgRTZs2VVtOHs35BBgxygqpMz8/H2fPnkWHDh1gb2+v9/bFKCvGmISUFWv8fMcl9etP+59p7n9ZWVnw8vLSKnhIetpq5syZePvttzWWcXNz4/6dlZWFgIAA+Pv748svv6yyfktLS1haWvLqi52dHerWrav2cXkmTT75YsQoK6ROOXt7e41jEqt9scoCuh2TkLJijgmoelxSv/60/5Uztf1PftpeG5IGjwYNGqBBgwa8yt67dw8BAQHo2LEjtmzZglq1jGKWMSGEmCSjuGCelZWFHj16oHnz5lixYgX+/fdf7rHGjRtL2DNCCKmZjCJ4HDlyBOnp6UhPT+duECRnhNf7CSHE6BnFuZ/x48eDMabyhxBCiP4ZRfAghBBiWCh4EEIIEYyCByGEEMGM4oK5PhQVFWlcqi9feMOHGGWF1Pns2TPut3zOtz7bF6OsGGMSUlas8fMdl9SvP+1/prn/FRUV8S5bWY0KHrGxsYiNjVXYVlpaKlFvCCHEeNWo4BEcHIzg4GCFbXl5eXBwcICVlRWvFZlCVtiKUZZPOXlCNGtra53WK6ScrsuKOSYhZXVdp9Bx0f5H+58u67SysuJdX2V0zYMQQohgFDwIIYQIRsGDEEKIYBQ8CCGECEbBgxBCiGAUPAghhAhGwYMQQohgFDwIIYQIVqMWCWpC6Ul01z6lh6D0JPLftP/ppn1KTyIxSk9CCCG6UaOCB6Un0U/7YpQ11fQQlJ5E9+2LUdZU9z9KT0IIIUSvKHgQQggRjIIHIYQQwSh4EEIIEYyCByGEEMEoeBBCCBGMggchhBDBKHgQQggRrEYtEtSE0pPorn1KD0HpSeS/af/TTfuUnkRilJ6EEEJ0o0YFD0pPop/2xShrqukhKD2J7tsXo6yp7n+UnoQQQoheUfAghBAiGAUPQgghglHwIIQQIhgFD0IIIYJR8CCEECIYBQ9CCCGCUfAghBAiGAUPQgghgtWoFeaaUG4r3bVPuYUot5X8N+1/ummfcltJjHJbEUKIbtSo4EG5rfTTvhhlTTW3EOW20n37YpQ11f2PclsRQgjRKwoehBBCBKPgQQghRDAKHoQQQgSj4EEIIUQwowkeb7zxBpo3bw4rKys0adIEY8aMQVZWltTdIoSQGslogkdAQAD27NmDq1ev4scff8SNGzfw5ptvSt0tQgipkYxmnUdISAj3b1dXV4SHhyMoKAjPnz+Hubm5hD0jhJCax2iCR0WPHj3CN998g86dO2sMHMXFxSguLtZYV25uLgAgOztbYzl5egJra+sq+ydGWSF1FhQU4OnTp8jKykJeXp7e2xejrBhjElJWrPHzHZfUrz/tf6a5/8k/9xhjVZZVwoxIaGgos7GxYQBYp06d2MOHDzWWj4yMZADoh37oh37oR8PPjRs3BH8eyxjTJuToxoIFC7Bw4UKNZZKSkuDj4wMAePjwIR49eoTbt29j4cKFcHBwwIEDByCTyVQ+l8+Rx5MnT+Dq6orMzEw4ODhoLOvr64ukpCSNZcQsy7dcXl4eXFxccOfOHdSpU0fv7YtRVqwxCSkrRp1CxkX7H+1/uq4zNzcXzZs3x+PHj1G3bl1edctJetpq5syZePvttzWWcXNz4/7doEEDNGjQAK1atULr1q3h4uKChIQE+Pv7q3yupaUlLC0tefXFwcGhyp3CzMyM144jVlkhdQJAnTp1dFqv1OMHdD8mIWXFGhPAb1xSv/60/5nu/lerlvC5U5IGD3kw0Ib8gKmqIwtCCCG6ZxQXzM+cOYMzZ87g9ddfR7169XDz5k1ERESgZcuWao86CCGEiMco1nlYW1tj79696NWrFzw9PTFhwgS88soriI+P531aihBCiO4YxZFH27ZtcezYMam7oXQvEH2XFVKnEGK0L1ZZKds3xTEJKUv7H39Sv/5ivVcVSTrbyhDIbwaVm5sr6AKTIaMxGQ9THBeNyXhUZ1xGcdqKEEKIYaHgQQghRLAaHzwsLS0RGRlpUhfeaUzGwxTHRWMyHtUZV42/5kEIIUS4Gn/kQQghRDgKHoQQQgSj4EEIIUQwCh6EEEIEo+BRgandJ/3WrVuYOHEi3N3dYW1tjZYtWyIyMhIlJSVSd63alixZgs6dO8PGxkZwKmlDsW7dOri7u8PKygodO3bEyZMnpe5Stfzxxx8IDAxE06ZNIZPJ8NNPP0ndpWpbunQpfH19YW9vDycnJwQFBeHq1atSd6ta1q9fD29vby5DsL+/P3799VfB9VDwqMDU7pN+5coVlJWVYePGjbh48SJiYmKwYcMGzJs3T+quVVtJSQlGjBiB6dOnS90VrezevRtz5szBxx9/jNTUVHTt2hUDBgxAZmam1F3TWmFhIdq1a4e1a9dK3RWdiY+PR3BwMBISEhAXF4cXL16gb9++KCwslLprWnN2dsayZcuQnJyM5ORk9OzZE0OGDMHFixeFVST49lE1yM8//8xkMhkrKSmRuis6Ex0dzdzd3aXuhs5s2bKFOTg4SN0NwV577TU2bdo0hW0vv/wyCw8Pl6hHugWA7du3T+pu6NyDBw8YABYfHy91V3SqXr167Ouvvxb0HDryUIPvfdKNTW5uLhwdHaXuRo1WUlKClJQU9O3bV2F73759cfr0aYl6RfjIzc0FAJP5P1RaWopdu3ahsLBQ8O0tKHhUEhYWBltbW9SvXx+ZmZn4+eefpe6Szty4cQNffPEFpk2bJnVXarSHDx+itLQUjRo1UtjeqFEjZGdnS9QrUhXGGObOnYvXX38dr7zyitTdqZYLFy7Azs4OlpaWmDZtGvbt24c2bdoIqsPkg8eCBQsgk8k0/iQnJ3PlP/zwQ6SmpuLIkSMwMzPD2LFjubsWGgqhYwKArKws9O/fHyNGjMCkSZMk6rlm2ozLmMlkMoW/GWNK24jhmDlzJtLS0vDdd99J3ZVq8/T0xLlz55CQkIDp06dj3LhxuHTpkqA6jOJ+HtUh9n3SpSB0TFlZWQgICIC/vz++/PJLkXunPaHjMlYNGjSAmZmZ0lHGgwcPlI5GiGGYNWsWfvnlF/zxxx9wdnaWujvVZmFhAQ8PDwCAj48PkpKSsHr1amzcuJF3HSYfPEzxPulCxnTv3j0EBASgY8eO2LJli1Y3uteX6rxXxsTCwgIdO3ZEXFwchg4dym2Pi4vDkCFDJOwZqYwxhlmzZmHfvn04ceIE3N3dpe6SKBhjgj/nTD548GWK90nPyspCjx490Lx5c6xYsQL//vsv91jjxo0l7Fn1ZWZm4tGjR8jMzERpaSnOnTsHAPDw8ICdnZ20neNh7ty5GDNmDHx8fLgjwszMTKO+HlVQUID09HTu74yMDJw7dw6Ojo5o3ry5hD3TXnBwML799lv8/PPPsLe3544WHRwcYG1tLXHvtDNv3jwMGDAALi4uyM/Px65du3DixAn89ttvwirS+ZwvI5WWlsYCAgKYo6Mjs7S0ZG5ubmzatGns7t27UndNa1u2bGEAVP4Yu3Hjxqkc1/Hjx6XuGm+xsbHM1dWVWVhYsA4dOhj99M/jx4+rfE/GjRsndde0pu7/z5YtW6TumtYmTJjA7XcNGzZkvXr1YkeOHBFcD6VkJ4QQIpjhngAnhBBisCh4EEIIEYyCByGEEMEoeBBCCBGMggchhBDBKHgQQggRjIIHIYQQwSh4EJ24desWZDIZt9JbzLpPnDgBmUyGJ0+e6LwtuQULFqB9+/ai1S+W8ePHIygoSJK23dzcuASWYr432nJzc8OqVat0Wuf48eO5MZvCnROFoOBhpB48eICpU6eiefPmsLS0ROPGjdGvXz/89ddfXBlT3aE7d+6M+/fvw8HBQbI+yAOa/KdevXro1q0b4uPjJesTAKxevRpbt27l/u7RowfmzJlT7XoLCwsRFhaGFi1awMrKCg0bNkSPHj1w4MABhXJRUVEK74080Ff++eSTT6rdJ3W2bt2q8tbESUlJmDJlik7bWr16Ne7fv6/TOo0F5bYyUsOHD8fz58+xbds2tGjRAv/88w+OHj2KR48eSd01rZWUlMDCwqLKchYWFgaTm+v333+Hl5cXHjx4gHnz5mHgwIH4+++/tUqgx3f8mogVUKdNm4YzZ85g7dq1aNOmDXJycnD69Gnk5OQolLO3t1f53ly9ehV16tTh/laVf6y0tBQymUy05J0NGzbUeZ0ODg6SfomRlM4TpxDRPX78mAFgJ06cUFvG1dVVIRePq6srY4yx9PR09sYbbzAnJydma2vLfHx8WFxcnNJzlyxZwt577z1mZ2fHXFxc2MaNGxXKJCYmsvbt2zNLS0vWsWNHtnfvXgaApaamMsYYe/HiBZswYQJzc3NjVlZWrFWrVmzVqlUKdYwbN44NGTKEffrpp6xJkyZcH6uqW55D6fHjx4wxxrp3764y/1BGRgZjjLEnT56wyZMns4YNGzJ7e3sWEBDAzp07p9CXpUuXMicnJ2ZnZ8cmTJjAwsLCWLt27dS+vhkZGQp9Yoyxu3fvMgBsw4YNjDHGTpw4wXx9fZmFhQVr3LgxCwsLY8+fP+fKd+/enQUHB7OQkBBWv3591q1bN17P+/7779krr7zCrKysmKOjI+vVqxcrKChQeE3l/678mty8eZO1bNmSffbZZwrjuXDhApPJZCw9PV3leB0cHNjWrVvVvh6Mle83MTExCtsqv1cVyW8hvH//fta6dWtmZmbGbt68yc6cOcN69+7N6tevz+rUqcO6devGUlJSFJ77+PFjNnnyZObk5MQsLS2Zl5cX279/v8r8WpGRkSr7d/v2bfbGG28wW1tbZm9vz0aMGMGys7O5xyMjI1m7du3Y9u3bmaurK6tTpw4bOXIky8vLUxoLTPS2u5rQaSsjZGdnBzs7O/z0009q0ygnJSUBALZs2YL79+9zfxcUFGDgwIH4/fffkZqain79+iEwMBCZmZkKz//888/h4+OD1NRUzJgxA9OnT8eVK1cAlJ/CGDx4MDw9PZGSkoIFCxbgf//7n8Lzy8rK4OzsjD179uDSpUuIiIjAvHnzsGfPHoVyR48exeXLlxEXF4cDBw7wqruyvXv34v79+9zPsGHD4OnpiUaNGoExhkGDBiE7OxuHDh1CSkoKOnTogF69enFHaXv27EFkZCSWLFmC5ORkNGnSBOvWreP5bvw/GxsbAMDz589x7949DBw4EL6+vjh//jzWr1+PTZs2YfHixQrP2bZtG2rXro0///wTGzdurPJ59+/fx6hRozBhwgRcvnwZJ06cwLBhw1TesGz16tXw9/fH5MmTudemefPmmDBhArZs2aJQdvPmzejatStatmypcmyNGzfGoUOHkJ+fL/h10eTp06dYunQpvv76a1y8eBFOTk7Iz8/HuHHjcPLkSSQkJOCll17CwIEDubbLysowYMAAnD59Gjt37sSlS5ewbNkymJmZoXPnzli1ahXq1KnDjVnV/sMYQ1BQEB49eoT4+HjExcXhxo0bGDlypEK5Gzdu4KeffsKBAwdw4MABxMfHY9myZTp9DYyWxMGLaOmHH35g9erVY1ZWVqxz587so48+YufPn1coA57fhtq0acO++OIL7m9XV1f27rvvcn+XlZUxJycntn79esYYYxs3bmSOjo6ssLCQK7N+/Xqlb+KVzZgxgw0fPpz7e9y4caxRo0asuLiY28anbk3fZleuXMnq1q3Lrl69yhhj7OjRo6xOnTqsqKhIoVzLli25oyl/f382bdo0hcf9/PwEHXkUFBSwqVOnMjMzM5aWlsbmzZvHPD09WVlZGfec2NhYZmdnx0pLSxlj5Uce7du3V6i3quelpKQwAOzWrVsq+1XxyEPexvvvv69QJisri5mZmbHExETGGGMlJSWsYcOGGo8s4uPjmbOzMzM3N2c+Pj5szpw57NSpUwplNB152NraKvw8fPiQy/pc+SiwshcvXjB7e3u2f/9+xhhjhw8fZrVq1eLe48rkRzSVVezfkSNHmJmZGcvMzOQev3jxIgPAzpw5wxgrP/KwsbFRONL48MMPmZ+fn1LdfP+vmRI68jBSw4cPR1ZWFn755Rf069cPJ06cQIcOHRQulqpSWFiI0NBQtGnTBnXr1oWdnR2uXLmidOTh7e3N/Vsmk6Fx48Z48OABAODy5cto164d900bgMp7nmzYsAE+Pj5o2LAh7Ozs8NVXXym107ZtW4Xz/HzrVuXXX39FeHg4du/ejVatWgEAUlJSUFBQgPr163NHbHZ2dsjIyMCNGze4Niu3wbfNzp07w87ODvb29ti/fz+2bt2Ktm3bcnVWvK1sly5dUFBQgLt373LbfHx8FOqr6nnt2rVDr1690LZtW4wYMQJfffUVHj9+zKuvck2aNMGgQYOwefNmAMCBAwdQVFSEESNGqH1Ot27dcPPmTRw9ehTDhw/HxYsX0bVrVyxatIhXmydPnsS5c+e4n3r16gEov35VcV8DyieDTJs2Da1ateKuKRQUFHD7zrlz5+Ds7My9x9q4fPkyXFxc4OLiwm2T/5+4fPkyt83NzQ329vbc302aNOH+H9R0dMHciFlZWaFPnz7o06cPIiIiMGnSJERGRmL8+PFqn/Phhx/i8OHDWLFiBTw8PGBtbY0333wTJSUlCuXMzc0V/pbJZCgrKwMAXvd037NnD0JCQvD555/D398f9vb2+Oyzz5CYmKhQztbWVuFvPnWrcunSJbz99ttYtmwZ+vbty20vKytDkyZNcOLECaXnqJqRI9Tu3bu5D5369etz25mK+5HLx1Zxu6rxa3qemZkZ4uLicPr0aRw5cgRffPEFPv74YyQmJgq6SD9p0iSMGTMGMTEx2LJlC0aOHKkQsFUxNzdH165d0bVrV4SHh2Px4sWIiopCWFhYlRf63d3dVb7e1tbWSuMdP348/v33X6xatQqurq6wtLSEv78/t4/q4iZMql5nVds1/T+o6ejIw4S0adMGhYWF3N/m5uYoLS1VKHPy5EmMHz8eQ4cORdu2bdG4cWPcunVLcDvnz5/Hs2fPuG0JCQlK7XTu3BkzZszAq6++Cg8PD+6bfnXrriwnJweBgYEYNmwYQkJCFB7r0KEDsrOzUbt2bXh4eCj8yG9527p1a6U2qmpTzsXFBS1btlQIHPJxnD59WiEYnj59Gvb29mjWrJna+vg8TyaToUuXLli4cCFSU1NhYWGBffv2qazPwsJCaR8AgIEDB8LW1hbr16/Hr7/+igkTJvAab+W+vnjxAkVFRYKfq8nJkycxe/ZsDBw4EF5eXrC0tMTDhw+5x729vXH37l1cu3ZN5fPVjbly3zMzM3Hnzh1u26VLl5Cbm4vWrVvrZiAmjoKHEcrJyUHPnj2xc+dOpKWlISMjA99//z2io6MV7oHt5uaGo0ePIjs7mzu14eHhgb179+LcuXM4f/48Ro8eLfib1OjRo1GrVi1MnDgRly5dwqFDh7BixQqFMh4eHkhOTsbhw4dx7do1zJ8/n7toX926Kxs2bBisra2xYMECZGdncz+lpaXo3bs3/P39ERQUhMOHD+PWrVs4ffo0PvnkEyQnJwMA3n//fWzevBmbN2/GtWvXEBkZiYsXLwp6TSqbMWMG7ty5g1mzZuHKlSv4+eefERkZiblz52qcilrV8xITE/Hpp58iOTkZmZmZ2Lt3L/7991+1H3hubm5ITEzErVu38PDhQ+69NjMzw/jx4/HRRx/Bw8OjytN0PXr0wMaNG5GSkoJbt27h0KFDmDdvHgICAhSm4OqCh4cHduzYgcuXLyMxMRHvvPOOwtFG9+7d0a1bNwwfPhxxcXHIyMjAr7/+yt1G1c3NDQUFBTh69CgePnyIp0+fKrXRu3dveHt745133sHZs2dx5swZjB07Ft27d1c6lUjUkOpiC9FeUVERCw8PZx06dGAODg7MxsaGeXp6sk8++YQ9ffqUK/fLL78wDw8PVrt2bW4abEZGBgsICGDW1tbMxcWFrV27VumiqqoLn+3ateOmPDLG2F9//cXatWvHLCwsWPv27dmPP/6ocAG5qKiIjR8/njk4OLC6deuy6dOns/DwcIWL0JUv7vKtu/IFc6i5Vah8qm5eXh6bNWsWa9q0KTM3N2cuLi7snXfeUbhYumTJEtagQQNmZ2fHxo0bx0JDQwVP1a2Mz1Tdyhezq3repUuXWL9+/VjDhg2ZpaUla9WqlcJkh8qv6dWrV1mnTp2YtbW1wmvCGGM3btxgAFh0dLTaMch9+umnzN/fnzk6OjIrKyvWokULNnv2bPbw4UOujLZTdSs7e/Ys8/HxYZaWluyll15i33//vVLdOTk57L333mP169dnVlZW7JVXXmEHDhzgHp82bRqrX7++TqbqVhQTE8P9X6oINfCCOd2GlpAa6s8//0SPHj1w9+5dNGrUqNr1ubm5Yc6cOTpZ0W5sZDIZ9u3bJ1lqGCnQaStCapji4mKkp6dj/vz5eOutt3QSOOTCwsJgZ2eH3NxcndVpyKZNm6ZytXxNQEcehNQwW7duxcSJE9G+fXv88ssvGi/gC3H79m08f/4cANCiRQvR0owYkgcPHiAvLw9A+TTeyrPnTBkFD0IIIYKZ/lcDQgghOkfBgxBCiGAUPAghhAhGwYMQQohgFDwIIYQIRsGDEEKIYBQ8CCGECEbBgxBCiGAUPAghhAj2f316eZetMj2GAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df=pd.read_csv(r\"https://raw.githubusercontent.com/GeostatsGuy/GeoDataSets/master/unconv_MV_v2.csv\")[[\"Por\",\"TOC\"]].iloc[0:100]\n",
    "\n",
    "x = StandardScaler().fit(df).transform(df)\n",
    "plt.scatter(x[:,0],x[:,1],color='darkorange',edgecolor='black',s=10); plt.xlim([-3,3]); plt.ylim([-3,3])\n",
    "add_grid(); plt.xlabel('Standardized Porosity S[Fraction]'); plt.ylabel('Standardized TOC S[fraction]'); plt.title('Two Correlated Spatial Features')\n",
    "\n",
    "plt.subplots_adjust(left=0.0, bottom=0.0, right=0.5, top=0.7, wspace=0.1, hspace=0.1); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interactive Feature Projection with Arbitrary Rotation\n",
    "\n",
    "#### Michael Pyrcz, Professor, The University of Texas at Austin \n",
    "\n",
    "Observed the partitioning of variance and corretion over 2 new features through orthogonal projection / rotation.\n",
    "\n",
    "##### The Inputs\n",
    "\n",
    "* **Angle**: data rotation angle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 323,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dashboard(Angle):\n",
    "    \n",
    "    fig = plt.figure(constrained_layout=False)\n",
    "    gs = GridSpec(2, 2, figure=fig)\n",
    "    \n",
    "    ax1 = fig.add_subplot(gs[:, 0])\n",
    "      \n",
    "    base = plt.gca().transData\n",
    "    #print(base)\n",
    "    rot = transforms.Affine2D().rotate_deg(int(Angle))\n",
    "    #line=ax16.plot(x[:,0],x[:,1], 'o', transform= rot + base, c = 'black', alpha = 0.3)\n",
    "    line=ax1.plot(norm[:,0],norm[:,1], 'o', c = 'black', alpha = 0.3)\n",
    "    \n",
    "    xdata=x[:,0]*math.cos(math.radians(int(Angle)))-x[:,1]*math.sin(math.radians(int(Angle)))\n",
    "    ydata=x[:,1]*math.cos(math.radians(int(Angle)))+x[:,0]*math.sin(math.radians(int(Angle)))\n",
    "    \n",
    "    eigen = np.zeros([2,2])\n",
    "    eigen[0,0] = math.cos(Angle*math.pi/180.0)\n",
    "    eigen[1,0] = math.sin(Angle*math.pi/180.0)\n",
    "    eigen[0,1] = -1*math.sin(Angle*math.pi/180.0)\n",
    "    eigen[1,1] = math.cos(Angle*math.pi/180.0)\n",
    "    \n",
    "    df2 = pd.DataFrame({'x':xdata, 'y':ydata})\n",
    "    data = df2.values\n",
    "    lists=[]\n",
    "    \n",
    "    ydataZeroed = np.zeros(len(ydata))\n",
    "\n",
    "    rotinv = transforms.Affine2D().rotate_deg(int(-Angle)) \n",
    "    ax1.plot(xdata, ydataZeroed,\"or\", c = 'red', alpha = 0.3,transform= rotinv + base,label=r'$C_1$')\n",
    "    ax1.plot(ydataZeroed, ydata,\"or\", c= 'blue', alpha = 0.3,transform= rotinv + base,label=r'$C_2$')\n",
    "    ax1.set_xlim(left=-3.5, right=3.5); ax1.set_ylim(bottom=-3.5, top=3.5)\n",
    "    ax1.set_title(\"Data and Arbitrary Feature Projection Components\"); ax1.set_xlabel(r'Standardized Porosity, $X_1$'); ax1.set_ylabel(r'Standardized TOC, $X_2$')\n",
    "    ax1.annotate(r'$C_1=X_1 \\cdot COS \\left(\\alpha \\cdot \\frac{180}{\\pi} \\right)-X_2 \\cdot SIN \\left(\\alpha \\cdot \\frac{180}{\\pi} \\right)$',(-3.0,-2.5)) \n",
    "    ax1.annotate(r'$C_2=X_1 \\cdot SIN \\left(\\alpha \\cdot \\frac{180}{\\pi} \\right)+X_2 \\cdot COS \\left(\\alpha \\cdot \\frac{180}{\\pi} \\right)$',(-3.0,-2.8)) \n",
    "    \n",
    "    add_grid2(ax1); ax1.legend(loc='lower right')\n",
    "    sizes = []\n",
    "    \n",
    "#     print('Your Estimated Principal Component/Eigen Vector #1 = ' + str(eigen[:,0]))\n",
    "#     print('Your Estimated Principal Component/Eigen Vector #2 = ' + str(eigen[:,1]))\n",
    "    \n",
    "    sumOfVariance=df2.var()['x']+df2.var()['y']\n",
    "    sizes.append(df2.var()['x']/sumOfVariance)\n",
    "    sizes.append(df2.var()['y']/sumOfVariance)\n",
    "    \n",
    "    ax2 = fig.add_subplot(gs[0, 1])\n",
    "    \n",
    "    n = ax2.pie(sizes, autopct='%1.1f%%',colors = ['lightcoral','royalblue'],shadow=True,startangle=90)\n",
    "    n[0][0].set_alpha(1.0); n[0][1].set_alpha(1.0)\n",
    "    ax2.axis('equal')\n",
    "    labels = [r'$\\frac{\\sigma_{C_1}^2}{\\sigma_{X_1+X_2}^2}$', r'$\\frac{\\sigma_{C_2}^2}{\\sigma_{X_1+X_2}^2}$']\n",
    "    ax2.legend(sizes, labels=labels,loc='upper left')\n",
    "    ax2.set_title('Components\\' Proportion of Variance')\n",
    "#    plt.tight_layout()\n",
    "\n",
    "\n",
    "    ax3 = fig.add_subplot(gs[1, 1])\n",
    "    nAngle = 30\n",
    "    var_pc1 = np.zeros(nAngle); var_pc2 = np.zeros(nAngle); corr = np.zeros(nAngle)\n",
    "    \n",
    "    for iAngle, lAngle in enumerate(np.linspace(0,180,nAngle)): \n",
    "        xdata=x[:,0]*math.cos(math.radians(int(lAngle)))-x[:,1]*math.sin(math.radians(int(lAngle)))\n",
    "        ydata=x[:,1]*math.cos(math.radians(int(lAngle)))+x[:,0]*math.sin(math.radians(int(lAngle)))\n",
    "        var_pc1[iAngle] = np.var(xdata); var_pc2[iAngle] = np.var(ydata) \n",
    "        corr[iAngle] = np.corrcoef(xdata,ydata)[0,1]\n",
    "\n",
    "    ax3.plot(np.linspace(0,180,nAngle),var_pc1/np.full((nAngle),2.0)-0.006,color='red',lw=2)\n",
    "    ax3.plot(np.linspace(0,180,nAngle),var_pc1/np.full((nAngle),2.0)+0.006,color='teal',lw=2)\n",
    "    #ax3.plot(np.linspace(0,180,nAngle),var_pc2/np.full((nAngle),2.0),color='blue',lw=2,label = r'Prop $\\sigma_{x_2}^2$')\n",
    "    ax3.fill_between(np.linspace(0,180,nAngle),var_pc1/np.full((nAngle),2.0),np.full((nAngle),0.0),color='red',alpha=0.4,label = r'$\\frac{\\sigma_{C_1}^2}{\\sigma_{X_1+X_2}^2}$',zorder=2)\n",
    "    ax3.fill_between(np.linspace(0,180,nAngle),np.full((nAngle),1.0),var_pc1/np.full((nAngle),2.0),color='blue',alpha=0.4,label = r'$\\frac{\\sigma_{C_2}^2}{\\sigma_{X_1+X_2}^2}$',zorder=2)\n",
    "    ax3.plot([Angle,Angle],[0,1.0],color='black',lw=3,ls='--',zorder=500)\n",
    "    \n",
    "    ax4 = ax3.twinx()\n",
    "    ax4.plot(np.linspace(0,180,nAngle),corr,color='black',lw=4,label = r'$\\rho_{C_1,C_2}$',zorder=101)\n",
    "    ax4.plot(np.linspace(0,180,nAngle),corr,color='white',lw=6,zorder=100)\n",
    "    ax4.plot([0,180],[0,0],color='black',lw=2,zorder=101)\n",
    "    ax4.plot([0,180],[0,0],color='white',lw=4,zorder=100)\n",
    "    add_grid2(ax3); plt.xlim([0,180]); ax3.set_ylabel('Proportion of Variance ($C_1$|$C_2$)')\n",
    "    ax3.set_ylim([0,1]); ax4.set_ylim([-1,1]); ax3.legend(loc='lower left'); ax4.legend(loc='lower right')\n",
    "    ax4.set_ylabel('Correlation'); ax3.set_title('Components\\' Variance Proportions and Correlation for all Angles')\n",
    "    ax3.set_xlabel(r'Rotation Angle ($\\alpha$)')\n",
    "    newYlabel = ['0.0|1.0','0.2|0.8','0.4|0.6','0.6|0.4','0.8|0.2','1.0|0.0']\n",
    "    ax3.set_yticklabels(newYlabel)\n",
    "    plt.subplots_adjust(left=0.0, bottom=0.0, right=1.8, top=1.1, wspace=0.2, hspace=0.1); plt.show()\n",
    "    \n",
    "title = widgets.Text(value='                           Understanding PCA with an Interactive Orthogonal Feature Projection, Professor Michael J. Pyrcz, The University of Texas at Austin',layout=Layout(width='950px', height='30px'))\n",
    "style = {'description_width': 'initial'}\n",
    "widget_angle = widgets.IntSlider(min=0, max = 180, value = 0, step = 5, description = r'Rotation Angle ($\\alpha$)',orientation='horizontal',continuous_update=False,layout=Layout(width='950px', height='30px'),style=style)\n",
    "uik2 = widgets.VBox([title,widget_angle],)\n",
    "interactive_plot = widgets.interactive_output(dashboard, {'Angle': widget_angle})\n",
    "interactive_plot.clear_output(wait = True)               # reduce flickering by delaying plot updating"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interactive Feature Projection with Arbitrary Rotation Angle to Explore Principal Component Analysis\n",
    "\n",
    "#### Michael Pyrcz, Professor, The University of Texas at Austin \n",
    "\n",
    "Observed the partitioning of variance over 2 new features through orthogonal projection / rotation.\n",
    "\n",
    "The Inputs: **Angle**: data rotation angle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 324,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "efcf7b33c6e74a1995cc48cd85e6ec06",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(Text(value='                           Understanding PCA with an Interactive Orthogonal Feature…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "15f7f96a8cfd4b9a90e91d43fbb45b35",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output(outputs=({'output_type': 'display_data', 'data': {'text/plain': '<Figure size 640x480 with 4 Axes>', 'i…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(uik2,interactive_plot)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comments\n",
    "\n",
    "This was a basic demonstration of PCA as a orthogonal transformation (a rotation as we preserve the inner product space, the pairwise distances and angles between all points).\n",
    "\n",
    "* by rotating the orthogonal components in 2D we can maximize the variance explained by a single component while reducing the absolute correlation between the components to 0.0.\n",
    "\n",
    "I have other demonstrations on the basics of working with DataFrames, ndarrays, univariate statistics, plotting data, declustering, data transformations, trend modeling, multivariate analysis, inferntial and predictive machine learning, deep learning and many other workflows available at https://github.com/GeostatsGuy/PythonNumericalDemos and https://github.com/GeostatsGuy/GeostatsPy. \n",
    "  \n",
    "I hope this was helpful,\n",
    "\n",
    "*Michael*\n",
    "\n",
    "#### The Author:\n",
    "\n",
    "### Michael Pyrcz, Professor, The University of Texas at Austin \n",
    "*Novel Data Analytics, Geostatistics and Machine Learning Subsurface Solutions*\n",
    "\n",
    "With over 17 years of experience in subsurface consulting, research and development, Michael has returned to academia driven by his passion for teaching and enthusiasm for enhancing engineers' and geoscientists' impact in subsurface resource development. \n",
    "\n",
    "For more about Michael check out these links:\n",
    "\n",
    "#### [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)\n",
    "\n",
    "#### Want to Work Together?\n",
    "\n",
    "I hope this content is helpful to those that want to learn more about subsurface modeling, data analytics and machine learning. Students and working professionals are welcome to participate.\n",
    "\n",
    "* Want to invite me to visit your company for training, mentoring, project review, workflow design and / or consulting? I'd be happy to drop by and work with you! \n",
    "\n",
    "* Interested in partnering, supporting my graduate student research or my Subsurface Data Analytics and Machine Learning consortium (co-PIs including Profs. Foster, Torres-Verdin and van Oort)? My research combines data analytics, stochastic modeling and machine learning theory with practice to develop novel methods and workflows to add value. We are solving challenging subsurface problems!\n",
    "\n",
    "* I can be reached at mpyrcz@austin.utexas.edu.\n",
    "\n",
    "I'm always happy to discuss,\n",
    "\n",
    "*Michael*\n",
    "\n",
    "Michael Pyrcz, Ph.D., P.Eng. Professor, Cockrell School of Engineering and The Jackson School of Geosciences, The University of Texas at Austin\n",
    "\n",
    "#### More Resources Available at: [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}