{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Criando uma Transformada de Karhunen-Loeve (TKL)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy.linalg as la"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inserindo dados"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "Data=np.array([[2.5, 2.4],\n",
    "      [0.5, 0.7],\n",
    "      [2.2, 2.9],\n",
    "      [1.9, 2.2],\n",
    "      [3.1, 3.0],\n",
    "      [2.3, 2.7],\n",
    "      [2.0, 1.6],\n",
    "      [1.1, 1.0],\n",
    "      [1.5, 1.6],\n",
    "      [1.1, 0.9]],dtype=float)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Guardar as linhas e coluna"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10, 2)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r,c = np.shape(Data)\n",
    "r,c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Média de cada atributo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.82, 1.9 ])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u=np.mean(Data,axis=0)\n",
    "u"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ajustando os dados "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.68  0.5 ]\n",
      " [-1.32 -1.2 ]\n",
      " [ 0.38  1.  ]\n",
      " [ 0.08  0.3 ]\n",
      " [ 1.28  1.1 ]\n",
      " [ 0.48  0.8 ]\n",
      " [ 0.18 -0.3 ]\n",
      " [-0.72 -0.9 ]\n",
      " [-0.32 -0.3 ]\n",
      " [-0.72 -1.  ]]\n"
     ]
    }
   ],
   "source": [
    "DataAdjust=Data-u #Data - medias\n",
    "print(DataAdjust)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Repare que DataAdjust está em formato de TABELA DE DADOS, o python só calcula matematicamente. Temos que usar a transposta.\n",
    "#### Matriz K de covariância da transposta da tabela DataAdjust"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.59955556, 0.61444444],\n",
       "       [0.61444444, 0.73555556]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "K=np.cov(DataAdjust.T)\n",
    "K"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Autovalores e autovetores de K\n",
    "#### Como K é simétrica, recomenda-se a função eigh()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.2857513 , 0.04935981])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "autoVal,autoVet=la.eigh(K) #,eigvals_only=True) #autovalores e autovetores da matriz de covariância de \n",
    "i=np.argsort(autoVal)[::-1] #índices dos autovalores ordenados em ordem decrescente ([::-1])\n",
    "autoVal=autoVal[i]    #autovalores ordenados por ordem decrescente\n",
    "O=autoVet[:,i]  #autovetores ortogonais associados a cada autovalor em ordem decrescente\n",
    "autoVal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.74498239,  0.66708413],\n",
       "       [ 0.66708413,  0.74498239]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "autoVet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.66708413, -0.74498239],\n",
       "       [ 0.74498239,  0.66708413]])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "O"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### O é o conjunto de autovalores ortogonais de K, que será a covariância da transformação de DataAdjust\n",
    "##### Direção dos autovalores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "eixoMax=np.array([u,u+autoVet[:,0]])\n",
    "eixoMin=np.array([u,u+autoVet[:,1]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fazendo a Transformada.\n",
    "### Dado um vetor de Variáveis Aleatórias X, achamos Y pela transformação Y = O'X\n",
    "#### Se colocarmos em formato de TABELA DE DADOS seria Y'=(O'X)'=X'O''=X'O"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Transformando Data em FinalData, ou dados X em dados Y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "FinalData=Data@O               #Y'=X'O e Y=O'X\n",
    "FinalDataAdjust=DataAdjust@O"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotando"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEVCAYAAAD6u3K7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XecVNX9//HXm6VYAQuKgoiV2LG71o3oV8QSNTF2JdiNiUa/xqhfo4kmxK/la+wSMSZKNP5ssWAFViyLwYIFEcWCYBddBTWLsJ/fH+cMDMPd3dndmbmzO5/n47GPnXLvPZ9bP+eee+8ZmRnOOedcri5pB+Ccc648eYJwzjmXyBOEc865RJ4gnHPOJfIE4ZxzLpEnCOecc4k8QRSBpIclHdPOaewiaXqB4qmVdFwhpuVKR1KNpNntnMYRkh4rVEyusniCyIOk9yR9J2mepE8k/VXSCk0Nb2Z7m9nf2lOmmT1lZoPaM418SdpQ0v+T9LmkryS9IukMSVWlKD9tkm6RdHE7p3GhJJO0XaHiKgQzG2Nm/1Xo6cb5/T7uE/WSnpVUnfX9GpJGS/pI0lxJb0j6naTls4aRpHckvV7o+FxheILI335mtgKwFbAt8D+5A8QNvkMtU0nrAc8Bs4DNzKwXcDCwDbBimrF1FJIEHAV8AbTrzLGD+WfcJ/oATwP3xH1gZaAOWBaoNrMVgT2B3sB6WePvCqwGrCtp29KG7vLRoQ5m5cDMPgAeBjaFRc03f5D0DPAtYWNf1KQjabikpyVdJulLSe9K2jszPUkrxzOSD+P398XPl2heiGcx50h6PQ73V0nLxO9WkvSgpM/idw9K6p/nLP0OeNbMzjCzj+I8Tjezw82sPk5/f0lTY02xVtJGOXGdFc86vom1xtVjM9tcSU9IWikOOzDWsk+I8/uRpDOzptVD0pXxuw/j6x7Zy0PSmZI+jeP+LGfcyyS9H8/ybpC0bEvjSjoBOAL4dawNPxA/P1vSB3Eepksa0swy3AVYEzgNOFRS96y4Wlr/P5M0LZbzjqQTkwqIy/junM+ulnRlVjnvxOm8K+mI7PLja0n6v7gMMmeKmzYzX3kxs++BvwF9gVWAM4C5wJFm9l4cZpaZnWZmr2SNegzwL2AslZVYOw4z878W/oD3gD3i67WAqcBF8X0t8D6wCdAV6BY/Oy5+Pxz4HjgeqAJOBj4EFL9/CPgnsFIcd7f4eQ0wOyeG12L5KwPPABfH71YBfgwsR6j1/z/gvqxxF8WTMG8fAz9rZt43BL4h1AC7Ab8GZgDds+KaBKwO9AM+BV4EtgR6AOOBC+KwAwEDbgeWBzYDPstatr+P01qNUCt9Nms51wAL4jDdgGGEhLxS/P5K4P64bFYEHgBG5jnuLZllGd8PIpxRrZkV93rNLKPRwJ1x2nOAg7K+a2n970OoVQvYLca1Ve42AKwR10Pv+L5rXNZbx2X5NTAoa9hNssp/Or7eC3iBUJMXsBGwRhv3iQuB2+LrHsClwKz4fhLwuxbGXy7GPIyw7X5O3Kb8r3z+Ug+gI/wRDoLzgHpgJnAdsGz8rhb4fc7wtSyZIGZkfbcc4SDZN+7IjZkDVc40Fh0csmI4Kev9MODtJuIdDHyZFE/CsN8DQ5uZ9/OBO7PedwE+AGqy4joi6/u7geuz3v+CmKxYnCB+kPX9/wKj4+u3gWFZ3+0FvJe1PL4DumZ9/ymwQzzYfUPWQRyoBt5tadz4+haWTBDrx+/3ALq1sG1kDnQHxPc3Av/K+r7J9d/E9O4DTmtiG3gYOD6+3hd4Pb5enrBt/pi4XeaUn0kQuwNvxmXWpZ37xIXA/Fjup4SKwNbxu7eyt9Umxj+SUDnoSkgw9cCBxdh//a/tf97ElL8DzKy3ma1tZqeY2XdZ381qYdyPMy/M7Nv4cgXC2cAXZvZlnjFklzOT0KyBpOUk3ShppqSvgYlAb+V3kXkOIVE1Zc1YVib+xhhHv6xhPsl6/V3C+9wL+onzkVtWzncAc8xsQdb7b+O0+xAOvC/EZrB64JH4eUvjLsXMZgCnEw6Cn0q6Q9KaScMCBxLOTsbG92OAvSVll93U+kfS3pImSfoixj0MWLWJsv5GOLAS/98ap/kNcAhwEvCRpIck/SBhvsYD1wDXAp9IGiWpZ+5wCnfQzYt/U5uIBULFobeZrWZmu5vZC/HzlrYpCE1Kd5rZAjNrAO7Bm5nKjieIwmhrl7izgJUl9c5z+LWyXg8gNFUAnEloFtnezHoSLv5BqFm35AlCzbMpHwJrZ95IUozjgzxjTtLUfCxRVs53zfmckIg2iQes3mbWy8IF1Hwstf7M7B9mtnOMx4BLmhj3GMLB/n1JHxOa97oBh7VUaLy+cjdwGbC6mfUmJJqm1tt9wObxusG+hGSUifdRM9uTcGB+A/hL0gTM7Coz25rQJLohcFbCME+Z2Qrxb5OW5iPBE8CBauKGjXh9bHfgSEkfx+X2E2CYpKaSo0uBJ4gUWbgo/DBwXbzQ3E3Srs2M8nNJ/RXuEjmXcO0CQpv7d0B9/O6CVoRxAbCjpEsl9QWQtL6k22LiuhPYR9IQSd0IyaiBcH2grc6PZz2bAD/Lmo/bgf+R1CceKH4L3NbSxOJZzV+A/5O0WpyHfpL2yjOeT4B1M28kDZK0ezyA/4ewbBfmjiSpHzCEcLAeHP+2ICSTfGrD3QnNK58BC+LF6yZvSTWz/wB3Af8A/m1m78c4Vle4kWB5wrqZ10S820raPq7Hb+K8LTVcAVwB9AT+JmntWHY/SVdI2pxwx9ebhEpNZrltCMwmj8TqSscTRPqOIlwHeIPQlnt6M8P+A3gMeCf+Ze7dv5JwS+HnhAuEj+RbuJm9TWivHwhMlfQVoVb7PDDXzKYTmjOujtPfj3DL7/x8y0jwJOFC9zjgMjPLPMh1cSz3FeBVwsXufJ9PODtOc1JsZnuCcADKx2hg49g8dR/hoP0nwvx+TLhofm7CeEcBU8zsMTP7OPMHXMXimn6TzGwu8EtCEv4SOJxwob05fyNc3L8167MuhMT9IeFW292AUxLG7UlIpF8Smu/mEM5eCsrMvgB2JGzXz0maS1jXXxHW0THAddnLLC63G/BmprKSuZPClTlJ7xEuND+RdixtJWkg8C7hwu+C5od2SSQNIFQm+prZ12nH4zo3P4NwroOIbfpnAHd4cnCl0DXtAJxzLYvXFz4hNA0NTTkcVyG8ick551wib2JyzjmXyBOEc865RJ4gmhA7eJspaZxCB3V53Z8dx7uwDeUNj53C1UqaIGn1VpS3bstDurSVeptqYlpDJe3TiuFHtGLY2jbGdLCkA9oybs50Fi0nSVe3d3qtKPdPyr9zzA7FE0TzbjWzIcDewBGStipyeZeaWQ3hXvXDs79o6qlUQn89niA6jlJvU0sws0fM7KFWjJJ3gmiHwwmdKxaMmf2ikNODZvfB2wjdnHQ6niDyEPtduhzYLz7tPE7SREl3Z/o7knSzpCcID08RPztb0jOSxksaoNC1d+YM4apmiuxJ6ACO2E/P9cBlCj9ss378vFahW+nhwOWSLpfUS6Gr74mZ6Us6UNK/YwzDirB4XBuUYpuSNFjSk3EbOjd+NlzScbG2fXHWZ8MVnqB/Nk7rXIWu0DeL099M0p/j9J6Kz2Og0HX7JEl/yip3j/jZJEl7xM/+FsedkH2gVXjyv4uZLVToDn68pDslvSTpIEmPKXSXvnwc/rcxnvEKz9U0tZwyXZyfE8t9TtKWOcvnhDitiZJmxM+Oi/P3VCZ5S3pZ0m2ELuG3iMt/kqQj47p8jdABYueTdm+B5fpHqJnndgF9PaGfnExPrhcTusHeDvhL/OxcQidvfYFH42c7x3H3AC6MnymnvOHAdMIPr8wk9vZJePK0f3x9C7B+fF0b/1/I4u6yzwKOiq9vArYn1G4GJpXpf51+m1o28xkwIb4fDhyXHUv8LPP58OxpEXuCja+Xi//3AP5AuE3+ufh/x6xt8mlCJacnoUuWbsC4JmLcDrg6vh4IvESouB5O7BU3zv9BhCfIb4yfbUToOXep5ZQdd1bM6wNjmlgvIwm/CbIq4Ul2EbqNz/RCPAdYPr6+P8bZLc57t/j5xLS3r2L8+RlE/voBHxG6Vh4t6UlCB2NrEpp4XorDZXq0HEjoMgJC9xHrE7qY6CLpHyzulTPbpRY6iNuXsAMCfGpmmR8Oyr4nOalDt/UI3VNkl3kxoX+jW+J7Vz6KvU2tA4yN092I0GVIRtK2dCehi5AxJD9r8WtJTxG2qTUJB9SZFp6KfyFrODOzry08zLfQ4g8KxVr4xTlNNbnb8esW+tb6kPD7J8TXK8V5qFG41nE9IQElLadsR0maSKgwLdUjr6SDCMl5TJzWFoRkeg/hdzMAplvoMRdC1/zvxXl6lyWXaafjCSIPCr/cdjqhnXQv4E0z243QZ5EIG8oWcfDMaex7WZ9tQ/itgyoz+62ZHU7oO6cp9YQaDITfi8j4ClhDoRO5DeJn3xN+iAZC/0xb55Q508yOA0YRnsJ1ZaBE29TJwCVxujNY8mD8FYu75N4s/v/ezM4gdKD4+/iZxXhXIfwGyC6E3wgRoa+qtWOTWHbzTRdJPRW6Eq+K399uZkcSumDP/nnRt1iyB19r4rUIZ9iPmVmNhWt1RzexnLKdQjhbOj5n/lHoEv1Y4L/jR+8Ck7Omv2f8PHsfrI9NYd0ICeXT+Hl7+iYrW/4kdfOOUvgh9ipglJm9JOkz4DxJ2xB2srfM7DlJJ0saR2geet/MPo7trc8SNp5jgO0k/ZFweprUp9JZCj8V2Z3kg/nfgb8CU1j8GwO1wB8lbU/oJO4fko4HXjGzSZIukbQDoUvq5pKSK41SblMPAddIep2lD2CvAGtKGktoQgHYX9KphN/WyPSiO0vhp07PB+ZJGh/HxcwWSPoroRnpyaxp/57QqaQIPfKuCNwfE8XXhI4YidP4QlKj8vjtEjN7WaF78FriLxOa2ajc5ZQz2r8Jv48yMWGSZxG6nn9C0sdmdqjCb2lMJPRyOx64KGec3xI6zawCrjWz7yVtFsvpdPxJaucqiMJvXs+LTSplQdLBwAIzuzftWNoiXqC/1sxa+uGwDscThHMVIt7FcxVwiJnl80NMrsLlnSDiKeDzwAdmtm9Ro3KuBCStRWi260toZx5lZn9ONyrnykdrrkGcBkwj3DngXGewADjTzF6UtCLhN60fN7PX0w7MuXKQ111MCo+R70O4Vcy5TsHMPjKzF+PruYQKUL90o3KufOR7BnEl8GvC3QiJ4lOXJwAss8wyWw8YMKD90bVSY2MjXbqU/s7dtMpNs+w05/nNN9/83Mz6FHKa8ancLQkPP2V/nvp23ZI010VzPK7WKcZ23W4tPUlHeGjrOlv8JOiDLY2z4YYbWhomTJhQUeWmWXaa8ww8bwV8WpRwC/ALwEHNDZfWdt2SNNdFczyu1in0dl2Iv3zS6E6E+6PfA+4Ado9PRDrX4cUHnu4mdMNwT9rxOFdOWkwQZnaOmfU3s4HAocB4C09EOtehSRIwGphmZlekHY9z5ab8GuKcK52dCD2A7i5pSvzzHm+di1rV1YaZ1RK6dnCuwzOzp0nu9NA5h59BOOeca4InCOecc4k8QTjnnEvkCcI551wiTxDOOecSeYJwzjmXyBOEc865RJ4gnHPOJfIE4ZxzLpEnCOecc4k8QTjnnEvkCcI551wiTxCuXerq6hg5ciR1dXVph+KcK7BW9ebqXLa6ujqGDBnC/Pnz6d69O+PGjaO6ujrtsJxzBeJnEJ1QqWr1tbW1zJ8/n4ULFzJ//nxqa2uLWp5zrrT8DKKTKWWtvqamhu7duy8qq6ampijlOJeGuro6amtrqampqdgzY08QnUxSrb5YG3d1dTXjxo3rsDuRpJuBfYFPzWzTtONx5cObTwNPEJ1MqWv11dXVHXnHuQW4Bvh7ynG4LOVQcy9lRauceYLoZHJr9QAjR47skDX8YjOziZIGph2HW6xcau7efBp4guiEMrX6ctnZOjJJJwAnAPTp06csL8TPmzev08Q1ZswYGhoaaGxspKGhgZtvvpmGhoZU4rr00kuZMmUKgwcPpqGhoSyXcbF5gujEmjxNbpgLD/8Gtjse1hycdphlzcxGAaMABg0aZOVYk8w+WywnbYmrR48ejBkzZlGlZsSIEQWv1OQbVzku01LzBFFAM+d8w4CVl0NS2qEATZwmf/AC3HUs1M+Eflt6gnBlpaPf+NDZeIIokM/nNbDvVU+zWf9eXLDfJgzqu2LaIS25s+22K9WN/4bRv4cV+sLwh2DtHdMO0bml5N74UA4XrSuVJ4gC6b1sN84aOojLH3uTYVc9xZHbD+BXe25I7+W6pxpXdXU11ZuuA/eeBO9MgI32h/2vgmVXSjWuciDpdqAGWFXSbOACMxudblQum19HS5c/SV0gXau6cHT1QGr/u4bDtxvArZNm8sPLarl10kwWNlp6gb35GFy/E7w/Cfa9En76d1h2Je9DCTCzw8xsDTPrZmb9PTmUn2af1jeDyaPh2y9Si6+z8zOIAltp+e5cdMCmHL79AH73wFTOv+81xkyayQX7bUL1equULpAFDfDEhTDpOlhtE/jJzbDaDwCvlbmOo8nbTed9BvedDDMeDzdd7Hx6qnF2Vn4GUSQbrdGT24/fgeuO2Iq5/1nAYX+ZxCljXmD2l98Wv/DP34KbhoTksN0JcPz4RckBvA+lSnT3C7Ope3tO2mG0WuY62kUXXbS4IvP2eLh+R3h3Igy7DHY6Le0wOy0/gygiSQzbbA12/8Fq3PjkO1z/5AzGTfuUE3dbj5N3W49lu1cVtkAzeOk2ePjX0HUZOPR2+MGwpQbzh4Aqy8JGY9TEd5j+yVyGbdaXc4dtRP+Vlks7rLwtumi9YD48dj48exX0+QEcfR+svkna4XVqniBKYJluVZy2xwb8ZJv+jBw7javGvcVdz8/inGEbse/ma7Tqttgm7+j4rh4e/BVMvQcG7gIHjYKeayZOw28lrCxVXcS/Tt2JURPf4braIldSimXO23D3sfDhS7DNCPivP0D3jpPkOipPECXUr/eyXHP4Vhxd/QUX3j+VX9z+ErdOmskF+23MJmv2anH8pGsHALz/HNx9HHz9Aex+Puz8K+jS/I7fwftQcq20TLcqfjlkA36ydX9GPvxGuyopJffyHfDQmdClK/z0Vth4/0Vf+S2wxeXXIFKw3Tor88AvduaPB27GjE/nsd/VT3Puva/yxTfzmx0v99rBkxPGM2DmnfDXvUHAiEdh1/9uMTm4yrVm72W5+rAtufPEanov151f3P4Sh9w4idc++Crt0Jb2n6/hnhPg3hOh7+Zw8jNLJYchQ4Zw/vnnM2TIkIq+I69YWkwQktaSNEHSNElTJfkVoQKo6iIO334AE86s4ZgdB/LPybOouXQCf33mXb5f2Jg4TubaQVVVFQNX7s7Jyz/Cuu+OgU0OgJOehrW2LfFcuI5qiUrKZ/PY75qnOeeeV5kzr7D9HrXZBy/AjbvCq/8Pas6F4Q9Cr/5LDOI3WxRfPmcQC4AzzWwjYAfg55I2Lm5YlaPXct24YL9NeOS0Xdhird787oHXGfbnp3j6rc+XGjZz7eD2C49m2mkr02ve27wx6Jfw49GwTMtNVM5ly66kDN9xIHc+P4sfXlbLzU83XUkpusZGePpKGP1f0LgAho+FmrMTz4qzK0x+s0VxtJggzOwjM3sxvp4LTAP6FTuwSrPB6ivy9xHbMeqorWlY0MiRo5/j+L8/z8w53ywe6PvvqP7ibg5ecDfdVlkHTpzIx2sMgXJuP3ZlL7eS8vsHQyXlqbc+K20gcz+G2w6EJy6AH+wDJz0Fazd9XSHxFlhXUK26SB37zt8SeC7hu9S7RU6r2+NCltsdOH8bePS9bjww/ROGTPuEoet045C+H7Pl9MtY4ZuZzOr/I95Z9yjstdmdYp5dechUUp6Y9ikXP/Q6R43+N3tuvDr/s89GrL3K8m2ebl4Xkt98DO47CeZ/C/tdBVsdnVfFx2+2KK68E4SkFYC7gdPN7Ovc78uhW+S0uj0uRrl7Amd9/R8uefgN5r18H9t+eC1dl+0JR9zNWhvswVpFLDsf5drFtGsfSey58ersuuGqjH76Xa4ZP4M9r5jIsbusw6k/XJ/le7TuxscWn9rPfuJ/9U3DE/99BhV2plyb5XUXk6RuhOQwxszuKW5ILmP1nstwxSGD+cXhBzF1+R1YcMJTsMEeaYflKkCPrlWcUrM+E/67hn03X4Pra9/mh5fVcs+Ls2lsRd9izV5Izn7if/uT4LhxnhzKTD53MQkYDUwzsyuKH5LLtdmmm7P1WQ+wzErJD745VyyZSso9p+zIGr2W4Yw7X+bHNzzLy7Pq8xo/8UKyGbx4a7hL6asP4LA7YO9LoNsyxZ0Z12r5nC/uBBwFvCppSvzsXDMbW7ywnHPlZKsBK3HvKTtx14uz+d9HpvOja5/h4K37c9bQQay2YtMH9qWe2h+8Edw1Ijzxv86ucOAo6LlGCefEtUaLCcLMniY8huWcq2BduoifbrMWe2/al2vGz+DmZ97l4dc+5pdD1mf4juvQvWtyg8SiC8nvPwc37BKe+B9yQehkzx/qLGv+JHXK2vO7DHV1dYwZM8afIG0HSUMlTZc0Q9Jv0o6nI1hxmW6cM2wjHj19V7ZbZ2X+OPYN9rpyIuPf+CR5hMaFMPHS+MS/4NjHYJczPDl0AN4XU4ra87sMmXEbGhoYM2aM3wfeBpKqgGsJN43NBiZLut/MXk83so5h3T4rcPPwbZkw/VMueuB1RtzyPDWD+nD+vhuzXp8VwkBffRC6ynjvKdjsYNjncn+oswPxM4gUtaergMy4jY2N3s1A220HzDCzd8xsPnAH8KOUY+pwfjhoNR45fVfOG7YRL7z3JXv930T+8NDrrPDxJLhhJ/jgRTjgBjjoL54cOhg/g0hRe36XITNuQ0ODdzPQdv2AWVnvZwPbpxRLh9a9axeO33VdDtiyH5c9Oh3VXc02Xf/B/NU2p/sht8Aq66UdomsDTxApas/vMmTGvfnmmxkxYoQ3L7VN0s0XS9zkXw49BLSk3J5q33tV+GKzHRj77lyWG3Q0vDqLJfNwuspteZUzTxApa09XAdXV1TQ0NHhyaLvZsOihdID+wIfZA5RDDwEtKden2mtr1y7TuMpzeZUjvwbhKtlkYANJ60jqDhwK3J9yTM6VDT+DcBXLzBZIOhV4FKgCbjazqSmH5VzZ8AThKlrsEcB7BXAugTcxOeda1J4HOos5LVdcfgbhnGtWex7oLOa0XPH5GYRzrlmF/O1n/x3pjsUTRAH5qbPrjAr528/+O9IdizcxFYifOrvOqj0PdBZzWq74PEEUSG1tLQ0NDTQ2NtLQ0EBtba1v/K7TKORvP/vvSHcc3sRUIKussgqNjY0ANDY2ssoqq6QckXPOtY8niAKZM2cOXbqExdmlSxfmzJmTckTOOdc+niAKpKamhh49elBVVUWPHj384ptzrsPzaxCEC8ztvWjmF9+cc51NxSeIQt595BffnHOdScU3MfmDO845l6ziE4Q/uOOcc8kqvonJrx0451yyik8Q4NcOnHMuScU3MTlXKbyvMNdafgbhKpKkg4ELgY2A7czs+XQjKi7vK8y1hZ9BuEr1GnAQMDHtQEoh01fYwoULF/UV5lxL/AzCVSQzmwYgKe1QSsL7CnNt4QnCuWZIOgE4AaBPnz5lWfOeN29ei3FNnjwZSZgZXbp0YfLkyWy44Yapx5WGco2rHHmCcJ2WpCeAvglfnWdm/8pnGmY2ChgFMGjQICvH52Qyt2g3p0ePHowZM2bRNYgRI0YU/RpEPnGloVzjKkd5JQhJQ4E/A1XATWb2p6JG5VwBmNkeacdQLvx5H9cWLSYISVXAtcCewGxgsqT7zez1YgfnXKUpRMeRTfHnfVxr5XMGsR0ww8zeAZB0B/AjwBOE67AkHQhcDfQBHpI0xcz2SjMmvxXVlZt8EkQ/YFbW+9nA9rkDlcPFvLQuPqV50asS57kQzOxe4N6048iW1HGkJwiXpnwSRNJ9gLbUB2VwMS+ti09pXvSqxHnurDIdR2bOIHz5urTlkyBmA2tlve8PfFiccJyrXH4h2ZWbfBLEZGADSesAHwCHAocXNSrnKpRfSHblpMUEYWYLJJ0KPEq4zfVmM5ta9Micc86lKq/nIMxsLDC2yLE455wrI95Zn3POuUSeIJxzziXyBOGccy6RJwjnnHOJPEE455xL5AnCOedcIk8QzjnnEnmCcM45l8gThHPOuUSeIJxzziXyBOGccy6RJwhXkSRdKukNSa9IuldS77Rjcq7ceIJwlepxYFMz2xx4Ezgn5XicKzueIFxFMrPHzGxBfDuJ8ENYzrkseXX33VpvvvnmPEnTizHtFqwKfF5B5aZZdprzPKjA0xsB/DPpi+zfWgcaJL1W4LILIc110RyPq3UKvV23m8yW+nnp9k9Uet7Mtin4hL3csim7I8yzpCeAvglfnWdm/4rDnAdsAxxkLewMac5zczyu1vG48leUMwjnyoGZ7dHc95KOAfYFhrSUHJyrRJ4gXEWSNBQ4G9jNzL5NOx7nylGxLlKPKtJ0vdzyKbujz/M1wIrA45KmSLqhROUWg8fVOh5XnopyDcI551zH57e5OuecS+QJwjnnXKJ2JwhJNZJmShonqVbSYZIOljRVUqOkpW7biuNc2IoypkuaIGmSpP2bGGaopPckzZL0m1ZM+0pJVfkOnzPuzZI+zb43XtItkga2ZXo5075Q0stxuT4saeus79aKy2NaXM6n5YxX08Yyd5FUL6l7M8OsLumtGNtUSb9rY1kDJe3ehvGqYvmNklaLn20ryfJZ7pnlI2mwpGPbUH7Ru+iQNKKVww+VNLK5fa6V02vV/tlMTNMlzcjdH+P0L27P9Jspd3j2vpLz3VL7azPTGSppnzbG0E/SJa0Yvsn9uRXTqJXUNW7fzd69F4c/TdJ2LQ1XqDOIW81sCLA3cASwEDgImFig6X9mZj8EfgicnvtlPMBfG/8uBg6TtHE+Ezaz081sYRvjugUY2sZx83FmXK7HA9dKWi5+viB+txGwA/DzfOe3BQcBdwFDmhlmWeAFM9vTu/nDAAAZBUlEQVQCGAwMlbRDG8oaCOSVICRlb6enATOBr4Efxc8OBJ5vTeFmNsXMRrdmnKjZLjpyYm1WM8O2KkGY2SPA3ynsPtdmWfvj3sDGtGJ/bC8zu8XMXmji61uI+2tL68nMHjGzh9oYxsnAba0YPnt/rqZw+3Nz/g6c0tJABW1iMrPvgMuBzYDrCAeQSzI19JjBnwCOyowj6WxJz0gaL2mApJVjNpwg6aqcIlYAusXxekl6UNJE4HZgBuFAcQZgwJFxGk9Lui6O00XSTZKelPRw/CyTeQfEGJ6RdHb87sI4/BOSboqfDY3jP084yH0BdFM4u7mP2GVDdnyZ+ZB0oKR/x3KGtWK5zgYeBbaVNBi4A7hO0rlmNjfO+z8kjY3LnDhPt8fyb4/vd5T0XCw/6SC0IfA74IA4jUU1vVgzG054qviHkmqBNeL83pi1/o6QdHIcZ3NJ1yq4Pg7zkKSV4nSOkjQuDntVjPXBuOwGxvV3FzBc0vaS6oDfAp8SnoQ9V9KThAQ6NU5n2Ti/4yX9U1K3uE1NyFk+2fN2Z1ynj0nqmb1AJF0Xt5FXJI0mJIirJY0nJNJ147Tul/QAsJekI+P28IykLXKmNzzG9RCwecL2fwKwWSxzM0l/jrE9JWlAnMYBcfoTJO0W18tOZjYdWB8YlbMcx0u6S9ILkvrHZTIuLu+7VaD9M87bncAzhDvEZgPfA98BD2Ste5qafs53F0r6i8L+d4Ok8+N8/zZ+v2hflHR01jh7KGH7JxzvLgEGAHtllZO0zQyXdJykLSX9Iw53m6RtlLBv56gxs1cl7aR4JhGX233x9W/jMhyvcNbbA7g8bus/AVYG/ifO1/Gx3Jcl7RXH3z6O/4yknyWUf6iWPGYtdWwzsy+BNSUpYfxFinEN4kPC06v7AlOA94DdFU5nFsaHl96OgfcFdjeznQg7/jnAVkBtPGPInGr1UUgEbxPOECAcYP5pZrsSFuh/CLeJXUpIUr2BPc1sZ6CnpA0INc5PzWw3IPf08WzgghjL7pLWjJ9PjTEPUGhOmBjH34HFXTCsCpxJWLlrJMS3nKTtgR8DPzWz3YGH44qrzfkb08JynU7YAHcA9pQ0iFDruMbMhhFq+BCS5eux/Kmx7GHA2bH8v2ZPXNJWwPNmNgtYXU3XsEYBT8TlOw1oiGcTmfX3AGHdw+Izkn2B92O51wAnxencamZDJG0LLB9jvSN+D7AacIiZ3Qz8HviEcFDONCH0IKy37DOB44D7Y1m1hHVyHHBTzvLJNjyu0zuBQ7K/MLNTgD2AWYRtK3teFgCNcdDuZrYf8BjwS2AXwtn0HxLKqzezfYCPydn+zWwU8KqZ1ZjZq/Gz3QiJ+8S4Xs4Dfhj3kacyE43LsYqw7WUvx5WAnwJXELaDBcC+cXlPo/37Z7avgMsI+/1BmeUFPMzidZ+JN2n6uV6PMa0LvBa3+0wzc9K+mJG0/QN0J6y/h7OGTdpmADCzl4B3Jd0IfGhmz5O8b2fLNNE+G2MjxvwvSZsB/cysBvh51jyvRtj2xhMqwtcBOwF/Ihxb9o7DQ9gX9gd2Bo7Q0k3Cucespo5tX7D4eJWozQkiZqjXCAeaEyS9Ft8fDHxE2GkHE3boNQkr+KU4euYUcCDwSnz9PKH28yTQJWbtI+N3n8WV8RMWN0usB7wYX78DLFHzIxw87lKo6e4cY9iQsNIws8ac4bOn9xKwTnydaa/8EOgFbB1rWeOATeJ33YGXYudvmfnJnl5m3i4m1AxuAdY3s/fjgSD77wiS9SMs13WAsQo1540JB+B/Z+Yrq8yk8q8DfirpVmDbnOkfBAyR9AihT5gdCWdiGdk1DTOzwYSE21XSppkyzOxrYL6kVQkHyYnARoRaTS3h4LZyTtlJsQK8nNX8t12M+XLCuq4C7gHGEnaujI2A02NZx8Tvsre9F7OGzTSHXBorIKcSalVPZLbnuE1/TDgrHpQ1L+8AGxAOftnT7QPMNLPvzew9wjaTq7ntP9evJT1F2HZ+RjjYbQhMjrG9AmwZh10PmJswvdfj9v4BIbEvD4yO29BPaP/+me0lwrYyJw6faQY9gKXXfT7zn73/ZV7Pi+staV/MaGqbSrr+kLTNZLsBOBb4cwvTzjAIOwnwiqQtiQkillUTy7qexcetlwmVl7sJLQLPmVkD8IaZfWJmHxISPcAWwP3ABEKlsU9O+bnHrKaObWLJfXwpbU4QZraHmW1K2GhHxdfbEHbkbwjts1MIK0/Au3HGYPEG/V7WZ9sQai5VZvZbMzuckDmzy3wM2Caepr7D4ppkX8LC/Z5w4OgfP7svZupnYgzTiRk9oYacPb0tY2yw9EHy14Qaxx6E2hLAfGCLuNFuljC9zLzNNLPjCLXnM/I9g5DUD9gTmExo37wklt8duI9Q68hdtknlfxlrxGcTaqTZtjWznc1sKKH2dWCcv0wNIzNfmWUMYUOcT2jXzZRBjOnXwFvxAD8d+HtMgDsD5+ZMJylWWFw7h1DbgXBQ6UFYx32Aewk14lXj99OB/41l7UBIiknbXsZgFp+9XEt4NmgPM9s0btMjCdv3AAv9N00HXieczawV5yU71s+AgbGZYiCLt5FsmWHfY+ntH+I2J2kVwtniLsD5hLOTTWIM28T4Nmfxgf0dQtNO4vQiEZpX3oy177sp0P4ZbUFoWlo3Dj8deBW4NmvdZzQ1/9msiddN7YsZ+WxTGUnbTLY/Ec6Wft/CtDO+z3p9F+EYWWVmX8SyHstUCIGjs+btbmAMofnUsj7PnmcI63ufOP6WZvZBTvm54zR1bFuZsB03qVBdbRwlqZqww48i1BrvB9YmtD9iZs9JOlmhzXkm4TTvY4W2zGcJB5pjgO0k/ZFwreGJhLLGEJMSod39eMIGuClhw/wDoaZ3EnChpAOyxr0f2C/WFucRmlwyLgH+Fk/XHjCzD5ponruXUBOYAnwZP5sDXElY2JkFnh3fK2Y2SdIlChd0VyBclHofqGlyqYZ2yc/jsjnVzL5TaLu+Jk5jHuFMbS5wt0I7bENWnLfFef0ozt/PJR0Ux820jV4dp/dpplAzeyOuz7MINeqxcR4hHIxXi+2lvySs82MJTQnHxGH+RagdZS4i3w9cpdBuT1xWTwIjJf3TzA6RdEysKc8FDifUdLMdRjhYdyHUpFYkJK8qwgXrTO+co4C/SDqFsHOcA9yUsHwypgPrxzOnWYRadrbfAR/GGt8jhPW+F6Em/2Ccl6+zlt1CSdcQmn4aWdwssJQmtn+AWZLuJiSFeXG5vRLHaZQ0EnhS0jdkJXoz+7ekRuAvhHWetBwBngPOU7jb6StCIi/U/rkK4UDan1BrfZew/FeVtGf28mpm/vOVtC9mf5e7/e/UxHSSthkgXO8hVOyuVbgetCcJ+3bO9CZK2tTMXiNsB2OI68nMXpb0cdyejHD99DFCM/H9ZnaFmrhTM8sFwP2xkvsFi5vPmpJ0bFuJ0GSWlDAXM7OC/xFqn7MJO+MnwKPFKCenzGGEs5a3Cb11FrW8WObthI3v+zi/x5ao3J3jxvUKYeeYAgwrUdmZGusrhDOI35ai3JwYaoAHS11uLHsGIZFklvsNacSREFfJ97mEGIYDx8XXJd8f84ivJPsrITle0orhS74/Eyp327c0nHe14ZwrCIW7qbqa2U1px+IKI+8EEdvXnwc+MLN9WxreuXInaS3C/eB9Cc1Bo8zsz82P5VzlaM01iNMIt8Tl3i3kXEeVeUDpRUkrAi9IetzMXk87MOfKQV53MUnqT3huwE8dXadhZh+Z2Yvx9VxCBahfulE5Vz7yPYO4knBL2YpNDaCs3+5dZpllth4wYEBTgxbMx9+EC/B9lw95rrGxkS5dSt//YFrlpll2mvP85ptvfm5mufd+t0u8JXVLwh0+2Z+XfLturTTXRVM+/qYRA9ZYvrzigvJcXlCc7brd8rjavS9wXXxdQx53j2y44YZWCj+94Vn76Q3PLno/YcKEkpSbK61y0yw7zXkmPPFdyDs6ViA8HHZQc8OVarturTTXRVN+esOz9l9/Gpt2GInKcXmZFX67LsRfPml0J2B/Se8RHt/fXVJrOqJyrmxJ6kZ8QMnM7kk7HufKSYsJwszOMbP+ZjYQOBQYb2ZJj9g716EoPAk5GphmZlekHY9z5ab8GuKcK52dCD2X7q7wu9RT1Ipedp3r7FrV1YaZ1RJ6O3SuwzOzp1myE0LnXBY/g3DOOZfIE4RzzrlEniCcc84l8gThnHMukScI55xziTxBOOecS+QJwjnnXCJPEM455xJ5gnDOOZfIE4RzzrlEniCcc84l8gThnHMukScI55xziTxBOOecS+QJwjnnXCJPEK5iSbpZ0qeSXks7FufKkScIV8luAYamHYRz5coThKtYZjYR+CLtOJwrV54gnHPOJWrVb1I7V2kknQCcANCnTx9qa2vTDSjBvHnzyi6u+vrvWLhwYdnFBeW5vMqVJwjnmmFmo4BRAIMGDbKampp0A0pQW1tLucV1/fQ66uvryy4uKM/lVa68ick551wiTxCuYkm6HagDBkmaLenYtGNyrpx4E5OrWGZ2WNoxOFfO/AzCOedcIk8QzjnnEnmCcM45l8gThHPOuUSeIJxzziXyBOGccy5RiwlC0lqSJkiaJmmqpNNKEZhzzrl05fMcxALgTDN7UdKKwAuSHjez14scm3POuRS1eAZhZh+Z2Yvx9VxgGtCv2IE551yh1dXVMWbMGOrq6tIOpUNo1ZPUkgYCWwLPJXxX8l4v6+u/A1hUVlq9NKbZO2QlzrNzbVFXV8eQIUNoaGhgzJgxjBs3jurq6rTDKmt5JwhJKwB3A6eb2de536fR6+X100MtoKYmrOS0emlMs3fISpxn59qitraW+fPn09jYyPz586mtrfUE0YK87mKS1I2QHMaY2T3FDck55wqvpqaG7t2706VLF7p37+4VnDzkcxeTgNHANDO7ovghOedc4VVXVzNu3DhGjBjhzUt5yqeJaSfgKOBVSVPiZ+ea2djiheVcx1NXV7eo6c0PPuUjd700NDT4+slTiwnCzJ4GVIJYnOuwMhdA58+fT/fu3b2G2grFTKxJ68Xlz5+kdq4AMhdAFy5cuOgCqGtZ5gB+/vnnM2TIkILffurrpX08QbiKJmmopOmSZkj6TVunk7kAWlVV1SkugNbV1TFy5MiiPy9Q7AN4Z1svpea/KOcqlqQq4FpgT2A2MFnS/W3pJSBzAbQzXIMoZXNZ5gCeKavQB/Ck9eJnEfnzBOEq2XbADDN7B0DSHcCPgMQEsdy3H8Bf92lyYtVAdV/gjafhjcIH25TB9fXwbu+CTW/N92fy8CFdMOuBJNZ87Hh4Y+1WTeO3c75iwYIF8Ndrmh2uGvho5FbU19fTu3dver1xccGXXe56KfTy6sw8QbhK1g+YlfV+NrB99gDZPQRs2rcH9fX1pYsuxzfffsO8efNYYYUVWH655Rd9vnDhwoLG1bVrVyRhBpLo2rVrq6e/YMFCMPIer1fPXlijlWT5Fnp5dWaeIFwlS7o7z5Z4k9NDQO9fPVOKuJbSXLNPoZ9q7w28n3VnUb82NC+deGMd9fX1PPqrvQsWV6GUbS8AZ5TfzaKeIFwlmw2slfW+P/BhSrEkytwC+v777y91MbeY1zmqq6s79HUUVxieIFwlmwxsIGkd4APgUODwdENaLPusoWvXrlRVVQH43Th58IcWC8MThKtYZrZA0qnAo0AVcLOZTU05rEWybwEFOP744xkwYIAf9FpQDg8tdpYE5QnCVbTYZUxZdhuTewvo0Ucf3aEPNqWS9GxFKZdbOSSoQvEE4VyZ6kzPVpRSsZ+taEnaCaqQPEE4V8Y6ysXicmpSSTuxpp2gCskThHOuXcqxQ7w0E2vaCaqQPEE459olsT+llWvSDitVHeXMryXeWZ9zrl28Q7zOyxOEc65dMk0qF110UYe+Y8ctzZuYnHPt1lmaVNyS/AzCOedcIk8QzjnnEnmCcM45l8gThHPOuUSeIJxzziXyBOGccy6RJwjnnHOJPEE455xL5AnCOedcIk8QLlV1dXWMHDmSurq6kpYr6WBJUyU1StqmpIU710F4VxsuNSn/8tZrwEHAjaUq0LmOxs8gXGoSu4kuETObZmbTS1agcx2Qn0G41HSEX96SdAJwAkCfPn1KmsTyNW/evLKLq77+OxYuXFh2cUF5Lq9y5QmilcrppxU7ssxyvPLKK5kzZ05RlqekJ4C+CV+dZ2b/ymcaZjYKGAUwaNAgK8ckltkey8n10+uor68vu7igPJdXucorQUgaCvwZqAJuMrM/FTWqMpVym3mqCpkYS7UczWyPgk/UuQrSYoKQVAVcC+wJzAYmS7rfzF4vdnDlJqnNvBISRKF/c7hSl6NzHU0+F6m3A2aY2TtmNh+4A/hRccMqTx3lpxULfetooS8ml8NylHSgpNlANfCQpEdLHoRzZS6fJqZ+wKys97OB7XMHSuNiXn39dwCLyirFxadLL72UKVOmMHjwYBoaGqitrU31oldu2VOnTuXMM8/k+++/p1u3blx++eVssskm7SqjZ8+edO3aFTOja9eu9OzZs93znLQcS8nM7gXuLWmhznUw+SQIJXxmS32QdTFvy7VWsJp3L21naC3rs/ArADZ5txcA9fX19O7du6hl1iwH7AjwIbw7tl3lfvX1V4vG7dWzV5viyS17nc9nsulhXTGrQhIDP7+etd9du03TzqhZDo68ZJvFsS43lvq32resk5ajc6685JMgZgNrZb3vD3xYnHAqw1dff8Unn3zCRx99DBiS2GKLLdqcJLL17t0bKeR0SQVLmL169ipIfM65jiOfBDEZ2EDSOsAHwKHA4c2N8O1y/eBnDxUgvOb9/sbQxv7Pn4ULnFNSun2tNeVmLvj+5z//wSyciFVVVXHRRftwzmnntLvsXkC3rDuOehXp4m9blnXB7oQakXRS65wrtBYThJktkHQq8CjhNtebzWxq0SPrpDIXfDPJQVLBL9RWV1eX3V1BlXyLsHMdVV5dbZjZWDPb0MzWM7M/FDuoziz7Dp4ePXpw4oknVsTBMs1uNZxzbeNPUpdYdXU148aNK4unsUv5VHhH6FbDObckTxApKIcmoFI3+ZRTYnTO5ccTRIVK42nmckiMzrn8eXffFaocnmZ2zpU3P4OoUN7k45xriSeICuZNPs655ngTk3POuUSeIJxzziXyBOGccy6RJwjnnHOJPEG4iiTpUklvSHpF0r2SittPvHMdkCcIV6keBzY1s82BN4HWd6XrXCfnCcJVJDN7zMwWxLeTCL9z4pzL4gnCORgBPJx2EM6VG2V+l6CgE5XmAtMLPuGWrQp8XkHlpll2mvM8yMxWbGkgSU8AfRO+Os/M/hWHOQ/YBjjIEnaG7N9aBzYFXmtz1MWT5rpojsfVOnlt16VUrATxvJltU/AJe7llU3ZnmGdJxwAnAUPM7NtSlVtoHlfreFz58642XEWSNBQ4G9gtn+TgXCXyaxCuUl0DrAg8LmmKpBvSDsi5clOsM4hRRZqul1s+ZXfoeTaz9dMot0g8rtbxuPJUlGsQzjnnOj5vYnLOOZeoKAlC0sGSpkpqlFSSq/KShkqaLmmGpN+UqMybJX0qqaS3PkpaS9IESdPicj6thGUvI+nfkl6OZf+uVGXH8qskvSTpwVKWm1V+WXbRkcY+10I8Jd8fW5LW/tqSNPfnlhTrDOI14CBgYpGmvwRJVcC1wN7AxsBhkjYuQdG3AENLUE6uBcCZZrYRsAPw8xLNL0ADsLuZbQEMBoZK2qFEZQOcBkwrYXm5yrWLjpLuc81JcX9syS2ks7+2JM39uVlFSRBmNs3MSvmg3HbADDN7x8zmA3cAPyp2oWY2Efii2OUklPuRmb0YX88lHDD7lahsM7N58W23+FeSC1mS+gP7ADeVorwk5dpFRwr7XHNS2R9bktb+2pI09+eWdJZrEP2AWVnvZ1MmC7jYJA0EtgSeK2GZVZKmAJ8Cj5tZqcq+Evg10Fii8lriXXQkq9j9sb3S2J+b0+bbXPPpxqCElPBZp789S9IKwN3A6Wb2danKNbOFwODY/n6vpE3NrKjtupL2BT41sxck1RS5rHy76FgAjClmLK2Nq0xU5P7YXmntz81pc4Iwsz0KGUg7zQbWynrfH/gwpVhKQlI3wsY0xszuSSMGM6uXVEto1y32hb+dgP0lDQOWAXpKus3Mjix0QS1t27GLjn0JXXSU7MBXZvtccypuf2yvctifk3SWJqbJwAaS1pHUHTgUuD/lmIpGkoDRwDQzu6LEZffJ3LkjaVlgD+CNYpdrZueYWX8zG0hYv+OLkRxaktVFx/7eRUeTKmp/bK809+eWFOs21wMlzQaqgYckPVqMcjLiRcNTgUcJF3juNLOpxSwTQNLtQB0wSNJsSccWu8xoJ+AoYPfYTcSUWLMuhTWACZJeIRwIHjezVG45TUlZdtFR6n2uOWntjy1JcX9tSZr7c7P8SWrnnHOJOksTk3POuQLzBOGccy6RJwjnnHOJPEE455xL5AnCOedcIk8QzjnnEnmCcM45l8gThHPOuUT/HzEUM65oAvkMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.subplot(2,2,1)\n",
    "plt.title('Dados Brutos',size=8)\n",
    "plt.plot(Data[:,0],Data[:,1],'k.')\n",
    "plt.axis([-1,4,-1,4])\n",
    "plt.plot(eixoMax[:,0],eixoMax[:,1]);\n",
    "plt.plot(eixoMin[:,0],eixoMin[:,1]);\n",
    "#plt.quiver(0,0,1.0,1.0)\n",
    "#plt.quiver(u[0],u[1],V[0,0],V[1,0])\n",
    "#plt.quiver(u[0],u[1],V[0,1],V[1,1])\n",
    "plt.grid(True)\n",
    "\n",
    "eixoMax=np.array([[0,0],autoVet[:,0]])\n",
    "eixoMin=np.array([[0,0],autoVet[:,1]])\n",
    "\n",
    "plt.subplot(2,2,2)\n",
    "plt.title('Dados ajustados (media zero)',size=8)\n",
    "plt.plot(DataAdjust[:,0],DataAdjust[:,1],'k.')\n",
    "plt.axis([-2,2,-2,2])\n",
    "plt.plot(eixoMax[:,0],eixoMax[:,1]);\n",
    "plt.plot(eixoMin[:,0],eixoMin[:,1]);\n",
    "plt.grid(True)\n",
    "\n",
    "\n",
    "udat=np.mean(FinalData,axis=0);\n",
    "xdat=np.array([ [-1, udat[1] ],[5, udat[1] ] ]); #reta horizontal\n",
    "ydat=np.array([ [ udat[1], -1 ],[ udat[1], 4 ] ]); #reta vertical\n",
    "plt.suptitle(\"Principal Components Analysis - PCA\")\n",
    "plt.subplot(2,2,3)\n",
    "plt.title(\"DadosRotacionados=Dados.Autovetores\",size=8)\n",
    "plt.plot(FinalData[:,0],FinalData[:,1],'k.')\n",
    "plt.plot(ydat[:,0],ydat[:,1])\n",
    "plt.plot(xdat[:,0],xdat[:,1])\n",
    "plt.axis([-1,4,-1,4])\n",
    "plt.grid(True)\n",
    "\n",
    "\n",
    "eixoMax=np.array([[-2,0],\n",
    "                  [ 2,0]]) #os pontos p1=(-2,0) e p2=(2,0)\n",
    "eixoMin=np.array([[0,-2],\n",
    "                  [0, 2]])\n",
    "plt.subplot(2,2,4)\n",
    "plt.title('Media zero rotacionados pelo maior eixo (vermelho)',size=8)\n",
    "plt.plot(FinalDataAdjust[:,0],FinalDataAdjust[:,1],'k.')\n",
    "plt.plot(eixoMin[:,0],eixoMin[:,1]);\n",
    "plt.plot(eixoMax[:,0],eixoMax[:,1]);\n",
    "\n",
    "plt.axis([-2,2,-2,2])\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matriz de Covariância dos dados transformados"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ky=np.cov(FinalDataAdjust.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.28575130e+00, 5.85951041e-17],\n",
       "       [5.85951041e-17, 4.93598137e-02]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Ky"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Observe que alguns números acima são da ordem de e-17, ou seja, praticamente zero"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "D=np.diag(autoVal)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.2857513 , 0.        ],\n",
       "       [0.        , 0.04935981]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Assim vemos que Ky é igual aos autovalores de K (ou Kx). Ou seja, as diagonais correspondem às variâncias e as não diagonais as covariâncias cov(xi,xj), que são praticamente zero. Ou seja, os novos dados estão descorrelacionados."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.2857513 , 0.        ],\n",
       "       [0.        , 0.04935981]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.round(Ky,8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}