26-PCA-使用矩阵奇异分解.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "本方法使用SVD分解，也可以使用书中的协方差分解"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "因为要最大化方差，也可以使用梯度上升的方法来进行求解X1求出第一主成分，然后去掉第一主成分的值为X2，再使用梯度上升的方法来进行求解X2求出第一主成分，依次迭代下去"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "\n",
    "def pca(X,K):\n",
    "    u,s,vT = np.linalg.svd(X,full_matrices=0)\n",
    "    return X.dot(vT.T[:, 0:K]),u,s,vT\n",
    "\n",
    "def recoverData(Z, v, K):\n",
    "    V_reduce = v[:,0:K]\n",
    "    X_rec = Z.dot(V_reduce.T)\n",
    "    return X_rec\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x1ec02de2978>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAGBJJREFUeJzt3X+MHOddx/HPJ3t2k0sLuOQosR3fpSKq6lSQNCfzI6g/sJu6pSRF8Iejawkl6BSbQKHiRyoLCoGKoiKKEEnUo0Qt+JQUGgom0KYhaQmBps25OD+cNKkT7OS4SHab/iAYtbX75Y+ZlffWu7Nzuzs3czvvlzTaneeZeeaZsfS98TzPfscRIQBAfZxVdgcAAKuLwA8ANUPgB4CaIfADQM0Q+AGgZgj8AFAzBH4AqBkCPwDUDIEfAGpmrOwOdHLeeefF1NRU2d0AgDXjwIEDX46IiTzbVjLwT01NaWFhoexuAMCaYfto3m151AMANUPgB4CaIfADQM0Q+AGgZgj8AFAzBH4AqBkCPwDkMD8vTU1JZ52VfM7Pl92j/lVyHj8AVMn8vDQ7K504kawfPZqsS9LMTHn96hd3/ADQw969p4N+04kTSflaROAHgB6eeWZl5VVH4AeAHrZsWVl51fUM/LZvtX3M9qNd6n/D9sF0edT2KdsvTeuO2H4krSP5DoA16b3vlcbHl5eNjyfla1GeO/4PS9rZrTIi3h8Rl0TEJZLeLelfI+L5lk1en9ZPD9ZVACjHzIw0NydNTkp28jk3tzYHdqUcs3oi4j7bUznbu1rSbYN0CACqaGZm7Qb6dkN7xm97XMn/DO5oKQ5Jn7J9wPbssI4FAOjfMOfx/5Skf297zHN5RCzZ/j5Jd9v+YkTc12nn9A/DrCRtWasjJgCwBgxzVs8utT3miYil9POYpI9L2tZt54iYi4jpiJiemMj1EhkAQB+GEvhtf7ek10r6h5ayc22/pPld0hWSOs4MAgCsnjzTOW+T9FlJr7C9aPta29fZvq5ls5+W9KmI+N+WspdJut/2Q5I+L+mfIuKTw+w8AKwFVcvzk2dWz9U5tvmwkmmfrWVPS/qhfjsGAKOginl++OUuABSoinl+CPwAUKAq5vkh8ANAgaqY54fADwAFqmKeHwI/ABSoinl+eAMXABSsanl+uOMHgJoh8ANAzRD4AaBmCPwARsog6RHKSK1QxjEZ3AUwMgZJj1BGaoWy0jk4IoprvU/T09OxsMAregGszNRUEjzbTU5KR44Ut2+/hnlM2wfyvuKWRz0ARsYg6RHKSK1QVjoHAj+AkTFIeoQyUiuUlc6BwA9gZAySHqGM1AplpXMg8AMYGYOkRygjtUJZ6RwY3AWAEcDgLgCgKwI/ANQMgR8AaqZn4Ld9q+1jth/tUv8621+3fTBdfqelbqftJ2wftn3DMDsOAOhPnjv+D0va2WObf4uIS9LlRkmy3ZB0k6Q3Sdoq6WrbWwfpLIDRUVSOmk2bkhkyzWXTpvzHHaRPe/ZIY2PJMcfGkvWq6pmrJyLusz3VR9vbJB2OiKclyfbtkq6S9FgfbQEYIUXlqNm0SVpaWl62tJSU//d/Zx9X6r9Pe/ZIt9xyev3UqdPrN9/c//kUJdd0zjTw3xkRr+pQ9zpJd0halLQk6dcj4pDtn5W0MyJ+Md3u7ZJ+OCKu73U8pnMCo62ovDh297qI7ONK/fdpbCwJ9u0aDenkyex9h2Ul0zmHkZ3zC5ImI+IF22+W9PeSLpLU6Z+g618Z27OSZiVpS5mvnwdQuLJy1PRz3Dx96hT0s8rLNvCsnoj4RkS8kH7/Z0nrbJ+n5H8AF7RsulnJ/wi6tTMXEdMRMT0xMTFotwBUWFk5arKOO0ifGo2VlZdt4MBv+/vt5D9YtrelbX5F0oOSLrJ9oe31knZJ2j/o8QCsfUXlqNm4Mbs867iD9Kl1nCBPeekiInORdJuk5yR9W8ld/LWSrpN0XVp/vaRDkh6S9ICkH2vZ982SnpT0lKS9vY7VXC677LIAMNr27YuYnIywk899+4bT7saNEckT/WTZuDH/cQfp0+7dEY1GcsxGI1lfTZIWImeMJVcPAIwAcvUAALoi8ANAzRD4AaBmCPwAMhWVWmEQg6RHyDqfrHareB36lncUeDUXZvUA1bBvX8T4+PJZMuPjw5uB04/du5f3p7nkmUWTdT5Z7VbxOrQTs3oADENRqRUGMUh6hKzzWVzs3u7mzdW7Du1WO2UDgBFVVmqFLIOkR8g6n273wKdOVfM6DIJn/AC6Kiu1QpZB0iNknU9Wu1W8DoMg8APoqqjUCoMYJD1C1vlktVvF6zCQvIMBq7kwuAtUR1GpFQYxSHqErPPJareK16GVGNwFgHohZQMAoCsCPwDUDIEfAGqGwA8ANUPgB0ZEGblkduxI8to0lx07ltdn5b7pt07KPtde12Gkcu70K+/0n9VcmM4JrEwZuWS2b++c22b79qQ+K/dNv3W9zrXXdVgLOXf6JaZzAvVSRk6d5E3bnUVk59SR+qs7eTL7XKXs61DF3EPDQq4eoGaqmEumn5w6eer6OddmXRWvUxl4xg+MgCrmksnKfdNvnZR9rr2uQxWvUxl6Bn7bt9o+ZvvRLvUzth9Ol/+w/UMtdUdsP2L7oG2e3QAFKSOXzPbt2eVZuW/6rZOyz7XXdRi5nDv96jUIIOk1kl4t6dEu9T8maUP6/U2SPtdSd0TSeXkHHJoLg7vAypWRS6Z9gLc5sNuUlfum37qI7HPtdR2qnnOnXxr24K7tKUl3RsSremy3If0DsSldPyJpOiK+vJI/RgzuAsDKlJmr51pJn2hZD0mfsn3AdmbSVNuzthdsLxw/fnzI3QIANA1tVo/t1ysJ/D/eUnx5RCzZ/j5Jd9v+YkTc12n/iJiTNCcld/zD6hcAYLmh3PHb/kFJH5J0VUR8pVkeEUvp5zFJH5e0bRjHAwD0b+DAb3uLpL+T9PaIeLKl/FzbL2l+l3SFpI4zgwBU1yDpEfptt8h9y2i3cnqN/kq6TdJzkr4taVHJ45zrJF2X1n9I0lclHUyXhbT85ZIeSpdDkvbmHXFmVg9QDYOkR+i33SL3LaPd1SJSNgAYhkHSI/TbbpH7ltHualnJrB4CP4Cuzjorufdt18zT063uO9/pv90i9y2j3dXCqxcBDMUg6RH6bbfIfctot4oI/AC6GiQ9Qr/tFrlvGe1WUt7BgNVcGNwFqmOQ9Aj9tlvkvmW0uxrE4C4A1AvP+AEAXRH4AaBmCPwAUDMEfgCoGQI/0EFtcrao97nW6VrUBS9bB9rMzyev+TtxIlk/evT0a/9mZsrrVxF6nWudrkWdMJ0TaLPWc7asRK9zrdO1WOuYzgkM4JlnVla+lvU61zpdizoh8ANt6pSzpde51ula1AmBH2hTp5wtvc61TteiTgj8QJuZGWluLnmObSefc3OjOZjZ61zrdC3qhMFdABgBDO4CALoi8ANAzeQK/LZvtX3M9qNd6m37z2wftv2w7Ve31F1j+0vpcs2wOg4A6E/eO/4PS9qZUf8mSRely6ykWyTJ9kslvUfSD0vaJuk9tjf021lgrduzRxobSwZKx8aS9dXYt6i0C6R7WKPyvrFF0pSkR7vUfVDS1S3rT0g6X9LVkj7YbbtuC2/gwijavTsieZ338mX37mL33bcvYnx8+X7j44O/XapXu0UdF52piDdw2Z6SdGdEvKpD3Z2S3hcR96fr90j6LUmvk3R2RPxBWv7bkv4vIv4461jM6sEoGhuTTp06s7zRkE6eLG7fotIukO6hWsqY1eMOZZFRfmYD9qztBdsLx48fH1K3gOroFLizyoe1b1FpF0j3sHYNK/AvSrqgZX2zpKWM8jNExFxETEfE9MTExJC6BVRHo7Gy8mHtW1TaBdI9rF3DCvz7Jf1cOrvnRyR9PSKek3SXpCtsb0gHda9Iy4DaaaYzzls+rH2LSrtAuoc1LM9AgKTbJD0n6dtK7uKvlXSdpOvSeku6SdJTkh6RNN2y7y9IOpwu78hzPAZ3Map2745oNJKBzkYj3+DsMPbdty9icjLCTj6HNcDaq92ijoszqYjB3dXE4C4ArAwpGwAAXRH4AaBmCPwAUDMEfgCoGQI/ANQMgR8AaobADwA1Q+AHgJoh8ANAzRD4AaBmCPwAUDMEfgCoGQI/ANQMgR8AaobADwA1Q+AHgJoh8ANAzRD4AaBmCPwAUDMEfgComVyB3/ZO20/YPmz7hg71H7B9MF2etP21lrpTLXX7h9l5AMDKjfXawHZD0k2S3iBpUdKDtvdHxGPNbSLi11q2/2VJl7Y08X8RccnwugwAGESeO/5tkg5HxNMR8S1Jt0u6KmP7qyXdNozOAQCGL0/g3yTp2Zb1xbTsDLYnJV0o6d6W4rNtL9h+wPZb++4pAGAoej7qkeQOZdFl212SPhYRp1rKtkTEku2XS7rX9iMR8dQZB7FnJc1K0pYtW3J0CwDQjzx3/IuSLmhZ3yxpqcu2u9T2mCciltLPpyV9Rsuf/7duNxcR0xExPTExkaNbAIB+5An8D0q6yPaFttcrCe5nzM6x/QpJGyR9tqVsg+0Xpd/Pk3S5pMfa9wUArJ6ej3oi4qTt6yXdJakh6daIOGT7RkkLEdH8I3C1pNsjovUx0CslfdD2d5T8kXlf62wgAMDq8/I4XQ3T09OxsLBQdjcAYM2wfSAipvNsyy93AaBmCPwAUDMEfgCoGQI/ANQMgR8AaobADwA1Q+AHgJoh8ANAzRD4AaBmCPxlmp+Xpqaks85KPufny+4RgBrIk5YZRZifl2ZnpRMnkvWjR5N1SZqZKa9fAEYed/xl2bv3dNBvOnEiKQeAAhH4y/LMMysrB4AhIfCXpdtbxnj7GICCEfjL8t73SuPjy8vGx5NyACgQgb8sMzPS3Jw0OSnZyefcHAO7AArHrJ4yzcwQ6AGsOu74AaBmCPwAUDMEfgComVyB3/ZO20/YPmz7hg71P2/7uO2D6fKLLXXX2P5SulwzzM6PNNI5AChIz8Fd2w1JN0l6g6RFSQ/a3h8Rj7Vt+tGIuL5t35dKeo+kaUkh6UC671eH0vtRRToHAAXKc8e/TdLhiHg6Ir4l6XZJV+Vs/42S7o6I59Ngf7eknf11tUZI5wCgQHkC/yZJz7asL6Zl7X7G9sO2P2b7ghXuK9uzthdsLxw/fjxHt0YY6RwAFChP4HeHsmhb/0dJUxHxg5L+RdJHVrBvUhgxFxHTETE9MTGRo1sjjHQOAAqUJ/AvSrqgZX2zpKXWDSLiKxHxzXT1LyRdlndfdEA6BwAFyhP4H5R0ke0Lba+XtEvS/tYNbJ/fsnqlpMfT73dJusL2BtsbJF2RliEL6RwAFKjnrJ6IOGn7eiUBuyHp1og4ZPtGSQsRsV/Sr9i+UtJJSc9L+vl03+dt/76SPx6SdGNEPF/AeYwe0jkAKIgjOj5yL9X09HQsLCyU3Q0AWDNsH4iI6Tzb8stdAKgZAj8A1AyBHwBqhsCfR795c3bsSGblNJcdO/K3O0iuHvL8AMgSEZVbLrvssqiMffsixscjpNPL+HhSnmX79uX7NJft23u32+8xB+kvgDVNySzLXDGWWT29TE0lSdLaTU5KR45038+dfrScishuV+rvmFL//QWwpq1kVg+vXuylqLw5/bSb55jk+QHQA8/4eykqb05Wu4Mckzw/AHog8PfSb96c7duzy7PaHSRXD3l+APSSdzBgNZdKDe5GJAOjk5MRdvKZd6C0fYC3ObCbp91+jznovgDWJDG4CwD1QsoGAEBXBH4AqBkCPwDUDIE/jz17pLGx5EdZY2PJelNWWgbSLgCoIH7A1cuePdItt5xeP3Xq9PqTT0r33LN8+3vuSYL/O94hzc5KJ04k5UePJutS7xeszM/3vy8A9MCsnl7GxpJg367R6FzeNDlJ2gUAq4ZZPcPULbhnBX1psNQJpF0AUCACfy+NxsrKm0i7AKCicgV+2zttP2H7sO0bOtS/y/Zjth+2fY/tyZa6U7YPpsv+YXZ+VTSfrXcqz0rLQNoFAFXV66e9khqSnpL0cknrJT0kaWvbNq+XNJ5+3y3poy11L+T9GXFzqVzKht27IxqNJO1Co5GsN2WlZSDtAoBVomGmbLD9o5J+NyLemK6/O/2D8Yddtr9U0p9HxOXp+gsR8eKV/DGq1OAuAKwBwx7c3STp2Zb1xbSsm2slfaJl/WzbC7YfsP3WPJ0CABQnzzz+Tq+S6vjfBNtvkzQt6bUtxVsiYsn2yyXda/uRiHiqw76zkmYlaQuDmABQmDx3/IuSLmhZ3yxpqX0j2zsk7ZV0ZUR8s1keEUvp59OSPiPp0k4HiYi5iJiOiOmJiYncJwAAWJk8gf9BSRfZvtD2ekm7JC2bnZM+1/+gkqB/rKV8g+0Xpd/Pk3S5pMeG1XkAwMr1DPwRcVLS9ZLukvS4pL+JiEO2b7R9ZbrZ+yW9WNLftk3bfKWkBdsPSfq0pPdFRDGBf5DcNlm5eCTp4ouX5+O5+OLTdevXL69bv/503fj48rr2KZqbNi2v39QydEKeHwBFyTv9ZzWXFU/n3LcvYnx8+bTK8fF8UyB3716+X3NpTtncurVz/datEevWda5bty7inHM6151zTtLuxo2d6zduHOx8BtkXwJql2r2Ba5DcNlm5eE6eTO7Ehy0iu13y/ABYofrl6hkkt02/uXiKRJ4fAAUajcA/SG6bfnPxFIk8PwAKNBqBf5DcNlm5eCRp69bO9Vu3SuvWda5bt04655zOdc3yjRs712/cSJ4fAMXKOxiwmktfuXoGyW2TlYsn4swB3q1bT9e1D/CuW3e6rn2Atzmw29Q+wLtx43DOhzw/QO2odoO7AFBz9RvcBQDkRuAHgJoh8ANAzdQn8BeVxiAr3UOvVBAAUII8aZnXvvn5ZHrmiRPJ+tGjp6drzsz03+6ePdItt5xeP3Vq+Xq3uptv7v+YADCgeszqKSqNQVa6Byk7FQQADNFKZvXU446/qDQG/aR7KDMVBACoLs/4i0pjkJXuoYqpIABAdQn8RaUxyEr30CsVBACUpB6PepoDuHv3Jo93tmxJgv4gA7vS6UHaubnkEU6jkQT21sHbrDoAKEE9BncBYMSRsgEA0BWBHwBqhsAPADWTK/Db3mn7CduHbd/Qof5Ftj+a1n/O9lRL3bvT8idsv3F4XQcA9KNn4LfdkHSTpDdJ2irpatvtr6W6VtJXI+IHJH1A0h+l+26VtEvSxZJ2Sro5bQ8AUJI8d/zbJB2OiKcj4luSbpd0Vds2V0n6SPr9Y5K223ZafntEfDMi/kvS4bQ9AEBJ8gT+TZKebVlfTMs6bhMRJyV9XdL35twXALCK8vyAyx3K2if/d9smz75JA/aspObPWl+w/USOvnVynqQv97lvnXCd8uE65cN1yqfI6zSZd8M8gX9R0gUt65slLXXZZtH2mKTvlvR8zn0lSRExJ2kuX7e7s72Q90cMdcZ1yofrlA/XKZ+qXKc8j3oelHSR7Qttr1cyWLu/bZv9kq5Jv/+spHvTt77vl7QrnfVzoaSLJH1+OF0HAPSj5x1/RJy0fb2kuyQ1JN0aEYds3yhpISL2S/pLSX9t+7CSO/1d6b6HbP+NpMcknZT0SxFBXmIAKFElc/UMwvZs+tgIGbhO+XCd8uE65VOV6zRygR8AkI2UDQBQMyMT+G3favuY7UfL7kuV2b7A9qdtP277kO13lt2nKrJ9tu3P234ovU6/V3afqsp2w/Z/2r6z7L5Ume0jth+xfdB2qXnnR+ZRj+3XSHpB0l9FxKvK7k9V2T5f0vkR8QXbL5F0QNJbI+KxkrtWKekvz8+NiBdsr5N0v6R3RsQDJXetcmy/S9K0pO+KiLeU3Z+qsn1E0nRElP57h5G544+I+5TMKEKGiHguIr6Qfv8fSY+LX1OfIRIvpKvr0mU07pKGyPZmST8p6UNl9wX5jUzgx8qlWVQvlfS5cntSTekjjIOSjkm6OyK4Tmf6U0m/Kek7ZXdkDQhJn7J9IM1UUBoCf03ZfrGkOyT9akR8o+z+VFFEnIqIS5T84nybbR4htrD9FknHIuJA2X1ZIy6PiFcryXT8S+nj6VIQ+GsofWZ9h6T5iPi7svtTdRHxNUmfUZJaHKddLunK9Nn17ZJ+wva+crtUXRGxlH4ek/RxlZipmMBfM+mg5V9Kejwi/qTs/lSV7Qnb35N+P0fSDklfLLdX1RIR746IzRExpeTX+vdGxNtK7lYl2T43nUwh2+dKukJSaTMQRybw275N0mclvcL2ou1ry+5TRV0u6e1K7s4Opsuby+5UBZ0v6dO2H1aSr+ruiGC6Ivr1Mkn3235ISb6yf4qIT5bVmZGZzgkAyGdk7vgBAPkQ+AGgZgj8AFAzBH4AqBkCPwDUDIEfAGqGwA8ANUPgB4Ca+X/CYex/cufiGwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "iris = datasets.load_iris()\n",
    "iris = datasets.load_iris()\n",
    "\n",
    "# 选取的维度为2维\n",
    "iris_x = iris.data[:, 2:]\n",
    "iris_y = iris.target\n",
    "# 我们选取两个来做二分类的为题就好了\n",
    "X = iris_x[iris_y<2]\n",
    "Y = iris_y[iris_y < 2]\n",
    "plt.scatter(X[Y==0, 0], X[Y==0, 1], color='r')\n",
    "plt.scatter(X[Y==1, 0], X[Y==1, 1], color='b')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(100, 2)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Z,u,s,vT = pca(X,1)\n",
    "X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(100, 1)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Z.shape   # 降到了一维"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_rec = recoverData(Z,vT,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X+UVOWd5/H3lwKkm2BQaR1padoxHHY1TCTTA7pkEo0/QGKQNbMbGVATE/sEdNfZzJLAtBtHxx5UJj/cjWiIYaOxgpOTGA5jjB2MuklcZWyUBCEhg8rPdoNRMWp3hC6++8e91VR3149b1dVdvz6vc/pU1/Pce+upivn05blPfa+5OyIiUjtGlXoAIiIyshT8IiI1RsEvIlJjFPwiIjVGwS8iUmMU/CIiNUbBLyJSYxT8IiI1RsEvIlJjRpd6AOlMmjTJm5ubSz0MEZGKsWXLlt+7e0OUbcsy+Jubm+ns7Cz1MEREKoaZ7Ym6raZ6RERqjIJfRKTGKPhFRGqMgl9EpMYo+EVEaoyCX0Skxij4RURqTFmu4xcRKQcbnj/A6o6ddB3qYfLEOpbPnc7CmY2lHtaQKfhFRNLY8PwBVj60jZ4jCQAOHOph5UPbACo+/DXVIyKSxuqOnX2hn9RzJMHqjp0lGlHxKPhFRNLoOtSTV3slUfCLiKQxeWJdXu2VJGfwm9k6MztoZi9k6F9uZlvDnxfMLGFmJ4Z9u81sW9inqmsiUjGWz51O3ZhYv7a6MTGWz51eohEVT5Qz/m8D8zJ1uvtqdz/b3c8GVgL/x91fT9nk/LC/ZWhDFREZOQtnNrLq8hk0TqzDgMaJday6fEbFX9iFCKt63P1nZtYc8XiLgPVDGZCISLlYOLOxKoJ+oKLN8ZtZPcG/DH6Q0uzAT8xsi5m1Fuu1RESkcMVcx/9x4KkB0zxz3L3LzE4GNpnZb9z9Z+l2Dv8wtAI0NTUVcVgiIpKqmKt6rmDANI+7d4WPB4EfArMy7ezua929xd1bGhoi3T1MREQKUJQzfjN7L/ARYElK23hglLu/Ff5+MXBLMV5PRKSclXuph5zBb2brgfOASWa2H7gJGAPg7veEm/1H4Cfu/k7KrqcAPzSz5Ot8190fLd7QRUTKTyWUeoiyqmdRhG2+TbDsM7XtJeADhQ5MRKQSZSv1UC7Br2/uiogUUSWUelDwi4gUUSWUelDwi4gUUSWUelA9fhGRIkrO41f0qh4REclPuZd60FSPiEiNUfCLiNQYBb+ISI1R8IuI1Bhd3BWRijeU2jilqKtT6lo+Cn4RqWhDqY1Tiro65VDLR1M9IlLRstXGGc59C1WK1xxIwS8iFW0otXFKUVenHGr5KPhFpKINpTZOKerqlEMtHwW/iFS0odTGKUVdnXKo5aOLuyJS0YZSG6cUdXXKoZaPufuIvVhULS0t3tnZWephiIhUDDPb4u4tUbbVVI+ISI1R8IuIlEI8Ds3NMGpU8BiPj9hL5wx+M1tnZgfN7IUM/eeZ2ZtmtjX8+VJK3zwz22lmu8xsRTEHLiJSseJxaG2FPXvAPXhsbR2x8I9ycffbwNeB+7Ns83N3vzS1wcxiwF3ARcB+4Fkz2+juOwocq4hUgeEoVzC7fRO/e+tw3/NTJoxlc9tFkV6z0PHcuGEb6zfvI+FOzIxFs6dw68IZmXeIx6GtDfbuDc7yE/2/xEV3d9C/eHF+b74AOYPf3X9mZs0FHHsWsMvdXwIwsweBywAFv0iNGo5yBQNDH+B3bx1mdvsmNrddlPU1gYLGc+OGbTzwzN6+5wn3vudpwz95ht/dHe6QGLwNBH8URkCx5vjPNbNfmtmPzeyssK0R2Jeyzf6wTURq1HCUKxgY+gPbs71moeNZv3lfXu20tR0L/WyamnJvUwTFWMf/HDDV3d82s/nABmAaYGm2zbh21MxagVaAphF68yIysiqlREKu8SQyLIPP1B7pTL6+Htrbc29XBEM+43f3P7j72+HvjwBjzGwSwRn+lJRNTwO6shxnrbu3uHtLQ0PDUIclImWo3EokFDqemKU7r83cnvFMPhYDM5g6FdauHZH5fShC8JvZn5gF79bMZoXHfA14FphmZqeb2VjgCmDjUF9PRCrXcJQrOGXC2Kzt2V6z0PEsmj2FBduf4IWvfIKXb7+Ul2+/lBdv/zgPbL0v/Q7t7cEZfar6erjvPjh6FHbvHrHQhwhTPWa2HjgPmGRm+4GbgDEA7n4P8FfAUjPrBXqAKzz4OnCvmV0PdAAxYJ27bx+WdyEiFWE4yhVsbrso66qeKK+Z73hufedXJH70VWJ+tK8thnPuo9+DZSfBmjX9d0iGenJVT1NT8MdgBMM+lUo2iIjkq7k5WHufTiwGvb0jOhxQyQYRkeGV7WJtpqWaZUTBLyKSTraSCtlWHsZimfvKhIJfRCRVPA6TJsGSJZlLKrS3Zw741taRG2uBFPwiItA/8F97bXB/sqQCBBdl77sPxo8/1j9qFCxdOvjCbhnSjVhEJK3hqKkzFHnXxkmR7b08fcknmfXo9xhF+m+dpvK9e49ts3hxyVblDJWCX0QGGY6aOkORd22cFNneyyl/93nOefR7OQM/qev4Bp59/kBJ/wAWg6Z6RGSQ4aipMxR518ZJkfpeFmx/gi13LmLHrZdw2QdPyyv0u0cfx+1/eWXJPoNi0hm/iAxSipo62eRdGydF16EeFmx/gr//6VpO6HkrctAnOfBG3QT+/oJWNp51Plaiz6CYFPwiMsjkiXUcSBNww1lTJ5uYWdqQz1gbJ2nZMnbdfXek+fuBHHh93ARuvjAI/KRSfQbFpOAXkUGWz53eb14chl5TZygWzZ7Sb44/tT2txkboCmpC5ruq3oF3x9Wz7cbbueqP7yubz6CYNMcvIoMsnNnIqstn0DixDgMaJ9ax6vIZJbuoeevCGSw5p6nvDD9mxpJzmtJf2E0J/ag8/Om1UTwz7z8zrucd/qLt+rL6DIpJtXpEpPKl3tYw30wbOxbWravYpZlJ+dTq0VSPiFS2gbc1zMdJJ8Gdd1Z86OdLwS8ilWXZsuCmJYlEUDZh3Lj8Qv8974F77qm5sE+l4BeRyhCPw+c+B2+/fawtkYB33ol+jAopqTDcFPwiVWCkyyss/ubTPPXi633P55xxIvFrz+17nq28Qt597/wKbrgBf+21vJdkJnU3/AkXfT4efD63PV7y8hOlpuAXqXAjXV5hYOgDPPXi6yz+5tPErz03a3kFIK++P3zrPg5vuoux7/4xa+g7/dfp94w5jhdu+if+ou36Y59P+L2EUpefKAcKfpEKl628wnAE28DQH9heSHmFZN/NHWtY/MtHiflREjaKP44ew9gj7+Yc01EzXpkwicl/+D1dx0/ijg9fxZbYmTzFyH8+lUDBL1LhqqG8QsKdmzvWcNXWR/rO3Ef7UcZHCH0HHvjAJdw0d1m/9mRphXL7fMqBvsAlUuEylRAoZXmFTO2Z+hbueLJf6Cflmt5JAA/NXjAo9OHY+y+3z6cc5Ax+M1tnZgfN7IUM/YvN7Ffhz/81sw+k9O02s21mttXM9I0skWGwfO506sb0L0wwnKUF5pxxYtb2TGUUFs2e0td3c8cadt2xgJdvv5QX7/g4X374yxlDfuC/Exx4bdwEli9czr88t5/Y3Wuyvv+R/nwqQZSpnm8DXwfuz9D/MvARd3/DzC4B1gKzU/rPd/ffD2mUIpJRcp56pFb1xK89N+uqnuQKnbQrd+JxbvzatRz3bk9f0MdyfNPWRo2CKVPwvXv53XsbWPWhK+mcM3/Qe8z0/kf686kEkUo2mFkz8LC7vz/HdicAL7h7Y/h8N9CSb/CrZINIFSr0G7Zaex9JPiUbij3H/xngxynPHfiJmW0xs/K/A7GIDJ+2NoV+mSha8JvZ+QTB/8WU5jnu/kHgEuA6M/twlv1bzazTzDpfffXVYg1LREpl2TIYPRrMgsc9e6LvW18PDzyg0B8mRQl+M/sz4F7gMnfvuz29u3eFjweBHwKzMh3D3de6e4u7tzQ0NBRjWCJSKsuWwd13ByUV4NhjNuPHB38kpk4NavHUcC2d4Tbkdfxm1gQ8BFzp7r9NaR8PjHL3t8LfLwZuGerricjIyFYGIrXv6pef4gs/v5/6/9cFTU3Q3h4EdwQOYIZ97nM5z+6HqyzFSJe7KAc5g9/M1gPnAZPMbD9wEzAGwN3vAb4EnASssWCNbm94geEU4Idh22jgu+7+6DC8BxEpsmxlIIC+vgXbn+B/PPzlY3e52rMHrrwyY018B7reezKnvvlq3zdsN519YXCDkwLHM5SQHulyF+VCN2IRkUHm3PZ42nvuNoZfejoQ3sD8aw9/Oa/54l4bxfu+sDHtcZ9a8dGCxpNtv1yG67iloBuxiMiQpCtnsGD7E/z9T9dyQs9bfW35Vsv87gfmRX69KP1DLbtQq+UcVLJBRAYZWM5gwfYn+Kcf38mJPW9h0PeT1dKlwY1SIHhcupRvfPJvI71e1P6hll2o1XIOCn4RGSRZ5mDB9if4xd2f5s6Hv8zYRG9+B1mzBnp7g/n+3l5Ys6bg8gnDVXahVss5aKpHRAZZuONJLvlf1zP2zUOF3fzkggvSH7fA8gnDVXahVss56OKuiPQ3lJuXQxD6jz1W3DFJTqUs2SAila6Q0gpTpwbftHVX6FcABb9IrYrHobkZRo0KHuPxoH3v3mx7HXPSScfCfvdufdO2gmiOX6QWDZzO2bMneA7Bt28z1dWpr1c5hSqgM36RWpRuOqe7O2hvbw8CfqCTTlLoVwmd8YtkUNU1XDJM5/jevVgY7N3Lv8i4V7roOn4S9877LGevuK563n+NU/CLpFHxNVzicbjhBngtLJZ70klw553HztYzTOd0Hd/As88fgDPPY+Vn7u17/wB1lfT+JStN9YiksbpjZ7/QA+g5kmB1x84SjSgP8Th8+tPHQh+C36+55tgF3PZ2esYc12+37tHHcftfXsnqjp2V/f4lJwW/SBoVXcOlrQ2OHBncfvhw0AeweDEr5l7P/uMbOIqx//gGVsy7no1nnU/XoZ7Kfv+Sk6Z6RNKYPLEubdXGiqjhkm05Zkpf55z5fOis8wdtMjmlAmemPqlsOuMXSaOia7g0NUXqy/YeK/r9S0464xdJo6JruLS3B3P8A6d7xo4N+kJR3mNFvn/JSbV6RCrBsmXBGvpEIihx3Nqa/VaFuVb1SNXRjVhEqknyxuVJicSx55nCf/FihbxkpDl+kXKX6cblEW9oLjKQgl+k3CUS+bWL5BBpqsfM1gGXAgfd/f1p+g24E5gPdAOfcvfnwr6rgRvDTW919/uKMXCRSnTjhm2s37yPhDsxMxbNnsKtC2dk3ykWSx/ysdjgtjSGo/RErmNWdbmLKhD1jP/bQPq7JAcuAaaFP63A3QBmdiJwEzAbmAXcZGYnFDpYkUp244ZtPPDMXhLhgoqEOw88s5cbN2zLut/TF32CgUswPGzPJVl64sChHpxjpSc2PH+gsDcR4ZjD8ZpSXJGC391/BryeZZPLgPs98Aww0cxOBeYCm9z9dXd/A9hE9j8gIlVr/eZ93Nyxhl13LODl2y9l1x0LuLljDes378u635Kzr+b+s+fTa6NwoNdGcf/Z81ly9tU5X3M4Si/kOqbKPZS/Yq3qaQRS/+vdH7Zlah/EzFoJ/rVAU7YvoIhUqC89ehdXbX2k7x62o/0oV219JHiyan7G/RLu3DR3GTfNXda/I8JS7OEovZDrmCr3UP6KdXE33f2YPUv74Eb3te7e4u4tDQ0NRRqWSPlY/MtHB/0fwsL2bGKW/nbnmdpTZSqxMJTSC7mOORyvKcVVrODfD0xJeX4a0JWlXaTmxPxoXu1Ji2ZPyas91XCUXsh1TJV7KH/FCv6NwFUWOAd4091fATqAi83shPCi7sVhm0jNsQyrcDK1J926cAZLzmnqO8OPmbHknKbcq4EIyjKsunwGjRPrMKBxYh2rLp8xpBU2uY45HK8pxRWpZIOZrQfOAyYBvyNYqTMGwN3vCZdzfp3gwm038Gl37wz3vQb4u/BQ7e7+v3O9nko2SFUa+A3cpKVLs5dfEImg6CUb3H1Rjn4HrsvQtw5YF+V1RKpaMtzzqbkjMgxUq0dkJK1Zo6CXklPJBhGRGqPgF4kqHofmZhg1KnhM3r9WpMJoqkcking8mI/v7g6e79kTPAeVP5aKozN+kSja2o6FflJ397Gbl4tUEAW/SBSZbmCe7cbmImVKwS8SRab6UaorJRVIwS8SRXs71Nf3b6uv73fzcpFKoeAXiWLx4uCLV1OnglnwuHatLuxKRdKqHpGodANzqRI64xcRqTEKfhGRGqPgFxGpMQp+EZEao+AXEakxCn4RkRqj4BcRqTEKfhGRGqPgFxGpMZGC38zmmdlOM9tlZivS9H/VzLaGP781s0MpfYmUvo3FHLyIiOQvZ8kGM4sBdwEXAfuBZ81so7vvSG7j7v8tZfv/AsxMOUSPu59dvCGLiMhQRDnjnwXscveX3P0w8CBwWZbtFwHrizE4EREpvijB3wjsS3m+P2wbxMymAqcDj6c0jzOzTjN7xswWFjxSEREpiijVOS1Nm2fY9grg++6eSGlrcvcuM/tT4HEz2+buLw56EbNWoBWgSTe3EBEZNlHO+PcDU1KenwZ0Zdj2CgZM87h7V/j4EvAk/ef/U7db6+4t7t7S0NAQYVgiIlKIKMH/LDDNzE43s7EE4T5odY6ZTQdOAJ5OaTvBzI4Lf58EzAF2DNxXRERGTs6pHnfvNbPrgQ4gBqxz9+1mdgvQ6e7JPwKLgAfdPXUa6N8D3zCzowR/ZG5LXQ0kIiIjz/rndHloaWnxzs7OUg9DRKRimNkWd2+Jsq2+uSsiUmMU/CIiNUbBLyJSYxT8IiI1RsEvIlJjFPwiIjVGwS8iUmMU/CIiNUbBLyJSY6JU55QRsOH5A6zu2EnXoR4mT6xj+dzpLJyZtvq1iMiQKPjLwIbnD7DyoW30HAmqWR841MPKh7YBKPxFpOg01VMGVnfs7Av9pJ4jCVZ37CzRiESkmin4y0DXoZ682kVEhkLBXwYmT6zLq11EZCgU/GVg+dzp1I2J9WurGxNj+dzpJRqRiFQzBX8ZWDizkVWXz6BxYh0GNE6sY9XlMwq7sBuPQ3MzjBoVPMbjRR6tiFQ6reopEwtnNg59BU88Dq2t0N0dPN+zJ3gOsHjx0I4tIlVDZ/yVLvUM/+qrj4V+Unc3tLWVZGgiUp50xl/JBp7hJxLpt9u7d+TGJCJlT2f8laytbfAZfjpNTcM/FhGpGJGC38zmmdlOM9tlZivS9H/KzF41s63hz2dT+q42s38Lf64u5uBrRqYLtlHO5Ovrob19OEcnIhUm51SPmcWAu4CLgP3As2a20d13DNj0n939+gH7ngjcBLQADmwJ932jKKOvdvE43HAD/tprWLJtzx56P3tt8D9cU1NwAXegWAyOHg3629t1YVdE+olyxj8L2OXuL7n7YeBB4LKIx58LbHL318Ow3wTMK2yoNSIeh0mTwAyWLIHU0A+N/mMP3cu/GIR6fX3/zvp6uO++IPh371boi8ggUYK/EdiX8nx/2DbQJ8zsV2b2fTObkue+ArBsWV/Y5zLula4g1NeuhalTgz8UU6cGzxX2IpJFlOAfeMIJwbRNqn8Bmt39z4DHgPvy2DfY0KzVzDrNrPPVV1+NMKwqkZy/N4O77468W9fxk4JfFi8Ozux1hi8iEUUJ/v3AlJTnpwFdqRu4+2vu/m749JvAn0fdN+UYa929xd1bGhoaooy98iWXY6abp8+ie/Rx3Dvvs7k3FBFJI8o6/meBaWZ2OnAAuAL469QNzOxUd38lfLoA+HX4ewfwj2Z2Qvj8YmDlkEddyRoboSvt376cHHijbgL/OPdzfGjFdcUdl4jUjJzB7+69ZnY9QYjHgHXuvt3MbgE63X0j8F/NbAHQC7wOfCrc93Uz+weCPx4At7j768PwPsrfhRfCT39a0K4OvFk3gZsuaKVzznzdnUtEhsTc0065l1RLS4t3dnaWehjFEY/DNdfA4cOF7T9+PHzjG5q7F5GszGyLu7dE2VYlG4ZLfT30FHAjFTNwD1boaA2+iAwDlWwotmXLgvDON/STyzG/850g+LVCR0SGic74i2EI8/cATJ4MBw4UbzwiIlko+PO04fkDrO7YSdehHn59x2Uc54m0X1bIJXllxcLQTz3u5Il1/S7gZuuLOtZ89hOR6qbgz8OG5w/wi1vupONH/5PxR4KvLeQb+snA/3nTB7hqUTtzzjiR//T8AVY+tI2eI0FZ5QOHelj50La+fTL1ZQvxDVmOqfAXqW1a1ZOHze/7c2a9+FzBZ/gO/M2lf8vGs87v19c4sY4DhwZfE2gMb7aeqe+pFR/N+Hpzbnu8oP1EpDJpVU+xnXUW7NjBLPI/w096JzaW9//3h9L2daUJ6GztufoKPaaI1Aat6sklDH0oIPRjMVi6FNwzhj7A5PDMPl17tr5sCt1PRKqfgj+XHQNvOxDBBRcESzJ7e2HNGgDmnHFi2k3nnHEiy+dOp25MrF973ZgYy+dOz9qXTaH7iUj1U/AnLVsGo0cH6+lHjw6e52vMmCDwH3tsUFf82nMHhf+cM04kfu25LJzZyKrLZ9A4sQ4jmIdfdfkMFs5szNqXTaH7iUj108VdCEI+XUnkpUujlUpeurTvzF5EpBTyubirM34Ibl6Sqf3MMzPvF87fK/RFpJJoVQ9AIpG5ffv2fhd4geCPwfbtIzM2EZEiq40z/uRdrkaNCh7j8f79sVi6vY61b98enNknfxT6IlLBqv+MP3mXq+7u4PmePcFzOFYErbU1/Vx+crsUN27YxvrN+0i4EzNj0ewp3LpwRnC4bz7NUy8eu91A8uJtkkoviEg5qP6Lu83N6W9tOHVqUAEzadmyYE4/kQjO9FtbB83d37hhGw88s3fQoZac08TLr77dL/STkuE/sIQCBMsrc620KXQ/EakturibOrWT6X62ewcE+Jo1wbr7AevvU63fvC/todZv3pc29IG+9tUdO/uFN0DPkQSrO3ZmfSuF7icikkn1TfUMnNrJpKkp70MnMvzrKFN7qkJLKKj0gogUW/Wd8be15Q79+vrg7lZ5iln6og2Z2lOp9IKIlIvqC/6BUzipkne5Wru2oLtbLZo9JWN7tpIMUHgJBZVeEJFiizTVY2bzgDuBGHCvu982oP/zwGeBXuBV4Bp33xP2JYBkcfm97r6gSGNPr6kp2sXcAiRX7xSyqid5ITbf1TmF7icikknOVT1mFgN+C1wE7AeeBRa5+46Ubc4HNrt7t5ktBc5z90+GfW+7+3vyGdSQVvWkm+Ovry/4LF9EpBIUe1XPLGCXu7/k7oeBB4HLUjdw9yfcPZm0zwCn5TPgolq8OAj5qVOHPLUjIlKNokz1NAKp6xj3A7OzbP8Z4Mcpz8eZWSfBNNBt7r4h3U5m1gq0AjQVsOKmn8WLFfQiIhlECf50S1bSzg+Z2RKgBfhISnOTu3eZ2Z8Cj5vZNnd/cdAB3dcCayGY6okwLhERKUCUqZ79QOpyltOAroEbmdmFQBuwwN3fTba7e1f4+BLwJDBzCOMVEZEhinLG/ywwzcxOBw4AVwB/nbqBmc0EvgHMc/eDKe0nAN3u/q6ZTQLmAHcUa/C5FFrjJls9nou+8iT/dvCdvm2nnTyeTZ8/D4D3rfwRvSn/VhltsGvVx/qe/7u2R/hj4tgG42LGb9rnAzC7fRO/e+twX98pE8ayue2iIb0P1fgRkXQi1eoxs/nA1wiWc65z93YzuwXodPeNZvYYMAN4Jdxlr7svMLP/QPAH4SjBvy6+5u7fyvV6xajVU2iNm2z1eDa/9Fq/0E+advJ4Xn71nX6hn5QM/4GhnzQuZry3fky/0E86ZcJYVs4/UzV+RCSnfFb1VG2Rtjm3Pc6BNGUNGifW8dSKj2bc74yVj6QtwRAzi1SaIZ3dt32M5hU/Kmjfxol1Bb2PQt+/iFQmFWmj8Bo3Q6nHMxxU40dEiq1qg7/QGjdDqcczHFTjR0SKrWqDv9AaN9nq8Uw7eXzavmknj2d0hr8LyfZxsfQbjIsZp0wYm7bvlAljVeNHRIqu+soyhwqtcZOrHk+hq3p+0z6/4FU9hbwP1fgRkUyq9uKuiEgt0cVdERHJSMEvIlJjFPwiIjWmai/u5jIc5QyylXqI0i8iMhJqMvgHljM4cKiHlQ8FNwkrNPwHlnpIuPc9v3XhjJz9IiIjpSanelZ37OxXwwag50iC1R07Cz7m+s37srbn6hcRGSk1GfzDUc4gV6mHcisFISK1qyaDfzjKGeQq9VBupSBEpHbVZPAPRzmDbKUeovSLiIyUmry4OxzlDHKVesjVLyIyUlSyQUSkCqhkg4iIZKTgFxGpMQp+EZEaEyn4zWyeme00s11mtiJN/3Fm9s9h/2Yza07pWxm27zSzucUbuoiIFCJn8JtZDLgLuAQ4E1hkZmcO2OwzwBvu/j7gq8Dt4b5nAlcAZwHzgDXh8UREpESinPHPAna5+0vufhh4ELhswDaXAfeFv38fuMDMLGx/0N3fdfeXgV3h8UREpESiBH8jkFpQZn/YlnYbd+8F3gROiriviIiMoChf4EpXU2Dg4v9M20TZNziAWSvQGj5928wKr5g22CTg90U8XrXR55OdPp/M9NlkN5Kfz9SoG0YJ/v1Aal2B04CuDNvsN7PRwHuB1yPuC4C7rwXWRht2fsysM+oXG2qRPp/s9Plkps8mu3L9fKJM9TwLTDOz081sLMHF2o0DttkIXB3+/lfA4x58JXgjcEW46ud0YBrwr8UZuoiIFCLnGb+795rZ9UAHEAPWuft2M7sF6HT3jcC3gO+Y2S6CM/0rwn23m9n3gB1AL3CduyfSvpCIiIyIsqzVU2xm1hpOJUka+nyy0+eTmT6b7Mr186mJ4BcRkWNUskFEpMZUdfCb2TozO2hmL5R6LOXGzKaY2RNm9msz225mN5R6TOXEzMaZ2b+a2S/Dz+fmUo+pHJlZzMyeN7OHSz2WcmNmu81sm5ltNbOyqjNf1VM9ZvZh4G3gfnd/f6nHU07M7FTgVHd/zswmAFuAhe6+o8RDKwvhN8/Hu/vbZjYG+AVwg7s/U+I0K5tZAAABrElEQVShlRUz+zzQAhzv7peWejzlxMx2Ay3uXnbfc6jqM353/xnBKiMZwN1fcffnwt/fAn6NvlXdxwNvh0/HhD/Ve5ZUADM7DfgYcG+pxyL5qergl2jCaqozgc2lHUl5CacxtgIHgU3urs+nv68BXwCOlnogZcqBn5jZlrAyQdlQ8Nc4M3sP8APgb9z9D6UeTzlx94S7n03wjfNZZqbpwpCZXQocdPctpR5LGZvj7h8kqGx8XTj1XBYU/DUsnLv+ARB394dKPZ5y5e6HgCcJSotLYA6wIJzHfhD4qJk9UNohlRd37wofDwI/pIwqEyv4a1R48fJbwK/d/SulHk+5MbMGM5sY/l4HXAj8prSjKh/uvtLdT3P3ZoJv6j/u7ktKPKyyYWbjw0UTmNl44GKgbFYXVnXwm9l64GlgupntN7PPlHpMZWQOcCXBmdrW8Gd+qQdVRk4FnjCzXxHUq9rk7lqyKFGdAvzCzH5JUJ/sR+7+aInH1Keql3OKiMhgVX3GLyIigyn4RURqjIJfRKTGKPhFRGqMgl9EpMYo+EVEaoyCX0Skxij4RURqzP8HtV23RToTlAYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(X[:,0],X[:,1])\n",
    "plt.scatter(X_rec[:, 0], X_rec[:, 1],c='red')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "END"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "实验一下波士顿房价"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(506, 13)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "boston = datasets.load_boston()\n",
    "X1 = boston.data\n",
    "X1.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(506, 8)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 我们降到8维\n",
    "Z1,u1,s1,vT1 = pca(X1,8)\n",
    "Z1.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(506, 13)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_rec1 = recoverData(Z1,vT1,8)\n",
    "X_rec1.shape"
   ]
  }
 ],
 "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}