week09d-monterey_bay_kelp.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Monterey Bay Kelp PCA\n",
    "\n",
    "This example uses data from the supplementary material of:\n",
    "\n",
    "Graham, M. H. (2003), Confronting multicollinearity in ecological multiple regression. _Ecology_, 84: 2809–2815. doi:10.1890/02-3114 \n",
    "http://onlinelibrary.wiley.com/doi/10.1890/02-3114/abstract\n",
    "\n",
    "The details of observational methods can be found in an earlier paper:\n",
    "\n",
    "Graham, M. H. (1997) Factors determining the upper limit of giant kelp, _Macrocystis pyrifera_ Agardh, along the Monterey Peninsula, central California, USA, _J. Exp. Mar. Bio. Ecol._, 218, 127-149."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Key concept\n",
    "\n",
    "The final exercise in this notebook is:\n",
    "\n",
    "Which principal component is most highly correlated with kelp extent? What is your interpretation of the relationship between kelp extent and environmental variability?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "from scipy import linalg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df = pd.read_csv('data/Graham2003_Kelp/Originaldata.txt',delimiter=' ')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Response</th>\n",
       "      <th>OD</th>\n",
       "      <th>BD</th>\n",
       "      <th>LTD</th>\n",
       "      <th>W</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>3.241</td>\n",
       "      <td>2.0176</td>\n",
       "      <td>4.87</td>\n",
       "      <td>-0.59</td>\n",
       "      <td>-4.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3.032</td>\n",
       "      <td>1.9553</td>\n",
       "      <td>4.78</td>\n",
       "      <td>-0.75</td>\n",
       "      <td>4.7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3.100</td>\n",
       "      <td>1.8131</td>\n",
       "      <td>3.14</td>\n",
       "      <td>-0.38</td>\n",
       "      <td>-4.9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3.156</td>\n",
       "      <td>2.5751</td>\n",
       "      <td>3.28</td>\n",
       "      <td>-0.16</td>\n",
       "      <td>-3.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>3.110</td>\n",
       "      <td>2.2589</td>\n",
       "      <td>3.28</td>\n",
       "      <td>0.01</td>\n",
       "      <td>5.6</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Response      OD    BD   LTD    W\n",
       "0     3.241  2.0176  4.87 -0.59 -4.1\n",
       "1     3.032  1.9553  4.78 -0.75  4.7\n",
       "2     3.100  1.8131  3.14 -0.38 -4.9\n",
       "3     3.156  2.5751  3.28 -0.16 -3.2\n",
       "4     3.110  2.2589  3.28  0.01  5.6"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Variables:\n",
    "\n",
    "* Response - shallow (upper) distributional limit of kelp Macrocystis pyrifera [m]\n",
    "* OD - wave orbital displacement [m]\n",
    "* BD - wave breaking depth [m]\n",
    "* LTD - minimum daily tidal height [m]\n",
    "* W - wind velocity [m/s], positive onshore"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'depth [m]')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEICAYAAABS0fM3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO2deXzbd33/n2/Z8m1Lvu04d+M0adIkTdPScrW0hR6c49iAAaUbK904N8YYY4xjMO7RMdivtOXaBoOWFii9y1HaUtqSpEmbO2mcy7djy7It27o+vz+kr6zYkiw7uvV+Ph5+xJK+kt5WJL2/7+v1FmMMiqIoSvFiy7YBiqIoSnZRR6AoilLkqCNQFEUpctQRKIqiFDnqCBRFUYocdQSKoihFjjoCJecQkWMictUi7vdpEfnfdNiUTaL/LhFZLiLjIlKyyMcaF5HVqbVQyXdKs22AoijJY4w5AdScxf0j9xWR7wOnjDH/nALTlDxGIwJFySFERE/OlIyjjkDJaURknYh0ichbw5eXiMhdIjIYvv6Dce63UkSMiNwoIj0i0isiH0nwPI+KyHuiLr9bRJ6IumxE5IMiclREhkTkKyJiizr29yLynyIyKiIHROTKqPs6ROQ7YRu6ReRzVmon6r5fF5Fh4NPzvB7W31UaZffnROTJcNrnlyLSKCI/FBG3iPxRRFbO+jvWiMiNwJ8D/2DdL9HzKoWNOgIlZxGRrcDDwAeMMT8Of/H+EtgNdABXAh8WkasTPMwrgE7gVcA/Lqb2EMWfANuArcDrgb+Iuu1FwFGgCfgUcLeINIRv+wHgB9YAF4RteU+M+7YAn1+EXW8F3knoNTkH+APwPaAB2B+25wyMMbcCPwS+bIypMca8dhHPqxQI6giUXOVlwD3A9caYe8PXXQQ0G2M+a4zxGmOOArcR+iKMx2eMMRPGmOcJfTm+7Sxs+pIxZjicp7951mMNADcbY3zGmJ8AB4FXi0grcC3w4bAdA8DXZ9ncY4z5T2OM3xgzuQi7vmeMecEYMwo8ALxgjPmVMcYP3EnI+ShKXDQfqeQqNwG/M8b8Nuq6FcASEXFFXVcCPJ7gcU5G/X4cOP8sbJr9WEuiLnebMxUcrdtXAHagV0Ss22yzHiv698XQH/X7ZIzLiy4uK8WBRgRKrnITsFxEvh513UmgyxjjjPqpNcZcl+BxlkX9vhzoiXPcBFAVdbltgY/VIVHf9FG3nwSmgaYom+uMMRuijs2WBLBKDyuAOgIldxkDrgFeLiJfDF/3DOAWkY+JSKWIlIjIRhG5KMHjfFJEqkRkA3AD8JM4x+0C3hg+dg3wlzGO+aiI1IvIMuBDsx6rBfigiNhF5C3AeuB+Y0wvoTrH10SkTkRsInKOiFyW5OuQTvoBnSlQ1BEouYsxxgW8ErhWRP7VGBMAXgtsAbqAIeB2wJHgYX4HHAF+DXzVGPNwnOO+DngJfTn+gFAhdTa/AHYQchr3Ad+Juu1pQkXpIUIF3zcbY06Hb3sXUAbsA0aAnwLtCWzOFN8BzhMRl4j8PNvGKNlDdDGNUoiEWya7AHu4aHq2j2eATmPMkRi3vRt4jzHmpWf7PIqSDTQiUBRFKXLUESiKohQ5mhpSFEUpcjQiUBRFKXLybqCsqanJrFy5MttmKIqi5BU7duwYMsY0x7ot7xzBypUr2b59e7bNUBRFyStE5Hi82zQ1pCiKUuSoI1AURSly1BEoiqIUOeoIFEVRihx1BIqiKEWOOgJFUZQiRx2BoihKkaOOQFHSxH3P9TI4Np1tMxRlXtQRKEoaGByb5n0/2skd2892C6WipB91BIqSBg71jwEwPOHNsiWKMj/qCBQlDViOwOXxZdkSRZmftDkCEakQkWdEZLeI7BWRz8Q4xiEiv4w65oZ02aMomeRQ/zgALo9GBEruk07RuWngCmPMuIjYgSdE5AFjzFNRx7wP2GeMea2INAMHReSHxhj99Ch5zWErIpjUiEDJfdLmCExo4814+KI9/DN7C44BakVEgBpgGDjr/bKKkk2MMZHU0IhGBEoekNYagYiUiMguYAB4xBjz9KxDvgmsB3qA54EPGWOCMR7nRhHZLiLbBwcH02myopw1A2PTuKf82EuEUa0RKHlAWh2BMSZgjNkCLAUuFpGNsw65GtgFLAG2AN8UkboYj3OrMWabMWZbc3PMvQqKkjNY0cCmpU5ckz50HayS62Ska8gY4wIeBa6ZddMNwN0mxBGgC1iXCZsUJV1YheKLVzUQCBrGpjXbqeQ26ewaahYRZ/j3SuAq4MCsw04AV4aPaQXOBY6myyZFyQSH+8dorC5jdVM1AK4JTQ8puU06u4bagR+ISAkhh3OHMeZeEbkJwBhzC/CvwPdF5HlAgI8ZY4bSaJOipJ2D/WN0ttZQX1UGhArGyxursmyVosQnnV1DzwEXxLj+lqjfe4BXpcsGRck0xhiO9I/zJ1s7cFbZAW0hVXIfnSxWlBTSOzrF2LSfztZanOGIQIfKlFxHHYGipBCrY2htS81MRKAtpEqOk84agaIUHYfDHUNrW2uprQh9vHSoTMl11BEoSgo51D9GU0059dWhtFBtealGBErOo6khRUkhhwbGWdtaE7nsrLZrjUDJedQRKEqKCHUMjbG2tTZynbOyTLuGlJxHHYGipIhu1yQT3gCd0RFBlZ0RTQ0pOY46AkVJEdGFYgtnVRmjmhpSchx1BIqSImZaR2ccQb1GBEoeoI5AUVLEof5xWmrLcYTnBwCclXbcUz4CQVUgVXIXdQSKkiIOD5xZKIZQasgYcGvBWMlh1BEoSgoIBg2H+8fPKBQDqjek5AXqCBQlBXS7Jpn0BeZEBNEKpIqSq6gjUJQUECkUz4oIrHqBrqxUchl1BIqSAqytZGtaNCJQ8g91BIqSAg73j9FWV4Gj0n7G9c5KVSBVch91BIqSAg4NjM0pFAPUVdoR0Z0ESm6jjkBRzpJg0HBkYHxOoRigxCbUVdi1a0jJadQRKMpZcnLEw5QvOKdQbKHTxUquo45AUc4Sq1DcGSMiAHBUlWlqSMlp1BEoyllitY52tsSPCLRYrOQy6ggUJQaT3gAvDI4ndezh/jGWOCqorbDHvN1Zacc1qRGBkruoI1CUGHz3911c/fXHODIwNu+xh/rH46aFIKQ35JrQiEDJXdQRKEoMDvSN4Q8a/u3+AwmPCwQNLwyOxy0UQ0hvaGzajy8QTLWZipIS1BEoSgy6hsYptQm/OTDA748MxT3uxLCHaX8wYURgTRePagupkqOoI1CUWRhj6Bqc4C3bltHhrORz9+2Pu09gRmMoUWpIp4uV3EYdgaLMYnB8mglvgHVttXzs2nXs73Vz185TMY89PE/HEIRqBKDTxUruoo5AUWbRNTgBwKqmal67qZ0ty5x89aGDeLz+Occe6h+nw1lJdXlp3MdTvSEl11FHoCiz6BqacQQiwidfs56BsWlufezonGMP9Y8lLBSDKpAquY86AkWZRdfpCcpKbSxxVgJw4YoGrju/jW//7ij97qnIcf5AkKODEwnrAxC1k0CLxUqOoo5AUWbRNTjBioYqSmwSue5j16zDHwzytYcPRq47PuzBG0jcMQRQV1FKiU00IlByFnUEijKLrqEJVjVVn3HdisZqrr90JXfuOMW+HjcwUyieLzUkIjgqVWZCyV3UEShKFIGg4fhpD6uaq+fc9oErOnFU2vm3+/djjInaSpbYEUCohVQdgZKrqCNQlCh6XJN4A0FWNc51BI4qOx+8opMnjgzx6MFBDvWPsayhkqqy+B1DFqo3pOQy6giUgmRofJopX2DB94vuGIrFOy5ZwcrGKj5//34O9I2xtiVxfcCivqqMEdUbUnIUdQRKQfK6/3yC//j14QXfL+IIYqSGAMpKbfzjtes5MjDOkYHEYnPROKrs2jWk5CzqCJSCw+sP0jM6xfZjwwu+b9fQBDXlpTTXlMc95uoNrVy8sgGYv1BsUV9Vpl1DSs6ijkApOKxc/L4eN8E4GkHx6BqaYGVTFSIS9xgR4V9eex6dLTVcvKohqcetr7Lj8QaY9i88XaUo6SZtjkBEKkTkGRHZLSJ7ReQzcY67XER2hY/5XbrsUYoHKxc/4Q1wfNizoPuGWkfnP8vf2OHgkb+7jKX1VUk9rsNSINXOISUHSWdEMA1cYYzZDGwBrhGRS6IPEBEn8F/A64wxG4C3pNEepUiITsHs6R5N+n5ef5BTI564heKzod5SINU6gZKDpM0RmBDWrj97+Gd2nP524G5jzInwfQbSZY9SPIxMzDiCveHhr2Q4MewhaGBVU3Jn+QvBWVk2xzZFyRXSWiMQkRIR2QUMAI8YY56edchaoF5EHhWRHSLyrjiPc6OIbBeR7YODg+k0WSkARsLpl6aaMvb2JB8RzLSOJlcAXghOjQiUHCatjsAYEzDGbAGWAheLyMZZh5QCFwKvBq4GPikia2M8zq3GmG3GmG3Nzc3pNFkpAKzU0IvPaWJfjxtjkisYdw2FAthYw2Rny8xyGo0IlNwjI11DxhgX8ChwzaybTgEPGmMmjDFDwGPA5kzYpBQuIxNeqspKuHBFPacnvPRFKYYmomvIQ0N1WUQtNJXUR5bTaESg5B7p7BpqDheDEZFK4Cpg9ibwXwAvE5FSEakCXgTsT5dNSnEw7PFSX1XGhiV1AOztTq5O0DU0npZCMUBVWQn2EomkrRQll0hnRNAO/FZEngP+SKhGcK+I3CQiNwEYY/YDDwLPAc8Atxtj9qTRJqUIcHl81FfbWd9eh0jyBeNYqqOpQkRwVpUxqnpDSg4yv1rWIjHGPAdcEOP6W2Zd/grwlXTZoRQfwxOhiKC6vJRVjdVJFYwnpv30u6fT5gggJDynekNKLqKTxUrB4QqnhgA2dDiSigiOnU4sNpcK6qvKVIFUyUnUESgFRygiCBV8Nyypo9s1OW///nyqo6nAoTsJlBxFHUGR8Ytd3bz11j9k24y04Q8EcU/5cVoRQbhgvK83cVTQNRhyBCvT0DpqUa+OQMlR1BEUGU8eOc1TR4cXpdWfD1gDWw3VliNwAMxbJ+g6PUG7o4LKspK02eZUBVIlR1FHUGRYPfWF+oVkDWxZA1wN1WUscVTMWydIZ8eQhbPKzrQ/WLBOWMlf1BEUGX2jYUdQoN0rwxNnRgQA5y1xzCs+lxFHYOkNFagTVvIXdQRFRu/oJFC4X0bW32V1DUGoTnB0aAKP1x/7PhNeXB5f2h1BRIFU6wRKjqGOoIjweP24p0JfhsMFqoJpdQfVV5/pCIyB/b1jMe/TlYHWUSAiXVGoTljJX9QRFBFWWggK98vIknCoj9IL2tARKhjvi1MwtjqG0h8R6HIaJTdRR1BERIuvFWqNYMTjpbzURqV9pvtniaOC+io7e+JoDnUNTVBiE5Y1pH4PQTTOSERQmK+9kr+oIygiiiIiCMtLRO8cFhE2LHGwtzdORDA0wbL6Suwl6f04RBRIdbpYyTHUERQRvWFH0FpXXrg1Ao/3jPqAxYYldRzqG8cXCM65LRMdQwAV9hLKS21aLFZyDnUERUS/ewpHpZ0lzsrCjQg8vjPqAxbnLanDGwhyuH/8jOuNMUkvrE8F9VVlupxGyTnUERQRvaNTtNVVUF/AE64jE7Ejgo3hgvGeWQXjfvc0k75AWvYUx8JZZdcawQL5/ZEh/nhsONtmZJ2+0am0DSOqIygi+kanaHOEHUEBF4tjRQSrGqupKith36wJ43TuKY6Fs8quXUML5F/v3cdXHzqYbTOyzuVf/S1fezg9r4M6giKizz1Fu6OChmp7QdYIAkGDa9JHQ9XciMBmE9a3183RHIo4gub01wggNF1cqNFYOjDGcGpkkqHx6WybklWmfAGmfMGImGKqUUdQJPgCQYbGp2mtq8BZVcakL1BwmjfuSR/GEPfDsmFJHft63ASDM8vsu4bGKS+10V5XkREb66vtEWE8ZX7ck37Gp/0MjRe383RF5mPUEShnwcDYNMYQjggKU/NmOPz3NMSoEUDIEUx4A5ElNBCKCFY2VmOzScz7pBpHZahYbIyZ/2CFUy4PAKOTPrz+uR1fxcLILDHFVBN3VWV41/B8DBpjrkyhPUqa6AtrDLU5KpjyhT5QwxNe2h2V2TQrpcxWHp3NjCS1m9XNoZpA19AEnS21mTGQ0MSzL2DweANUl6dtU2zB0D0yGfl9eMJLmyMzkVuukTVHAJQA1yW4XYB7UmuOki6sGYI2RwXuyZDeUKEVjEdiKI9Gs7a1FnuJsLfHzWs3L8EfCHJi2MOrNrRlzEZnlN6QOoL56XbNOIKh8emidQTpTg0leie+1xhzPNGdReRvUmyPkiasqeL2ukpKJPT7cIGmhuJ9WMpKbXS21EYKxt2uSXwBw6o0biWbjVW/cHl8LK3P2NPmLdERQTEXjLNWIzDGPDHfnZM5RskN+kanqLDbqKssjfTZF9pgk/X3xJojsLAKxtYgGWSuYwjAWalS1Auh2zVJhT30NVXMBeN0p4bmLRaLyGtE5FkRGRYRt4iMiUjidU9KzhFqHa1ERCJfRoXWQjo84cNeIlQnWDe5scPB6Qkvfe6pjCysn03ECaveUFJ0uyYjtZ3TRR0ReKmw26iwp2eVajJdQzcD1wONxpg6Y0ytMaYuLdYoaaMvPFUMUFpiw1Fpj2j3Fwouz1zBudlYy+z3drvpGpqgtqKUxgQRRKqxnLBOFydH98gka1trqLDbijo1FJJOSd/7NBlHcBLYY7TfLa/pDU8VW9RX2RkusC+j4bDyaCLWt9chEuocssTmEjmOVGMtpxktsLRcOpj0Bjg94aXDWUlTTXlRp4ZcHl/ahskgcbHY4h+A+0Xkd0DEJRtj/j1tVikpJRg0DIzNcgTVhSd+5vL4qK9OnEOtLi9lVWM1e3tG6Rqa4MIVma3YlpeWUFVWohFBElgdQx31lTTWlBd1RODyeCPRZDpIJiL4POABKoDaqB8lTzg94cUXMJHUEEBDVVnh1Qg880cEEFIiffaki27XZEbrAxYhBVJ1BPMRcQTOKppryoo6IgjJq6fPESQTETQYY16VNguUtNMXNUNgUV9dxv7ewqr5uzzepMLnjR0O7n2uF8hsodjCUWkvuGgsHVitox31lTRWl7P7VOzFQsVAulNDyUQEvxIRdQQ5gDGGy77yW/77D8cWdD9rRWX7nBpB4XwZGWMY8fhoSOKsySoYQ3YcgeoNJUe3y0OJTWitLaepNhTBRutEFQvGhMQUY6nqpopkHMH7gAdFZFLbR7PL6KSP46c9PHX09ILuF5GXqDszIpjyBZn0FobwnHvKTyBokkoNWe2IACuz4AhUgTQ5ukcmaauroLTERlNNeURdttgYmw69t52VWSwWG2O0HpAjWGf2s7dsJXO/UpvQWFMeuc6Sah7xeKksy3+9IasVNhlH0FBdRrujAl8gSF1F+s6y4qE7CZKj2zVJR33ovWm9d0+PT8eVEClUXGHplHQNk0GCiEBE5hVgSeYYJXVYuf6uoYmYu3fj0Ts6RUttOSVRCpvWYFOhFIxHIlPFyX1YrlzfwovPaUqnSXFxVoVSQ9qRnZjukUmWOkOOoKkm9H4dLMLOoZF5pFNSQaKI4H5g6zz3T+YYJUX0hyMCf9Bw/PQEa5JUzeybNUMAM2+qQklRLPTD8rk3nJ9OcxJSX1VGIGgYm/ZnJSLJB3yBIH3uqUhE0ByOCIqxc8hKh2UlIgA2h2sC8X7GgNa0WabMoW905mzoyEDy6SFLXiIaq6haKP3slvJoOs+aUoXD0hsqMPXXVNI3OkXQQIdzbmqo2JiRV89C15AxpiQsKRHvp9YY05E2y5Q59LmnqAlLFydbJzDG0Dc6RWtdnIig4FJDue8IrNde9YbiEz1MBiFpjhKbFOVQ2Uz9K7tzBEqO0O+eYnlDFaOTPg4nGRGMTfvxeANntI5C6KxUpLBqBCU2oa4i99/SMzsJNCKIR2SGIBwR2GxCQ3UZp4swNWS9TxxpnCzO/U+NEsHK9bfUlSftCKwCc+ssR1BaYqOuwl4wNYLhiVCfdSZ1gxbLzE6Cwnjt04EVESxxzqQ0m4pUZmJ00kdtRSmlJenbLKw7i/OIfncoxdPZUsPRwXECSQzXWJvJZkcEEGqjLJSz0mSninMBKyJQmYn4dI9M0lRTfobsclNNGYNFGREkJ51yNiQVEYhICaHCcOR4Y8yJdBmlzGXaH1JibKuroN1RwbQ/yKkRDyvm2a7Vb8lL1M11BPVVhSNFPTzhjcxG5Dq6nGZ+omcILJpqyiM7JIqJkAR1ervLkllM8wGgH3gEuC/8c28S96sQkWdEZLeI7BWRzyQ49iIRCYjImxdge1Ex4A6FxG2Ocs5pCS1eT6ZgbEUEs4vFECpaFkqNIKTFkh+tmKUlNmrLSwsmLZcOul0zMwQWTTVlDI1PF938RSai3WRSQx8CzjXGbDDGnB/+2ZTE/aaBK4wxm4EtwDUicsnsg8LRxpeAhxZieLFhzRC01lWwxnIESdQJ+tyTNNWUUVY697+6vrpwpA6GPd68mjh1VtsZLUK5hGQIBk3MiKCxppwpX5CJApFFSZZMnOQku5hmwbJ/JoT1TWUP/8Ry5R8A7gIGFvocxUR/JCKowFFpp7WunMMDY/PeL9YwmUVDgTgCY0xe1QhA9YYSMTQxjdcfjHQMWTQV6SxBVmsEIvJ34V+PAo+KyH0scDFN+Gx/B7AG+JYx5ulZt3cAfwJcAVyU4HFuBG4EWL58+XxPW5BYOkNWrr+zpZYXkogIekenWFofW0uovmpGeK4ywZ7fXGd82o8vYJJSHs0VnFV2rRHEYXbrqIUlMzE0Pj1vbaxQ8AeCjE35sxoRWAtoThCqD5RFXVeTzIMbYwLGmC3AUuBiEdk465CbgY8ZYxLGesaYW40x24wx25qbm5N56oKj3z1Feakt0ku8pqWGwwPj8+ZL+93xIwKrAJXvctTWF2o+RQSh5TT5/bqni9nDZBZNRSgzYclLZC0iMMZ8BkBE3mKMuTP6NhF5y0KexBjjEpFHgWuAPVE3bQN+HO79bgKuExG/MebnC3n8YsCaDrb65Dtba/B4A/SMTs05c7KY8gUY8fhidgzBzBTuSHgvbL5iFbzzpWsIZoTnlLlEL6SJZsYRFE9qaOYkJ/s1go8ned0ZiEiziDjDv1cCVwEHoo8xxqwyxqw0xqwEfgr8jTqB2PS5p874Ql/TbHUOxa8TzGwmi/0lbxVX8z1XvVDl0VzAWVXG6KQvqVmQYqPbNUltRekcQT7r/To0lt/v14WQCZ0hSFwjuBa4DugQkW9E3VQH+JN47HbgB+E6gQ24wxhzr4jcBGCMuWXxZhcf/e4pNi11Ri53toaUR48MjHP5uS0x7xNrM1k0VriZ7y2kmZDpTTXOSjvGwNhUelcQ5iPdI5MxI9SycGr09ETxRATWwGe65wgSDZT1ANuB1xEq+FqMAX873wMbY54DLohxfUwHYIx593yPWaxYwnGvOi9qsUx1GY3VZQlVSPsSzBDAzJsr34fK8kl51KI+Sv1VHcGZdLsm4zY4WLMExUKmTnIS1Qh2A7tF5EeAAOsItX8eNMbk9zdHnjE66WPaH5zzhW4VjOMR6TSKExFEhOdyrHvFGIN70o8jybOgEY8Xm0BdGkW5Uo21djAU+hdHB0yydI9M8qJVDTFva6wpL6rUkLXJLtnPwmJJpkbwSuAF4BvAN4Ej4bSRkiHifaF3ttZwuH8sbudQ3+gUteWlEenq2ZSWhELtXOte+eVzvVz8b79K+sxvxOPFEZYpzhdUbyg2o5M+xqb9cwrFFs015QwVVWrIS6lNqI3zGU4VyTiCfwdeYYy53BhzGfAK4OtptUo5g744ekFrmmtwT/kZHIv9wegdnYwbDVg05KDMxGOHBpn2B9nb407q+JEJX17sIYjGqTsJYjIzQ1AV8/ammjKG4rzfC5GR8FRxulV1k3EEA8aYI1GXj6JTwBklWl4iGqtgHC891OeentcROKtyT4p65/ERAA72JekIMjB5mWpm6jMaEUQTb4bAorGmHPeUH68/+Z3d+UymJuaTcQR7ReR+EXm3iFwP/BL4o4i8UUTemGb7FGZWVM5xBGHNoXgF477RybgzBBYN1WUM59CX0fCEl6NhhcmDfcntXAipM+aXI6itCNVndJbgTLpHPMDcqWKLiMxEkaSHXB5fRK02nSTjCCoIqY9eBlwODAINwGuB16TNMiVCn3uKxuq5wnHNteXUVZTG1BzyB4IMjk3HbR21yLUJVysacFTaOdifbGrIm/b2ulRTYpOcrM9km27XJOWltoicxGwiMhNFUjAeyVBEMG8FwhhzQ9qtUBJiLaSZjYjQ2VobU456cHyaoJm7mWw2oYjAizEmJ7Z77TgxQqlNeO3mdu7cfopA0CQsAhtjQqmhPKsRQGiWoFAWA6WKbldohiDee9FaYl8sBWOXx8f5HTkQEYjIWhH5tYjsCV/eJCL/nHbLlAiJ9II6W2pipoYSbSaLxllVxrQ/yKQvN6R9dxwfYUOHg81LnUz7gxw7nXgRyaQvwLQ/mHepIQitYTwxz99XbHS7puLWByDUNQQUTcE4Uyc5yaSGbiMkKeGDyKDYW9NplHIm8SICCM0SnJ7wzun8mdlMllhDyFLszIXOIV8gyO6TLi5cXs+6tjoADvUlltrO1ORlOtjY4WB/3xi+QHEUPpMh3lSxRWM4NXQ6B96v6WYqfJKTzqX1Fsk4gipjzDOzrktGYkJJAV5/kKFxb9yi75o4BePe0cTDZBb1kUXq2U9R7OtxM+0PcuGKejpba7AJHJjPEUxYOkP5FxFsWFKH1x9MOB1eTEz5AgyNTyd0BNXlpVTaS4oiIsikdEoyjmBIRM4hvFQmvE6yN61WKREGxqwv9PKYt8+0kJ75hdnnnqKs1DbvmbL1BZoLEcGOcKF46wonFfYSVjZWc3DeiCD/dIYsNnY4ANjTveC9TwVJzzytoxZNtcUhMzEjnZL+iCCZcbX3AbcC60SkG+gC3pFWq5QI8WYILJY4KqguK5lTMO4bDamVzlcAtr5Ac2GWYMeJETqclbSH1VLXttZyMIG6KkRJUOeR8qjFqsZqqstK2NvjZkG67gVKZIZgHkn0xuryokgNWR1l6ZaXgOS6ho4CV16g2scAACAASURBVIlINWAzxsy/H1FJGfFmCCxEhDUxCsaJVlRG05BDEcHO4yNsWzmjMXNuWy0P7etjyhegwh57g1o+LqWxsNmE85bU8bxGBED8PQSzaaop51R43qCQydRSGkhuVeXs64HkVlUqZ8/sFZWxWNNSyxNHBufc74Llzjj3mMESnst2G2OPa5Le0SkujLJ5XVstxsDh/nHOX+qIeT/LgWVi6CYdbFji4Cd/PDlvm2wx0O2apMQm8w5BNtWUseukK0NWZY9cqRFYaym3AX8NdIR/bgLOS7tlChBKDZWV2hJuKFrTUkO/exr3VOjL3JKtnu8DBTODTdmWot55IlQfuHDFmREBwIEEUhMuj5e6ilJKS5Ipd+UeGzscTPoCdA1pwbh7JDQJP9//ZVNNOcMT0wQLfKlPpraTQQJHYIz5THhdZROw1RjzEWPMR4ALCe0gVjJAMrn+2VITwxNevIFgUqkhCAvPZblGsOP4CBV2G+vaayPXrWisprzUlrBgPOzxRdJb+cj5kYJxclPUhcwpV+LWUYummjKCJjfqWunE5fFSYbfFTYumkmROo5YD0a+4F1iZFmuUOcxeURmLztawIwgXjOfbTDab+uqy7EcEx0fYvNSJPepssMQmdLbWJCwYZ0qUK12c0xxydto5FJ4hmKc+ADPTxYVeMM6khlYyjuB/gGdE5NMi8ingaeAH6TVLseh3T80rE7G0voryUlukhXS+zWSzqa8qy2qNYNIbYG+PmwtX1M+57dzWusQRwYQ3ryOC0hIb69vr2NNT3I7AHwjS555KMiIojuniTJ7kzOsIjDGfB24ARgAXcIMx5gvpNkyJzvXHniGwKLEJ5zTPbCubiQjm/1BBqE85mxHBc6dc+IMmpiNY11bLwNh0XPtcYb32fGZjRx17u90Fn/NORP/YNIGgSSoiaK4NfTkOFvgsQSgiyMx7O6kKmzFmpzHmP8I/z6bbKCVEvBWVsehsnWkh7RudwibEVXCcTUN1qEYQb9NZutkRLhRfsHyuI1gbKRjHjgqGJ7w05HFqCGDjEgdj035OFkFLZDxmFtIkkRqqDqeGxgs7NRSKCHLIESjZYb6dw9Gsaa7h1MgkHq+f3tEpWmrn776wqK8uw+sP4vFmR3hu5/ERVjdXx0zxrAs7glhLaqZ8ASZ9gbyUl4hmoxaM6XaF9xAkERE4Ku2U2qTgp4tD0W6OpIaU7BFvRWUsrILxCwMTCdVKY9GQxeliYww7jo9wYYxoAKClthxnlZ2DMaS281leIprO1hrsJVLUdYKFRAQ2m9BQXdgyE8YYXJM5lhpSssOAO/FUcTRrWmY0h3qTnCGwcGZxbWLX0AQjHl/M+gCEBhjPba2NGRFkUoslnZSXlrC2tbYgO4em/QF++PTxyIxLPLpdkzTVlCXdKtlUU17QqSH3lJ9A0ORU15CSJfrm0RmKZkVjFfYS4fDAeNLyEhYRmYksRASW0Fw8RwCh9NCh/vE5NYxIRJDnqSEI1Qn2dI9mrU6TLn62s5tP/GwPH71zd8K/7dQ88tOzaaotL+iIYDTcxZcJCWpQR5DTxFtRGQt7iY1VTdXsOuFifNq/IEdgfZFmo3No54kR6ipKOae5Ju4xa9tqGZ/2cyqcPrAolNQQhDqHRjw+esLpwELh7p3dlJXaeGhvP9/7/bG4x3W7kpshsGiqLmOogCOCTL+31RHkMP2j8RfSxKKzpTZyhp3sMBlkt0aw4/gIW1fUY0ugszNTMD6zc2hmF0F+p4YANhSgJPWJ0x6eOTbMh67s5JXntfKFB/bH1AgyxtCT5FSxhRURFFoEZTET7WpEUPT0LbDoe05LDd7wtquF1AjqLOG5DEcEo5M+DvWPxy0UW6wN71yYPWFsDcE5K/M/IljfVodNYG8BOYK7nz2FCLzhgg6++ubNtNRW8L4f7ozIK1ucnvAy5QsuyBE0VodWrI5PF+aOrEyr6qojyGESraiMhaU5BMm1nFqU2ARnpT3jNYJnT1iLaBI7gtoKOx3OyjkRwfCEl9ry0qRSZ7lOZVkJa1pq2NNTGC2kxhju3tnNpasb6XBW4qiy860/38rA2BR/f+dzZ5zJz8hPVyX9+NZ0caEWjC1nmSlV3fz/BBUo862ojIXVQgrJy0tYhPSGMts1tPP4CDaBzcvml8s+t612jiNwebw4CyAtZLGxw1EwqaEdx0c4MezhTVtn9Cm3LHPy8WvX86v9/Xznia7I9ckupImmqTYsM1GgBeMRLRYrMLOisnUeeYloVjVVY5NQO+VCFQsbqsoyXiPYcWKEdW111JTPvyjv3LZaXhgcx+ufWfQ+7PHl/VRxNBuXOBgYm2bAnf8F47t2nqLSXsI1G9vOuP6Gl6zk6g2tfPGBAxHp8WQX0kTTGG5wKFRHkGl59aJyBFO+QN7ouURWVC4gxVNeGtrz25akxlA0zqqyjG4pCwQNu064EraNRrOurRZ/0NA1NBG5Lt+VR2djTRjvzfP00JQvwL3P9XLtxjaqZzl5EeHLb95Mu7OC94frBd2uSWrLSxd09tsciQgKMzU04vFltC26aBzBL3f3cN6/PDinBTFXsVZULiQ1BPDul6zk7RcvW/DzNVTbMxoRHOwbY8IbSNoRxFpSMzzhzfthsmjOW1IH5H/n0K/29zM25edNF8ZeW+KotPOtt29lcHyaj9yxm1MjngVFAzAz+1KwEcGkL6Nb95JZXl8QLHFWEjRwqH+M5Y3JF6WyRTIrKmPxrktXLur5rBqBMWbehfepYMeJ+QfJolndVEOpTc6oE7gyfNaUbmrKS1ndVJ33UhN37ThFu6OCS1Y3xj1m01Inn7huPZ/+5T5KbMLla5sX9Bz2ktDWvkIuFmdyPqZoIoI14Y6awwP5sRIwmRWVqaS+qgxvIHPCczuPj9BcW87SJM8Ey0ptrG6ujjgCb7h1sBCGyaLZ0OHIa/G5gbEpHjs8xBsu6Jh3B/P1L17Jdee3JS0/PZummsKdLh7xZDbaLRpH4Ki001ZXweEE265yiWRWVKYSq+iaqTqBJTS3kL/v3La6iBy1q4DkJaLZuKSObtdk1jfGLZZ7dvUQCBretLVj3mNFhC++aRMvXdPEZQuMCCBUMC5UR+CayJzyKBSRI4BQe2W+RATJrKhMJRGZiQzUCQbGpjgx7Ek6LWSxrq2WbtckY1O+yMxDIdUIIEqSOk/TQ3ft7GbzUkdEBHE+6irs/O97XsSV61sX/FxNtYUpPOcLBBmb9md04VJxOYKWWo4MjOdF59BAEisqU0lDuB8/ExHBzuMhmYH5Bslmc254wvhQ/3hk5qGQ2kcBNkQKxvmXHtrX42Z/r5s3bo1dJE41zTXlBbmlbHTSUtXN3Hu7aIrFEIoIJn0Bul2TLGvI3YKxMYY+9xSvXMAMwdlihaGuFO8unvYHODnsoWvIw7GhCbpOT/D00dOUldjY2FG3oMc6N0pzyDpbKqT2UQj9PUvrK/MyIvjZs6ewlwiv3bwkI8/XWF3G2JSfaX+A8tKFzc3kMpGp4gxGBGlzBCJSATwGlIef56fGmE/NOubPgY+FL44Df22M2Z0um9a2WgXjsZx2BO5JP1O+5FZUpopU1Qh8gSC3P97Fky8M0TU0QY9rkugAzFllZ2VjNX/7ymUL/vB2OCupLivhYJ87ssIynxfXx2PjEkfeaQ75A0F+9mwPrzi3JWP/J9Z08elxL0sWMJWc64xkWGcI0hsRTANXGGPGRcQOPCEiDxhjnoo6pgu4zBgzIiLXArcCL0qXQWuaZ1ILV6xbeE4yUyxkRWWqqKu0Y5OzqxEc7h/jb+/YxZ5uNxuW1LF1eT1v3LqUVU1VrGysZlVT9Vm9uW02YW1bLQf7x2gJO8l8X1wfi40ddTy4tw/3lI+6ivz4+x4/MsTQ+HTG0kJw5nRxITkCKyrPZP0rbY7AhFSlrMqsPfxjZh3zZNTFp4C0voscVXZaass5HGPtYS6x2BmCs6HEJoueLg4GDd9/8hhffPAANeWl3PKOC+dIC6SKdW21PLinjw1LHFSVlSxYSiMfsCSp9/W4E/biB4KGfvdUTnwJ3r2zG2eVnVesW3j3z2KJjggKiWzs2UhrsVhESkRkFzAAPGKMeTrB4X8JPBDncW4Uke0isn1wcPCsbFrbWsvhgdxuIe0fTX4zWSpxVtkXXCPocU3yzu8+zWfv3cfL1jTx4IdfljYnAKH/vxGPj0P9YwU3Q2Cxccn8uwmmfAFu+t8dvOzLv+VYlOxGNnBP+Xh4bx+v27wko7n65rACaaEVjLNRI0irIzDGBIwxWwid6V8sIhtjHSciryDkCD4W63ZjzK3GmG3GmG3NzWd3xrGmpSbnO4cWsqIylTQsICIwxvDzZ7u5+ubHePaEiy+88Xxuv34bLbXptdkqGG8/NlIQC2li0VxbTltdRVzNobEpH9d/9xke2ddPIGjYHl5GlC3uf66XaX8wo2khgMaawpSZGPH4KLVJUmKMqSIj7aPGGBfwKHDN7NtEZBNwO/B6Y8zpdNuytrUWjzdAz2juag4tZEVlKqmvTk6BdGTCy/t/9Cwf/sku1rbW8sCHXsbbLl6ekeG3dW2hTqNJX6BgIwII1QliRQRD49O87ban2HF8hJv/bAs15aXsOpldR3D3zm5WN1ezeakjo89bVVZKVVlJwaWGXB4fzip7xoZJIY2OQESaRcQZ/r0SuAo4MOuY5cDdwDuNMYfSZUs0lmZ/LtcJFrqiMlUkExGMTvp49Tce56G9fXz06nO5472XsqKxOkMWhrqELOXJQnYEG5Y4eGFwHI93ZgPXqREPf3rLHzgyMM5t79rGGy7oYNNSR8z1j5nCWkf5pq1LM/rFZVGIMhPZUNVN5ylnO/BbEXkO+COhGsG9InKTiNwUPuZfgEbgv0Rkl4hsT6M9wMwWr1yuE/S5pxa0hyBVOKtDNYJEe2C///tj9IxO8aO/uoT3vWLNvHoy6cAaLCu0qeJoNnY4CBrY3xtKDx0ZGOMtt/yBwfFp/ucvX8Qr1rUAoWUvB3rHmPJlRiNqNg/u7QVC6yizQWNNWcFFBJnWGYL0dg09B1wQ4/pbon5/D/CedNkQC2dV6IzyUC5HBO4pNmU4zIZQROANBJnwBmLmJ91TPr7zxFFeeV4rF69qyLh9Fue21fLEkaGC0xmKxhq229PtpsRm44bvPUOJzcYd772U9e0zg3hbljnxBw17e0a5cEXm/0+e73bT4axc0HaxVNJUU87JYU9WnjtduDw+li5gbWcqKCqJCYvOltzVHLJWVGYjNRTRG4qTHvrB74/hnvLzoSs7M2nWHKyCcSGnhtrqKmisLuOunad4+21PUVNRyl1/faYTgJAjAHj2RHbSQ3t7RiOyGNmgqabwhOdcHl/GI4KidARrW2s50j+WMAWSLawVlZmcIbBINF08Pu3n9ie6uHJdS0QYLVtYXzzZSJ9lChFhQ4eD506Nsqy+ip/e9OKYtZiWugqWOCqyUieYmPbTNTTBhiXZez801ZQzPOElkMNdgAtlxOPNeLRbVFpDFmtaapjwBugZncpaSBuPxayoTBVWO2aszqEfPHmM0UkfH7oqu9EAhAqpP77xErYtULQu33jrRctwVNr53Os34khwhrhluZPdpzLvCA70uTFmZrNaNmiqKSdoQu/Zppr8PzGY9AaY9gczPjFftBEBkJO7CRa7ojIVWKmW2Y5gYtrP7Y8f5RXnNrNpqTPjdsXiktWNGVvsnS2uO7+d/3zbBQmdAITSQyeHJzmd4RTJvvCcQzZTQ4U2S2B99pyVhdM1lLNEOodSXDD+8TMn+Ny9+84q5dSfBXkJC0ssbHjizOni/3nqOCMeHx/Mcm1Aic3msHPOdHpob4+b+io77VmIXi2sKOBsOoeCQcN3n+iKpGWzSTZ0hqBIHUF9dRlNNWUpbSE9MjDOJ3+xh9uf6OI7T3Qt+nEyvaIymrqKsPBcVI3A4/Vz22NHefnaZi5YXtipmHzl/KUOSmzC7iw4gvOW1GVlfsDCcgRnExH8/oUhPnvvPr7z+OI/t6liRl5CI4KM0NlSm7LOIWMM//zz56m0l3DZ2ma++MABdp5Y3LSntZksGx8uW1h4Ljo19MOnTnB6wpv1TiElPlVlpaxtreXZDDoCXyDIwf6xrBaKIdQ1BDB0FhHBT3ecAuD+Pb1ZbyCxJKgzLZ9SvI6gtYYj/eMp+Y//2bPdPHV0mI9du45vvO0C2p0VvP+HOxe1d9baVZwt6qvsEUcw6Q3w7cde4KVrmha8VlLJLFuWOdh90pUxDa0XBsfx+oNZrQ9AaBd5qU0WHRGMTvp4cE8fLbXlnByejKvvlCm0RpBhOltrGZv2RwTeFovL4+Xz9+3nguVO3nbRchyVdr719q0Mjk/zkTt3L/iD2Z/hFZWzaaiekZn44dPHGRr35kSnkJKYLcucuKf8dJ3OjBLp3u7sF4oh1GbbWFPG0NjiHMF9YcG8L795EyU24YE9vSm2cGFYayq1ayhDWAXjs50w/tKDB3BN+vj8G87HFpZb2LTUySeuW89vDgxw2+NHk34sa0VlWxb7451VZYxM+JjyBfj2Y0e5dHUjF63M3hSxkhxbloUitkzVCfb1uqmw21jVVJOR50tEU005pxe5We+nO06ytrWGy9Y2c8nqBh54vi+r6aGRCS+V9szv2ShaR5CKFtIdx4f5v2dO8hcvWTmnl/r6F6/kuvPb+PJDB9l+bDipx8vGisrZNIRrBP/3zAkGx6Y1GsgT1rTUUF1WkrHOob09o6xrq8uK1tRsFis8d2RgnJ0nXLz5wpBg3rUb2zk6NJFV+ZmRLEwVQxE7gobqMhqryxbdQuoLBPnEz/awxFHBh69aO+d2EeGLb9rE0vpKPvB/z86r6ume8vHtx14AMruicjaWFPUtv3uBF61qSLghS8kdSmzC+RlSIjXGsK/HnfW0kMVihefu2nmKEpvwhi0hwbyrN7QhAvc/n730kMvjxZEF6ZSidQQQOotabAvp937fxYG+MT79ug1Ux1kgUVcRqhecHvfyd3fsilkvGPX4+Pojh3jpF3/Dfz36Aletb+GKsLJkNmiotuMLGPrd09oplGdsWVbP/l532pVIT41M4p7yZ3WiOJpVjdX0jE5yoC/5Qm8gaLh75ykuW9sc2X/dXFvORSsbsloncE1qRJBxQmsrF945dGrEw9cfOcxV61t51YbEaxk3djj45GvW8+jBQW4Jn/FDSM/nKw8d4CVf+g3/8evDXLK6kXs/8FJuv/4iqsqyp/xh9S9ftLKeS8/RaCCf2LLMiS9g2Neb3s6XvZGJ4uy2jlq889IV1JaX8vn79id9n8cPD9LvnuYtF565Ve26jW0c6h/nSJZEKUMS1BoRZJTO1hrGpvz0uxeWX/z0PftC/77uvKSOf8clK3jNpna+9vAhHtzTxxce2M9LvxSKAC5b28wDH3oZt75rW9bF3ABWNFQhAn971dqsDgopC8dSIt2VZiXSfT2jlNiEdWEV2GzjrCrjg1d28vjhIX53KLmd5j/dcQpnlZ0r1p8ZfV+zsR2AB7MUFbg8vnklRdJBcTuClnDBeAHpoYf39vGr/f387Ss7k9YMFxG+8MbzWd5QxU3/u4PbHgtp+j/84ZfzrT/fOkdaOJtcvKqBpz5+JS9e05RtU5QF0uaooK0u/Uqke3vcnNNcnfHOlkS889IVLG+o4t/u2z+vEumox8fD+/p5/eYllJee+Te0OSrYutzJ/c/3pdPcmASDBlcWltJAsTuC1oW1kE5M+/n0PXtZ11bLDS9ZtaDnqq2wc9u7LuTGl6/mkb+7jP946wV0tubGGVU0IpLVriXl7NiyLP1KpPt63ZyXQycvAOWlJfzjtes42D/GHdtPJjz2nud68PqDvGXbspi3X3d+O/t63RzP0EyGxdi0n6DJzp6NonYETTXlNFSXcSTJiODmXx2iZ3SKz//JRuyLUL5c01LLP123nnOas997rRQmm5c5OX7aM2+X2mIZnvDSOzqVM/WBaK7d2Ma2FfV87eFDjE/74x730x2nWNdWG7fr6epw3e+BPZmNCrKlMwRF7ggg3DmURERwctjD9588xp9tW5aVlYCKkgxWnSBdg2V7e0aB7E8Ux0JE+MSr1zM0Ps23f/dCzGMO94+x++TM7EAsljVUsWmpI+OOwNIZclZqaijjdLbUcCiJbWXf+u0RRIQPv1JbKpXcZdNSBzZJnyS1tYMgV1pHZ3PB8npeu3kJtz1+lN7RyTm3/3THKUptwhsu6Ej4ONdubGf3SRenRjK3D9nSGcq04ByoI2Btay3uKT+DCbRKjp+e4M4dp3j7xctpd+TWRjNFiaa6vJTOltq0OYK9PaFl9dlIXyTLP1x9LsEgfOWhg2dc7w8EufvZbi4/t2XebWbXbgylhx7MYFQwakUEmhrKPMloDn3j10cotQl/c/k5mTJLURaNVTBOh2bO3p7RnI0GLJY1VHHDS1Zy985u9nSPRq5//PAQg2PTvHnW7EAsVjZVs769LmlHMOULRHL8iyUSEagjyDxW5068FtKjg+P87NlTvPOSFZEJREXJZbYsd+Ly+Dh2OrVpDY/Xz9GhiZzrGIrF37xiDfVVdj5338zGwDt3nKShuizpyf3rNrax/fgIfaOJFYp7Rye57huPc/XNjzE25Ut4bCKsGkFdReYHSoveETTVlOGssseNCP7j14cpLy3hJo0GlDwhXQXjA31jGJObheLZOCrtfPiqtTx1dJhf7R9gZMLLr/YN8PotSygrTe5r79rzQ+mhh/bGjwpOnPbwllv+QP/oFP3uab75myOLttnl8VJXUZqVXdxF7whEhLUttTFbSA/3j3HP7h6uf/HKeXOKipIrdLbUUGlPvRJpRFoiBybgk+HtL1rO6uZqvnD/fu7aeQpvIJhUWshiTUstnS01cbWHjgyM85ZvP8n4tJ//u/ES/nTbUr77+y5eGFycPIXL46O+Oju1l6J3BABrWms4FGNb2c2/PkyVvYQbX746S5YpysIpLbFx/lJHyldX7usZxVllZ0kW1XEXgr3ExsevXc/RoQm+/NBBzmuvW/D8w7Ub23ima3hOM8nenlH+7Nt/IBCEH994CZuWOvno1euoKC3hs7/ct6j6zIjHm7UivDoCQmdQo5M+BqM0zff3urnvuV7+4qWraMiSl1aUxXLBMif7e9xM+1OnRLqvJzRRnE8aVFetb+GS1Q14/QuLBiyuPb+doIGH982kh549McLbbn2K8lIbd7z3Eta1hVJlzbXlfOiqTn53aJBf7x9Y8HO5PL6szBCAOgJgZknNkag6wc2/OkRteSnvealGA0r+sXmZE28gyP7exS9eisYfCHKgbywv6gPRiAifff1Grlrfwhu3Jp4diMW6tlpWNVVHuof+8MJp3nH709RXl3HHTZeyepZKwPUvXsmalho+e+++BcuBuyazozME6giA6BbS0IdmT/coD+3t5y9ftiorSoCKcrbMKJGOpOTxXhicYNofzElpiflY21rL7ddftKi0i4hwzcY2nnzhND9/tpt3f+8ZljgrufO9l8YUnbSX2PjUa8/jxLCH7zzRtaDnck34NDWUTZpry3FU2jkc1iD/+iOHcFTa+YuXLkxYTlFyhXZHBS215SkrGO/rDfXj5/oMQTq4bmM7gaDhwz/ZxZqWGn7y3ksTtpK/rLOZqze08s3fHIk53RwLXyDI2LQ/KzMEoI4ACHn9zpYaDg+Ms+uki18fGODGl6+mrkKjASU/ERE2L3Oy+9To/Acnwd5uN+WlNlY3Vafk8fKJjR11rG+v46KV9fzory5Jqmb4z68+j6AxfOH+A0k9hysyVaypoazS2VrL4f4x/v2RQ9RX2bn+xSuzbZKinBVbljnpGpo464lXCLWOrmuvy0qPe7YREX7+vhdzx3svxZFkMXdZQxXvvewc7tndwzNdw/Mev+N46Bh1BFmms6WGEY+Pxw4NctNl51ATZw+xouQLF4TrBE8n8UWUCGNMSFoiDyaK00V5acmCu6X++rJzWOKo4FP37I27LOdw/xjv+cEfuel/d9LuqOCS1dlZD6uOIIy1pKappox3Xroiy9YoytmzdUU9S+sr+ewv951VVNDtCi2rz7eOoWxTWVbCJ159Hvt73fzomRNn3DYwNsXH736eq29+jKePDvMP15zLb//+8qwthVJHEOa89jrKSm184IrOrC6PV5RUUWEv4Ztv38rA2BR/f+dzixahm1lWr45goVx3fhuXrm7kaw8fZGTCy8S0n68/cojLv/Iod24/ybsuXcmjH72cv7l8TVZXf+o3XpjGmnL++E9XUVepL4lSOGxZ5uTj167ns/fu4/bHu/irRUzJ7+1xYxMig1NK8ogIn3rdebz6G0/w1z/cwQuDEwyOTfPq89v56NXnsjJHiu/6rReFzgwohcgNL1nJM13DfOnBA2xdUc+FK+oXdP99PW5WN9dQWZY7y+rziXVtdbzzkhV8/8ljXLSynm+/80K2Ll/Y/0G60dSQohQ4IsKX3ryJdmcFH/jRTkYWuM94X8+opoXOkk+8ej33vP8l3PHeS3POCYA6AkUpChyVdr719q0MjXv5yJ27CcbpYpnNyISXntEpdQRnib3ExqalzpzVaUqbIxCRChF5RkR2i8heEflMjGNERL4hIkdE5DkR2ZouexSl2Nm01MknXr2e3xwY4LbHjyZ1n3294R3F7fknLaEkTzojgmngCmPMZmALcI2IXDLrmGuBzvDPjcD/S6M9ilL0vOvSFVx3fhtffugg24/NP1+wtyc0mawRQWGTNkdgQlhynvbwz+x49PXAf4ePfQpwikh7umxSlGJHRPjimzaxtL6S9//oWYZj1AtOj09z33O9fOJnz3Pb410scVRkbWGKkhnS2jUkIiXADmAN8C1jzNOzDukATkZdPhW+7oyVQCJyI6GIgeXLl6fNXkUpBuoqQvWCN/7Xk/zdHbv4xtsuYPuxYX5/5DRPvnCaF0d+zQAABRRJREFU/eF0UE15KS9a1cBbL9bPXKGTVkdgjAkAW0TECfxMRDYaY/ZEHRKrcjKnimWMuRW4FWDbtm2Lm4pRFCXCxg4Hn3zNej75i71s/szDGANlpTa2rajn71+1lhevaWJTh6MotYWKkYzMERhjXCLyKHANEO0ITgHLoi4vBXoyYZOiFDvvuGQFo5M+Jn0BXnJOE1tX1Gd1ulXJHmlzBCLSDPjCTqASuAr40qzD7gHeLyI/Bl4EjBpjYm+KVhQlpYgI77+iM9tmKDlAOiOCduAH4TqBDbjDGHOviNwEYIy5BbgfuA44AniAG9Joj6IoihKDtDkCY8xzwAUxrr8l6ncDvC9dNiiKoijzo5UgRVGUIkcdgaIoSpGjjkBRFKXIUUegKIpS5KgjUBRFKXLUESiKohQ5stg9ptlCRAaB44u8exMwlEJz0kU+2Kk2pga1MTWojfOzwhjTHOuGvHMEZ4OIbDfGbMu2HfORD3aqjalBbUwNauPZoakhRVGUIkcdgaIoSpFTbI7g1mwbkCT5YKfamBrUxtSgNp4FRVUjUBRFUeZSbBGBoiiKMgt1BIqiKEVO0TgCEblGRA6KyBER+cds2xMLETkmIs+LyC4R2Z5tewBE5LsiMiAie6KuaxCRR0TkcPjf+mzaGLYplp2fFpHu8Ou5S0Suy6J9y0TktyKyX0T2isiHwtfnzGuZwMaceR3D9lSIyDMisjts52fC1+fSaxnPxpx6LS2KokYQXo5zCHglofWYfwTeZozZl1XDZiEix4BtxpicGYwRkZcD48B/G2M2hq/7MjBsjPli2KnWG2M+loN2fhoYN8Z8NZu2hW1pB9qNMTtFpBbYAbwBeDc58lomsPFPyZHXEUBEBKg2xoyLiB14AvgQ8EZy57WMZ+M15NBraVEsEcHFwBFjzFFjjBf4MfD6LNuUFxhjHgOGZ139euAH4d9/QOjLIqvEsTNnMMb0GmN2hn8fA/YDHeTQa5nAxpzChBgPX7SHfwy59VrGszEnKRZH0AGcjLp8ihx8gxN6ozwsIjtE5MZsG5OAVmu3dPjflizbk4j3i8hz4dRR1lNYACKyktD2vqfJ0ddylo2QY6+jiJSIyC5gAHjEGJNzr2UcGyHHXksoHkcgMa7LRe/8EmPMVuBa4H3hdIeyeP4fcA6wBegFvpZdc0BEaoC7gA8bY9zZticWMWzMudfRGBMwxmwBlgIXi8jGbNs0mzg25txrCcXjCE4By6IuLwV6smRLXIwxPeF/B4CfEUpp5SL94XyylVceyLI9MTHG9Ic/jEHgNrL8eoZzxXcBPzTG3B2+Oqdey1g25trrGI0xxgU8Sij3nlOvpUW0jbn6WhaLI/gj0Ckiq0SkDHgrcE+WbToDEakOF+gQkWrgVcCexPfKGvcA14d/vx74RRZtiYv1pRDmT8ji6xkuHn4H2G+M+feom3LmtYxnYy69jgAi0iwizvDvlcBVwAFy67WMaWOuvZYWRdE1BBBu07oZKAG+a4z5fJZNOgMRWU0oCgAoBX6UCzaKyP8BlxOS0O0HPgX8HLgDWA6cAN5ijMlqoTaOnZcTCsENcAx4r5VDzoJ9LwUeB54HguGr/4lQDj4nXssENr6NHHkdAURkE6FicAmhk9k7jDGfFZFGcue1jGfj/5BDr6VF0TgCRVEUJTbFkhpSFEVR4qCOQFEUpchRR6AoilLkqCNQFEUpctQRKIqiFDnqCBRFUYocdQSKoihFzv8HfGpkowVlqogAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(df.index,df['Response'])\n",
    "plt.title('kelp upper limit')\n",
    "plt.ylabel('depth [m]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "a) Become familiar with data set, make plot each variable vs. the response (kelp) to see which variables have strongest correlation with response (kelp)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "# insert code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Goal of analysis:\n",
    "\n",
    "* Instead of correlating with one variable, we will use PCA to create an \"index\" of common environmental variability.\n",
    "\n",
    "#### Key concept (again)\n",
    "\n",
    "The final exercise in this notebook is:\n",
    "\n",
    "Which principal component is most highly correlated with kelp extent? What is your interpretation of the relationship between kelp extent and environmental variability?\n",
    "\n",
    "### PCA Step 1: data matrix\n",
    "\n",
    "Put our environmental variables into matrix form. This matrix (which excludes the response variable) has four columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# insert code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Standardize data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "b) After creating an array for the environmental data:\n",
    "\n",
    "Standardize the environmental data by transforming to z-scores.\n",
    "\n",
    "For each column, subtract the mean and divide by the standard deviation. This should transform each column into a standardized variable with a mean of 0 and a standard deviation of 1.\n",
    "\n",
    "*Hint:* Use the `ddof=1` option with `np.std`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# insert code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PCA Step 2: correlation (or covariance) matrix\n",
    "Next create the correlation matrix. We can also think of this as the covariance matrix of the standardized variables. If we don't standardize the variables, we end up weighting them differently simply based on how big the numbers are."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# insert code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PCA Step 3: extract eigenvalues and eigenvectors\n",
    "\n",
    "$$RV = V\\Lambda$$\n",
    "\n",
    "The eigenvalues are in the matrix $\\Lambda$ \n",
    "\n",
    "The eigenvectors are in the matrix $V$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# insert code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "c) After extracting the eigenvectors and eigenvalues of the correlation matrix:\n",
    "\n",
    "Find the percentage of total variance explained by each principal component."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# insert code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Factor loading matrix\n",
    "\n",
    "Eigenvectors times square root of eigenvectors\n",
    "\n",
    "$$A = V\\sqrt{\\Lambda}$$\n",
    "\n",
    "The columns of this new A matrix are called the factor loadings for each principal component. We can think of new axes for the data in a space defined by waves, wind and tides. The factor loadings define those axes - both the orientation (eigenvector) and the spread of the data across each axis (square root of eigenvalue - which describes the standard deviation)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# insert code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A common way of visualizing the results of PCA is to plot the factor loadings for the first two principal components (the first two columns of the A matrix). If variables are clumped together, it means that they covary (similar direction with respect to origin), with similar magnitudes (similar distance from origin). \n",
    "\n",
    "This plot can be thought of as viewing the original data axes (which exist in four-dimensional space) from a perspective that is perpendicular to the plane defined by the first two PC's. It is the _projection_ of the original axes onto the 2D plane defined by the first two PCs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAEKCAYAAAAmfuNnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxU9b3/8dcnE8KegOz7orgACsWIIG2FWlsXvFRvF7letdaWimArUJWKdb322rrhzkXr0k2rv2ql1latiNpCWxZRFrGiIoagiEDCDkk+vz/OSRhCJkzIzJlJ8n4+HvOYmbN+Mgx553vO93yPuTsiIiLplpPpAkREpGlQ4IiISCQUOCIiEgkFjoiIREKBIyIikVDgiIhIJDIaOGb2sJltMLPlCeabmd1tZqvN7C0zGxY37zQzeyecNz26qkVE5FBkuoXzKHBaLfNPBwaEjwnAAwBmFgPuC+cPBMab2cC0VioiIvWS0cBx99eATbUsMg74pQf+AbQzs27AcGC1u7/v7nuAJ8JlRUQkS+VmuoCD6AF8FPe+KJxW0/QTa9qAmU0gaB3RunXr448++uj0VCoi0kgtXrx4o7t3qu92sj1wrIZpXsv0Aye6zwZmAxQWFvqiRYtSV52ISBNgZh+mYjvZHjhFQK+49z2BYiAvwXQREclSme40cDBzgAvC3mojgBJ3Xw8sBAaYWT8zywPODZcVEZEsldEWjpk9DowGOppZEXAd0AzA3WcBzwNnAKuBHcBF4bwyM5sMvADEgIfdfUXkP4CIiCQto4Hj7uMPMt+BSQnmPU8QSCIi0gBk+yE1ERFpJBQ4IiISCQWOiIhEQoEjIiKRUOCIiEgkFDgiIhIJBY6IiERCgSMiIpFQ4IiISCQUOCIiEgkFjoiIREKBIyIikVDgiIhIJBQ4IiISCQWOiIhEQoEjIiKRUOCIiEgkFDgiIhKJjAaOmZ1mZu+Y2Wozm17D/CvMbGn4WG5m5WZ2WDhvjZktC+ctir56ERGpi9xM7djMYsB9wKlAEbDQzOa4+8rKZdz9VuDWcPmzgCnuviluM2PcfWOEZYuIyCHKZAtnOLDa3d939z3AE8C4WpYfDzweSWUiIpJymQycHsBHce+LwmkHMLNWwGnA7+MmO/CimS02swlpq1JERFIiY4fUAKthmidY9izg79UOp41y92Iz6wy8ZGar3P21A3YShNEEgN69e9e3ZhEROUSZbOEUAb3i3vcEihMsey7VDqe5e3H4vAF4huAQ3QHcfba7F7p7YadOnepdtIiIHJpMBs5CYICZ9TOzPIJQmVN9ITMrAE4Gno2b1trM2la+Br4CLI+kahEROSQZO6Tm7mVmNhl4AYgBD7v7CjO7JJw/K1z0bOBFd98et3oX4Bkzg+Bn+K27/yW66kVEpK7MPdFpk8ansLDQFy3SJTsiInVhZovdvbC+29FIAyIiEgkFjoiIREKBIyIikVDgiIhIJBQ4IiISCQWOiIhEQoEjIiKRUOCIiEgkFDgiIhIJBY6IiERCgSMiIpFQ4IiISCQUOCIiEgkFjoiIREKBIyIikVDgiIhIJBQ4IiISCQWOiIhEQoEjIiKRyGjgmNlpZvaOma02s+k1zB9tZiVmtjR8XJvsuiIikl1yM7VjM4sB9wGnAkXAQjOb4+4rqy36uruPPcR1RUQkS2SyhTMcWO3u77v7HuAJYFwE64qISAZkMnB6AB/FvS8Kp1U30szeNLM/m9mgOq6LmU0ws0VmtujTTz9NRd0iInIIMhk4VsM0r/Z+CdDH3YcA9wB/qMO6wUT32e5e6O6FnTp1OuRiRUSkfjIZOEVAr7j3PYHi+AXcvdTdt4WvnweamVnHZNYVEZHsksnAWQgMMLN+ZpYHnAvMiV/AzLqamYWvhxPU+1ky64qISHbJWC81dy8zs8nAC0AMeNjdV5jZJeH8WcDXgYlmVgbsBM51dwdqXDcjP4iIiCTFgt/fTUNhYaEvWrQo02WIiDQoZrbY3Qvrux2NNCAiIpFQ4IiISCQUOCIiEgkFjoiIREKBIyIikVDgiIhIJBQ4IiISCQWOiIhEQoEjIiKRUOCIiEgkFDgiIhIJBY6IiERCgSMiIpFQ4IiISCQUOCIiEgkFjoiIREKBIyIikVDgiIhIJDIaOGZ2mpm9Y2arzWx6DfPPM7O3wsd8MxsSN2+NmS0zs6VmpvtGi4hkudxM7djMYsB9wKlAEbDQzOa4+8q4xT4ATnb3zWZ2OjAbODFu/hh33xhZ0SIicsgy2cIZDqx29/fdfQ/wBDAufgF3n+/um8O3/wB6RlyjiIikyEFbOGY2tYbJJcBid19aj333AD6Ke1/E/q2X6i4G/hz33oEXzcyB/3P32TWtZGYTgAkAvXv3rke5IiJSH8kcUisMH38M358JLAQuMbOn3P3nh7hvq2Ga17ig2RiCwPl83ORR7l5sZp2Bl8xslbu/dsAGgyCaDVBYWFjj9kVEJP2SOaTWARjm7tPcfRpB+HQCvgh8ux77LgJ6xb3vCRRXX8jMjgMeAsa5+2eV0929OHzeADxDcIhORESyVDKB0xvYE/d+L9DH3XcCu+ux74XAADPrZ2Z5wLnAnPgFzKw38DRwvrv/O256azNrW/ka+AqwvB61iIhImiVzSO23wD/M7Nnw/VnA4+Ev+pWJV6udu5eZ2WTgBSAGPOzuK8zsknD+LOBaghbW/WYGUObuhUAX4JlwWi7wW3f/y6HWIiIi6WfuBz+tYWaFwCiC8y5/c/cGed1LYWGhL1rUIEsXEckYM1sc/rFfL8leh/MGwfmV3HDnvd19bX13LiIiTUcy3aIvA64DPgHKCVo5DhyX3tJERKQxSaaF80PgqPgeYiIiInWVTC+1jwgu9BQRETlkybRw3gfmmdmfiOsG7e53pK0qERFpdJIJnLXhIy98iIiI1NlBA8fdb4iiEBERadwSBo6ZzXT3y83sj9Qwxpm7/0daKxMRkUalthbOr8Ln26IoREREGreEgePui8PnV6MrR0REGqvaDqktI8HtAgDcXRd+iohI0mo7pDY2fJ4UPlceYjsP2JG2ikREpFGq7ZDahwBmNsrdR8XNmm5mfwduTHdxIiLSeCQz0kBrM6u606aZnQS0Tl9JIiLSGCVz4efFwMNmVhC+3wJ8J30liYhIY5TMhZ+LgSFmlk9w/xyNqyYiInWW1P1wzOxMYBDQIrzLJu6uczgiIpK0g57DMbNZwLeAywjuhfMNoE+a6xIRkUYmmU4DJ7n7BcDmcFy1kUCvVOzczE4zs3fMbLWZTa9hvpnZ3eH8t8xsWLLriohIdkkmcHaGzzvMrDuwF+hX3x2bWQy4DzgdGAiMN7OB1RY7HRgQPiYAD9RhXRERySLJnMN5zszaAbcCSwhGH3gwBfseDqx29/cBzOwJYBywMm6ZccAv3d2Bf5hZOzPrBvRNYl2RrOBlZZSXlFC+ZQvln22i/JMPKVv/AaVFH7Jj/Xo2b17Php2ljHpoDp269c90uSJpk0wvtZvCl783s+eAFinqqdaD4G6ilYqAE5NYpkeS6wJgZhMIWkf07t27fhVLk1exc2cQHFu2UL55M+VbtlBW+X7LFso3baJ848eUb9pI+ZYSyrdup2Ln3oTbsxynbXNna5scrnpxCneP/w1t8tpE+BOJROeggWNmzYCJwBfDSfPM7P/cPfH/ouRYDdOqj92WaJlk1g0mus8GZgMUFhYmHBtOmhavqKBi69Z9QbFlC2VhgOx7lOwXLOVbtuC7dyfcZk4zJ5ZXTqx5BbG8CvJaVpBT4OzOy6MkrzWf5eXzWfP27DmsB7Hu/Sk4fCC9jz6OAb27s6F4Hm+8+iMu+eslzPryLIWONErJHFJ7AGgG3B++Pz+c9t167ruI/Tsf9ASKk1wmL4l1pYnwvXv3D46qkCipFiBx4VFSAhUVNW8wx4i1bkmsdTNiLYxmeeW06LCbWIdtxHJ3VQVKbvMKrEWMnW07sqFlJ94v78CqnQUUVXRkHR0pyetCh279OKZnBwZ1L2BQ93xO7tSGWM6Bfy+d0ucUbj35Vq549Qom/nUis06dRetmGtBDGhcLTo/UsoDZm+4+5GDT6rxjs1zg38ApwDpgIfBf7r4ibpkzgcnAGQSHzO529+HJrFuTwsJCX7RoUX3KljRyd3zHjgMPU23ecmBwxAVIxfbtCbdpzZsTa9cueBTkE2udR6ylEcurIBbbTSxnGzEvIbf8U2JlG4g120NOM8cqM6FFO2jXi92te/BprBMflnVg5Y4CFm1uzeKStmwkHzA6t23O4B5BqASPAnq2b4lZTY3xxF768CWuePUKju14rEJHsoaZLXb3wvpuJ5kWTrmZHe7u74U77g+U13fH7l5mZpOBF4AY8LC7rzCzS8L5s4DnCcJmNcEI1RfVtm59a6qvKVOm0KdPHy6//HIAvvrVr9KrVy8eeughAKZNm0aPHj2YOnVqJsuMhFdU7DtRHn+I6oDDVvuHh+9NfKQ2p23bfeHRvj15/fvFhUk7clvlEmu2h1jOdmJeSqxiIzk718OWj6BkOWz/dP8NWg607QYFvaDdCCjoiRf0YmOsM+/sLGBJaRve+LiMFcWlbFiz71Banw6tGNQ7n4u6F1SFS6e2zVPyuZ3a59T9WjoPfPkBhY40Gsm0cE4BHgHeJzh30ge4yN1fSX95qZXuFs5TTz3FU089xZNPPklFRQUnnHACeXl5LFiwAICRI0cyc+ZMTjyxxv4NWatiz57ErYxEAVJSAom+W7HYvqCoCpB25FafFhcusTatsV0bw/AogpK1ca/D5z3b9t9Pbkso6AntegXPBb3jXveivE033t+0mxXFpSxfV8KK4lJWFJdQuqssKDPHOKJTmyBUwtbLwO755LdoluZPHF5c8yJXvnYlQzoN4f4v36/QkYxKVQvnoIET7qw5cBRB4Kxy98RnTrNYugOnuLiY4cOHU1RUxLJly7jttttYv349v/vd72jVqhVdunRhw4YN5OXlpa2G2rg7Fdt3UL5lc80BkiA8KnYkvv2RtWwZFw4FVa9z27dPGB45bdoceKhpz459QVJSFIbJR/tel64Dr9awbnlYGCDhIy5MaNcbWnWg8tjYrr3l/PuTrVWhsqK4lLfXl7Jrb3Aep3luDkd3bVsVLIO6F3B017a0aBZL6b9BXbyw5gWueu0qhnQawgNffoBWzVplrBZp2tJ+SM3Mzkkw63Azw92fru/OG5vu3buTm5vL2rVrmT9/PiNHjmTdunUsWLCAgoICjjvuuJSFjZeXU15aun9I7Bcim/dNK6k8J1ICtR2yKijYFxqdOtJ8wBH7WhgJWh85LVokUazDjk1BmBS9GRcmH+17veOz/dexGOR3D8Kj94iagyWv5r/6t+0uY2VxKSuK17B8XRAwqzdso6wi+OOqbfNcBnbP57+G92FwjyBcDu/UmtxYMtdBR+erfb+K40x/bXrV4TWFjjRktZ3DOauWeQ4ocGowatQo5s+fz/z585k6dSrr1q1j/vz5FBQUcNJJJ9W4TsXu3YlbGJsre17FzyuhorQ08SGr3Nx9h6gK2pHXty8tq0KifdUhrP3CIz8fy01qLNcDlZfB1uLaD3ftrdZKatYqDJCe0H1otTDpFZxbiR28ns+2hYfEwlbLyuJSPti4rxNBxzbNGdQ9n1OO6VzVU6xX+1bk1NBTLBud1vc0cJj++nQufflS7j/lfoWONFi13fHzoigLaQzcnRHDhvH6Cy/w5r/+Rd9NmyjYvJlbnnmGNrEY5x57LOumTg2v99h3Qt137ky4TWvVar9zG3k9eh54fiPuPEisXTtyWreuc++oWu3ZXnOrpPJw19Zi8GpdjFt1DMKj01FwxKkHHu5q2R7qUKO7s27LzvCQWCkri0tYvq6Uj0t3VS3Ts31LBncv4JzP9WBQj3wGdy+gc34SLbAsd1q/0wC46vWrFDrSoB3in7QSLxaLcWSLlnh5OXsqytlWUcHhzZuz8LvfZewHHwBQ5s6a997jk379+cbgwTTr0oUWRx2VODzCAMlJ9/ked9i+seZWyZa1weudm/dfJyd33+Guvp8/MEzye0Deof9CLK9wPti4vepcS+Xzlh3B4cAcg8M7tWFE/8OCVkuPfAZ1K6CgVfpP5mfKaf1OCw6vvT6dSS9P4r5T7lPoSIOjwEmBli1aMHfaNGLt2vHK2rWcd8/dXHj22fQYP57Dp07l+OOP55+LF/Pcn//MOeecQ8evfoWLLoqoAVm+NzjhHh8mW9buHyxlu/ZfJ6/NvsNdPQtrONzVFXJSczJ9T1lFeDK/pKr18vb6UnbsCToI5MVyOKprW04b1LXqhP4xXfNpmZe5k/mZcnq/03F3fvy3HzN57mTu/dK9Ch1pUJLqpdZYpKuXWps2bdi2LeiS+9RTT/Gb3/yGP/zhD6xZs4axY8eyfPnyqmXnzp3LtGnTeOONN1Kz891bazjcVbTv9db1HDDqT+vOcd2Fq4VJQc86H+5K1vbdZby9vrSq1bJ8XSnvbtjK3vKgvtZ5MQZ1L2BgePHk4B4FHNG5Dc2y7GR+pv3p/T9x9d+uprBLIfeeci8tc1tmuiRp5CK58DO8rXSnyos+46Yf5+5v1XfnjcXOnTsZOnQou3btYv369cydOzfhssOGDWPVqlXJbbiiIrhYMdG5k5K1sKvaOKo5zaCgRxAe/Ufvf7iroFcwr1n6f0Ft3r5nX7CEzx9s3F7Vz6FD6zwGds/ni0f2r+op1uewhnMyP5PO7H8mjjPjbzOY/PJkhY40GLV1i/4mMBPYEA7g+W13XxjOfhQYlmjdpqZly5YsXboUgAULFnDBBRfs16qJt1+LsmwPlBYd2CqpDJiSdVBe7ZKnvLb7WiO9T9z/3ElBT2jTJWWHu5Lh7nxcuquq+3FlT7F1W/Z1hOjRriUDu+czbkiP8CLKfLrmt0htx4YmZmz/sbg71/z9Gi57+TLuOeUehY5kvdpaOFcDx7v7ejMbDvzKzK4Or7/Rb4oERo4cycaNG/n0009hdymU7YZ3/lx1Ev6Nvy/kmM7N4LajYNsnHHC4q03XIDi6DYGjxx54QWPLdhn5uQAqKpw1n22vOtdSGTCbtu8BgqNw/Tq25vg+7blgZJ+qbsjtW2fmQtfG7qzDgysXZvxtBpfNvYx7vqTQkexWW+DE3H09gLv/y8zGENyMrScJbgXQ5Gz7FDZ/ABXl8LeZUPIRq95eSfn2TXR4cBg7Pi2BTTvg8XMBWFMa40dP7eKyrx4JR4yqdh6lZ/DITc2YXPW1t7yCdz/Ztl9PsZXFpWwPT+Y3ixlHdmnLl8PrWwb3yOforvm0bq5+KFGKD50fzP0B93zpHlrkNvyu4NI41fbbYWv8oJ1hS2c08AdgUBTFZb0XZ8Bbv2Pnrl0M/fo0yInhOc14bOJIYkOHw85WvDfrZ3zu/3Vl194K2ha057IbJ0bXQy1JO/eU8/bHpaxYt6+n2Dsfb2VPeXBtTau8GMd0y+frx/esOql/ZJe25OXqZH42OOvws3Cca/52TVVLR6Ej2ShhLzUzGwJsd/fV1aY3A77p7r+JoL6USnkvtaJFwZAslYe7WuSnbttpUrJj736tluXFpbz/6TbCUV9o16oZg7vvG6hycI8C+nZoXeM9XCS7PLv6WX7y958wotsI7v7S3QodSZkoeqltB7oQ3Bog3gjgH/XdcaPQs96ff9q4Oxu27q7qflwZMkWb953M71bQgkHd8znj2G4MDkdE7l6gk/kN1bgjxuE41/79Wn74yg+5a8xdCh3JKrUFzkyCjgPV7Qzn1TbWmkTI3Vm7acd+wbKiuISN2/ZULdOvY2uG9GrHeSf2qbpJWIc22XG+SFLna0d8DXfnuvnXcfkrl3PXl+6ieUz/zpIdagucvjVda+Pui8ysb9oqaoCKioqYNGkSK1eupKKigrFjx3Lrrbcyf/58xo0bR//+/dmxYwddunThyiuvZOzYsYe8r7LyClZ/uo0V60qrBq18u7iUrbuDe7jk5hgDurRl9FGdq4bZP6ZbW9pGcA8XyQ5nDzgbgOvmX8cP5/5QoSNZo7bAqa0trr6XIXfnnHPOYeLEiTz77LOUl5czYcIEZsyYwZlnnskXvvAFnnvuOQCWLl3K1772NVq2bMkpp5xy0G3v2lvOqo+3Vh0WW1lcwqqPt7K7LDiZ36JZDsd0y2fc57qH510KGNClTUbv4SLZoTJ0rp2/7/CaQkcyrbbAWWhm33P3B+MnmtnFwOL0ltVwzJ07lxYtWlT1PIvFYtx5553069ePMWPG7Lfs0KFDufbaa7n33nsPCJzSXXtZGd55cmXYU2z1p9soD8/m57fIZVD3gv2ub+nfqY1O5ktCZw84G2ff4bWZY2YqdCSjagucy4FnzOw89gVMIZAHnF2fnZrZYcDvgL7AGoJeb5urLdML+CXQFagAZrv7XeG864HvAZU3qb/a3Z+vT02HasWKFRx//PH7TcvPz6d3796sXl29v0UwtM0tP/s5r7yzoSpgVhSXsnbTvvvFdG7bnME9CvjKoC5Vh8V6tm+pk/lSZ+cMOAd35/oF1zPllSnMHDOTvJguxJXMqO1+OJ8AJ4UXfA4OJ//J3RMPFJa86cDL7n6LmU0P319VbZkyYJq7LzGztsBiM3vJ3VeG8+9099tSUEu9uHuNQVDZ3Xzn3nL+snx91Qn9fy1awgcbt3PRI8EoQX06tGJwj3y+dUKvqnDp1FZ/hUrq/OeR/4nj3LDghqqWjkJHMqG2sdRaAJcARwDLgF+4e1mK9jsOGB2+fgyYR7XACUc5qBzpYKuZvQ30AFaSRQYNGsTvf/97yiuc9z/dxvLiEha/W8zbqz/gjn+WsOG9z7jk10uI5RhHdGpDj4qPaT54IPdPGMHA7vnk62S+RODrR34dx7lxwY1MmTeFO0ffqdCRyNV2SO0xYC/wOnA6cAzBYbZU6BI3bM56M+tc28Jhr7jPAf+MmzzZzC4AFhG0hDbXsCpmNgGYANC7d+/6Vx4q3rKTV//9Kcu3duatNRvoefYVND9mDF5RzpaX7qfniacz4ugerPywLY9OGsVRXdvy77dXMO7ux3jooYcY0b9DymoRScY3jvwGADcuuJGp86Zyx+g7FDoSqdpGGljm7seGr3OBf7l70iNEm9lfCc6/VDcDeMzd28Utu9nd2yfYThvgVeDmcOBQzKwLsJFgTLebgG7u/p2D1ZTKkQb+vGw9E3+zhLbNc+nbYif/fmYmpR9/SLMcOGvsmdxx++0sWLBgv27RnTt35sorr+Sss3QJk2TOk+88yU3/uImTe56s0JGkpGqkgdoCZ0l8wFR/X6+dmr0DjA5bN92Aee5+VA3LNQOeA15w9zsSbKsv8Jy7D65pfrxUBs623WV8tm03vdrrHi7S8FSGzuieo7l99O0KHalVqgKnttEXh5hZafjYChxX+drMSuu53znAheHrC4Fnqy9gwZn4XwBvVw+bMKQqnQ3UfPOZNGrTPJc+HVorbKRB+uZR3+SaE69hXtE8ps2bxt7yvZkuSZqAhIHj7jF3zw8fbd09N+51fUepvAU41czeBU4N32Nm3c2ssnvzKOB84EtmtjR8nBHO+7mZLTOzt4AxwJR61iPS5Hzr6G8x48QZzCuax9RXpyp0JO0SHlJrjFI+WrRII/D4qsf56T9/ypheY7j95NtpFlPPSdlfFIfURKQJGH/0eK4+8Wpe+egVfvTqj9TSkbRR4IgI448ez4+H/5i5H83liteuYG+FQkdST4EjIgD81zH/xfTh03l57ctc8apCR1JPgSMiVc475ryq0Lny1SsVOpJSChwR2c95x5zHVSdcxV/X/pWrXrtKoSMpU9vQNiLSRP33wP/GcX6+8OfwGvzsiz+jWY56r0n9KHBEpEbnDzwfd+fWRbdirxm3fPEWhY7UiwJHRBK6YNAFOM5ti26raunk5ujXhhwafXNEpFYXDgpGobpt0W3Y68YtX7hFoSOHRN8aETmoCwddiLtz++LbMYz//cL/KnSkzvSNEZGkfHvwt3GcOxbfgWH89As/VehInejbIiJJu2jwRTjOnYvvBIOffl6hI8nTN0VE6uQ7g7+DuzNzyUwM4+bP36zQkaToWyIidXbxsRfjOHctuQsIWjqxnFiGq5Jsp8ARkUPy3WO/C8BdS+7CzLh51M0KHamVAkdEDtl3j/0u7s7db9yNYfzPqP9R6EhCChwRqZfvHfc9HOeeN+7BMG4adZNCR2qkwBGReptw3ATcnXuX3ouZceNJNyp05AAZCRwzOwz4HdAXWAN8090317DcGmArUA6UVd7iNNn1RSQ63x/yfRznvqX3ASh05ACZuj3BdOBldx8AvBy+T2SMuw+tdj/tuqwvIhG5ZMglXDr0Uua8N4fr5l9HeUV5pkuSLJKpQ2rjgNHh68eAecBVEa4vImkycchEcLj/zfsxM2446QZyTLfekswFThd3Xw/g7uvNrHOC5Rx40cwc+D93n13H9TGzCcAEgN69e6fsBxCRxCYOnYjjPPDmAwAKHQHSGDhm9legaw2zZtRhM6PcvTgMlJfMbJW7v1aXOsKQmg1QWFjodVlXRA7dpUMvxXFmvTkLw7j+pOsVOk1c2gLH3b+caJ6ZfWJm3cLWSTdgQ4JtFIfPG8zsGWA48BqQ1PoiklmXDrkUIAgdM64beZ1CpwnL1L/8HODC8PWFwLPVFzCz1mbWtvI18BVgebLri0jmmRmXDrmU7x/3fZ5+92luWHADFV6R6bIkQzJ1DucW4EkzuxhYC3wDwMy6Aw+5+xlAF+AZM6us87fu/pfa1heR7GNmTBo6CceZ/dZsDOPakdeqpdMEZSRw3P0z4JQaphcDZ4Sv3weG1GV9EclOZsbkoZNxdx5c9iCAQqcJ0kgDIhIJM+Oyz10GwIPLHsTM+MmInyh0mhAFjohEpjJ0HOehZQ9hGNeMuEah00QocEQkUmbGDz73A9ydXyz/BYYxY8QMhU4ToMARkciZGT8c9kMc5+HlD2NmzDhxBmEnIWmkFDgikhFmxuXDLsdxHln+CIBCp5FT4IhIxpgZU4ZNAYdHVih0GjsFjohklJkx5fgpOM6jKx7FMK4+8WqFTiOkwBGRjDMzph4/FXfnsZWPYWb8ePiPFTqNjAJHRBeRU5MAAA0SSURBVLKCmTGtcBqO88uVv8Qwpg+frtBpRBQ4IpI1zIwfFf4Ix/nVyl9hZlx1wlUKnUZCgSMiWcXMuKLwCtydX7/9awzjyhOuVOg0AgocEck6ZkHIAPz67V8DKHQaAQWOiGSl6qFT2fJR6DRcChwRyVqVoVN1TofgHI9Cp2FS4IhIVqvsOOC+r/fatMJpCp0GSIEjIlnPLOgi7ey7Tmfq8VMVOg2MAkdEGoTKi0Hd941IMOX4KQqdBkSBIyINhlkw7I3jwdhrBlOGKXQaiozcgMLMDjOzl8zs3fC5fQ3LHGVmS+MepWZ2eTjvejNbFzfvjOh/ChHJhMpbGXzrqG/xyPJHmLlkJu6e6bIkCZlq4UwHXnb3W8xsevj+qvgF3P0dYCiAmcWAdcAzcYvc6e63RVSviGSRqpaOh/fTIbi/jlo62S1TgTMOGB2+fgyYR7XAqeYU4D13/zC9ZYlIQ5FjOcwYMQMnvHNoeCdRhU72ylTgdHH39QDuvt7MOh9k+XOBx6tNm2xmFwCLgGnuvjkNdYpIFsuxHK4ZcQ2O89CyhzCMyz53mUInS6UtcMzsr0DXGmbNqON28oD/AH4cN/kB4CbAw+fbge8kWH8CMAGgd+/eddm1iDQAOZbDT0b8BHfnwWUPAih0slTaAsfdv5xonpl9YmbdwtZNN2BDLZs6HVji7p/EbbvqtZk9CDxXSx2zgdkAhYWFOrMo0gjlWA7XjrwWgAeXPYiZMXnoZIVOlsnUIbU5wIXALeHzs7UsO55qh9Mqwyp8ezawPB1FikjDER86s9+ajWFMGjqpxtBp06YN27Ztq3p/44038vTTTwOwbNkyjj32WAC+973vsX79eh555BE6derE9u3bOe6447j55ps5+uijI/ipGhfLRHdCM+sAPAn0BtYC33D3TWbWHXjI3c8Il2sFfAT0d/eSuPV/RdCDzYE1wPfjAiihwsJCX7RoUap/HBHJIhVewQ0LbuDpd5/mkiGXMGnopAOWqR44lcrKyujYsSNbtmypmnbNNdfQsWNHLr/8cgAef/xxpk6dyvLly+nQoUP6fpAsYmaL3b2wvtvJSAvH3T8j6HlWfXoxcEbc+x3AAf+i7n5+WgsUkQYrx3K4buR1AMx6cxaGcenQS1O2/fHjx/PHP/6RJ554gkmTDgwzSUwjDYhIo1MZOu7OA28+AJDS0Bk2bBirVq1K2faaCgWOiDRKOZbD9SddjxOEjmFMHDoxJdvWyAaHRoEjIo1WjuVww0k34O7c/+b9YDBxSP1D54033uDzn/98CipsWhQ4ItKoVYUOzv1L78eoX1fpJ598kldeeYV77rknRRU2HQocEWn0YjkxbjzpRgBeX/c6O3bsoGfPnlXzp06dytSpUxOuf+utt/Loo4+yfft2jj32WF555ZUm00MtlTLSLTpT1C1apGkrryhnd/luWjVrlelSGpQG3S1aRCQTYjkxWuUobDIlI/fDERGRpkeBIyIikVDgiIhIJBQ4IiISCQWOiIhEQoEjIiKRUOCIiEgkFDgiIhIJBY6IiERCgSMiIpFQ4IiISCQUOCIiEomMBI6ZfcPMVphZhZklHIHUzE4zs3fMbLWZTY+bfpiZvWRm74bP7aOpXEREDlWmWjjLgXOA1xItYGYx4D7gdGAgMN7MBoazpwMvu/sA4OXwvYiIZLGMBI67v+3u7xxkseHAand/3933AE8A48J544DHwtePAV9LT6UiIpIq2Xw/nB7AR3Hvi4ATw9dd3H09gLuvN7POiTZiZhOACeHb3Wa2PB3FplhHYGOmi0iC6kydhlAjqM5Uayh1HpWKjaQtcMzsr0DXGmbNcPdnk9lEDdPqfHtSd58NzA5rWpSKu9alm+pMrYZQZ0OoEVRnqjWkOlOxnbQFjrt/uZ6bKAJ6xb3vCRSHrz8xs25h66YbsKGe+xIRkTTL5m7RC4EBZtbPzPKAc4E54bw5wIXh6wuBZFpMIiKSQZnqFn22mRUBI4E/mdkL4fTuZvY8gLuXAZOBF4C3gSfdfUW4iVuAU83sXeDU8H0yZqfwx0gn1ZlaDaHOhlAjqM5Ua1J1mnudT4uIiIjUWTYfUhMRkUZEgSMiIpFodIHTUIbNSWY/ZnaUmS2Ne5Sa2eXhvOvNbF3cvDMyUWO43BozWxbWsaiu60dRp5n1MrNXzOzt8Pvxw7h5af0sE33X4uabmd0dzn/LzIYlu27EdZ4X1veWmc03syFx82r8DmSgxtFmVhL3b3ltsutGXOcVcTUuN7NyMzssnBfJZxnu62Ez22AJrk9M+XfT3RvVAziG4CKleUBhgmViwHtAfyAPeBMYGM77OTA9fD0d+Fma6qzTfsKaPwb6hO+vB36U5s8yqRqBNUDH+v6M6awT6AYMC1+3Bf4d92+ets+ytu9a3DJnAH8muPZsBPDPZNeNuM6TgPbh69Mr66ztO5CBGkcDzx3KulHWWW35s4C5UX6Wcfv6IjAMWJ5gfkq/m42uheMNZ9icuu7nFOA9d/8wTfXUpL6fRdZ8lu6+3t2XhK+3EvR87JGmeuLV9l2rNA74pQf+AbSz4PqyZNaNrE53n+/um8O3/yC4Ni5K9fk8suqzrGY88HiaaqmVu78GbKplkZR+Nxtd4CSppmFzKn/57DdsDpBw2Jx6qut+zuXAL+XksJn7cJoOVyVbowMvmtliC4YSquv6UdUJgJn1BT4H/DNucro+y9q+awdbJpl1U6Wu+7qY4C/fSom+A6mUbI0jzexNM/uzmQ2q47qpkPS+zKwVcBrw+7jJUXyWyUrpdzObx1JLyLJk2JyD7qSWOuu4nTzgP4Afx01+ALiJoO6bgNuB72SoxlHuXmzBmHYvmdmq8C+nlEnhZ9mG4D/35e5eGk5OyWeZaJc1TKv+XUu0TCTf04PUcOCCZmMIAufzcZPT/h1IssYlBIedt4Xn4v4ADEhy3VSpy77OAv7u7vGtjCg+y2Sl9LvZIAPHG8iwObXVaWZ12c/pwBJ3/yRu21WvzexB4LlM1ejuxeHzBjN7hqC5/RpZ9lmaWTOCsPmNuz8dt+2UfJYJ1PZdO9gyeUmsmyrJ1ImZHQc8BJzu7p9VTq/lOxBpjXF/RODuz5vZ/WbWMZl1o6wzzgFHLiL6LJOV0u9mUz2klg3D5tRlPwcc4w1/sVY6m+AeQ6l20BrNrLWZta18DXwlrpas+SzNzIBfAG+7+x3V5qXzs6ztu1ZpDnBB2CNoBFASHhpMZt3I6jSz3sDTwPnu/u+46bV9B6KusWv4b42ZDSf4HfdZMutGWWdYXwFwMnHf1wg/y2Sl9rsZRU+IKB8EvzCKgN3AJ8AL4fTuwPNxy51B0FPpPYJDcZXTOxDc1O3d8PmwNNVZ435qqLMVwX+Ygmrr/wpYBrwV/kN3y0SNBL1U3gwfK7L1syQ4/OPh57U0fJwRxWdZ03cNuAS4JHxtBDcbfC+so7C2ddP4f+dgdT4EbI77/BYd7DuQgRonhzW8SdCx4aRs/CzD998Gnqi2XmSfZbi/x4H1wF6C35sXp/O7qaFtREQkEk31kJqIiERMgSMiIpFQ4IiISCQUOCIiEgkFjoiIREKBI1IDC0bvrRzJ96lwCJLK6zyeMLP3zGylmT1vZkeG8/5iZlvMLOGFo2b2qJl9PQX19bVwhF8zKzSzu+u7TZF0U+CI1Gynuw9198HAHuCS8ILCZ4B57n64uw8Erga6hOvcCpwfdaHuvsjdfxD1fkXqSoEjcnCvA0cAY4C97j6rcoa7L3X318PXLwNbk92omZ1iZm9YcO+Th82seTj9WjNbGLauZsddOX+8BYNSLgAmxW1ndGWryoJ7+zxsZvPM7H0z+0Hccj8xs1UW3DPocTP7Ub0+FZE6UuCI1MLMcgnGslsGDAYWp2i7LYBHgW+5+7EE4xpODGff6+4nhK2rlsDYcPojwA/cfeRBNn808FWCMbiuM7NmFtyM8D8JRsk+B0h4c0KRdFHgiNSspZktBRYBawnGYUulo4APfN+YZI8R3AwLYIyZ/dPMlgFfAgaF4261c/dXw2V+Vcu2/+Tuu919I8FApl0IhvZ51t13enA/oD+m+OcROagGOVq0SAR2uvvQ+AlmtgKo9wn/ys3VODFo+dxPMGbVR2Z2PdAiXD7Zcah2x70uJ/h/XuP+RKKkFo5I8uYCzc3se5UTzOwEMzv5ELa1CuhrZkeE788HXiUIF4CNFty75+sA7r4FKDGzynvQnFfH/f0NOMvMWoTbPfMQahapF7VwRJLk7m5mZwMzzWw6sIvg/vOXA5jZ6wTnT9qYWRFwsbu/kGBbu8zsIuCp8DzRQmCWu++24J48y8JtL4xb7SLgYTPbAdS43VpqX2hmcwhGIf6Q4FBhSV22IVJfGi1apIkwszYe3AmzFcENvSa4+5JM1yVNh1o4Ik3HbDMbSHDY7jGFjURNLRwREYmEOg2IiEgkFDgiIhIJBY6IiERCgSMiIpFQ4IiISCT+P2SACCPn5HkKAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.xlim([-1,1])\n",
    "plt.ylim([-1,1])\n",
    "\n",
    "plt.xlabel('PC1 loading')\n",
    "plt.ylabel('PC2 loading')\n",
    "\n",
    "variables = ['OD','BD','LTD','W']\n",
    "for i,txt in enumerate(variables):\n",
    "    plt.plot([0,A[i,0]],[0,A[i,1]],'-')\n",
    "    plt.text(A[i,0],A[i,1],txt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Principal component score\n",
    "This is the projection of each data vector onto new component axes. It is the new variables created from the original data, based on rotating axes. Each principal component score can be thought of as an \"index\" of envrionmental variability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# insert code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "source": [
    "d) Key concept\n",
    "\n",
    "Which principal component is most highly correlated with kelp extent? What is your interpretation of the relationship between kelp extent and environmental variability?"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}