Fisherface+Eigenface.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import scipy\n",
    "import numpy as np\n",
    "import numpy.linalg as LA\n",
    "\n",
    "%matplotlib inline\n",
    "import pylab\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn import neighbors\n",
    "import pdb\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 数据加载"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 202,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import scipy.io as sio\n",
    "\n",
    "def load_Yale_data(randomsplit):\n",
    "    #matlab文件名\n",
    "    matfn=u'face_data/Yale/Yale_32x32.mat'\n",
    "    data=sio.loadmat(matfn)\n",
    "    face = data['fea']\n",
    "    label = data['gnd'] - 1 ### python 从0 开始\n",
    "    face = face  / 255.0\n",
    "\n",
    "    #### 显示第一张图像\n",
    "#     f1 = face[5].reshape((32,32)).T\n",
    "#     pylab.imshow(f1, cmap=plt.cm.gray), pylab.show()\n",
    "\n",
    "    ###  p = 6, randomsplit = 1\n",
    "    matfn = u'face_data/Yale/6Train/'+str(randomsplit)+'.mat'\n",
    "    index = sio.loadmat(matfn)\n",
    "    testIdx = index['testIdx'] -1    ### python 从0 开始\n",
    "    trainIdx = index['trainIdx'] -1  ### python 从0 开始\n",
    "\n",
    "    train_face = face[trainIdx]\n",
    "    test_face = face[testIdx]\n",
    "    train_face = np.squeeze(train_face)  ## 将维度为1去掉\n",
    "    test_face = np.squeeze(test_face)    ## 将维度为1去掉\n",
    "    train_face = train_face.T ### 将数据维度变为: 样本维度 x 样本个数\n",
    "    test_face = test_face.T ### 将数据维度变为: 样本维度 x 样本个数\n",
    "    train_label = np.squeeze(label[trainIdx])\n",
    "    test_label = np.squeeze(label[testIdx])\n",
    "    return train_face, test_face, train_label, test_label\n",
    "\n",
    "def load_ORL_data(randomsplit):\n",
    "    #matlab文件名\n",
    "    matfn=u'face_data/ORL/ORL_32x32.mat'\n",
    "    data=sio.loadmat(matfn)\n",
    "    face = data['fea']\n",
    "    label = data['gnd'] - 1 ### python 从0 开始\n",
    "    face = face  / 255.0\n",
    "\n",
    "    #### 显示第一张图像\n",
    "#     f1 = face[5].reshape((32,32)).T\n",
    "#     pylab.imshow(f1, cmap=plt.cm.gray), pylab.show()\n",
    "\n",
    "    ###  p = 6, randomsplit = 1\n",
    "    matfn = u'face_data/ORL/6Train/'+str(randomsplit)+'.mat'\n",
    "    index = sio.loadmat(matfn)\n",
    "    testIdx = index['testIdx'] -1    ### python 从0 开始\n",
    "    trainIdx = index['trainIdx'] -1  ### python 从0 开始\n",
    "\n",
    "    train_face = face[trainIdx]\n",
    "    test_face = face[testIdx]\n",
    "    train_face = np.squeeze(train_face)  ## 将维度为1去掉\n",
    "    test_face = np.squeeze(test_face)    ## 将维度为1去掉\n",
    "    train_face = train_face.T ### 将数据维度变为: 样本维度 x 样本个数\n",
    "    test_face = test_face.T ### 将数据维度变为: 样本维度 x 样本个数\n",
    "    train_label = np.squeeze(label[trainIdx])\n",
    "    test_label = np.squeeze(label[testIdx])\n",
    "    return train_face, test_face, train_label, test_label"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# BaseLine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 203,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "39.28 4.03724878819\n"
     ]
    }
   ],
   "source": [
    "acc_list = []\n",
    "for i in xrange(1,51):\n",
    "    train_face, test_face, train_label, test_label = load_Yale_data(randomsplit=i)\n",
    "    mn = np.mean(train_face, keepdims=True, axis=1)\n",
    "    train_face = train_face - mn\n",
    "    test_face = test_face - mn\n",
    "    clf = neighbors.KNeighborsClassifier(algorithm='auto',n_neighbors=1, weights= 'distance')  \n",
    "    clf.fit(train_face.T, train_label) \n",
    "    \"\"\"测试准确率\"\"\"\n",
    "    acc = np.float(sum(test_label==clf.predict(test_face.T))) / len(test_label)\n",
    "    acc_list.append(acc)\n",
    "print 100.0-np.mean(acc_list)*100, np.std(acc_list)*100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# EigenFace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 204,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def PCA(data, K=2):  \n",
    "    \"\"\"\n",
    "    data : M x N matrix (M dimensions, N samples)\n",
    "    K: The final reduced dimension\n",
    "    \"\"\"\n",
    "#     M, N = data.shape\n",
    "#     mn = np.mean(data, keepdims=True, axis=1)\n",
    "# #     data = data - mn\n",
    "#     covariance = np.matmul(data, np.transpose(data)) / (N-1)\n",
    "#     eigenValues, eigenVectors = LA.eig(covariance)\n",
    "#     eigValInd=np.argsort(eigenValues)\n",
    "#     eigValInd=eigValInd[::-1][0:K]\n",
    "#     principalEigVects=eigenVectors[:,eigValInd]\n",
    "#     lowDimData=np.matmul(np.transpose(data), principalEigVects)\n",
    "    M, N = data.shape\n",
    "    mn = np.mean(data, keepdims=True, axis=1)\n",
    "#     data = data - mn\n",
    "#     print data.shape\n",
    "    U, S, V = LA.svd(data, full_matrices=True)\n",
    "    principalEigVects= np.transpose(U)[0:K,:]\n",
    "    \n",
    "    return principalEigVects.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 205,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "39.92 4.5208455453\n"
     ]
    }
   ],
   "source": [
    "acc_list = []\n",
    "for i in xrange(1,51):\n",
    "    train_face, test_face, train_label, test_label = load_Yale_data(randomsplit=i)\n",
    "    mn = np.mean(train_face, keepdims=True, axis=1)\n",
    "    train_face = train_face - mn\n",
    "    test_face = test_face - mn\n",
    "    Eigen_faces = PCA(train_face, K=32)\n",
    "    ##train数据\n",
    "    lowDim_TrainData = np.matmul(Eigen_faces.T, train_face)\n",
    "    ##test数据\n",
    "    lowDim_TestData = np.matmul(Eigen_faces.T, test_face)\n",
    "    clf = neighbors.KNeighborsClassifier(algorithm='auto',n_neighbors=1, weights= 'distance')  \n",
    "    clf.fit(lowDim_TrainData.T, train_label) \n",
    "    \"\"\"测试准确率\"\"\"\n",
    "    acc = np.float(sum(test_label==clf.predict(lowDim_TestData.T))) / len(test_label)\n",
    "    acc_list.append(acc)\n",
    "print 100.0-np.mean(acc_list)*100, np.std(acc_list)*100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 206,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAA9CAYAAAC0qDI5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztncmOXdlxrv9z9un7LltmFimxGkswBEiCRx7ZI3mghzDg\noR/CgN/AIz+DAQGG38IwJNmDkkpFVrHLYubJPH3fnzsgvsg4yb1p2gCBi3t3AAmyWJl51l4rVsQf\nf/xrbSm22GKLLbbYYosttthiiy222GKLLbbYYosttthiiy222GKLLbbYYosttthiiy222GKLLbbY\nYosttthiiy222GKLLbbYYosttthiiy222P6fssSn+sX/+I//uJ/NZhoOh+r3++r1eppOp5rNZprP\n59rv98rlcioUCjo6OtLx8bHK5bLS6bRWq5W2262SyaSy2azy+bxyuZzS6bRyuZx95fN55fN5SdJ8\nPtdwONTf/u3fvvdMv/nNb/ZBEGi5XGo0Gmk8HmuxWGi9Xmu322m/3yuVSmm73Wq5XKrdbmuxWGi/\n3yuZTCqVStk4yuXye2PJZDLa7/fK5/MqFAoKgkC//vWv3xvHr371q306nValUlE+n1c2m1W5XFal\nUlEQBFqv19put9put9rv96rX60qlUkok3v2qIAgUBIF9niSt12ttNhut12tNp1ONRiNNJhNNp1P1\nej3927/923vj+PWvf73/j//4D3W7XQVBoGKxaM+RTCYlSZlMRtlsVpVKRYlEQrvdzsaQTCaVTCaV\nyWSUTqdtPIx7u93aWufzeVWrVf3rv/7re+P4+7//+/2jR49ULBZVrVYVBIFSqZR9JZNJ7XY7+5x8\nPn/wOel02sa72+203W4P/vRzdXd3p36/r7/7u797bxx/8zd/sy8WiwqCQIlEQqlUyn53KpVSvV5X\nq9VSOp2WJI3HY00mEy2XS+12O202Gy0WC00mE/V6PS0WC83nc61WKwVBcPBZi8VChUJBv/3tb98b\nx1//9V/vHz16ZPuiVCrZ/OKfi8VCw+FQvV5Pk8lEq9VKiUTCxhYEgbbbrRKJhPnNfr83f0kkEubv\njUZD//zP/xzqH7/61a/05MkTNRoN7XY7LZdLJZNJVSoVJZNJTadT3dzcqNvtajabaTAY2GfNZjNt\nt1vNZjNbv+12a2Nlfvn+ZDKpf/qnf3pvHH/1V3+1z+Vy5g/MRSaTUaVSUaFQUD6fV61WUzabtTma\nzWb2jOl0Wvv9XqPRSNPpVJPJROv1Wsvl0vbbbrfTaDRSNpvVv/zLv7w3jr/8y7/c53I5XVxc6Be/\n+IWePn2qfD5v8SOXy6lYLGq73arf72s4HGq1Wimfz6tUKimbzdrzTyYTLRYLTadTrddrrVYrrVYr\n21/stX/4h3/4qNic+phv+t8YgSeVSmm/32u329kgWbxkMqn9fm8PJEm5XM6c0Dv/er22iZ/NZioU\nCtpsNkqlUspkMrb5wmw2mymXy2m5XGqxWGi73SqbzapYLCqVSmm321mw6Ha7KhQKkqTtdqt0Oq1S\nqaRisahCoSA2Os4UBIFtYhyCzRI2J5K02WxsQxBwpHfOjFOkUilLZP7n/WcS0JlXAk4ymbSNGmbT\n6dSerVwuq16vq1qtKpfLWWAmmG63W5t71oQNVa/Xlc1mLUCQfBi/JBtbmOXzeUsGjJs/5/O51uu1\nxuOxBQGSZxAE2u/3FqQ2m41Wq5XNAV9s/kKhoEKhoNFoFDmWXC6nbDZrG2q1WimVSimXy2kwGGg+\nn1tQ32639pnr9fpgQ/o12m63tg/wPeYpzDxgKBQKqtfrqlQqllyDINBms9FsNtN4PFa/39dgMNB2\nu7Xxs99Ys/1+r9VqZfvGz1WUfzQaDQMnm81Gb9++1fX1tQVg1nu9XiuZTJoPTiYTDYdDA0DMRyKR\nMEBQKpWUz+cNLGw2m8g1YZ1JxIVCQeVyWZlMRo1GQxcXF3r06JEqlYqWy6X2+72Wy6Umk4l2u50K\nhYI9w3A41N3dnSaTiWazmQHJxWKh1WplQCDMUqmUgiBQpVJRo9FQvV7XdrvVcDi0ee12uxbPNpuN\nBWLWYDqd6vb2VqPRSIlEQslk0pI0a8JcfGhO3hvbR3/n/9Bwqvl8bosMcvULyobFOUAT6XTagvJq\ntbINvd/vbTMvl0tDZx8K4GxGsi7Bu1wuK5vNKpPJaL1eazabmaOMx2Ot12vL9qVSyQL4er2WJBvv\nfr/XZDLRZrOxpBJmbGqQWTKZ1Hq9Vq/X02azseBPdTKdTi1Is7F9QMhkMhYgCSoEAb4/zPisfD6v\nRqOhk5MTtVotmwvWaT6fG5rnOfP5vAW2YrFoSESSBafdbqd0Oq31eq35fB7pI1Q/yWRS6XTa1sgj\n2NVqZQiq3++rUCgok8nYz7O+y+XSAAMoLZvN2nqXSiXlcrnQcRDs0+m0ptPpwbwvl0t7fjY580TQ\nGo1GWi6XVlVls1nzRcZIAF4sFpEbFKBTLpctUABi8LVSqaRqtap6vW7+u1gsrJpjbhaLhZbLpZbL\npYGLYrGo2Wymfr9/AKge2nw+N9R/d3enr7/+Wq9fv9ZisVAikbC9Vq/XbU8QjHu9nt68eaN+v69c\nLmeVVTqd1na7VSaT0XA4VK1Ws+o5KpHMZjNDsHwOCQkfC4JA2WxWuVxOi8XC/r1QKKjRaFhFXygU\nDJzsdjt7FmIQyS3MxuOxoelGo6FsNqvBYGDVayKR0HQ6teqqWCxaDAOk8bXf762aJ5ER9Gezmfnx\nx9onC+BsAJ9B2QSeBvFlz3a71Xw+t0UgWBEYfVm02+0ssOMIUQ8OcifT+qDMGFOplE1qs9m0BU0k\nEjbRBAfGyrPxTDhQlENCk/Bs8/ncqCWSSyKRsCDlUT6OKknL5VLD4VC73U6JRMICaTqdVjab1Waz\nsaQZZlQc+XxelUrFNhPJM5/PW7AhcLAJCYRsGqgM1leSrSPBI2o+QLM4+mAwUK/XM4QIZUCVttvt\njNbhOUF6JFtoOCqZRCKhSqWicrlsfvLQCPZQBaw58zCfzy0xMX8873A4NJ+lQuG58AXmh2BCtfnQ\nlsulgiBQoVBQOp22RAa6lN5VCrVaTblczvYT8wW4APkRuNlHfP9sNtN0Oo1cl/F4bM/5ww8/6O7u\nTpvNRpVKRZlMxgIoFRdUZ6/XU6fT0WKxULPZVKvV0vHxsQXv/X6vcrlsiRqLGgcUA0GcuSPRBUGg\nfr+vUqlkPgFNwlgBZv1+3+gNaBQCZrVaValUMhrooU2nUyWTSVWrVVWrVaOJZrOZVUjSO2BEhdpo\nNFStVrXZbDQajbRareznK5WKsRAk9GKxqGKxaPTcx9onReCgAQIcyJuStlgsKpvNGp0wGo0OMiOL\nvd1ubQMWi0Xj1NhQbLqoDLpYLCTdl6gkFBDqdDrVarWywFSpVHRxcWFcH9kTh/WomGTCM4AYw8yj\nheVyqfF4rHa7bT/DBgPx8f2UcARM+P5EImGJg6TE3H+oEsD56S3w+1OplIrFojk/Y4AKqNfrqtfr\nhmL4IhH6tWeTen7voVGJBEGg6XSqdrtt6NcnAygnz+3CK5JoSey+quP55/O5bfAwgyLr9XqG2HO5\nnK35YrHQbDazsrdQKFhgGo1Glkh8vySRSFgQhbYDCEQFcCgJ1phy3ycqElK9XlcikdBgMLAAO5/P\nVS6XFQSBZrOZ8a3ZbNZAA2DF9zUe2mazsX4KlWipVFK5XDaun3U6Pj62BJlOpzUcDiVJFxcXajab\nNk4/B6VSSbvdTv1+X5KMsnxoJC5PwyYSCVWrVdVqNVtb6AeQMfRIKpUyZAtQKBQKhnj9PGez2UjA\ns9lstNvtVK1Wlc1mNZlM1Ol0NBgMlEgktN1u1e12tdlslMvl1Gq1VK/XdXp6qtFoZAEcH4U6ZE2H\nw+EBcPwQnfPQPlkAh4Pi755CoWlVKBQsgEvvEPTV1ZUhFxwe/pHNQ2MCBCzdl6hhtlwubcODPilt\nPMJk43l+t9frGc8I173ZbAxV5fN5qwa2263xb2G23W4t2IEKcDSCJ5uAUtH3EQj0vvFaq9UMca7X\nawt+8P5hxuYlCUEJEKxJiCQ61u38/Fynp6f2nOPx2ILiarVSNpu1BLndbg2tgOge2mQyMVQ8HA41\nHA5tE0j33CPri0/VajU1Gg21223b5PCf4/HY1ptgTiKJ8g9+v3TfR4C3ZC5Zh1QqpWq1qouLCwVB\noLdv3xpCBz0TIPE5mlpw6VEIazgc2hh98mNtfA8I3+12u+p0OsbHP3r0yBAv60IiymQyFozYB1H+\nQXXLM9Ogg36pVCr2uzabjdFx9JTwy9VqpfF4bPvZV3mg9yg/JejhB3Dgx8fHury8NEqI9RoMBhoM\nBua/k8nEKm6ovl6vZxQKPSBf4YQZlR/xYTgcajwea7vdajAYqN/vazweq1qtajKZqFar2X6cz+e6\nvb3V3d2dAbDb21vrKVEl7fd7FQoFq+g+1j5ZACfIwvlKMtRXLpctm4GScrmcVquV+v2+UqmULi4u\nbLFByDwsKAVuzDdKwozJopyksQFFMp1ObSxMPBuZ4MzvyeVyuru7M+fy9AulflQlQPOU8aNiIJgS\nJKT7ZhyZGYSOumM6napYLNpn8f2MR1LkfOx2O83nc83nc202G1POMBZKdxIDgZyqyCdOSTY/GEnW\nK3vCzPsHa0QC9JQbz8dGrFQqOj09lST1+311u12rnkg8BDxPpURVJCBiQIJvfhJgms2mdrud8vm8\nWq2WKpWKJW2qORp0zB20Aeu0WCy0WCwiS/Vut2vfgy95bh9fpCEL8iOA87n0a0qlkqR7/py5oPcS\nZZvNRo1Gw0p9/AX0HASBLi8vTU3G3mu1WqpWq+p2u4baUWsBTiqViprNps0PSSDMQL31et3+jX2M\nsm29XqtWqxlVsdvtrL+FSqdUKllzs9fraT6fW5ICxHglT5ifsm9JiCQwnoPGL/Pe6/V0d3dnsQuh\ngKdfoG+JlcPh0Kq9j7VPFsApm3zXGy6Lck+SdfAnk4m63e5BI282m5k8Dn4J2Q5I1POdUSUyqhUQ\nMEoSxumbLQRCAs9isVAQBCqVSkZVeDmfl/lBf0QtAAiS8VLGQQ+waRkTXOxDREqZn0qlNB6PLdjA\nz5PUogIWaN43Y9nUD5vHJDkCKgqTYrEo6R2KJgFIOqAxQK8fSiQEKGgXPptki78QmOhd1Ot1a27i\nE15y5yV0BPUPzQf0DGoj0CJzXSgUDFWTKFinZrNpCBFgQLAloTInBPIoP0WeuN1u1el0NJlMzA/o\nKRSLRVPUsDfW67UGg4EymYzNDxUAdBCJgWeOQp13d3fWECWwUUHTtMPf2J/I5kj0CAUqlYrOz891\ndHRkAAblxWw2M3ohzE5PT3VycmJIfrlcqt/vKwgCa/xDqdCjYZyFQkGdTkd3d3eGbjebjSaTyXu0\ni1cthRkBH8SNgCCVStnak+BqtZrm87k6nY5evHihVCqlx48fq1qtqlwua7VaqVgsqt/va7PZ2ByO\nx2NLJFGAJ8w+WQAnu6F4IFiAGD3Cmc/nGo/HVrKBunyjgUlkM8HreT4xqkk1HA5NN14ul42vA42A\nkiRZAwh5ERxirVYz6RQIwm9WnvVDHCcO51UiXtIFp+4bkXwm5S+deJqeJJ2HhiQwzNi8UDGr1cq6\n5nwewZpE4hU3IAdoMBAuz+Q5eN98DRsH6AtqAuUHjUH8qFAoqNlsHihWkJMxVhKnr4qghz7E+RK8\nqXgI4KzvarVSJpMxBEXVQGO1XC4fzLVv3jM3zAO8eJh5VAqwAeVSMZHgfUOfwEmAozdAJYSv88W/\nRxmJhPUB8PBFpRQEgSaTifr9vhaLhaFQkiBAi6RM89VXRiiNwuz8/FyNRkOSrMdFNcD+ZN+1Wi3b\nT5PJxOZjs9no9vbWKhh+zqtYmJ9yuRw6Dt8o55yFV6VB5WSzWePGR6ORXr58afMhySo41gwfJ2Hz\nPR+icx7aJwvgoAyQHt105HZMIM7ps7p0H+wIMDRz0um08bulUskQAg3GMIM7pMTl4ArIkg3nETfO\njnIlnU6rXq9b44LglE6nbQE8eg8zggE0DpIkAiYLK8k4d3jwZDJp3C4cGx15NjLow6OKMCNg8cxo\nU2muSDpAtVRMrFk+n7fPrtVqB3p0r8MFgUZxzyBCEjNjZzOQENLptNEWcLDSu6BP8PQ6dB+8vYwr\nKnD6QIn6BKQJNVYqlZRKpYyu8ZJJtOn4EeMhULDurCNV40MjacAbE2xownl/ZK18UxcwNJ1ODdkh\n22W+CGQol8Jsv9+r2+0eIHCCHIG33+8rmUxaIxKESqPTJ3Aan/gD4MNLRcMM7TWInnUH7FFRNxoN\ntVotS3IkFnyhUqloMpnYvz1M1IlEwprRYdZoNCyZI2dlnX21fXd3p2+++cao1dlspsVioU6nY1z9\n559/bp9L1Us8JE59KLk+tE8WwG9vb01LDAL3pSVBjxIVdOXLaHSWLDALD3ImyFMaRzWHPOIksFFS\nM1kcbgH54pCgoeVyaaoZEDbdbLhQNkoU0qMBhaKGkt0/kySNRqMDrt+jdlD38fGxyaw8h+edIyqA\nk6QkWfkGsvKSSSoBkqdPFslk0ho5IDQ2JZuVRBYl70QtgI4YmgR6hHkEzbFG/E6UAzSS/bhJjATm\nD6lhOp2OOp2OVT8kQ18NopJaLpemZGBd8EGUERzq8YdFfAKJWhcSgKfp/DkJEiiUlaeKPEU0Ho/N\nf9kTJLVms2nAJYpSWiwW6na7evLkie1B5IM0mj31eXR0ZEGrXC6r1+tZFdBqtbRYLKxxx1io+j7U\nTKVPhm+0Wi1TJOGPmUxGZ2dn1rT0SJY5aTab1gCF9mKdqfb2+31k/Pjyyy+VyWTMr5G60mimcvzD\nH/6gdrut09NTS740f1GzbTYbnZ2dGS00GAxMuUMQj0qsYfbJAjgbgknlIA6bir+zEdk8/iQjjlYq\nlazhJungoAabiIcPM1AJAYnP8EiJgOEPxEj3aJXmQq1WUyqVsmaTP/bv1R1hBidM5QHCw6FJaEgS\nPdXCiUef3NA9S3rvsJRXcISNA56Q4EBVExY4SJwEJNAyyhH0y6AOAh3zEVUSZjKZAyULG5USH9+g\nSoPjRZ9NtQPFRKOTZMRaeLQXZpTd9XrdyuSH8rv5fG6IinVEs7terw0pkngx5Hj4DKgwyk8p132l\nxdwSsHK5nMrlsvl7Mpm0ahSu3mu50aB7NA4lF2bZbFbtdltv3rw5+HnQOweE2FdffPGFNXFbrZb6\n/b4lweFwaPuTxit+WS6XraoOM+gHEjW0q288kqw7nc7B8zBu5u3s7MySFp/Ps4CQo7jnn/70pxqN\nRoa2fd8BGeBoNNLV1ZWq1apWq5Wurq6Mp0+n07q8vFS/3z9IuIACJI70l6J6RmH2yQJ4v9/X0dGR\nlUhsZgICaJwAViwW7f8TSFn4QqGgk5MTa6xR/hC0POURZu122/TdBF4mkECVSCRsA3Gqzjfm4M8Y\nN4Ge/6Y0pCoIM36OzA2HR+kP957P57VYLA4OazAG0C/UBoHcN2UJZFFNma+++spKNeaP8Un3aI5N\nihb6+vpamUxGJycnRjH4asr3DyTZ3EY5ZL1et8BIAvfH0X2gQ2cLxeIrpul0qkKhYFw1QXK329lm\n2+12H5TvUc1kMhkbP78HtASNhJrh5uZG3377rX1PpVLR2dmZnj9/rn6/b4dagiCwZh20QJgSBbDj\nqQ+SG5RXOp02VRJBOQgClctla/pXq1W9efPm4HShD34AoCiq7+LiwmSirDN7EYUOwfHJkydqNpsH\na3p6emoSYmSqNPmGw6FVoV5uGWYPARrPSbLabrdG4XU6He12O9XrdTUaDWsUAr5KpZKazaZx0NAh\nt7e31t+I4sB/8pOf6N///d9tHrmHBZ/rdruaTqem+vnjH/+o58+f682bNwZGvvrqK2WzWY3HY93c\n3BxQnkgZoYKj5iPMPikH7mkTgq1HAgTAfD5vl8nAsXLijzJnv99bIxIEyGKwwWgWPDQOvrRaLY1G\nI+sce4QHsvInHHFSeMbZbKZer2fO7zW1IHuyepiBdClHvaoFVE7godzzh2ZyuZxVBAQzKgnpPgCD\n3KOQzcXFhV6/fm1OQ+AFXTOnHA+eTqcaj8dW6pFEK5WKJS8QLtUNPG8ikTjY4N6KxaKpLgiuXPQF\n4oNzJDjB1zNf+NDV1ZWNxzeY/Zhub29Dx7FerzWZTCw4ESxBcaw3ldo333yj0Wiku7s7vXjxwjTN\n/iRtvV63NYa39/xrmGWzWd3e3urJkycHvpDNZu0yrVqtZv7YarWMn6cvwX0fHMXHH0nooFkkqWEG\nuhyNRtbnIBHw+zabjT777DNJOqgOjo+PjZLgZ1CpoKsnwMJrcwDtodH49YkGwQE/y4GYWq2mzz//\nXH/+53+ui4sLLZdLffvtt3r79u3BHTjMAVWipy+9XNHb+fm5SqXSwWljT9tSSaM7f/nypTqdjl1J\nMJ1OdX19rSdPnlily76BEqJSRPH2sfZJT2JCexBYPIpkU0FrwJ+RlSiJPbL08jL+W5IFjKjMxQ1u\n/n6Ih6iVZpX/Av2jKyVocMESvwsnIIBHHVxh0/ggQSMPZ/fNMYIB6NDfLYMOlfn1fK1XjIRZtVrV\n2dmZbm9vrRKhUUng9AcoWAfuBeHAxGbz7ng1lYv/PK/siEI2NAbhbPk51BRQOHx5qgxnp6FIsEFK\nCO2AXGyxWOj169eR46DaIHmxFlBHVDO9Xk9v377VbDbT3d2dfS6USrvdtkTPAQ9PyUE1hFkul9N4\nPFa321Wr1TJpG2iPvYL/4n/M2WQysXXxSBlwAfVGFRdF5Zyenmo6ndppX1A4gR0UTNMQcAFwqVQq\nlswIlIgAfEVJcEZp8tDOzs7sd7IW9GzYn+ylZrOpv/iLvzCwsFqt9Gd/9mc6Pj7Wixcv9OrVK0nv\n+ktUvfwOqv6oy84ePXqkk5MTlUolo+ugVFGeVatV7XY7vXnzRpL05MkTPXnyRMlk0g74pFIpnZyc\nWM8JhRS9LaiUD2n0H9onC+CbzcZOoPlACapNJpOmOkkmk2o0Gjo6OjJhPotL0ADRSzpA8qAUmm9h\nBgqjfEPlgDP5AyCU9VALVAM0HuGoQeDeEaAaopANzkLJD3L3WuUguD8GTVBlzCwuPKdvokEjsSk/\nFMA5VACVQ8IiOHAsHL55u91a4qFvgKbf9wdIOlAKUGNRXLwPsAQ639jmOQAD0Dqe/oJDRxnBmPEz\nNka73dazZ88i/QPqx99y+fDAC+tH4uFSpmq1andfLBYL48bp9fi+CHMUtS6oT7z6w292aJ3ZbKZ2\nu631em1NYAI4clnWgTMQniPOZDKRASubfXdB1nfffafHjx+bYoUKgoBOVdLv9/Xdd98plUrZ8X8a\npgR6EKen+EDjX331Veg4/AVx8/lctVrtgLqg8ubKgJubGxUKBVWrVaVSKf3www+6uroyWnK5XFqi\noTrjy/tcmD19+lQ3NzcHB9t8hZtKpexyr/Pzc5trrj3IZDKq1WoaDAYGevB/6f1rmj/WPmkAh1dm\nwQlQ/B20yUYbj8fG0202G93d3R2oJfxDEjBBz9VqNbI0pYnUbDat8ehPUZFQ4OOpDNh40DkEJzSc\nODQcNnr2Xq8XOg5QE0HFB3CPoBkDqLff72symZhKAyTGIR/f1MFofIZZuVxWs9k0LS/PAsoF1UBn\n0SiCIy6Xy6ZHJ6AQXPk7G53TeWEGSkPRIt0f86eRC93m77/Z7XYmy/InbKG5kFuS1MbjsV68eGGS\nt7B1AdX5BioBnAQHlYcy4vz8XK1Wy+6hWSwW1uxig+LnoFT47TC7uLjQd999ZydhCcze5/EVgs5D\nyo4196jbrwlUJIkhzJrNpjKZjO7u7iyxSvdnFbh/hyZkp9PR27dvLVGA0pvNpt0nz/rSXwEFc4Iz\nzEqlkq2Zl9lBxfgDY6PRSG/fvtV4PDauu9frqdfrqVqtinvWvRYcWq7b7VqfLsyoeHwfgPEDpIbD\noQqFgp4+faqzszOrHqmCFouFXr16ZSeoAYPEQ+ldnCqXy5F+Gmaf9CQmXVucFkToFSAEjOVyaY0d\neCOCKiUxzsimIABw0COKS+Po7e3trc7Pz61k8lI1/pvmGdI0PtPTNaA7n0QkGVKIUjt4ORblIJvD\nl+7ZbNYoHRC1pwwIqJIMDTxUe9D4jRoHpRwb1OvuCT78PhAgjV0SH2W7L81ZH9/Yi2qm+oNPNPtY\nW+aDRMV4QH7+INNDOSnNZpCMR6thVq1WreHF7yaAz2Yzq+Bo2ELxIU0jCOAPXjHlm+v0NqIoJWgF\n5gOUSZBAwQFtRsCU7jXeJAnAAEnZV5wg+Kh1oeqcTqd2MAVp6na7tVODfCb3cUOx8CydTscqRoK9\nV2tBBR4fH4eO4/PPP9ef/vQnawATwKERCao8d6lUMjDAlRzNZtManfQJqIKhgEjSUVTOo0ePrIHv\nnwHdPeuNbHq5XKper9t4FouFXrx4oZubG/N3YgF7DB/OZrO6u7sLHUeYfbIAXiqVDvTAHkVI9w7H\nZveHUWgeEbA4Vei5RB+w0D5HcYtwtG/evNGjR49MZcCY2LQgJhxFulcisAk8f+gbdQTkD8mR4GdB\n4J5+8Q1IKgHuSvFKGUp9r5d/KC+kifajH/0odBzJZFI//vGPlUwm9f3335sTe408G9LfkQ6qJUGh\n9SXJsqE8xwnqCTOQNeU5SIfE7vls37sAebEpfUBnDvzZg8FgoB9++CGyIqHx5hMzARitLrfyoX4i\nySClwzc910zvgjVMp9O6uLiITKwcPwcwMH7+TtDwggCSFj8D6vM/74ESyR5aI8x4C0+xWFS73dZg\nMDCtN1UG6wNdQz+H5ElfYjKZqNVqWWPaJ+J8Pq+jo6MPUga+ImNv+VOyJP1Go2GodjQa2eFA7mBJ\npVJGEUITUjFDQz1+/Dh0DLvdTsfHx+p0OrbGXoJLslwul6YR7/V65m/0R+gBctWFP7SIoGG73ert\n27eR8/HQPmkAl2RNFRxP0gGy5ZBGsVi0bjolHkED9AACgoPGgVi8KBWK1zJfX1+bqB5kKemgoQjS\nAfVgoF9JwUDGAAAPA0lEQVQ+i8QDn8Zpr6iN4bXRlGwEcjhkEhr/TgIkeBJY5/O5KXHYDATNSqWi\narUaWaovFgudnp7qs88+M94UR0QiBarg8x82NZknScYfoophbVELRFUkJFJoKg7AkCRBs1694ZVI\n/CxUgj+tS8JbLBZqt9vq9XqRSK9Wq6ndbltA5lmZf4KSV0hJsjUCFRNUMU+VZTIZFYtFXVxcRM5H\noVDQ2dnZAe9LYiQIgxY5iesbeSQ8r4jxKiovl6WyCDPQcqPR0PPnz9XtdnV+fm6HzuiXgBipVkmo\noGT6Vsypp5NQt6BACzMoFtRXgBfGvdlsLEZAsRBcqXS9EgiqjINVq9XKklG1Wo30Dyg5qFQqAuJA\nqVSyqoRqnwvGvELHxzTAp1c8hfUq/jv7ZAEcxQQT5XlaHiqVSllHm83uD0Gw0D6A0GCiBPNvJ4ni\n9GhccLUjF8mAlEC7OICkAx065TwIk1KRTQE/ygVdUUiPn2MDev7cOyDP7REoKB0nYZNLMqeglORG\ntKhLtVgXXuZQLBbNqbnbAXRAgGKjkihIRHCKBAovGYWrjyrVC4WC6d0Zr0+YNFMTicPTc75Ryxh8\nKcrnowN/8eKF5vN5pJwRXpfrXD2Xi78dHx8f9ANQZaDZx8dRpPjzAdK76qtcLluyCDNAgU+enidF\n8+617SBrEryvAjw6J7kEQWCNxijOF3TdbDZVLpfVbrd1eXlpPSYqQig01pp/x6+l+z1M0GePsner\n1WqkaouqBekgB6kI4JwKJbl6X/B7kmY9754l2ZHITk5O9JOf/CSS2kIBVqvVdH19bfuc5242m/rx\nj3+sTqejXq9n1JunWR8yDT6RAeZYGx8D/zv7pNfJYiBMMqM/4AGCYrOCnvg50BWBD8RBYKhUKioW\ni6rVapHqj/Pzc2UyGXW7XdPvwrHj/KBcAoTnk/1hGRzQl4JeAw7SCjNoFj7Dl8g4FM00/o0DJHxe\nNps1SSOJDEcAAYFIog6udDodffbZZ3a4oV6v22bkOXxTlYDom6z+Lg6OzzMPJBLmMerCfm6lA6FQ\nTUiHyQ6EixTr4UsoSIgENPxkvV6r0+no6upK6XQ6sllGc5jETjOW38d9GyQ1PkM6fDmGV61QJUIT\noeUGEYYZlSd9COg2KjrUQEjZCAa+gQ3wYP95kEKQ4yBUVNOfF2BUq1Wdnp7q9vbWThkGQWB/gjC9\nP0gytEmwBOwwHtabpPz8+fPQcRAzaPrd3t5aBURzmt4F8+ZBFY1KLubyPTTQeBAEOjs705dffhkJ\nAJFTlkol1Wo1XV1daTAYWDL3vSniGX0S1hHahBjC+RZUayS3169f6+XLl6HjCLNPfpDHy8Ao5XgI\n+EVJFph4GwxlYRj1IskmBf6brnmYffXVVwdHcLl/AK6aUpAylbKVcYBqQBU4oVdFSPcvbIhqhgTB\nu2sw4eUoY0H4lIx8jqcrvI4YZyHASPc38IG+h8Nh5Ek7AnWtVtOPfvSjA+UEiRSU6WVfXtbmn5nv\nY5xIRCWZYiXM8vm80VsgfH93DsiEv4NG0RX7stwfEGNtZrOZOp2OttutHj9+rEePHoWOo1qtHrzB\nBj6SNUDyhg/gkwQOEKdPJD6Y5PN5nZycmDohaj787+XZoCAIvoyHz/R3vPh9QtDGZ5gnkq2/wuCh\nce9PLpfTl19+qZcvX+rt27eq1+um5vC3Yi6Xy/eAB/0PUCVrSUJmPNfX17q5uQkdB5JNxppOv3vZ\nSaPRsD4HqBbZJqiY9SJZ+GY//a79fm9VwG6307fffhs6jj/84Q86Pj42Kmu1WqnT6dh1B5yN4JyJ\nv7iPYM4JWa5TYP+STFOplG5ubvSf//mf/3dQKJx+9FcmPpS2+cDuF4PA5MsjX1bA+cJJ0Q2OUhl8\n+eWXFugp2ZlskJ5H15KM7yQoeP23NxAkpbY/gv/QQJEcJUYnz0k30DufBYr0iY+sD+KAdiKI0Ueg\nsRdmJIN6va6joyPTUPPmIdA0Vxb49zAyH2xc6bBMlu5VNJ63DDOadl7pg3zwYcP6YUJ/KK2jSvCB\nFklnrVbT5eVlJMd5enpqQe/u7s4CEOUvARO/8NQW6wq1QaONTZhKvbvk6unTp7q8vLTbLKPMV6al\nUsnWlWqEvYCPhD27N5A3nL4/9fqhSpEg/vTpU/3sZz/T73//e3W7Xdt3PlGzn7xsjkSED3i1Vjb7\n7rVk9CY6nU7oOAhwkuxdkldXV3YimOeDTvI9G0lWIXjJIX4KLcPVt8+ePbNDOA/t5uZGtVpN/X5f\nz58/17NnzzQej1Uul+2KAMZAYgJZ+xiXTqftxTXQMAT7zWajb775RplMRhcXF/rd734X6SMH6/tR\n3/W/MMoX35WGT/Z0BJNJebjdbg/eakLTymtZQZncEV2v1zUcDiO72ZVKRUdHR4aQuJENSoPMyub0\n6MFTJaAeEI9XLJB0er1eJOc7Ho+tcUICAWUQwAnMvnFGScyJUMoz37yEDqpUKjo5OVEQBLq+vg4d\nB/riIAisrH/16pV++OEH6+DDJaOLphpgTkgqUEg+kCKFIwlFBc4gCHR0dGSfR1W13W51dXVlFJNX\nAkn3J00xv158/2Zz/4YWXvUVhWx4uw4JCz/inmvmVpL5JAmFoMHJRa46BXFmMu8uY/vss890fHxs\n6xpm/nevVit77+Z6vTZu1PP/BEmvymEvUZGxHplMxio+brKM8lOQOsn15ORE9Xpd0+nUmsEoT6gI\nCUz4iW/G+aY+gXw0Gunbb79VrVaL5HwJ3svlUu122wQCnA8hgYKu6T2QNBANMD6ejTHxHKvV6qAi\nf2ipVMp07qvVSre3t1b5zGYz69/xtiKf0KnIvM8ANKk+8/m8XdT285//PLJyDh3bR3/n/9BAcaBj\nUColKFmaksRfeIU0DQUDmZ1mF2ia4Esg87y1NzaNl2RRFczncwsaXmrlN4Pnmdm0BG+kSuVyWScn\nJ0omk5GIgiDMNaxwnByKAFXu93vbICQUAkG1Wj043eh5UJ6BE2tRgeL4+NhKUO5+qVQqajQa1khE\nlcLnEER9MOf5fQOYBg7vHvTa4IdWLBbtrmXpvmQeDAYH90fzp0/60n0Qk3SArLweH1qJyiHMoLGC\nILC3yPi3n3vpqz/YA5LygRs1Ed+fz+f19OlTO+BRr9f129/+NnI+giBQrVZTq9XS0dGRlsulbm5u\nDFlCO/lj2F7CSdUmydaOQ2jQMY8fP9bR0ZG63W7oOPDBarVqL+XlMIp/IQLBirUAnEmyPQawwRc4\nqfrq1St9/fXXSqXevT4xzFCojMdjff/99+p0OioUCur3+8pms6rX6weyykQiYTSYV26xtqwb88Pt\niTc3N/rTn/4UmdD6/b5evnxpr7LDl7gmgBcRU52kUimTinp57mLx7nV6vgldrVbVarU0GAx0dnam\nX/7yl5GJJMw+OQdO1icDkYE9l0qQBDmAsECknm/0x9rT6bQhxtPT08imzN3dnaksQNloWqfT6YHC\nAdUMaMKrQ3xz1TeLkFyBjqO4RYIqTUM2uHSPIiQZDULwBj2SLHhXIcEEbpEGGZv66OgodBzValXz\n+VxXV1dqNBoHh4gISl4zSxMHCsy/EYnPZz6CINBgMFC73baAG1UZPX36VI8fP9ZyubQ70H0n3isp\npHsk59UMoBoSLMF7MBgol8vp8vLSfjbqhGyv19N+/+6OG65x8AEcBEx/AoTHl39NF34D4KjVanZL\n4RdffKHb29vIS7V4ycfl5aUuLy+VyWTsbhXeBO+f2dNq+C6Umz84588UNBoNPX36VMfHxx+8g4T7\nTLixjybwfr9Xp9NRIvHuDTscbvHrIN3flc+hJPwJDv7777/X69evDYCFGQGfuX39+rXOzs7Mx4gD\nnu+nauTZ8UmkiAAUgjXV02AwiJTd/vGPf9RsNrOXetfr9QNNOO+3JFHR1KY3RxIjXrA+HPunQbpY\nLOydoh9rnzyAcw8BwQr07O9T8HJB3zVn88P7cucEfGwQBBqPx8pkMrq9vY28TYzTUtxHQEkGh+wb\nmKA9FoMvxoqTgugp84MgMCnjh07aceR6Op3aW2agSvb7/cELaZkjvzkpiwnYDxUGvP3koSb5oaFb\nB4m/efPGmi6SrHE4GAwOgiXNZppclP1+zdhwyWTSpF5hRvlPw5O7algvNiGJgP4Im9VvDF8dkaxp\niHMgKUqfjyzt/PxciUTCqAKSIXPCxVXQJr7/4dU7vvm72WzU7XatRD4+Po48cNZqtdRsNnVxcaEv\nvvhCvV7PEhE+5wOBP8zlVR5+DF6VQYXBlatRiZXf+ezZMz179kzdblfp9P1lVFy/y50+vvnq0bC/\noRK/LZfLNr/j8Vi5XC7yjhpOyHa7Xb148ULX19cGaODh8aHd7vDN8X7fQr+ClDmlXCwW7QZHkHWY\nQd1Qaf/0pz+18dF0B0gS53zfjqoeygS/pU+DrwOkeDHIx9gnbWJy8xcZGo6WjQWHJskeGB0uqALa\nhCPIZF3Q/PX1tR08iLpzA0fLZDLq9/u2WGx2uLbtdmvddQIGm8Wjv4dfBACcPEo2Rxmbyby7vezk\n5MRKYlQWjIX5IJl4hAOi4TAT3zMej3V9fa3RaKR6vR4ZOGezmUqlkqEa1BpIxEhU9CG8ygR6xXPS\nVFcEFE6CDgYDO5kWZv5agHa7bW8ugbv0lMVD/le6V+l4uZqXZnK9Ak3jKF18o9FQo9HQq1evNBgM\njEsluflGN5uQioMg4uWlvnpk4/7Xf/2XKpWKnjx58kGVUqVSUavVsnnx/DZVqlcIASx8VeI5VJ9w\nAEP4TZT+mteAff311/Z6tkajoePjY1sfki1CAt+4p5L26+yljQAzElPUujCXvg+DL1Bpsef97/MJ\nlcTiaSXojnK5fNBwj0rwACqC+OXlpbrdrjqdjl2twJioglkv6F3fbEVyud/vjSeX7u8GiqJyYost\ntthiiy222GKLLbbYYosttthiiy222GKLLbbYYosttthiiy222GKLLbbYYosttthiiy222GKLLbbY\nYosttthiiy222GKLLbbYYosttthiiy222GKLLbbYYosttthiiy222GKLLbbYYosttthiiy222GKL\nLbbYYosttthiiy222GKLLbbYYosttthiiy222GKLLbb/j+3/AMGOWThP4SsrAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa7a8148110>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "###保存Eigenfaces的图像\n",
    "n_row = 10  ###显示多少张图像\n",
    "fig = pylab.figure()\n",
    "for i in range(n_row):\n",
    "    p = Eigen_faces[:,i].reshape((32,32)).T\n",
    "    fig.add_subplot(1, n_row, i+1)\n",
    "    pylab.imshow(p, cmap=plt.cm.gray)\n",
    "    pylab.axis('off') \n",
    "pylab.show()   \n",
    "fig.savefig('Eigenfaces.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# FisherFace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 209,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_SB(data, label):  ### Get the between-class scatter matrix\n",
    "    \"\"\"\n",
    "    data : M x N matrix (M dimensions, N samples)\n",
    "    label: The label\n",
    "    \"\"\"\n",
    "    M, N = data.shape\n",
    "    mu = np.mean(data, keepdims=True, axis=1)  ### 全部样本数据的均值\n",
    "    C = label.max() + 1  ### label的个数\n",
    "    S_B = np.zeros((M, M)) ### 类间散度矩阵\n",
    "    for i in xrange(C):\n",
    "        i_index_list = [ind for ind, x in enumerate(label) if x == i]  ### 第i类label对应的下标\n",
    "        data_i = data[:,i_index_list]  ### 第i类样本数据\n",
    "        mu_i = np.mean(data_i, keepdims=True, axis=1)  ### 第i类样本数据的均值\n",
    "\n",
    "        N_i= len(i_index_list)  ### 第i类样本数据的数量\n",
    "        Sb_i = N_i * np.matmul((mu_i-mu), (mu_i-mu).T)\n",
    "        S_B = S_B + Sb_i\n",
    "    return S_B\n",
    "\n",
    "def get_SW(data, label): ### Get the within-class scatter matrix\n",
    "    \"\"\"\n",
    "    data : M x N matrix (M dimensions, N samples)\n",
    "    label: The label\n",
    "    \"\"\"\n",
    "    M, N = data.shape\n",
    "    C = label.max() + 1  ### 类的个数\n",
    "    S_W = np.zeros((M, M)) ### 类内散度矩阵\n",
    "    for i in xrange(C):\n",
    "        i_index_list = [ind for ind, x in enumerate(label) if x == i]  ### 第i类label对应的下标\n",
    "        data_i = data[:,i_index_list]  ### 第i类样本数据\n",
    "        mu_i = np.mean(data_i, keepdims=True, axis=1)  ### 第i类样本数据的均值\n",
    "        Sw_i = np.matmul((data_i - mu_i), (data_i - mu_i).T)\n",
    "        S_W = S_W + Sw_i\n",
    "    return S_W\n",
    "\n",
    "def LDA(data, label, K=14):\n",
    "    S_B = get_SB(data, label)\n",
    "    S_W = get_SW(data, label)\n",
    "    N = len(label) ### 样本数量\n",
    "    C = label.max() + 1  ### 类的个数\n",
    "    \n",
    "    #### 先执行PCA将原始数据的维度降到(N-C)维，用以解决S_W的奇异性问题####\n",
    "    S_T = S_W + S_B\n",
    "    eigenValues, eigenVectors = LA.eig(S_T)\n",
    "    eigValInd = np.argsort(eigenValues)\n",
    "    eigValInd = eigValInd[::-1][0:(N-C)]\n",
    "    W_pca = eigenVectors[:,eigValInd]\n",
    "    \n",
    "    ##### 再执行FLD过程#######\n",
    "    S_B_bar = np.matmul(np.matmul(W_pca.T, S_B), W_pca)\n",
    "    S_W_bar = np.matmul(np.matmul(W_pca.T, S_W), W_pca)\n",
    "    eigenValues, eigenVectors = LA.eig(np.matmul(LA.inv(S_W_bar), S_B_bar))\n",
    "\n",
    "    eigValInd = np.argsort(eigenValues)\n",
    "    eigValInd = eigValInd[::-1][0:K]\n",
    "    W_fld = eigenVectors[:,eigValInd]\n",
    "    \n",
    "    #### 得到最优 W ###\n",
    "    W_opt = np.matmul(W_pca, W_fld)\n",
    "    return W_opt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 201,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "train_face, test_face, train_label, test_label = load_Yale_data(randomsplit=5)\n",
    "Fisherfaces = LDA(train_face, train_label, K=14).real"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 197,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1024, 14)"
      ]
     },
     "execution_count": 197,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Fisherfaces.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 210,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAA9CAYAAAC0qDI5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztvcmPnml5PXzeeZ7nsYa3qmyXXbaxu6EH6A4xoERswiIR\nYpEViyirCLFJpPwFEavsImWRAMkCIQJiSgSEDqHbbtx2u1yueXzneZ7Hb+HvXCm33zfi95O8+b7n\nbBDd7vJdz3M/131d5zrXuQEFChQoUKBAgQIFChQoUKBAgQIFChQoUKBAgQIFChQoUKBAgQIFChQo\nUKBAgQIFChQoUKBAgQIFChQoUKBAgQIFChQoUKDg/1NQvaof/B//8R8zv9+Pi4sLVKtVnJ2dYTQa\nIRaLwWw2I5PJ4PT0FKlUCrPZDOvr63A4HLh69So2NjYwnU4xHo8xHA7RarVQq9VgNBoxGo0wmUxg\nt9thNpuRTqeRyWTg9/tx584dvPXWWy/9Tnfv3p3dvn0bFosFRqMRs9kMxWIRx8fHmM1mcLvdcDgc\nSCQSiMViAIBisYhyuYxmswmVSgWVSoVWq4VgMAiTyYRIJIJIJILpdIp6vY5CoYB+v4/xeIzt7W38\n+Mc/fmkd3/3ud2dutxvf+9738PHHH2MymSCRSGBtbQ1erxdutxsAkE6nUSwWodfr4Xa7YbFYMBwO\nMZ1O0Ww2MR6P0e12YbPZMBwOUalUkM/n0e12odFocOXKFaytreHrX/86PB7PS+v4+te/Ptvd3YXd\nbsedO3fQarXQarXQbDahVqvh8XhgtVphNpuRy+XQ6/XQ7XYxGAywvLwMj8cDn8+HYDCIWq2GdruN\nyWSCyWQCjUYDs9mMUCiEer2O4+Nj2Gw2/M3f/M1L6/jGN74xc7lcqNfr0Gg0iEajiMfjAIB+vw+/\n349er4fT01M4HA7odDpotVqUSiXk83kYjUa0Wi2kUil4vV7Y7XY0Gg30ej1Z93g8hsPhgMvlwhtv\nvIG/+Iu/eGkdX/va12ZWqxXlchkGgwHdbhdqtRpWqxV6vR4ulwvD4RAAcPXqVajVatTrdTSbTWg0\nGthsNqjVarTbbQwGA3S7XYzHY+j1egyHQ/R6PZTLZYxGI/h8PnQ6HXz/+99/aR3f/va3Zx999BEm\nkwnK5TJyuRz8fj9isRiGwyHa7TYMBgPW19dhNpvRbrfR7/dhNBphNBoBAMPhEFqtFtVqFZPJBGaz\nGb1eD5PJBAaDAV6vFzqdDpPJBJ1OB9/85jdfWsdf/uVfziKRCHQ6HVwuF6xWKzqdDjKZDEqlkjxX\nfk9cA/enTqeDWq3GZDKB0WiE1WqFx+OB2+2GRqNBp9NBo9FAKpWCRqOBz+fD3/7t3760jnv37s0i\nkQhMJhP8fj+CwSAMBgP6/T76/T4GgwHOz88xm82wsrKCSCQCl8sFr9eLyWSCXq+H2WyGfr8v31M2\nm0Wr1ZL11et1jEYjGAwGdDod/NM//dNL6/i7v/u7mcFgwNnZGY6OjjAYDAAAiUQCN27cQCQSwdbW\nFgqFAn7xi1/g+PgYLpcLsVgMKpUK9XodBoMBN27cgMfjAQAcHh5ie3sbyWQS6XQaJpMJgUAAPp8P\n8Xgcf/3Xf/17xWbt7/OH/m9QKpVwcHCA8XiMarWKTqcDs9ksH7rFYoHT6YRGo0E4HJZfjJvP6/XC\nZDKhVCqh0+lArVaj1+uhVquh0WjAZDLBarWiVCqh2WzCYrHAYrHMXUsgEIDD4YDRaIRKpcJ0OkW/\n30c2m0Wz2YRer4fdbodKpYLT6YTP54PD4UC/38dkMkG73UY2m8VoNEK1WoXP50Ov14NWq4Ver8ds\nNoNer0ez2cRsNoPZbJ67DoPBgOPjYwwGA3i9XthsNmg0GrTbbcxmM5yenmIymUClUkGv12MwGKDT\n6cDtdstHZLFYUK/X4XQ60el0YLfb0e/34XK50Gw2UalU4PV6ZbPOQ6VSQalUwrVr1+B0OjEajWA0\nGhGPxzEej6HT6WA0GuHz+XDz5k0cHR2h0WhgPB5jaWkJZrMZg8EA29vbcqDo9XpYrVZotVoJNBqN\nRj6gebDZbNDpdIhEIjAYDAgEAlCpVCgWi3JoTyYTJJNJqNVqeL1euFwuqNVqOBwOtNttJBIJvP76\n6wCA6XSKXq+Hfr+PXq+HYrGITCaD8XiMZrOJarU6dx21Wg2VSgXBYBBarRY2mw1+v19+DgB4PB7o\ndDpYLBa4XC7YbDYkk0kMh0MJ9olEAvV6HYPBALlcDgDQbrehUqngcDiQSqVQLBbhcrnmrkOr1UKr\n1aLb7aLT6SAYDMLpdGI2m6Hb7UqgTKfTcLlc6PV6aLfb0Ov10Gg00Gg0smf9fr88D4vFAr1ej1qt\nhmKxCLvdDp1Oh2QyuXCf2mw2CfpMUAqFAgaDAcxmMzQajRwERqMRk8kEAGCxWNDv9zGdTmE2mzGb\nzaDRaGStOp1ODuLRaIR+v4+NjY2562AyYTAY5OButVoYjUYYj8eYTqdIp9OYTCbw+Xyo1+swm80Y\nDocYjUbo9XpIp9M4Pj7GyckJwuEwtFotDAYDnE4nHA4HrFYrer0e1Go1nE7n3HVsbm6i3W7jwYMH\nGA6HaDabsNls8Pl8soczmQzS6TQ6nQ6KxSJqtRo0Gg0sFgum0ylCoRC63S4MBgOMRiNcLpccOru7\nu8hkMnC73QiHw+h0OnPXMXfP/N5/8v8QjUYDmUwGZrMZJycn0Ov1cDqdEpyWlpZgs9lQKpUQCASw\nvLyMwWCA4XAIlUoFo9GIbreLXq8ngVOtVkOr1cLr9UoQ6ff78rIWBXCbzSZZq8FgQLPZxGAwwOrq\nqmSeLpcLrVYL7XYbdrsds9kMg8EAs9kMarUa0WgU/X4fZrNZMp5CoQCbzQan04npdIrJZCIbdh4O\nDw9hMpng9XolUzOZTLKucrmMhw8fQq1WIxAIyCbR6/UIBAIAgNPTUxQKBezv70um63Q6YbVa4fP5\nMBqNJBNiEPkkjEYjtFottra2ZP0mkwnT6RStVgulUgnZbBadTgeRSAR2ux2tVkueNZ8Pg1cmk0Gl\nUoHH48HS0hIMBgOGwyEGgwEGgwF0Ot3cdTSbTfj9flitVgwGAxQKBYxGI+h0OgQCAYxGI+TzeZhM\nJpTLZYzHY/j9fskK+c8ZxBio+PsEAgEsLS3h/PwcR0dHCw80l8uFZDL5QnDhzxuNRlCr1ZI9TqdT\ndDod6HQ6mEwmyTz57/x+P2azGTweD0ajkezRs7MzpFIpVKvVhYGiUqmg0+lgPB7D6/UiEAggGo1i\nOByCFUKhUMBkMkGxWITP54PRaIRGo4FarYZOp4PZbIbFYkG320W1WkW1WoXNZkMikYDVakWr1UK9\nXodOp1sYKGw2GxqNhhy+pVIJtVpNDhCn04loNAqdTodCoYByuSzf7HQ6lW+VAWk6ncLlcknFNBwO\nodfrYTKZ5J3PAzN8HsCsgv1+P3w+H5rNJkwmE9RqtVQkrVYLKpUK5XIZ/X4fnU4Ho9EILpcLKpUK\ns9kMFotF9iSrA36/81CtVjEcDtHtduFyuWCxWBCLxTAajZDL5dBsNtFsNrG/vy8xazKZSJzR6/VS\nAVxcXGA6nWI6nUKlUsn36PF45LkuOljn4ZUFcKvVitlshnw+j2w2i6WlJajVahgMBvh8Pjmhw+Gw\nBNjRaIRsNotSqSRZGMHgqdPp0O/3odfrsbKyArPZjO3tbSmb54Hlp8/nQ7fblQrAZDLh1q1bkoGT\nnhmPx9BoNJJRMiN1Op0YDAay4ZgtptNpCWzMYBeBQXA4HCISicDpdOL69eu4efMmjo+Pcfv2bWQy\nGezt7aHdbkvlUa/XpRSvVqvQarVwOp0IBAIIh8Po9/sS1BqNBobD4cKA5ff7EY/H4ff75QOx2WyY\nTqc4Pz+XLIfVCcvRcrmMRCIBlrU6nQ56vR7BYBD1eh2VSkWyDZVKhUePHsk/mwdWYcPhELPZDCaT\nCXq9Hp1OR2gJBtFf/vKXAIButysHVr/fx/HxMXZ2dmAymTAYDBCPx/GlL30JbrcbJpMJLpcLg8EA\n7XZ74Trsdjui0ahklmq1GoPBQJ6DxWJBq9WSoMGKglWYyWR6Yf+Vy2XMZjPYbDbEYjEsLy/DaDTi\n7Ozsfw1YzWZT6AmVSgW1Wo1+vy8VR7VaFVrJbDbD5/PJdzUejzGbzSSDb7Va8kzr9TouLi6wtraG\nUCiE4+NjtNttOByOuevgnuOaZrOZ7HW73S4HnlqtlgOnXq/D7/fD6XTC6/XKukejERqNBhqNBvx+\nP3Q6HabTKUajEUKhELRa7cLAGQgE5EDVarVQqVTyjVWrVXnHbrdb9vvJyYkcVO12GxqNBiqVCn6/\nHw6HA6PRCFqtFmazGQ6HQyoKu92Oer0+dx2sxK5cuQK73Q69Xo9qtYrt7W243W4YDAa43W74fD5J\n+Lh3YrGY0MGlUknii1arleTM6XTCaDQiEonIPvx98coCeCaTQSqVQigUwmc/+1mo1Wo5yXO5HPR6\nPfR6PQCgUCgglUqh1WqhUqmg3+8jFArB5XJhMplAq9UKT80SMZ/PQ6/XY319HWdnZ+h2u3j//ffn\nroUBjiehSqUSasRut8PhcLxwAjebTQnmwPMyGAAmkwkCgQA6nQ5qtRpKpRLq9Tq0Wi36/T6GwyEM\nBgPUavXcdezv78PtdmM0GiEajUrWvbu7C71eD4PBgFgsJhn8xcWFBB+uW6VSQavVSmBYWlqCz+fD\ncDiUiuDBgwcol8s4Ojqauw6Px4NgMCg/OxgMQq1Ww2g0Ynl5GaVSCf1+HzqdDs1mE4FAAEaj8QUO\ndTgcyqFCukuj0WAymSCXy0Gn06Hb7aLdbsNkMs1dB2mkVCqFyWQif67ZbGJ3d1cy4kqlAo1Gg7W1\nNcmuWZr6/X643W60222USiXJ2p1OJyqVimSfpHMW7Q+PxwOHwwGz2Yxut4tmswmtVovBYIDpdAqD\nwSAfn9FoRLvdRrVahd1uh9frhUajQa1WQz6fx7Nnz1Cv17G6ugqVSiXZeigUkux8Hra3t6Wa0Ov1\nqFQqSKVSGI/HMJlMKBaL8qyZHbManEwmGI1GyGQyyOVyKBQKcDgcWF9fR6fTQafTQblcxtLSknw/\nkUhk7jri8bjQAPzdY7EYEokEzGYz+v0+UqmU0Jk+n0+CYiAQwHQ6RaFQQKVSQbvdRq1WQ7fbRTAY\nxNLSktBCVqtVOPJF+9RqtQp9FwqFoNfr0e/30Wq14PF4MJvNJCj7fD6pmmezmVBwdrsda2tr0nvj\n4cS/2+12w2w2LzxY2Utg72cwGCCVSsFoNMrBwGCsUqmkygcAk8mEyWSC8/NzVCoVGI1GxGIxBAIB\n6Z+wag6FQkgkElheXp67jnl4ZQF8Op3CYrFgY2MDsVgMFxcXyGQyQtSHQiHhyAqFAprNJvL5PDqd\nDpxOJ3K5HBqNBkKhEFqtlpSi1WoVzWZTuChmvolEAuVyee5a+JH1ej00Gg35qJxOJ2w2mwQy8ll6\nvV5OUIfDgVarJXzfdDqFTqdDOByWn1+v16WB2Ov15GD6JAqFAlZXV6FWq4WTLBQKaLfb6HQ6uHbt\nmmzsUCgk2TT/fLfblSwklUrJ5u50OtLcYTlpt9sXBnCWma1WC+PxGLVaDfV6XZqw4/EYxWIRDodD\nMuBer4fBYIAPP/wQ1WoVlUoFer0eR0dH8Hq9AJ5zyfwIyK3zUJsHUgvD4RA6nQ6j0Ug+SB5eBoMB\nBoNBnv9wOEStVoPP55O9tLy8jMPDQ7z11ltCt41GI6lWdDodcrkc0un0wr06Go0QDAbRaDSk8jEY\nDFK+u91u2O12GAwGSThyuRy8Xq9UPqenp5jNZmg0GkLzlUol2Q9sOi46SEqlEnQ6Hfx+v1BIzD7b\n7TbC4TD8fj+63a5QS4PBAFevXkWr1cKDBw+kma7VaqHRaOQwZnZOSiEej0ti8kkwuWGywKpZr9cj\nGo2i3W6j2WxKckVKjYmN2WyGWq1GsVjEdDqFVqtFNBqFSqUSfno8HqPVauHw8FAa1/P2KanUUCgE\nh8OBbreL4XAIm82GdruNRqMhB73P54PBYJDGOr8hNteZeKhUKqFbQqEQvF4vqtUqxuPx3HWwuW8w\nGGCxWIRyDQQC0Gq1WFpaglarxfn5OUajEZxOJ27fvi09rouLCxEs9Ho9EStUKhWYTCYEg0EYjUZ5\nXsVice465uGVBfDNzU04HA4JAna7Xcr0i4sLtFotKa1CoRACgQAePHiA2WwmWWG/34fT6UQ+n0cm\nkxEqg8HGZrNJh5pZ6Dxw8zcaDfT7feF72Rxh55yKiMsfr8vlkq63Xq9HJpMBAAns5O3JawHPOcR5\nuH79OlZXV7G3twetVovhcIh+vy9ZeT6flyYZG7cajUYahCaTSXi/dDoth1i5XIZerxcKgFkZ1/pJ\nkKMkd3t4eChVj0qlkmyEdIHD4UAul4PNZoPL5ZLsigqaXC4nvQ1mr/1+H5VKBdlsFna7fe461Go1\nGo2GZCOk3PR6PRKJhPCosVgMP/7xj5FOp+Hz+USlwfLUZrNhbW0N6+vraLVaOD09RTqdlkZZq9XC\n0dERTk5O5q7j4uJCuNRSqQSNRiNNQuA5R67T6bC9vQ2fz4darSZZZTabhclkwurqqjQeXS4Xrl+/\nLg21Wq0me4oUxzxEIhGo1WpsbGwgkUjgyZMnEjSsVqsETCpIBoOBZMHkXS0WC87OzvCpT30KxWIR\nxWIRVqtVVCPc65PJBLVabe46Go2GCAlItTCLLJVKmE6nqFQqKBaLuHHjBjQaDcbjMS4uLlAqlRCP\nx9HpdGAymTAej2EwGKQJP5lMUKlUAACdTgfD4XDhwXpwcCCqsdlshl6vh2q1Ks3UcDiM8/NzmM1m\nXL16VYQJPLDYL/F6vUilUtJEZZ+KmW+9XpfnNA8GgwHtdhvr6+sYjUZS9dVqNVSrVTx9+hTlchk6\nnU4SDvZ1Dg8PAQDRaFSSovPzcxQKBWi1WkynU0QiEdhsNmm8L2q2z8MrC+DlchkWiwUajQb5fB6T\nyQSDwUBeHk/ser2OeDwuD9Vms+Hi4gJqtRo2m01eoEqlQqfTwXQ6xdLSEtxuN9xuN7xeLxwOB0wm\n08KT3OfziZyn1Wqh1+sJH95oNKDVauFwOF5ovF3uwofDYTx79gy1Wg27u7uSffOkNxqNogbR6/Ui\n6fokHA4HvF4vxuMxBoOBcISUm52cnMDlcsFgMMButyOfz6Pdbktzzmw2S2Pqj/7oj6DRaODxePD0\n6VMAEOrG7/eLOmAe2Dg+Pj6GVquVLjwzlp2dHSlvU6mUNOacTieCwSA++ugjGI1GJJNJfOYzn5GM\nS6PR4OTkRJozbL5RdvVJ9Ho9WK1WWK1WGAwGFItFjEYjrK6uCqfZ7/eRy+VgNBrh9XqhVqtRLpfR\n6/WwurqK2WyGSCSCt9566wXeNZ/PC9fIwLmIymm1WnC73S80QPk+qTjQ6/UiSyP1Q6miyWSCzWaT\n/ReLxZDL5RCLxTCbzaSaI6e/iGIrFovY3NyEy+USJQ+VUcw+Gczi8bgc5DxwqG74yU9+guXlZahU\nKqEHNBqN0DBmsxn5fH7he2EWa7VaodPphJdvNps4PDyUA6RareLWrVvC2ZIbBp5XV1SkaDQapFIp\nqcpYkUwmE8RisRd6XZdB/vyyRLPRaEhVQwkiq/XJZIJSqQSVSiWJg81mw/LyMn73u99J8jYYDNBs\nNmEwGGAymSSZXNS7slqtUKvVcLlcOD8/R7/fR6PRQKFQQKvVQrfbFcrq6tWrcDqd2NnZQavVwpMn\nTxCLxRCNRoUGGgwGCIVC0oejIGM4HOL8/Hxhb2IeXlkA393dxac+9SnE43H87ne/g9/vFw6NpSQb\nmNSGsnnidDqRyWQQi8UwmUxw+/ZtjMdjpFIplMtl2Gw2hMNhOBwOBAIB5HI5xOPxheS/2+0WPo6Z\nL3W5/X4fwWAQmUwG0+kUwWBQAmW9XofJZBK5UaVSkYYeSzYG13q9Lo0klkifxPHxMd59912Ew2Gc\nnZ3B4/FgbW1Nqg5mRNPpFMPhEPv7+9IkZKMxkUhgb28Pu7u7mE6nL+imbTYbPB4P8vm8qDPmgWWs\nw+FAr9fDcDiUzJJddB4YkUgEo9EIy8vLUr4fHBwgkUig2WzK70EaiR99s9kUvrbRaMxdx8XFBTY3\nN1/QuOt0OtGV9/t9FAoFPHv2TLJ1q9WK09NTuN1ujMdjbG1t4Ze//KU0+JaWlhCJRHB8fCzZ9NbW\nFjqdDk5PT+euY2VlBVqtVvTBs9lMVAJ6vR7T6RRqtRrhcFhkYKzUWBl1u114PB6RY+ZyOVFF8bkG\nAgHZL/Pw+PFjLC8vy+9ttVoxnU6l+VWtVhEKheB2u9FsNhEMBjEcDmE2m2XPcJ80Gg3h5XU6HWKx\nmGThVJgsqkjY6yHV02g0RFZKinE2m+Hw8BD1el2kwe12W/ZOIBBAMpmUZjsPRR6kg8EAbrcbk8lk\nYU+AChybzSbVA+lCnU4n37DRaBS9vcPhgEqlwsHBAQqFAsLhsDQ8j46OhFriOrhvzGazJB6fBClH\nJmij0Ui+H7VajXw+j36/Lwf7dDpFJpOB1WqFxWKBw+GAwWDA1atXMR6PpY8wGo1kLWxas+fy++KV\nqlB0Oh08Hg80Go1sSIvFIlwmu7WUM5nNZpydncFms0Gr1Up3OxQKYXV1FcFgULIxloB8WOvr6wt/\n8WQyKae9wWCQDcwg0Wg04PF4ZFigVquh1+tJM5GKlGq1KpvqcllM6oal1v+2Echd8kX7/X7Rt/Kl\nlstlHB4evqC4GA6HcLvd0snf3d1Fo9HA7u6u0CGUt7GEXCRXY2lLLp0fxe7uLkwmE5xOJ1wulyhw\n6vU6VlZWMJvNhP5iY/ng4ACvvfaaBLzJZAK1Wg2fzweTyYThcLiwJCS32Gw2RZ45Go2k6jAYDMhm\ns8KzslEHQIL+5uYmjo6O8MEHH+Djjz+Gz+dDv99HrVaTw49Vj9/vX7hXg8EgVCrVCyqUwWAg6iRW\nE4FAAIVCQbhltVoNlUqFWq0mgZVBnr0Daoy9Xi+KxeLC9xKJRJDL5RAOh0UBRIkadd7dbhdnZ2cA\nnqt4+PFT+thoNHDr1i35NliBsIdBKV+3211IXWQyGaEkO50ONBqNVANUl1ksFly7dk2CORvszGop\ny2VVQFUPG4IckhoMBiKR/STIL1cqFanw7Ha7fI+z2Uz+HQ9dyjdtNhsqlQq63S5OT09lLzAYu91u\noSr5HBdVRpFIRJIVHqDcv91uF7FYTA6ZTCYj/LzT6UQ4HJak1efzCTdPOSVjGfsuFotl4fOY+4x+\n7z/5fwhmJsBzCoM8Y7FYRKFQkD/HiTHSF9PpVF4MRfqkX3hqNptNNBoN2Gw2ZLNZxONxUVHMA8tG\nrkmn00GlUiESiUCv18uAxtWrV4WmMRqN0Ov1ktWy6cgSvVqtSnOTPDV/zqJS7OjoCMlkUgYJKI3U\naDQYDAYiH0smkzg8PJTDj531TCaDwWCAGzduyCQgGzRarRbBYFAmIbvdLjY3N/HDH/7wpXWsrKxI\nNktVBQCMx2MprVlSzmYz3LlzBwCElnC73RgMBojFYjg6OsLBwQHW19fhdDqlkRwMBlEsFtHpdBZm\nnOQQ+T7G4zF6vZ4MmnA/UJam0WjkoLNYLOj1eggGg3j99dext7cnygObzYbV1VVpcmq1Wlit1oUD\nNGazWaYtKSGl2onvl4NEGo1GuH3SFsDzrDWdTkuCwj3E9/vJjH4eSAFwH8xmM5GtMYin02m0Wi2s\nra1JM242m8FgMAgt53Q6hSqkRI3qFg4qAc8pvXlS0+l0Khr7brcrhxh7PlTVcFaBE8ta7fNwQnXX\njRs3pOq9TDOMRiN0Oh1R0iyqnAOBAJxOJ9LptAgXyC+zaT8YDISDpwKMw2VspNpsNpH68YC5nIWz\nwZtKpeaugyKCSqWC2Wwm8sNyuSw0qkqlEoFEJpORdVQqFTncAMhcCzl0Mg+UKnq93oVN7nl4ZQF8\neXkZrVYLFxcXMhJOvS1lXxxk4QdSq9VEa8lNkcvlMB6P4fF4sLm5KdrpcrmMer0ufC7LyHmgUuLt\nt9+G1WpFu92WF8xsxOl0Ynt7G06nE/F4HNlsVj5UNkM4ks0XdJnvZvbGKb558Pv9otRgw4QBUa/X\nw+fzSekbDAbh8XjQ6XSkSQhANPKdTgerq6soFAoydMIpVh4+i7CxsQGdTofz83N5DuRmPR4P+v0+\n8vk8XC4X7t69C7/fj2q1KiqEZDKJ8XiMZDKJeDwOq9WKs7MzyTbYAAYgw1vzEIvFkEwmcfv2bWg0\nGjgcDlQqFRQKBUynU5lmY4XT7/eRyWSws7ODeDyOO3fuwG63o9lsIhwOQ6VSodfryVSn3W6Xad9q\ntYp8Pj93HQ6HA6enp1hfX5dJXaoVmCl5PB6hMvihUX/MvVCtVtHr9WCz2aQiZMJAWaXL5ZIM+pPY\n3NyEVqtFoVCQaUy32y2VHr+JO3fuyAHBoMTBMPYjWJIbjUaZTiX9EY/HRWo5D8zeqWgC8ALFSUkk\nDyXaXVBKC0D+GZv7DodDqEYmc0x4FokPisUirly5Ij0y/kwGbap0MpkMEokEHA6HfE9er1feTS6X\ng9vthkqlgs1mkwSBtGe9XofX613Y1DUajQgEAiiVSjg9PZWJZcpZOSsRDofx4YcfYmdnR2TCnKId\nDAZC0wKATqeTv1uj0cDr9SIUCmE8Hi9U083DKwvgP/nJT3Dv3j1pQBiNRmxvb6PVagl3RB5xdXUV\n3W5XTqtsNotbt25heXkZxWJRypterydNzMlkgsPDQxHOUwc9D7dv35asTa/XY3l5GV6vF91uVwLi\n/v6+ZBPMBlUqFSqVCk5PT2VDl8tlOBwOJJNJ1Ot1ybLIGU+n04XNIZfLhZ2dHQQCAeHl6AVjMBhE\nytbpdPDHf/zHKBQKePDgAWq1GmKxGCKRCBqNBn7wgx8gFoshHA7jjTfewHg8Fv01P7CVlZWFI+xU\nljDDOzg4kMEonv6hUAj37t1Dp9PB8fGxHGIc5aYaIZPJiOql3+8Ln8kGIjnHeXjjjTfw7//+7xgM\nBhIYg8EUqkk9AAAgAElEQVSgZJzUHzNrzefzKBQK8Hq9iMfj2NvbQ6lUgtVqlWaf3+9HPp9Ho9EQ\njpgfSjAYnLsOjtGTEhmPxzJibTQaRfdLKRqnfpeWlmSa8PT0VKbzCoUCnj59KtQDFTac2F3UtLNa\nreh2u5Ix12o1JJNJuN1uBAIBSWLIMZfL5ReyZGbHzWYTsVgM29vbEripwuK0JlUd88CgwgpMpVIh\nFAqJ7pxNTp1Oh729PQyHQzQaDane+F2wX2O1WqXRWK1WxQ5iNBohkUgslO+xUiBnTl8bWgSwoful\nL30JyWRSDjHSbz6fD2trawgGg/jVr36F119/XZRDarUadrsdJpNJFFNra2tz15HNZjEYDBCNRmUY\niwcjn9fdu3dF7ki5ZyaTkWrkzp07+MUvfoFnz54hEAhIxafVauHxeGQq3Wg0Lmy2z8MrC+BskLEx\ndXZ2JsoSbmCj0YhqtSpZFbOp1dVVxONxjEYjPHnyRLJ3q9UqL54ZKCV+lUpl4WACG5VsWOj1ehmt\nbzQaqNVq8Hg8UhKS36xWq8JFr62t4f79+9jY2MDKygrS6TRyuZxMlDKjIH86Dy6XS2iYdruNZDKJ\naDQKl8uFbrcrWQrlTaenp9jf3xfKhdOqbrcbv/3tb/GpT31KVBwcVzYYDCiVSlLmLVrH2dkZisUi\nAoGAyDSpyAkGg7BYLPjHf/xHeL1efP7zn8fe3p4od3Q6HWq1GuLxON577z3cv38ff/iHfwin0ynN\nSwCitlkkz/r0pz+Nn/3sZzKlx/FvBpCVlRX4/X6YTCZ8/PHHcDqdePToEabTqfRTut0ubty4ga2t\nLVxcXMDhcGBlZUVKdI6OU/s/DxaLRcr0jY0N2Gw2nJyciJKKE5fRaBQPHz6UgbLj42OYzWbhPMvl\nMv70T/8UFxcX+O1vfyt03HQ6lV4GpwLngUoFnU6HRqOBZDKJpaUlLC0tyQF5fHwMr9eLw8NDJJNJ\nOUx8Pp9MX7IiqtVqSKfTcljTW4U0zaJSvdVq4datW0IZmEwm7O3toVarIRQKIRaLwefz4f3333+B\nChmPx+h0OpIYDQYDfOYzn4Fer8eVK1eQyWRwfn4uk6AMwItkphxF51yDzWbD8fExjo+P5XldvXoV\nmUwGo9FIBAfr6+tScW1sbKDVauELX/gCzs/P4fF4EAqFXvj5NLZaZMXRaDQQi8XkkGRF02w2Ybfb\nsbW1JUNw7FWMx2NcvXpVpJ9arRZXrlyB0WjEz3/+c9y7dw+3bt1Cv99HOp1GPp+XBnA0Gp27jnl4\npTrwYrEo3hjkLwFI0+3k5EQm2vx+P1ZXV9Fut0XVwDJRp9NJ1ms2m7G8vIxer4dsNovhcCiThItO\ncrPZLNwmuUC6zlEmZjKZxFuCPi4HBweo1WpYXl5GPp+HSqXC8fExwuEwbt++LYoPgpt9kf46FotJ\nQFOpVBiNRuKmyMnT3d1dVCoVOBwOhMNh/NVf/ZVIpKgz5fj7s2fP0Ov18M4778gocCqVEo57EXXx\nJ3/yJ/iHf/gHeQakPex2uxxirVZLhpyOjo7wz//8z3A4HOL4Rwe3t99+Wxp21LDTk+Xs7AwajQb3\n7t3Dd77znZfWMRwOsbq6KlkQADmIdDodEokEfD4fnj59ijfeeANnZ2f4gz/4A3nPbAxNJhPs7u7K\nhC/VLxxI4oexuro693l89rOfxY9//GOoVCqcn5/DbreL7w4z762tLRnW4qAQM0QmB71eDx999BF8\nPh/C4TBOTk4kQNBZU6vV4q233pq7jlu3biGTyaDVaokOvVwuC70EAKurq0in0/jtb38rsxasQDkf\nQa4ewAsBUK1WI5VKifXB9evX8eTJk5fW8ezZM1y/fl36F+PxGNFoFLFYDBaLBR6PB41GA3q9Hhsb\nG/jwww+lGTgcDvHkyROsr6+jXq8jHA5LH4dcMRvSXq8XTqdzYZM7lUrJZCKnotfX1xEIBMS4qtPp\n4MmTJ0gkEjg5OZF40Ww2cX5+jpOTE7z55pvodDoyWcsRe3oX1Wo1oV/m4fvf/74EZ9o9GAwGXLt2\nDYFAACaTCdvb22KTcHFxgStXrmA0GmF3dxenp6f43Oc+Jwexz+fD2dmZ9GcoFtBoNCgWi/jggw/m\nrmMeXlkAZ7ZWKBSEuiAVQj22RqPB0tISNBoNNjc35SNIp9OIx+OSYWWzWVEpdDodHBwcyERkrVbD\n+fm5NHXmwW634/3338drr70mLnYAxMWQCpN6vY5AIACDwYDl5WXY7XZRaPAwok6dMiKVSiXrrtVq\nMtI7D/V6HRsbG/j444+h1+tFBsmOPOmbeDwuQ0yXrTMpdaMaIBAIYDKZ4MmTJyKF4mASm2Dz8OjR\nI2SzWZGKUSni8Xhk6m86nWJ5eRndbhdPnjwRjtzpdOLWrVtotVo4OztDNBrF5uamNP1MJhNUKhUK\nhYLw14vey3/+53/inXfewZMnT6QjT/5+a2tLJij39vaws7MjnCv9N2azmThdXqbQarWaOAuyN/G/\nNVMLhQLW1taQzWbRaDTkIMzlcjCZTHj33Xelf7GysoJmsymGUGtra/LvCoWCWK7ygONzqVarQqUs\nyvR2dnZw584dnJyc4OrVqyJdJRXELNNqtcr30el0XvDR6Pf7iEQiwunzG2HQoqcLnfTm4dOf/jQu\nLi7EvQ/AC9TgaDRCvV6XAS8Of1FOZzabhZunrcJsNpP3e9n+lsF2HhhkWZmzyU01y2QyQbfbhd/v\nl0ar2WzG/v4+2u22zCDs7u5ic3NTPIX432q12heED4u8lCqVCgaDAR4+fCiJIA/UarWKjz76CLVa\nDbdu3UIsFsPa2ppU2PR3+q//+i98+ctfRiAQwOc+9zkkk0nh5llRspe1iGKbh1cWwGOxmEz8XX6B\nvV4PKpVKLF6B5w+InuH86LvdrgR9n88nJSj9R/gClpaWRPK2qFSPxWLY2NhAr9eDy+XCaDRCsVgU\no57T01OMRiNxTnM4HLDZbPB6vaIfplHPbDYTFza/3y/VBPm3/837w2w2Y2dnB263W/h/+j2MRiN4\nPB7httnEOzw8FO7PYDCIjJD83WQyQTAYxMrKimRgxWIR7XZ7YVe9VCrJUIXBYJBRXup1qbfl5Gqv\n18PVq1dx584dCeT5fF5KdA45sLKgzS3Nmb797W/PXcf5+bl04mn+RDkWGznHx8cIBALyvshx0oud\n/hipVEoqClJO5M3ZgFuU6Z2enuLatWv47//+bwDPgxWnJelCyKnOo6Mj+b3W19dF8klfDPZKWFlW\nKhUxTuL+I8U0b39QJhkKhaQqYpCx2+1wu91YXV3F0tIS/u3f/g29Xg+3b98Ww6tIJILhcIhwOIzJ\nZILNzU00m00Zf+dA1MrKysKmbj6fh9frlXep1+uFGiwWi0gmk2LP3Ol0EI1GMZ1OJaFotVqoVqvY\n2NgQqwaaSlHWyGqO1gTzwCqKlQslx/Q4YbOcB/1bb72FaDQq9AktImhLsLq6imQyKZkyKScqTN59\n99256+h0Onj06BGsVqskbUxEnzx5gn6/j7t37+LatWuSIHKQqtvt4itf+YrszfX1dUwmE+j1euRy\nORlsoosk+yC/L15ZAL927Rru37+P6XSKaDSKcDgsWk1+EMDzj4W+Fna7HRaLBaenp+LxsbKyIpkA\ntbmXPYZpB9lsNhfqWrvdLj772c/i4OBA+OHRaASVSiWBl0b7BwcH0vxgac8PiJamzIyp1aVMrN/v\nS5Y+Dxwcol7V5XIJx8sLGjhdSUtRSiXZyCI9QT0v+U+HwyHSQGqkI5GITGlexs7Ojoy708yHBkG0\nQGU2RWMlKmXYKa/VagiHw3C73SgWiyIrpMad2Z7T6VzY3afcjjKu4XAo/4xVxGQyEctbSgs5zkzv\nC07pUTXBhAGAfOTT6RQrKytz10F98cbGhjTleKDG43HY7XZ5v1tbWyiXy9BqtfIsgOfBNxqNiqyz\nVquJ1p77nZPJiyq0nZ0doegYvPmdsHrN5XJQq9UyOHRxcYHDw0Osrq4iFotJBQFAhkioZqJkkBTP\nD37wg7nruHLliiiByONTYcIeDJMxKo/8fr/IGpk40L2TvSbyzAy6lNwtOtA2NzdFdcJxfhpbkfbj\nPuJFK+xlUTNPEzxWX/TNZ5BkhUQ6Zx4ikQgePnyIL37xi/I+KGgYj8dIJBKi3Savfrm35nK5hL59\n/Pix3BnAZ8H9QYnyo0eP5q5jHl5ZAD8/P4dGo5Hpt2vXrokPMEdkKSPki2THvFAoyEQkGx0EeSi6\neFF3/fDhw4UuXhsbGzKIcHZ2JvQNG0x8yYVCASaTSXxAqGAgD8jSmTw9TXkYxFkCLWpSMQNi5hKJ\nRKTbb7VaUSgUkE6nUalUZATaZDLBZDLJoAbLLXpkcIKM3X/6adPDex6y2aw0hDk5yk3kcDikMuGk\nIQc0Go2GqII4fGOz2RAMBpHP50Unzkwlk8kszPKA5y6APAg4SXdZw8tATEqL1AMbd+FwWOR5LKuH\nw6HwwZyYm06nWF9fx/7+/tx1qNVqPHz4UMzNSCfRq5zDSaurq/KeuQ80Gg3cbjeCwaCYOjGwsXph\n4C6VSlhaWlpYqs9mM1ElsR9x2UMln8/jww8/FBVGPp8X24DZbIZyuQyj0QibzSZyQgYLfkd8Z3QR\nnIfpdCp7md8IPXYuqybobU+HQJpFqdVq0ZEzoI1GI7np6PINRv8bZUAaFoAooKjjZo+Cl8Twu6Jd\nNKWxVGhdXFxI45uyTh76lAO/9957c9dBjpoTxTy8KCf1+XyikCIdxCr2srT3cnXCi2M47MXKczwe\n480338TPfvazBV/NJ9b2e/2p/wswANKlLRKJyFgtPyq60NEch+PuhUJBMvBIJILJZCKG7G63WyYB\nAQjH+vjxY1y9enX+L/n/3uYBQG50ob8xnQBrtRoMBgPi8bi40zmdTng8Hvmz5DJTqZQYFbEEo6aV\n/iHzwJHkSCSCcrksagByzxy57nQ6cLlc4nfM20weP34sTTrSHxw2YibBrJkfxzwsLS3JdU7Xrl2T\n6TiWcPRpv+x3TT06AzizCtqw8mYRVhmVSgWNRgOHh4cLJ8voscLf8/DwEMViEZVKRT4y6sx5cQH1\nxrlcDqlU6oU+SCqVglqtFpqJFRvL5b29vYX7lWUt5WfU9HNwh43jwWCAZDIpzU2r1fpC4sDbeDi2\nzgZ8s9kUM6NF/tfMqGmaptFohApiUsFDlb40zARbrZZc8uDxeOSyCe4Jfmv8/8fHxwupC7vdLu+V\njT7qvdl7oeb+4uJCEgjOc3AAjbYM1NJTY0/HSSYBiw60Gzdu4OTkRKZYg8Eg2u32C779TBh4Q9Vl\nuTCFAaRKKH3kQUDnU5VKBZfLhYODg7nroDzQYrHIoB0Hl9gzY+JB614epLzykFPEPOAnk4mo5gaD\ngWjUx+PxQguMeXhlAbzT6SCfzyMcDkOv12N3dxcrKysyesvLEuj10O/34fP5ZOqOIn3eEMObePgx\n84qpZrMpPrqL7FP5sjnNxxKU2YLZbBapUCwWEykU8D/NG47Zh0Ih/OIXv0A2m5VszOVyCX+cSqVe\nsJq9DKoYzGYzgsGg3Ml5ubm2tbUlrnidTkcyTw7q5HI5+QBok0vZG7NWtVqNpaUl3L9/f+466PhH\n6oTeEsyKWJ6Xy2WcnZ0J70c6gxkoqSz6pJDzValUwhVfvXoVXq93bkZRqVTkZhfqxfm7cfQ7kUgg\nk8ng6OgIbrcboVBIPkzq0Nl8pDKCa+VHQ23tolK91+shFAqJRw/pCzbBaEvL8Xy73Y6VlRU5ZGw2\nmzRWDQYDqtWq9HAof6MF6eULIT4JvnM29ur1ujxrSvM4RMJeALPtUqkkvZdut4tcLodgMIhIJCKV\nDYMYs85Fo+Mej0cC5t7enmTOvGjDYrGgVqsJJ81GO51HZ7OZSFr9fj8KhQIODg7kebE5ySB+fn4+\ndx2kXeiHwr3HZ03bB3LITHT4HpjR0tOI//1lz3dWXJPJBB999NHcddAniAcFr9LjfbVmsxk3btyQ\nA61cLqPRaLxwOxF7dvwZjB2TyUR6A8PhEF6vFz/96U/nrmMeXlkAp4sfByGy2aw04lgWOhwOlEol\n0VTzI2Zwpv47GAwiGo2KYoWlCpsjKysr6PV6CzcCM91GoyG+BAzebAZSTWCxWGSog1ku7SsByEun\nYoR3abJUGwwGku1/EolEQu6sdDqdcvsQlQm87WR/fx+np6fodDp4+PChZOg0wWETj5addGwk9aBS\nqRAOhxGJRPD48eOX1sFBBDZmbTabZEd8rpe9Zlg+k57QarUSwEkj0BmQo9t8jhqNBteuXZv7PJiZ\n8CJpjjUzuFABwxuIaLxVq9WEZtDpdLh+/Tpmsxnee+89TCYT8Zzgcw2FQrh58yasViu+973vvbSO\nx48fy6FLPTMPZGbZfO+DwUAya04QVqtVRCIRuYyj3W4LrQFA/pdulIvgcrkQDoelSmXGSBUNDdTO\nz8+FnqCxGpuddA4k/US6gsHb6XTKBRc+nw//8i//8tI62HOi8ov8LBv2Z2dnQiMGg0HMZrMXPPHp\ne0IXP9o1MGli05CGdouudqP6rFwui+si+zuX79SllzytJTqdjlwCQh8bGs+x+Ul6jRXG48ePF1bO\nX/nKV/DrX/8aAESnTR9/AEJX8u7Nvb098bKnSgb4n14L7awvU4U8tPX659dN/r54ZQH8/v37ohop\nFAo4PDyUl7q+vi6X46rVapkAZKOSpR6bZPylGXR5n18mk5GR/VKptPByVPKjvEGGpylvWeF0Zb1e\nFzOZbrcrU3Gj0UiMcTj5xktuuW66vCUSiYV32rHsPD09xdrampjf0CyKG4HND6fTiW63i5OTE7mN\nJRqNwuPxSJAtl8tCAZFr5fVdi+x1TSaT8M+saC7z5xwaYhBj85mSSmq1jUYjstkszs7ORCYXCoWk\nTGVwWUQZ8KNXq9VyTRu5fE77kUNPJBJiT9psNtFqtaRMDofDOD09RblclvF2cs904ms0GgsP+Eql\nglarhXQ6LUGI/DWbgtTfVyoVHB0diZ6Xz2c2m8lEKTXiLKk5Gr+0tITV1VX85je/WbiOd955R2gh\njvIbDAaZZeAhxjkJ3rnI4ErveKpd2HNyu90yG1AqlZDL5RaqttxuN3Z3dyXDvnzpAHtIvKSc/DAr\nagCSkVMvTzqJlCP/WSKRQKFQwPXr1+c222mfO5lMxOSOE6q0SSA1w8qGtBc5eWbwVBYx2aKfOnsZ\nH374IdbW1mRI6JPfC4M/L2TIZrNSkTFGhMNhRKNRrKysSCVbKpVekNKaTCZEIhGkUimUSiWhqjio\ntbe3hzfffHPue5mHVxbAKeYvFovIZrN49uwZut0utra2hDphwFhaWpJT/HJQ5MQlJxWpQCEn3e12\nodVqxSZya2tr7lqoT240GtIQ48ugJznvblSr1Tg5ORF1ATlhlqEcr6Z0qtfrIRKJIBQK4ezsDLPZ\nDDdu3Ji7DpZtBwcHwnmtra3h6OgIuVwODodDslkA0s2mLtvn84l3MEs0TmeOx2PE43E0Gg2RSi4a\nlQ6Hw9ja2sKjR4/Q7XblVnY24agE4bPhEBEVNLTt5Kj4ZDIRBzuuhdaiGxsbCwM4b28hd02ukxcU\nGI1GaXqxaQZAmq1UOezv72N3d1cymGKxKJlzsVgU18ZFJTIrKh5QvCKv2+0KXcHLCOr1ulygATxX\nKJDj5uXdHGphL4L/H3juy7PovfAdsHqgXA2AlNmdTgfhcFguOWbPgc8vFovJIBTnF5hg6PV65PN5\n/PCHP4RKpVqYcfLwZVlPmS7pAk4jXlxcyJ6x2+2iCuIAEiss9lY4pTgajRCJRBAIBPDo0SN88Ytf\nXLg/WF0cHx8jFAoJf00hALNYVoS0tmAfoNFoiJSz3+/LFYikWkqlkiRf9+7dw89//vOX1vH06VO0\n223EYjEcHBygUqnA7/fLlDn/zvX1dYRCIZkhYGOfvDgrelZEzNKpNhoMBrhy5crcZ7EIryyA12o1\nXFxcwOPxwG634+7duzIWzIfJW7gZoJhdkBIgD86TnlwSg3kkEkE6ncbdu3fFcGYeeDD0ej3hCrkZ\nJ5MJ3G63mCBdbrrQEpdBhA0XGrGz8RUIBOD1evHkyRPUajW8/fbbC9dhsVgQjUah1WqxtrYGl8sl\nPQI2ccnZ88YeNoQZXJk5Xw6UbJgMBgNsbm7KtVnzYLFYxEC+Wq2KgRgbqSxvmSHTTY3To6wSeCky\nrXfpf8JM0OfzwWq1LlTlqNVqsQg+ODjAnTt3EI1GhcZh+U7VEgAJcFS6cFiG49jMWp1OpxwczWYT\nOzs7SCQSC9fh8XhQLBaRyWRgsVjkA/f7/XLlFiVtNCmipG84HMqe4OwCD8B2uw2z2YylpSXs7+8j\nk8ng+vXrc9cRDAbFw5xUhNFolGYiezhUqHCe4Pz8XCgWOmTa7Xa5fpABqtFoIJ1OY21tTS4MmQd6\n+vCaQcowL0sI2QSkKygbceyr0DSONCeVXLQV4OXEvLh4HtgQJuXJ4TPGBSYLtG3lgXn5+ZHy4nvm\nN0QaYzKZyBVoi3z8ScHRNnk4HMLhcAjlWy6XJZkolUoyFMigDfyPtQgTITbJuWdo3OX1ehdeRj4P\nryyAc+AmFArhc5/7nExvqVQqnJyciLEQuahWqyXuZYPBQG5cX11dhdPpFH9wBp3Lpj1XrlxBsVhc\naAzPU9ftduP4+Bg6nU78C8LhsDSP2DR1u91imQk8H9mmcxk9zrkeyu+Oj49l+nFRM3VlZUWaj51O\nB2dnZ1CpVLh7965wu3xO3LQ07L+sIY5Go6hUKsjlcuJZ3e8/v2g2EokIzbEo0ysWiwiHw1Cr1bh+\n/TouLi5kcs5msyEQCOD8/Fw05myAkW/k5mfmyTKatBNvtllbW5PAPA+vvfaaaHCB59UOaTG73Y5s\nNiuTojwkNBqNXLLMadjXXnsN29vbcvCRRuIzODs7Q7PZXOiVw8Od1SD9T2hoFovFRPtPJ0pWcPw7\nmBVS9ki/bXKuJycnSKfTci3dom+GkteNjQ3cv38fGo1GZhI49ERpHw23KGOl/pwBO5/PvxC0SKFs\nbGzg6dOn+PznP48f/ehHL60jnU5jb28PNpsNDocDoVBIBs8mk4mofC7fQMU7M5ks8e9lsKcZHZMS\nVnBLS0sLL/yg9erlqWBWfBcXF4jFYi/QOGyG80Z6DkbRJIpqMXqmkOJhnFkkPvjqV7/6wkUupVJJ\n7IV58XS1WsXFxYXcIctr99hEJ7VjNBqlsmBcGwwG8Hg84snz8ccfz13HPLyyAB6JRNBqtRCPx+Xk\n+s1vfoNcLicn0OUbtgnylpRQcRKM2XKj0Xhh5HkymWB1dRXZbHYufwVATuzLEh1OP+r1ekQiEcm4\nqA+/3BwsFAo4OTkRxQFP8Mvl4mw2e0FyOA8ul0sUDQaDQbxMbt++jWg0KnQJDaqYQZFTDYVCiEaj\ncuAB/5NV0IbS6/Xi/v37ODo6Wujt8Oabb4op15MnT/DTn/4UKysruH37NkwmE6LRqDjIMfNhkIjF\nYqKJJe/ND2c8Hstt4S6XC8lkUqyB56FSqUhZGQ6Hcf/+fZEOMhPhgcpsnRULG6UrKysiOyTPSHkl\nlU3pdFo+knlot9uitlhZWYHVakW/30e5XEaxWBQnSP6dPMwY4PlsaBHK5i7vSQyHw9jf3xcd9iLf\n+tXVVZFG2mw28cq4zOeymdvpdKSfo1KpEAwGhe4zmUxIJpNCGTKQAs/ps3w+LwF5Ht577z34fD7U\n63XRitM9k+sgdUJ1h81mE/Mw0nIMrqykuObJZIJCoYCNjQ1otdqF9rq8sYcunJdvxWGDk1m32WyW\nS7rJt1MRwiqbMxVcF3sbTMIWHayhUAgff/wxvv/974u3DqteDi4xAw+FQnI5OvcIaRYAL8y/0Cue\nfQ4qoR48eDB3HfPwygI4Lwh++vQpAoEAMpkMHA6HcHgAJFvgVCTNnkhP8JRjk4IlCtUkRqMRBwcH\neP/992UQaB7UajWOjo7En+PBgwdot9tYW1vDcDiURlitVpOyi/IobjqqMXiLyGUfhWw2i5s3b4pI\n/8tf/jK++93vvrSOZDIpIn6aalFbTlULDw6Px4ObN2+K3welWbPZ7AVKhRQCOdnt7W14PB6k0+mF\nlNIXvvAF6PV6/OpXv8JsNsOnP/1pbG5uwu/34+zsTLTm9JXweDwyLh+Px+HxeNBut5FOp+Uwovf7\nbDZDIpFAPp+XhuCi3sRljxwenrlcDhsbGzL11263EY1Gxdyfl3IwMBUKBZTLZdk3NO/nu+NU49e+\n9rWFgYIKqE6nI8NeyWQS7777Lm7evCkKB06t1mo1Caps3nHaj7Qgpw+pM47FYnj06BE2NzcXHmjk\nbOmtw0sp8vn8C3a2VFFwqIZyOipOGCSoGR+NRtIopjLn6dOnC9UO9ISnZfKzZ88wGAzwla98BVtb\nW3IfJKkWrVYrd1dy8pJ/j9lslpkJDqLxkNnf38frr7++cIKaDosc9kkmk8hkMjKVnM1mpaHe7XZf\nmPKmnYLL5XqBfmI/h15KlGPqdLq5jVQA+Na3viUH9bVr1+B2u2UfXm58Mygzm2YzWq/Xiy8UABmM\n4yFI6ovqpXkGY4vwygL4N7/5Tfz93/+9uPONx2M8fvwY77zzDu7evYvT01McHR3J5a2cUgLwQobD\n4QpeoUQZD4ctrl27hl//+tfSwZ0Haq4dDgeWlpbEy4F+I1S+aLVa8fBgpx94nq1zTPpy5568H7nj\nr371q/jOd76zcGR7fX0dP/rRj2T4YGdnR07z8Xgs8kUO5TDjYkleKpXEB/nyEBTtZOPxOEqlksgD\nyb99EvV6HbVaDV6vF2+++SY+/PBD7O3tSXXT7/fhcrmEr2cQ5MYzGo2S9Xk8HvGu4S0tsVgMarUa\nFxcXSCQSC6VzbH663W6srKxIxkk3RA5vseKhuoBNKPpbM/i73W7JNhnQeRDfu3cPNpsN3/rWt15a\nx6BQd38AAAGkSURBVOXfgVOEVqsVf/Znf4ZwOIwf/ehH8vfy0FSpVNJUZiOZz5v8N/eTyWQSFVOp\nVFpotXDz5k3867/+K65fvy6eGGwYUnlBywMGcp1OJwGFwZvctd/vl++IthDk+NkAX/S9eL1eBAIB\nVCoV7O/vw2w248qVK3jnnXfw6NEj8fzgtWuUVHLsnhLdy9OclUpFrphjFpzL5fD222/PfS83btyQ\n+1DL5bLc/nTr1i3x8z84OEA0GkUwGBQunWZiAERRRZrMZDLB7XaLfUYikRCDr0XUxQcffIBCoYC7\nd+8iGAxie3sbf/7nf44PPvgAzWZTPOgpHCgUChgOhyiVSjKZSTmjWq2WYM/Reap8vvGNb2B5eRkb\nGxtz5b8KFChQoECBAgUKFChQoECBAgUKFChQoECBAgUKFChQoECBAgUKFChQoECBAgUKFChQoECB\nAgUKFChQoECBAgUKFChQoECBAgUKFChQoECBAgUKFChQoECBAgUKFChQoECBAgUKFChQoECBAgUK\nFChQoECBAgUKFChQoECBAgUKFChQoECBgv8f4/8BoI4g4+ZEuT0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa7a8486b10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "###保存Fisherfaces的图像\n",
    "n_row = 10  ###显示多少张图像\n",
    "fig = pylab.figure()\n",
    "for i in range(n_row):\n",
    "    p = Fisherfaces[:,i].reshape((32,32)).T\n",
    "    fig.add_subplot(1, n_row, i+1)\n",
    "    pylab.imshow(p, cmap=plt.cm.gray)\n",
    "    pylab.axis('off') \n",
    "pylab.show()   \n",
    "fig.savefig('Fisherfaces.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 211,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "22.4533333333 4.09105256763\n"
     ]
    }
   ],
   "source": [
    "acc_list = []\n",
    "for i in xrange(1,51):\n",
    "    train_face, test_face, train_label, test_label = load_Yale_data(randomsplit=i)\n",
    "    mn = np.mean(train_face, keepdims=True, axis=1)\n",
    "    train_face = train_face - mn\n",
    "    test_face = test_face - mn\n",
    "    Fisherfaces = LDA(train_face, train_label, K=14)\n",
    "    ##对train数据进行LDA降维\n",
    "    lowDim_TrainData = np.matmul(Fisherfaces.T, train_face)\n",
    "    ##对test数据进行LDA降维\n",
    "    lowDim_TestData = np.matmul(Fisherfaces.T, test_face)\n",
    "    clf = neighbors.KNeighborsClassifier(algorithm='auto',n_neighbors=1, weights= 'distance')  \n",
    "    clf.fit(lowDim_TrainData.T, train_label) \n",
    "    \"\"\"测试准确率\"\"\"\n",
    "    acc = np.float(sum(test_label==clf.predict(lowDim_TestData.T))) / len(test_label)\n",
    "    acc_list.append(acc)\n",
    "print 100.0-np.mean(acc_list)*100, np.std(acc_list)*100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 187,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8\n"
     ]
    }
   ],
   "source": [
    "train_face, test_face, train_label, test_label = load_Yale_data(randomsplit=5)\n",
    "mn = np.mean(train_face, keepdims=True, axis=1)\n",
    "train_face = train_face - mn\n",
    "test_face = test_face - mn\n",
    "Fisherfaces = LDA(train_face, train_label, K=14)\n",
    "##对train数据进行LDA降维\n",
    "lowDim_TrainData = np.matmul(Fisherfaces.T, train_face)\n",
    "##对test数据进行LDA降维\n",
    "lowDim_TestData = np.matmul(Fisherfaces.T, test_face)\n",
    "clf = neighbors.KNeighborsClassifier(algorithm='auto',n_neighbors=1, weights= 'distance')  \n",
    "clf.fit(lowDim_TrainData.T, train_label) \n",
    "\"\"\"测试准确率\"\"\"\n",
    "acc = np.float(sum(test_label==clf.predict(lowDim_TestData.T))) / len(test_label)\n",
    "print acc"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}