第04章 线性分类器.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 第四章 线性分类器\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 二轮总结笔记\n",
    "\n",
    "\n",
    "### 一、线性判别函数的基本概念\n",
    "\n",
    "#### 1. 线性判别函数表达式\n",
    "\n",
    "$$\n",
    "g(\\mathbf{x}) = \\mathbf{w}^\\text{T}\\mathbf{x} + w_0\n",
    "$$\n",
    "\n",
    "#### 2. 线性判别函数的性质\n",
    "\n",
    "1. 决策面是超平面，决策面的法向量为$\\mathbf{w}$。\n",
    "\n",
    "2. 原点到决策面的距离为\n",
    "\n",
    "$$\n",
    "\\frac{|w_0|}{||\\mathbf{w}||}\n",
    "$$\n",
    "\n",
    "3. 任意一点$\\mathbf{x}$到决策面的距离为\n",
    "\n",
    "$$\n",
    "\\frac{|g(\\mathbf{x})|}{||\\mathbf{w}||}\n",
    "$$\n",
    "\n",
    "\n",
    "### 二、Fisher线性判别\n",
    "\n",
    "#### 1. 思想\n",
    "\n",
    "将高维样本投影成一维，使类内离散度最小，类间离散度最大。\n",
    "\n",
    "#### 2. 准则函数\n",
    "\n",
    "$$\n",
    "J_\\text{F}(\\mathbf{w}) = \\frac{\\mathbf{w}^\\text{T}\\boldsymbol{S}_\\text{b}\\mathbf{w}}{\\mathbf{w}^\\text{T}\\boldsymbol{S}_\\text{w}\\mathbf{w}}\n",
    "$$\n",
    "\n",
    "其中\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\mathbf{m}_1 &= \\frac{1}{N_1} \\sum_{\\mathbf{x}_i\\in\\mathscr{X}_1} \\mathbf{x}_i \\\\\n",
    "\\mathbf{m}_2 &= \\frac{1}{N_2} \\sum_{\\mathbf{x}_i\\in\\mathscr{X}_2} \\mathbf{x}_i \\\\\n",
    "\\mathbf{S}_1 &= \\sum_{\\mathbf{x}_i\\in \\mathscr{X}_1}(\\mathbf{x}_i-\\mathbf{m}_1)(\\mathbf{x}_i-\\mathbf{m}_1)^\\text{T}\\\\\n",
    "\\mathbf{S}_2 &= \\sum_{\\mathbf{x}_i\\in \\mathscr{X}_2}(\\mathbf{x}_i-\\mathbf{m}_2)(\\mathbf{x}_i-\\mathbf{m}_2)^\\text{T}\\\\\n",
    "\\mathbf{S}_\\text{w} &= \\mathbf{S}_1 + \\mathbf{S}_2 \\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "#### 3. 判别函数和判别规则\n",
    "\n",
    "##### (1) 一般情况\n",
    "\n",
    "$$\n",
    "f(\\mathbf{x}) = \\mathbf{w}^\\text{T}\\mathbf{x} + w_0\n",
    "$$\n",
    "\n",
    "##### (2) 样本为正态分布且两类协方差矩阵相同时\n",
    "\n",
    "$$\n",
    "若g(\\mathbf{x}) = \\mathbf{w}^\\text{T}(\\mathbf{x} - \\frac{1}{2}(\\mathbf{m}_1 + \\mathbf{m}_2)） \\gtrless \\ln\\frac{P(\\omega_2)}{P(\\omega_1)}，则\\mathbf{x}\\in\\begin{cases}\\omega_1\\\\ \\omega_2\\end{cases}\n",
    "$$\n",
    "\n",
    "#### 4. 投影方向$\\mathbf{w}$的计算\n",
    "\n",
    "$$\n",
    "\\mathbf{w}^* = \\mathbf{S}_\\text{w}^{-1}(\\mathbf{m}_1 - \\mathbf{m}_2)\n",
    "$$\n",
    "\n",
    "#### 5. 阈值$w_0$的计算\n",
    "\n",
    "##### (1) 样本为正态分布且两类协方差矩阵相同时\n",
    "\n",
    "$$\n",
    "w_0 = -\\frac{1}{2}(\\mathbf{m}_1+\\mathbf{m}_2)\\mathbf{S}_\\text{w}^{-1}(\\mathbf{m}_1-\\mathbf{m}_2)-\\frac{P(\\omega_2)}{P(\\omega_1)}\n",
    "$$\n",
    "\n",
    "如果样本不是正态分布，虽然不能保证是最优的，但是通常效果也很好。\n",
    "\n",
    "##### (2) 其他情况\n",
    "\n",
    "$$\n",
    "w_0 = -\\frac{1}{2}(\\mathbf{m}_1+\\mathbf{m}_2)\n",
    "$$\n",
    "\n",
    "或\n",
    "\n",
    "$$\n",
    "w_0 = -\\mathbf{m}\n",
    "$$\n",
    "\n",
    "### 三、感知器\n",
    "\n",
    "#### 1. 思想\n",
    "\n",
    "用梯度下降算法使线性判别面收敛。\n",
    "\n",
    "#### 2. 准则函数\n",
    "\n",
    "$$\n",
    "J_\\text{P}(\\boldsymbol{\\alpha}) = \\sum_{\\boldsymbol{\\alpha}^\\text{T}\\mathbf{y}_k \\leqslant 0}-\\boldsymbol{\\alpha}^\\text{T}\\mathbf{y}_k\n",
    "$$\n",
    "\n",
    "#### 3. 判别函数和判别规则\n",
    "\n",
    "$$\n",
    "若g(\\mathbf{y}) = \\mathbf{\\boldsymbol{\\alpha}}^\\text{T}\\mathbf{y} \\gtrless 0，则\\mathbf{x}\\in\\begin{cases}\\omega_1\\\\ \\omega_2\\end{cases}\n",
    "$$\n",
    "\n",
    "#### 4. 性质\n",
    "\n",
    "只能解决一阶谓词逻辑问题（AND和OR）。\n",
    "\n",
    "#### 5. 算法\n",
    "\n",
    "每次考察一个样本，若$\\mathbf{\\boldsymbol{\\alpha}}^\\text{T}\\mathbf{y}_j \\leqslant 0$，则按梯度$\\nabla J_\\text{P}(\\boldsymbol{\\alpha}) = -\\mathbf{y}_j$更新权值$\\boldsymbol{\\alpha}$，即$\\boldsymbol{\\alpha} = \\boldsymbol{\\alpha} + \\mathbf{y}_j$，直到收敛。\n",
    "\n",
    "\n",
    "### 四、最小平方误差判别（MSE）\n",
    "\n",
    "#### 1. 思想\n",
    "\n",
    "通过最小平方误差逼近来最小化错分样本。\n",
    "\n",
    "#### 2. 准则函数\n",
    "\n",
    "$$\n",
    "J_\\text{S}(\\boldsymbol{\\alpha}) = ||\\mathbf{Y}\\boldsymbol{\\alpha} - \\mathbf{b}||^2 = \\sum_{i=1}^N (\\boldsymbol{\\alpha}^\\text{T}\\mathbf{y}_i - b_i)^2\n",
    "$$\n",
    "\n",
    "#### 3. 算法\n",
    "\n",
    "##### (1) 伪逆法\n",
    "\n",
    "用$\\mathbf{Y}$的伪逆求$\\boldsymbol{\\alpha}$\n",
    "\n",
    "$$\n",
    "\\boldsymbol{\\alpha}^* = (\\mathbf{Y}^\\text{T}\\mathbf{Y})^{-1}\\mathbf{Y}^\\text{T}\\mathbf{b} = \\mathbf{Y}^+\\mathbf{b}\n",
    "$$\n",
    "\n",
    "##### (2) LMS算法\n",
    "\n",
    "即根据准则函数的梯度$\\nabla J_\\text{S}(\\boldsymbol{\\alpha}) = 2\\mathbf{Y}^\\text{T}(\\mathbf{Y}\\boldsymbol{\\alpha} - \\mathbf{b})$进行梯度下降直到收敛。\n",
    "\n",
    "#### 4. $\\mathbf{b}$的取值\n",
    "\n",
    "##### (1) 同一类样本的$b_i$是相同值\n",
    "\n",
    "等价于Fisher线性判别，即将样本和权向量都还原成增广以前的形式后有\n",
    "\n",
    "$$\n",
    "\\mathbf{w}^* \\propto \\mathbf{S}_\\text{w}^{-1}(\\mathbf{m}_1 - \\mathbf{m}_2)\n",
    "$$\n",
    "\n",
    "##### (2) 第一类样本对应的$b_i=N/N_1$，第二类样本对应的$b_i=N/N_2$\n",
    "\n",
    "在等价于Fisher线性判别的基础上，阈值$w_0^*$是全体样本均值在判别函数方向的投影\n",
    "\n",
    "$$\n",
    "w_0^* = -m^\\text{T}\\mathbf{w}^*\n",
    "$$\n",
    "\n",
    "##### (3) $b_i$全为1，样本数量$N\\rightarrow\\infty$\n",
    "\n",
    "MSE是贝叶斯估计的最小平方误差逼近。\n",
    "\n",
    "贝叶斯判别函数为\n",
    "\n",
    "$$\n",
    "g_0(\\mathbf{x}) = P(\\omega_1|\\mathbf{x}) - P(\\omega_2|\\mathbf{x})\n",
    "$$\n",
    "\n",
    "最小平方误差逼近是指平方误差的期望最小，即\n",
    "\n",
    "$$\n",
    "\\int [\\boldsymbol{\\alpha}^\\text{T}\\mathbf{y} - g_0(\\mathbf{x})]^2p(\\mathbf{x})\\text{d}\\mathbf{x}\n",
    "$$\n",
    "\n",
    "在$\\boldsymbol{\\alpha} = \\mathbf{Y}^+\\mathbf{1}_N$时取得最小值。\n",
    "\n",
    "\n",
    "### 五、线性支持向量机（Linear SVM）\n",
    "\n",
    "#### 1. 决策函数\n",
    "\n",
    "$$\n",
    "f(\\mathbf{x}) = \\text{sgn}(g(\\mathbf{x})) = \\text{sgn}(\\mathbf{w} \\cdot \\mathbf{x} + b)\n",
    "$$\n",
    "\n",
    "#### 2. 原问题（primary problem）\n",
    "\n",
    "##### (1) 线性可分情况下\n",
    "\n",
    "$$\\begin{aligned}\n",
    "& \\min\\limits_{\\mathbf{w},b}\\frac{1}{2}||\\mathbf{w}||^2 \\\\\n",
    "&\\text{s.t.}\\quad y_i(\\mathbf{w}\\cdot\\mathbf{x}_i+b)-1 \\geqslant 0, \\quad i = 1,2,\\cdots,N\n",
    "\\end{aligned}$$\n",
    "\n",
    "##### (2) 线性不可分情况下\n",
    "\n",
    "$$\\begin{aligned}\n",
    "& \\min\\limits_{\\mathbf{w},b}\\frac{1}{2}||\\mathbf{w}||^2 + C\\sum_{i=1}^N \\xi_i \\\\\n",
    "&\\text{s.t.}\\quad y_i(\\mathbf{w}\\cdot\\mathbf{x}_i+b)-1 + \\xi_i \\geqslant 0,  \\quad \\xi_i\\geqslant 0,\\quad i = 1,2,\\cdots,N\n",
    "\\end{aligned}$$\n",
    "\n",
    "#### 3. 对偶问题（dual problem）\n",
    "\n",
    "##### (1) 线性可分情况下\n",
    "\n",
    "$$\\begin{aligned}\n",
    "&\\max\\limits_{\\boldsymbol{\\alpha}} Q(\\boldsymbol{\\alpha}) = \\sum_{i=1}^N\\alpha_i - \\frac{1}{2}\\sum_{i,j=1}^N\\alpha_i\\alpha_jy_iy_j(\\mathbf{x}_i\\cdot\\mathbf{x}_j)\\\\\n",
    "&\\text{s.t.} \\sum_{i=1}^Ny_i\\alpha_i = 0 \\\\\n",
    "&\\alpha_i \\geqslant 0,\\quad i=1,2,\\cdots,N\n",
    "\\end{aligned}$$\n",
    "\n",
    "##### (2) 线性不可分情况下\n",
    "\n",
    "$$\\begin{aligned}\n",
    "&\\max\\limits_{\\boldsymbol{\\alpha}} Q(\\boldsymbol{\\alpha}) = \\sum_{i=1}^N\\alpha_i - \\frac{1}{2}\\sum_{i,j=1}^N\\alpha_i\\alpha_jy_iy_j(\\mathbf{x}_i\\cdot\\mathbf{x}_j)\\\\\n",
    "&\\text{s.t.} \\sum_{i=1}^Ny_i\\alpha_i = 0 \\\\\n",
    "&0 \\leqslant \\alpha_i \\leqslant C,\\quad i=1,2,\\cdots,N\n",
    "\\end{aligned}$$\n",
    "\n",
    "其中，$C$较小表示对错误比较容忍而更强调对于正确分类的样本的分类间隔，$C$较大表示更强调对分类错误的惩罚。\n",
    "\n",
    "#### 4. 支持向量\n",
    "\n",
    "##### (1) 线性可分情况下\n",
    "\n",
    "满足$\\alpha_i > 0$的样本点是支持向量，支持向量满足$y_i(\\mathbf{w}^*\\cdot\\mathbf{x}_i + b^*) - 1 = 0$。\n",
    "\n",
    "##### (2) 线性不可分情况下\n",
    "\n",
    "满足$0 < \\alpha_i < C$的样本点是边界支持向量，边界支持向量满足$\\xi_i = 0$，$y_i(\\mathbf{w}^*\\cdot\\mathbf{x}_i + b^*) - 1 = 0$\n",
    "\n",
    "满足$\\alpha_i = C$的样本点是错分支持向量，错分支持向量满足$\\xi_i > 0$，$y_i(\\mathbf{w}^*\\cdot\\mathbf{x}_i + b^*) - 1 < 0$\n",
    "\n",
    "支持向量具体位置：\n",
    "\n",
    "$$\n",
    "\\begin{cases}\n",
    "0 < \\alpha_i < C, \\xi_i=0:支持向量落在间隔面上\\\\\n",
    "\\alpha_i = C ,\\begin{cases}\n",
    "           0<\\xi_i<1:支持向量落在间隔面和分类面之间 \\\\\n",
    "           \\xi_i=1:支持向量落在分类面上\\\\\n",
    "           \\xi_i>1:支持向量落在分类面错分的对侧\n",
    "         \\end{cases}\\\\\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "#### 5. 权向量$\\mathbf{w}$解法\n",
    "\n",
    "$$\n",
    "\\mathbf{w}^* = \\sum_{i=1}^N\\alpha_i^*y_i\\mathbf{x}_i\n",
    "$$\n",
    "\n",
    "#### 6. 阈值$b$解法\n",
    "\n",
    "$$\n",
    "y_i(\\mathbf{w}^*\\cdot\\mathbf{x}_i+b^*) - 1 = 0\n",
    "$$\n",
    "\n",
    "##### (1) 线性可分情况下\n",
    "\n",
    "用所有的**支持向量**用上面的公式求$b$，再平均。\n",
    "\n",
    "\n",
    "##### (2) 线性不可分情况下\n",
    "\n",
    "用所有的**边界支持向量**用上面的公式求$b$，再平均。\n",
    "\n",
    "#### 7. SVM的性质\n",
    "\n",
    "1. SVM所得的分类超平面使VC维最小。\n",
    "2. 最大化分类间隔提高推广能力。\n",
    "3. 支持向量代表样本集全部分类信息。\n",
    "\n",
    "\n",
    "### 六、多类线性分类器\n",
    "\n",
    "#### 1. 多个两类分类器的组合\n",
    "\n",
    "##### (1) 一对多\n",
    "\n",
    "1. 用$c-1$个两类分类器实现分类\n",
    "2. 优点：分类器少\n",
    "3. 缺点：样本数量不均衡、歧义区域大\n",
    "\n",
    "##### (2) 逐对分类\n",
    "\n",
    "1. 用$c(c-1)/2$个两类分类器实现分类\n",
    "2. 优点：样本数量均衡、歧义区域小\n",
    "3. 缺点：分类器多\n",
    "\n",
    "##### (3) 比较线性判别函数的阈值\n",
    "\n",
    "1. 用$c$个两类分类器实现分类\n",
    "2. 思想：比较每个分类器判别函数的输出值\n",
    "2. 问题：不同分类器训练后的权值和阈值尺度不相同，线性判别函数的输出值可能没有绝对值可比性\n",
    "\n",
    "#### 2. 多类线性判别函数\n",
    "\n",
    "##### (1) 判别函数\n",
    "\n",
    "$$\n",
    "g_i(\\mathbf{x}) = \\mathbf{w}_i^{\\text{T}}\\mathbf{x} + w_{i0}, \\quad i=1,2,\\cdots,N\n",
    "$$\n",
    "\n",
    "##### (2) 判别规则\n",
    "\n",
    "$$\n",
    "若g_i(\\mathbf{x})>g_j(\\mathbf{x}),\\forall j\\neq i, \\quad 则\\mathbf{x}\\in \\omega_i\n",
    "$$\n",
    "\n",
    "##### (3) 实现算法\n",
    "\n",
    "类似于感知器的梯度下降（逐步修正法）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## 一轮学习笔记（包含代码实现）\n",
    "\n",
    "线性分类器的决策面是超平面，决策函数是线性函数。\n",
    "\n",
    "有两类情况需要分开考虑：线性可分和线性不可分。\n",
    "\n",
    "首先考虑线性可分的情况，生成数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.svm import SVC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "X是样本矩阵，每**行**是一个样本。\n",
    "\n",
    "理论推导的样本是列向量，写代码的时候为了方便一般用行向量。\n",
    "\n",
    "real_class是类别，分为1和-1两类。前二十个是1类，后二十个是-1类。\n",
    "\n",
    "生成的这组数据是线性可分的。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.99720813 -1.55464404]\n",
      " [-2.15332905 -0.79988732]\n",
      " [-1.6271845  -1.89607447]\n",
      " [-2.88804683 -0.72722164]\n",
      " [-0.79871052 -2.86639942]\n",
      " [-1.9454445  -1.2659932 ]\n",
      " [-2.63452428 -1.21045949]\n",
      " [-2.28814449 -3.60491927]\n",
      " [-0.14831176 -1.88885453]\n",
      " [-1.31782927 -1.09599919]\n",
      " [-3.50049996 -1.19430091]\n",
      " [-2.6730074  -3.10028377]\n",
      " [-2.34421754 -4.356574  ]\n",
      " [-1.17966467  0.04991699]\n",
      " [-2.38234616 -1.87415735]\n",
      " [-3.07994777 -1.70599279]\n",
      " [-1.41432354 -2.58620074]\n",
      " [-1.96619849 -2.00863597]\n",
      " [-2.14275402 -3.29307151]\n",
      " [-2.44862645 -1.26654069]\n",
      " [ 3.17028293  2.87021525]\n",
      " [ 1.72407696  2.27430904]\n",
      " [ 1.94459064  2.02676625]\n",
      " [ 2.95914397  0.92932175]\n",
      " [ 2.70985559  1.76582962]\n",
      " [ 2.71996217  1.25905597]\n",
      " [ 0.40617109  1.47900559]\n",
      " [ 1.718164    1.57826132]\n",
      " [ 0.52405299  1.09889933]\n",
      " [ 2.56521233  2.62989855]\n",
      " [ 1.1198243   1.04418024]\n",
      " [ 1.86167209  1.17372623]\n",
      " [ 0.09299459  2.37263639]\n",
      " [ 1.42149269  2.73232963]\n",
      " [ 1.44622845  3.0853286 ]\n",
      " [ 1.63555418  0.77169289]\n",
      " [ 1.62201784  2.85085454]\n",
      " [ 1.97566001  1.33620458]\n",
      " [ 1.4990023   2.388804  ]\n",
      " [ 3.14544988  0.89199308]]\n",
      "[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1]\n"
     ]
    }
   ],
   "source": [
    "X = np.r_[np.random.randn(20, 2) - [2.0, 2.0], np.random.randn(20, 2) + [2.0, 2.0]]\n",
    "real_class = 20 * [1] + 20 * [-1]\n",
    "\n",
    "print(X)\n",
    "print(real_class)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面画的图中蓝色是-1样本，红色是1样本。\n",
    "\n",
    "由图像可以明显看出来是线性可分的。\n",
    "\n",
    "（这个数据是随机生成的，如果你发现你生成的是线性不可分的，就再运行几次重新生成数据，因为下面的实验必须保证线性可分）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGdCAYAAAAvwBgXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAlBUlEQVR4nO3dX2xUZf7H8c8w2gJCZ8WtWJipBUzcGFfcoChq99cqq7sXWjKpuxGj6BKzGNhA2LhSEyVemLorWTGEKIkueCH+SR1k/4QVgy1043+UrMsuJLAlLW2RCnEGezE10/O7OE7p0Omfmc6Z55wz71cymczpgXlm6u75cJ7v830ClmVZAgAAMGCK6QEAAIDSRRABAADGEEQAAIAxBBEAAGAMQQQAABhDEAEAAMYQRAAAgDEEEQAAYMxFpgcwlsHBQfX09GjmzJkKBAKmhwMAACbAsiydO3dOc+bM0ZQpY9/zcHUQ6enpUSQSMT0MAACQh66uLoXD4THPcXUQmTlzpiT7g1RUVBgeDQAAmIhEIqFIJDJ0HR+Lq4NIejqmoqKCIAIAgMdMpKyCYlUAAGAMQQQAABhDEAEAAMYQRAAAgDEEEQAAYAxBBAAAGEMQAQAAxhBEAACAMa5uaAYAcIdUSmpvl3p7paoqqbZWCgZNjwp+QBABAIwpFpPWrpVOnjx/LByWXnhBikbNjQv+wNQMAGBUsZjU2JgZQiSpu9s+HouZGRf8gyACAMgqlbLvhFjWyJ+lj61bZ58H5IsgAgDIqr195J2Q4SxL6uqyzwPy5WgQefHFF3XdddcN7Z67ZMkS7dmzx8m3BAAUSG9vYc8DsnE0iITDYT377LM6ePCgPvvsM91+++1qaGjQ4cOHnXxbAEABVFUV9jwgm4BlZZv9c86sWbP03HPPaeXKleOem0gkFAqFFI/HVVFRUYTRAQDSUimppsYuTM12pQgE7NUzHR0s5UWmXK7fRasRSaVSeuONN9Tf368lS5ZkPSeZTCqRSGQ8AABmBIP2El3JDh3DpV9v3kwIweQ4HkS+/PJLzZgxQ+Xl5Vq1apV27dqla665Juu5zc3NCoVCQ49IJOL08AAAY4hGpZYWae7czOPhsH2cPiKYLMenZgYGBtTZ2al4PK6Wlha9/PLL2r9/f9YwkkwmlUwmh14nEglFIhGmZgDAMDqr+o+Tv9NcpmaKXiOydOlSLViwQNu2bRv3XGpEAMCbCC7u5nS3XFfWiKQNDg5m3PUAAPhLLGYXudbXS8uX2881NXRhdQu3dct1NIg0NTXpwIEDOnHihL788ks1NTWpra1N999/v5NvCwAwxG0XOWRyY7dcR4PI6dOn9eCDD+rqq6/WHXfcoU8//VTvvvuufvaznzn5tgAAA9x4kUMmN3bLdXT33VdeecXJvx4A4CK5XOTq6oo2LAzjxm65jgYRAIC35VJ06saLHDK5sVsum94BALLKtejUjRc5ZKqttVfHXNigLi0QkCIR+7xiIYgAAEbIp+jUjRc5ZHJjt1yCCAAgQ75Fp268yGEkt3XLJYgAADJMZmWF2y5yyC4alU6ckFpbpZ077eeODjO/H4pVAQAZJlt0Go1KDQ2l3VnVC51lg0F3rF4iiAAAMhSi6NQtFzkTnG6f7jdMzQAAMlB0mj86y+aOIAIAyEDRaX7oLJsfgggAYASKTnPnxvbpXkCNCAAgK4pOc0Nn2fwQRAAAoyrlotNc0Vk2P0zNAABQABT55ocgAgBAAVDkmx+CCAAABUKRb+6oEQEAoIAo8s0NQQQAkDMvtDA3iSLfiSOIAPAsLoZm0MIchUSNCABPisWkmhqpvl5avtx+rqmhhbbTaGGOQiOIAPAcLoZm0MIcTiCIAPAULobm0MIcTiCIAPAULobmFKqFeSoltbVJr79uPxMaSxvFqgA8hf08zClEC3MKXXEh7ogA8BT28zBnsi3Mqe1BNgQRAJ7Cfh7mTKaFObU9GA1BBICnsJ+HWfm2MKe2B6MhiADwHPbzMCsalU6ckFpbpZ077eeOjrG/d2p7MBqKVQF4Evt5mJVrC3NqezAagggAz2I/D+9I1/Z0d2evEwkE7J9T21N6mJoBADiO2h6MhiACACgKanuQDVMzAICiobYHFyKIAACKajK1PakUIcZvCCIAAE+gPbw/OVoj0tzcrBtvvFEzZ87U5ZdfrmXLluno0aNOviUAwIdoD+9fjgaR/fv3a/Xq1froo4/03nvv6bvvvtOdd96p/v5+J98WAOAjtIf3t4BlZfvVOqOvr0+XX3659u/fr5/+9Kfjnp9IJBQKhRSPx1VRUVGEEQIA3KatTaqvH/+81lb6yrhFLtfvotaIxONxSdKsWbOy/jyZTCqZTA69TiQSRRkXAMC9aA/vb0XrIzI4OKh169bp1ltv1bXXXpv1nObmZoVCoaFHJBIp1vAAAC5Fe3h/K9rUzKOPPqo9e/bon//8p8LhcNZzst0RiUQiTM0AQAlLpaSamvHbw3d0sJTXLXKZminKHZE1a9bob3/7m1pbW0cNIZJUXl6uioqKjAcAoLTRHt7fHA0ilmVpzZo12rVrl95//33NmzfPybcDAPgU7eH9y9Fi1dWrV2vnzp3avXu3Zs6cqVOnTkmSQqGQpk2b5uRbAwB8hvbw/uRojUjgwnto39u+fbseeuihcf88y3cBAPAe1yzfLWKLEgAAPIe9c9hrBgAAI9g7x1a0PiIAAMDG3jnnEUQAACgi9s7JRBABAKCI2ttH3gkZzrKkri77vFJAEAEAoIjYOycTQQQAgCJi75xMBBEAAIqottZeHTNKqy0FAlIkYp9XCggiAAAUEXvnZCKIAACKKpWS2tqk11+3n0tldchw7J1zHg3NAABFQxOv89g7x+boXjOTxV4zAOAf6SZeF1510tMRpXYnwM9yuX4zNQMAcBxNvDAagggADEP9gjNo4oXRUCMCAN+jfsE5NPHCaLgjAgBiEzKn0cQLoyGIACh51C84jyZeGA1BBEDJo37BeTTxwmgIIgBKHvULxUETL2RDsSqAkkf9QvEUqolXKkUjML8giAAoeen6he7u7HUigYD9c+oXCiMYlOrq8v/zrG7yF6ZmAJQ86he8g9VN/kMQAQBRv+AFrG7yJ6ZmAOB7bELmbrmsbprM1E8xUetCEAGADJOtX4Bz/La6iVoXG1MzAABP8NPqJmpdziOIAAA8wS/dWal1yUQQAQC4TrZdkP2yuolOvpkIIgAAV4nFpJoaqb5eWr7cfq6psY/7YXWT32pdJotiVQCAa6RrJy6ctkjXTqTDhpdXN/mp1qUQApaVbZbKHRKJhEKhkOLxuCoqKkwPBwBYbumgVMq+8zHatEW6w21Hh7e/8/TnHK+Tr5c/Zy7Xb6ZmAGCCxpoywOSVSu2EX2pdCoUgAgATwHJL55VS7YQfal0KhRoRABjHeMstAwF7uWVDg/l/xXp56qjUaie8XutSKAQRABiHV1qLe71TZynugkwnX4enZg4cOKC7775bc+bMUSAQ0DvvvOPk2wGAI7wwZeCHqSNqJ0qTo0Gkv79fCxcu1NatW518GwBwlNunDPzUqZPaidJTtOW7gUBAu3bt0rJlyyb8Z1i+C8AN3L7csq3NXsEzntZW70wDeLnWBbldv6kRAYBxpKcMGhvt0DE8jLhhysALU0e5onZiYvwQ2Fy1fDeZTCqRSGQ8AMAN3Dxl4PapIzjDL31tXBVEmpubFQqFhh6RSMT0kABgSDQqnThhT3Hs3Gk/d3SYr1vwy660mLjJFCdn21DQJFfViCSTSSWTyaHXiURCkUiEGhEAGEf6wiRlnzoyfdcGhTOZVvjFWuLt2Rbv5eXlqqioyHgAAMbn5qkjFFa+rfDdusTb0WLVb7/9VseOHRt63dHRoUOHDmnWrFmqrq528q0BoOTQqbM05FOc7ObuwI4Gkc8++0z1w9aUrV+/XpK0YsUK7dixw8m3BoCSxGoT/8unONnN3YEdDSJ1dXUqUgkKAMBl/LC01I3yaYXv5iXerqoRAQD4g1+WlrpRPq3w3bzEmyACACgotxZF+kmuxcluXuJdtOW7+aDFOwB4y2SWliJ3uUx/FXOJt2eX7wIAvC3fpaXIT7o4+b777Oexwp1bl3iz1wwAoGDcXBQJdy7xJogAAArGzUWRsLltiTdTMwCAgnFzUSTciSACACiYfJaWorQRRAAABeXWoki4EzUiAICCc2NRJNyJIAIAcITbiiLhTkzNAAAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMuMj0AoCSkUlJ7u9TbK1VVSbW1UjBoelQAYBxBBHBaLCatXSudPHn+WDgsvfCCFI2aGxcAuABTM4CTYjGpsTEzhEhSd7d9PBYzMy4AcImiBJGtW7eqpqZGU6dO1U033aRPPvmkGG8LmJVK2XdCLGvkz9LH1q2zzwOAEuV4EHnzzTe1fv16bdy4UZ9//rkWLlyou+66S6dPn3b6rQGz2ttH3gkZzrKkri77PAAoUY4HkT/96U965JFH9PDDD+uaa67RSy+9pOnTp+vPf/6z028NmNXbW9jzAMCHHA0iAwMDOnjwoJYuXXr+DadM0dKlS/Xhhx+OOD+ZTCqRSGQ8AM+qqirseQDgQ44Gka+//lqpVEqzZ8/OOD579mydOnVqxPnNzc0KhUJDj0gk4uTwAGfV1tqrYwKB7D8PBKRIxD4PAEqUq1bNNDU1KR6PDz26urpMDwnIXzBoL9GVRoaR9OvNm+knAqCkORpEfvjDHyoYDOqrr77KOP7VV1/piiuuGHF+eXm5KioqMh6Ap0WjUkuLNHdu5vFw2D5OHxEAJc7RIFJWVqZFixZp3759Q8cGBwe1b98+LVmyxMm3BtwjGpVOnJBaW6WdO+3njg5CCACoCJ1V169frxUrVuiGG27Q4sWLtXnzZvX39+vhhx92+q0B9wgGpbo606MAANdxPIj86le/Ul9fn5566imdOnVK119/vf7xj3+MKGAFAAClJ2BZ2do+ukMikVAoFFI8HqdeBAAAj8jl+u2qVTMAAKC0EEQAAIAxBBEAAGAMQQQAABhDEAEAAMYQRAAAgDEEEQAAYAxBBAAAGEMQAQAAxhBEAACAMQQRAABgDEEEAAAYQxABAADGEEQAAIAxBBEAAGAMQQQAABhDEAEAAMYQRAAAgDEEEQAAYAxBBAAAGEMQAQAAxlxkegAYRSoltbdLvb1SVZVUWysFg6ZHBQBAQRFE3CgWk9aulU6ePH8sHJZeeEGKRs2NCwCAAmNqxm1iMamxMTOESFJ3t308FjMzLgAAHEAQcZNUyr4TYlkjf5Y+tm6dfZ5XpFJSW5v0+uv2s5fGDgBwHEHETdrbR94JGc6ypK4u+zwviMWkmhqpvl5avtx+rqnhrg4AYAhBxE16ewt7nklMMQEAJoAg4iZVVYU9zxQ/TjEBABxBEHGT2lp7dUwgkP3ngYAUidjnuZnfppgAAI4hiLhJMGgv0ZVGhpH0682b3d9PxE9TTAAARxFE3CYalVpapLlzM4+Hw/ZxL/QR8csUEwDAcQHLyjaR7w6JREKhUEjxeFwVFRWmh1NcXu6smkrZq2O6u7PXiQQCdrDq6PDOZwIATFgu1286q7pVMCjV1ZkeRX7SU0yNjXboGB5GvDTFBABwHFMzcIYfppgAAI5z7I7IM888o7///e86dOiQysrK9M033zj1VnCraFRqaPDuFBMAwHGOBZGBgQHde++9WrJkiV555RWn3gZu5+UpJgCA4xwLIk8//bQkaceOHU69BQAA8DhXFasmk0klk8mh14lEwuBogAnw8uomAHABVxWrNjc3KxQKDT0ikYjpIQGjY1M/AJi0nILIhg0bFAgExnwcOXIk78E0NTUpHo8PPbq6uvL+uwBHsakfABRETlMzv/vd7/TQQw+Nec78+fPzHkx5ebnKy8vz/vNAUYy3qV8gYG/q19DANA0AjCOnIFJZWanKykqnxgJ4Qy6b+rFiCADG5Fixamdnp86ePavOzk6lUikdOnRIknTVVVdpxowZTr3txJRKgWGpfM5iY1M/ACgYx4LIU089pVdffXXo9U9+8hNJUmtrq+pM/isxFrNvqw//F204bLck91O3z1L5nCawqR8AFExpbXqXLjC88COn9z/xS+vxUvmcprCpHwCMKZfrd+kEkfTFY7S5fb9cPCb7Od0yneOWcYwmHfak7Jv6EfYAlLBcrt+u6iPiqFwKDL1sMp/TLX0x3DKOsbCpHwAUROkEkVIpMMz3c7qlL4ZbxjER0ah04oTU2irt3Gk/d3QQQgAgB6UTREqlwDCfzzleXwzJ7ouRSk16eGNyyzhykd7U77777Gc3TR8BgAeUThCprbVvm6fn8C8UCEiRiH2el+XzOd0ybeWWcQAAiqZ0gkgwaC9dlUZepNOvN2/2/r9o8/mcbpm2css4AABFUzpBRCqdAsNcP6dbpq3cMg4AQNGUzvLd4dy+NLRQJvo53dIXwy3jAABMSi7Xb8c6q7pausDQ7yb6OdPTOY2N9sU+W1+MYkxbuWUcAICiKa2pGYzOLdNWbhkHAKAoSnNqBqNzy7SVW8YBAMgZUzPIn1umrdwyDgCAo5iaAQAAxnBHBHAjpqYAlAiCCOA2sZjd6n54l9lw2F5RRLEuAJ9hagZwEy9t+gcABUAQAdzCi5v+AcAkEUQAt2DTPwAliCACuAWb/gEoQQQRwC3Y9A9ACWLVTDGwFBMTUVtrr44Zb9O/2trijw0AHMIdEafFYvaOsvX10vLl9nNNDasfMFJ60z/p/CZ/aWz6B8CnCCJOYikmcsWmfwBKDJveOSWVsu98jLYKIn2bvaODf+FiJKbzAHgYm965QS5LMdncDRdi0z8AJYKpGaewFBMAgHERRJzCUkwAAMbF1IxT/LIUk1oFAICDCCJOSS/FbGy0Q8fwMOKVpZh+2wWWUOUu/D4AiKkZZ3l5Kabflh7Tz8Vd+H0A+B7Ld4vBa//y89vS43SouvA/9fSdKbeHQr/h9wH4Xi7Xb4IIRmprs/+FOp7WVvcvMfVbqPI6fh9AScjl+s3UDEby09LjXPq5wHn8PgBcgGJVjGRy6XGhp7H8FKr8gN8HgAs4dkfkxIkTWrlypebNm6dp06ZpwYIF2rhxowYGBpx6SxRKeunxhRuvpQUCUiRS+KXHThQw0s/FXfh9ALiAY0HkyJEjGhwc1LZt23T48GE9//zzeumll/TEE0849ZYoFBO7wDq1SsdUqEJ2/D4AXKCoxarPPfecXnzxRf3vf/+b0PkUqxqWrY9IJGKHkEKuanC6gDEdcqTs/VxYpVFc/D4A33NtsWo8HtesWbNG/XkymVQikch4wKBoVDpxwl4ds3On/dzRUfiLhNMFjF7u5+JH/D4ADFO0YtVjx45py5Yt2rRp06jnNDc36+mnny7WkDARxdgFthgFjNGo1NDgrX4ufsbvA8D3cp6a2bBhg/7whz+Mec5///tf/ehHPxp63d3drf/7v/9TXV2dXn755VH/XDKZVDKZHHqdSCQUiUSYmvE7P/UtAQA429Csr69PZ86cGfOc+fPnq6ysTJLU09Ojuro63XzzzdqxY4emTJn4bBA1IiUiXSMy3gaBNLkCAE/I5fqd89RMZWWlKisrJ3Rud3e36uvrtWjRIm3fvj2nEIIS4ocNAgEAeXEsGXR3d6uurk7V1dXatGmT+vr6dOrUKZ06dcqpt4SXUcAIACXJsWLV9957T8eOHdOxY8cUDoczfubi7W1gEgWMAFBy2PQOAAAUlGv7iAAAAAxHEAEAAMYQRAAAgDEEEQAAYAxBBAAAGEMQAQAAxhBEAACAMQQRAABgDEEEAAAYQxABAADGEEQAAIAxBBEAAGAMQQQAABhDEAEAAMYQRAAAgDEXmR4A4AmplNTeLvX2SlVVUm2tFAyaHhUAeB5BBBhPLCatXSudPHn+WDgsvfCCFI2aGxcA+ABTM8BYYjGpsTEzhEhSd7d9PBYzMy4A8AmCCDCaVMq+E2JZI3+WPrZunX0eACAvBBFgNO3tI++EDGdZUleXfR4AIC8EEWA0vb2FPQ8AMAJBBBhNVVVhzwMAjEAQAUZTW2uvjgkEsv88EJAiEfs8AEBeCCLAaIJBe4muNDKMpF9v3kw/EQCYBIIIMJZoVGppkebOzTweDtvH6SMCAJNCQzNgPNGo1NBAZ1UAcABBBJiIYFCqqzM9CgDwHYII4BfshwPAgwgigB+wHw4Aj6JYFfA69sMB4GEEEcDL2A8HgMcRRAAvYz8cAB5HEAG8jP1wAHgcQQTwMvbDAeBxjgaRe+65R9XV1Zo6daqqqqr0wAMPqKenx8m3BEoL++EA8DhHg0h9fb3eeustHT16VG+//baOHz+uxsZGJ98SKC3shwPA4wKWla3c3hl/+ctftGzZMiWTSV188cXjnp9IJBQKhRSPx1VRUVGEEQIela2PSCRihxD6iAAoslyu30VraHb27Fm99tpruuWWW0YNIclkUslkcuh1IpEo1vAAb2M/HAAe5Xix6uOPP65LLrlEl112mTo7O7V79+5Rz21ublYoFBp6RCIRp4fnT6mU1NYmvf66/UwPidKQ3g/nvvvsZ0IIAA/IOYhs2LBBgUBgzMeRI0eGzn/sscf0xRdfaO/evQoGg3rwwQc12mxQU1OT4vH40KOrqyv/T1aqYjGppkaqr5eWL7efa2rorgkAcKWca0T6+vp05syZMc+ZP3++ysrKRhw/efKkIpGIPvjgAy1ZsmTc96JGJEfpVt8X/krTRYstLdQLAAAc52iNSGVlpSorK/Ma2ODgoCRl1IGgQMZr9R0I2K2+Gxq4ZQ8AcA3HilU//vhjffrpp7rtttt06aWX6vjx43ryySe1YMGCCd0NQY5yafVdV1e0YQEAMBbHilWnT5+uWCymO+64Q1dffbVWrlyp6667Tvv371d5eblTb1u6aPUNAPAgx+6I/PjHP9b777/v1F+PC9HqGwDgQew14xe0+gYAeBBBxC/c3OqbviYAgFEQRPwkGrWX6M6dm3k8HDa3dJe+JgCAMRR1r5lc0UckT6mUO1p909cEAEpSLtdvggickUrZdz5GW1IcCNh3ajo66GsCAD6Ty/WbqRk4I5e+JgCAkkUQgTPoawIAmACCCJxBXxMAwAQQROAM+poAACaAIAJnuLmvCQDANQgicI4b+5oAAFzFsb1mAEl22GhocEdfEwCA6xBE4LxgUKqrMz0KAIALMTUDAACMIYgAAABjCCIAAMAYakSQP7dsrgcA8CyCCPITi0lr12buJxMO271DWJYLAJggpmaQu1hMamwcualdd7d9PBYzMy4AgOcQRJCbVMq+E2JZI3+WPrZunX0eAADjIIggN+3tI++EDGdZUleXfR4AAOMgiCA3vb2FPQ8AUNIIIshNVVVhzwMAlDSCCHJTW2uvjrlwR920QECKROzzAAAYB0EEuQkG7SW60sgwkn69eTP9RAAAE0IQQe6iUamlRZo7N/N4OGwfp48IAGCCaGiG/ESjUkMDnVUBAJNCEEH+gkGprs70KAAAHsbUDAAAMIYgAgAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMKUoQSSaTuv766xUIBHTo0KFivCUAAPCAogSR3//+95ozZ04x3goAAHiI40Fkz5492rt3rzZt2uT0WwEAAI9xtLPqV199pUceeUTvvPOOpk+fPu75yWRSyWRy6HUikXByeAAAwDDH7ohYlqWHHnpIq1at0g033DChP9Pc3KxQKDT0iEQiTg0PAAC4QM5BZMOGDQoEAmM+jhw5oi1btujcuXNqamqa8N/d1NSkeDw+9Ojq6sp1eAAAwEMClmVZufyBvr4+nTlzZsxz5s+fr1/+8pf661//qkAgMHQ8lUopGAzq/vvv16uvvjrueyUSCYVCIcXjcVVUVOQyTAAAYEgu1++cg8hEdXZ2ZtR49PT06K677lJLS4tuuukmhcPhcf8OgggAAN6Ty/XbsWLV6urqjNczZsyQJC1YsGBCIQQAAPgfnVUBAIAxji7fHa6mpkYOzQIBAACP4o4IAAAwhiACAACMIYgAAABjCCIAAMCYohWrwgdSKam9XertlaqqpNpaKRg0PSoAgIcRRDAxsZi0dq108uT5Y+Gw9MILUjRqblwAAE9jagbji8WkxsbMECJJ3d328VjMzLgAAJ5HEMHYUin7Tki2HjDpY+vW2ecBAJAjggjG1t4+8k7IcJYldXXZ5wEAkCOCCMbW21vY8wAAGIYggrFVVRX2PAAAhiGIYGy1tfbqmEAg+88DASkSsc8DACBHBBGMLRi0l+hKI8NI+vXmzfQTAQDkhSCC8UWjUkuLNHdu5vFw2D5OHxEAQJ5oaIaJiUalhgY6qwIACooggokLBqW6OtOjAAD4CFMzAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBhXd1a1LEuSlEgkDI8EAABMVPq6nb6Oj8XVQeTcuXOSpEgkYngkAAAgV+fOnVMoFBrznIA1kbhiyODgoHp6ejRz5kwFLtyC3icSiYQikYi6urpUUVFhejhG8V1k4vs4j+8iE99HJr6P89zyXViWpXPnzmnOnDmaMmXsKhBX3xGZMmWKwuGw6WEURUVFRcn/DyiN7yIT38d5fBeZ+D4y8X2c54bvYrw7IWkUqwIAAGMIIgAAwBiCiGHl5eXauHGjysvLTQ/FOL6LTHwf5/FdZOL7yMT3cZ4XvwtXF6sCAAB/444IAAAwhiACAACMIYgAAABjCCIAAMAYgoiL3HPPPaqurtbUqVNVVVWlBx54QD09PaaHVXQnTpzQypUrNW/ePE2bNk0LFizQxo0bNTAwYHpoxjzzzDO65ZZbNH36dP3gBz8wPZyi27p1q2pqajR16lTddNNN+uSTT0wPyYgDBw7o7rvv1pw5cxQIBPTOO++YHpIxzc3NuvHGGzVz5kxdfvnlWrZsmY4ePWp6WMa8+OKLuu6664YamS1ZskR79uwxPawJIYi4SH19vd566y0dPXpUb7/9to4fP67GxkbTwyq6I0eOaHBwUNu2bdPhw4f1/PPP66WXXtITTzxhemjGDAwM6N5779Wjjz5qeihF9+abb2r9+vXauHGjPv/8cy1cuFB33XWXTp8+bXpoRdff36+FCxdq69atpodi3P79+7V69Wp99NFHeu+99/Tdd9/pzjvvVH9/v+mhGREOh/Xss8/q4MGD+uyzz3T77beroaFBhw8fNj208Vlwrd27d1uBQMAaGBgwPRTj/vjHP1rz5s0zPQzjtm/fboVCIdPDKKrFixdbq1evHnqdSqWsOXPmWM3NzQZHZZ4ka9euXaaH4RqnT5+2JFn79+83PRTXuPTSS62XX37Z9DDGxR0Rlzp79qxee+013XLLLbr44otND8e4eDyuWbNmmR4GimxgYEAHDx7U0qVLh45NmTJFS5cu1YcffmhwZHCbeDwuSfz/hKRUKqU33nhD/f39WrJkienhjIsg4jKPP/64LrnkEl122WXq7OzU7t27TQ/JuGPHjmnLli36zW9+Y3ooKLKvv/5aqVRKs2fPzjg+e/ZsnTp1ytCo4DaDg4Nat26dbr31Vl177bWmh2PMl19+qRkzZqi8vFyrVq3Srl27dM0115ge1rgIIg7bsGGDAoHAmI8jR44Mnf/YY4/piy++0N69exUMBvXggw/K8knz21y/C0nq7u7Wz3/+c91777165JFHDI3cGfl8HwBGWr16tf7973/rjTfeMD0Uo66++modOnRIH3/8sR599FGtWLFC//nPf0wPa1y0eHdYX1+fzpw5M+Y58+fPV1lZ2YjjJ0+eVCQS0QcffOCJ22vjyfW76OnpUV1dnW6++Wbt2LFDU6b4Kzfn89/Gjh07tG7dOn3zzTcOj84dBgYGNH36dLW0tGjZsmVDx1esWKFvvvmmpO8YBgIB7dq1K+N7KUVr1qzR7t27deDAAc2bN8/0cFxl6dKlWrBggbZt22Z6KGO6yPQA/K6yslKVlZV5/dnBwUFJUjKZLOSQjMnlu+ju7lZ9fb0WLVqk7du3+y6ESJP7b6NUlJWVadGiRdq3b9/QBXdwcFD79u3TmjVrzA4ORlmWpd/+9rfatWuX2traCCFZDA4OeuL6QRBxiY8//liffvqpbrvtNl166aU6fvy4nnzySS1YsMAXd0Ny0d3drbq6Ol155ZXatGmT+vr6hn52xRVXGByZOZ2dnTp79qw6OzuVSqV06NAhSdJVV12lGTNmmB2cw9avX68VK1bohhtu0OLFi7V582b19/fr4YcfNj20ovv222917NixodcdHR06dOiQZs2aperqaoMjK77Vq1dr586d2r17t2bOnDlUMxQKhTRt2jTDoyu+pqYm/eIXv1B1dbXOnTunnTt3qq2tTe+++67poY3P7KIdpP3rX/+y6uvrrVmzZlnl5eVWTU2NtWrVKuvkyZOmh1Z027dvtyRlfZSqFStWZP0+WltbTQ+tKLZs2WJVV1dbZWVl1uLFi62PPvrI9JCMaG1tzfrfwYoVK0wPrehG+/+I7du3mx6aEb/+9a+tK6+80iorK7MqKyutO+64w9q7d6/pYU0INSIAAMAY/028AwAAzyCIAAAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMOb/ARwQYacflUlzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i in range(40):\n",
    "    if real_class[i] == -1:\n",
    "        plt.scatter(X[i][0], X[i][1], c='b')\n",
    "    else:\n",
    "        plt.scatter(X[i][0], X[i][1], c='r')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 一、Fisher线性判别\n",
    "\n",
    "Fisher线性判别解决的是线性可分的问题。\n",
    "\n",
    "Fisher线性判别的核心思想是投影，使投影后的样本比较好的分开。\n",
    "\n",
    "比较好的分开是指：让投影后的每一类样本都尽可能聚集在这类的均值附近，并且不同类之间尽可能分开的远一点。\n",
    "\n",
    "即寻找一个$\\mathbf{w}$，使$y_i = \\mathbf{w}^\\text{T}\\mathbf{x}_i$，让$y_i$尽可能好的分开。\n",
    "\n",
    "#### 1. $\\mathbf{w}$的求法：\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\mathbf{m}_1 &= \\frac{1}{N_1} \\sum_{\\mathbf{x}_i\\in\\mathscr{X}_1} \\mathbf{x}_i \\\\\n",
    "\\mathbf{m}_2 &= \\frac{1}{N_2} \\sum_{\\mathbf{x}_i\\in\\mathscr{X}_2} \\mathbf{x}_i \\\\\n",
    "\\mathbf{S}_1 &= \\sum_{\\mathbf{x}_i\\in \\mathscr{X}_1}(\\mathbf{x}_i-\\mathbf{m}_1)(\\mathbf{x}_i-\\mathbf{m}_1)^\\text{T}\\\\\n",
    "\\mathbf{S}_2 &= \\sum_{\\mathbf{x}_i\\in \\mathscr{X}_2}(\\mathbf{x}_i-\\mathbf{m}_2)(\\mathbf{x}_i-\\mathbf{m}_2)^\\text{T}\\\\\n",
    "\\mathbf{S}_\\text{w} &= \\mathbf{S}_1 + \\mathbf{S}_2 \\\\\n",
    "\\mathbf{w} &= \\mathbf{S}_\\text{w}^{-1}(\\mathbf{m}_1 - \\mathbf{m}_2)\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "$\\mathbf{m}_1$和$\\mathbf{m}_2$是两类样本的均值，$\\mathbf{S}_\\text{w}$是pooled with-in class scatter matrix（总类内离散度矩阵，可以理解为协方差矩阵）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_mean(_X: np.ndarray) -> np.ndarray:\n",
    "    \"\"\"\n",
    "    计算m_i（样本均值）\n",
    "    :param _X: 样本矩阵，每一行是一个样本\n",
    "    :return: 均值（一维向量）\n",
    "    \"\"\"\n",
    "    return np.mean(_X, 0)\n",
    "\n",
    "\n",
    "def get_within_class_scatter_matrix(_X: np.ndarray) -> np.ndarray:\n",
    "    \"\"\"\n",
    "    计算S_i（类内离散度矩阵）\n",
    "    :param _X: 样本矩阵，每一行是一个样本\n",
    "    :return: 类内离散度矩阵（d * d）\n",
    "    \"\"\"\n",
    "    ret = np.zeros((_X.shape[1], _X.shape[1]))\n",
    "    m = get_mean(_X)\n",
    "    print(\"m_i = \")\n",
    "    print(m)\n",
    "    for row in _X:\n",
    "        ret += (row - m).reshape(m.shape[0], 1) @ (row - m).reshape(m.shape[0], 1).T\n",
    "    print(\"S_i = \")\n",
    "    print(ret)\n",
    "    return ret\n",
    "\n",
    "\n",
    "def get_pooled_within_class_scatter_matrix(_X1: np.ndarray, _X2: np.ndarray) -> np.ndarray:\n",
    "    \"\"\"\n",
    "    计算S_w（总类内离散度矩阵）\n",
    "    :param _X1: 第一类样本矩阵，每一行是一个样本\n",
    "    :param _X2: 第二类样本矩阵，每一行是一个样本\n",
    "    :return: 总类内离散度矩阵（d * d）\n",
    "    \"\"\"\n",
    "    return get_within_class_scatter_matrix(_X1) + get_within_class_scatter_matrix(_X2)\n",
    "\n",
    "\n",
    "def get_w(_X1: np.ndarray, _X2: np.ndarray) -> np.ndarray:\n",
    "    \"\"\"\n",
    "    计算最终需要的投影方向w\n",
    "    :param _X1: 第一类样本矩阵，每一行是一个样本\n",
    "    :param _X2: 第二类样本矩阵，每一行是一个样本\n",
    "    :return: 投影方向w（一维向量）\n",
    "    \"\"\"\n",
    "    return np.linalg.inv(get_pooled_within_class_scatter_matrix(_X1, _X2)) @ (get_mean(_X1) - get_mean(_X2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "m_i = \n",
      "[-1.99651597 -1.91231467]\n",
      "S_i = \n",
      "[[13.15645887  0.11877799]\n",
      " [ 0.11877799 22.70508313]]\n",
      "m_i = \n",
      "[1.81307045 1.82796564]\n",
      "S_i = \n",
      "[[14.72090202 -0.89478405]\n",
      " [-0.89478405 11.15114561]]\n",
      "w = \n",
      "[-0.13981965 -0.11368015]\n",
      "最终投影后的y = \n",
      "[ 0.31616146  0.39200903  0.4430584   0.48647637  0.43752815  0.41592967\n",
      "  0.50596349  0.72973533  0.23546217  0.30885178  0.62520699  0.72617969\n",
      "  0.82302367  0.15926573  0.5461533   0.62457474  0.49174992  0.50325523\n",
      "  0.67395599  0.48634664 -0.76955436 -0.49960364 -0.50229508 -0.51939192\n",
      " -0.57963085 -0.52343384 -0.22492428 -0.41965008 -0.19819595 -0.65763436\n",
      " -0.27527601 -0.39372772 -0.28272413 -0.50936426 -0.55295178 -0.31640878\n",
      " -0.55087554 -0.42813603 -0.48114958 -0.54119762]\n"
     ]
    }
   ],
   "source": [
    "w = get_w(X[:20], X[20:])\n",
    "print(\"w = \")\n",
    "print(w)\n",
    "y = X @ w\n",
    "print(\"最终投影后的y = \")\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面是投影后的图像，可以明显的看出，投影后两类样本已经很好的分开"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGdCAYAAADuR1K7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAnl0lEQVR4nO3df3RU5Z3H8c8kIUkVJik/kiEQRBD5oShtaGJoqR7JbhDOCkc8sMgisqzYFbBLPBxBo3Hr2rhKC7iAHN31sByhQaS4FFldDFpRYsSALb8VRPnlBClNgiJJSJ79IycjA5NkZpI7kzy8X+fc0+bOc5/7/c6dzP1wc2d0GWOMAAAALBET7QIAAADaEuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGCVuGgXEA319fU6efKkunTpIpfLFe1yAABAEIwxOnv2rNLS0hQT0/T1mSsy3Jw8eVLp6enRLgMAAITh2LFj6t27d5OPX5HhpkuXLpIanhy32x3lagAAQDCqqqqUnp7uO4835YoMN41/inK73YQbAAA6mJZuKeGGYgAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFglIuFm2bJl6tu3rxITE5WVlaWPPvqo2fHr1q3ToEGDlJiYqKFDh2rz5s1Njv3FL34hl8ulxYsXt3HVAACgI3I83Kxdu1Z5eXkqKCjQzp07dfPNNys3N1enTp0KOH779u2aPHmyZsyYoV27dmn8+PEaP3689uzZc9nYDRs26MMPP1RaWprTbQAAgA7C8XDz29/+Vvfff7+mT5+uIUOGaMWKFbrqqqv08ssvBxy/ZMkSjR49WvPmzdPgwYP11FNP6cc//rGWLl3qN+7EiROaM2eOVq9erU6dOjndBgAA6CAcDTc1NTUqKytTTk7O9zuMiVFOTo5KSkoCblNSUuI3XpJyc3P9xtfX12vq1KmaN2+ebrjhhhbrqK6uVlVVld8CAADs5Gi4OX36tOrq6pSamuq3PjU1VV6vN+A2Xq+3xfH//u//rri4OD300ENB1VFYWKikpCTfkp6eHmInAACgo+hwn5YqKyvTkiVLtHLlSrlcrqC2WbBggSorK33LsWPHHK4SAABEi6Phpnv37oqNjVV5ebnf+vLycnk8noDbeDyeZsdv27ZNp06dUp8+fRQXF6e4uDh9+eWXevjhh9W3b9+AcyYkJMjtdvstAADATo6Gm/j4eGVkZKi4uNi3rr6+XsXFxcrOzg64TXZ2tt94SdqyZYtv/NSpU/XnP/9Zn3zyiW9JS0vTvHnz9NZbbznXDAAA6BDinN5BXl6epk2bpuHDhyszM1OLFy/Wt99+q+nTp0uS7r33XvXq1UuFhYWSpF/+8pe69dZb9Zvf/EZjx45VUVGRPv74Y7344ouSpG7duqlbt25+++jUqZM8Ho8GDhzodDsAAKCdczzcTJo0SV9//bWeeOIJeb1eDRs2TG+++abvpuGjR48qJub7C0gjRozQmjVrlJ+fr0cffVQDBgzQ66+/rhtvvNHpUgEAgAVcxhgT7SIiraqqSklJSaqsrOT+GwAAOohgz98d7tNSAAAAzSHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsEpFws2zZMvXt21eJiYnKysrSRx991Oz4devWadCgQUpMTNTQoUO1efNm32O1tbV65JFHNHToUF199dVKS0vTvffeq5MnTzrdBgAA6AAcDzdr165VXl6eCgoKtHPnTt18883Kzc3VqVOnAo7fvn27Jk+erBkzZmjXrl0aP368xo8frz179kiSzp07p507d+rxxx/Xzp079fvf/14HDx7UnXfe6XQrAACgA3AZY4yTO8jKytJPfvITLV26VJJUX1+v9PR0zZkzR/Pnz79s/KRJk/Ttt99q06ZNvnW33HKLhg0bphUrVgTcx44dO5SZmakvv/xSffr0abGmqqoqJSUlqbKyUm63O8zOAABAJAV7/nb0yk1NTY3KysqUk5Pz/Q5jYpSTk6OSkpKA25SUlPiNl6Tc3Nwmx0tSZWWlXC6XkpOTAz5eXV2tqqoqvwUAANjJ0XBz+vRp1dXVKTU11W99amqqvF5vwG28Xm9I48+fP69HHnlEkydPbjLFFRYWKikpybekp6eH0Q0AAOgIOvSnpWprazVx4kQZY/TCCy80OW7BggWqrKz0LceOHYtglQAAIJLinJy8e/fuio2NVXl5ud/68vJyeTyegNt4PJ6gxjcGmy+//FJbt25t9m9vCQkJSkhICLMLAADQkTh65SY+Pl4ZGRkqLi72rauvr1dxcbGys7MDbpOdne03XpK2bNniN74x2Hz22Wd6++231a1bN2caAAAAHY6jV24kKS8vT9OmTdPw4cOVmZmpxYsX69tvv9X06dMlSffee6969eqlwsJCSdIvf/lL3XrrrfrNb36jsWPHqqioSB9//LFefPFFSQ3B5u6779bOnTu1adMm1dXV+e7H6dq1q+Lj451uCQAAtGOOh5tJkybp66+/1hNPPCGv16thw4bpzTff9N00fPToUcXEfH8BacSIEVqzZo3y8/P16KOPasCAAXr99dd14403SpJOnDihjRs3SpKGDRvmt6933nlHt912m9MtAQCAdszx77lpj/ieGwAAOp528T03AAAAkUa4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwSly0C7BKXZ20bZv01VdSz57SyJFSbGzH3ndL89bUSMuXS4cPS9deKw0dKp0+/f1YSXr33YZF+n57r7dhOXNGiomRbrutYZHC7+PiWlNSGtadOnX5PNF6ri4ec+KE9PXXUo8eUq9e4fcZqd4Qtmi+3MKdszUvz3DqbeqxcOtpfFv67DPJ5ZKysqS0NKm+XnrvvYYxjW85Tc3j1HG7+C2zf3/pwQel+PiWt2tNPXV1/m/DLfXeGu3mLchEwNKlS80111xjEhISTGZmpiktLW12/KuvvmoGDhxoEhISzI033mjeeOMNv8fr6+vN448/bjwej0lMTDSjRo0yn376adD1VFZWGkmmsrIyrH4CWr/emN69jZG+X3r3bljvNKf23dK88+YZExvr//jFS7duxnTu3PTjly6dOzdsE04fgWoNNE+0nquWamxNn073hrBF8+XWFnM6OXdLL9t588Krp6W3pUvfogLN49RxC1RbbGzD+nCfv5asX3/522pzvbdGJN6Cgj1/Ox5uioqKTHx8vHn55ZfN3r17zf3332+Sk5NNeXl5wPEffPCBiY2NNc8++6zZt2+fyc/PN506dTK7d+/2jXnmmWdMUlKSef31182f/vQnc+edd5prr73WfPfdd0HV1ObhZv16Y1yuy189LlfD4uTJxal9tzTvuHHBh5bWLMH00VStl87Tmn205rlqfAcPpsZw+nSyN4QtWr+a4czbFi/PcOpt7VvDpfXMmxfeXJf+G8SJ49ZSbU0FnNbUs359aL23RqROg+0m3GRmZppZs2b5fq6rqzNpaWmmsLAw4PiJEyeasWPH+q3LysoyDzzwgDGm4aqNx+Mxzz33nO/xiooKk5CQYH73u98FVVObhpsLF5q/YuByGZOe3jCurTm175bmjfTSXB9tVatTz5XL1fB4r17B1eFEn06+BhFQtH41w5k3lJeWk3OHs1xcT3V18FdsLl16926Yx6njFkxtsbEN40J5/trirbGx99aI5Gkw2PO3ozcU19TUqKysTDk5Ob51MTExysnJUUlJScBtSkpK/MZLUm5urm/8kSNH5PV6/cYkJSUpKyuryTmrq6tVVVXlt7SZbduk48ebftwY6dixhnFtzal9tzRvpDXXR1vV6tRzZUzD4ydOBDefE306+RpEQNH61Qxn3lBeWk7OHY6L61m+vOF+j3AcP94wj1PHLZja6uoaxl2sNfUE+9w39t4a0TwNNsXRcHP69GnV1dUpNTXVb31qaqq8Xm/Abbxeb7PjG/83lDkLCwuVlJTkW9LT08PqJ6Cvvmrbce1h307U2hYC1dXWtbaH58qpPtvrcbVQtH81Q5nXyZd8JF5yjfs4fLj18zh13IKt7dJxraknkscpmqfBplwRHwVfsGCBKisrfcuxY8fabvKePdt2XHvYtxO1toVAdbV1re3huXKqz/Z6XC0U7V/NUOZ18iUfiZdc4z7692/9PE4dt2Bru3Rca+qJ5HGK5mmwSa3/C1jTqqurTWxsrNmwYYPf+nvvvdfceeedAbdJT083ixYt8lv3xBNPmJtuuskYY8zhw4eNJLNr1y6/MT//+c/NQw89FFRdjtxz09TdcZG456at993SvJFegvnDclvcnejEc3XxPTfB1OhEn9xzE3HR+tVszT03rXl5tsXc4SxO3XPT1settffchFNPNO65icRpsF3ccxMfH6+MjAwVFxf71tXX16u4uFjZ2dkBt8nOzvYbL0lbtmzxjb/22mvl8Xj8xlRVVam0tLTJOR0VGystWdLw/10u/8caf1682JkP+ju175bmdbmkceNCLjcsLfXRXK2B5gk0zsnnSmp4/PnnW57L5QqvT6d6Q9ii9asZzrwXz9mc5l6ewcwdyss2GJfWEx8v5eWFPo/UUGNsrHPHLZja8vIu/76b1tQT7HFt7L01onkabFLrc1TzioqKTEJCglm5cqXZt2+fmTlzpklOTjZer9cYY8zUqVPN/PnzfeM/+OADExcXZxYuXGj2799vCgoKAn4UPDk52fzP//yP+fOf/2zGjRsX3Y+CGxP4A/7p6dH7npu22HdL80bie26C7aOl77lpnCdaz1VLNbamT6d7Q9ii+XJrizmdnLull21z33PTXD1Ofc9NWxy3tvyem1DeMqL5PTdt/RYU7PnbZYwxTgeopUuX6rnnnpPX69WwYcP0/PPPKysrS5J02223qW/fvlq5cqVv/Lp165Sfn68vvvhCAwYM0LPPPqsxY8ZcHMhUUFCgF198URUVFfrZz36m5cuX6/rrrw+qnqqqKiUlJamyslJut7vtGuUbivmGYr6hGJfgG4pbrpdvKG55O76huEGw5++IhJv2xrFwAwAAHBPs+fuK+LQUAAC4chBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWcSzcnDlzRlOmTJHb7VZycrJmzJihb775ptltzp8/r1mzZqlbt27q3LmzJkyYoPLyct/jf/rTnzR58mSlp6frBz/4gQYPHqwlS5Y41QIAAOiAHAs3U6ZM0d69e7VlyxZt2rRJ7733nmbOnNnsNnPnztUf/vAHrVu3Tn/84x918uRJ3XXXXb7Hy8rKlJKSoldeeUV79+7VY489pgULFmjp0qVOtQEAADoYlzHGtPWk+/fv15AhQ7Rjxw4NHz5ckvTmm29qzJgxOn78uNLS0i7bprKyUj169NCaNWt09913S5IOHDigwYMHq6SkRLfcckvAfc2aNUv79+/X1q1bg66vqqpKSUlJqqyslNvtDqNDAAAQacGevx25clNSUqLk5GRfsJGknJwcxcTEqLS0NOA2ZWVlqq2tVU5Ojm/doEGD1KdPH5WUlDS5r8rKSnXt2rXtigcAAB1anBOTer1epaSk+O8oLk5du3aV1+ttcpv4+HglJyf7rU9NTW1ym+3bt2vt2rV64403mq2nurpa1dXVvp+rqqqC6AIAAHREIV25mT9/vlwuV7PLgQMHnKrVz549ezRu3DgVFBTob//2b5sdW1hYqKSkJN+Snp4ekRoBAEDkhXTl5uGHH9Z9993X7Jh+/frJ4/Ho1KlTfusvXLigM2fOyOPxBNzO4/GopqZGFRUVfldvysvLL9tm3759GjVqlGbOnKn8/PwW616wYIHy8vJ8P1dVVRFwAACwVEjhpkePHurRo0eL47Kzs1VRUaGysjJlZGRIkrZu3ar6+nplZWUF3CYjI0OdOnVScXGxJkyYIEk6ePCgjh49quzsbN+4vXv36vbbb9e0adP09NNPB1V3QkKCEhISghoLAAA6Nkc+LSVJd9xxh8rLy7VixQrV1tZq+vTpGj58uNasWSNJOnHihEaNGqVVq1YpMzNTkvTP//zP2rx5s1auXCm32605c+ZIari3Rmr4U9Ttt9+u3NxcPffcc759xcbGBhW6GvFpKQAAOp5gz9+O3FAsSatXr9bs2bM1atQoxcTEaMKECXr++ed9j9fW1urgwYM6d+6cb92iRYt8Y6urq5Wbm6vly5f7Hn/ttdf09ddf65VXXtErr7ziW3/NNdfoiy++cKoVAADQgTh25aY948oNAAAdT1S/5wYAACBaCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFUcCzdnzpzRlClT5Ha7lZycrBkzZuibb75pdpvz589r1qxZ6tatmzp37qwJEyaovLw84Ni//OUv6t27t1wulyoqKhzoAAAAdESOhZspU6Zo79692rJlizZt2qT33ntPM2fObHabuXPn6g9/+IPWrVunP/7xjzp58qTuuuuugGNnzJihm266yYnSAQBAB+Yyxpi2nnT//v0aMmSIduzYoeHDh0uS3nzzTY0ZM0bHjx9XWlraZdtUVlaqR48eWrNmje6++25J0oEDBzR48GCVlJTolltu8Y194YUXtHbtWj3xxBMaNWqU/vrXvyo5OTno+qqqqpSUlKTKykq53e7WNQsAACIi2PO3I1duSkpKlJyc7As2kpSTk6OYmBiVlpYG3KasrEy1tbXKycnxrRs0aJD69OmjkpIS37p9+/bpV7/6lVatWqWYmODKr66uVlVVld8CAADs5Ei48Xq9SklJ8VsXFxenrl27yuv1NrlNfHz8ZVdgUlNTfdtUV1dr8uTJeu6559SnT5+g6yksLFRSUpJvSU9PD60hAADQYYQUbubPny+Xy9XscuDAAadq1YIFCzR48GD9wz/8Q8jbVVZW+pZjx445VCEAAIi2uFAGP/zww7rvvvuaHdOvXz95PB6dOnXKb/2FCxd05swZeTyegNt5PB7V1NSooqLC7+pNeXm5b5utW7dq9+7deu211yRJjbcLde/eXY899pj+9V//NeDcCQkJSkhICKZFAADQwYUUbnr06KEePXq0OC47O1sVFRUqKytTRkaGpIZgUl9fr6ysrIDbZGRkqFOnTiouLtaECRMkSQcPHtTRo0eVnZ0tSVq/fr2+++473zY7duzQP/7jP2rbtm3q379/KK0AAABLhRRugjV48GCNHj1a999/v1asWKHa2lrNnj1bf//3f+/7pNSJEyc0atQorVq1SpmZmUpKStKMGTOUl5enrl27yu12a86cOcrOzvZ9UurSAHP69Gnf/kL5tBQAALCXI+FGklavXq3Zs2dr1KhRiomJ0YQJE/T888/7Hq+trdXBgwd17tw537pFixb5xlZXVys3N1fLly93qkQAAGAhR77npr3je24AAOh4ovo9NwAAANFCuAEAAFYh3AAAAKsQbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWiYt2AdFgjJEkVVVVRbkSAAAQrMbzduN5vClXZLg5e/asJCk9PT3KlQAAgFCdPXtWSUlJTT7uMi3FHwvV19fr5MmT6tKli1wuV9DbVVVVKT09XceOHZPb7XawwuixvUfb+5Ps79H2/iR6tIHt/UnR6dEYo7NnzyotLU0xMU3fWXNFXrmJiYlR7969w97e7XZb+2JtZHuPtvcn2d+j7f1J9GgD2/uTIt9jc1dsGnFDMQAAsArhBgAAWIVwE4KEhAQVFBQoISEh2qU4xvYebe9Psr9H2/uT6NEGtvcnte8er8gbigEAgL24cgMAAKxCuAEAAFYh3AAAAKsQbgAAgFUINy04c+aMpkyZIrfbreTkZM2YMUPffPNNs9t4vV5NnTpVHo9HV199tX784x9r/fr1Eao4NOH0J0klJSW6/fbbdfXVV8vtduvnP/+5vvvuuwhUHLpwe5Qavg3zjjvukMvl0uuvv+5soWEKtb8zZ85ozpw5GjhwoH7wgx+oT58+euihh1RZWRnBqpu3bNky9e3bV4mJicrKytJHH33U7Ph169Zp0KBBSkxM1NChQ7V58+YIVRq+UHp86aWXNHLkSP3whz/UD3/4Q+Xk5LT4nLQHoR7HRkVFRXK5XBo/fryzBbZSqP1VVFRo1qxZ6tmzpxISEnT99de3+9dqqD0uXrzY996Snp6uuXPn6vz58xGq9iIGzRo9erS5+eabzYcffmi2bdtmrrvuOjN58uRmt/mbv/kb85Of/MSUlpaaw4cPm6eeesrExMSYnTt3Rqjq4IXT3/bt243b7TaFhYVmz5495sCBA2bt2rXm/PnzEao6NOH02Oi3v/2tueOOO4wks2HDBmcLDVOo/e3evdvcddddZuPGjebQoUOmuLjYDBgwwEyYMCGCVTetqKjIxMfHm5dfftns3bvX3H///SY5OdmUl5cHHP/BBx+Y2NhY8+yzz5p9+/aZ/Px806lTJ7N79+4IVx68UHu85557zLJly8yuXbvM/v37zX333WeSkpLM8ePHI1x58ELtsdGRI0dMr169zMiRI824ceMiU2wYQu2vurraDB8+3IwZM8a8//775siRI+bdd981n3zySYQrD16oPa5evdokJCSY1atXmyNHjpi33nrL9OzZ08ydOzfClRtDuGnGvn37jCSzY8cO37r//d//NS6Xy5w4caLJ7a6++mqzatUqv3Vdu3Y1L730kmO1hiPc/rKyskx+fn4kSmy1cHs0xphdu3aZXr16ma+++qrdhpvW9HexV1991cTHx5va2lonygxJZmammTVrlu/nuro6k5aWZgoLCwOOnzhxohk7dqzfuqysLPPAAw84WmdrhNrjpS5cuGC6dOli/vu//9upElstnB4vXLhgRowYYf7zP//TTJs2rV2Hm1D7e+GFF0y/fv1MTU1NpEpstVB7nDVrlrn99tv91uXl5Zmf/vSnjtYZCH+WakZJSYmSk5M1fPhw37qcnBzFxMSotLS0ye1GjBihtWvX6syZM6qvr1dRUZHOnz+v2267LQJVBy+c/k6dOqXS0lKlpKRoxIgRSk1N1a233qr3338/UmWHJNxjeO7cOd1zzz1atmyZPB5PJEoNS7j9XaqyslJut1txcdH9z83V1NSorKxMOTk5vnUxMTHKyclRSUlJwG1KSkr8xktSbm5uk+OjLZweL3Xu3DnV1taqa9euTpXZKuH2+Ktf/UopKSmaMWNGJMoMWzj9bdy4UdnZ2Zo1a5ZSU1N144036te//rXq6uoiVXZIwulxxIgRKisr8/3p6vPPP9fmzZs1ZsyYiNR8sSvyP5wZLK/Xq5SUFL91cXFx6tq1q7xeb5Pbvfrqq5o0aZK6deumuLg4XXXVVdqwYYOuu+46p0sOSTj9ff7555KkJ598UgsXLtSwYcO0atUqjRo1Snv27NGAAQMcrzsU4R7DuXPnasSIERo3bpzTJbZKuP1d7PTp03rqqac0c+ZMJ0oMyenTp1VXV6fU1FS/9ampqTpw4EDAbbxeb8DxwfYfaeH0eKlHHnlEaWlpl4W69iKcHt9//33913/9lz755JMIVNg64fT3+eefa+vWrZoyZYo2b96sQ4cO6cEHH1Rtba0KCgoiUXZIwunxnnvu0enTp/Wzn/1MxhhduHBBv/jFL/Too49GomQ/V+SVm/nz58vlcjW7BPsmE8jjjz+uiooKvf322/r444+Vl5eniRMnavfu3W3YRdOc7K++vl6S9MADD2j69On60Y9+pEWLFmngwIF6+eWX27KNZjnZ48aNG7V161YtXry4bYsOgdOv0UZVVVUaO3ashgwZoieffLL1hcNxzzzzjIqKirRhwwYlJiZGu5w2cfbsWU2dOlUvvfSSunfvHu1yHFFfX6+UlBS9+OKLysjI0KRJk/TYY49pxYoV0S6tzbz77rv69a9/reXLl2vnzp36/e9/rzfeeENPPfVUxGu5Iq/cPPzww7rvvvuaHdOvXz95PB6dOnXKb/2FCxd05syZJv9UcfjwYS1dulR79uzRDTfcIEm6+eabtW3bNi1btiwiL2Qn++vZs6ckaciQIX7rBw8erKNHj4ZfdIic7HHr1q06fPiwkpOT/dZPmDBBI0eO1LvvvtuKyoPjZH+Nzp49q9GjR6tLly7asGGDOnXq1NqyW6179+6KjY1VeXm53/ry8vIm+/F4PCGNj7Zwemy0cOFCPfPMM3r77bd10003OVlmq4Ta4+HDh/XFF1/o7/7u73zrGv8hFRcXp4MHD6p///7OFh2CcI5hz5491alTJ8XGxvrWDR48WF6vVzU1NYqPj3e05lCF0+Pjjz+uqVOn6p/+6Z8kSUOHDtW3336rmTNn6rHHHlNMTOSup1yR4aZHjx7q0aNHi+Oys7NVUVGhsrIyZWRkSGo48dXX1ysrKyvgNufOnZOkyw5ibGys75fVaU7217dvX6WlpengwYN+6z/99FPdcccdrS8+SE72OH/+fN8vZ6OhQ4dq0aJFfm++TnKyP6nhik1ubq4SEhK0cePGdnMFID4+XhkZGSouLvZ9DLi+vl7FxcWaPXt2wG2ys7NVXFysf/mXf/Gt27Jli7KzsyNQcejC6VGSnn32WT399NN66623/O6xao9C7XHQoEGXXdnOz8/X2bNntWTJEqWnp0ei7KCFcwx/+tOfas2aNaqvr/edHz799FP17Nmz3QUbKbwez507F/DcJzV8rUZERfwW5g5m9OjR5kc/+pEpLS0177//vhkwYIDfx2yPHz9uBg4caEpLS40xxtTU1JjrrrvOjBw50pSWlppDhw6ZhQsXGpfLZd54441otdGkUPszxphFixYZt9tt1q1bZz777DOTn59vEhMTzaFDh6LRQovC6fFSaqefljIm9P4qKytNVlaWGTp0qDl06JD56quvfMuFCxei1YZPUVGRSUhIMCtXrjT79u0zM2fONMnJycbr9RpjjJk6daqZP3++b/wHH3xg4uLizMKFC83+/ftNQUFBh/goeCg9PvPMMyY+Pt689tprfsfr7Nmz0WqhRaH2eKn2/mmpUPs7evSo6dKli5k9e7Y5ePCg2bRpk0lJSTH/9m//Fq0WWhRqjwUFBaZLly7md7/7nfn888/N//3f/5n+/fubiRMnRrx2wk0L/vKXv5jJkyebzp07G7fbbaZPn+73hnLkyBEjybzzzju+dZ9++qm56667TEpKirnqqqvMTTfddNlHw9uLcPozxpjCwkLTu3dvc9VVV5ns7Gyzbdu2CFcevHB7vFh7Djeh9vfOO+8YSQGXI0eORKeJS/zHf/yH6dOnj4mPjzeZmZnmww8/9D126623mmnTpvmNf/XVV831119v4uPjzQ033NAu/yFxqVB6vOaaawIer4KCgsgXHoJQj+PF2nu4MSb0/rZv326ysrJMQkKC6devn3n66afbxT8omhNKj7W1tebJJ580/fv3N4mJiSY9Pd08+OCD5q9//WvE63YZE+lrRQAAAM65Ij8tBQAA7EW4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBV/h/HypPZOkipKgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.clf()\n",
    "for y_i in y[:20]:\n",
    "    plt.scatter(y_i, 0, c='b')\n",
    "for y_i in y[20:]:\n",
    "    plt.scatter(y_i, 0, c='r')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2. 分类\n",
    "\n",
    "假设样本服从协方差矩阵相同的正态分布，可以得到分类规则：\n",
    "\n",
    "$$若 g(\\mathbf{x}) = \\mathbf{w}^\\text{T}(\\mathbf{x} - \\frac{1}{2}(\\mathbf{m}_1 + \\mathbf{m}_2)) \\gtrless \\ln\\frac{P(\\omega_2)}{P(\\omega_1)}， 则 \\mathbf{x}\\in\\begin{cases}\\omega_1\\\\\\omega_2\\end{cases}$$\n",
    "\n",
    "先验概率可以通过样本类别的数量进行估算，因为上面的样本两类数量一样，所以这里$P(\\omega_1)=P(\\omega_2)$\n",
    "\n",
    "如果样本不符合正态分布，这样的分类阈值不一定为最优解，但是通常也可以取得较好效果（也就是说分界点歪了点，但是也能用）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def classify_Fisher(_X: np.ndarray, _w: np.ndarray, _m_1: np.ndarray, _m_2: np.ndarray) -> np.ndarray:\n",
    "    \"\"\"\n",
    "    Fisher判别器分类\n",
    "    :param _X: 样本矩阵（N * d）\n",
    "    :param _w: 上面计算出来的w（投影方向）\n",
    "    :param _m_1: 第一类样本均值\n",
    "    :param _m_2: 第二类样本均值\n",
    "    :return: 分类结果（一维向量）\n",
    "    \"\"\"\n",
    "    g = ((_w @ (_X - 1 / 2 * (_m_1 + _m_2)).T) >= 0.0).astype(int)\n",
    "    for i, g_i in enumerate(g):\n",
    "        if g_i == 0:\n",
    "            g[i] = -1\n",
    "    return g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ans_Fisher = \n",
      "[ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1 -1 -1 -1 -1\n",
      " -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1]\n",
      "正确率： 1.0\n"
     ]
    }
   ],
   "source": [
    "ans_Fisher = classify_Fisher(X, w, get_mean(X[:20]), get_mean(X[20:]))\n",
    "print(\"ans_Fisher = \")\n",
    "print(ans_Fisher)\n",
    "print(\"正确率：\", 1 - np.sum(abs(ans_Fisher - real_class)) / 2 / X.shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 二、感知器\n",
    "\n",
    "感知器的思想是通过迭代梯度下降，找到合适的线性判别函数。如果样本是线性可分的，梯度下降最终就会收敛。\n",
    "\n",
    "感知机解决的也是线性可分的情况。\n",
    "\n",
    "#### 1. 预处理数据\n",
    "\n",
    "1. 将$\\mathbf{x}$增广为$\\mathbf{y} = [1, \\mathbf{x^\\text{T}}]^\\text{T}$\n",
    "2. 再次处理$\\mathbf{y}=\\begin{cases} \\mathbf{y}, &若\\mathbf{y}\\in \\omega_1 \\\\ -\\mathbf{y}, &若\\mathbf{y}\\in \\omega_2\\end{cases}$\n",
    "3. 设权向量$\\boldsymbol{\\alpha}$是$d+1$维向量\n",
    "4. 线性判别函数变为$g(\\mathbf{y}) = \\boldsymbol{\\alpha}^\\text{T}\\mathbf{y}$\n",
    "5. 可以设置余量$b$来避免噪声误差，也可以使$b=0$\n",
    "\n",
    "预处理后感知器的核心问题就是求出权向量$\\boldsymbol{\\alpha}$的值，因为求出了权向量，就相当于找到了线性判别函数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.         -0.99720813 -1.55464404]\n",
      " [ 1.         -2.15332905 -0.79988732]\n",
      " [ 1.         -1.6271845  -1.89607447]\n",
      " [ 1.         -2.88804683 -0.72722164]\n",
      " [ 1.         -0.79871052 -2.86639942]\n",
      " [ 1.         -1.9454445  -1.2659932 ]\n",
      " [ 1.         -2.63452428 -1.21045949]\n",
      " [ 1.         -2.28814449 -3.60491927]\n",
      " [ 1.         -0.14831176 -1.88885453]\n",
      " [ 1.         -1.31782927 -1.09599919]\n",
      " [ 1.         -3.50049996 -1.19430091]\n",
      " [ 1.         -2.6730074  -3.10028377]\n",
      " [ 1.         -2.34421754 -4.356574  ]\n",
      " [ 1.         -1.17966467  0.04991699]\n",
      " [ 1.         -2.38234616 -1.87415735]\n",
      " [ 1.         -3.07994777 -1.70599279]\n",
      " [ 1.         -1.41432354 -2.58620074]\n",
      " [ 1.         -1.96619849 -2.00863597]\n",
      " [ 1.         -2.14275402 -3.29307151]\n",
      " [ 1.         -2.44862645 -1.26654069]\n",
      " [-1.         -3.17028293 -2.87021525]\n",
      " [-1.         -1.72407696 -2.27430904]\n",
      " [-1.         -1.94459064 -2.02676625]\n",
      " [-1.         -2.95914397 -0.92932175]\n",
      " [-1.         -2.70985559 -1.76582962]\n",
      " [-1.         -2.71996217 -1.25905597]\n",
      " [-1.         -0.40617109 -1.47900559]\n",
      " [-1.         -1.718164   -1.57826132]\n",
      " [-1.         -0.52405299 -1.09889933]\n",
      " [-1.         -2.56521233 -2.62989855]\n",
      " [-1.         -1.1198243  -1.04418024]\n",
      " [-1.         -1.86167209 -1.17372623]\n",
      " [-1.         -0.09299459 -2.37263639]\n",
      " [-1.         -1.42149269 -2.73232963]\n",
      " [-1.         -1.44622845 -3.0853286 ]\n",
      " [-1.         -1.63555418 -0.77169289]\n",
      " [-1.         -1.62201784 -2.85085454]\n",
      " [-1.         -1.97566001 -1.33620458]\n",
      " [-1.         -1.4990023  -2.388804  ]\n",
      " [-1.         -3.14544988 -0.89199308]]\n"
     ]
    }
   ],
   "source": [
    "Y = np.concatenate((np.ones((X.shape[0], 1)), X), 1)\n",
    "\n",
    "for i in range(20, 40):\n",
    "    Y[i] = -Y[i]\n",
    "\n",
    "print(Y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2. 权向量$\\boldsymbol{\\alpha}$的计算\n",
    "\n",
    "采用迭代法，每次考察一个样本$\\mathbf{y}_i$，判断$\\boldsymbol{\\alpha}^\\text{T}\\mathbf{y}_i \\leq b$是否成立。\n",
    "\n",
    "如果不成立，令$\\boldsymbol{\\alpha} = \\boldsymbol{\\alpha} + \\rho\\mathbf{y}_i$\n",
    "\n",
    "其中，$\\rho$可以设置为$1$，也可以按照**绝对修正法**设置为\n",
    "\n",
    "$$\\rho = \\frac{|\\boldsymbol{\\alpha^\\text{T}}\\mathbf{y}_i|}{||\\mathbf{y}_i||^2}$$\n",
    "\n",
    "这样收敛更快。\n",
    "\n",
    "接着考察下一个样本，如果当前一轮所有样本都考察过并且都满足条件，则判断为收敛。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_alpha(_Y: np.ndarray, _b: float, _is_rho_enable: bool) -> (np.ndarray, int):\n",
    "    \"\"\"\n",
    "    迭代计算最终的alpha值\n",
    "    :param _Y: 预处理后的样本矩阵（N * d+1）\n",
    "    :param _b: 余量b\n",
    "    :param _is_rho_enable: False: rho设置为1, True: rho设置为绝对修正法\n",
    "    :return: 最终的alpha, 迭代次数\n",
    "    \"\"\"\n",
    "    i, n = 0, 0\n",
    "    iters = 0\n",
    "    alpha = np.zeros(_Y.shape[1])  # 初始化 alpha = [0. 0. 0.]\n",
    "    while n < _Y.shape[0]:\n",
    "        if np.sum(alpha * _Y[i]) <= _b:\n",
    "            rho = 1\n",
    "            if _is_rho_enable:\n",
    "                rho = abs(np.sum(alpha * _Y[i])) / np.sum(_Y[i] * _Y[i])\n",
    "                if rho == 0:  # 这里要特别注意，书上没说。步长显然是一个正数，不能为0，如果是0会陷入死循环。\n",
    "                    rho = 1\n",
    "            n = -1\n",
    "            alpha = alpha + rho * _Y[i]\n",
    "            iters += 1\n",
    "        i += 1\n",
    "        if i >= _Y.shape[0]:\n",
    "            i = 0\n",
    "        n += 1\n",
    "        if iters > 100000:\n",
    "            raise Exception('迭代100000次还不收敛，数据可能不是线性可分的')\n",
    "    return alpha, iters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "将rho设置为1：\n",
      "alpha = [ 1.         -0.99720813 -1.55464404] \n",
      "迭代次数 = 1\n",
      "将rho采用绝对修正法：\n",
      "alpha = [ 1.         -0.99720813 -1.55464404] \n",
      "迭代次数 = 1\n"
     ]
    }
   ],
   "source": [
    "alpha_without_rho, iters_without_rho = get_alpha(Y, 0, False)\n",
    "print(\"将rho设置为1：\\nalpha =\", alpha_without_rho, \"\\n迭代次数 =\", iters_without_rho)\n",
    "alpha_with_rho, iters_with_rho = get_alpha(Y, 0, True)\n",
    "print(\"将rho采用绝对修正法：\\nalpha =\", alpha_without_rho, \"\\n迭代次数 =\", iters_without_rho)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "因为数据维度太小了，收敛的太快了，一次迭代就收敛了，所以没显出绝对修正法的速度优势"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3. 分类\n",
    "\n",
    "判别规则如下：\n",
    "\n",
    "$$若g(\\mathbf{y}) = \\boldsymbol{\\alpha}^\\text{T}\\mathbf{y} \\gtrless 0, 则\\mathbf{y}\\in\\begin{cases}\\omega_1\\\\ \\omega_2\\end{cases}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def classify_perceptron(_Y: np.ndarray, _alpha: np.ndarray) -> np.ndarray:\n",
    "    \"\"\"\n",
    "    感知机方法分类\n",
    "    :param _Y: 增广后的样本矩阵（每一行为一个样本），*****注意这里只增广，不改变正负号*****\n",
    "    :param _alpha: 前面算出的alpha值\n",
    "    :return: 返回分类好的矩阵（一维向量）\n",
    "    \"\"\"\n",
    "    g = ((_alpha @ _Y.T) >= 0).astype(int)\n",
    "    for i, g_i in enumerate(g):\n",
    "        if g_i == 0:\n",
    "            g[i] = -1\n",
    "    return g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ans_perceptron = \n",
      "[ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1 -1 -1 -1 -1\n",
      " -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1]\n",
      "正确率： 1.0\n"
     ]
    }
   ],
   "source": [
    "# 注意我们前面把Y预处理过，把原来的第二类样本y变成了-y，这里要变回来才能分类\n",
    "ans_perceptron = classify_perceptron(np.concatenate([Y[:20], -Y[20:]]), alpha_without_rho)\n",
    "print(\"ans_perceptron = \")\n",
    "print(ans_perceptron)\n",
    "print(\"正确率：\", 1 - np.sum(abs(ans_perceptron - real_class)) / 2 / X.shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 三、最小平方误差判别\n",
    "\n",
    "这个算法用于解决线性不可分，解决的思想是：引入松弛变量，使投影后的样本线性可分。并且让松弛变量尽可能小。具体使松弛变量变小的方法是用最小平方误差使投影后的样本离人为设定的阈值最近。\n",
    "\n",
    "先生成一下新的线性不可分的数据。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-1.34299662 -1.04252449]\n",
      " [-2.79854016 -2.59819657]\n",
      " [-0.77424817 -2.18673294]\n",
      " [-1.34896419  0.1280363 ]\n",
      " [-2.39101614 -1.86949986]\n",
      " [-1.05532767 -0.6941443 ]\n",
      " [-0.01296092 -0.14118512]\n",
      " [ 0.18570744 -0.42184708]\n",
      " [-1.17572825 -0.48215859]\n",
      " [-0.6787651  -0.74512697]\n",
      " [ 0.07235573  0.22487309]\n",
      " [-1.45626639 -1.29698314]\n",
      " [-3.03224089 -0.73369765]\n",
      " [-1.63819157  0.04351657]\n",
      " [-0.75748584 -2.12514056]\n",
      " [-0.62754058 -1.70106403]\n",
      " [-0.14226139 -0.44542636]\n",
      " [-0.66372502 -1.22693198]\n",
      " [-2.32776608 -0.59538424]\n",
      " [-1.17255684 -1.10848414]\n",
      " [ 1.01004954  1.96209479]\n",
      " [-0.46661866 -0.02766581]\n",
      " [ 0.85500022  1.15566666]\n",
      " [ 1.47413048  1.28779157]\n",
      " [ 0.4647262   2.94319899]\n",
      " [ 0.79042834  1.69130797]\n",
      " [ 1.91553629  0.52362328]\n",
      " [ 0.84999404  1.36781803]\n",
      " [ 1.11544251  0.23309943]\n",
      " [ 0.46791749  1.96889589]\n",
      " [ 1.17958567 -0.29744786]\n",
      " [ 0.57513099  1.41559357]\n",
      " [ 1.53657285  1.13196758]\n",
      " [ 2.1828135  -0.06019617]\n",
      " [ 0.59310113  0.97142505]\n",
      " [ 0.30757167  0.34764882]\n",
      " [ 1.67442632  0.9752627 ]\n",
      " [ 0.03075503  1.00573105]\n",
      " [ 1.45106052  1.98296215]\n",
      " [ 0.80973014  0.20798233]]\n",
      "[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1]\n"
     ]
    }
   ],
   "source": [
    "X = np.r_[np.random.randn(20, 2) - [1, 1], np.random.randn(20, 2) + [1, 1]]\n",
    "real_class = 20 * [1] + 20 * [-1]\n",
    "\n",
    "print(X)\n",
    "print(real_class)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "再生成一下散点图，发现确实是线性不可分的情况（如果你发现是线性可分的，就重新跑几次上面生成数据的代码）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGdCAYAAAAvwBgXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAivElEQVR4nO3da2xb9f3H8c+pWUMZiVnXC2ntklI2aWiCSQW6dgtKpg7YA0hlpZqoNApDFaCCGnWa1k7aqj1AmQTaUjHERdvaPVgLW+bChLQxVCUlk+jY2KoNpFb/QqqkbloCCLvkgYNO/H9wcFonTuLbOb9zeb8ky+T4l/jX05Lz8e/yPVahUCgIAADAgEWmOwAAAKKLIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAmCtMd2A+U1NTOnfunJqbm2VZlunuAACAChQKBV28eFGrVq3SokXzj3n4OoicO3dOyWTSdDcAAEANRkdHlUgk5m3j6yDS3NwsyfmDtLS0GO4NAACoRC6XUzKZnL6Oz8fXQaQ4HdPS0kIQAQAgYCpZVsFiVQAAYAxBBAAAGEMQAQAAxhBEAACAMQQRAABgDEEEAAAY42oQeeaZZ3TTTTdNb7/duHGj/vKXv7j5lgAAIEBcDSKJREI///nP9dZbb+lf//qXvvWtb6mrq0vvvPOOm28LAAACwioUCgUv33Dp0qV64okn9OCDDy7YNpfLKR6PK5vNUtAMQKTZtjQ0JI2NSa2tUnu7FIuZ7hVQXjXXb88qq9q2rT/+8Y+amJjQxo0by7bJ5/PK5/PTX+dyOa+6BwC+lU5Lu3ZJZ89eOpZISPv3S6mUuX4BjeD6YtX//e9/uvrqq9XU1KSHH35YR44c0Y033li2bW9vr+Lx+PSDG94BiLp0WuruLg0hkpTJOMfTaTP9AhrF9amZyclJjYyMKJvNqr+/X7/+9a917NixsmGk3IhIMplkagZAJNm21NY2O4QUWZYzMjI8zDQN/KWaqRnP14hs3rxZ69at03PPPbdgW9aIAIiywUGps3PhdgMDUkeH270BKlfN9dvzOiJTU1Mlox4AgPLGxhrbDvAjVxer7t27V9/5zne0Zs0aXbx4UYcOHdLg4KBeffVVN98WAEKhtbWx7QA/cjWIvP/++7rvvvs0NjameDyum266Sa+++qq+/e1vu/m2ABAK7e3OGpBMRio3iV5cI9Le7n3fgEZxNYj85je/cfPHA0CoxWLOFt3ubid0XB5GLMt57utjoSqCjXvNAICPpVJSf7+0enXp8UTCOU4dEQSdZwXNAAC1SaWkri4qqyKcCCIAEACxGFt0EU5MzQAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYV4NIb2+vbr31VjU3N2vFihXasmWLTp065eZbAgCAAHE1iBw7dkw7d+7U8ePH9dprr+nTTz/VHXfcoYmJCTffFgCABdm2NDgoHT7sPNu26R5Fk1UoFApevdn4+LhWrFihY8eO6fbbb1+wfS6XUzweVzabVUtLiwc9BAB/sm1paEgaG5NaW6X2dikWM92ryvix7+m0tGuXdPbspWOJhLR/v5RKmetXWFRz/fZ0jUg2m5UkLV261Mu3BYBAS6eltjaps1Pats15bmtzjvudH/ueTkvd3aUhRJIyGed4EM5rmHg2IjI1NaV77rlHH3/8sf7+97+XbZPP55XP56e/zuVySiaTjIgAiKziRXPmb2rLcp77+/37Cd6PfbdtJwjNDCFFluWMjAwPmx+1CTJfjojs3LlTb7/9tl544YU52/T29ioej08/ksmkV90DAN+xbWf6oNzHxeKxnh5/rm3wa9+HhuYOIZLTt9FRpx284UkQefTRR/XKK69oYGBAiURiznZ79+5VNpudfoyOjnrRPQDwpSBfNP3a97GxxrZD/a5w84cXCgU99thjOnLkiAYHB7V27dp52zc1NampqcnNLgFAYAT5ounXvre2NrYd6udqENm5c6cOHTqkl19+Wc3NzTp//rwkKR6Pa8mSJW6+NQAEXpAvmn7te3u7swYkkyk/bVRcI9Le7m2/oszVxapWcUXSDAcOHND999+/4PezfRdAlBUXVi500fTjwko/9724iFYq7VsQFgAHhW8WqxYKhbKPSkIIAERdLObUtZAuXSSLil/39fkvhEj+7nsq5YSN1atLjycShBATuNcMAPhYkC+afu57KiWdOSMNDEiHDjnPw8P+Pp9h5Wll1WoxNQMADj9WJ61UkPuO2lRz/XZ1sSoAoDFiMamjw3QvahPkvsN9TM0AAABjCCIAAMAYgggAADCGIAIAAIxhsSoAYBZ2usArBBEAQIl02rlz7uU3rUsknAJl1NlAozE1AwCYVix/PvPOuZmMczydNtMvhBdBBAAgyZmO2bWr/L1hisd6epx2QKMQRAAAkpw1ITNHQi5XKEijo047oFEIIgAASc7C1Ea2AypBEAEASHJ2xzSyHVAJgggAQJKzRTeRkCyr/OuWJSWTTjugUQgiAABJTp2Q/fud/54ZRopf9/VRTwSNRRABAExLpaT+fmn16tLjiYRznDoiaDQKmgEASqRSUlcXlVXhDYIIAGCWWEzq6DDdC0QBUzMAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBh2zQBAiNg2224RLAQRAAiJdFratav0DrqJhFMtlUJk8CumZgAgBNJpqbu7NIRIUibjHE+nzfQLWAhBBAACzradkZBCYfZrxWM9PU67et9ncFA6fNh5rvfnARJBBAACb2ho9kjI5QoFaXTUaVerdFpqa5M6O6Vt25zntjZGWlA/gggABNzYWGPbzcS0D9xEEAGAgGttbWy7y3k17YPoIogAQMC1tzu7Yyyr/OuWJSWTTrtqeTHtg2gjiABAwMVizhZdaXYYKX7d11dbPRG3p338gEW4ZhFEACAEUimpv19avbr0eCLhHK+1joib0z5+wCJc86xCodzMnz/kcjnF43Fls1m1tLSY7g4A+F6jK6vatnNhzmTKrxOxLCfsDA8Hr4JrcRHuzD9XcRSpngAXddVcvwkiAIB5FS/YUulFO8gX7GLAmmv9S5ADlh9Uc/1magYAMC+3pn1MYhGuf3CvGQDAglIpqasrPDfUi8Ii3KAgiAAAKhKLSR0dpnvRGGFfhBskTM0AACLHzdorqA5BBABQsbDU3HCz9gqqQxABAFQkbDU3FlqE29UVjtDldwQRAMCCwnrju1RKOnNGGhiQDh1ynoeHndfCFLr8jDoiAIB5Ra3mBoXO6kcdEQBAw0Sp5gZ3G/YeQQQAMK8o1dyIUujyC4IIAGBeUaq5EaXQ5RcEEQDAvKJUcyNKocsvCCIAgHlFqeZGlEKXXxBEAGCGsBTtaqQw3viunCiFLr9wNYi8/vrruvvuu7Vq1SpZlqWXXnrJzbcDgLqFrWhXI81VcyMsIaQoKqHLL1y96d3ExIRuvvlmff/731eKvzkAPjdX/Yhi0a6gXIRs27275IbpxnfzqeVuw26e9zDzrKCZZVk6cuSItmzZUvH3UNAMgFfCUrQrnXbqYFz+50gknOmGIISooOK8lwpsQbN8Pq9cLlfyAAAvhKF+RFjLsPsd570+vgoivb29isfj049kMmm6SwAiIuj1I6gIagbnvX6+CiJ79+5VNpudfoyOjpruEoCICHr9iDCM6AQR571+ri5WrVZTU5OamppMdwNABBXrR2Qy5T/dFteI+LV+RNBHdIKK814/X42IAIApQa8fEfQRnaDivNfP1SDyySef6MSJEzpx4oQkaXh4WCdOnNDIyIibbwsANQly/QgqgprBea+fq9t3BwcH1dnZOev49u3bdfDgwQW/n+27AEwIaj2I4u4NqXR6qXiR9HuYCirO+2zVXL89qyNSC4IIAFSnXD2LZNKZVoraxdBLnPdSBBEAiLCgjugEHef9kmqu377aNQMAqF9UyrD7Dee9NuyaAQAAxhBEAACAMQQRAABgDEEEAAAYQxABAADGEEQAAIAxBBEAAGAMQQQAABhDEAEAAMYQRAAAgDEEEQAAYAxBBAAAGEMQAQAAxhBEAACAMQQRAABgzBWmOwAAAMqzbWloSBobk1pbpfZ2KRYz3avGIogAAOBD6bS0a5d09uylY4mEtH+/lEqZ61ejMTUDAIDPpNNSd3dpCJGkTMY5nk6b6ZcbCCIAAPiIbTsjIYXC7NeKx3p6nHZhQBABAMBHhoZmj4RcrlCQRkeddmFAEAEAwEfGxhrbzu8IIgAA+Ehra2Pb+R1BBAAAH2lvd3bHWFb51y1LSiaddmFAEAEAwEdiMWeLrjQ7jBS/7usLTz0RgggAAD6TSkn9/dLq1aXHEwnneJjqiFDQDAAAH0qlpK4uKqsCABosCmW70RixmNTRYboX7iKIAICHolK2G6gUa0QAwCNRKtsNVIogAgAeiFrZbqBSBBEA8EDUynYDlSKIAIAHola2G6gUQQQAPBC1st1Apdg1AwDzadBe22LZ7kym/DoRy3JeD0vZ7mqxpTm6GBEBgLmk01Jbm9TZKW3b5jy3tdW0vSWoZbttWxoclA4fdp7dWEzbwNOMACKIAEA5Luy1DVrZbi8CAluaYRUK5QYJ/SGXyykejyubzaqlpcV0dwBEhW07V9y5trkU51GGh2sawgjCNEQxIMy8QhRHbxoRnFw+zTComus3IyIAMJPLe22LZbvvvdd59ttF1quaJ2xphkQQAYDZIr7X1quAEPHTjM+wawaoRRDG1lG7iO+19SogRPw04zOMiADVYol/+BX32s7c3lJkWVIyGdq9tl4FhIifZnyGIAJUgyX+0RDUvbYN4lVAiPhpxmcIIkCluGtZtARtr20DeRkQInya8Rm27wKVGhx0pmEWMjDgbIVAOER4PVA67WTvywcAk0knhDQ6IET4NIdSNddvFqsifNz6jcYS/2gq7rWNoFRK6uryJiBE+DRHHkEE4VLuI1wi4Ywz1/sRjiX+iCACAtzGGhGEh9sLSVniDwANRxBBOHixkJQl/gDQcJ4EkaefflptbW268sortWHDBr355ptevC2ixKtSkCzxB4CGcn2NyIsvvqjdu3fr2Wef1YYNG9TX16c777xTp06d0ooVK9x+e0SFlwtJvVzBFyLsigBQjusjIr/4xS+0Y8cOPfDAA7rxxhv17LPP6qqrrtJvf/tbt98aUeL1QlK/37XMZyhGC2AurgaRyclJvfXWW9q8efOlN1y0SJs3b9Ybb7wxq30+n1culyt5ABVhIalvUYwW8CfbdsojHT7sPJuqxehqEPnggw9k27ZWrlxZcnzlypU6f/78rPa9vb2Kx+PTj2Qy6Wb3ECYsJPUlitEC/uSnUUpf7ZrZu3evstns9GN0dNR0lxAkLCT1Ha/WEAOonN9GKV1drLps2TLFYjFduHCh5PiFCxd07bXXzmrf1NSkpqYmN7uEsGMhqa9QjBbwl4VGKS3LGaXs6vLu16arIyKLFy/W+vXrdfTo0eljU1NTOnr0qDZu3OjmWyPKWEjqG4EpRuuXyXLAZX4cpXR9++7u3bu1fft23XLLLbrtttvU19eniYkJPfDAA26/NQDDimuIM5nyn8Asy3nd6BpiN28LAPiMH0cpXQ8i3/3udzU+Pq6f/vSnOn/+vL72ta/pr3/966wFrADCp7iGuLvbCR2XhxFfrCEuTpbPTEnFyXLWFiFk/DhKaRUK5T6n+EM1txEG4F9e3k6+YrbtbBOYa5y6OFwzPMz0HkKj+M9+oVHKev/ZV3P99tWuGQDhlEpJZ85IAwPSoUPO8/Cw4cEGP06WAy7zY6UD16dmAEDy4e3k/ThZDnigWOmg3NIoE6OUBBEA0eTHyXLAI36qdEAQARBNgdjSA7jHL6OUrBEBEE1+nCwHIoggAiC6uC0AYBxTMwCizU+T5ZjFtvmrCTuCCAD4ZbIcJSh6Gw1MzQAAfMdvd4iFewgiAABfWegOsZJzh1juTRgOBBEAgK9Q9DZaCCIAAF+h6G20EEQAAL5C0dtoIYgAAHylWPR2Zp25Isty7t5M0dtwIIgAQCPYtjQ4KB0+7DyzkrJmFL2NFoIIgHDyMhik01Jbm9TZKW3b5jy3tbHHtA4UvY0Oq1Aot0HKH3K5nOLxuLLZrFpaWkx3B0BQeFkJq1jwYuav0uJHd66adaGyajBVc/0miABe4zeru7wMBrbtjHzMtde0eAff4WH+jhEp1Vy/mZoBvMQQvru8roRFwQugbgQRwCvUrHaf18GAghdA3QgigBeoWe0Nr4MBBS+AuhFEAC8whO8Nr4MBBS+AuhFEAC8whO8Nr4MBBS+AuhFEAC8whO8NE8GAghdAXdi+GxVsGTWruM0zkym/ToRtno1Vro5IMumEELeCAf+PAdOoI4JSXhZ3wtyKu2ak0jBC4St3EAwAYwgiuISqj/5i4pM6AHiMIAIHVR/9iU/qAEKumuv3FR71CSZUs2W0o8OzbkVeLMb5BoDPsGsmzNgyCgDwOUZEwowto8HHNA6AkGNEJMyo+hhs3CAPQAREM4jYtjQ4KB0+7DyH9f4eVH0MLm6QByAiohdEovYpk6qPweOnG+RFJbQDMCZa23ejXFODtQbBMTjoBOSFDAy4u/uGQngAasT23XIW+pRpWc6nzK6ucF6g2TIaHH7Y7TRXaC9ODYU5tAPwVHSmZrgNO4LC9G4nP00NAQi96AQRP3zKBCphercToR2Ah6ITREx/ygQqZXq3E6EdgIeiE0RMf8oEqmFytxOhHYCHorlrRuI27AgGE7udijdLzGTKrxPhZokAFlDN9Ts6IyISNTUQPMXdTvfe6zx7ceE3PTUEIFKiNSJSRE0NYGHl6ogkk04IIbQDmEc11+9oBhEAlSG0A6gBBc0ANAaF8AC4LFprRAAAgK8QRAAAgDEEEQAAYAxBBAAAGONaEHn88ce1adMmXXXVVbrmmmvcehsAABBgrgWRyclJbd26VY888ohbbwEAAALOte27P/vZzyRJBw8edOstALOosQEAdaOOCFCLclVHEwmnNDpVRwGgYr5arJrP55XL5UoegO8Ub554eQiRnJvEdXc7rwMAKlJVENmzZ48sy5r3cfLkyZo709vbq3g8Pv1IJpM1/yzAFbbtjISUuzNC8VhPj9MOALCgqu41Mz4+rg8//HDeNtdff70WL148/fXBgwfV09Ojjz/+eMGfn8/nlc/np7/O5XJKJpPcawb+MTgodXYu3G5gYP7S6KwvARBirt1rZvny5Vq+fHldnZtPU1OTmpqaXPv5QN3Gxupvx/oSAJjm2mLVkZERffTRRxoZGZFt2zpx4oQk6YYbbtDVV1/t1tsC7mptra9dcX3JzIHI4vqS/n7CCIBIqWpqphr333+/fve73806PjAwoI4K7+ZZzdAO4AnbltranOBQ7n8dy3JGN4aHZ0+1FL935iLXSr4XAAKkmuu3a7tmDh48qEKhMOtRaQgBfCkWc6ZQJCc4XK74dV9f+SAxNDR3CJGcYDM66rQDgIjw1fZdIBBSKWcKZfXq0uOJxPxTK41YXwIAIUNBM6AWqZTU1VXdzpd615cAQAgRRIBaxWLzb9Gdqb3dGTVZaH1Je3vDutgwbDcG4BKmZgCv1LO+xKR02llk29kpbdvmPLe1UUEWQEMQRAAv1bq+xBTK2QNwmWvbdxuB7bsIrSBMdbDdGECNXKusCqBBql1fYkI12439/mcB4FsEEWCmIIxWeIHtxgA8QBABLsd9YC5huzEAD7BYFShiYWap4nbjmTt8iixLSib9ud0YQGAQRADJmY7Ztat8fY/isZ4ep11UBHW7MYBAIYgAEveBmUvQthsDCBzWiAASCzPnU0s5ewCoEEEEkFiYuZAgbDcGEEhMzQASCzMBwBCCCCCxMBMADCGIAEUszAQAz7FGBLgcCzMBwFMEEWAmFmYCgGeYmgEAAMYQRAAAgDEEEQAAYAxBBAAAGEMQAQAAxhBEAACAMQQRAABgDEEEAAAYQxABAADGEEQAAIAxBBEAAGAMQQQAABjDTe8AVMe2uTsxgIYhiACoXDot7dolnT176VgiIe3fL6VS5voFILCYmgFQmXRa6u4uDSGSlMk4x9NpM/0CEGgEEQALs21nJKRQmP1a8VhPj9MOAKpAEAGwsKGh2SMhlysUpNFRpx0AVIEgAmBhY2ONbQcAnyGIAFhYa2tj2wHAZwgiABbW3u7sjrGs8q9blpRMOu0AoAoEEQALi8WcLbrS7DBS/Lqvj3oiAKpGEAFQmVRK6u+XVq8uPZ5IOMepIwKgBhQ0CxuqXsJNqZTU1cW/MQANQxAJE6pewguxmNTRYboXAEKCqZmwoOolACCACCJhQNVLAEBAEUTCgKqXAICAIoiEAVUvAQABRRAJA6peAgACiiASBlS9BAAEFEEkDKh6CQAIKNeCyJkzZ/Tggw9q7dq1WrJkidatW6d9+/ZpcnLSrbeMNqpeAgACyLWCZidPntTU1JSee+453XDDDXr77be1Y8cOTUxM6Mknn3TrbaONqpcAgICxCoVyxSfc8cQTT+iZZ57Re++9V1H7XC6neDyubDarlpYWl3sHAAAaoZrrt6cl3rPZrJYuXTrn6/l8Xvl8fvrrXC7nRbcAAIAhni1WPX36tJ566ik99NBDc7bp7e1VPB6ffiSTSa+6BwAADKg6iOzZs0eWZc37OHnyZMn3ZDIZ3XXXXdq6dat27Ngx58/eu3evstns9GN0dLT6PxEAAAiMqteIjI+P68MPP5y3zfXXX6/FixdLks6dO6eOjg59/etf18GDB7VoUeXZhzUiAAAEj6trRJYvX67ly5dX1DaTyaizs1Pr16/XgQMHqgohADxi2+y0AmCMa4tVM5mMOjo6dN111+nJJ5/U+Pj49GvXXnutW28LoBrptHPn5stvmphIOAXyqD0DwAOuBZHXXntNp0+f1unTp5VIJEpe83DHMIC5pNNSd7dzd+bLZTLOcQrhAfCAp3VEqsUaEcAlti21tZWOhFzOspyRkeFhpmkAVK2a6zeLNoAoGhqaO4RIzijJ6KjTDgBcRBABomhsrLHtAKBGBBEgiv7v/ypr19rqbj8ARB5BBIga25aef37hdomEs5UXAFxEEAGiZmjI2RmzkB07WKgKwHUEESBqKl338aUvudsPABBBBIieStd9sD4EgAcIIkDUtLc76z8sq/zrliUlk6wPAeAJgggQNbGYU8Jdmh1Gil/39bE+BIAnCCJAFKVSTgn31atLjycSlHYH4CnX7jUDwOdSKamrizvvAjCKIAJEWSwmdXSY7gWACGNqBgAAGEMQAQAAxhBEAACAMQQRAABgDEEEAAAYQxABAADGEEQAAIAxBBEAAGAMQQQAABhDEAEAAMYQRAAAgDEEEQAAYAw3vWsE2+YOpgAA1IAgUq90Wtq1Szp79tKxRELav9+5zToAAJgTUzP1SKel7u7SECJJmYxzPJ020y8AAAKCIFIr23ZGQgqF2a8Vj/X0OO0AAEBZBJFaDQ3NHgm5XKEgjY467QAAQFkEkVqNjTW2HQAAEUQQqVVra2PbAQAQQQSRWrW3O7tjLKv865YlJZNOOwAAUBZBpFaxmLNFV5odRopf9/VRTwQAgHkQROqRSkn9/dLq1aXHEwnnOHVEAACYFwXN6pVKSV1dVFYFAKAGBJFGiMWkjg7TvQAAIHCYmgEAAMYQRAAAgDEEEQAAYAxBBAAAGEMQAQAAxhBEAACAMQQRAABgDEEEAAAYQxABAADG+LqyaqFQkCTlcjnDPQEAAJUqXreL1/H5+DqIXLx4UZKUTCYN9wQAAFTr4sWLisfj87axCpXEFUOmpqZ07tw5NTc3y7Is194nl8spmUxqdHRULS0trr1P1HBe3cF5dQ/n1h2cV/f49dwWCgVdvHhRq1at0qJF868C8fWIyKJFi5RIJDx7v5aWFl/9RYYF59UdnFf3cG7dwXl1jx/P7UIjIUUsVgUAAMYQRAAAgDEEEUlNTU3at2+fmpqaTHclVDiv7uC8uodz6w7Oq3vCcG59vVgVAACEGyMiAADAGIIIAAAwhiACAACMIYgAAABjCCIz3HPPPVqzZo2uvPJKtba26nvf+57OnTtnuluBdubMGT344INau3atlixZonXr1mnfvn2anJw03bVQePzxx7Vp0yZdddVVuuaaa0x3J7CefvpptbW16corr9SGDRv05ptvmu5SKLz++uu6++67tWrVKlmWpZdeesl0lwKvt7dXt956q5qbm7VixQpt2bJFp06dMt2tmhFEZujs7NQf/vAHnTp1Sn/605/07rvvqru723S3Au3kyZOamprSc889p3feeUe//OUv9eyzz+rHP/6x6a6FwuTkpLZu3apHHnnEdFcC68UXX9Tu3bu1b98+/fvf/9bNN9+sO++8U++//77prgXexMSEbr75Zj399NOmuxIax44d086dO3X8+HG99tpr+vTTT3XHHXdoYmLCdNdqwvbdBfz5z3/Wli1blM/n9bnPfc50d0LjiSee0DPPPKP33nvPdFdC4+DBg+rp6dHHH39suiuBs2HDBt1666361a9+Jcm5z1UymdRjjz2mPXv2GO5deFiWpSNHjmjLli2muxIq4+PjWrFihY4dO6bbb7/ddHeqxojIPD766CP9/ve/16ZNmwghDZbNZrV06VLT3QA0OTmpt956S5s3b54+tmjRIm3evFlvvPGGwZ4Blclms5IU2N+pBJEyfvSjH+nzn/+8vvjFL2pkZEQvv/yy6S6FyunTp/XUU0/poYceMt0VQB988IFs29bKlStLjq9cuVLnz5831CugMlNTU+rp6dE3vvENffWrXzXdnZpEIojs2bNHlmXN+zh58uR0+x/+8If6z3/+o7/97W+KxWK67777xAzWbNWeV0nKZDK66667tHXrVu3YscNQz/2vlnMLIHp27typt99+Wy+88ILprtTsCtMd8MIPfvAD3X///fO2uf7666f/e9myZVq2bJm+/OUv6ytf+YqSyaSOHz+ujRs3utzTYKn2vJ47d06dnZ3atGmTnn/+eZd7F2zVnlvUbtmyZYrFYrpw4ULJ8QsXLujaa6811CtgYY8++qheeeUVvf7660okEqa7U7NIBJHly5dr+fLlNX3v1NSUJCmfzzeyS6FQzXnNZDLq7OzU+vXrdeDAAS1aFInBuJrV828W1Vm8eLHWr1+vo0ePTi+inJqa0tGjR/Xoo4+a7RxQRqFQ0GOPPaYjR45ocHBQa9euNd2lukQiiFTqH//4h/75z3/qm9/8pr7whS/o3Xff1U9+8hOtW7eO0ZA6ZDIZdXR06LrrrtOTTz6p8fHx6df4xFm/kZERffTRRxoZGZFt2zpx4oQk6YYbbtDVV19ttnMBsXv3bm3fvl233HKLbrvtNvX19WliYkIPPPCA6a4F3ieffKLTp09Pfz08PKwTJ05o6dKlWrNmjcGeBdfOnTt16NAhvfzyy2pubp5eyxSPx7VkyRLDvatBAdP++9//Fjo7OwtLly4tNDU1Fdra2goPP/xw4ezZs6a7FmgHDhwoSCr7QP22b99e9twODAyY7lqgPPXUU4U1a9YUFi9eXLjtttsKx48fN92lUBgYGCj773P79u2muxZYc/0+PXDggOmu1YQ6IgAAwBgm6gEAgDEEEQAAYAxBBAAAGEMQAQAAxhBEAACAMQQRAABgDEEEAAAYQxABAADGEEQAAIAxBBEAAGAMQQQAABhDEAEAAMb8P4DhGF4a4vyAAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i in range(40):\n",
    "    if real_class[i] == -1:\n",
    "        plt.scatter(X[i][0], X[i][1], c='b')\n",
    "    else:\n",
    "        plt.scatter(X[i][0], X[i][1], c='r')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最小平方误差判别有两种实现方式：伪逆法求解和梯度下降法求解。\n",
    "\n",
    "梯度下降法和上面的感知器实现的类似，只不过迭代公式有变化，就不实现了。梯度下降的实现也叫LMS算法。\n",
    "\n",
    "下面实现伪逆法：\n",
    "\n",
    "#### 1. 计算权向量$\\boldsymbol{\\alpha}$\n",
    "\n",
    "最小平方误差的判别函数和感知器类似，但是引入了新的常数向量$\\mathbf{b}$（这个向量是人为设置的）\n",
    "\n",
    "如果$\\mathbf{b}$把同一类样本的$b_i$设置为相同的值，就和Fisher线性判别等价；\n",
    "\n",
    "如果$\\mathbf{b}$是全为1的向量，最后MSE算法的解是贝叶斯判别函数的最小平方误差逼近，也就是平方误差的期望最小。\n",
    "\n",
    "判别函数为$g(\\mathbf{y}) = \\boldsymbol{\\alpha}^\\text{T}\\mathbf{y} - \\mathbf{b}$\n",
    "\n",
    "$\\boldsymbol{\\alpha}$的算法为：\n",
    "\n",
    "$$\\begin{align}\\boldsymbol{\\alpha} &= \\mathbf{Y}^+\\mathbf{b} \\\\ \\mathbf{Y}^+&=(\\mathbf{Y}^\\text{T}\\mathbf{Y})^{-1}\\mathbf{Y}^\\text{T}\\end{align}$$\n",
    "\n",
    "$\\mathbf{Y}^+$也叫$\\mathbf{Y}$的伪逆，可以用`np.linalg.pinv()`计算\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.02712775 -0.38992454 -0.33913745]\n"
     ]
    }
   ],
   "source": [
    "# 预处理数据和上面感知器是一样的\n",
    "Y = np.concatenate((np.ones((X.shape[0], 1)), X), 1)\n",
    "\n",
    "for i in range(20, 40):\n",
    "    Y[i] = -Y[i]\n",
    "\n",
    "# 这里我把b设为全1的向量\n",
    "b = np.ones(Y.shape[0])\n",
    "\n",
    "alpha = b @ np.linalg.pinv(Y).T\n",
    "print(alpha)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2. 分类\n",
    "\n",
    "分类和感知器一样，决策规则：\n",
    "\n",
    "$$若g(\\mathbf{y}) = \\boldsymbol{\\alpha}^\\text{T}\\mathbf{y} \\gtrless 0, 则\\mathbf{y}\\in\\begin{cases}\\omega_1\\\\ \\omega_2\\end{cases}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def classify_MSE(_Y: np.ndarray, _alpha: np.ndarray) -> np.ndarray:\n",
    "    \"\"\"\n",
    "    感知机方法分类\n",
    "    :param _Y: 增广后的样本矩阵（每一行为一个样本），*****注意这里只增广，不改变正负号*****\n",
    "    :param _alpha: 前面算出的alpha值\n",
    "    :return: 返回分类好的矩阵（一维向量）\n",
    "    \"\"\"\n",
    "    g = ((_alpha @ _Y.T) >= 0).astype(int)\n",
    "    for i, g_i in enumerate(g):\n",
    "        if g_i == 0:\n",
    "            g[i] = -1\n",
    "    return g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ans_MSE = \n",
      "[ 1  1  1  1  1  1  1  1  1  1 -1  1  1  1  1  1  1  1  1  1 -1  1 -1 -1\n",
      " -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1]\n",
      "正确率： 0.95\n"
     ]
    }
   ],
   "source": [
    "# 注意我们前面把Y预处理过，把原来的第二类样本y变成了-y，这里要变回来才能分类\n",
    "ans_MSE = classify_MSE(np.concatenate([Y[:20], -Y[20:]]), alpha)\n",
    "print(\"ans_MSE = \")\n",
    "print(ans_MSE)\n",
    "print(\"正确率：\", 1 - np.sum(abs(ans_MSE - real_class)) / 2 / X.shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "为了证明如果$\\mathbf{b}$把同一类样本的$b_i$设置为相同的值，就和Fisher线性判别等价，\n",
    "\n",
    "我又用Fisher线性分类算了一下结果。显然和MSE完全一致"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "m_i = \n",
      "[-1.15692593 -0.9509051 ]\n",
      "S_i = \n",
      "[[16.21931747  4.68690613]\n",
      " [ 4.68690613 12.48347856]]\n",
      "m_i = \n",
      "[0.94086771 1.039338  ]\n",
      "S_i = \n",
      "[[ 7.99324042 -1.16075032]\n",
      " [-1.16075032 13.21892985]]\n",
      "ans_Fisher = \n",
      "[ 1  1  1  1  1  1  1  1  1  1 -1  1  1  1  1  1  1  1  1  1 -1  1 -1 -1\n",
      " -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1]\n",
      "正确率： 0.95\n"
     ]
    }
   ],
   "source": [
    "w = get_w(X[:20], X[20:])\n",
    "ans_Fisher = classify_Fisher(X, w, get_mean(X[:20]), get_mean(X[20:]))\n",
    "print(\"ans_Fisher = \")\n",
    "print(ans_Fisher)\n",
    "print(\"正确率：\", 1 - np.sum(abs(ans_Fisher - real_class)) / 2 / X.shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 四、线性支持向量机\n",
    "\n",
    "线性支持向量机的思想是：\n",
    "\n",
    "如果样本是线性可分的，则找出两类样本的分界面，使边界样本到分界面的距离尽可能大。边界样本又叫支持向量。\n",
    "\n",
    "推广到线性不可分的样本，引入松驰变量$\\xi_i \\quad i = 1, 2, \\cdots, N$\n",
    "\n",
    "这一节的符号表示和前面不同，用别的教材比较常用的符号来表示。样本向量为$\\mathbf{x}_i$，类别为$y_i = 1$或$-1$\n",
    "\n",
    "线性支持向量机的算法如下\n",
    "\n",
    "#### 1. 求分类面参数\n",
    "\n",
    "这个算法的数学推导非常之复杂，书上说了跟没说一样，这篇文章讲的很详细，可以参考一下：[支持向量机的数学原理详细推导和一些优化技巧](https://zhuanlan.zhihu.com/p/413953469)\n",
    "\n",
    "具体算法如下：\n",
    "\n",
    "(1) 解决最优分类超平面的对偶问题（the dual problem）即通过下面的公式求出$\\boldsymbol{\\alpha}^*$（具体的求法书中没给出，我直接用sklearn调包实现了）\n",
    "\n",
    "$$\\begin{align}\n",
    "&\\max\\limits_{\\boldsymbol{\\alpha}}Q(\\boldsymbol{\\alpha}) = \\sum_{i=1}^{N}\\alpha_i - \\frac{1}{2}\\sum_{i,j=1}^{N}\\alpha_i\\alpha_jy_iy_j(\\mathbf{x}_i\\cdot\\mathbf{x}_j) \\\\\n",
    "&\\text{s.t.}\\sum_{i=1}^{N}y_i\\alpha_i = 0 \\\\\n",
    "&0 \\leqslant \\alpha_i \\leqslant C,\\quad i=1,\\cdots,N\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "其中$C$是指定常数，强调对分类错误的惩罚。\n",
    "\n",
    "(2) 通过下面的公式求出$\\mathbf{w}^*$\n",
    "\n",
    "$$\\mathbf{w}^* = \\sum_{i=1}^N\\alpha_i^*y_i\\mathbf{x}_i$$\n",
    "\n",
    "(3)通过下面的公式求出b\n",
    "\n",
    "$$y_i[(\\mathbf{w^*}\\cdot\\mathbf{x}_i) + b] - 1 + \\xi_i = 0$$\n",
    "\n",
    "需要注意这里的$\\mathbf{x}_i, y_i$是一个样本，求b的时候把所有的**边界支持向量**样本都用上，每个支持向量算一个b出来，然后算平均值，作为最终的b。\n",
    "\n",
    "支持向量是$\\alpha_i>0$的样本，分为两种：\n",
    "\n",
    "$0<\\alpha_i<C，\\xi_i=0$ 的叫边界支持向量（对应分类正确但处在分类边界面上的样本）\n",
    "\n",
    "$\\alpha_i=C， \\xi_i>0$ 的叫错分支持向量（对应分类错误的样本）。\n",
    "\n",
    "代码实现如下（我直接用sklearn了）："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 用SVC而不是LinearSVC的原因是LinearSVC没有提供查看支持向量的方法\n",
    "svc = SVC(kernel=\"linear\")\n",
    "_ = svc.fit(X, real_class)  # _没有意义，是为了防止输出"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "支持向量为：\n",
      "[[-0.46661866 -0.02766581]\n",
      " [ 1.17958567 -0.29744786]\n",
      " [ 0.30757167  0.34764882]\n",
      " [ 0.03075503  1.00573105]\n",
      " [ 0.80973014  0.20798233]\n",
      " [-0.01296092 -0.14118512]\n",
      " [ 0.18570744 -0.42184708]\n",
      " [ 0.07235573  0.22487309]\n",
      " [-0.14226139 -0.44542636]]\n",
      "对偶问题中的alpha系数中大于零的部分为：\n",
      "[[-1.         -1.         -1.         -0.40243771 -0.43370871  1.\n",
      "   1.          1.          0.83614642]]\n"
     ]
    }
   ],
   "source": [
    "alpha = svc.dual_coef_  # alpha是对偶问题的alpha系数中大于0的部分（即支持向量的alpha系数）\n",
    "support_vectors = svc.support_vectors_  # support_vectors是支持向量\n",
    "\n",
    "print(\"支持向量为：\")\n",
    "print(support_vectors)\n",
    "print(\"对偶问题中的alpha系数中大于零的部分为：\")\n",
    "print(alpha)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对偶问题的系数alpha大于零的部分分为两类：等于C（sklearn的默认C值是1）时是错分支持向量，在0和C之间是边界支持向量。\n",
    "\n",
    "注意sklearn实现时alpha带正负号，实际上是$y_i\\alpha_i$，和书上写的有点区别。\n",
    "\n",
    "下面可视化一下，更好的理解SVM。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGdCAYAAAAvwBgXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB760lEQVR4nO3dd3hUZdrH8e9Meg+hBpIQOoTeEkB67yCruzZERVQEBFkpwUJZJTQRKXZfcVlRUFroIBJAwBCqkITQIQQSQktvM3PeP44gKCXATM7k5P5cVy73TCYzv8yGzJ2n3I9BURQFIYQQQggNGLUOIIQQQoiSSwoRIYQQQmhGChEhhBBCaEYKESGEEEJoRgoRIYQQQmhGChEhhBBCaEYKESGEEEJoRgoRIYQQQmjGUesA92KxWLhw4QJeXl4YDAat4wghhBCiEBRFISMjg4oVK2I03nvMw64LkQsXLhAYGKh1DCGEEEI8hMTERAICAu55H7suRLy8vAD1G/H29tY4jRBCCCEKIz09ncDAwJvv4/di14XIjekYb29vKUSEEEKIYqYwyypksaoQQgghNCOFiBBCCCE0I4WIEEIIITQjhYgQQgghNCOFiBBCCCE0I4WIEEIIITQjhYgQQgghNCOFiBBCCCE0I4WIEEIIITQjhYgQQgghNCOFiBBCCCE0I4WIEEIIITRT4guRRXGL+HDvh2TmZ2odRQghhChxSnQhcj33OvMPzGdh7EL6rOzD6pOrURRF61hCCCFEiVFkhci0adMwGAyMGjWqqJ7yvnxdfZnZbiZBXkFczrnMhF8nMGjDII5ePap1NCGEEKJEKJJCJCYmhs8//5wGDRoUxdM9kLYBbVnRbwUjm4zEzdGNA5cO8K81/+L9394nLS9N63hCCCGErtm8EMnMzOTZZ5/lyy+/pFSpUrZ+uofi7ODMy/VfJrJ/JN2Du2NRLCxJWELvFb358diPmC1mrSMKIYQQumTzQmTYsGH06tWLzp073/e+eXl5pKen3/ZRlCp4VGBmu5l83fVrqvtW53redabsnsIz657h4KWDRZpFCCGEKAlsWoj88MMP7N+/n4iIiELdPyIiAh8fn5sfgYGBtox3V6H+ofzY50fGh47H08mTuCtxDFw/kHd+fYfLOZc1ySSEEELokc0KkcTEREaOHMl3332Hq6trob4mPDyctLS0mx+JiYm2indfjkZHnq3zLKsfX03/6v0BWHVyFX1W9OG/sf+lwFKgWTYhhBBCLwyKjfarrly5kscffxwHB4ebt5nNZgwGA0ajkby8vNs+dyfp6en4+PiQlpaGt7e3LWIW2qHUQ0RERxB7JRaAaj7VCA8LJ8w/TNNcQgghhL15kPdvmxUiGRkZnD179rbbXnzxRWrXrs24ceOoV6/efR/DngoRAItiYcXxFczZP4fredcB6Fq5K2Oaj6GCRwVtwwkhhBB2wi4KkTtp3749jRo1Ys6cOYW6v70VIjek5aUx/8B8lh5bikWx4OboxpD6QxhUdxDODs5axxNCCCE09SDv3yW6s+rD8nHx4e0Wb7O091KalGtCjimHuQfm0n9Vf7af3651PCGEEKLYKNIRkQdlryMit1IUhXWn1/Hh3g9JzUkF1CZp45qPI8g7SON0QgghRNGTEZEiZDAY6FW1F6sfX82LdV/E0eDI9vPb6b+qP3P3zyW7IFvriEIIIYTdkkLESjycPBjdbDTL+i2jpX9LCiwFfHn4S/qt6sfGMxvlMD0hhBDiDmRqxgYUReGXc78wI2YGF7IuABBWIYzwsHCq+VbTOJ0QQghhWzI1ozGDwUCnyp1Y2X8lQxsOxcXBhejkaJ6IfIIZMTPIyM/QOqIQQghhF6QQsSE3Rzdeb/Q6K/utpGNgR0yKiUVxi+izog+rTqzColi0jiiEEEJoSqZmitDOpJ1M2zONM+lnAGhYtiETwiYQUjpE22BCCCGEFdltQ7MHpbdCBKDAXMCi+EV8dugzckw5GDDwRM0neKPxG/i6+modTwghhHhkskbEjjk5OPFSvZdY3X81Pav0REHhx2M/0ntlb5YmLMVsMWsdUQghhCgyMiKisZjkGCL2RHD82nEA6vjVYULYBBqVa6RtMCGEEOIhydRMMWOymFiSsIQFBxfc3FHTt1pf3mz6JmXcymicTgghhHgwMjVTzDgaHXm2zrOs7r+aATUGYMBA5MlIeq/ozbex31JgKdA6ohBCCGETMiJihw6nHmZq9FSOXDkCQFWfqowPHU/Lii01TiaEEELcn0zN6IBFsbDyxErm7JvDtbxrAHSp3IUxzcbg7+mvcTohhBDi7mRqRgeMBiMDagxg9eOreab2MxgNRjaf3UzflX35/NDn5JnztI4ohBBCPDIpROycj4sP4WHhLO29lKblm5JrzmX+wfn0X9mfbYnbtI4nhBBCPBKZmilGFEVh/en1fLj3Qy7lXAKgTaU2jAsdR2XvyhqnE0IIIVQyNaNTBoOBnlV7svrx1bxU7yUcjY7sSNrB46se5+P9H5NdkK11RCGEEOKBSCFSDLk7ufNm0zdZ3nc5j1V8jAJLAV8d/oq+K/uy4fQG7HiQSwghhLiNTM0Uc4qisDVxKzNiZpCUmQRA8wrNCQ8Np0apGhqnE0IIURLJ1EwJYjAY6BjUkZX9VvJ6o9dxcXAhJjmGJ1c/yfQ900nPT9c6ohBCCHFXUojohKujK0MbDmVV/1V0DuqMWTHzv/j/0WdFH1YcX4FFsWgdUQghhPgbmZrRqV0XdjFtzzROp50GoEGZBkwIm0DdMnU1TiaEEELvpLPqA9gcl0J2vom+DStiMBhs8hxaKTAX8F38d3x66FOyTdkYMDCgxgDeaPIGfq5+WscTQgihU1KIFFJWnolOH24jOT2X0GA/JverSx1//Y28XMq+xEf7PmLNqTUAeDt7M7zxcJ6s+SSORkeN0wkhhNAbKUQKKc9k5qsdp5n3y3FyCywYDTCwRWVGd6mFj7uT1Z9Pa/tT9jM1eioJ1xIAqFWqFuFh4TQt31TjZEIIIfRECpEHlHQ9h6lr41l7+CIAfh7OjO1Wi382C8Ro1Nd0jdli5sdjPzLvwLybO2p6Ve3F6KajKedeTuN0Qggh9EAKkYe068RlJkbGcvxSJgANA3yY3K8ejQJ9bf7cRe1a7jXmHpjLsmPLUFBwd3TntYav8Vyd53By0N9okBBCiKIjhcgjKDBb+HbXGT7++TgZeSYA/tksgLHda1PG06VIMhSl2MuxTI2eyu+Xfwcg2DuY8NBwWlVqpXEyIYQQxZUUIlZwKSOX6esTWLb/PABero6M7lKTgS0q4+igr/YrFsVC5MlIPtr3EVdzrwLQKagTY5qPoZJnJY3TCSGEKG6kELGifWevMTHyCEeS1PUUtSt4MalvXVpULa1JHlvKyM/gk4Of8P3R7zErZlwcXBhcfzAv1n0RV0dXreMJIYQoJqQQsTKzRWFJTCIzNh7lenYBAH0aVmRCz9r4+7hplstWjl87zrQ909iTvAeASp6VGNt8LB0CO+iu14oQQgjrk0LERq5l5fPh5gQWR5/DooC7swPDO1ZncOsquDg6aB3PqhRFYePZjcyKmUVKdgoAj1V8jHGh46jiU0XjdEIIIeyZ3Rx69+mnn9KgQQO8vb3x9vamZcuWrF+/3pZPaVOlPJx5v399Ioe3plnlUmTnm5mxIYHuc3YQlXBJ63hWZTAY6B7cncj+kQypPwQnoxM7L+xkQOQAPtr3EdkF2VpHFEIIoQM2HRFZvXo1Dg4O1KhRA0VR+Pbbb5k5cyYHDhygbt37n3libyMit1IUhRUHkohYf5TUjDwAOtcpz3u9Qwgq7a5xOus7l36OaXumsSNpBwDl3Mrx72b/pkeVHjJdI4QQ4jZ2PTXj5+fHzJkzGTx48H3va8+FyA0ZuQXM3XKcb3aewWRRcHY08lrbqgxtXx03Z31N1wBsS9zG9JjpJGYkAtC0fFPCQ8Op5VdL42RCPDyzGXbsgIsXwd8f2rQBB/398xWiyNhlIWI2m/nxxx8ZNGgQBw4cICQk5G/3ycvLIy8v7+Z1eno6gYGBdl2I3HDiUgaTIuP49cRlACr5uvFu7zp0q1tBdyMGeeY8vo39li9//5Jccy5Gg5Gnaj3F641ex8fFR+t4QjyQ5cth5Eg4f/7P2wIC4OOPYcAA7XIJUZzZVSFy+PBhWrZsSW5uLp6enixevJiePXve8b6TJk1i8uTJf7u9OBQioE7XbDiSzPtr40m6ngNA6+plmNQ3hOrlvDROZ30XMy8ya+8sNp3dBEApl1KMajqK/tX7YzToq9eK0Kfly+GJJ+CvvwVv/O3w009SjAjxMOyqEMnPz+fcuXOkpaXx008/8dVXX7Ft2zbdjYjcKiffzKfbTvLZtpPkmyw4Gg28+Fgwb3SqgZer/tqn/3bxN6ZFT+Nk2kkA6pWux4SwCdQvW1/jZELcndkMwcG3j4TcymBQR0ZOn5ZpGiEelF0VIn/VuXNnqlWrxueff37f+xaHNSL3cu5KNlPWxPFzvLr9tayXCxN61qZ/o0q6m64psBTwffz3fHLoE7IKsgAYUGMAbzR+g9Ju+mv+Joq/qCjo0OH+99u6Fdq3t3UaIfTFbrbv3onFYrlt1EPPgkq789WgZnzzYnOqlPEgNSOPN5cc4snPdnMkKU3reFblZHTi+brPs+bxNfSt1heA5ceX02dFH76L/w6TxaRxQiFud/Gide8nhHg4Ni1EwsPD2b59O2fOnOHw4cOEh4cTFRXFs88+a8untTsdapVjw6g2jO1eCzcnB/aevUbf+b/yzsrDXM/O1zqeVZVxK8MHrT9gUY9F1PGrQ0ZBBtP2TOOfa/5JTHKM1vGEuMnf37r3E0I8HJtOzQwePJgtW7Zw8eJFfHx8aNCgAePGjaNLly6F+vriPjVzJxfTcpi67iirD10AoJS7E291q8VTzYNwMOprusZsMbPs+DLmHphLWp46AtQjuAejm42mgkcFjdOJku7GGpGkpL8vVgVZIyLEo7DrNSIPQo+FyA27T15hUmQsCSkZANSr5M3kvvVoWrmUxsms73rudeYdmMePx35EQcHN0Y1XG7zKwJCBODs4ax1PlGA3ds3A7cWI7JoR4tFIIVJMmMwWFv12ltmbj5GRq66h+EeTAMb1qEU5L/2ddht3JY6p0VM5lHoIgGDvYMaFjqN1pdYaJxMl2Z36iAQGwpw5UoQI8bCkEClmLmfmMWPDUZbuVX8Terk4MrJzDQa1CsbJQV/9OCyKhTWn1jB772yu5F4BoENgB8Y0H0OgV6DG6URJJZ1VhbAuKUSKqQPnrjExMpbfz6vrKWqU82Ry37q0ql5G42TWl5mfyaeHPmVx/GJMiglnozMv1X+Jl+q9hJujm9bxhBBCPAIpRIoxi0Vh6d5Epm84yrXsAgB61ffn7V51qOirvzfok9dPErEnguiL0QD4e/gztvlYOgV10l2vFSGEKCmkENGB69n5zN58jP/9dhaLAm5ODgzvWJ2X21TBxVFfY8aKorD57GZm7p1JclYyAC39WzI+bDxVfapqnE4IIcSDkkJER+IupDMx8ggxZ64BEFzanff6hNCxdnmNk1lfdkE2Xx/5mm+OfEOBpQBHgyPPhTzHqw1exdPZU+t4QgghCkkKEZ1RFIVVBy8wdV08lzLUrrQda5fjvd4hBJfx0Did9SWmJzIjZgZR56MAtUna6Kaj6V21t0zXCCFEMSCFiE5l5pmYt+U4X/96GpNFwdnByCttq/J6h2q4OztqHc/qtp/fzvQ90zmXcQ6AJuWaEB4WTm2/2honE0IIcS9SiOjciUuZTF4dy47jlwGo6OPK271C6Fm/gu5GDPLN+fw37r988fsX5JhyMBqMPFnzSUY0HoGPi4/W8YQQQtyBFCIlgKIobIxN4T9r4ki6ngNAq2qlmdS3LjXLe2mczvqSs5KZtXcWG89sBMDXxZc3mrzBgOoDcDDqa/GuEEIUd1KIlCA5+WY+23aSz7adJM9kwcFo4IVWwYzsXANvVyet41ndnot7iNgTwYnrJwAIKR3ChLAJNCzbUONkQgghbpBCpARKvJrNf9bEsSkuBYAyni6M71GbAY0rYdTZYXoFlgKWHF3CgoMLyCzIBKBftX6MajqKMm76a/4mhBDFjRQiJdi2Y6lMjozl1OUsAJoE+TKlXz3qVdLfeorLOZf5eP/HrDyxEgBPJ0+GNRrGv2r/Cyej/kaDhBCiuJBCpITLN1n4v52nmbvlONn5ZgwGeDo0iDFda1HKQ3+n3R5KPcTU6KnEXYkDoLpvdcJDwwn1D9U4mRBClExSiDyIi4fA6Ajl69rm8TWUnJZLxPp4Vh28AICPmxNvdavFM6FBOOhsusZsMbP8xHLm7p/L9bzrAHQL7sZbzd6igkcFbcMJIUQJI4VIYVnM8Hk7uBQHzV+GDhPAzdf6z6Ox6FNXmBgZy9HkDABC/L2Z0q8uzYL9NE5mfWl5acw/MJ+lx5ZiUSy4OboxpP4QBtUdhLOD/kaDhBDCHkkhUlg512H1GxC3Sr12LwOdJ0GjZ8FotP7zachktrB4zzlmbUwgPdcEwIDGlRjfozblvF01Tmd9R68eZWr0VA5cOgBAkFcQ40LH0TagrcbJhBBC/6QQeVAnt8L6sXD5mHpdqSn0nKn+V2euZOYxc2MCS/Ymoijg6eLIyE41eOGxYJwc9FV8KYrCmlNr+GjfR6TmpALQPqA9Y5uPJdA7UON0QgihX1KIPAxTPuz5HKKmQ34GYIAmA6HTRPDQ35bQQ4nXeS8ylkOJ1wGoXs6TSX3q0rqG/r7XrIIsPj/0OYviFmFSTDgbnXmh3gu8XP9l3BzdtI4nhBC6I4XIo8hIhp8nwaHv1WtXH+jwDjR7CRz0dZ6LxaLw475EZmxI4EpWPgA96lXg7V51CCjlrnE66zuVdopp0dPYfXE3ABU8KjCm2Ri6VO6iu9b4QgihJSlErOHcb7DuLUg+rF6Xrwc9ZkDwY0Wbowik5RTw0eZjLPrtLGaLgquTkWHtqzOkbVVcnfTVPl1RFH459wszYmZwIUvdTRTmH0Z4aDjVfKtpnE4IIfRBChFrsZhh30L45T+Qc029rf6T0GUKeFcs+jw2djQ5nYmrYok+fRWAID933u0dQuc65XQ3YpBjyuGbI9/w9eGvybfk42hw5Ok6TzO04VC8nPV3Vo8QQhQlKUSsLfuqWozs/QZQwMkD2o2FFq+Do762hCqKwurfL/LB2jhS0vMAaF+rLBP71KVKGQ+N01nf+YzzzIyZyS+JvwBQ2rU0o5uNpnfV3hgN+lq8K4QQRUUKEVu5cBDWjYHze9Tr0tWhx3So3lnTWLaQlWdi/tYTfLXjFAVmBWcHI4PbVGF4h+p4uOhrrQzAzqSdTNszjTPpZwBoWLYhE8ImEFI6RNtgQghRDEkhYksWC/y+BDa/B1mX1Ntq94ZuH0CpYE2j2cKp1Ewmr45j2zF1+2sFb1fe7lWH3g38dTddU2AuYFH8Ij479Bk5phwMGHiy5pOMaDwCX1dfreMJIUSxIYVIUchNU7f6Rn8GihkcXeGxUdB6FDjpa0uooij8HH+JKWtiSbyaA0CLqn5M6luX2hXs7P8XK0jJSuHDfR+y/vR6AHxcfHij8Rv8o8Y/cDDqa/GuEELYghQiRelSvNoM7fR29do3CLpNVUdJdDZikFtg5ovtp1iw9QR5JgsORgMDW1TmzS418XHT32m3MckxROyJ4Pi14wDU8avDhLAJNCrXSNtgQghh56QQKWqKoraJ3/g2pJ9Xb6vWEbpPh7I1tc1mA+evZfPB2njWH0kGoLSHM+N61OaJJgEYdXaYnsliYknCEhYcWEBGgXpWT99qfXmz6ZuUcdNf8zchhLAGKUS0kp8Fv34EOz8Gc756qm+L19UdNi762xL66/HLTIw8wsnULAAaBfoyuW9dGgb6ahvMBq7kXGHugbksP74cAA8nD4Y2HMozdZ7Byai/0SAhhHgUUoho7cpJ2DgBjm1Qrz0rQNf/qD1IdDZdk2+y8O2uM3y85TiZeSYMBvhXs0DGdKtFaU8XreNZ3eHUw0yNnsqRK0cAqOpTlfCwcFr4t9A4mRBC2A8pROzFsY2wfhxcO61eB7WCnjOgQn1tc9nApfRcpq0/yvIDSQB4uzryVrdaPBMahKPODtOzKBZWnljJnH1zuJanNrrrUrkLY5qNwd/TX+N0QgihPSlE7ElBLuyeDzs+hIJsMBih2WDoMAHc/bROZ3V7z1zlvVWxxF1MB6B2BS+m9KtHaBX9fa9peWl8cvATfkj4AYtiwdXBlZfrv8wL9V7AxUF/o0FCCFFYUojYo7TzsOkdiF2hXrv5QeeJ0Hgg6GxLqNmisHjPOWZtTCAtpwCAfo0qMqFnHcp7u2qczvoSriYQsSeCfSn7AAjwDGBc6DjaBbTTXa8VUbKZzbBjB1y8CP7+0KYNOOjr15ewErspRCIiIli+fDlHjx7Fzc2NVq1aMX36dGrVqlWor9dVIXLDqW3qdE1qvHpdsTH0nAUBzbTNZQNXs/KZtSmB7/ecQ1HAw9mBNzrV4MXHquDsqK/pGkVRWH96PR/u/ZBLOWqjuzaV2jAudByVvStrnE7cT3F4g9U64/LlMHIknD//520BAfDxxzBgQNHlEMWD3RQi3bt356mnnqJ58+aYTCYmTJjAkSNHiIuLw8Pj/ueW6LIQATAXwJ4vISoC8tQpDBo9p46QeJbTNpsNHD6fxnuRRzhw7joAVct6MKlPXdrWLKttMBvILsjm898/579x/8VkMeFkdGJQ3UEMqT8Edyd3reOJOygOb7BaZ1y+HJ54Qu1UcKsbA34//WQ/r5WwD3ZTiPxVamoq5cqVY9u2bbRt2/a+99dtIXJD5iX4eRIc/E69dvFW1440HwIO+jrPxWJRWLb/PNM3HOVyZj4AXUPK827vEAL99PcGfTrtNNNjprMzaScA5dzLMabZGLoFd5PpGjtSHN5gtc5oNkNw8O1F0F9zBATA6dP2N4oktGO3hciJEyeoUaMGhw8fpl69en/7fF5eHnl5eTev09PTCQwM1G8hckPiHlj3Flw8pF6XC4EeM6BKG21z2UB6bgEf/3ychbvOYLYouDgaGdq+Gq+1q4ark75+iymKQlRiFNNjppOUqe4mal6hOeGh4dQoVUPbcKJYvMHaQ8aoKOjQ4f7327oV2re3TQZR/DxIIVJkE/UWi4VRo0bx2GOP3bEIAXVNiY+Pz82PwMDAooqnrcBQGLIVes8Bt1JwKQ6+7Q0/vghpSVqnsypvVyfe7R3C+pFtaFm1NHkmC3N+Pk7n2dvYGJuMHa+dfmAGg4EOQR1Y2W8lrzd6HRcHF2KSY3hy9ZNM3zOd9Px0rSOWaDt23P0NHtQRiMRE9X5asYeMFy9a935C/FWRFSLDhg3jyJEj/PDDD3e9T3h4OGlpaTc/EhMTiyqe9owO0OxFGLFf3d5rMELscpjfDHbMBlPe/R+jGKlZ3ovFQ8KY/0xj/H1cOX8th1cX7eP5/9vDydRMreNZlaujK0MbDmVV/1V0DuqMWTHzv/j/0WdFH1YcX4FFsWgdsUQqDm+w9pDRv5CtcQp7PyH+qkimZoYPH86qVavYvn07VapUKfTX6X6NyL1cPATrxkLib+q1XzXoMR1qdNE2lw1k55tYsPUEX24/Tb7ZgpODgZceq8KITjXwdNHXWhmAXRd2MW3PNE6nqY3uGpRpwISwCdQtU1fjZCVLcZhysIeMN6aHkpL+vk4F7GMKS9gfu1kjoigKI0aMYMWKFURFRVGjxoPNi5foQgTUf/W/L4XN70JminpbzR7QfSr4VdU2mw2cuZzFlDVx/HJU3f5azsuFt3vVoW/Dirpb4FlgLuC7+O/49NCnZJuyMWBgQI0BjGwyklKupbSOVyIUhzdYe8l4Y8Es3J7Dnhb1CvtiN4XI66+/zuLFi1m1atVtvUN8fHxwc3O779eX+ELkhtx02DYdoj8DiwkcXOCxkdD6TXDW346TLfEpTFkTx9kr2QCEBvsxqW9dQirq72cgNTuV2ftms+bUGgC8nL0Y0XgET9Z8Ekej/kaD7E1xeIO1l4x32kIcGAhz5mj/Ggn7YzeFyN3+iv3mm2944YUX7vv1Uoj8RWoCrB8Lp6LUa59A6PYB1Omru8P0cgvMfLXjFPO3niC3wILRAM+1qMy/u9TCx11/p93uT9nP1OipJFxLAKBWqVqEh4XTtHxTjZPpX3F4g7WXjFo3VRPFh90UIo9KCpE7UBSIj4SNb0PaH4t5q7RTt/uWq61tNhtIup7D1LXxrD2srsbz83BmbLda/LNZIEajvoovs8XMj8d+ZN6BeTd31PSq2ovRTUdTzl1/je7sSXF4gy0OGYW4QQqRkiA/G379CHZ+DOY8MDpC2GvQbhy46u+12nXiMhMjYzl+Sd1R0yDAh8l969I4SH/rKa7lXmPugbksO7YMBQV3R3dea/gaz9V5DicH/Y0GCSH0RwqRkuTqadg4ARLWqdee5aHLFGjwL91N1xSYLXy76wxzfj5OZp4JgCebBjCuR23KeOrvtNvYK7FMjZ7K76m/AxDsHUx4aDitKrXSOJl4GDKiIUoSKURKouOb1cP0rp5UrwPDoOdM8G+obS4buJSRy/T1CSzbr06Ye7k6MrpLTQa2qIyjg74O07MoFiJPRvLRvo+4mnsVgE5BnRjTfAyVPCtpnE4UltZnxQhR1KQQKalMebB7AWyfBQVZgEFtktbxXXD30zqd1e07e42JkUc4kqSup6hV3otJfevSslppjZNZX0Z+Bp8c/ITvj36PWTHj4uDC4HqDebHei7g6umodT9yD1mfFCKEFKURKurQktffIkWXqtVsp6PQeNBmkdnDVEbNF4YeYc8zcmMD17AIAejfw5+1edfD3uf8W8eLm+LXjTNszjT3JewCo5FmJsc3H0iGwg+56reiBPZwVI4QWpBARqtM71O2+l+LUa/+G0GMmBIVpm8sGrmXl8+HmBL6LPoeigJuTAyM6VWdw6yq4OOrrN7yiKGw8u5FZMbNIyVYb3T1W8THGhY6jik/hOxcL27OHzqhCaMEuD70TGqjSBl7dAd2ng4uP2jb+/7rCitcgI0XrdFZVysOZ9/vXZ/Xw1jStXIqcAjMzNiTQfc4OtiZc0jqeVRkMBroHdyeyfyRD6g/ByejEzgs7GRA5gNn7ZpNVkKV1RPEHezgrRgh7JyMiJUVmKmyZBAf+p147e0GHcAh9BXS2JVRRFFYcSCJi/VFSM9TDAjvXKce7vUOoXNpD43TWdy79HNP2TGNHknoEazm3coxuNpqeVXrKdI3GZERElFQyNSPu7vxeWPcWXDigXpetrR6mV7W9prFsISO3gLlbjvPNzjOYLArOjkZea1uVoe2r4+asr+kagG2J25i2ZxrnM9UFCU3LNyU8NJxafrXu85XCVuzlrBghipoUIuLeLBY4sAi2TIbsK+ptIf2h6/vgG6hpNFs4npLBpNWx7Dyhfq+VfN14p1cduteroLsRgzxzHt/GfsuXv39JrjkXo8HIU7We4vVGr+Pj4qN1vBLJXs6KEaIoSSEiCifnGmydCjFfgWIBRzdo+29oOQKc9LUlVFEUNhxJ5v218SRdzwGgdfUyTOobQvVyXhqns76LmReZtXcWm85uAqCUSylGNR1F/+r9MRpkaVhRs5ezYoQoKlKIiAeTfBjWjYVzu9TrUsHqAtda3TWNZQs5+WY+jTrBZ9tPkW+y4Gg08OJjwbzRqQZervpaKwPw28XfiIiO4FTaKQDqla7HhLAJ1C9bX+NkJY90VhUliRQi4sEpChz+Se0/kvHHEv4a3aB7BJSupm02Gzh3JZspa+L4OV7dPVTWy4XwHrV5vHEl3U3XFFgK+D7+ez459MnNHTUDagzgjcZvUNpNf83fhBDak0JEPLy8DNg+E3Z/ApYCcHCGViOgzb/BWX87TrYevcTk1bGcuZINQLPKpZjUty71KulvPcXlnMt8tO8jIk9GAuDl5MWwxsP4V61/4Wh01DidfZPRDCEejBQi4tFdPq42Qzv5i3rtXQm6faAuatXZiEGeyczXv55m3pYT5BSYMRrgmbAg3upaC193Z63jWd3BSweZGj2V+KvxANQoVYMJoRNoVqGZxsnsk5wTI8SDk0JEWIeiwNG1sDEcrp9Tbwtuox6mV66Ottls4GJaDlPXHWX1oQsAlHJ34q1utXiqeRAORn0VX2aLmWXHlzH3wFzS8tIA6FGlB/9u+m/Ke5TXOJ39sOU5MTLKIvRMChFhXQU5sPNj+PUjMOWCwQHCXoX248FVf1MYu09eYVJkLAkpGQDUq+TN5L71aFq5lMbJrO967nXmH5zP0oSlKCi4ObrxaoNXGRgyEGcH/Y0GPQhbnhMjoyxC76QQEbZx7SxsnABH16jXHmWhyxRo8BQY9bUl1GS2sOi3s8zefIyMXBMA/2gSwLgetSjnpa+tzQDxV+KZGj2Vg6kHAQj2DmZc6DhaV2qtbTAN2aorqpzGK0oCKUSEbZ34GdaPgysn1OuA5up0TcXG2uaygcuZeczYcJSle9U/Xb1cHBnZuQaDWgXj5KCv4ktRFNacWsOHez/kSq7a/K1DYAfGNh9LgFeAxumK3vffwzPP3P9+ixfD008X7jH1chqvTCuJ+5FD74RtVe8MQ3dD58ng5AHnY+CLDrB6FGRf1TqdVZXxdGHGEw1Z8XorGgT4kJFn4v218fT8eAe7TlzWOp5VGQwG+lTrw5rH1zAoZBCOBke2Jm6l38p+LDi4gBxTjtYRi5S/v3XvB+qb992KEFBHSRIT1fvZq+XL1WKqQwe1UOvQQb1evlzrZKK4khER8WjSL8Dm9+Dwj+q1qy90eheavghGff2JZLEoLN2byIyNCVzNygegV31/3u5Vh4q+bhqns76T108SsSeC6IvRAPh7+DO2+Vg6BXXSXa+VO7HFOTG2GGUpSjKtJApLpmZE0TuzU93um3JEva5QH3rOgqAW2uaygbTsAmZvTmDRb2exKODm5MDwjtV5uU0VXBz1VXwpisLP535mRswMkrOSAWjp35LxYeOp6lNV43S2Z+1zYorzabx6mVYSRUMKEaENswn2/h9sfR9y1S2hNHgKukwGrwraZrOBuAvpTIw8QsyZawAEl3bnvT4hdKytv+2vOaYcvj78Nd8c+YZ8Sz6OBkeeC3mOVxu8iqezp9bxbMqa58QU59N4i3MRJYqeFCJCW1mX1ZN99y8CFHD2gvbjIOw1cNDXeS6KohB56AIfrI3nUkYeAJ1ql+O9PiFULq2/TrSJGYnMiJlBVGIUAGXcyjC66Wh6V+2t6+kaay7OLK6n8Rb3aSVRtKQQEfYhaR+sG6P+F6BMTegxHap11DaXDWTmmZj3y3H+79fTFJgVnB2MvNK2Kq93qIa7s/7ap28/v53pe6ZzLkNtdNe4XGMmhE2gtl9tjZMVD8Wxj4iMiIgHIbtmhH2o1BQG/wz9FoB7Gbh8DBY9DksG/tmpVSc8XRwJ71GH9SPb0qZGGfLNFuZvPUHnD7ex9veL2HG9/1DaBrRlRb8VjGwyEjdHNw5cOsC/1vyL9397/2anVnFvf/2RsPcfkTZt1GLpXgNffn7q6JHZXHS5RPEnIyKiaORch6gI2PMlKGZwdIM2o6HVG+CkrwZhiqKwKS6F/6yJ4/w1dctrq2qlmdy3LjXKe2mczvqSs5KZvXc268+sB8DXxZc3mrzBgOoDcNDZzilrKM47T+42rfRX9j66I2xPpmaE/UqJhXVj4eyv6nWpYOgWAbV66O4wvdwCM59tO8mnUSfJM1lwNBoY1CqYkZ1r4O2qr7UyADHJMUyNnsqJ62qju5DSIUwIm0DDsg01TmY/9LDz5E7TSn9VHIoqYVtSiAj7pigQuxw2vgMZ6gFzVO8C3adBmeraZrOBxKvZvL82jo2xKYDaJG18j9oMaFwJo84O0zNZTCxJWML8A/PJLMgEoF+1foxqOooybmU0Tqc9vayzMJvV7+Wf/4Srd+lhWByKKmE7skZE2DeDAer9A4bHQOs3wegEJzbDJy3g50mQl6l1QqsK9HPn84HN+PalUKqW8eByZh5v/XiIJz7bxZEkfa2ncDQ68mydZ1n9+Gr6V+8PwKqTq+izog//i/sfBZYCbQNq7OJF695PKw4O6sfdihAoHl1ihX2QQkRox8UTOk+C139TR0QsBeoJv/Obw+Gf7H/13gNqV7MsG0a1JbxHbTycHdh/7jp95v/KhBWHufZHp1a9KONWhv889h++6/kddUvXJbMgk+kx0/nn6n+y5+Kev90/KysLg8GAwWAgKytLg8RFwxZt47Wil6JKaE+mZoR9UBQ4tkE9TO/6WfW2yq2h5wwoX1fbbDaQkp5LxLp4Vh5Up6Z83Jx4q1stngkNwkFn0zUWxcLy48v5eP/HXM+7DkC34G681ewtKnioje7S07Pw8VEbo61bl0nXrh66HM4vzg3N/kov00zCNuxmjcj27duZOXMm+/bt4+LFi6xYsYL+/fsX+uulECmBCnJh11zYMRtMOWBwgOYvQ4dwcCuldTqriz51hYmRsRxNzgAgxN+bKf3q0izYT+Nk1peWl8b8A/NZemwpFsWCm6MbQ+oPwef4IN58o4ALF250aM0kIMDDLnddWKOxWXFtaPZXD1tUycm9JYPdrBHJysqiYcOGLFiwwJZPI/TEyRXajYXhe6BOX3Wr757PYV4ztVOrxaJ1QqsKq1qaNSNaM7lvXbxdHYm7mM4Tn+1m9JKDXErP1TqeVfm4+PB2i7dZ0nsJjcs1JseUw9wDc3nn5OOkl9l1232TktQ3a3s60dVap84OGKAWG5Uq3X57QEDxKUJALR4+/lj933/d8Hbjes6c24sMOblX3EmRTc0YDAYZEREP7uRWdbrmcoJ6Xakp9Jyp/ldnrmTmMXNjAkv2JqIoapO0kZ1q8MJjwTg56Gs5V2ZmJuvPbGDS1vkYva8AkH4gneTFyeSnJgJqe/xKlSAuDry9tW2Xb4veH3oZGSjsWTzFuX+KeHB2MzVz2xMVohDJy8sjLy/v5nV6ejqBgYFSiJR0pnx1VCRqOuRnAAZoMhA6TQQP/W0JPZR4nfciYzmUeB2A6uU8mdSnLq1r6Od7vXEujdHVSNm+ZSnTtQwGRwOWAguX118mdU0qSv6fv5q0XMqmh94ftna/okpew5Kn2BYikyZNYvLkyX+7XQoRAUBGMmyeCL//oF67+kCHd6DZS+Cgr/NcLBaFn/afZ/r6o1z5Y0dNj3oVeLtXHQJKuWuc7tH99YA8Z39nKj5bEc966jqR/Cv5JH+fTPredEDbQkQWZT46eQ1LHrtZI/KgwsPDSUtLu/mRmJiodSRhT7wqwIDP4aWNUKE+5KbB+jHweVs4s1PrdFZlNBr4Z7NAfnmrPS+0CsZogPVHkuk8extztxwnt6B4H+aRmZnJunWZgPqRfzGRM7POcHbuWfIvl8e5tDNBw4MIHvMC30Qe1jSrbFN9dPIainuxq0LExcUFb2/v2z6E+JugFvDKNuj1Ibj6wqVYWNgTlr0M6Re0TmdVPm5OTOpbl3Uj2xBWxY/cAguzNx+j60fb2RyXUmwP0/Pw8KBrVw8CAjwwGDy4sSYkY38Gxycs4dLKoVgKnPGsu5ePrw1kRswMMvIzNMmqp94fWpHXUNyLXU3N/JUsVhX3lXUFfvkP7FsIKODkoe66afE6ODprnc6qFEVhze8X+WBtPMl/7KhpX6ssE/vUpUoZbRdzPqw/FzBmAX9u3zUYPHAqm0jf2TM5atoKQGnX0oxuNpreVXtjNBTd31B66v2hFXkNSx67mZrJzMzk4MGDHDx4EIDTp09z8OBBzp3T1xHwQkMepaHPHHhlKwQ0h4Is+HkifNoSTvysdTqrMhgM9GlYkS3/bsfQ9tVwcjAQlZBKt4+2M33DUbLyTFpHfGA3trJWrHj77QEB8P2ngfz47Fw+7fwpwd7BXMm9wtu/vs3z658n7kpckWV8mG2q4nbyGop7semISFRUFB3usEJp0KBBLFy48L5fLyMi4oFYLOpC1s0TIeuSelvt3tDtA/WUX505lZrJ5NVxbDuWCkAFb1cm9KpDnwb+f1sMau/u11m1wFzAovhFfHboM3JMORgw8ETNJ3ij8Rv4uvoWScbCblMVdyevYclhl7tmHoYUIuKh5KapW32jP1Mbojm6wmOjoPUocHLTOp1VKYrCz/GXmLImlsSrOQCEVfFjcr+61K6gv38zKVkpzN43m3Wn1wFqk7Q3Gr/BP2r8Awej7f+c1kvvDy3Ja1gySCEiBMCleFg/Fk5vV699g6DbVHWUpJiNGNxPboGZL7afYsHWE+SZLDgYDQxsUZk3u9TEx81J63hWtzd5LxF7Ijh27RgAdfzqMCFsAo3KNdI2mBACkEJEiD8pCsSthI3vQPof48HVOkKPGVCmhqbRbOH8tWw+WBvP+iPJAJT2cGZc99o80TQAo84O0zNZTCxNWMr8g/Nv7qjpU7UPbzZ9k7LuZTVOJ0TJJoWIEH+Vn6UepLdrLpjzwegELYaqO2xcvLROZ3U7jqcyKTKWk6lZADQM9GVK37o0DPTVNpgNXM29ytz9c1l+fDkKCh5OHgxtOJRn6jyDk1F/o0FCFAdSiAhxN1dOwsYJcGyDeu1ZAbq+D/Wf0N10Tb7Jwre7zjDn52Nk5ZsxGOBfzQIZ060WpT1dtI5ndUcuH2Fq9FQOX1YboFX1qcr40PG0rNhS42RClDxSiAhxP8c2qofpXTutXge1gp4z1I6tOnMpPZeI9UdZcSAJAG9XR97qVotnQoNw1NlhehbFwqoTq5izfw5Xc68C0KVyF8Y0G4O/p3TLEqKoSCEiRGEU5MLu+bDjQyjIBoMRmg2Gjm+DWymt01ldzJmrvLcqlviL6vkttSt4MaVfPUKr+GmczPrS89P55OAn/HD0B8yKGVcHV16u/zIv1HsBFwf9jQYJYW+kEBHiQaSdh03vQOwK9dq9NHR6DxoPhCLYElqUzBaFxdFnmbXpGGk5BQD0b1SR8J51KO/tqnE66zt27RgR0RHsTdkLQIBnAONCx9EuoF2x67UiRHEihYgQD+PUNnW7b+pR9bpiY+g5CwKaaZvLBq5m5TNzYwI/xJxDUcDD2YE3OtXgxceq4Oyor+kaRVHYcGYDs2JmcSlHbXTXplIbxoWOo7J3ZY3TCaFPUogI8bDMBbDnS4iKgDx1CoNGz0HnieBZTttsNnD4fBrvRR7hwLnrAFQt68GkPnVpW1N/21+zC7L54vcv+DbuW0wWE05GJ54PeZ5XGryCu5O71vHEA5LGaPZNChEhHlXmJbVV/KHF6rWLN3SYAM2HgIOjttmszGJRWLb/PNM3HOVyZj4A3eqW551eIQT66e8N+kzaGabFTGNn0k4AyrmXY0yzMXQL7ibTNcXEnVrFBwSo59lIq3j7IIWIENaSuAfWvQUXD6nX5ULUZmhV2mibywbScwuYs/k43+4+g9mi4OJoZGj7arzWrhquTvr6U1NRFKISo5geM52kTHU3UfMKzQkPDadGKf01utOTP09svv32GzXkTz9JMWIPpBARwposZtj/X9gyGXKuqbfVHaD2H/GppG02G0hIzmBSZCy7T10BIKCUG+/2DqFrSHndjRjkmnJZGLuQrw5/RZ45DweDA0/XfpqhjYbi7Sy/c+yN2QzBwbePhNzKYFBHRk6flmkarUkhIoQtZF+FX96Hfd+AYgEnd2g7BloOA0d9bQlVFIW1hy/ywdp4LqblAtC2Zlkm9gmhWllPjdNZX1JmErNiZvHzuZ8B8HP1Y1STUfSr3g+jQV+Ld4uzqCi4w4Huf7N1K7Rvb+s04l6kEBHCli4egnVjIDFavfarBt2nQc2u2uaygex8Ewu2nuDL7afJN1twcjDwUusqjOhYA08Xfa2VAdh1YRfT9kzjdJra6K5BmQZMCJtA3TJ1NU4mAL7/Hp555v73W7wYnn7a9nnE3UkhIoStKQr8vgQ2vweZKeptNXtA96ngV1XbbDZw5nIWU9bE8ctRdftrOS8X3u5Vh74NK+puuqbAXMDio4v55OAnZJuyMWBgQI0BjGwyklKu+mt0V5zIiEjxIYWIEEUlNx22TYfoz8BiAgcXeGwktH4TnPW342RLfApT1sRx9ko2AKHBfkzuV5c6/vbx7zM/Px9nZ+eH/vytUrNTmb1vNmtOrQHA29mb4Y2H82TNJ3E06m80qDi4sUYkKenvi1VB1ojYkwd5/5bJTyEehas3dPsAhu6Cqu3BnAfbZ8CCUIhbdefflsVYpzrl2TiqLW91rYmrk5E9Z67Sa+4O3lt1hLTsAk2zLVmyhPr165OYmIjZrP71/P336n/NZkhMTKR+/fosWbKkUI9X1r0sEW0i+Lb7t9QqVYv0/HSmRk/lqTVPsS9ln02/F3FnDg7qFl34+xmVN67nzJEipLiREREhrEVRID4SNr4NaYnqbVXaqdt9y9XWNpsNJF3PYeraeNYevgiAn4czY7vV4p/NAjEai3a6Jj8/n/r163Ps2DHKl6+KwRBFcnLgzc9XqJCIorQnJeUUNWvW5PDhw4UeGQEwW8z8eOxH5h2YR3q+2uiuV9VejG46mnLu+mt0Z+/u1EckMFAtQmTrrn2QqRkhtJSfDb9+BDs/VkdIjI4Q9hq0G6eOoOjMrhOXmRgZy/FLmQA0DPBhcr96NAr0LdIciYmJNG+uFhtQFYgCAoFEoD1wivLlqxITE0VgYODdH+geruVeY+6BuSw7tgwFBXdHd15r+BrP1XkOJwcnK30nojCks6p9k0JECHtw9TRsnAAJ69Rrz/LQZQo0+Nffx5WLuQKzhW93nWHOz8fJzDMB8M9mAYztXpsynkWztdlshoCARJKT2wM3ipFFwMCb1/7+USQmBj7yG1bslVimRk/l99TfAQj2DiY8NJxWlVo92gMLoRNSiAhhT45tgg3j4epJ9TqwBfScAf4Ntc1lA5cycpm+PoFl+9Uxcy9XR0Z3qcnAFpVxdLj3krRHXWj6546KP0dA/vTnCIm1dlRYFAuRJyP5aN9HXM29CkCnoE6MaT6GSp76a3QnxIOQxapC2JOaXeH13dBpotoELfE3+KI9rBmtNknTkXJernz4z4YsG9qSepW8ycg1MXl1HL3n/cpvf3RqvZNbF5reaaVpYRaaXrx4438Foo6E3GrRH7ffer9HYzQY6V+9P2seX8NzdZ7DweDAlnNb6LeyH58e/JRcU+7fviYrKwuDwYDBYCArK8s6Qf5wpwW6QhQHMiIiRFFKS4LN78KRZeq1Wyno9B40GQRGfU1wmy0KP8ScY+bGBK7/saOmT8OKTOhZG38ft5v3u3WhadXy5YkyGAhMTr75+cQKFWivKJxKSbnnQtOiHhH5q+PXjjNtzzT2JO8BoJJnJcY2H0uHwA43e62kp2fh46N2pl23LpOuXT2ssq5BDoET9kamZoSwd6d3wPqxcClOvfZvCD1nQWCotrls4FpWPh9uTuC76HMoCrg7OzC8Y3UGt66Ci6P6LpyYmEj75s05lZJyl2WmqEVKTMxdF5oW5RqRu1EUhY1nNzIrZhYp2Wqju8cqPsa40HEc2FKFESOyuHDhRov8TAICPB65WJBD4IQ9kkJEiOLAbIKYr2DrVMhLU29r+Ax0ngRe5TWNZgtHktKYGBnLvrPqwYFVyngwsU8I7WuVU6dfAgJon5x8lxICovz9CUxMvOvWiKLYNVNY2QXZfHX4KxbGLqTAUoADjqSsG8ilyOew5N74/zYTg8EDePhiQQ6BE/ZKChEhipPMVNgyCQ78T7129oIO4RD6CuhsS6iiKKw4kMTUdUe5nJkHQOc65XmvbDpBvTrdZ1KFu/buvl8fEX//RCyWh+8j8rCOphxl9oGP2J2yC4CCa6VJXnKQtN/SgBRALUQqVYK4OPD29nigx5eW58JeyWJVIe7EXlfzeZaFfgvg5S1QsTHkZ6jbfj9rDaeitE5nVQaDgQFNAtj6Vjtebl0FR6OBn+NT6Lwjm9mtn6WMo8s9lply15Wmzs7OTJkyhZo1axITE8X58+pakMWL1TfhxMRAYmKiqFmzJlOmTCmSIgSgToU6fNnjC85+dJb8S/k4lbpC4GuBVBlfBZeAIMAT8CQpyfPm2pEHUdiFt9ZaoCuELciIiCgZistqPosFDiyCLZMh+49dJiH9oev74Gvb6QQtHE/JYNLqWHaeUL/X8mmXSP3lS04f233zPoUZEbnBmmfNWMOtBwIanAyU6V6Gsr3LYnQxolgUrm65SsqKFCzZFkAdMXoQMiIi7JVMzQhxq+K4mi/nmrp2JOYrUCzg6AZt/w0tR4CTq9bprEpRFDb8foGJX27hkmdp9cYzB5j58+eMv3K+0GtE7FFWVhbbt0PPnjdvwckvgApPV8CnuQ8ApnRfUn56nUXjetOrp9cDPb4cAifslRQionizZu/m4r6aL/mIurvm7E71ulQwdJ8OtbprGsvaEhMTadfiMa5WbYtv2D/A0RlHs4kn9kWydOf3nM7Pue+uGXt1e7GQhTodAx51tuL/3Ee4VlJXxNQtXY+3wyZQv2z9B3r8G3U23F6M2HOdLfRPChFRfFl7CkUPY9eKovYd2fQOZPwx2V+jG3SPgNLVtM1mBX/tI7LYqyxfNujPzzVaAOCXfZ1re3/kzO5VRbrQ1Jr+HJT7sxCBTAyOzpTu9D2Vn/6EPNQGZwNqDOCNxm9Q2q30Az2+HAIn7IksVhXF043f1n8dvUhKUm9fvvzBH1MPq/kMBqj/BAyPgcdGgtEJjm+ET1rAlimQX7gOnfn5+Y/0eVu5daFpVEwMYUcP8tWITnxTJZsq7gauuvuitB1C8MtzeW38+8WuCAG1GPjpJ6hY8fbbA/yd+PyV59nwzzX0rdYXgOXHl9NnRR++i/8Ok8VU6Mc/c4bbFuiePi1FiCgeimREZMGCBcycOZPk5GQaNmzIvHnzCA29f+MmGREpQWw1haKHEZG/unxcna45+Yt67V1JXcxa9/G7Hqa3ZMkS3nvvPX7++WcqVgz828zXhQuJdO7cmSlTpvCvf/2rCL+ZP91pIWmeyczXv55m3pYT5BSYMRrgmbAg3upaC1/34leQ3K+z6sFLB5kaPZX4q/EA1ChVgwmhE2hWoZkWcYV4aHY1NbNkyRKef/55PvvsM8LCwpgzZw4//vgjCQkJlCtX7p5fK4VICWKrgkGvq/kUBY6uhY3hcP2celtwG+g5E8rVue2u9+uxUaFCIopS9D02HsTFtBymrjvK6kMXACjl7sRb3WrxVPMgHIz6OsnYbDGz7Pgy5h6YS9ofje56VOnBv5v+m/Ie+mt0J/TJrgqRsLAwmjdvzvz58wGwWCwEBgYyYsQIxo8ff8+vlUKkBPn+e3jmmfvfb/FiePrpB3tsPa/mK8iBnXPh19lgygWDA4S9Cu3Hg6vPzbvZU9fRR7H75BUmRcaSkJIBQL1K3kzuW4+mlUtpnMz6rudeZ96Befx47EcUFNwc3Xi1was8H/I8TjprdGcPrLlGXtjRGpH8/Hz27dtH586d/3xCo5HOnTuze/fuv90/Ly+P9PT02z5ECeHvb9373erGBH2lvxzNHhBQvIsQACc3aD8Ohu2B2r1BMcNvn8C8pnBwsdqXBKhYMRCDIQq1CDmFWnzs4paTXDAao6hY0X6LEICW1Uqz9o3WTOwTgperI0eS0vnHp7v499JDXMr4+2m3xZmvqy/vtnyXH3r/QMOyDckx5TBn/xwGRA7g16RftY6nK8uXqwOnHTqofw916KBeP8yyNPHgbDoicuHCBSpVqsSuXbto2bLlzdvHjh3Ltm3biI6Ovu3+kyZNYvLkyX97HBkRKQGKYgqlJPzJc2ILrB8HV46r1wHNoedMoo411vRk2hus2XDscmYeMzYcZeledV2Rl4sjIzvXYFCrYJwc9LUO36JYWHNqDbP3zuZKrtr8rUNgB8Y2H0uAV4DG6Yq34thmqDiwmxGRBxUeHk5aWtrNj8TERK0jiaLi4KBu0YW/L7i8cT1nzqMVDg4O6rvs00+r/9VbEQJQvRMM3QVdpoCzJ5yPgS86EHBgJKXdrqBOx9y9ibotNw8tWbKE+vXrq/+u79BuPzExkfr167NkyZJCPV4ZTxdmPNGQFa+3okGADxl5Jt5fG0/Pj3ew68Rl230jGjAajPSt1pfVj6/m+ZDncTA4sDVxK/1W9mPBwQXkmHK0jlgsmc3qtuc7/e1z47ZRo+znNAi9sumISH5+Pu7u7vz000/079//5u2DBg3i+vXrrFq16p5fL2tESiBpiGA96Rdh83tweCkAV3N8eeeX4Xy+71Msyulb7mj7EZG/9gqJMhgITE6++fnEChVoryicSkl5qAWzFovC0r2JzNiYwNUsdRtyr/r+vN2rDhV93az+/Wjt5PWTROyJIPqiOqpc0aMiY5qPoVNQp9vayot70+OmOnthd4tVQ0NDmTdvHqAuVg0KCmL48OGyWFXcWUmYQilKZ3ehrB2D4dIRAA5cNDN8vR+7En8ABnJjjYi/fxSJiYE2e6kTExNp37w5p1JS7rJclkfunpqWXcDszQks+u0sFgXcnBwY1qEaL7epiquTvn6GFEVh89nNzNw7k+Qstahr6d+S8WHjqepTVeN0xYMt18iXdHY1NTN69Gi+/PJLvv32W+Lj4xk6dChZWVm8+OKLtn5qUVyVhCmUolS5Fed7/Y/wHR5cy1Fo7O/AzpfS+G//hVTw/JEbC1gtlvZcuGC76dDAihWJMhjusVwWooxGAv/a9esB+Lg7MblfPdaMaEPz4FLkFJiZtekY3eZs55ejKY/+TdgRg8FA1+CurOq3ilcavIKz0ZndF3fzj1X/4MO9H5KZn6l1RLtnyzXyovCKpKHZ/PnzbzY0a9SoEXPnziUsLOy+XycjIkI8ulunRUIqBxMe1pFnai/DaFDIyPPkowOv8uXvCzl/8bRt+4j8MQ5+7+WyWG0cXFEUIg9d4IO18VzKyAOgY+1yvNc7hOAyHo/8+PYmMT2RGTEziDofBUAZtzKMbjqa3lV7y3TNXei1zZA9sKupmUchhYgQ1vHXzqoH1+6jctxYyuTuBaDApwqvrLxO99en266z6i3j4LuAx2751E6g1Y0LK4+DZ+aZmPfLcf7v19MUmBWcHYwMaVuFYR2q4+7saLXnsRfbz29n+p7pnMtQG901KdeE8LBwavvV1jiZfdJzmyEtSSEihPibv22NtVjg0GLYPBGy/9hlUqcPdJsKvkHWD1DEIyJ/deJSJpNXx7LjuPq9+vu48k6vEHrWr6C7EYPMnEx+OPEDX/z+BTmmHIwGI0/WfJIRjUfg4+LzQNukSwJZI299UogIIQov5zpETYM9X6gN0RzdoM1oaPUGOLla73nMZhIDAmifnHxzTcgibl0uC1H+/gQmJtpsHFxRFDbFpTBldRxJ19Utry2rlmZyv7rULO9lk+csareOfhl9nJiw8UP2Z28AwNfFl4HBA/nwhQ+ZMlm7c4XskayRty4pRIQQDy4lFtaNhbN/dO30rQzdp0GtHnc9TO9BFMWumcLKLTDz2baTfBp1kjyTBQejgRdaBTOycw28XYtv+/S7nSvkUTsG/2en4hp4AoCc0zk4bnPk942/y8iIsAkpRIQQD0dRIHY5bHwHMtQD5qjeGbpPhzLVH/ph79tHxN+f9hbLQ/cReViJV7P5z5o4NsWpO2rKeLowvkdtBjSuhLGYHqZ313OFjGfw69iN8gMccHBX/9TvV60fo5qOooxbGQ0TCz2SQkQI8WjyMmHHLNg1HywFYHSCVsOhzVvg4vlQD3nrlEFgxYp/GwdPvHCBzp07M2VK0U8ZbDuWyuTIWE5dzgKgSZAvU/rVo14ln/t8pf0xmyEgIJHk5PZwh0kwB++qVHm+H67NNgPg6eTJ641e56naT+FkLL6jQcK+SCEihLCOyydgw3g4ob5p4VURuv4H6v3joaZrHuWsGWueU3PHrzdZ+L+dp5m75TjZ+WYMBng6NIgxXWtRyqP4TF/82S303suC/2/DITYUTCXuShwA1X2rEx4aTqh/aJHmFfpkVw3NhBDFWJnq8OyP8PQP6pqRjAuwbDAs7K2uKbmDrKwsDAYDBoOBrKys2z53v0Lhbp+39jk1d3xuRyOvtavGL/9uT79GFVEUWBx9jvazolj021nMFrv9m+02f54XdO9zhVyvNmRxz8VMbDkRXxdfTlw/weBNg3lr21s3O7UKURRkREQIUTgFubBrLuyYDaYcMDhA6BBoHw5uvjfvlpWejqePOqWRuW4dHl27PtL2A1ufU3M30aeuMDEylqPJGQCE+HszpV9dmgX7PfJj21JhR0Ru3SWdlpfG/APzWXpsKRbFgpujG0PqD2FQ3UE4OxSf0SBhP2RqRghhO9fPwca3IT5SvXYvA50nQaNnYeVKskaMwPOCutA1E/AICFBPVn6Ehgxa7bgxmS18F32ODzclkJ5rAmBA40qM71Gbct5W3NpsRfdbI3Kvc4USriYwNXoq+y/tByDIK4hxoeNoG9C2KL8FoQNSiAhhLdJc4O5OboX1Y+HyMfXarSrMO0xWkpkby1kzAQ9rtKjUuAfJlcw8Zm5MYMneRBQFPF0cGdW5BoNaBePkYF8z3HfdNXNL2Va+fFViYqLuWLQpisLa02uZvXc2qTmpALQLaMe45uMI9LbttmqhH1KICGENd2q3aIW/7nXFlE/er/Nw3jUbQ34miqKQe6CAoC15XM5WSAFunupSqRLExeHxMP+WNe7KesOhxOu8FxnLocTrAFQv58mkPnVpXcM+tr/erY/IDf7+iVgsapFyv2msrIIsPj/0OYviFmFSTDgZnXih7gu8XP9l3J3ci+pbEsWULFYV4lHdOIDi1iIE1NOxnnhC/bwAR2dcO7xFxYiL/PdQPgaDAbcmzhwb7smw5k5UNIAnf3wkJd1cO/LA/liBee/ll9y6UtMmGgb6smJoK2b8owGlPZw5cSmT576OZuj/9nH+WrZNn7swnJ2dmTJlCjVr1iQmJorz59W1IIsXqzVaYmIgMTFR1KxZkylTptxzLY2Hkwejm41mWb9ltPRvSYGlgC8Pf0m/Vf3YdGYTdvw3rChmZEREiL+6cSTnX4uQG+RIztvcek5Lq0AH5vdwpbG/+rocSjYzfH0uv54z37zPQ/3KsZMRkVul5RTw0eZjN3fUuDoZGda+OkPaVsXVSdufC2tvdVYUhV/O/cKMmBlcyFLX/4T5hxEeGk4132qPnFfoj0zNCPEo/tx2cG9F+KZnz7KysmD7dujZU702wMSmTnzQ0RU/N7VIMf1eQP7PuShL1uLRo8eDP4kdnFNzN0eT05m4Kpbo01cBCPJz573eIXSqU053h+nlmHL45sg3fH34a/It+TgaHHm6ztMMbTgUL2d9nNUjrEOmZkqaO/RVEI+gsMP7Np4GKC48PDzw6NoVj4AAPAwGPBT4bG8BNedlUrAvHxQFxwZOuI/wwsPzKJjyH/g5Ei9cULfo8ucISKs//lsVtRhpb7GQ+MdunaJUu4I3P7zSgrlPN6a8twvnrmbz8n/38uLCGE5fzrr/AxQjbo5uvN7odVb1X0XHwI6YFBOL4hbRZ0UfIk9GYlEsWkcUxZCMiBR3sqDS+mRE5OH8sa4mS1Fu3zVT0QF6uELAHyMVpatDj+nqGTaFYK/n1NxJVp6J+VtP8NWOUxSYFZwdjAxuU4XhHarj4eKoSSZb2pm0k2l7pnEm/QwADcs2ZELYBEJKh2gbTGhOpmZKihsLKv/6f6E1tkuWZDfWiCQl/f21BVkjci/Ll/+9j0hgIHw0G6rlw+b3IEvdEkrt3tDtAygVfN+Htedzau7kVGomk1fHse2Y+r1W8Hbl7V516N3AX3fTNQXmAhbFL+KzQ5+RY8rBgIEnaz7JiMYj8HX11Tqe0IgUIiWBLKi0rRtFHtxejEiRd1/37KyamwZR0yH6M1DM4OgKrd+Ex0aCk9s9H9fWZ81Ym6Io/Bx/iSlrYkm8mgNAi6p+TO5bj1oV9LeeIiUrhdn7ZrPu9DoAfFx8eKPxG/yjxj9wMMrvoJJGCpGSQKYPbO9O016BgTBnzu1FiDQ9e3CX4tVmaKe3q9e+QdAtAmr3eqjD9OxZboGZL7afYsHWE+SZLDgYDTzfsjKjOtfEx01/p93uTd5LxJ4Ijl1TG93V8avDhLAJNCrXSNtgokhJIVISfP89PPPM/e+3eDE8/bTt8+jV/YoMWaPz8BQF4laq7eLTk9TbqnWEHjOgTA1No9nC+WvZvL8mng2x6vqW0h7OjOtRmyeaBGA06qv4MllMLE1YyvyD88nIV8/q6VutL282fZMybvbR/E3YlhQi96OHv2BlRER7Wq3R0cPP763ys9SD9HbNBXM+GJ2gxVBoNxZc9DeFseN4KpMiYzmZqu6oaRToy5R+dWkQ4KttMBu4mnuVufvnsvz4chQUPJw8GNpwKM/UeQYno/5Gg8SfpBC5F738BSsLKrWl1Rodvfz83smVk7BxAhzboF57VoCu70P9J3Q3XZNvsvDtrjPM+fkYWflmDAZ4qnkgY7rVxs/Dfta5WMuRy0eYGj2Vw5cPA1DVpyrhYeG08G+hcTJhK1KI3I3edpnIgkrtaDEipbef37tJ2AAbxsO10+p1UCvoORMq1NM2lw1cSs8lYv1RVhxQp6a8XR15q1stngkNwtHODtN7VBbFwqoTq5izfw5Xc9Xmb10qd2FMszH4e/prnE5YmxQid6LXXSaFXVAprKuo1+jo9ef3bgpyYfd82PEhFGSDwQjNX4YOE8CtlNbprC7mzFXeWxVL/MV0AGpX8GJKv3qEVvHTOJn1peen88nBT/j+6PdYFAuuDq68XP9lXqj3Ai4OLlrHE1Yihcid6HlNhd7WDBQHRf3zpOef33tJOw+b3oHYFeq1e2noNBEaDwSjvkYMzBaFxdFnmbXpGGk5BQD0b1SR8J51KO/tqnE660u4mkDEngj2pewDIMAzgHGh42gX0E53vVZKImnxfid6btvt4KC++Tz9tPpfKUJsr00bdQTibr8wDQZ1ZKpNG+s8n55/fu/FJwCeXAjPR0LZOpB9BVa/AV91gvP7tE5nVQ5GAwNbBrP1rfY8HRqEwQArD16g46woPt92knyTvtqn1/KrxTfdvmF6m+mUcyvH+czzjPhlBMO2DONs+lmt44kiVHIKEf9CzkEW9n6iZHNwUBeIwt+LkRvXc+ZYrygs6T+/VdvBazvUXiMu3nBhP3zVEVYNg8xUrdNZlZ+HMxED6hM5rDWNg3zJyjcTsf4o3T/ezvZj+vpeDQYDPav2ZPXjq3mp3ks4Gh3ZkbSDx1c9zsf7Pya7IFvriKIIlJypGdllImyhqNboyM/vnzIvwc+T4OB36rWLj7p2pPnL4KCv81wsFoVl+88zfcNRLmeqhwV2q1ued3qFEOjnrnE66zuTdoZpMdPYmbQTgHLu5RjTbAzdgrvJdE0xI2tE7kZ2mQhbKKo1OvLze7vEGFj3Flw8qF6Xqws9Z0Bwa01j2UJ6bgEfbT7Gf3efxWxRcHE0MrR9NV5rVw1XJ30VnoqiEJUYxfSY6SRlqruJmldoTnhoODVK6a/RnV5JIXIvsstEFGfy83s7ixn2/xe2TIEcdUso9f4BXf4DPpW0zWYDCckZTIqMZfepKwAElHLj3d4hdA0pr7sRg1xTLt/EfsPXh78mz5yHg8GBp2s/zdBGQ/F2lk7b9k4KkfuRXSaiOJOf37/LvgpbP4C9/weKBZw8oO1b0HIYOOprS6iiKKw9fJEP1sZzMS0XgLY1yzKxTwjVynpqnM76kjKTmBUzi5/P/QyAn6sfo5qMol/1fhgNJWeZY3EjhYgQomS6+DusGwOJv6nXftWgx3So0UXbXDaQnW9iwdYTfLn9NPlmC04OBl5qXYURHWvg6aKvtTIAu5J2EbEngjPpZwBoUKYBE8ImULdMXW2DiTuyi0Lkgw8+YO3atRw8eBBnZ2euX7/+wI8hhYgQ4oEpCvy+FDa/C5kp6m01e0D3qeBXVdtsNnDmchZT1sTxy9FLAJT3dmFCzzr0bVhRd9M1BeYCvov/jk8PfUq2KRsDBgbUGMDIJiMp5aq/RnfFmV0UIhMnTsTX15fz58/z9ddfSyEitCdTGiVLbjpsnwG/fQoWEzi4wGMjofWb4Ky/HSdb4lOYsiaOs1fULa+hwX5M7leXOv76+915KfsSH+37iDWn1gDg7ezN8MbDebLmkzga9TcaVBzZRSFyw8KFCxk1apQUIkJbej4sTtxbagKsHwunotRrn0D1ML2Qfro7TC+3wMzXv55m3i/HyS2wYDTAwBaVGd2lFj7u+jvtdl/KPiKiI0i4lgBArVK1mBA2gSblm2icTBTbQiQvL4+8vLyb1+np6QQGBkohIh5NSTksTtydokD8avV037RE9bYq7aDHDChXW9tsNpB0PYepa+NZe1jttOvn4czYbrX4Z7NAjEZ9FV8mi4mfjv3E3ANzycjPAKBX1V6Mbjqacu7lNE5XchXbQmTSpElMnjz5b7dLISIe2qMeFifTOfqSnw0758Cvc8CcB0ZHCHsN2o0DV/39jtl14jITI2M5fikTgIYBPkzuV49Ggb7aBrOBq7lXmbt/LsuPL0dBwd3RnaENh/JsnWdxctDfaJC9s9lZM+PHj8dgMNzz4+jRow8dPDw8nLS0tJsfiYmJD/1YQgBqEXG3IgTUv5QTE9X7/dXy5WoR06GDetJuhw7q9fLltkorbM3ZXe3COiwaavVU147sng/zm8GhH+7ctbYYa1W9DOtGtuGdXnXwcnHk0Pk0+i/YydifDnE5M+/+D1CM+Ln6ManVJL7v9T0NyjQg25TNh/s+ZEDkAHYl7dI6nriHBxoRSU1N5cqVK/e8T9WqVXF2dr55LWtEhKa+/14tIu5n8WL10MAbZDqnZDi+GdaPg6sn1evAMOg5E/wbapvLBi5l5DJ9fQLL9quFuZerI6O71GRgi8o4OuirH4dFsbDqxCrm7J/D1Vy10V2noE6MaT6GSp76a3Rnj4rt1MxfSSEiHllUlDqScT9bt6onF8OjT+eI4sWUB799AttmQkEWYIBmL0LHd8HdT+t0Vrfv7DUmRh7hSFI6ALXKezGpb11aViutcTLrS89P59ODn/L90e8xK2ZcHFwYXH8wL9Z9EVdHV63j6ZpdFCLnzp3j6tWrREZGMnPmTHb8MfRdvXp1PD0L1/1PChHxyB7msLiHKV5E8ZeWpPYeObJMvXYrBZ3egyaDwKivgtNsUVgSk8iMjUe5nl0AQJ+GFZnQszb+Pm4ap7O+49eOE7EngpjkGAAqeVZibPOxdAjsoLteK/bCLgqRF154gW+//fZvt2/dupX2hfzlLYWIsIoHPSzuYadzipIsorWd0zvU7b6X4tRr/4bQcxYEhmqbywauZeXz4eYEFkefw6KAu7MDwztWZ3DrKrg46uvnSVEUNp7dyKyYWaRkq43uHqv4GONCx1HFp4rG6fTHLgoRa5BCRFjNgxwWZ+8jItITxfbMJoj5CrZOhbw09baGz0DnSeBVXtNotnAkKY2JkbHsO3sNgCplPHivTwgdaulv+2t2QTZfHf6KhbELKbAU4Gh05PmQ53m1wau4O+mv0Z1WpBAR4k4KO4rwMNM5RUUW0RatzFTYMgkO/E+9dvGG9uMh9BXQ2ZZQRVFYcSCJiPVHSc1Qd9R0rlOe93qHEFRaf2/Q59LPMW3PNHYkqcsGyrmVY3Sz0fSs0lOma6xAChFRMllzuuJBp3OKgiyi1c75vbDuLbhwQL0uW1tthla1nba5bCAjt4C5W47zzc4zmCwKzo5GXmtblaHtq+PmrL+fq22J25geM53EDLVdRNPyTQkPDaeWXy2NkxVvUoiIkscW0xUPMp1TFOx9ykjvLBY4sAi2TIbsP9oYhPRX28X7BmoazRZOXMpgUmQcv564DEAlXzfe7V2HbnUr6G7EIM+cx7ex3/Ll71+Sa87FaDDyVK2neL3R6/i4+Ggdr1iSQkSULLacrrCnRaHFYRFtSZBzTV07EvMVKBZwcoc2o6HlCHDS15ZQRVHYcCSZ99fGk3Q9B4DW1cswqW8I1ct5aZzO+i5mXmTW3llsOrsJgFIupRjZZCSP13gco0FfvVZsTQoRUXKUpOkKGRGxL8mHYd0YOLdbvS5VBbpPg1rdtc1lAzn5Zj6NOsFn20+Rb7LgaDTw4mPBvNGpBl6u+lorAxB9MZqI6AhOpqmN7uqVrseEsAnUL1tf42TFhxQiouQoSW/O9ryItqRSFDj8E2x6BzKT1dtqdIPuEVC6mrbZbODclWymrInj53h1+2tZLxcm9KxN/0aVdDddU2Ap4Pv47/nk0CdkFWQBMKDGAN5o/Aal3fTX/M3abHbWjBB25+JF697Pnjk4qGte4O/H19+4njNHipCiZDBAgydhxF54bCQYneD4RvikBWyZAvlZWie0qqDS7nw1qBnfvNCc4NLupGbk8eaSQzz52W5iL6RpHc+qnIxOPF/3edY8voa+1foCsPz4cvqs6MN38d9hspg0TqgfMiIiireSNCJyg70tohV/unxcbYZ28hf12ruSupi17uN/Lx6LuTyTma92nGb+LyfIKTBjNMCzYZX5d9ea+Lo73/8BipmDlw4yNXoq8VfjAahRqgbhoeE0r9Bc42T2SaZmRMlRUqcr7GkRrbidosDRtbAxHK6fU28LbqMepleujrbZbODC9Rymrotnze/qqGMpdyfGdKvNv5oH4mDUV/FltphZdnwZcw/MJe2PRnc9gnswutloKnhU0DidfZFCRJQs9tjzQ4iCHNj5Mfz6EZhyweAAYa+qDdFc9bcldPfJK0yKjCUhJQOA+pV8mNS3Lk0rl9I4mfVdz73OvAPz+PHYjygouDm68WqDVxkYMhBnB/2NBj0MKUREySPTFcJeXTsLGyfA0TXqtUdZ6DwZGj4NRn0t0zOZLSz67SyzNx0jI09dQ/GPJgGM71Gbsl4uGqezvrgrcURER3Aw9SAAlb0rMz50PK0rtdY2mB2QQkSUTDJdIezZiS2wfhxcOa5eBzRXp2sqNtY2lw2kZuQxY8NRftyn/mHg5eLIqC41eb5lZZwc9FV8KYrCmlNr+HDvh1zJVRvddQjswJjmYwj00l+ju8KSQkQIIeyRKR+iP4VtMyA/EzBA00HQ8T3w0N+W0P3nrjFxVSyHk9T1FDXLezKpT11aVS+jcTLry8zP5NNDn7I4fjEmxYSz0ZmX6r/ES/Vews3RTet4RU4KESGEsGfpF2Hze3B4qXrt6gsd34FmL4FRX6N4ZovC0r2JzNhwlGvZBQD0qu/P273qUNFXf2/QJ6+fJGJPBNEXowHw9/BnbPOxdArqpLteK/cihYgQQhQHZ3aq231TjqjXFepDz1kQ1ELbXDZwPTuf2ZuP8b/fzmJRwM3JgeEdq/Nymyq4OOqr+FIUhc1nNzNz70ySs9RGdy39WzI+bDxVfapqnK5oSCEihBDFhdkE+76BX/4DuX80BWvwFHSZDF762xIadyGdiZFHiDlzDYDKpd2Z2CeEjrXLa5zM+nJMOXx1+CsWHllIviUfR4Mjz4U8x2sNX8PDyUPreDYlhYgQQhQ3WZfVbqz7/wso4OwF7cZC2GvgqK8toYqisOrgBaaui+dSRh4AHWuX473eIQSX0d8bdGJ6IjNiZhB1PgqAsm5lGd1sNL2q9NLtdI0UIkIIUVwl7VMP00vap16XqQk9pkO1jtrmsoHMPBPzthzn619PY7IoODsYeaVtVV7vUA13Z0et41nd9vPbmb5nOucy1EZ3Tco1ITwsnNp+tTVOZn1SiAghRHFmscChxbB5ImRfVm+r0xe6fQC+Qdpms4ETlzKZvDqWHcfV77Wijytv9wqhZ/0KuhsxyDfn89+4//LF71+QY8rBaDDyZM0nGdF4BD4u+ml0J4WIECWJ9E/Rr5zrEBUBe74ExQyObtBmNLR6A5xctU5nVYqisDE2hf+siSPpeg4ALauWZnK/utQs76VxOutLzkrmw70fsuHMBgB8XXwZ2WQkj1d/HAcd7JySQkSIkuJOHWUDAtRTeqWjrH6kxMK6sXD2V/W6VDB0i4BaPXR3mF5ugZlPo07y2baT5JksOBgNvNAqmJGda+Dt6qR1PKuLSY5havRUTlw/AUBI6RAmhE2gYdmGGid7NFKICFES3Dhj56//hOWMHX1SFDiyDDa9CxkX1Nuqd1HXj5Supm02G0i8ms1/1sSxKS4FgDKeLozvUZsBjSth1NlheiaLiR+O/sCCgwvILMgEoF+1foxqOooybsWz+ZsUIkLo3Y1Th28dCbmVXk8dFpCXCTtmwa75YCkAB2doOQzavAUunlqns7ptx1KZHBnLqctZADQJ8mVKv3rUq6Sf9RQ3XM65zMf7P2bliZUAeDp5MqzRMP5V+184GYvXaJAUIkLoXVQUdOhw//tt3Qrt29s6jdDC5ROwYTyc2Kxee1WErv+Bev/Q3XRNvsnC/+08zdwtx8nON2MwwNOhQYzpWotSHvra2gxwKPUQU6OnEnclDoDqvtUJDw0n1D9U42SF9yDv3/o6fUiIkuLiReveTxQ/ZarDsz/C0z+Ab2V1umbZYFjYW11ToiPOjkZea1eNX/7dnn6NKqIosDj6HO1nRbHot7OYLXb79/RDaVi2IYt7LmZiy4n4uvhy4voJBm8azFvb3rrZqVVPZEREiOJIRkTErQpyYddc2DEbTDlgcIDQIdA+HNx8tU5nddGnrjAxMpajyRkAhPh7M6VfXZoF+2mczPrS8tKYf2A+S48txaJYcHN0Y0j9IQyqOwhnB/sdDZKpGb2T7ZrixhqRpKS/L1YFWSNSUl0/BxvfhvhI9dq9jNoqvuEzYNTXALjJbOG76HN8uCmB9FwTAAMaV2J8j9qU89bX1maAhKsJTI2eyv5L+wEI8gpiXOg42ga01TjZnUkhomeyXVPccGPXDNxejMiuGXFyK6wfB5cT1OtKzaDnDKjUVNtcNnA5M4+ZGxJYui8RRQFPF0dGda7BoFbBODnoq/hSFIW1p9cye+9sUnNSAWgX0I5xzccR6B2ocbrbSSGiV7JdU/zVnQrTwECYM0d+Fko6cwFEfw5R0yA/AzBAk4HQaSJ4FM8tofdyMPE6E1cd4dB59eDA6uU8mdSnLq1r6O97zSrI4vNDn7MobhEmxYST0YkX6r7Ay/Vfxt3JXet4gBQi+iTbNcXdyFSduJeMZLVV/O8/qNeuPtDxXWj6Ijjo6zwXi0Xhx32JzNiQwJWsfAB61KvA273qEFDKPt6grelU2immRU9j98XdAFTwqMCYZmPoUrmL5q3xpRDRI1mcKIR4FOeiYd1bkPy7el2+HvScCZVbaZvLBtJyCvho8zH+u/sMFgVcnYwMa1+dIW2r4uqkryJdURS2nNvCzJiZXMhSG92FVQgjPCycar7aNbqT7bt6JNs1hRCPIigMXomCXrPBrRSkHIFvesCylyFdX783fNycmNS3LmvfaENoFT9yCyx8uPkYXT/azua4FOz47+8HZjAY6Fy5Myv7r2Row6E4G52JTo7micgnmBEzg4z8DK0j3pfNCpEzZ84wePBgqlSpgpubG9WqVWPixInk5+fb6in1zd/fuvcTQpQ8RgdoPhhG7IdmLwEGOPwjzG8Gv84Bk75+P9fx92bJKy2Y+3Rjynu7cO5qNkP+u5cXF8Zw+o9OrXrh5ujG641eZ2X/lXQI7IBJMbEobhF9VvQh8mQkFsWidcS7stnUzIYNG1iyZAlPP/001atX58iRIwwZMoSBAwcya9asQj2GTM3cQrZrCiGs7cJBWDcGzu9Rr0tXV8+uqd5Z01i2kJVnYv7WE3y14xQFZgVnByOD21RheIfqeLjoa60MwK9JvzJtzzTOpp8F1CZpE8ImEFI6pEie327XiMycOZNPP/2UU6dOFer+Uoj8hWzXFEJYm8UCvy+Bze9B1iX1ttq9odsH6im/OnMqNZPJq+PYdkzd/lrB25UJverQp4G/5gs8rS3fnM+iuEV8/vvn5JhyMGDgyZpPMqLxCHxdfW363Ha7RiQtLQ0/P/11visyAwaoxUalSrffHhAgRYgQ4uEYjdDoaRixF1oMU7uyHl0DC8JgawQU5Gid0KqqlvVk4YvN+fL5ZgT6uZGcnssb3x/gqS9+42hyutbxrMrZwZnB9Qezuv9qelTpgYLC0mNL6b2yN0sTlmK2mLWOCBThiMiJEydo2rQps2bNYsiQIXe8T15eHnl5eTev09PTCQwMlBGRv5LtmkIIW7kUD+vHwunt6rVPEHSfqo6S6GzEILfAzBfbT7Fg6wnyTBYcjAYGtqjMm11q4uNWvE67LYyY5Bgi9kRw/NpxAOr41WFC2AQalWtk9eey6dTM+PHjmT59+j3vEx8fT+3atW9eJyUl0a5dO9q3b89XX31116+bNGkSkydP/tvtUoiIEkkKTqEVRYG4VWq7+PQ/ehdV7QA9ZkDZmtpms4Hz17L5YG0864+oB8qV9nBmXPfaPNE0AKNRX8WXyWJiScISFhxYQEaBuqOmb7W+TG41GUej9dbK2LQQSU1N5cqVK/e8T9WqVXF2Vg/juXDhAu3bt6dFixYsXLgQ4z3OO5ARESH+IK38hT3Iz1IP0ts1F8z5YHSEFq9Du7Hg4qV1Oqv79fhlJkYe4WSquqOmYaAvU/rWpWGgr7bBbOBKzhXmHpjL8uPL6VGlBzPazrDq49vNYtWkpCQ6dOhA06ZN+d///ofDA/41J4tVRYkkrfyFvblyEjZOgGMb1GvPCtD1P1D/Sd1N1+SbLHy76wwfbzlOZp4JgwH+1SyQMd1qUdrTRet4Vnc49TDl3MtR3qO8VR/XLgqRpKQk2rdvT+XKlfn2229vK0IqVKhQqMeQQkSUONLKX9izhA2wYTxcO61eB7VSD9OrUF/bXDZwKT2XaeuPsvxAEgDero681a0Wz4QG4aizw/RswS4KkYULF/Liiy/e8XOFfUopRESJI638hb0ryIXd82H7LDDlgMEIzQZDx7fVjq06s/fMVd5bFUvcRXVHTe0KXkzpV4/QKrID9F7sYvvuCy+8gKIod/wQQtzFqlWFu5+08hdacXKFtm/B8BgI6Q+KBWK+hHlNYd9CsJMtodbSLNiP1SNa85/+9fBxc+Jocgb//Hw3I384QEp6rtbxdEEOvRPCXpjNUKECXL58//vKiIiwF6e2qdt9U4+q1xUbQ89ZENBM21w2cDUrn1mbEvh+zzkUBTycHXijUw1efKwKzo4yXXMru5iasQYpRESJUthpmbJl1RERWSMi7IW5APZ8CVERkPdHU7BGz0HnSeBZVtNotnD4fBrvRR7hwLnrAFQt68HEPnVpV1N/3+vDsoupGSHEAyrsdMuzz0oRIuyLgxO0fB2G74WGz6i3HfyfOl3z22dgNmmbz8rqB/iw7LVWzHqyIWU8nTmVmsWg/9vDK//dS+LVbK3jFTtSiAhhLwp7cnK/frbNIcTD8ioPj38KgzeDf0PIS4MN4+DzNnB6h9bprMpoNPBE0wB+eas9g1tXwcFoYFNcCp1nb2POz8fILdDXWhlbkqkZIezF/U5YBggMlK27oniwmGH/t7BlCuRcU2+rOwC6vg8+le79tcXQsZQMJq6KZfcpteFnQCk33u0dQteQ8ro7TK8wZI2IEMWVnLAs9Cb7KvzyPuz7Rt1h4+QObcdAy2HgqK8GYYqisPbwRT5YG8/FNHVHTduaZZnYJ4RqZT01Tle0pBARoji7U3v3wECYM0eKEFF8XTwE68ZAYrR67VcNekyHGl20zWUD2fkmPtl6ki+2nyLfbMHJwcBLraswomMNPF2sd56LPZNCRIjiTg68E3qkKPD7Etj8HmSmqLfV7KGe7utXVdtsNnDmchZT1sTxy9FLAJT3dmFCzzr0bVhR99M1UogIIYSwX7npsG06RH8GFhM4uMBjI6H1m+DsrnU6q9sSn8KUNXGcvaLuqAkN9mNyv7rU8dfv+5oUIkIIIexfaoLaDO1UlHrtE6guZg3pp7vD9HILzHz962nm/XKc3AILRgMMbFGZ0V1q4ePupHU8q5NCRAghRPGgKBAfCRvfhrRE9bYq7aDHDChXW9tsNpB0PYepa+NZe1jtG+Tn4czYbrX4Z7NAjEb9FF9SiAghhChe8rNh5xz4dQ6Y88DoCGGvQbtx4Kq/3/+7TlxmYmQsxy9lAtAwwIfJ/erRKNBX22BWIoWIEEKI4unqadg4ARLWqdee5aHLFGjwL91N1xSYLfx391nmbD5GRp7affafzQIY2702ZTyL99ZmKUSEEEIUb8c3w/pxcPWkeh3YAnrOUDu26syljFymr09g2X51y76XqyOju9RkYIvKODoUzwboUogUBdleKYQQtmXKg98+gW0zoSALDEZo+iJ0fAfc/bROZ3X7zl5jYuQRjiSpBwfWruDFpL51aVG1tMbJHpwUIrZ2p4ZTAQHw8cfScEoIIawtLQk2vwtHlqnXbqWg03vQZBAY9fUHoNmisCQmkRkbj3I9uwCAPg0rMqFnbfx93DROV3hSiNjSjRbcf33ZpAW3EELY1ukd6nbfS3HqtX9D6DkLAkO1zWUD17Ly+XBzAoujz2FRwN3ZgeEdqzO4dRVcHO2/+JJCxFZuHEp260jIrQwGdWREDiUTQgjbMJsg5ivYOlU93Reg4TPQeZJ6+q/OHElKY2JkLPvOqgcHVinjwcQ+IbSvVU7jZPcmhYitREVBhw73v9/WrdC+va3TCCFEyZWZClsmwYH/qdcu3tB+PIS+Ag76ahCmKAorDiQRsf4oqRl5AHSuU573eocQVNo+O9E+yPt38VyOq5WLF617PyGEEA/Hsyz0WwAvb4GKjSEvXd32+1lrOLVN63RWZTAYGNAkgF/+3Y4hbargaDTwc3wKnT/axuxNCeTkm7WO+EikEHkQ/v7WvZ8QQohHE9AMXv4F+swF99KQehT+2xeWDoLriVqnsyovVyfe7hXChlFtaF29DPkmC3N/OUHn2dtYf/gidjzBcU8yNfMgbqwRSUr6+2JVkDUiQgihpZxr6tqRmK9AsYCTO7QZDS1HgJOr1umsSlEUNhxJ5v218SRdzwGgdfUyTOobQvVyXhqnkzUitnVj1wzcXozIrhkhhLAPyYdh3Vg4t0u9LlUFekyHmt20zWUDOflmPok6wefbT5FvsuBoNPDiY8G80akGXq7arZWRQsTW7tRHJDAQ5syRIkQIIeyBosDhn2DTO5CZrN5Woxt0j4DS1bTNZgPnrmQzZU0cP8enAFDWy4UJPWvTv1ElDBq0xpdCpChIZ1UhhLB/eRmwfSbs/gQsBeDgDK1GQJt/g7OH1umsbuvRS0xeHcuZK9kANKtcisn96lK3ok+R5pBCRAghhLjV5eNqM7STv6jX3pWg6/tQ93HdHaaXZzLz9a+nmbflBDkFZowGeDasMv/uWhNfd+ciySCFiBBCCPFXigJH18LGcLh+Tr0tuA30nAnl6mibzQYupuUwdd1RVh+6AEApdyfe6laLp5oH4WC0bfElhYgQQghxNwU5sPNj+PUjMOWCwQHCXlUborkW7RRGUdh98gqTImNJSMkAoF4lbyb3rUfTyqVs9pxSiAghhBD3c+2s2gTt6Br12qMsdJ4MDZ8Go77abJnMFhb9dpbZm46RkWcC4ImmAYzrXpuyXi5Wfz4pRIQQQojCOrEF1o+DK8fV64Dm6nRNxcba5rKB1Iw8Zmw4yo/71F2fXi6OjOpSkxdaBVt1ukZavAshhBCFVb0TDN0FXf4Dzp5wPga+6ACrR0LWFa3TWVVZLxdmPtmQ5a+3okGADxl5JjYcuYiNl4zck4yICCGEEDekX4TN78Hhpeq1qy90fAeavQRGfbVosFgUlu5NpEGALyEVrfseK1MzQgghxKM4uwvWjYGUI+p1hfrQYyZUbqltrmLCbqZm+vbtS1BQEK6urvj7+zNw4EAuXLhgy6cUQgghHl3lVvDKNug5S91Jk3wYvukOy1+BjGSt0+mKTQuRDh06sHTpUhISEli2bBknT57kiRvntAghhBD2zMERQofAiP3QZBBggN+XwLymsHMumPK1TqgLRTo1ExkZSf/+/cnLy8PJ6f6H8cjUjBBCCLuRtF+drknaq16Xqakepleto7a57JDdTM3c6urVq3z33Xe0atXqrkVIXl4e6enpt30IIYQQdqFSExi8GfotAPcycPkYLHocljz3Z6dW8cBsXoiMGzcODw8PSpcuzblz51i1atVd7xsREYGPj8/Nj8DAQFvHE0IIIQrPaITGz8GIfRA2VO3KGr8a5jeHqOlq11bxQB54amb8+PFMnz79nveJj4+ndu3aAFy+fJmrV69y9uxZJk+ejI+PD2vWrLnjscR5eXnk5eXdvE5PTycwMFCmZoQQQtinlFhYNxbO/qpe+1aG7tOgVg/dHab3IGy6fTc1NZUrV+7d4KVq1ao4O//9hL/z588TGBjIrl27aNny/lugZI2IEEIIu6coELscNr4DGX/sDK3eGbpPhzLVtc2mkQd5/3Z80AcvW7YsZcuWfahgFosF4LZRDyGEEKJYMxig3j+gRjfY8SHsmgcnfoZPWkDLYdB2DLh4ap3Sbtls10x0dDQxMTG0bt2aUqVKcfLkSd59911SUlKIjY3FxeX+h+zIiIgQQohi5/IJ2DAeTmxWr70qQtf/qMVKCZmusYtdM+7u7ixfvpxOnTpRq1YtBg8eTIMGDdi2bVuhihAhhBCiWCpTHZ79EZ7+QV0zknEBlg2Ghb3VNSXiNtLiXQghhLCVglzYNRd2zAZTjrrLJnQItA8HN1+t09mMXYyICCGEECWekyu0GwvD90CdvqCYIfoztTvr/kXwx9rJkkwKESGEEMLWfIPgX4tg4Eq1I2v2ZYgcDl93hqR9WqfTlBQiQgghRFGp1gFe2wld3wdnL7UI+bITRI6ArMtap9OEFCJCCCFEUXJ0hlYjYMReaPAUoMD+/8K8JhD9BZhNWicsUlKICCGEEFrwqgADPoeXNkKF+pCbBuvHwBft4OwurdMVGSlEhBBCCC0FtYBXtkGvD8HVF1KOwDc9YNnLkH5B63Q2J4WIEEIIoTWjAzR/GUbsh6YvAgY4/CPMawa/zgFTvtYJbUYKESGEEMJeeJSGPnPgla0Q0BwKsuDnifBpS7VtvA5JISKEEELYm4qN4aVN0P9T8CgLV07A//4BPzwL185onc6qpBARQggh7JHRCI2egRH7oMUwtSvr0TWwIAy2RkBBjtYJrUIKESGEEMKeufpA96kwdCdUaQumXNg2DRaEQvxqsN+TWgpFChEhhBCiOChXB56PhCcXgncAXD8HS56D/w2A1GNap3toUogIIYQQxYXBAHUfV8+uafMWODjDyV/Uxayb3oW8DK0TPjApRIQQQojixtkDOr0Lr/8GNbqBxaSe8juvGfy+tFhN10ghIoQQQhRXpavBs0vh6SVQqgpkJsPyIfBNT0g+rHW6QpFCRAghhCjuanVXR0c6vgtO7nBuF3zeFta+BdlXtU53T1KICCGEEHrg5Apt34LhMeo6EsUCMV/C/GawbyFYzFonvCMpRIQQQgg98QlQd9Y8Hwlla0P2FVg9Er7qBOf3ap3ub6QQEUIIIfSoajt47VfoFgEu3nDhgFqMrBwGmalap7tJChEhhBBCrxycoOXrMHwvNHxGve3g/2BeU/jtMzCbtM2HFCJCCCGE/nmVh8c/hcGbwb8h5KXBhnHweRs4vUPTaFKICCGEECVFYCgM2Qq954BbKbgUp64f0XBkxFGzZxZCCCFE0TM6QLMXIaQfbP0AavYAB+3KASlEhBBCiJLI3Q96fah1CpmaEUIIIYR2pBARQgghhGakEBFCCCGEZqQQEUIIIYRmpBARQgghhGakEBFCCCGEZqQQEUIIIYRmpBARQgghhGaKpBDJy8ujUaNGGAwGDh48WBRPKYQQQohioEgKkbFjx1KxYsWieCohhBBCFCM2L0TWr1/Ppk2bmDVrlq2fSgghhBDFjE3PmklJSWHIkCGsXLkSd3f3+94/Ly+PvLy8m9fp6em2jCeEEEIIjdlsRERRFF544QVee+01mjVrVqiviYiIwMfH5+ZHYGCgreIJIYQQwg488IjI+PHjmT59+j3vEx8fz6ZNm8jIyCA8PLzQjx0eHs7o0aNvXqelpREUFCQjI0IIIUQxcuN9W1GU+97XoBTmXrdITU3lypUr97xP1apV+ec//8nq1asxGAw3bzebzTg4OPDss8/y7bff3ve5zp8/L6MiQgghRDGVmJhIQEDAPe/zwIVIYZ07d+62kYwLFy7QrVs3fvrpJ8LCwu4bDMBisXDhwgW8vLxuK2isKT09ncDAQBITE/H29rbJc+iFvFaFJ69V4cjrVHjyWhWevFaFZ6vXSlEUMjIyqFixIkbjvVeB2GyxalBQ0G3Xnp6eAFSrVq1QRQiA0Wgs9H0flbe3t/zAFpK8VoUnr1XhyOtUePJaFZ68VoVni9fKx8enUPeTzqpCCCGE0IxNt+/eKjg4uFCLVoQQQghRcpT4EREXFxcmTpyIi4uL1lHsnrxWhSevVeHI61R48loVnrxWhWcPr5XNFqsKIYQQQtxPiR8REUIIIYR2pBARQgghhGakEBFCCCGEZqQQEUIIIYRmpBC5Rd++fQkKCsLV1RV/f38GDhzIhQsXtI5ld86cOcPgwYOpUqUKbm5uVKtWjYkTJ5Kfn691NLv0wQcf0KpVK9zd3fH19dU6jl1ZsGABwcHBuLq6EhYWxp49e7SOZHe2b99Onz59qFixIgaDgZUrV2odyW5FRETQvHlzvLy8KFeuHP379ychIUHrWHbn008/pUGDBjebmLVs2ZL169drlkcKkVt06NCBpUuXkpCQwLJlyzh58iRPPPGE1rHsztGjR7FYLHz++efExsby0Ucf8dlnnzFhwgSto9ml/Px8nnzySYYOHap1FLuyZMkSRo8ezcSJE9m/fz8NGzakW7duXLp0SetodiUrK4uGDRuyYMECraPYvW3btjFs2DB+++03Nm/eTEFBAV27diUrK0vraHYlICCAadOmsW/fPvbu3UvHjh3p168fsbGx2gRSxF2tWrVKMRgMSn5+vtZR7N6MGTOUKlWqaB3Drn3zzTeKj4+P1jHsRmhoqDJs2LCb12azWalYsaISERGhYSr7BigrVqzQOkaxcenSJQVQtm3bpnUUu1eqVCnlq6++0uS5ZUTkLq5evcp3331Hq1atcHJy0jqO3UtLS8PPz0/rGKKYyM/PZ9++fXTu3PnmbUajkc6dO7N7924Nkwk9SUtLA5DfTfdgNpv54YcfyMrKomXLlppkkELkL8aNG4eHhwelS5fm3LlzrFq1SutIdu/EiRPMmzePV199Vesoopi4fPkyZrOZ8uXL33Z7+fLlSU5O1iiV0BOLxcKoUaN47LHHqFevntZx7M7hw4fx9PTExcWF1157jRUrVhASEqJJFt0XIuPHj8dgMNzz4+jRozfvP2bMGA4cOMCmTZtwcHDg+eefLzFn5DzoawWQlJRE9+7defLJJxkyZIhGyYvew7xWQoiiM2zYMI4cOcIPP/ygdRS7VKtWLQ4ePEh0dDRDhw5l0KBBxMXFaZJF9y3eU1NTuXLlyj3vU7VqVZydnf92+/nz5wkMDGTXrl2aDVkVpQd9rS5cuED79u1p0aIFCxcuxGjUfV1708P8XC1cuJBRo0Zx/fp1G6ezf/n5+bi7u/PTTz/Rv3//m7cPGjSI69evy0jkXRgMBlasWHHbayb+bvjw4axatYrt27dTpUoVreMUC507d6ZatWp8/vnnRf7cRXb6rlbKli1L2bJlH+prLRYLAHl5edaMZLce5LVKSkqiQ4cONG3alG+++aZEFSHwaD9XApydnWnatClbtmy5+aZqsVjYsmULw4cP1zacKLYURWHEiBGsWLGCqKgoKUIegMVi0ey9TveFSGFFR0cTExND69atKVWqFCdPnuTdd9+lWrVqJWI05EEkJSXRvn17KleuzKxZs0hNTb35uQoVKmiYzD6dO3eOq1evcu7cOcxmMwcPHgSgevXqeHp6ahtOQ6NHj2bQoEE0a9aM0NBQ5syZQ1ZWFi+++KLW0exKZmYmJ06cuHl9+vRpDh48iJ+fH0FBQRomsz/Dhg1j8eLFrFq1Ci8vr5vrjXx8fHBzc9M4nf0IDw+nR48eBAUFkZGRweLFi4mKimLjxo3aBNJkr44d+v3335UOHToofn5+iouLixIcHKy89tpryvnz57WOZne++eYbBbjjh/i7QYMG3fG12rp1q9bRNDdv3jwlKChIcXZ2VkJDQ5XffvtN60h2Z+vWrXf8+Rk0aJDW0ezO3X4vffPNN1pHsysvvfSSUrlyZcXZ2VkpW7as0qlTJ2XTpk2a5dH9GhEhhBBC2K+SNbEvhBBCCLsihYgQQgghNCOFiBBCCCE0I4WIEEIIITQjhYgQQgghNCOFiBBCCCE0I4WIEEIIITQjhYgQQgghNCOFiBBCCCE0I4WIEEIIITQjhYgQQgghNCOFiBBCCCE08//z85YKHC8GAwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 画出判别面\n",
    "w = np.sum(alpha.reshape(alpha.shape[1], 1) * support_vectors, 0)\n",
    "b = svc.intercept_\n",
    "x1 = np.arange(-3, 3.1, 0.2)\n",
    "plt.plot(x1, (-w[0] * x1 - b) / w[1])\n",
    "plt.plot(x1, (-w[0] * x1 - b + 1) / w[1])\n",
    "plt.plot(x1, (-w[0] * x1 - b - 1) / w[1])\n",
    "\n",
    "# 画出样本点\n",
    "for i in range(40):\n",
    "    if real_class[i] == -1:\n",
    "        plt.scatter(X[i][0], X[i][1], c='b')\n",
    "    else:\n",
    "        plt.scatter(X[i][0], X[i][1], c='r')\n",
    "\n",
    "# 画出支持向量\n",
    "for i, x in enumerate(support_vectors):\n",
    "    if abs(alpha[0][i]) == 1:\n",
    "        plt.scatter(x[0], x[1], marker='x', s=50, c='black')\n",
    "    else:\n",
    "        plt.scatter(x[0], x[1], marker='+', s=70, c='black')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "图中的三条直线分别是判别函数$g(\\mathbf{x}) = 0$，$g(\\mathbf{x}) = 1$和$g(\\mathbf{x}) = -1$\n",
    "\n",
    "标记为+号的样本点是边界支持向量，标记为x号的样本点是错分支持向量。\n",
    "\n",
    "可以看到，边界支持向量在$g(\\mathbf{x}) = 1$和$g(\\mathbf{x}) = -1$上并且分类正确；但是错分支持向量不是只有分类错误，即在判别面对侧的样本，也包含了和判别面距离小于1的样本。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. 分类\n",
    "\n",
    "广义最优分类面的判别函数是：\n",
    "\n",
    "$$f(\\mathbf{x}) = \\text{sgn}\\{(\\mathbf{w}^*\\cdot\\mathbf{x}) + b\\}=\\text{sgn}\\{\\sum_{i=1}^{N}\\alpha_i^*y_i(\\mathbf{x}_i\\cdot\\mathbf{x}) + b^*\\}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 五、多类线性分类器\n",
    "\n",
    "#### 1. 多个两类分类器的组合\n",
    "\n",
    "书中介绍了三种实现方式：$c-1$个一对多、$\\frac{c(c-1)}{2}$个逐对分类、$c$个一对多分类器。\n",
    "\n",
    "注意c个一对多分类器中**不一定**可以直接比较线性判别函数的输出值，用来确定类别。因为$\\mathbf{w}$的尺度可能不相同，不同分类器的输出值不一定有可比性\n",
    "\n",
    "#### 2. 多类线性判别函数\n",
    "\n",
    "实际上是感知器的变形，也和感知器一样通过迭代来运算。\n",
    "\n",
    "书中对这个内容介绍篇幅非常少，应该不重要，而且和感知器非常类似，这里就不用代码实现了。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}