第09章 非监督模式识别.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 第九章 非监督模式识别\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 二轮总结笔记\n",
    "\n",
    "### 一、基于模型的方法（单峰子集分离）\n",
    "\n",
    "#### 1. 思想\n",
    "\n",
    "将样本投影到一维，在一维上估计概率密度，寻找单峰进行划分。\n",
    "\n",
    "#### 2. 算法\n",
    "\n",
    "1. 对**所有样本**进行主成分分析，选择第一主成分$\\mathbf{u}_j$方向对样本进行投影。\n",
    "2. 用**非参数方法**估计投影后的样本的概率密度函数，如直方图法。\n",
    "3. 在波谷（概率密度极小值点）做垂直于$\\mathbf{u}_j$的分类超平面，将样本划分为若干子集。如果没有波谷，则选择下一主成分进行投影。\n",
    "4. 对每个子集重复上述步骤，直到得到预想的的聚类数量或者每个子集都是单峰分布。注意每次进行投影使用的主成分都是**所有样本**的主成分，包括子集的投影。\n",
    "\n",
    "#### 3. 性质\n",
    "\n",
    "1. 把方差最大的方向作为投影方向。\n",
    "2. 不一定适用于所有情况。\n",
    "\n",
    "\n",
    "### 二、混合模型的估计（非监督参数估计问题）\n",
    "\n",
    "#### 1. 非监督最大似然估计\n",
    "\n",
    "##### (1) 假设条件\n",
    "\n",
    "假设除了几个参数外，已经完全知道了问题的概率结构，即\n",
    "\n",
    "1. 类别数量已知，但是不知道样本来自哪个类。\n",
    "2. 每类先验概率$P(\\omega_j)$已知。\n",
    "3. 每类概率密度函数形式$p(\\mathbf{x}|\\omega_j,\\boldsymbol\\theta_j)$已知。\n",
    "4. 未知的只有$c$个参数向量$\\boldsymbol\\theta_j$。\n",
    "\n",
    "##### (2) 似然函数\n",
    "\n",
    "$$\\begin{aligned}\n",
    "p(\\mathbf{x}|\\boldsymbol\\theta) &= \\sum_{i=1}^cp(\\mathbf{x}|\\omega_j,\\boldsymbol\\theta_j)P(\\omega_j)\\\\\n",
    "l(\\boldsymbol\\theta) &= p(\\mathscr{X}|\\boldsymbol\\theta)= \\prod_{k=1}^Np(\\mathbf{x}_k|\\boldsymbol\\theta) \\\\\n",
    "H(\\boldsymbol\\theta) &= \\ln[l(\\boldsymbol\\theta)] = \\sum_{k=1}^N \\ln p(\\mathbf{x}_k|\\boldsymbol\\theta)\n",
    "\\end{aligned}$$\n",
    "\n",
    "其中，$P(\\omega_j)$称为**混合参数**，$p(\\mathbf{x}|\\omega_j,\\boldsymbol\\theta_j)$称为**分量密度**。\n",
    "\n",
    "##### (3) 可识别性问题\n",
    "\n",
    "1. 常见**连续**随机变量的分布密度函数都是**可识别**的。\n",
    "2. **离散**随机变量的混合概率函数往往使**不可识别**的。\n",
    "\n",
    "##### (4) 计算问题\n",
    "\n",
    "1. 先验概率已知：\n",
    "\n",
    "最大似然估计解$\\hat{\\boldsymbol\\theta}$满足\n",
    "\n",
    "$$\n",
    "\\nabla_{\\boldsymbol\\theta_i}H(\\hat{\\boldsymbol\\theta}) = \\sum_{k=1}^N P(\\omega_i|\\mathbf{x}_k,\\hat{\\boldsymbol\\theta}_i)\\nabla_{\\boldsymbol\\theta_i}\\ln p(\\mathbf{x}_k|\\omega_i,\\hat{\\boldsymbol\\theta}_i) = 0,\\quad i=1,\\cdots,c\n",
    "$$\n",
    "\n",
    "2. 先验概率未知：\n",
    "\n",
    "最大似然估计解$\\hat{\\boldsymbol\\theta}$和$\\hat{P}(\\omega_i)$满足\n",
    "\n",
    "$$\\begin{aligned}\n",
    "&\\sum_{k=1}^N \\hat{P}(\\omega_i|\\mathbf{x}_k,\\hat{\\boldsymbol\\theta}_i)\\nabla_{\\boldsymbol\\theta_i}\\ln p(\\mathbf{x}_k|\\omega_i,\\hat{\\boldsymbol\\theta}_i) = 0,\\\\\n",
    "&\\hat{P}(\\omega_i|\\mathbf{x}_k,\\hat{\\boldsymbol\\theta}_i) = \\frac{p(\\mathbf{x}_k|\\omega_i,\\hat{\\boldsymbol\\theta}_i)\\hat{P}(\\omega_i)}{\\displaystyle \\sum_{j=1}^cp(\\mathbf{x}_k|\\omega_i,\\hat{\\boldsymbol\\theta}_i)\\hat{P}(\\omega_i)}, \n",
    "\\end{aligned}\\quad i=1,\\cdots,c$$\n",
    "\n",
    "3. 性质\n",
    "\n",
    "求出封闭解很困难，通常使用迭代法求解。\n",
    "\n",
    "#### 2. 正态分布情况下的非监督参数估计\n",
    "\n",
    "##### (1) 只有均值向量$\\boldsymbol\\mu_i$未知\n",
    "\n",
    "1. 解的方程：\n",
    "\n",
    "$$\\begin{aligned}\n",
    "\\hat{\\boldsymbol\\mu}_i &= \\frac{\\displaystyle\\sum_{k=1}^NP(\\omega_i|\\mathbf{x}_k,\\hat{\\boldsymbol\\mu}_i)\\mathbf{x}_k}{\\displaystyle\\sum_{k=1}^NP(\\omega_i|\\mathbf{x}_k,\\hat{\\boldsymbol\\mu}_i)} \\\\\n",
    "P(\\omega_i|\\mathbf{x}_k,\\hat{\\boldsymbol\\mu}_i) &= \\frac{p(\\mathbf{x}_k|\\omega_i,\\hat{\\boldsymbol\\mu}_i)P(\\omega_i)}{\\displaystyle\\sum_{j=1}^cp(\\mathbf{x}_k|\\omega_j,\\hat{\\boldsymbol\\mu}_j)P(\\omega_j)}\n",
    "\\end{aligned}$$\n",
    "\n",
    "2. 算法：\n",
    "\n",
    "+ 迭代公式：\n",
    "\n",
    "$$\n",
    "\\hat{\\boldsymbol\\mu}_i(j+1) = \\frac{\\displaystyle\\sum_{k=1}^NP(\\omega_i|\\mathbf{x}_k,\\hat{\\boldsymbol\\mu}_i(j))\\mathbf{x}_k}{\\displaystyle\\sum_{k=1}^NP(\\omega_i|\\mathbf{x}_k,\\hat{\\boldsymbol\\mu}_i(j))}\n",
    "$$\n",
    "\n",
    "+ 缺点：一般梯度法的缺点：得到的可能是局部最优解，受$\\hat{\\boldsymbol\\mu}_i(0)$初值影响。\n",
    "\n",
    "3. 性质：$\\boldsymbol\\mu_i$的最大似然估计是样本加权平均。\n",
    "\n",
    "##### (2) 均值向量$\\boldsymbol\\mu_i$、协方差矩阵$\\boldsymbol\\Sigma_i$、先验概率$P(\\omega_i)$未知，但是类别数$c$已知\n",
    "\n",
    "1. 解的方程：\n",
    "\n",
    "$$\\begin{aligned}\n",
    "\\hat{P}(\\omega_i) &= \\frac{1}{N}\\sum_{k=1}^N\\hat{P}(\\omega_i|\\mathbf{x}_k,\\hat{\\boldsymbol\\theta}_i)\\\\\n",
    "\\hat{\\boldsymbol\\mu}_i &= \\frac{\\displaystyle\\sum_{k=1}^NP(\\omega_i|\\mathbf{x}_k,\\hat{\\boldsymbol\\theta}_i)\\mathbf{x}_k}{\\displaystyle\\sum_{k=1}^NP(\\omega_i|\\mathbf{x}_k,\\hat{\\boldsymbol\\theta}_i)} \\\\\n",
    "\\hat{\\Sigma}_i &= \\frac{\\displaystyle\\sum_{k=1}^NP(\\omega_i|\\mathbf{x}_k,\\hat{\\boldsymbol\\theta}_i)(\\mathbf{x}_k-\\hat{\\boldsymbol\\mu}_i)(\\mathbf{x}_k-\\hat{\\boldsymbol\\mu}_i)^\\text{T}}{\\displaystyle\\sum_{k=1}^NP(\\omega_i|\\mathbf{x}_k,\\hat{\\boldsymbol\\theta}_i)}\\\\\n",
    "P(\\omega_i|\\mathbf{x}_k,\\hat{\\boldsymbol\\theta}_i) &= \\frac{p(\\mathbf{x}_k|\\omega_i,\\hat{\\boldsymbol\\theta}_i)P(\\omega_i)}{\\displaystyle\\sum_{j=1}^cp(\\mathbf{x}_k|\\omega_j,\\hat{\\boldsymbol\\theta}_j)P(\\omega_j)}\n",
    "\\end{aligned}$$\n",
    "\n",
    "2. 算法：迭代法。\n",
    "3. 性质： 所有样本都对估计值起某种作用，这些估计仍然是加权的频数比、样本均值和样本协方差矩阵。\n",
    "\n",
    "##### (3) 均值向量$\\boldsymbol\\mu_i$、协方差矩阵$\\boldsymbol\\Sigma_i$、先验概率$P(\\omega_i)$、类别数$c$全都未知\n",
    "\n",
    "无法用最大似然估计求解。\n",
    "\n",
    "\n",
    "### 三、动态聚类算法\n",
    "\n",
    "#### 1. C均值（C-means）/k均值（k-means）/广义Llogd（GLA）算法\n",
    "\n",
    "##### (1) 原理\n",
    "\n",
    "最小平方划分：使误差平方和准则$J_\\text{e}$最小\n",
    "\n",
    "$$\n",
    "J_\\text{e} = \\sum_{i=c}^c\\sum_{\\mathbf{y}\\in\\Gamma_i}||\\mathbf{y}-\\mathbf{m}_i||^2\n",
    "$$\n",
    "\n",
    "##### (2) 算法\n",
    "\n",
    "1. 初始划分$c$个聚类，计算均值$\\mathbf{m}_i$。\n",
    "2. 选取一个样本$\\mathbf{y}$，设$\\mathbf{y}\\in\\Gamma_i$。\n",
    "3. 计算误差平方和的增加量和减少量\n",
    "\n",
    "$$\\begin{aligned}\n",
    "\\rho_j &= \\frac{N_j}{N_j+1}||\\mathbf{y}-\\mathbf{m}_j||^2,\\quad j\\neq i\\\\\n",
    "\\rho_i &= \\frac{N_i}{N_i-1}||\\mathbf{y}-\\mathbf{m}_i||^2\n",
    "\\end{aligned}$$\n",
    "\n",
    "4. 考察增加量$\\rho_j$的最小值$\\rho_k$，若$\\rho_k<\\rho_i$，则把$\\mathbf{y}$从$\\Gamma_i$移动到$\\Gamma_k$中。\n",
    "5. 重新计算均值$\\mathbf{m}_i$，另选样本继续迭代，直到全部样本在一轮迭代中均不移动，则收敛。\n",
    "\n",
    "##### (3) 样本初始划分方法\n",
    "\n",
    "1. 代表样本点的选择：\n",
    "\n",
    "+ 凭经验选。\n",
    "+ 将数据随机分成$c$类，计算每类重心，重心作为代表点。\n",
    "+ 用密度法选择代表点，找密度最大的点作为代表点。\n",
    "+ 按照样本天然顺序或随机排序后选择前$c$个点作为代表点。\n",
    "+ 从$(c-1)$聚类划分问题的解中产生$c$聚类划分问题的代表点。\n",
    "\n",
    "2. 初始划分分类：\n",
    "\n",
    "+ 归入最近的代表点。\n",
    "+ 归入最近的代表点同时更新代表点为该类的重心（均值）。\n",
    "+ 距离阈值划分法，小于阈值划入，大于阈值划分新类。\n",
    "+ 标准化后归类，$\\displaystyle\\frac{(c-1)[\\text{SUM}(i)-\\text{MI}]}{\\text{MA}-\\text{MI}} + 1$最接近的整数就是类别序号。\n",
    "\n",
    "##### (4) 性质\n",
    "\n",
    "1. 使局部搜索算法，不能保证收敛到全局最优解。算法的结果收到初始化分和样本调整顺序的影响。\n",
    "2. 误差平方和随聚类数目$c$增加单调下降，在拐点处的聚类数量$c$为最优。\n",
    "\n",
    "#### 2. ISODATA方法（Iterative Self-Organization Data Analysis Technique）\n",
    "\n",
    "##### (1) 特点\n",
    "\n",
    "1. 分裂：一类中样本**特征**最大的**方差**大于阈值，则进行分裂。\n",
    "2. 合并：不同类**均值**之间的**距离**小于阈值，则进行合并。\n",
    "3. 最终类别数可以不等于期望类别数。\n",
    "\n",
    "#### 3. 基于样本与核的相似性度量的动态聚类算法\n",
    "\n",
    "##### (1) 准则函数\n",
    "\n",
    "$$\n",
    "J_K = \\sum_{j=1}^c\\sum_{\\mathbf{y}\\in\\Gamma_j}\\Delta(\\mathbf{y},K_j)\n",
    "$$\n",
    "\n",
    "##### (2) 算法\n",
    "\n",
    "1. 初始化分为$c$类，初始化每类的初始核$K_j$。\n",
    "2. 把每个样本移动到准则函数$J_K$最小的类中。\n",
    "3. 修正核$K_j$，迭代直至收敛。\n",
    "\n",
    "##### (3) 核与相似性度量\n",
    "\n",
    "1. 正态核函数：\n",
    "\n",
    "$$\\begin{aligned}\n",
    "K_j(\\mathbf{y},V_j) &= \\frac{1}{(2\\pi)^{\\frac{d}{2}}|\\hat{\\boldsymbol\\Sigma}_j|^{\\frac{1}{2}}}\\exp\\left\\{-\\frac{1}{2}(\\mathbf{y}-\\mathbf{m}_j)^\\text{T}\\hat{\\boldsymbol\\Sigma}{_j^{-1}}(\\mathbf{y}-\\mathbf{m}_j)\\right\\} \\\\\n",
    "\\Delta(\\mathbf{y},K_j) &= \\frac{1}{2}(\\mathbf{y}-\\mathbf{m}_j)^\\text{T}\\hat{\\boldsymbol\\Sigma}{_j^{-1}}(\\mathbf{y}-\\mathbf{m}_j)+\\frac{1}{2}\\log |\\hat{\\boldsymbol\\Sigma}_j|\n",
    "\\end{aligned}$$\n",
    "\n",
    "2. 主轴核函数：\n",
    "\n",
    "$$\\begin{aligned}\n",
    "K_j(\\mathbf{y},V_j) &= \\boldsymbol{U}_j^\\text{T}\\mathbf{y}\\\\\n",
    "\\Delta(\\mathbf{y},K_j) &= [(\\mathbf{y}-\\mathbf{m}_j)-\\mathbf{U}_j\\mathbf{U}_j^T(\\mathbf{y}-\\mathbf{m}_j)]^\\text{T}[(\\mathbf{y}-\\mathbf{m}_j)-\\mathbf{U}_j\\mathbf{U}_j^T(\\mathbf{y}-\\mathbf{m}_j)]\n",
    "\\end{aligned}$$\n",
    "$\\Delta(\\mathbf{y},K_j)$为样本与主轴之间欧氏距离的平方。\n",
    "\n",
    "##### (4) 性质\n",
    "\n",
    "C均值是基于样本与核的相似性度量的动态聚类算法的一个特例。\n",
    "\n",
    "\n",
    "### 四、模糊聚类方法\n",
    "\n",
    "#### 1. 模糊集\n",
    "\n",
    "##### (1) 隶属度函数\n",
    "\n",
    "$\\mu_A(\\mathbf{x})$表示一个对象$\\mathbf{x}$隶属于集合$A$的程度，$\\mu_A(\\mathbf{x})=1$表示$\\mathbf{x}$完全属于$A$，$\\mu_A(\\mathbf{x})=0$表示$\\mathbf{x}$完全不属于$A$。\n",
    "\n",
    "##### (2) 支持集和支持点\n",
    "\n",
    "$\\mu_A(\\mathbf{x})>0$的对象的集合叫做模糊集$A$的支持集$S(A)$。支持集中的元素称作模糊集$A$的支持点，或模糊集$A$的元素。\n",
    "\n",
    "#### 2. 模糊C均值算法（Fuzzy C-Means，FCM）\n",
    "\n",
    "##### (1) 损失函数\n",
    "\n",
    "$$\\begin{aligned}\n",
    "&J_\\text{f} = \\sum_{j=1}^c\\sum_{i=1}^n[\\mu_j(\\mathbf{x}_i)]^b||\\mathbf{x}_i-\\mathbf{m}_j||^2 \\\\\n",
    "&\\sum_{j=1}^c\\mu_j(\\mathbf{x}_i)=1\n",
    "\\end{aligned}$$\n",
    "\n",
    "##### (2) 算法\n",
    "\n",
    "迭代以下公式直至收敛\n",
    "\n",
    "$$\\begin{aligned}\n",
    "\\mathbf{m}_i &= \\frac{\\displaystyle\\sum_{i=1}^n[\\mu_j(\\mathbf{x}_i)]^b\\mathbf{x}_i}{\\displaystyle\\sum_{i=1}^n[\\mu_j(\\mathbf{x}_i)]^b}\\\\\n",
    "\\mu_j(\\mathbf{x}_i) &= \\frac{1}{||\\mathbf{x}_i-\\mathbf{m}_j||^{\\frac{1}{b-1}}\\displaystyle \\sum_{k=1}^c\\frac{1}{||\\mathbf{x}_i-\\mathbf{m}_k||^{\\frac{1}{b-1}}}}\n",
    "\\end{aligned}\\quad i=1,\\cdots,n,\\quad j=1,\\cdots,c$$\n",
    "\n",
    "##### (3) $b$的选择\n",
    "\n",
    "1. $b>1$。\n",
    "2. 如果$b\\rightarrow 1$，则FCM等同于C均值。\n",
    "3. 如果$b = \\infty$，则算法得到完全模糊的解。\n",
    "\n",
    "#### 3. 改进的模糊C均值算法\n",
    "\n",
    "##### (1) 损失函数\n",
    "\n",
    "$$\\begin{aligned}\n",
    "&J_\\text{f} = \\sum_{j=1}^c\\sum_{i=1}^n[\\mu_j(\\mathbf{x}_i)]^b||\\mathbf{x}_i-\\mathbf{m}_j||^2 \\\\\n",
    "&\\sum_{j=1}^c\\sum_{i=1}^n\\mu_j(\\mathbf{x}_i)=n\n",
    "\\end{aligned}$$\n",
    "\n",
    "##### (2) 算法\n",
    "\n",
    "迭代以下公式直至收敛\n",
    "\n",
    "$$\\begin{aligned}\n",
    "\\mathbf{m}_i &= \\frac{\\displaystyle\\sum_{i=1}^n[\\mu_j(\\mathbf{x}_i)]^b\\mathbf{x}_i}{\\displaystyle\\sum_{i=1}^n[\\mu_j(\\mathbf{x}_i)]^b}\\\\\n",
    "\\mu_j(\\mathbf{x}_i) &= \\frac{n}{||\\mathbf{x}_i-\\mathbf{m}_j||^{\\frac{1}{b-1}}\\displaystyle \\sum_{k=1}^c\\sum_{l=1}^n\\frac{1}{||\\mathbf{x}_l-\\mathbf{m}_k||^{\\frac{1}{b-1}}}}\n",
    "\\end{aligned}\\quad i=1,\\cdots,n,\\quad j=1,\\cdots,c$$\n",
    "\n",
    "##### (3) 性质\n",
    "\n",
    "1. 得到的隶属度函数值可能大于$1$。\n",
    "2. 改进后的FCM比普通FCM有更好的鲁棒性，有**野值**的情况下效果好。\n",
    "3. 对预先确定的**聚类数目不敏感**。\n",
    "4. 对聚类中心**初值敏感**。\n",
    "5. 如果迭代过程中出现某个聚类中心距离某个样本非常近，可能会得到只包含这**一个样本**的聚类。\n",
    "\n",
    "#### 4. C均值、普通FCM和改进的FCM对比\n",
    "\n",
    "##### (1) 预先确定的类别数比实际多（不同类有重叠）时\n",
    "\n",
    "C均值和普通FCM不会反应重叠，改进的FCM会反应重叠。改进的FCM能合理的给出比预先确定的数目少的聚类，前两个不能。\n",
    "\n",
    "##### (2) 预先确定的类别数比实际少时\n",
    "\n",
    "C均值和普通FCM会把不同的类混在一起，改进的FCM会正确给出聚类中心。\n",
    "\n",
    "\n",
    "### 五、分级聚类方法（hierarchical clustering）\n",
    "\n",
    "#### 1. 思想\n",
    "\n",
    "自底向上。\n",
    "\n",
    "#### 2. 算法\n",
    "\n",
    "##### (1) 自底向上算法\n",
    "\n",
    "1. 初始化：每个样本形成一个类。\n",
    "2. 合并：计算两两类别之间的距离或相似性，合并最小的。迭代直至合并为两类。\n",
    "\n",
    "##### (2) 距离或相似性度量\n",
    "\n",
    "1. 样本之间的距离或相似性度量：**欧氏距离**或**相关系数**。\n",
    "2. 两类之间的距离或相似性度量（类间的连接）：\n",
    "+ 最近距离（single linkage或single-link）：$\\displaystyle \\Delta(\\Gamma_i,\\Gamma_j)=\\min\\limits_{\\mathbf{y}\\in\\Gamma_i,\\tilde{\\mathbf{y}}\\in\\Gamma_j}\\delta(\\mathbf{y},\\tilde{\\mathbf{y}})$\n",
    "+ 最远距离（complete linkage或complete-link）：$\\displaystyle \\Delta(\\Gamma_i,\\Gamma_j)=\\max\\limits_{\\mathbf{y}\\in\\Gamma_i,\\tilde{\\mathbf{y}}\\in\\Gamma_j}\\delta(\\mathbf{y},\\tilde{\\mathbf{y}})$\n",
    "+ 均值距离（average linkage或average-link）：$\\displaystyle \\Delta(\\Gamma_i,\\Gamma_j)=\\delta(\\mathbf{m}_i,\\mathbf{m}_j)$\n",
    "\n",
    "#### 3. 性质\n",
    "\n",
    "1. 分级聚类是局部搜素算法，算法对噪声比较敏感，个别样本变动可能会对聚类结果影响很大，样本数量较多时，影响比较微弱。\n",
    "2. 聚类树树枝长度表示两类间相似度或距离。\n",
    "3. 聚类树画法不唯一，可以对调位置。\n",
    "\n",
    "#### 4. 双向聚类与双聚类\n",
    "\n",
    "##### (1) 双向聚类（two-way clustering）\n",
    "\n",
    "分别对数据的两种排列索引做聚类。\n",
    "\n",
    "##### (2) 双向聚类（biclustering）\n",
    "\n",
    "寻找一组特征和一组样本，使这些样本在找出的特征上呈现为很强的聚类。\n",
    "\n",
    "\n",
    "### 六、自组织映射神经网络（self-organizing map，SOM）\n",
    "\n",
    "#### 1. 思想\n",
    "\n",
    "1. 属于竞争学习型神经网络。\n",
    "2. 神经元节点在同一层上，在一个平面上呈规则排列，常见的有**方形网格**排列和**蜂窝状**排列。\n",
    "3. 神经元节点之间没有直接联系，但是相邻的节点构成邻域相互作用，作用范围随训练次数增加而缩小。\n",
    "4. 相应强度是神经元节点的权值向量和输入向量的匹配程度，可以用欧氏距离或内积表示，欧氏距离最小或者内积最大的是**获胜节点**（winnner）。\n",
    "\n",
    "#### 2. 算法\n",
    "\n",
    "##### (1) 算法内容\n",
    "\n",
    "1. 用**小随机数**初始化神经元的权值向量，初始值不能相等。\n",
    "2. 按随即顺序或指定顺序遍历样本训练。\n",
    "3. 计算样本的神经元响应，找到获胜节点c。\n",
    "4. 权值竞争学习。对所有神经元节点，更新各自权值$\\displaystyle\\mathbf{m}_i = \\begin{cases}\\mathbf{m}_i+\\alpha(t)||\\mathbf{x}(t)-\\mathbf{m}_i(t)||&i\\in N_c(t)\\\\\\mathbf{m}_i(t)&i\\notin N_c(t)\\end{cases}$，$\\alpha(t)$是学习步长。\n",
    "5. 更新步长$\\alpha(t)$和邻域$N_c(t)$，继续迭代直到收敛。\n",
    "\n",
    "##### (2) 注意事项\n",
    "\n",
    "1. SOM是自学习过程，没有训练误差之类的收敛目标，终止条件是确定的迭代次数。\n",
    "2. 步长$\\alpha(t)$和邻域$N_c(t)$随时间增加单调减小，到算法终止时邻域缩小到只包含获胜节点。\n",
    "3. 理论上无法严格证明，确定参数需要试算。\n",
    "4. 邻域函数还可以用高斯函数。\n",
    "\n",
    "#### 3. 性质\n",
    "\n",
    "##### (1) 像、原像和像密度\n",
    "\n",
    "1. 每个输入样本$\\mathbf{x}$的获胜节点会逐渐趋于固定，固定的获胜节点叫样本$\\mathbf{x}$的**像**，样本$\\mathbf{x}$称为获胜节点的**原像**。一个样本只能由一个像，一个神经元可能有多个原像。\n",
    "2. 每个神经元上有多少个原像叫做**像密度**，把各个节点的像密度按照神经元本来的排列图示出来，就得到**像密度图**。\n",
    "\n",
    "##### (2) 自组织现象\n",
    "\n",
    "1. 原样本空间内距离相近的样本趋向于映射到**同一个**像节点或在神经元平面**排列相近**的像节点。节点的像密度与原空间的样本密度形成近似的**单调关系**。\n",
    "2. SOM完成的是从原样本空间到二维平面上神经网络的映射，这种映射时**拓扑保持**的。\n",
    "3. SOM是一种映射空间高度网格化的**非线性特征变换**。\n",
    "\n",
    "##### (3) 与C均值的关系\n",
    "\n",
    "1. 如果邻域函数固定（不随学习时间改变），则SOM是通过梯度下降最小化势函数$\\displaystyle \\sum_k\\sum_ih_{ci}||\\mathbf{x}_k-\\mathbf{m}_i||^2$。\n",
    "2. 如果取消邻域作用，**只修正获胜节点的权值**，则SOM退化为**C均值**算法的一种随机迭代实现，**聚类数**就是**神经元节点数目**。\n",
    "\n",
    "#### 4. SOM用于模式识别\n",
    "\n",
    "##### (1) 用于监督学习\n",
    "\n",
    "1. 思想：\n",
    "\n",
    "+ 先用样本训练SOM网络，训练后每个神经元的原像会趋于一类，相邻的神经元的原像也来自相同或相近的类别。\n",
    "+ 识别新样本时，考察获胜节点的原像类别即可。\n",
    "\n",
    "2. 特点：\n",
    "\n",
    "+ 在非监督学习的基础上进行监督模式识别。\n",
    "+ 相邻聚类之间有相似性，即使识别错误也倾向于识别为相同或类似的类别，**接近人的识别特性**。\n",
    "\n",
    "##### (2) 用于非监督学习（自组织映射/自组织分析/SOMA）\n",
    "\n",
    "1. 思想：\n",
    "\n",
    "画出像密度图，在像密度较高且集中的地方分为一类。\n",
    "\n",
    "2. 特点：\n",
    "\n",
    "+ **无须事先确定聚类数目**，能更好的适应不同的分布情况，是有效的聚类方法。\n",
    "+ 能反映出样本集中**无明显聚类**的分布特性。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## 一轮学习笔记（包含代码实现）\n",
    "\n",
    "这一章介绍的是非监督学习中的聚类。主要分为两个思路：一个是根据概率模型进行聚类，另一个是根据数据设计算法进行聚类。\n",
    "\n",
    "书中在介绍概率模型聚类用到了大篇幅的概率公式推导，可以说完全看不懂，我尝试了一个小时之后彻底放弃，只能赌一手考研不考。至于到底考不考，等我今年考完研就知道了。\n",
    "\n",
    "下面是书上介绍的基于数据的算法。\n",
    "\n",
    "### 一、C均值算法\n",
    "\n",
    "C均值算法的思想是最小化误差平方和，即最小化$\\displaystyle J_e = \\sum_{i=1}^c\\sum_{\\mathbf{y}\\in\\Gamma_i}||\\mathbf{y}-\\mathbf{m}_i||^2$。\n",
    "\n",
    "具体算法如下：\n",
    "\n",
    "1. 初始化划分$c$个聚类，具体初始化划分书上给了好多方法，不过都是启发式的算法，具体怎么初始化划分都可以。\n",
    "2. 遍历所有样本，如果当前样本所在的类的样本数量$N_i > 1$，则比较把这个样本移动到别的类造成的误差平方和变化，即比较：\n",
    "\n",
    "$$\\begin{align}\n",
    "\\rho_j &= \\frac{N_j}{N_j + 1}||\\mathbf{y}-\\mathbf{m}_j||^2, \\quad j\\neq i \\\\\n",
    "\\rho_i &= \\frac{N_i}{N_i - 1}||\\mathbf{y}-\\mathbf{m}_i||^2\n",
    "\\end{align}$$\n",
    "\n",
    "其中$\\mathbf{m}_i$是第$i$类的均值。\n",
    "\n",
    "3. 考察$\\rho_j$中的最小值$\\rho_k$，如果$\\rho_k<\\rho_i$，则把当前样本从$\\Gamma_i$类移动到$\\Gamma_k$类中。移动后重新计算$\\mathbf{m}_i$是第$i$。\n",
    "4. 如果当前一轮遍历把所有样本都考察过了，并且没有移动任何样本，说明已经达到了局部的最小误差平方和，认为收敛。\n",
    "\n",
    "下面实现的代码中，我采用了寻找代表点的方法(2)和初始化分类的方法(1)。其实这些方法实际效果区别不大，用简单的就可以了。\n",
    "\n",
    "实现如下："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.datasets import load_iris\n",
    "from sklearn.preprocessing import StandardScaler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "class CMeans:\n",
    "    \"\"\"\n",
    "    实现C均值算法的类\n",
    "    \"\"\"\n",
    "    def __init__(self, c: int):\n",
    "        \"\"\"\n",
    "        构造函数\n",
    "        :param c: 指定的聚类数量\n",
    "        \"\"\"\n",
    "        self._c: int = c\n",
    "        # self._c表示聚类数\n",
    "\n",
    "        self._d: int = 0\n",
    "        # self._d表示样本维数\n",
    "\n",
    "        self._N: int = 0\n",
    "        # self._N表示样本总数\n",
    "\n",
    "        self._X: np.ndarray = np.zeros(0)\n",
    "        # self._X存储所有样本\n",
    "\n",
    "        self._clusters: list[set[int]] = []\n",
    "        # self._clusters[i] = {1, 3} 表示第i个聚类中包含1号样本和2号样本\n",
    "\n",
    "        self.m: list[np.ndarray] = [np.zeros(0)] * c\n",
    "        # self._m[i]表示第i个聚类的均值\n",
    "\n",
    "        self.costs: list[float] = []\n",
    "        # self.costs记录每次迭代后的误差平方和\n",
    "\n",
    "    def _update_mean(self):\n",
    "        \"\"\"\n",
    "        更新每一类的均值\n",
    "        \"\"\"\n",
    "        for i in range(self._c):\n",
    "            cluster_sum = np.zeros(self._d)\n",
    "            for j in self._clusters[i]:\n",
    "                cluster_sum += self._X[j]\n",
    "            self.m[i] = cluster_sum / len(self._clusters[i])\n",
    "\n",
    "    def initialize(self, _X: np.ndarray):\n",
    "        \"\"\"\n",
    "        初始化聚类\n",
    "        :param _X: 样本矩阵，每一行是一个样本（N * d）\n",
    "        \"\"\"\n",
    "        if _X.shape[0] <= self._c:\n",
    "            raise Exception(\"样本数量必须比分类数量大！\")\n",
    "\n",
    "        self._X = _X.copy()\n",
    "        self._N = _X.shape[0]\n",
    "        self._d = _X.shape[1]\n",
    "\n",
    "        np.random.shuffle(_X)\n",
    "        for i in range(self._c):\n",
    "            self._clusters.append({i})\n",
    "        self._update_mean()\n",
    "\n",
    "        for i in range(self._c, _X.shape[0]):\n",
    "            min_distance = 999999\n",
    "            min_index = 0\n",
    "            for j in range(self._c):\n",
    "                distance = np.sum((_X[i] - self.m[j]) ** 2)\n",
    "                if min_distance > distance:\n",
    "                    min_distance = distance\n",
    "                    min_index = j\n",
    "            self._clusters[min_index].add(i)\n",
    "            self._update_mean()\n",
    "\n",
    "    def _compute_cost(self) -> float:\n",
    "        \"\"\"\n",
    "        计算误差平方和\n",
    "        :return: 误差平方和\n",
    "        \"\"\"\n",
    "        ret = 0\n",
    "        for i in range(self._c):\n",
    "            for j in self._clusters[i]:\n",
    "                ret += np.sum((self._X[j] - self.m[i]) ** 2)\n",
    "        return ret\n",
    "\n",
    "    def cluster(self):\n",
    "        \"\"\"\n",
    "        开始运行迭代聚类算法\n",
    "        \"\"\"\n",
    "        clusters_order = np.random.permutation(self._c)\n",
    "        # 随机生成遍历的顺序\n",
    "\n",
    "        iters = 0\n",
    "        confirmed_samples = 0\n",
    "        while confirmed_samples < self._N:\n",
    "            flag = False\n",
    "            for c in clusters_order:\n",
    "                if flag:\n",
    "                    break\n",
    "                cluster = list(self._clusters[c])\n",
    "                N_i = len(self._clusters[c])\n",
    "                samples_order = np.random.permutation(N_i)\n",
    "                # 随机生成遍历的顺序\n",
    "                for i in samples_order:\n",
    "                    if flag:\n",
    "                        break\n",
    "                    rho_i = np.sum(N_i / (N_i - 1) * (self._X[cluster[i]] - self.m[c]) ** 2)\n",
    "                    min_value_of_rho_j = 999999\n",
    "                    min_cluster_of_rho_j = 0\n",
    "                    for j in range(self._c):\n",
    "                        if c == j:\n",
    "                            continue\n",
    "                        N_j = len(self._clusters[j])\n",
    "                        rho_j = np.sum(N_j / (N_j + 1) * (self._X[cluster[i]] - self.m[j]) ** 2)\n",
    "                        if min_value_of_rho_j > rho_j:\n",
    "                            min_value_of_rho_j = rho_j\n",
    "                            min_cluster_of_rho_j = j\n",
    "                    iters += 1\n",
    "                    if min_value_of_rho_j < rho_i:\n",
    "                        self._clusters[c].remove(cluster[i])\n",
    "                        self._clusters[min_cluster_of_rho_j].add(cluster[i])\n",
    "                        self._update_mean()\n",
    "                        self.costs.append(self._compute_cost())\n",
    "                        confirmed_samples = 0\n",
    "                        flag = True\n",
    "                        break\n",
    "                    confirmed_samples += 1\n",
    "                    if iters >= 1000000:\n",
    "                        raise Exception(\"迭代了100000次还不收敛，数据可能有问题\")\n",
    "\n",
    "    def result(self) -> np.ndarray:\n",
    "        \"\"\"\n",
    "        获取聚类结果\n",
    "        :return: 聚类结果向量\n",
    "        \"\"\"\n",
    "        ret = np.zeros(self._N, dtype=int)\n",
    "        for i, cluster in enumerate(self._clusters):\n",
    "            for j in cluster:\n",
    "                ret[j] = i\n",
    "        return ret"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下来用水仙花的数据来测试一下聚类效果。\n",
    "\n",
    "因为C均值算法考虑的是样本误差平方和，所以需要将各维度大小进行统一，防止某一个维度的影响过大。所以下面用`StandardScaler`对样本进行预处理。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "iris = load_iris()\n",
    "X = iris.data\n",
    "y = iris.target\n",
    "X = StandardScaler().fit_transform(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "原始数据类别为：\n",
      "[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0 0 0 0 0 0 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 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2\n",
      " 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2\n",
      " 2 2]\n",
      "聚类获得的类别为：\n",
      "[2 1 2 2 2 0 2 2 0 2 1 0 1 1 1 1 0 2 0 0 0 0 2 0 2 2 0 1 0 0 0 2 1 0 2 2 1\n",
      " 0 2 1 0 0 1 1 2 0 2 1 0 2 2 2 0 0 2 1 0 1 1 1 1 1 0 1 0 0 1 0 1 0 1 0 1 0\n",
      " 2 1 2 0 0 1 2 1 1 0 1 0 0 0 0 2 1 1 2 0 1 2 1 2 0 2 2 1 0 1 0 0 2 1 0 2 2\n",
      " 1 2 2 2 0 2 2 0 1 0 0 0 1 0 1 2 1 2 2 1 2 1 0 2 1 2 2 2 2 2 1 0 2 0 1 0 2\n",
      " 1 0]\n"
     ]
    }
   ],
   "source": [
    "cmeans = CMeans(3)\n",
    "cmeans.initialize(X)\n",
    "cmeans.cluster()\n",
    "y_clustered = cmeans.result()\n",
    "print(\"原始数据类别为：\")\n",
    "print(y)\n",
    "print(\"聚类获得的类别为：\")\n",
    "print(y_clustered)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可以看到聚类获得类别和原始类别非常不同，这是非常正常的现象。\n",
    "\n",
    "接下来可视化的观察一下随着迭代，误差平方和的变化。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkEAAAGwCAYAAACuIrGMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABIWklEQVR4nO3deXxU1f3/8ddkJQGSEAgJIQmoIItK1KBINQqSsmgRDPzqQhUtX2krKFu1pf2q1S5o64pSbast+vgqWCyguFAohEUJFCIRVGQzmABZwEhCEgghc39/3GTIkG0mmcmdZN7Px2MeM3PvmTufmVDn03M+5xybYRgGIiIiIn4mwOoARERERKygJEhERET8kpIgERER8UtKgkRERMQvKQkSERERv6QkSERERPySkiARERHxS0FWB+AL7HY7R48epWvXrthsNqvDERERERcYhsHJkyeJj48nIMD9fh0lQcDRo0dJTEy0OgwRERFpgby8PBISEtx+nZIgoGvXroD5JUZERFgcjYiIiLiitLSUxMREx++4u5QEgWMILCIiQkmQiIhIO9PSUhYVRouIiIhfUhIkIiIifklJkIiIiPglJUEiIiLil5QEiYiIiF9SEiQiIiJ+SUmQiIiI+CUlQSIiIuKXlASJiIiIX9KK0d5ir4Zjm+FUPoT1gphUCAi0OioRERGpoSTIG/KWQ9YsqDh87lh4AqS8AInp1sUlIiIiDhoO87S85bB5snMCBFBxxDyet9yauERERMSJkiBPslebPUAYDZysOZY122wnIiIillIS5EnHNtfvAXJiQEWe2U5EREQspSTIk07le7adiIiIeI2SIE8K6+XZdiIiIuI1SoI8KSbVnAWGrfE2QV2gx7VtFpKIiIg0zNIkaNOmTYwfP574+HhsNhsrV650Ol9YWMg999xDfHw84eHhjB07lv379zu1OX36NDNmzKB79+506dKFSZMmUVhY2Iafoo6AQHMaPNBoInS2DLb/VMXRIiIiFrM0CSovLyc5OZlFixbVO2cYBhMnTuTrr7/m3XffZefOnfTp04e0tDTKy8sd7ebMmcOqVatYtmwZGzdu5OjRo6SnW7gWT2I6pL4D4b2dj4cnwsUPgC0Avv47bLkDqs+YyVDhBji0xLxXciQiItImbIZhNDSfu83ZbDZWrFjBxIkTAdi3bx8DBgzg888/55JLLgHAbrcTFxfHH/7wB/7nf/6HkpISYmJieOutt5g8eTIAX331FYMGDSIzM5NrrrnGpfcuLS0lMjKSkpISIiIiPPOBGlsxOm85fHI72Ksg6gqoLIJTR869rrFFFbUCtYiIiJPW/n77bE1QZWUlAJ06dXIcCwgIIDQ0lI8//hiArKwsqqqqSEtLc7QZOHAgSUlJZGZmNnnt0tJSp5vHBQRC7Ajoe4d5X5uwJKbD9asgIARO7HROgKDhRRXzlsN7fWHdSNhyp3n/Xl8tvCgiItIKPpsE1SYz8+fP57vvvuPMmTM89dRTHD58mPx8c4p5QUEBISEhREVFOb02NjaWgoKCRq+9YMECIiMjHbfExERvfpT64tIgOLKRk+ctqqgVqEVERLzCZ5Og4OBgli9fzr59+4iOjiY8PJyMjAzGjRtHQEDrwp4/fz4lJSWOW15enoeidtGxzVB5rIkGNYsqrv8+ZE5FK1CLiIh4nk9voJqSkkJ2djYlJSWcOXOGmJgYhg0bxtChQwGIi4vjzJkznDhxwqk3qLCwkLi4uEavGxoaSmhoqLfDb5yriyUWZTTToM4K1LEjWhuViIiIX/HZnqC6IiMjiYmJYf/+/ezYsYMJEyYAZpIUHBzMunXrHG337t1Lbm4uw4cPtyrc5rm6WGJsWvNtQCtQi4iItIClPUFlZWUcOHDA8TwnJ4fs7Gyio6NJSkpi2bJlxMTEkJSUxO7du5k1axYTJ05k9OjRgJkcTZs2jblz5xIdHU1ERAQPPPAAw4cPd3lmmCVqF1WsOELDQ1028/zgX0Lhf5q/XqeaXi/NIBMREXGZpUnQjh07GDlypOP53LlzAZg6dSqLFy8mPz+fuXPnUlhYSK9evbj77rt55JFHnK7x3HPPERAQwKRJk6isrGTMmDH8+c9/btPP4bbaRRU3T8ZcVLFuIlSzyGLK8+YQV5PJUo0vn4LSvfDl750LqBubbi8iIiK+s06QlbyyTpAr8pZD1qzzEpdEMwGqTVxqZ4cB9ZMlA2xBYJxt5A1qEqrUd5QIiYhIh9Pa328lQViYBIFrQ1hNJUsRg+Cjy8F+ppE3qBlauyVHQ2MiItKhtPb326dnh/mF2kUVm5KYDr0nNJwsFW5oIgECzSATERFpmJKg9qKxZMnVmWGaQSYiIuKkXUyRlya4Ot3e1XYiIiJ+QklQe1c73b62CLoem1k/FJPallGJiIj4PCVB7V3tdHug4UTIMAuoVRQtIiLiRElQR5CYbk6DD+/dwMkA6HJRm4ckIiLi65QEdRSJ6XDLIRiVAd97y7xPmATYYfv9YNitjlBERMSnaHZYR3L+DLKu/aBgNRzfAl+/Dhfda1loIiIivkY9QR1ZeAJc+pj5OPthqCy2Nh4REREfoiSooxs4GyIHQ+Vx+OzXVkcjIiLiM5QEdXQBwTC0ZkPZA3+Bb7dbG4+IiIiPUBLkD2JvgL4/AgzY/jNzvzIRERE/p8Jof3HFn+DIKijOggOvQOQlTW/aKiIi0sEpCfIXYXEw5HeQ9QDseAAwzp0LTzAXXExMtyw8ERGRtqbhMH/SKa7mgeF8vOIIbJ4MecvbPCQRERGrKAnyF/Zq2DmnkZM1SVHWbNULiYiI31AS5C+ObYaKw000MKAiz2wnIiLiB5QE+YtT+Z5tJyIi0s4pCfIXYb08205ERKSdUxLkL2JSzVlg2BppYIPwRLOdiIiIH1AS5C8CAs1p8EDDiZBhriWk9YJERMRPKAnyJ4npkPoOhPc+70RNUnR8W5uHJCIiYhUtluhvEtOh9wRzFljtitFnTsDmW2HvcxA3CnrfbHWUIiIiXqckyB8FBELsCOdjFz8I+xbC1ntgXHYDvUUiIiIdi4bDxHTFH6HbFVB5HLb8SIsmiohIh6ckSEyBoXDtUgjqDEUb4IvfQ+EGOLTEvFdSJCIiHYyGw+SciIth6J9h61TY/ZjzOW2yKiIiHYx6gsRZcJeGj2uTVRER6WCUBMk59mrImtXISW2yKiIiHYuSIDlHm6yKiIgfURIk52iTVRER8SNKguQcbbIqIiJ+REmQnNPsJqtAWLw2WRURkQ5BSZCc0+wmq0BgFzCq2iwkERERb1ESJM4a22S1Uy9zIcWyfbB1GhiGNfGJiIh4iBZLlPoa2mQ1JhWKNkLGGPjmLYgcBJf+r9WRioiItJiSIGlYQ5usxt0IQ1+C7T+FXY9AxEBIuLV+shQQaEnIIiIi7lASJO7p/xMo3QN7X4AtUyA4CiqLzp3X9hoiItJOqCZI3HfF0xB1BdjPOCdAoO01RESk3VASJC1gg8rCRs5pew0REWkflASJ+45thlNHm2ig7TVERMT3KQkS92l7DRER6QCUBIn7tL2GiIh0AEqCxH3Nbq9hg/BEba8hIiI+TUmQuK/J7TVqnqc8r/WCRETEpykJkpZpbHuNkGjzuNYJEhERH6ckSFouMR1uOQSjMiD+ZvNYr9FKgEREpF3QitHSOrXbawQEw9EP4OiHYK8yn4uIiPgw9QSJZ3S/Bjr1hKoSc6NVERERH6ckSDwjIBB632I+zltpaSgiIiKuUBIknpMw0bw/vBIMw8pIREREmqUkSDwnbhQEdYFTR6A4y+poREREmqQkSDwnsBPEjzMfH15hbSwiIiLNUBIknlV3SExERMSHWZoEbdq0ifHjxxMfH4/NZmPlypVO58vKypg5cyYJCQmEhYUxePBgXnnlFac2p0+fZsaMGXTv3p0uXbowadIkCgsL2/BTiJP4m8AWBCVfQuk+q6MRERFplKVJUHl5OcnJySxatKjB83PnzmX16tX83//9H3v27GH27NnMnDmT9957z9Fmzpw5rFq1imXLlrFx40aOHj1KeroW67NMSBTEjjQfH37X0lBERESaYjMM35jGY7PZWLFiBRMnTnQcu/TSS7ntttt45JFHHMdSUlIYN24cv/vd7ygpKSEmJoa33nqLyZMnA/DVV18xaNAgMjMzueaaaxp8r8rKSiorKx3PS0tLSUxMpKSkhIiICO98QH+y78+wYwb0+B6M/sTqaEREpIMqLS0lMjKyxb/fPl0T9L3vfY/33nuPI0eOYBgGGRkZ7Nu3j9GjRwOQlZVFVVUVaWlpjtcMHDiQpKQkMjMzG73uggULiIyMdNwSExO9/ln8SkLNekHHM+FUgbWxiIiINMKnk6AXX3yRwYMHk5CQQEhICGPHjmXRokVcf/31ABQUFBASEkJUVJTT62JjYykoaPzHd/78+ZSUlDhueXl53vwY/ic8AaKvAgw4ssrqaERERBrk03uHvfjii2zdupX33nuPPn36sGnTJmbMmEF8fLxT74+7QkNDCQ0N9WCkUk/irVC8HfJWQL/7rI5GRESkHp9Ngk6dOsWvfvUrVqxYwc03mzuUDxkyhOzsbJ5++mnS0tKIi4vjzJkznDhxwqk3qLCwkLi4OIsiF8CcKv/Zr6BwHVSVQrBqrURExLf47HBYVVUVVVVVBAQ4hxgYGIjdbgfMIung4GDWrVvnOL93715yc3MZPnx4m8Yr54kYCF0vBvsZOLra6mhERETqsbQnqKysjAMHDjie5+TkkJ2dTXR0NElJSdxwww089NBDhIWF0adPHzZu3Mgbb7zBs88+C0BkZCTTpk1j7ty5REdHExERwQMPPMDw4cMbnRkmbcRmM3uD9vzRXDixzw+tjkhERMSJpVPkN2zYwMiRI+sdnzp1KosXL6agoID58+ezZs0aiouL6dOnD9OnT2fOnDnYbDbAXCxx3rx5LFmyhMrKSsaMGcOf//xnt4bDWjvFThpxLBPWfs8cCks/BoEhVkckIiIdSGt/v31mnSArKQnyEsMOK+LhdCGM/Df0Gm11RCIi0oF06HWCpJ2zBUDCBPNxnjZUFRER36IkSLwr4Vbz/si7Zs+QiIiIj1ASJN4VOxICu8CpfPj891C4AezVVkclIiLiu+sESQdx9AOgJunZ/ah5H54AKS9Aoja6FRER66gnSLwnbzlsngzVp5yPVxwxj+cttyYuERERlASJt9irIWsW0NDkw5pjWbM1NCYiIpZREiTecWwzVBxuooEBFXlmOxEREQsoCRLvOJXv2XYiIiIepiRIvCOsl2fbiYiIeJiSIPGOmFRzFhi2RhrYIDzRbCciImIBJUHiHQGB5jR4oNFEKOV5s52IiIgFlASJ9ySmQ+o7EN77vBOBcN0yrRMkIiKW0mKJ4l2J6dB7gjkLrPwb+O99YK+CiIutjkxERPyceoLE+wICIXYEXDgVYtPMY0c/sjQkERERJUHStuLHmff5q62NQ0RE/J6SIGlbtUnQsY+h6qS1sYiIiF9TEiRtq2s/6NLPrAsqWGd1NCIi4seUBEnbix9r3uerLkhERKyjJEjaXq+aIbGjH4HR0AarIiIi3qckSNpe7AgICDU3UC3dY3U0IiLip5QESdsLCjcTIdBUeRERsYySILFG3SExERERCygJEmvUFkcf2wxVZdbGIiIifklJkFij68XQ+QKwn4HC9VZHIyIifkhJkFjDZtPq0SIiYiklQWKdeE2VFxER6ygJEuvEjoSAECg/BKV7rY5GRET8jJIgsU5QZ+h5vflYq0eLiEgbUxIk1tJUeRERsYiSILFWbV1Q0UY4W25tLCIi4leUBIm1IgZC5z41U+U3WB2NiIj4ESVBYi2bDXrVLJyoITEREWlDSoLEeo71gjRVXkRE2o6SILFe7I0QEAxlX8PJ/VZHIyIifiLI3RecPn2aF198kYyMDIqKirDb7U7nP/30U48FJ34iuCvEpJrbZ+SvhoiLrY5IRET8gNtJ0LRp01izZg2TJ0/m6quvxmazeSMu8Tfx48wk6OhHMOBBq6MRERE/4HYS9P777/Phhx9y7bXXeiMe8Ve9xsLOh6BwHXz9ujljLCYVAgKtjkxERDoot5Og3r1707VrV2/EIv6sdC/YAsFeBVvvMY+FJ0DKC5CYbmloIiLSMbldGP3MM8/wi1/8gm+++cYb8Yg/ylsOH/8/MKqdj1ccgc2TzfMiIiIe5nZP0NChQzl9+jQXXngh4eHhBAcHO50vLi72WHDiB+zVkDULaGhqvAHYIGs29J6goTEREfEot5OgO+64gyNHjvCHP/yB2NhYFUZL6xzbDBWHm2hgQEWe2S52RFtFJSIifsDtJGjLli1kZmaSnJzsjXjE35zK92w7ERERF7ldEzRw4EBOnTrljVjEH4X18mw7ERERF7mdBD355JPMmzePDRs28O2331JaWup0E3FLTKo5C4zGhlVtEJ5othMREfEgm2G4t1lTQICZN51fC2QYBjabjerq6oZe5tNKS0uJjIykpKSEiIgIq8PxP3nLzVlggHOBdM2/sdR3NE1eRETqae3vt9s1QRkZGW6/iUiTEtPNRCdrlnORdEg0DPurEiAREfEKt5OgG264wRtxiL9LTDenwR/bDF/+0dxRPjFdCZCIiHiN20nQpk2bmjx//fXXtzgY8XMBgeY0+KoSMwk69rHVEYmISAfmdhI0YsSIesfq1ge1x5og8TEx15n3pXvgdBF06mltPCIi0iG5PTvsu+++c7oVFRWxevVqrrrqKtasWeONGMXfhHaHqMvMx0WbrY1FREQ6LLd7giIjI+sd+/73v09ISAhz584lKyvLI4GJn4u5Hk7shqKNkDTJ6mhERKQDcrsnqDGxsbHs3bvXU5cTfxdbU4Bf1HQNmoiISEu53RO0a9cup+eGYZCfn8+TTz7J5Zdf7qm4xN/VLo54Yhec+Q5Culkbj4iIdDhuJ0GXX345NpuN89dYvOaaa/j73//uscDEz4XFQcQAKN0LRR9DwnirIxIRkQ7G7eGwnJwcvv76a3JycsjJyeGbb76hoqKCLVu2MHDgQLeutWnTJsaPH098fDw2m42VK1c6nbfZbA3e/vSnPznaFBcXM2XKFCIiIoiKimLatGmUlZW5+7HEF8XULLdQtNHaOEREpENyOwnq06eP0y0xMZHTp0+36M3Ly8tJTk5m0aJFDZ7Pz893uv3973/HZrMxadK5QtkpU6bwxRdfsHbtWt5//302bdrE9OnTWxSP+JieqgsSERHvcXvvsKeeeoq+ffty2223AfDDH/6Qd955h169evHhhx+SnJzcskBsNlasWMHEiRMbbTNx4kROnjzJunXrANizZw+DBw9m+/btDB06FIDVq1dz0003cfjwYeLj4116b+0d5qPK8+DdJLAFwuTvILir1RGJiIgPae3vt9s9Qa+88gqJiYkArF27lrVr17J69WrGjRvHQw895HYAriosLOSDDz5g2rRpjmOZmZlERUU5EiCAtLQ0AgIC2LZtW6PXqqysdOx6X3sTH9Q5ETpfAEY1HNtidTQiItLBuJ0EFRQUOJKg999/nx/+8IeMHj2ahx9+mO3bt3s8wFqvv/46Xbt2JT393F5SBQUF9OzpvJpwUFAQ0dHRFBQUNHqtBQsWEBkZ6bjVfh7xQT1VFyQiIt7hdhLUrVs38vLyAHPoKS0tDTCnyntzy4y///3vTJkyhU6dOrX6WvPnz6ekpMRxq/084oNq64KOqS5IREQ8y+0p8unp6dx5553079+fb7/9lnHjxgGwc+dO+vXr5/EAATZv3szevXt5++23nY7HxcVRVFTkdOzs2bMUFxcTFxfX6PVCQ0MJDQ31SqziYbU9Qd/+F85WQFC4tfGIiEiH4XZP0HPPPcfMmTMZPHgwa9eupUuXLoA5k+v+++/3eIAAr732GikpKfWKrocPH86JEyectupYv349drudYcOGeSUWaWNdLoSw3mCvgm8br/MSERFxl9s9QcHBwfz85z+vd3zOnDluv3lZWRkHDhxwPM/JySE7O5vo6GiSkpIAs/J72bJlPPPMM/VeP2jQIMaOHct9993HK6+8QlVVFTNnzuT22293eWaY+DibzewN+mYJFG6E2JFWRyQiIh2Ex/YOa4kdO3ZwxRVXcMUVVwAwd+5crrjiCh599FFHm6VLl2IYBnfccUeD13jzzTcZOHAgo0aN4qabbuK6667jr3/9a5vEL21EdUEiIuIFbq8T1BFpnSAfV7IHPhgMgZ1g8gkIVD2XiIhYsE6QSJuLGAihMVB9Gop3WB2NiIh0EEqCxPfV1gWB1gsSERGPaVESdOLECV599VXmz59PcXExAJ9++ilHjhzxaHAiDtpHTEREPMzt2WG7du0iLS2NyMhIDh06xH333Ud0dDTLly8nNzeXN954wxtxir+r7Qk69gnYz0KA2/90RUREnLjdEzR37lzuuece9u/f77R680033cSmTfp/6eIlUZdBcBScLYPvdlodjYiIdABuJ0Hbt2/nJz/5Sb3jvXv3bnK/LpFWsQVAz1TzsYbERETEA9xOgkJDQxvcdX3fvn3ExMR4JCiRBjnqglQcLSIired2EnTLLbfwxBNPUFVVBYDNZiM3N5df/OIXTJo0yeMBijg4ZohtBsNubSwiItLuuZ0EPfPMM5SVldGzZ09OnTrFDTfcQL9+/ejatSu///3vvRGjiKnbFRDUBapOwIndVkcjIiLtnNtTbCIjI1m7di2ffPIJn332GWVlZVx55ZWkpaV5Iz6RcwKCIOZayP+3WRfULbn514iIiDTCrSSoqqqKsLAwsrOzufbaa7n22mu9FZdIw3reYCZBue9AaA8I6wUxqRAQaHVkIiLSzriVBAUHB5OUlER1dbW34hFpWu1Wd8c2ndtQNTwBUl6AxHTr4hIRkXbH7ZqgX//61/zqV79yrBQt0mbylsOu/61/vOIIbJ5snhcREXGR27vIX3HFFRw4cICqqir69OlD586dnc5/+umnHg2wLWgX+XbAXg3v9YWKw400sJk9QrfkaGhMRMRPtPb32+3C6IkTJ7r9JiKtdmxzEwkQgAEVeWa72BFtFZWIiLRjbidBjz32mDfiEGnaqXzPthMREb/Xol3kRdpcWC/PthMREb/ndk9QdXU1zz33HP/85z/Jzc3lzJkzTudVMC1eEZNq1vxUHAEaKWMLTzTbiYiIuMDtnqDHH3+cZ599lttuu42SkhLmzp1Leno6AQEB/OY3v/FCiCKYxc4pL9Q8sTXc5oKpKooWERGXuZ0Evfnmm/ztb39j3rx5BAUFcccdd/Dqq6/y6KOPsnXrVm/EKGJKTIfUdyC8t/PxwHDzfu/zULyzzcMSEZH2ye0kqKCggMsuuwyALl26UFJSAsAPfvADPvjgA89GJ3K+xHS45RCMyoDvvWXepx+H2FFwtgw23ATl31gdpYiItANuJ0EJCQnk55szcC666CLWrFkDwPbt2wkNDfVsdCINCQg0p8H3vcO8Dw6D1H9B1GVwugAyxkHlcSjcAIeWmPf2RlY5t1e71k5ERDoctwujb731VtatW8ewYcN44IEH+NGPfsRrr71Gbm4uc+bM8UaMIs0LiYQRH8Ka4VC6B1b0Bnudov2GttbIWw5Zs5zXH9IWHCIifsPtFaPPl5mZSWZmJv3792f8+PGeiqtNacXoDuSr5+HThpLxmmLq1HfMBCdvubnVRr2ZZue1ExERn9Xa3+9WJ0EdgZKgDqLZrTUwd56/+q+wbTqcOd5II23BISLSHrT5thlvvPFGk+fvvvtut4MQ8Yhmt9bArBXa3FwPz3lbcNirzcen8s3FGGNSlRyJiHQAbidBs2bNcnpeVVVFRUUFISEhhIeHKwkS67i6ZUZoDFQea75d8U44U6y6IRGRDsrt2WHfffed062srIy9e/dy3XXXsWTJEm/EKOIaV7fMuORXrrXbORc2T6rfu1RxxKwnylvuXnwiIuJTPLJ3WP/+/XnyySfr9RKJtKnarTUaW1Eam7m1Rr/7m2kH2Jpa7qGmjC5rtqbUi4i0Yx7bQDUoKIijR4966nIi7mtya42a5ynPQ1BIM+1scGlzvUV16oZERKRdcrsm6L333nN6bhgG+fn5vPTSS1x77bUeC0ykRWq31miwjuf5c3U8zbWrrnTt/VytQxIREZ/jdhI0ceJEp+c2m42YmBhuvPFGnnnmGU/FJdJyienQe0LzM7qaale4wbX3crUOSUREfI7bSZDdbvdGHCKeVbu1Rkvb1dYXVRyh/oKKdWifMhGRdstjNUEiHYor9UUAW++BL58Cw9A+ZCIi7YzbPUFz5851ue2zzz7r7uVFfEdTdUNXPgvfboM9T0P2L6FwI5zYDae0npCISHvhdhK0c+dOdu7cSVVVFQMGDABg3759BAYGcuWVVzra2WxNTD8WaS+aqhtKmgxhCfDpbMj/qP5ra9cT0j5kIiI+ye0kaPz48XTt2pXXX3+dbt26AeYCivfeey+pqanMmzfP40GKWKqp+qKLZ8LnT5grS9djADZzPaHeE7TVhoiIj3F7A9XevXuzZs0aLrnkEqfjn3/+OaNHj26XawVpA1VpscINsG5k8+1GZbhWqC0iIi5r7e+324XRpaWlHDtWf9+lY8eOcfLkSbcDEGnXXF0nSOsJiYj4HLeToFtvvZV7772X5cuXc/jwYQ4fPsy//vUvpk2bRnq66h7Ez7i6TpDWExIR8Tlu1wS98sor/PznP+fOO++kqqrKvEhQENOmTeNPf/qTxwMU8WmurCcUGAbdLm/LqERExAVu1wTVKi8v5+DBgwBcdNFFdO7c2aOBtSXVBEmr5C03Z4EBjSZCkYMhdSVE9DfXD2puNWsREWlWa3+/W5wE1Q1g/fr1DBw4kIEDB7bmUpZREiStlre8gfWEEqH//bDvRTh1FIIjzdlkOa83sF+Z1hMSEXFXmydBP/zhD7n++uuZOXMmp06dIjk5mUOHDmEYBkuXLmXSpEluB2E1JUHiEY318JzKN3uKjm9p5IU1a2ppPSEREbe0+eywTZs2kZqaCsCKFSswDIMTJ06wcOFCfve737kdgEiHUbueUN87zPvaIa6wXjDyPxDY2JBxzf8PyZqtrTZERNqQ20lQSUkJ0dHRAKxevZpJkyYRHh7OzTffzP79+z0eoEiHULwNqsubaGBARZ7ZkyQiIm3C7SQoMTGRzMxMysvLWb16NaNHjwbMVaM7derk8QBFOgStJyQi4nPcniI/e/ZspkyZQpcuXejTpw8jRowAzGGyyy67zNPxiXQMWk9IRMTnuJ0E3X///QwbNozc3Fy+//3vExBgdiZdeOGFqgkSaYwr6wmFJ5rtRESkTbR6inxHoNlh0iaaW09oyO/h0l+1aUgiIu1Zm88OE5EWSkw3p8GH93Y+Hhhm3u99Dk4ebPu4RET8lNvDYSLSConp0HuC83pC3a6E9TdCcRZsvBlGZ0JIN6sjFRHp8JQEibS12vWE6rphFfx7GJTuhU3pcMOH5rR6ba0hIuI1SoJEfEFYL7jhfVh7HRRtgOU9oLri3HltrSEi4nFu1wStXr2ajz/+2PF80aJFXH755dx555189913Hg1OxK90GwIDHjQf102AwJxVtnmyWVwtIiIe4XYS9NBDD1FaWgrA7t27mTdvHjfddBM5OTnMnTvXrWtt2rSJ8ePHEx8fj81mY+XKlfXa7Nmzh1tuuYXIyEg6d+7MVVddRW5uruP86dOnmTFjBt27d6dLly5MmjSJwsJCdz+WiPXs1ebmqg3S1hoiIp7mdhKUk5PD4MGDAfjXv/7FD37wA/7whz+waNEiPvroI7euVV5eTnJyMosWLWrw/MGDB7nuuusYOHAgGzZsYNeuXTzyyCNOK1PPmTOHVatWsWzZMjZu3MjRo0dJT9eQgbRDxzY77y5fj7bWEBHxJLdrgkJCQqioMLvq//Of/3D33XcDEB0d7eghctW4ceMYN25co+d//etfc9NNN/HHP/7Rceyiiy5yPC4pKeG1117jrbfe4sYbbwTgH//4B4MGDWLr1q1cc801bsUjYiltrSEi0qbc7gm67rrrmDt3Lr/97W/573//y8033wzAvn37SEhI8FhgdrudDz74gIsvvpgxY8bQs2dPhg0b5jRklpWVRVVVFWlpaY5jAwcOJCkpiczMzEavXVlZSWlpqdNNxHLaWkNEpE25nQS99NJLBAUF8c477/Dyyy/Tu7e58NtHH33E2LFjPRZYUVERZWVlPPnkk4wdO5Y1a9Zw6623kp6ezsaNGwEoKCggJCSEqKgop9fGxsZSUFDQ6LUXLFhAZGSk45aYmOixuEVarHZrDWyNNLBpaw0REQ9yezgsKSmJ999/v97x5557ziMB1bLb7QBMmDCBOXPmAHD55ZezZcsWXnnlFW644YYWX3v+/PlORdylpaVKhMR6AYHmNPjNkzETofO31jAg5XmtFyQi4iGt2jbj9OnTXhtW6tGjB0FBQY4i7FqDBg1yzA6Li4vjzJkznDhxwqlNYWEhcXFxjV47NDSUiIgIp5uIT2hsaw2AiMGQcGvbxyQi0kG5nQSVl5czc+ZMevbsSefOnenWrZvTzVNCQkK46qqr2Lt3r9Pxffv20adPHwBSUlIIDg5m3bp1jvN79+4lNzeX4cOHeywWkTaVmA63HIJRGfC9t+Dat8EWDKVfQsFaq6MTEekw3B4Oe/jhh8nIyODll1/mrrvuYtGiRRw5coS//OUvPPnkk25dq6ysjAMHDjie5+TkkJ2dTXR0NElJSTz00EPcdtttXH/99YwcOZLVq1ezatUqNmzYAEBkZCTTpk1j7ty5REdHExERwQMPPMDw4cM1M0zat/O31ji+Bfa+AJ/9GuK+D7bG6oZERMRlhpsSExONjIwMwzAMo2vXrsb+/fsNwzCMN954wxg3bpxb18rIyDAwCx+cblOnTnW0ee2114x+/foZnTp1MpKTk42VK1c6XePUqVPG/fffb3Tr1s0IDw83br31ViM/P9+tOEpKSgzAKCkpcet1Im2mosAwloYbxpsYRu4Kq6MREfEJrf39thmGcX71ZZO6dOnCl19+SVJSEgkJCSxfvpyrr76anJwcLrvsMsrKyjyeqHlbaWkpkZGRlJSUqD5IfFf2r+DLBRB5KYzLVoG0iPi91v5+u10TdOGFF5KTkwOYa/L885//BGDVqlX1pqqLiAcNfgiCI6Hkc8h92+poRETaPbeToHvvvZfPPvsMgF/+8pcsWrSITp06MWfOHB566CGPBygiNUK6waCfm493PQb2KmvjERFp59weDjvfN998Q1ZWFv369WPIkCGeiqtNaThM2o2qk/DehVB5HK7+G/T7H6sjEhGxTGt/v1udBHUESoKkXdnzLOycZ64ePX4/BIZaHZGIiCVa+/vt0hT5hQsXMn36dDp16sTChQubbPvggw+6HYSIuKH/z+CrZ80d5fe/DN0uNzdVDetlbqmhgmkREZe41BN0wQUXsGPHDrp3784FF1zQ+MVsNr7++muPBtgW1BMk7c7+v8D2n2KW9dnPHQ9PMLfeSEy3KjIRkTbTJj1BtbPBzn8sIhYJqV2d3e58vOKIufdY6jtKhEREmtGqvcMMw0AlRSJtzF5t1gQ1qOZ/j1mzzXYiItKoFiVBr732GpdeeimdOnWiU6dOXHrppbz66quejk1EGnJsM1QcbqKBYdYLHdtsPrVXQ+EGOLTEvFdyJCICtGDvsEcffZRnn33WsUcXQGZmJnPmzCE3N5cnnnjC40GKSB2n8l1vl7ccsmY5J02N1Q3Zq83ESUXWIuIn3J4iHxMTw8KFC7njjjucji9ZsoQHHniA48ePezTAtqDCaGlXCjfAupHNt+s5Eoo24Bgic6jZfLVu3ZA7yZKIiI9o820zqqqqGDp0aL3jKSkpnD171u0ARMRNMalmgkIzO8kXZVA/AYJ6dUN5y81i6vOH2GqLrPOWtz5mEREf5HYSdNddd/Hyyy/XO/7Xv/6VKVOmeCQoEWlCQKDZQwPUT4Rs5i2hud6bmrqhDT+AzKm4lCyJiHQwLtUEzZ071/HYZrPx6quvsmbNGq655hoAtm3bRm5uLnfffbd3ohQRZ4np5nBWg0NYz0N1JRx2oQenYHUzDeoUWceOaEXAIiK+x6UkaOfOnU7PU1JSADh48CAAPXr0oEePHnzxxRceDk9EGpWYDr0nNFzMXLjBtWv0HFFTN9QMV4uxRUTaEZeSoIyMDG/HISItERDYcA9Nbd1QxREaHuqymecv+V/XkqCwXq2LU0TEB7VqsUQR8VHN1g1hDpvFjmimyNpmbtQak+qNKEVELKUkSKSjqq0bCu/tfDw84dz0+CaTpRopz2u9IBHpkNxeLFFE2pGm6obqtmmoyBpgyO+0TpCIdFhKgkQ6usbqhuo6P1k69H9w9EM4kd0WEYqIWEJJkIiY6iZLkZeYSVDecijPg86JloYmIuINqgkSkfq6DTGnzxvVsH+R1dGIiHiFkiARadiAWeb9gb/C2QprYxER8QIlQSLSsN7jofMFcOY7s0ZIRKSDURIkIg0LCISLZ5qP9y4Eo6FFF0VE2i8lQSLSuIt+DEGdoeQLKFxndTQiIh6lJEhEGhcSBRfcYz7+6oWmWoqItDtKgkSkaQMeMO+PfgAnD1gbi4iIBykJEpGmRQyAXuMAA/a+aHU0IiIeoyRIRJpXO13+639AVam1sYiIeIiSIBFpXq/REDEQzp6ErxdbHY2IiEcoCRKR5tlsMOBB8/FXC6FgPRxaAoUbwF5taWgiIi2lvcNExDUX3A2f/hzKD8L6UeeOhydAygvabV5E2h31BImIa/L/DdUNbJ9RcQQ2TzY3WxURaUeUBIlI8+zVkDWrkZM1K0lnzdbQmIi0K0qCRKR5xzZDxeEmGhhQkWe2ExFpJ1QTJCLNO5XvXjt7tZkQncqHsF4Qk2ruRSYi4kOUBIlI88J6udau/BDk/gs+ne3cc6TiaRHxQTbD0NbQpaWlREZGUlJSQkREhNXhiPgeezW819csgqYl/8mwmXep7zgnQuoxEpFWaO3vt2qCRKR5AYFmTw7gSGgcbOat94QmLtBA8XTecjOxWjcSttxp3r/XV7PMRKTNKAkSEdckpps9OeG9nY+HJ5jHB85u5gI1xdOH3zUTnc2T6xdba7q9iLQh1QSJiOsS080en4aGsA4tce0aH08CWyAND6sZgM3sMeo9QUNjIuJVSoJExD0BgRA7ov5xV4unAYym1hOqM92+ofcREfEQDYeJiGfEpJpDY/VqhmrZIDwRrv6ra9dzdVq+iEgLKQkSEc9otngaSHkeuvZ37Xru9CyJiLSAkiAR8ZzmiqcT013vMYpJ9Xa0IuLnVBMkIp7VVPE0nOsx2jwZMxFqoEA65XkVRYuI1ykJEhHPa6x4ulZtj1HWrPrT5KNTtLK0iLQJJUEiYo3ze4zsZ2DrNCjeAUdXQ/xYqyMUkQ5OSZCIWOf8HqMTu+CrZ80eotjdEBhiWWgi0vGpMFpEfMdlj0GnWDi5D/Y+b3U0ItLBKQkSEd8RHAGXP2U+/vy3UHHU2nhEpENTEiQivuWCu6D7NXC2DLIftjoaEenAlASJiG+xBcBVLwE2OPQmFH1sdUQi0kEpCRIR3xOdAhf9j/l4x0zIX2du0Fq4AexN7TsmIuI6S5OgTZs2MX78eOLj47HZbKxcudLp/D333IPNZnO6jR3rPG22uLiYKVOmEBERQVRUFNOmTaOsrKwNP4WIeEXy7yEwHE58BhlpsOVOWDcS3usLecutjk5EOgBLk6Dy8nKSk5NZtGhRo23Gjh1Lfn6+47ZkyRKn81OmTOGLL75g7dq1vP/++2zatInp06d7O3QR8bZjm6G6ov7xiiPmatNKhESklSxdJ2jcuHGMGzeuyTahoaHExcU1eG7Pnj2sXr2a7du3M3ToUABefPFFbrrpJp5++mni4+MbfF1lZSWVlZWO56WlpS38BCLiFfZqc62gBhmADbJmm4stansNEWkhn68J2rBhAz179mTAgAH87Gc/49tvv3Wcy8zMJCoqypEAAaSlpREQEMC2bdsaveaCBQuIjIx03BITE736GUTETcc2199Ow4kBFXlmOxGRFvLpJGjs2LG88cYbrFu3jqeeeoqNGzcybtw4qqvNwsiCggJ69uzp9JqgoCCio6MpKCho9Lrz58+npKTEccvLy/Pq5xARN53K92w7EZEG+PS2Gbfffrvj8WWXXcaQIUO46KKL2LBhA6NGjWrxdUNDQwkNDfVEiCLiDWG9PNtORKQBPt0TdL4LL7yQHj16cODAAQDi4uIoKipyanP27FmKi4sbrSMSkXYgJhXCEwBb422CI6HHdW0Wkoh0PO0qCTp8+DDffvstvXqZ/+9v+PDhnDhxgqysLEeb9evXY7fbGTZsmFVhikhrBQRCygs1TxpJhKpK4NNZWjdIRFrM0iSorKyM7OxssrOzAcjJySE7O5vc3FzKysp46KGH2Lp1K4cOHWLdunVMmDCBfv36MWbMGAAGDRrE2LFjue+++/jvf//LJ598wsyZM7n99tsbnRkmIu1EYjqkvgPhvZ2PhyfChdMAG+z/M3xyO1SfNpOhwg3NL6roajsR6fBshmEYVr35hg0bGDlyZL3jU6dO5eWXX2bixIns3LmTEydOEB8fz+jRo/ntb39LbGyso21xcTEzZ85k1apVBAQEMGnSJBYuXEiXLl1cjqO0tJTIyEhKSkqIiIjwyGcTEQ+xV5uzwE7lmzVAMalmT9E3/4TMH4G9CiIGQ9UJOFVnw9XwBLM3KTH93LG85ebU+7ozzxpqJyLtQmt/vy1NgnyFkiCRdqpgPWy4GeynGzhZM4yW+o6Z4OQtNxdZ5Pz/5J3XTkTajdb+frermiARESc9b4Dgxv7DV5Ps7HgQKovN+3oJUJ12WbM1NCbiZ5QEiUj7dWwzVBY10cCAU0fgX93N+6baafFFEb+jJEhE2i9PL5aoxRdF/IqSIBFpv1xdLDH5Kc9eT0Q6BCVBItJ+Nbuoos2cUj9gdvOLL4Z0N68nIn5DSZCItF9NLqpY8zzleQgKaX7xxTPfQvbDYD+rtYRE/ISmyKMp8iLtXoPr/ySaCZAr6wRFXw2Hl5vPIy6BqmLn+iCtJSTik7ROkAcoCRLpABpbVNHVdnnL4ZMprq05JCI+obW/3z69i7yIiMsCAiF2RMvb9Z4AIZFwuqEkyABs5lpCvSecS65cTbxExCcpCRIRATOZOV3YRIM6awnFjtAWHCIdgAqjRUTA9TWCMu+GzZPMW90ECKDiiLk1R95yz8cnIh6nJEhEBFxfI6gir4kkR1twiLQnSoJERMC1NYfCesPgXzVzofO24NB0exGfpZogERE4t+bQ5smYiVDdibM1idHQhVBd6dr1Ko6obkjEx6knSESkVmK6OQ0+vLfz8fCEc9PjXR02+3SO6oZEfJzWCULrBInIeZqa+m6vhvf6mskMLf3Pp81MrG7J0ZR6kVZo7e+3eoJERM5Xu5ZQ3zvM+7qJSrNbddhg8Pxm3uC8uiERsYSSIBERdzU3bBZ1mWvXcXVavoh4hQqjRURaIjHdXD26oWGzwg2uXcPV+iIR8QolQSIiLdXYFhy10+2bqhuyBUJgmDejE5FmaDhMRMTTmq0bAoxq+M/1sP8vYBhaT0jEAuoJEhHxhtq6oYbWCUr+PeS+A0feg+0/NROfsgNw6ohzu4bWE9KmrSIeoynyaIq8iHhRY0mLYcCeP0H2fMDewAtreoxq1ycCLb4ocp7W/n4rCUJJkIhYxF4NK3pB5bFGGtRZT+jIuzWrWZ//n+wGkiX1FomfaO3vt4bDRESscmxzEwkQONYT+jAZynNouMjaAGzmpq29J5jJknqLRFyiwmgREau4uk5Q6RdQXdFEg5pk6Yvfm71F2qpDxCVKgkRErOLqOkEJt7rWbvfjNN5bhNlbVHfWmWakiZ/TcJiIiFWaXU+opiao/ww4vMKFCzZUYF2rzlYdsSNcL7JWfZF0YEqCRESsUrue0ObJmAXOdROhmoLnlOfNpKW5ZCk4Cqq+a/49j2dC5bfw8f+rf63aYbPaImvNRpMOTsNhIiJWam4fssR01xZfHDjbtff77Ffwye00O2z2zTvu1RdpaE3aIU2RR1PkRcQHuDLs1GDPTKLZW9R7ArzXt+mtOgLDoLoKONt8PCE94MzxRk7WmbofEKgeI7GM1gnyACVBItJuNJUs5S2vGVqDBofWUt+BMydg2zTPxDIqA84Uu75+kYiHaZ0gERF/0timrdD0Vh0pz5vnXd3h3hXb7quZ5u/C+kUqphYfpCRIRKQjSUw3k47GeotcmZEW2qOZRRxrlB1opsF5M9JEfIwKo0VEOpra3qK+d5j3dXthXCmyvmqRmSjVO1+nXac46D/TtXjqLgqpAmrxIeoJEhHxN64Mm9kCm566f9UiCImG/S81/361SZAKqMXHqDAaFUaLiJ9qbkZaU7PREtPN1zc3I61W5KVQ8nkDJ1RALS2n2WEeoCRIRKQRriRKTc1I6zUWCv4DRlUTb3LelHsRF7X291s1QSIi0rim6oug+cUeR34I33ujmTepU0At0oZUEyQiIq3T3Iw0Vwcc6hZQi7QBJUEiItJ6Ta1fFNbLtWvUttOmrdJGlASJiIh3Nbs2EYANvtsFp4/DzjmaQSZtQjVBIiLiXa6sTYQBn86CT/6f65u2irSSkiAREfG+pgqor1sGKS/S+OKMdXa31+KK4kEaDhMRkbbRVAF14QaaXmuogS04VDskraQkSERE2k5jBdSuzgxzd/VpJUrSBCVBIiJiPVdnkB38u1kjlP0w9XqOamuHalefdmebDiVLfklJkIiIWM+lGWRA4X/MW4MMwFZTO2SHT35Y/1rnJ0qgPc38mAqjRUTEes3OILPB5U9B7KhmLlRTO7R1Kg0nU+cVWddu++HKjDR7tVm7dGiJed9Ykbar7cRy6gkSERHf4Mru9uGJULiu+WtVVzRxsiZR2jgejm+h8WSpplep9wQ48q5rvUXqVWpXtIEq2kBVRMSnNFWfU7gB1o1s23guexx2/4b6yVJNj1XdGqTNk5tvV8uVOiTVKjVJu8h7gJIgEZF2wl4N7/VtonbIBqE9oPJY89eKSXVt01ZbEBhnGztp9vT84AC8f1H9YbXz292SYyYxrvQYqVepWdpFXkRE/Icrq09ftchMFhpdfNFmDqtd+phr79loAgSOobVVTSVAddoVrnetDsmdWiVwrQ5JtUr1qCcI9QSJiLQ7DfaSJJ6rHXIMTYFzj1GdoaneE5rvVQqOgKoSz8ZuCwSjiQQkJAZsNqgsauwC3u1VakdDcO26J2jTpk2MHz+e+Ph4bDYbK1eubLTtT3/6U2w2G88//7zT8eLiYqZMmUJERARRUVFMmzaNsrIy7wYuIiLWSkyHWw7BqAz43lvm/S05537Qm9qmo7Y2x5VepYFzXYvnovtcj72pBAjgzLEmEiBwWj3b071KecvNxHDdSNhyp3n/Xt+G923rAL1Pls4OKy8vJzk5mR//+Mekpzc+vrlixQq2bt1KfHx8vXNTpkwhPz+ftWvXUlVVxb333sv06dN56623vBm6iIhYrbHVp2s1tU1H3TZNzUjrPQEO/q3p3qLwBEh5CfI/ar7d4F/Cjhkt+bT1fXw7nD3ZyPvVzG7bMavO80ba1J0B11Bhd0vXVmoHNU0+Mxxms9lYsWIFEydOdDp+5MgRhg0bxr///W9uvvlmZs+ezezZswHYs2cPgwcPZvv27QwdOhSA1atXc9NNN3H48OEGk6aGaDhMRMTPNTUE5MrQmqtDcCHRbT+7zRU9R8K325pYWqDOEFxjyVLdzwnuzZRroXY9HNYcu93OXXfdxUMPPcQll1xS73xmZiZRUVGOBAggLS2NgIAAtm3b1uh1KysrKS0tdbqJiIgfq+1V6nuHed9Qb1FTQ2uutqtdGbupou2wBPPWZJveMOjn7n7KxhVluLa20ofJsOUumlyI8r8zYPuMptvULlZpMZ9eLPGpp54iKCiIBx98sMHzBQUF9OzZ0+lYUFAQ0dHRFBQUNHrdBQsW8Pjjj3s0VhER6cBcGVpzpV1tHdLmyZhJTgM9RkNr6pSabLPQ7FXa87RnPl/c96FgbfPtSr9opoEBlY3//jra1NY0NTWc2QZ8ticoKyuLF154gcWLF2OzNZYNt8z8+fMpKSlx3PLy8jx6fRER6YCa6i1yp50rPUZt2asUngiDftH4564r4VbX2rniVL7nrtVCPtsTtHnzZoqKikhKSnIcq66uZt68eTz//PMcOnSIuLg4ioqcK+jPnj1LcXExcXFxjV47NDSU0NBQr8UuIiLSJFeLttuiVynleTNZa3ID25qaoP4z4PCK1n56U1gvz1ynFXw2CbrrrrtIS0tzOjZmzBjuuusu7r33XgCGDx/OiRMnyMrKIiUlBYD169djt9sZNmxYm8csIiLisuZmt7nSxpX91sC1Ns0lVK4mS2E1vVenmkmoYlIb/1xtxNIkqKysjAMHDjie5+TkkJ2dTXR0NElJSXTv3t2pfXBwMHFxcQwYMACAQYMGMXbsWO677z5eeeUVqqqqmDlzJrfffrvLM8NERETaNU/0KtW28USy5Grvkw8swGhpErRjxw5Gjjw3VXDuXHNRqqlTp7J48WKXrvHmm28yc+ZMRo0aRUBAAJMmTWLhwoXeCFdERMQ3eaJXCTybLLnSxmI+s06QlbROkIiIiJtc2V7Dy1twtPb322drgkRERMSHear3yUI+O0VeRERExJuUBImIiIhfUhIkIiIifklJkIiIiPglJUEiIiLil5QEiYiIiF9SEiQiIiJ+SUmQiIiI+CUlQSIiIuKXtGI0ULtzSGlpqcWRiIiIiKtqf7dbugOYkiDg5MmTACQmJlociYiIiLjr5MmTREZGuv06baAK2O12jh49SteuXbHZbC69prS0lMTERPLy8rTpqgX0/VtL37+19P1bR9+9tc7//g3D4OTJk8THxxMQ4H6Fj3qCgICAABISElr02oiICP0PwUL6/q2l799a+v6to+/eWnW//5b0ANVSYbSIiIj4JSVBIiIi4peUBLVQaGgojz32GKGhoVaH4pf0/VtL37+19P1bR9+9tTz9/aswWkRERPySeoJERETELykJEhEREb+kJEhERET8kpIgERER8UtKglpg0aJF9O3bl06dOjFs2DD++9//Wh1Sh7Rp0ybGjx9PfHw8NpuNlStXOp03DINHH32UXr16ERYWRlpaGvv377cm2A5owYIFXHXVVXTt2pWePXsyceJE9u7d69Tm9OnTzJgxg+7du9OlSxcmTZpEYWGhRRF3LC+//DJDhgxxLAo3fPhwPvroI8d5ffdt58knn8RmszF79mzHMX3/3vWb3/wGm83mdBs4cKDjvKe+fyVBbnr77beZO3cujz32GJ9++inJycmMGTOGoqIiq0PrcMrLy0lOTmbRokUNnv/jH//IwoULeeWVV9i2bRudO3dmzJgxnD59uo0j7Zg2btzIjBkz2Lp1K2vXrqWqqorRo0dTXl7uaDNnzhxWrVrFsmXL2LhxI0ePHiU9Pd3CqDuOhIQEnnzySbKystixYwc33ngjEyZM4IsvvgD03beV7du385e//IUhQ4Y4Hdf3732XXHIJ+fn5jtvHH3/sOOex798Qt1x99dXGjBkzHM+rq6uN+Ph4Y8GCBRZG1fEBxooVKxzP7Xa7ERcXZ/zpT39yHDtx4oQRGhpqLFmyxIIIO76ioiIDMDZu3GgYhvl9BwcHG8uWLXO02bNnjwEYmZmZVoXZoXXr1s149dVX9d23kZMnTxr9+/c31q5da9xwww3GrFmzDMPQv/228NhjjxnJyckNnvPk96+eIDecOXOGrKws0tLSHMcCAgJIS0sjMzPTwsj8T05ODgUFBU5/i8jISIYNG6a/hZeUlJQAEB0dDUBWVhZVVVVOf4OBAweSlJSkv4GHVVdXs3TpUsrLyxk+fLi++zYyY8YMbr75ZqfvGfRvv63s37+f+Ph4LrzwQqZMmUJubi7g2e9fG6i64fjx41RXVxMbG+t0PDY2lq+++sqiqPxTQUEBQIN/i9pz4jl2u53Zs2dz7bXXcumllwLm3yAkJISoqCintvobeM7u3bsZPnw4p0+fpkuXLqxYsYLBgweTnZ2t797Lli5dyqeffsr27dvrndO/fe8bNmwYixcvZsCAAeTn5/P444+TmprK559/7tHvX0mQiDRrxowZfP75505j8uJ9AwYMIDs7m5KSEt555x2mTp3Kxo0brQ6rw8vLy2PWrFmsXbuWTp06WR2OXxo3bpzj8ZAhQxg2bBh9+vThn//8J2FhYR57Hw2HuaFHjx4EBgbWq0AvLCwkLi7Ooqj8U+33rb+F982cOZP333+fjIwMEhISHMfj4uI4c+YMJ06ccGqvv4HnhISE0K9fP1JSUliwYAHJycm88MIL+u69LCsri6KiIq688kqCgoIICgpi48aNLFy4kKCgIGJjY/X9t7GoqCguvvhiDhw44NF//0qC3BASEkJKSgrr1q1zHLPb7axbt47hw4dbGJn/ueCCC4iLi3P6W5SWlrJt2zb9LTzEMAxmzpzJihUrWL9+PRdccIHT+ZSUFIKDg53+Bnv37iU3N1d/Ay+x2+1UVlbqu/eyUaNGsXv3brKzsx23oUOHMmXKFMdjff9tq6ysjIMHD9KrVy/P/vtvRfG2X1q6dKkRGhpqLF682Pjyyy+N6dOnG1FRUUZBQYHVoXU4J0+eNHbu3Gns3LnTAIxnn33W2Llzp/HNN98YhmEYTz75pBEVFWW8++67xq5du4wJEyYYF1xwgXHq1CmLI+8YfvaznxmRkZHGhg0bjPz8fMetoqLC0eanP/2pkZSUZKxfv97YsWOHMXz4cGP48OEWRt1x/PKXvzQ2btxo5OTkGLt27TJ++ctfGjabzVizZo1hGPru21rd2WGGoe/f2+bNm2ds2LDByMnJMT755BMjLS3N6NGjh1FUVGQYhue+fyVBLfDiiy8aSUlJRkhIiHH11VcbW7dutTqkDikjI8MA6t2mTp1qGIY5Tf6RRx4xYmNjjdDQUGPUqFHG3r17rQ26A2nouweMf/zjH442p06dMu6//36jW7duRnh4uHHrrbca+fn51gXdgfz4xz82+vTpY4SEhBgxMTHGqFGjHAmQYei7b2vnJ0H6/r3rtttuM3r16mWEhIQYvXv3Nm677TbjwIEDjvOe+v5thmEYHuipEhEREWlXVBMkIiIifklJkIiIiPglJUEiIiLil5QEiYiIiF9SEiQiIiJ+SUmQiIiI+CUlQSIiIuKXlASJiIiIX1ISJCJtbsSIEcyePdvqMJzYbDZWrlxpdRgi0oa0YrSItLni4mKCg4Pp2rUrffv2Zfbs2W2WFP3mN79h5cqVZGdnOx0vKCigW7duhIaGtkkcImK9IKsDEBH/Ex0d7fFrnjlzhpCQkBa/Pi4uzoPRiEh7oOEwEWlztcNhI0aM4JtvvmHOnDnYbDZsNpujzccff0xqaiphYWEkJiby4IMPUl5e7jjft29ffvvb33L33XcTERHB9OnTAfjFL37BxRdfTHh4OBdeeCGPPPIIVVVVACxevJjHH3+czz77zPF+ixcvBuoPh+3evZsbb7yRsLAwunfvzvTp0ykrK3Ocv+eee5g4cSJPP/00vXr1onv37syYMcPxXiLi+5QEiYhlli9fTkJCAk888QT5+fnk5+cDcPDgQcaOHcukSZPYtWsXb7/9Nh9//DEzZ850ev3TTz9NcnIyO3fu5JFHHgGga9euLF68mC+//JIXXniBv/3tbzz33HMA3HbbbcybN49LLrnE8X633XZbvbjKy8sZM2YM3bp1Y/v27Sxbtoz//Oc/9d4/IyODgwcPkpGRweuvv87ixYsdSZWI+D4Nh4mIZaKjowkMDKRr165Ow1ELFixgypQpjjqh/v37s3DhQm644QZefvllOnXqBMCNN97IvHnznK75v//7v47Hffv25ec//zlLly7l4YcfJiwsjC5duhAUFNTk8Ndbb73F6dOneeONN+jcuTMAL730EuPHj+epp54iNjYWgG7duvHSSy8RGBjIwIEDufnmm1m3bh333XefR74fEfEuJUEi4nM+++wzdu3axZtvvuk4ZhgGdrudnJwcBg0aBMDQoUPrvfbtt99m4cKFHDx4kLKyMs6ePUtERIRb779nzx6Sk5MdCRDAtddei91uZ+/evY4k6JJLLiEwMNDRplevXuzevdut9xIR6ygJEhGfU1ZWxk9+8hMefPDBeueSkpIcj+smKQCZmZlMmTKFxx9/nDFjxhAZGcnSpUt55plnvBJncHCw03ObzYbdbvfKe4mI5ykJEhFLhYSEUF1d7XTsyiuv5Msvv6Rfv35uXWvLli306dOHX//6145j33zzTbPvd75BgwaxePFiysvLHYnWJ598QkBAAAMGDHArJhHxXSqMFhFL9e3bl02bNnHkyBGOHz8OmDO8tmzZwsyZM8nOzmb//v28++679QqTz9e/f39yc3NZunQpBw8eZOHChaxYsaLe++Xk5JCdnc3x48eprKysd50pU6bQqVMnpk6dyueff05GRgYPPPAAd911l2MoTETaPyVBImKpJ554gkOHDnHRRRcRExMDwJAhQ9i4cSP79u0jNTWVK664gkcffZT4+Pgmr3XLLbcwZ84cZs6cyeWXX86WLVscs8ZqTZo0ibFjxzJy5EhiYmJYsmRJveuEh4fz73//m+LiYq666iomT57MqFGjeOmllzz3wUXEcloxWkRERPySeoJERETELykJEhEREb+kJEhERET8kpIgERER8UtKgkRERMQvKQkSERERv6QkSERERPySkiARERHxS0qCRERExC8pCRIRERG/pCRIRERE/NL/B0MbdX/d2VodAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "costs = cmeans.costs\n",
    "iters = np.arange(1, len(costs) + 1)\n",
    "plt.plot(iters, costs, color='orange', marker='o')\n",
    "plt.xlabel(\"iteration\")\n",
    "plt.ylabel(\"bias square sum\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "从图像可以看出每次迭代误差平方和都在下降。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 二、ISODATA方法\n",
    "\n",
    "ISODATA方法能一定程度上自动优化聚类的数量，是C均值的改进版，比C均值要灵活一些。\n",
    "\n",
    "这个算法的问题在于参数非常多，而且都是需要按照经验给出的。书上也没给出实际应用中可行的参数大小，实现这个算法需要大量时间来试错调参，因此这个算法我就不实现了。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 三、基于样本与核的相似性度量的动态聚类算法\n",
    "\n",
    "这个算法是C均值算法的衍生，定义了两个函数：类核函数和距离函数，实现算法就是通过类和函数计算出每个类的核心，再通过距离函数计算每个样本到每个类核心的距离，把所有样本都放到距离最小的类，然后更新每个类的核心。反复这样迭代直到收敛。\n",
    "\n",
    "C均值算法实际上就是这个方法的一个特例，C均值算法的类核函数是类均值，距离函数是欧氏距离的平方。\n",
    "\n",
    "这个方法要解决的问题是样本在各个维度上的方差不统一，也就是各个维度度量大小不统一，但是其实我们在前面实现C均值算法时通过归一化处理已经把各个维度度量统一了，所以这个算法其实意义不大，这个算法也不实现了。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 四、模糊C均值算法\n",
    "\n",
    "模糊C均值算法用$\\displaystyle J_e = \\sum_{j=1}^c\\sum_{i=1}^n[\\mu_j(\\mathbf{x}_i)]^b||\\mathbf{x}_i-\\mathbf{m}_j||^2$替代C均值算法中的误差平方和函数，其中$b$通常取$2$，$\\mu_j(\\mathbf{x})$是$\\mathbf{x}$对$j$类的隶属度。\n",
    "\n",
    "下面直接用改进的模糊C均值算法的约束条件：$\\displaystyle\\sum_{j=1}^c\\sum_{i=1}^n\\mu_j(\\mathbf{x}_i)=n$。\n",
    "\n",
    "这个条件下推导出的算法如下：\n",
    "\n",
    "1. 获得初始的聚类中心，计算方法可以和C均值初始化方式一样，也可以用C均值计算的均值结果。\n",
    "2. 用下面的公式计算每个样本的隶属度函数：\n",
    "\n",
    "$$\n",
    "\\mu_j(\\mathbf{x}_i) = \\frac{\\displaystyle n\\left(\\frac{1}{||\\mathbf{x}_i-\\mathbf{m}_j||^2}\\right)^{\\frac{1}{b-1}}}{\\displaystyle\\sum_{k=1}^c\\sum_{l=1}^n\\left(\\frac{1}{||\\mathbf{x}_l-\\mathbf{m}_k||^2}\\right)^{\\frac{1}{b-1}}},\\quad i=1,2,\\cdots,n,\\quad j=1,2,\\cdots,c\n",
    "$$\n",
    "\n",
    "3. 用下面的公式更新各类聚类中心：\n",
    "\n",
    "$$\n",
    "\\mathbf{m}_j = \\frac{\\displaystyle\\sum_{i=1}^n[\\mu_j(\\mathbf{x}_i)]^b\\mathbf{x}_i}{\\displaystyle\\sum_{i=1}^n[\\mu_j(\\mathbf{x})_i]^b},\\quad j=1,2,\\cdots,c\n",
    "$$\n",
    "\n",
    "反复迭代2、3步直到收敛。\n",
    "\n",
    "实现代码如下："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "class FuzzyCMeans:\n",
    "    \"\"\"\n",
    "    这个类用于实现FCM算法\n",
    "    \"\"\"\n",
    "    def __init__(self, c: int, b: int):\n",
    "        \"\"\"\n",
    "        构造函数\n",
    "        :param c: 要聚类的数量\n",
    "        :param b: 公式中的b\n",
    "        \"\"\"\n",
    "        self._c: int = c\n",
    "        # self._c是聚类的数量\n",
    "\n",
    "        self._b: int = b\n",
    "        # self._b记录公式中的b\n",
    "\n",
    "        self._X: np.ndarray = np.zeros(0)\n",
    "        # self._X存储样本矩阵\n",
    "\n",
    "        self._N: int = 0\n",
    "        # self._N存储样本数量\n",
    "\n",
    "        self._d: int = 0\n",
    "        # self._d存储样本的维度\n",
    "\n",
    "        self.m: np.ndarray = np.zeros(0)\n",
    "        # self._m存储各类的聚类中心\n",
    "\n",
    "        self.mu: np.ndarray = np.zeros(0)\n",
    "        # self._mu存储各个样本的隶属度\n",
    "\n",
    "        self.costs: list[float] = []\n",
    "        # self.cost记录损失函数的值\n",
    "\n",
    "    def initialize(self, _X: np.ndarray):\n",
    "        self._X = _X.copy()\n",
    "        self._N = _X.shape[0]\n",
    "        self._d = _X.shape[1]\n",
    "        self.m = np.zeros((self._c, self._d))\n",
    "        self.mu = np.zeros((self._N, self._c))\n",
    "        # 用C均值算法确定初始值，用来加速收敛\n",
    "        cmeans = CMeans(self._c)\n",
    "        cmeans.initialize(self._X)\n",
    "        cmeans.cluster()\n",
    "        for i in range(self._c):\n",
    "            self.m[i] = cmeans.m[i]\n",
    "\n",
    "    def _compute_cost(self) -> float:\n",
    "        \"\"\"\n",
    "        计算损失函数\n",
    "        :return: 损失函数的值\n",
    "        \"\"\"\n",
    "        ret = 0\n",
    "        for i in range(self._N):\n",
    "            for j in range(self._c):\n",
    "                ret += self.mu[i][j] ** self._b * np.sum((self._X[i] - self.m[j]) ** 2)\n",
    "        return ret\n",
    "\n",
    "    def _update_mu(self):\n",
    "        \"\"\"\n",
    "        更新隶属度mu\n",
    "        \"\"\"\n",
    "        inf = 999999\n",
    "        for i in range(self._N):\n",
    "            for j in range(self._c):\n",
    "                _sum = 0\n",
    "                for l in range(self._N):\n",
    "                    for k in range(self._c):\n",
    "                        square_bias = np.sum((self._X[l] - self.m[k]) ** 2)\n",
    "                        if square_bias != 0:\n",
    "                            _sum += (1 / square_bias) ** (1 / (self._b - 1))\n",
    "                        else:\n",
    "                            _sum += inf\n",
    "                square_bias = np.sum((self._X[i] - self.m[j]) ** 2)\n",
    "                if square_bias != 0:\n",
    "                    self.mu[i][j] = self._N * ((1 / square_bias) ** (1 / (self._b - 1))) / _sum\n",
    "                else:\n",
    "                    self.mu[i][j] = self._N * inf / _sum\n",
    "\n",
    "    def _update_m(self):\n",
    "        \"\"\"\n",
    "        更新各类聚类中心m\n",
    "        \"\"\"\n",
    "        for j in range(self._c):\n",
    "            self.m[j] = np.sum(((self.mu[:, j].reshape(self._N, 1) ** self._b) * self._X), 0) / np.sum(self.mu[:, j] ** self._b)\n",
    "\n",
    "    def cluster(self, _iters: int):\n",
    "        for _iter in range(_iters):\n",
    "            self._update_mu()\n",
    "            self._update_m()\n",
    "            _cost = self._compute_cost()\n",
    "            self.costs.append(_cost)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "从书中第203页的图片可以发现，改进的模糊C均值算法的优势在于就算实际的类别数量和给定的类别数量不符，算法也能得到较好的结果，我们在下面测试一下。\n",
    "\n",
    "首先生成数据，显然生成的数据实际上有四个聚类。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGdCAYAAAAvwBgXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA0PUlEQVR4nO3df3Bb1Zn/8Y+sNCJQy4RggomUOEC3wNICm0DoUG/tklnS6R9mhZNuwjKE4ct36YStXXamC7TTtLNtXRpa7DIdoO0szXSbQOqKZdsu2TIhZr1dmLKkMJCGtPmSTIxiQkqDlWU6Csj6/nF7HcuW9cO6V+f+eL9mNIqurn2PhQY9Oud5nhMpFAoFAQAAGNBkegAAACC8CEQAAIAxBCIAAMAYAhEAAGAMgQgAADCGQAQAABhDIAIAAIwhEAEAAMbMMz2AciYmJnTkyBE1NzcrEomYHg4AAKhCoVDQiRMndN5556mpqfych6cDkSNHjiiZTJoeBgAAmIPR0VElEomy53g6EGlubpZk/SHxeNzwaAAAQDWy2aySyeTk53g5ng5E7OWYeDxOIAIAgM9Uk1ZBsioAADCGQAQAABhDIAIAAIwhEAEAAMYQiAAAAGMIRAAAgDEEIgAAwBgCEQAAYIynG5oBAFBWPi+NjEhjY1Jbm9TRIUWjpkeFGhCIAAD8KZ2Wenul118/dSyRkAYHpVTK3LhQE5ZmAAD+k05LPT3FQYgkZTLW8XTazLhQMwIRAIC/5PPWTEihMPM5+1hfn3UePI9ABADgLyMjM2dCpioUpNFR6zx4HoEIAMBfxsacPQ9GEYgAAPylrc3Z82AUgQgAwF86OqzqmEik9PORiJRMWufB8whEAAD+Eo1aJbrSzGDEfjwwQD8RnyAQAQD4TyolDQ1JS5YUH08krOP0EfENGpoBAPwplZK6u+ms6nMEIgAA/4pGpc5O06NAHVxfmslkMvrbv/1bLVq0SAsWLNCHPvQh/c///I/blwUAAD7g6ozI8ePHdc0116irq0tPPvmkWltb9bvf/U4LFy5087IAAMAnXA1E7r33XiWTST3yyCOTx5YvX+7mJQEAgI+4ujTzb//2b1q5cqXWrl2rc845R1dccYW+973vzXp+LpdTNpstugEAgOByNRB57bXX9OCDD+oDH/iA/uM//kOf/vSn9ZnPfEZbt24teX5/f79aWlomb8lk0s3hAQAAwyKFQqntC50xf/58rVy5Uv/93/89eewzn/mMnn/+eT377LMzzs/lcsrlcpOPs9msksmkxsfHFY/H3RomAABwUDabVUtLS1Wf367OiLS1temSSy4pOnbxxRfr8OHDJc+PxWKKx+NFNwAAEFyuBiLXXHON9u/fX3Tst7/9rZYtW+bmZQEAgE+4Goh89rOf1XPPPaevfe1rOnDggLZt26bvfve72rRpk5uXBQAAPuFqIHLllVfq8ccf1/bt23XppZfqn/7pnzQwMKAbb7zRzcsCAACfcDVZtV61JLsAAABv8EyyKgAAQDkEIgAAwBgCEQAAYAyBCAAAMIZABAAAGEMgAgAAjCEQAQAAxhCIAAAAYwhEAACAMQQiAADAGAIRAABgDIEIAAAwhkAEAAAYQyACAACMIRABAADGEIgAAABj5pkeAAAAnpfPSyMj0tiY1NYmdXRI0ajpUQUCgQgAAOWk01Jvr/T666eOJRLS4KCUSpkbV0CwNAMAwGzSaamnpzgIkaRMxjqeTpsZV4AQiAAAUEo+b82EFAozn7OP9fVZ52HOCEQAAChlZGTmTMhUhYI0OmqdhzkjEAEAoJSxMWfPQ0kEIgAAlNLW5ux5KIlABACAUjo6rOqYSKT085GIlExa52HOCEQAACglGrVKdKWZwYj9eGCAfiJ1IhABAGA2qZQ0NCQtWVJ8PJGwjtNHpG40NAMAoJxUSuruprOqSwhEAACoJBqVOjtNjyKQWJoBAADGMCMCBFyhkNfbb4/o5MkxzZ/fpjPP7FAkwpQyAG8gEAEC7NixtA4c6FUud6o7ZCyW0IUXDqq1lSQ7AOaxNAME1LFjae3d21MUhEhSLpfR3r09OnaMzboAmEcgAgRQoZDXgQO9kkps1vWnYwcO9KlQYLMuAGYRiAAB9PbbIzNmQooVlMuN6u232awLaKh8XhoelrZvt+7ZuZccESCITp6sbhOuas8D4IB0WurtLd7RN5GwureGuDEaMyJAAM2fX90mXNWeB6BO6bTU01MchEhSJmMdT4c3Z4tAxAGFQl7Hjw/r6NHtOn58mHV3GHfmmR2KxRKSZtmsSxHFYkmdeSabdQGuy+etmZBCiZwt+1hfX2iXaViaqRPlkfCiSCSqCy8c1N69PbKCkan/A7SCkwsvHKCfCNAIIyMzZ0KmKhSk0VHrvBB2b2VGpA6UR8LLWltT+vM/H1IsVrxZVyyW0J//+RCBMtAoY1XmYlV7XsAwIzJHlcsjIzpwoE9nn93Nt04Y09qa0tlnd9NZFTCprcpcrGrPCxgCkTmqpTxy4cLOBo0KmCkSifIeBEzq6LCqYzKZ0nkikYj1fEc4c7ZYmpkjyiMBAFWJRq0SXckKOqayHw8MWOeFEIHIHIW1PJIKIQCYg1RKGhqSlhTnbCmRsI6HuI8ISzNzZJdH5nIZlc4TiSgWSwSqPJIKIQCoQyoldXdb1TFjY1ZOSEdHaGdCbA2bEfn617+uSCSivr6+Rl3SVXZ55J8eTX9WUrDKI6kQAgAHRKNWie769dZ9yIMQqUGByPPPP6+HH35YH/7whxtxuYYJS3kkG6gBANzi+tLM//7v/+rGG2/U9773PX3lK19x+3INF4bySCqEalMo5AP9fgAAJ7keiGzatEmf/OQntXr16oqBSC6XUy6Xm3yczWbdHp4jgl4eSYVQ9cijAeCYfD4U+SSuLs08+uij2rNnj/r7+6s6v7+/Xy0tLZO3ZDLp5vBQpbBWCNWKPBoAjkmnpfZ2qatL2rDBum9vD+TmeK4FIqOjo+rt7dWPfvQjnXbaaVX9zN13363x8fHJ2+joqFvDQw3YQK0y8mgAOCZkO/W6Foi88MILevPNN/UXf/EXmjdvnubNm6dnnnlG3/72tzVv3jzlS+wyGIvFFI/Hi24wL2wVQnNRSx4NAMwqhDv1uhaIXHvttXr55Zf14osvTt5WrlypG2+8US+++KKiAVznCrKwVAjNFXk0ABxRy069AeFasmpzc7MuvfTSomNnnHGGFi1aNOM4/CEMFUJzRR4NAEeEcKdeOquiJkGvEJqrMHbaBeCCuezU6/PqmoYGIsPDw428HNAwdh7N3r09svJmpgYj5NEAqFKtO/Wm01ZOydTlnETC2mTPJ/vXsOkd4BDyaADUrZadegNSXRMpFEqFXN6QzWbV0tKi8fFxKmjgG3RWBVC3UjMdyaQVhKRS1nJMe/vsia32zMnBg0aWaWr5/CYQAQDAi8rlfgwPW03OKtm929pcr8Fq+fwmWTVk+LYOAD5h79RbyhNPVPc7fFBdQyASIuyDAgABkE5bSzTVqLYKxyCSVUOCfVAAIADszquVRCJWTkmH91sGEIiEAPugAEBAVOq8aisUTlXXeByBSAiwDwqA0MrnrcTO7dute7/v0VJtzkdfn2/6iJAjEgLsgwLAl+rtGBqAZl8zVJvz0d3t7jgcxIxICLAPCgDfSaetPhldXdKGDdZ9e3v1TbqcbPblpVkVu/Pq9GZnNh/lhtgIRELA3gfFbjU+U0SxWJJ9UAB4Q71BhJ3QWapNln2sr6+6gKLegMhptXRe9QkCkRCw90H506Ppz0ry9j4ohUJex48P6+jR7Tp+fJikWiDInAgiKiV0FgrS6Kh1XjlebaGeSklDQ9KS4u0klEhYx3227ESOSEjY+6CU7iMy4Nk+IvQ+AUKmliBitmZf1SZ0ljuvUkAUiVgBUXe3mdmHVMq6to933bURiBjWyE6nra0pnX12t286q9q9T6aXHdu9T9hIDgggJ4KIahM6y53nREDktnKdV32EQMQgE9/2I5GoFi7sdOV3O6ly75OIDhzo09lnd3s2kAIwB04EEXZCZyZTekbD3hCuXEKnEwERqkKOiCF0Oi2P3idASDlRFeJEQqcTARGqQiBiAJ1OK6P3CRBSTlWF1JvQGcAyWa8iEHFZqYoPvu1XRu8TIMScqgpJpaRDh6Tdu6Vt26z7gwer+/kAlsl6FTkiLpotB6S1taeqnw/zt32790kul1HpmaOIYrEEvU+AoHKqKqSehE47ICrVnXVgwHdlsl4VKRRKZfJ4QzabVUtLi8bHxxWPx00PpyazVXxYfTuqe8kvu2y3LxJL3XLqNZSKXzPr2whVMwAaot5W8yFUy+c3MyIuqCYHRIpKmpjlHL7tS/7tfQIgYJwqkyWgKYlAxAWVc0AkyU5EnT5D4v1Op43kt94nQdTIXjdAYAVxAz6HEIi4oNrcjiVL+vT73w/xbb8Cv/Q+CSI62wIOsFvFT8+EsFvF+7Atu5PIEXHB8ePDeumlrornXXbZbp15ZgffNuFJ5fOcyNEBqpLPW5vkzdal1W6udvBgoJZpavn8pnzXBbXsdmt/21+8eL0WLuwkCIEn0OsGcIhTG/AFGIGIC/y+2y1ArxuETj4vDQ9L27db9+V2960FreIrIhBxiV3xEYsVN+SJxRJMaYdUqeZ2XkVnW4RKOm0tn3R1SRs2WPft7dbxetEqviKSVV1ExUe4lKsu8VvSJ51tERpuJ5I6sQFfwJGsCjigXKAhyXdJn4VCXs89116xs+3VVx8ksIZ/NSqR1A52pOJgxG4VH8CqGZJV0VB+WnJwQ/mdlG/Q/v3/V35L+iTPCaHQqERSp/bOCSiWZlAXvy05OK2a6pL33nur3G+YTPr0Wq8UOtuiLn7oItrIRFKn9s4JIAIRzNlsfSasmYAeTy45OK26LrqVeTXpkzwnzIlfuog2OpHUqVbxAcPSDOaEPhMWpwIILyd90usGNbHzIaYvedjJn05UojjFTiSNzNLzKRKRkslQJ5I2AoEI5oQ+E5b6A4hTze0A38vnrZmQUjUQ9rG+Pud6dNQrGrVmaaSZwYj9eGCA5ROXEYhgTugzYammi+68eYv+9DxJnwg4P3YRJZHUOHJEMCfVzgS8733nuDwSs+zqEitXpvROyh/84HcliaRPBJ9fu4iGLZHUY4nEBCKYE3smYPY+E5ZXX92oD3wg2BU01VaXkPSJwPNzF9GwJJJ6MJGYhmaYs1NVM9LswYh3m3Y5rVxnVSAU7AZhlbqIBmynWd+YrYusC43VaGiGhrBnAubPP6/MWeGpoKG6BKFH8qd3eTiRmEAEdWltTenii7dWOCscFTQARPKnV3k4kZgcEdTt5Mk3qzzPYwlqANwRtuRPP/BwIjGBCOrGTq0AZghL8qdfeDiRmEAEdatcQWPt1ErTLgAVlSotlZhdqZfdRbZSIrGBLrIEIqhbNb00aNoFoKJSpaWLFln3b03ZPNKL+9Z4nZ1I3NNjBR1TgxHDicSuJqv29/fryiuvVHNzs8455xxdf/312r9/v5uXhCF2BU0sVpygFoslQlG6C6BOs+1R89ZbxUGI5M19a/zAo4nErvYRWbNmjf7mb/5GV155pd577z3dc889euWVV/Sb3/xGZ5xxRsWfp4+I/9BLA0DN7P4j5ao6pqulJ4nHOoka14DXo5bP74Y2NDt27JjOOeccPfPMM/rLv/zLiucTiABACAwPS11dc/vZ3bvLJ8V6sJNoGHi2odn4+Lgk6ayzzir5fC6XUzabLboBAAKunpLRcj8723IPSzue0rBAZGJiQn19fbrmmmt06aWXljynv79fLS0tk7dkMtmo4cHHCoW8jh8f1tGj23X8+HBDOriauCYQWPWUjM72sx7uJIpiDVua+fSnP60nn3xS//Vf/6VEIlHynFwup1wuN/k4m80qmUyyNINZHTuWnmWzOfc22jNxTSDQKu1RU0qlHJFdu6TVqyv/nkpLO5gTzy3N3HHHHfrZz36m3bt3zxqESFIsFlM8Hi+6AbOxN92bGhBIUi6X0d69PTp2zPlpVxPXBAKv3B41pVQqN02npXXrqru2gU6iKOZqIFIoFHTHHXfo8ccf19NPP63ly5e7eTmESKGQ14EDvSrdQK0gqaDf/vZ2TUycbOA1w7G5H+CK2UpLFy061UvEVq7c1M4L+cMfqruugU6iKOZqQ7NNmzZp27ZteuKJJ9Tc3Kw33nhDktTS0qIFCxa4eWkE3Ntvj8yYlZju3XeP6dlnE/qzP3vIkSWTytc8tbnfwoWddV8PCJ3Z9qiRqis3LZcXMp3BTqIo5mog8uCDD0qSOqetvz3yyCPauHGjm5dGwFW7gd677x7T3r09jjRVq/aabO4H1GG2PWqqyeOotMPsdIY6iaKYq4FIA1uUIGRq3UDvwIE+nX12d13N1djcD/Cgqc25fvOb6n5m0SLpu9+lj4hHhHKvGbp/+l/ljfamcmbJ5NQ1Ky8JVcJ7EHBAqWZl1XjsMenaa90ZE2oWukCE0stgKN5orzr1LplEIlFdcMH9+s1v1pY97//9v39Qa2tq1sCC9yDgADsptZaZdzsvhHJdT2loZ1XTKL0MFnujvfe97+yqzndiyaSaa9mzL6XwHgQcUEtSqs3wDrNl5fNWm/vt2637kDVZC00gQullMLW2pvSRj2T0vve1ljkrolgsqTPPrD87vp6EVd6DgENqTUqVjO8wO6t02mrm1tUlbdhg3be3h6r9fGgCkVpKL+EvTU3z9Wd/9pCkyJ9uU1mPL7xwwJEcjHoSVnkPAg6ptgnZF74gbdtmdU89eNCbQQh74YQnEKH0MtjsZZpYrLgZUiyWcKR012YnrM4MeGyzz77wHgQcUm0Tsmuvldavt3JCvLgcw144kkKUrErpZfC1tqZ09tndrlajFCfJRlS8zFJ+9oX3IOCQjg5rqWW2vWn80Kys0vJSoSCNjlrnBTy5NjQzIvV8k4VZtex0G4lEtXBhpxYvXq+FCztdKYmd6+wL70HAIeX2pvFyUupU1S4vhWAvnNDMiNTzTRbmeLXUdS6zL7wHAQfZe9NM7yOSSFhBiNfyQaardnkpBHvhRAoebn9ayzbC1Sr9wZbUhRcO0MPBY+xS15lVJtaHtpO5H43EexCBNLXDabn9YIJy3Xrl81Z1TKXlpYMH/fH3TFPL53foAhGJrpZ+UCjk9dxz7WWqTCKKxRK6+uqDvvxvx3sQgVKqw2kiYS2feH1mwiS7akYqDkbs5SUvlhtXiUAEvnf8+LBeeqmr4nmXXbabnW4Bk2brcBqAD9OGKBXEJZP+WF4qo5bP79DkiMBfKHUFfKBSCWokYpWgdnf7cnmhIVIp6/Xx4/KSQwhE4EmUugI+QAmqM6LRUL8+oSnfhb9Q6gr4ACWocACBCDzJLnX906Ppz0qi1BUwjhJUOIBABJ7VqLbtAObI7nA6vamYLRKxEi+93OEUxpEjAk9rRNt2AHNkdzjt6bGCjlIlqF7vcArjmBGB5zWibTuAObI7nC4pnrlUIkHpLqrCjAgAoD6UoKIOBCIAgPqFvAQVc8fSDAAAMIZABAAAGEMgAgAAjCEQAQAAxhCIAAAAYwhEAACAMQQiAADAGAIRAABgDIEIAAAwhkAEAAAYQyACAACMIRABAADGEIgAAABjCEQAAIAxBCIAAMAYAhEAAGAMgQgAADCGQAQAABhDIAIAAIwhEAEAAMYQiAAAAGMIRAAAgDEEIgAAwJiGBCLf+c531N7ertNOO02rVq3Sr371q0ZcFgAAeJzrgchjjz2mO++8U5s3b9aePXt02WWX6brrrtObb77p9qUBAIDHuR6IfOtb39Jtt92mW265RZdccokeeughnX766frnf/5nty8NAAA8ztVA5OTJk3rhhRe0evXqUxdsatLq1av17LPPzjg/l8spm80W3QAAQHC5Goj8/ve/Vz6f1+LFi4uOL168WG+88caM8/v7+9XS0jJ5SyaTbg4PAAAY5qmqmbvvvlvj4+OTt9HRUdNDAgAALprn5i8/++yzFY1GdfTo0aLjR48e1bnnnjvj/Fgsplgs5uaQAACAh7g6IzJ//nytWLFCu3btmjw2MTGhXbt26SMf+YiblwYAAD7g6oyIJN155526+eabtXLlSl111VUaGBjQO++8o1tuucXtS8Pj8nlpZEQaG5Pa2qSODikaNT0qAEAjuR6IfOpTn9KxY8f0xS9+UW+88YYuv/xy7dy5c0YCK8IlnZZ6e6XXXz91LJGQBgelVMrcuAAAjRUpFAoF04OYTTabVUtLi8bHxxWPx00PBw5Jp6WeHmn6Oy8Sse6HhghGAMDPavn89lTVDIIvn7dmQkqFv/axvj7rPABA8BGIoKFGRoqXY6YrFKTRUes8AEDwEYigocbGnD0PAOBvrierAlO1tTl7HoC5y0/kNXJ4RGMnxtTW3KaOpR2KNlG6hsYiEEFDdXRY1TGZTOk8kUjEer6jo/FjA8IkvS+t3p29ej17aq00EU9ocM2gUheTLY7GYWkGDRWNWiW60qkqGZv9eGCAfiKAm9L70urZ0VMUhEhSJptRz44epfelDY0MYUQggoZLpawS3SVLio8nEpTuAm7LT+TVu7NXBc2ckrSP9e3sU36C0jU0BkszMCKVkrq76awKNNrI4ZEZMyFTFVTQaHZUI4dH1Nne2biBIbQIRGBMNCp1dpoeBRAuYyeqK0mr9jygXizNAECItDVXV5JW7XlAvQhEACBEOpZ2KBFPKKJIyecjiigZT6pjKaVraAwCEQAIkWhTVINrrNK16cGI/XhgzQD9RNAwBCIAEDKpi1MaWjekJfHi0rVEPKGhdUP0EUFDsfsuAIQUnVXhllo+v6maAYCQijZFKdGFcSzNAAAAYwhEAACAMQQiAADAGAIRAABgDIEIAAAwhkAEAAAYQyACAACMIRABAADGEIgAAABjCEQAAIAxBCIAAMAYAhEAAGAMgQgAADCG3XcBD8vnpZERaWxMamuTOjqkKLu0AwgQAhHAo9JpqbdXev31U8cSCWlwUEqlzI0LAJzE0gw8K5+Xhoel7dut+3ze9IgaJ52WenqKgxBJymSs4+m0mXEBgNMIROBJ6bTU3i51dUkbNlj37e3h+ADO562ZkEJh5nP2sb6+cAVmTshP5DV8aFjbX96u4UPDyk/wAgJewNIMPMeeDZj+QWzPBgwNBXtpYmRk5kzIVIWCNDpqndfZ2bBh+Vp6X1q9O3v1evbUC5uIJzS4ZlCpiwP8ZjIoP5HXyOERjZ0YU1tzmzqWdijaRIITZmJGBJ7CbICVmOrkeWGX3pdWz46eoiBEkjLZjHp29Ci9LwTTbA2W3pdW+2C7urZ2aUN6g7q2dql9sJ3XGiURiMBTapkNCKq2NmfPC7P8RF69O3tV0MzI1j7Wt7OPZRoHEfihVgQi8BRmA6wS3URCikRKPx+JSMmkdR7KGzk8MuMDcaqCChrNjmrkcIAj2wYi8MNcEIj4QJiqR5gNsPqEDA5a/54ejNiPBwboJ1KNsRPVRazVnofyghb4keDcGCSrelzYeknYswGZTOk8kUjEej7oswGplJWUW+q//cBAMP/bu6GtubqItdrzUF6QAj8SnBuHGREPC2MvCWYDTkmlpEOHpN27pW3brPuDBwlCatGxtEOJeEIRlV7niiiiZDypjqUBj2wbJCiBX6PyXJhxsUQKhVLfO70hm82qpaVF4+PjisfjpofTUPm81TdjtsRNe2bg4MFgfiiXmglKJpkNQO3sDxVJRbkLdnAytG6Ib7gOyU/k1T7Yrkw2UzJPJKKIEvGEDvYe9Gwpr/03zLbE5NTfEPQZl1o+v5kR8aiwV48wGwCnpC5OaWjdkJbElxQdT8QTBCEOizZFNbjGmtKcPgtlPx5YM+DZIERqTJ4LlUXFyBHxKKpHrJkeGnbBCamLU+r+YDcNthrADvxKfdsfWDPg+cDP7TyXSpVFEUXUt7NP3R/sDs37k0DEo6geAZwVbYqqs73T9DBCwc+Bn9t5LrXMuITl/Uog4lFUjwDwM78GfnaCc6U8l7kmOAepssgpruWIHDp0SLfeequWL1+uBQsW6IILLtDmzZt18uRJty4ZKFSPAEDjuZ3nEpTKIie5Foi8+uqrmpiY0MMPP6y9e/fq/vvv10MPPaR77rnHrUsGjt1LYklxjp0SieBv/AYApriZ4ExJ+UwNLd/dsmWLHnzwQb322mtVnR/m8t2p8nmrOmZszMoJ6ehgJgQA3ObWDsJhKCmv5fO7oTki4+PjOuuss2Z9PpfLKZfLTT7OZrONGJbnUT0CAI3nVp6L3yuLnNawQOTAgQN64IEHdN999816Tn9/v7785S83akgAABjh58oip9W8NHPXXXfp3nvvLXvOvn37dNFFF00+zmQy+tjHPqbOzk59//vfn/XnSs2IJJNJXyzNsHyCRuL9Bpjn1tJNENSyNFNzIHLs2DG99dZbZc85//zzNX/+fEnSkSNH1NnZqauvvlo/+MEP1NRUfX6sX3JEwrYxHczi/QaYF/QW7fVyNRCpRSaTUVdXl1asWKF/+Zd/UbTGr2x+CETsjemmv4p2iS3VLXAS7zfAPDvZdHqfkSAlm9bLE4FIJpNRZ2enli1bpq1btxYFIeeee25Vv8PrgUjYN6YLG9PLIbzfAPMatSme33li07unnnpKBw4c0K5du5RIJNTW1jZ5C4qwb0wXJum0FQR0dUkbNlj37e3W8Ubh/QaY14hN8cLGtUBk48aNKhQKJW9BwcZ04WAvh0wPAjIZ63ijghHeb4B5tGh3nmuBSBiwMV3w5fNWYmip+Nk+1tdnnec23m+AebRodx6BSB3sjemm7wVji0SkZJKN6fzMS8shvN8A82jR7jwCkTqwMd1M+bw0PCxt327dN2KmwE1eWg7h/QaYV25TPFs9m+KFEYFIndiY7hQvJHQ6zWvLIbzfAPPsFu1nLZi5ZUmpYyivoZve1crr5btTmS7tNG22/ha2H//Yer6USq+dydfWLpnNZEr/baZKZsP+fgNMS+9L64YdN8w4Ti8Riyf6iDjBT4FImFXqbyFZH5Lbt0tr1xYfr9Ql1AtdRO0gSyoORmgiBoQTvUQq80QfEYRHpYROyQpW1q0rXqapVBb7uc/NrWzW6TwVlkMATEUvEWcRiKButSRq2qWu1ZTFfutbtZfNupWnkkpJhw5Ju3dL27ZZ9wcPEoQAYUQvEWfNMz0A+F8tiZpTS10rlcWWm8mYWjbb2Wkdmy1PxZ5BqXf2Iho9dS0A4UUvEWcxI4K62f0tqjU25ly5q/17vNR4DECwOdVLJD+R1/ChYW1/ebuGDw0rPxHO/0ERiKBuU/tbVKOtzblyV/v3eKnxGIBgK9dLxH5cqZdIel9a7YPt6trapQ3pDera2qX2wXal9/m438EcEYjAEamUVaJbroR0aufParqEVvu7JG81HgMQfHYvkSXx4iz2RDxRsXQ3vS+tnh09MxJeM9mMenb0zCkY8fPsCjkicExPj1Wpsm7dzOdKdf4cHLR+JhIpXRZ7553SffdZ/y71/NTf5bXGYwCCL3VxSt0f7NbI4RGNnRhTW3ObOpZ2lJ0JyU/k1buzVwXNXEcuqKCIIurb2afuD3ZXXfqb3pdW787eosAmEU9ocM2gL3qZMCMCR61dK/3kJzNzRkqVulYqi/3GN6ovmzW1D0vQWtoDqE20KarO9k6t/9B6dbZ3VgwenC79dWN2pdFoaAZX1NL506nOqo1uPOZUszW6pMLr8hP5mr71Y3bbX96uDekNFc/bltqm9R9aX/YcLzdWq+Xzm6UZuKKWUtdK51b7u+wZllLBwcCA80GIE6XCXugcC5Tj92n/ckwEWE6W/tYyu9LZ3lntEBuOGREEjtszDJVa2le7/8xswUypGRxmTWCCPe0/PZ8hCPupmAqw7FmMTDZTMk+kllkMJ2dXnEaLd4SaPYOyfr117/QHthOlwrX0PQnirsbwvkpJlZLUt7PPV9UZNpN5FU6U/tqC0liNQASokROlwtUGM1/96tz22wHqFdT9VLwQYNVT+juVU43VTCNHBKiRE6XC1QYzg4Ozz5pEItasSXc3yzRwXlD3U/FKXsVcSn+ns2dXenb0KKJIUXBV6+yKScyIADVyolS42mDmD3+Y/Tm6xcJNQZn2n85LAVatpb+lODW7YhIzIkCN7Jb25ZqxTW22VoodzGQypWc8IhFp4cLygYiNbrFwgz3tXymp0uvT/tMFMcByYnbFJGZEgDmo1IytUunt1P15ps+s2I97e6sbC91i4QYnkyq9JCh5FdM5MbtiCoEIMEeplHTokLR7t7Rtm3V/8GD1/T8qBTOf/7yZbrGALQjT/tMFNcDyM/qIAIaV6xHS6G6xQClB7Kxaqo9IMp7UwJoBXwZYXlPL5zeBiE/Q0Cq8SnVfTSad7xYLhE0QAyyvIBAJGNqAg0AUteADFqYRiARILW3A4W0EE2gEJ1uXE9BgrghEAsKpPU1gHrNaaAQn94YJ8mZ39SA4qw6BSEAMD1v7ilSye3f1O906jW/5lTGrhUZwckv4IG92Vw+Cs+qx6V1AOLGniZvYjK2yWja3A+rh1N4wXtiLxYtMbpQXdAQiHubEniZusb/lsxlbeU7s1AtUw6nW5W5tdpefyGv40LC2v7xdw4eGfRXIEJy5ixbvHlZNG/BEovENrSp9y2cztlO8PquF4HCqdXktAU21+RJ+X9LwykZ5QcWMiIdV0wa80p4mbuBbfvW8PKuFYHGqdXm1Ac3v/vA7tQ+2q2trlzakN6hra5faB9tnLFEEYUnDSxvlBRGBiMfVu6eJG/iWXz0nduoFquFU6/JqAppFCxbpS8Nfqhhc1LKk4eWlmyBulOclBCI+UO+eJk7jW371vDqrhWByYm+YSgGNHUBUE1xUu6Tx1ZGvVjW7YkpQN8rzCsp3UTO7v0ml3BX6m5xCm3Y0khO9Lmbbi+X//MX/0ebhzRV/fvfNuzV2Ykwb0htqHr/kvVJhe4lJKg7CvDZOr6CPCFzHZmy1o+cK/KZUQLNj746qgottqW1qa25T19YqmiHNopbeJ43ARnnVq+Xzm6oZzImdu1KqWyjf8kuLRs01ngPmItoUnVEFUku+hL2kkclmSi7lVOK1apTUxSl1f7CbzqoOIxDBnKVSVoku3/KB8KgUXNizGPYH9OCaQfXs6CnKL7HPqzY48VI1SqngzK+80q6eZFXUxf6Wv369dU8QAgRbrdU55RJov9z55aquSTWK89L70p5JECZHBABQs1rzJUp9+5ak9sH2irMrXskRCYpG7CVEsioAwHVOVedQjdI4Tm6OWA6b3gEAXGfnS6z/0Hp1tnfO6YPLid4nqJ5bewnVoyHJqrlcTqtWrdJLL72kX//617r88ssbcVkAgA9QjdI4XmxX35BA5HOf+5zOO+88vfTSS424HADAZ4JUjeJlXmxX7/rSzJNPPqlf/OIXuu+++9y+FAAAKMOL7epdDUSOHj2q2267TT/84Q91+umnVzw/l8spm80W3QAAgDOc2hzRSa4FIoVCQRs3btTtt9+ulStXVvUz/f39amlpmbwlk0m3hgcAQCh5LUG45vLdu+66S/fee2/Zc/bt26df/OIX2rFjh5555hlFo1EdOnRIy5cvL5usmsvllMvlJh9ns1klk0nKdwEAcJibnVVd7SNy7NgxvfXWW2XPOf/887Vu3Tr99Kc/VWTK3uf5fF7RaFQ33nijtm7dWvFa9BEBAMB/PNHQ7PDhw0U5HkeOHNF1112noaEhrVq1SolEouLvIBABAMB/PLH77tKlS4sev//975ckXXDBBVUFIQAAIPjorAoAAIxpSEMzSWpvb5eHt7UBAAAGNCwQARotn5dGRqSxMamtTerokKJ0jAYATyEQQSCl01Jvr/T6lL2dEglpcFBKsYcWAHgGOSIInHRa6ukpDkIkKZOxjqfTZsYFAJiJQASBks9bMyGl0pHsY3191nkAAPNYmkGgjIzMnAmZqlCQRkelBx6QFi8mdwQATCMQQaCMjVV33mc/e+rf5I4AgDkszSBQ2tpq/xlyRwDAHAIRBEpHhzXDEYlUPtdG7ggAmEMggkCJRq1lFqn2YGR01MoxAQA0DoEIAieVkoaGpCVLav/ZanNMAADOIBBBIKVS0qFD0u7d0rZt0v33V/dzc8kxAQDMHVUzCKxoVOrstP6dz0vf/KaVmFqqx0gkYuWWdHQ0dIgAEHrMiCAUyuWO2I8HBugnAgCNRiCC0JgtdySRsI7TRwQAGo+lGYRKKiV1d7MrLwB4BYEIQmdq7ggAwCyWZgAAgDHMiAAAECL5ibxGDo9o7MSY2prb1LG0Q9Emc+vTBCIAAIREel9avTt79Xr21DbliXhCg2sGlbrYTMY+SzMAAIRAel9aPTt6ioIQScpkM+rZ0aP0PjM7fxKIAAAQcPmJvHp39qqgmR0d7WN9O/uUn2j8zp8EIgAABNzI4ZEZMyFTFVTQaHZUI4cbv/MngQgAAAE3dqK6HT2rPc9JBCIAAARcW3N1O3pWe56TCEQAAAi4jqUdSsQTiihS8vmIIkrGk+pY2vidPwlEAAAIuGhTVINrrJ0/pwcj9uOBNQNG+okQiAAAEAKpi1MaWjekJfHinT8T8YSG1g0Z6yMSKRQKM2t5PCKbzaqlpUXj4+OKx+OmhwMAgO81orNqLZ/fdFYFACBEok1RdbZ3mh7GJJZmAACAMQQiAADAGAIRAABgDIEIAAAwhkAEAAAYQyACAACMIRABAADGEIgAAABjCEQAAIAxdFYFAPhKI1qUo3EIRAAAvpHel1bvzl69nn198lgintDgmkFjm7ahPizNAAB8Ib0vrZ4dPUVBiCRlshn17OhRel/a0MhQDwIRAIDn5Sfy6t3Zq4JmbhhvH+vb2af8RL7RQ0OdCEQAAJ43cnhkxkzIVAUVNJod1cjhkQaOCk4gEAEAeN7YiTFHz4N3uBqI/PznP9eqVau0YMECLVy4UNdff72blwMABFRbc5uj58E7XKua+clPfqLbbrtNX/va1/Txj39c7733nl555RW3LgcACLCOpR1KxBPKZDMl80QiiigRT6hjaYeB0aEergQi7733nnp7e7Vlyxbdeuutk8cvueQSNy4HAAi4aFNUg2sG1bOjRxFFioKRiCKSpIE1A/QT8SFXlmb27NmjTCajpqYmXXHFFWpra9MnPvGJijMiuVxO2Wy26AYAgCSlLk5paN2QlsSXFB1PxBMaWjdEHxGfcmVG5LXXXpMkfelLX9K3vvUttbe365vf/KY6Ozv129/+VmeddVbJn+vv79eXv/xlN4YEAAiA1MUpdX+wm86qAVLTjMhdd92lSCRS9vbqq69qYmJCkvT5z39eN9xwg1asWKFHHnlEkUhEP/7xj2f9/XfffbfGx8cnb6Ojo/X9dQCAwIk2RdXZ3qn1H1qvzvZOghCfq2lG5B/+4R+0cePGsuecf/75Ghuzyqem5oTEYjGdf/75Onz48Kw/G4vFFIvFahkSAADwsZoCkdbWVrW2tlY8b8WKFYrFYtq/f78++tGPSpLeffddHTp0SMuWLZvbSAEAQOC4kiMSj8d1++23a/PmzUomk1q2bJm2bNkiSVq7dq0blwQAAD7kWh+RLVu2aN68ebrpppv0xz/+UatWrdLTTz+thQsXunVJAADgM5FCoTCzM4xHZLNZtbS0aHx8XPF43PRwAABAFWr5/GavGQAAYAyBCAAAMIZABAAAGEMgAgAAjHGtasYJdh4te84AAOAf9ud2NfUwng5ETpw4IUlKJpOGRwIAAGp14sQJtbS0lD3H0+W7ExMTOnLkiJqbmxWJWNs8Z7NZJZNJjY6OUtLrAl5fd/H6uovX1128vu4K0utbKBR04sQJnXfeeWpqKp8F4ukZkaamJiUSiZLPxeNx3/+H8jJeX3fx+rqL19ddvL7uCsrrW2kmxEayKgAAMIZABAAAGOO7QCQWi2nz5s2KxWKmhxJIvL7u4vV1F6+vu3h93RXW19fTyaoAACDYfDcjAgAAgoNABAAAGEMgAgAAjCEQAQAAxgQiEPn5z3+uVatWacGCBVq4cKGuv/5600MKnFwup8svv1yRSEQvvvii6eEEwqFDh3Trrbdq+fLlWrBggS644AJt3rxZJ0+eND003/rOd76j9vZ2nXbaaVq1apV+9atfmR5SIPT39+vKK69Uc3OzzjnnHF1//fXav3+/6WEF0te//nVFIhH19fWZHkrD+D4Q+clPfqKbbrpJt9xyi1566SX98pe/1IYNG0wPK3A+97nP6bzzzjM9jEB59dVXNTExoYcfflh79+7V/fffr4ceekj33HOP6aH50mOPPaY777xTmzdv1p49e3TZZZfpuuuu05tvvml6aL73zDPPaNOmTXruuef01FNP6d1339Vf/dVf6Z133jE9tEB5/vnn9fDDD+vDH/6w6aE0VsHH3n333cKSJUsK3//+900PJdD+/d//vXDRRRcV9u7dW5BU+PWvf216SIH1jW98o7B8+XLTw/Clq666qrBp06bJx/l8vnDeeecV+vv7DY4qmN58882CpMIzzzxjeiiBceLEicIHPvCBwlNPPVX42Mc+Vujt7TU9pIbx9YzInj17lMlk1NTUpCuuuEJtbW36xCc+oVdeecX00ALj6NGjuu222/TDH/5Qp59+uunhBN74+LjOOuss08PwnZMnT+qFF17Q6tWrJ481NTVp9erVevbZZw2OLJjGx8clifeqgzZt2qRPfvKTRe/hsPB1IPLaa69Jkr70pS/pC1/4gn72s59p4cKF6uzs1B/+8AfDo/O/QqGgjRs36vbbb9fKlStNDyfwDhw4oAceeEB/93d/Z3oovvP73/9e+XxeixcvLjq+ePFivfHGG4ZGFUwTExPq6+vTNddco0svvdT0cALh0Ucf1Z49e9Tf3296KEZ4MhC56667FIlEyt7s9XVJ+vznP68bbrhBK1as0COPPKJIJKIf//jHhv8K76r29X3ggQd04sQJ3X333aaH7CvVvr5TZTIZrVmzRmvXrtVtt91maORAZZs2bdIrr7yiRx991PRQAmF0dFS9vb360Y9+pNNOO830cIzwZIv3Y8eO6a233ip7zvnnn69f/vKX+vjHP66RkRF99KMfnXxu1apVWr16tb761a+6PVRfqvb1XbdunX76058qEolMHs/n84pGo7rxxhu1detWt4fqS9W+vvPnz5ckHTlyRJ2dnbr66qv1gx/8QE1Nnvx+4GknT57U6aefrqGhoaKquZtvvllvv/22nnjiCXODC5A77rhDTzzxhP7zP/9Ty5cvNz2cQPjXf/1X/fVf/7Wi0ejksXw+r0gkoqamJuVyuaLngmie6QGU0traqtbW1ornrVixQrFYTPv3758MRN59910dOnRIy5Ytc3uYvlXt6/vtb39bX/nKVyYfHzlyRNddd50ee+wxrVq1ys0h+lq1r69kzYR0dXVNzuYRhMzN/PnztWLFCu3atWsyEJmYmNCuXbt0xx13mB1cABQKBf393/+9Hn/8cQ0PDxOEOOjaa6/Vyy+/XHTslltu0UUXXaR//Md/DHwQInk0EKlWPB7X7bffrs2bNyuZTGrZsmXasmWLJGnt2rWGR+d/S5cuLXr8/ve/X5J0wQUXKJFImBhSoGQyGXV2dmrZsmW67777dOzYscnnzj33XIMj86c777xTN998s1auXKmrrrpKAwMDeuedd3TLLbeYHprvbdq0Sdu2bdMTTzyh5ubmybyblpYWLViwwPDo/K25uXlGrs0ZZ5yhRYsWhSYHx9eBiCRt2bJF8+bN00033aQ//vGPWrVqlZ5++mktXLjQ9NCAsp566ikdOHBABw4cmBHYeXDF1PM+9alP6dixY/riF7+oN954Q5dffrl27tw5I4EVtXvwwQclSZ2dnUXHH3nkEW3cuLHxA0KgeDJHBAAAhAML0gAAwBgCEQAAYAyBCAAAMIZABAAAGEMgAgAAjCEQAQAAxhCIAAAAYwhEAACAMQQiAADAGAIRAABgDIEIAAAwhkAEAAAY8/8BOT4KnWWMikIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "X = np.r_[\n",
    "    np.random.randn(20, 2) + [-3.0, -3.0],\n",
    "    np.random.randn(20, 2) + [3.0, 3.0],\n",
    "    np.random.randn(20, 2) + [-3.0, 3.0],\n",
    "    np.random.randn(60, 2) + [3.0, -3.0],\n",
    "]\n",
    "y = 20 * [0] + 20 * [1] + 20 * [2] + 20 * [3]\n",
    "\n",
    "for i in range(80):\n",
    "    if y[i] == 0:\n",
    "        plt.scatter(X[i][0], X[i][1], c='b')\n",
    "    elif y[i] == 1:\n",
    "        plt.scatter(X[i][0], X[i][1], c='r')\n",
    "    elif y[i] == 2:\n",
    "        plt.scatter(X[i][0], X[i][1], c='y')\n",
    "    else:\n",
    "        plt.scatter(X[i][0], X[i][1], c='g')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "fcm = FuzzyCMeans(3, 2)\n",
    "fcm.initialize(X)\n",
    "fcm.cluster(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "运行上面的代码我们可以发现这个算法极其的慢，分析一下每次迭代的时间复杂度高达$\\Theta(n^2c^2d)$，虽然我们只有80个样本，而且只迭代了十次，但是上面的算法还是跑了好几秒。接下来我们看一下损失函数的变化曲线。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAGwCAYAAACzXI8XAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA21klEQVR4nO3de3xT9f3H8Xd6v5AEitBS20KRclFuAg5rZTLsRFSU0Q0vzNt4zD0cILR4YxtuOBV0QxiCMJEfyOMnQ1Hxtp8oqwxEgSEXxTlLxWqr0HoBmrbQUpr8/gjNCDebNslJcl7PxyOPnJyTnHyOdea97/mc77G4XC6XAAAAwlCU0QUAAAC0FkEGAACELYIMAAAIWwQZAAAQtggyAAAgbBFkAABA2CLIAACAsBVjdAGB5nQ6tW/fPlmtVlksFqPLAQAALeByuVRTU6P09HRFRZ153CXig8y+ffuUmZlpdBkAAKAVKioqlJGRccbtER9krFarJPc/CJvNZnA1AACgJRwOhzIzMz2/42cS8UGm+XSSzWYjyAAAEGa+ry2EZl8AABC2CDIAACBsEWQAAEDYIsgAAICwRZABAABhiyADAADCFkEGAACELYIMAAAIWwQZAAAQtiJ+Zt+AcDZJ37wjHdkvJXaROg2ToqKNrgoAANMhyPiq4iVp+xTp8Jf/XZeUIQ3+i5Q51ri6AAAwIU4t+aLiJemdn3qHGEk6/JV7fcVLxtQFAIBJEWRaytnkHomR6zQbj6/bPtX9PgAAEBQEmZb65p1TR2K8uKTDFe73AQCAoCDItNSR/f59HwAAaDOCTEsldvHv+wAAQJsRZFqq0zD31UmynOENFikp0/0+AAAQFASZloqKdl9iLenUMHP89eB5zCcDAEAQEWR8kTlWGvaClHSu9/qEzu71zCMDAEBQEWR8lTlWuvZz6fL1kjXHvY7J8AAAMARBpjWioqXU4dI5l7hf13xqaDkAAJgVQaYtbD3dzzV7jK0DAACTIsi0hfV4kHEQZAAAMAJBpi1svdzPNSWS63S3LgAAAIFEkGmLdj3cz0cPSg3fGVsLAAAmRJBpi5hEKSnLvUyfDAAAQUeQaSsafgEAMAxBpq1o+AUAwDCGB5mvvvpKP//5z9WxY0clJiaqX79+ev/99z3bXS6XHnjgAXXp0kWJiYnKz89XaWmpgRWf5MSGXwAAEFSGBpmDBw8qLy9PsbGxeuONN/Txxx9rzpw56tChg+c9jz32mObPn6/Fixdr69atSk5O1siRI1VfX29g5SdgRAYAAMPEGPnljz76qDIzM7Vs2TLPuuzsbM+yy+XSvHnz9Lvf/U7XXXedJGnFihVKTU3Vyy+/rBtuuCHoNZ/C0yNTKrmcksXwQS4AAEzD0F/dV199VUOGDNHPfvYzde7cWRdeeKGWLFni2V5WVqbKykrl5+d71tntdg0dOlSbN28+7T4bGhrkcDi8HgGV1FWKipWcDdLhisB+FwAA8GJokPnss8+0aNEi5eTk6M0339Sdd96pu+66S88884wkqbKyUpKUmprq9bnU1FTPtpPNmjVLdrvd88jMzAzsQURF/3c+GU4vAQAQVIYGGafTqUGDBumRRx7RhRdeqDvuuEO//OUvtXjx4lbvc/r06aqurvY8KiqCMErS3PDroOEXAIBgMjTIdOnSReeff77Xuj59+qi8vFySlJaWJkmqqqryek9VVZVn28ni4+Nls9m8HgFnZS4ZAACMYGiQycvLU0mJ9yjGnj171LVrV0nuxt+0tDQVFxd7tjscDm3dulW5ublBrfWsmBQPAABDGHrVUmFhoS655BI98sgjGjdunP71r3/pqaee0lNPPSVJslgsmjp1qh566CHl5OQoOztbM2bMUHp6usaMGWNk6d64BBsAAEMYGmQuuugirVmzRtOnT9eDDz6o7OxszZs3T+PHj/e8595771VdXZ3uuOMOHTp0SJdeeqnWrl2rhIQEAys/SXOQqftcamqQouMNLQcAALOwuFwul9FFBJLD4ZDdbld1dXXg+mVcLumF9lKjQ7rqI6n9BYH5HgAATKKlv9/M3uYPFotkbb5VAaeXAAAIFoKMv9DwCwBA0BFk/IWGXwAAgo4g4y/MJQMAQNARZPyFU0sAAAQdQcZfmkdk6r+Wjh4ytBQAAMyCIOMvse2kxHT3Mn0yAAAEBUHGn+iTAQAgqAgy/kSfDAAAQUWQ8ScuwQYAIKgIMv5ka57dt+Ts7wMAAH5BkPGnE0dkIvsWVgAAhASCjD+1y5Ys0VLTYenIPqOrAQAg4hFk/CkqVmrX3b1Mwy8AAAFHkPE3Gn4BAAgagoy/NTf8Omj4BQAg0Agy/sakeAAABA1Bxt+YFA8AgKAhyPhb84hM7WeSs9HYWgAAiHAEGX9LTJdikiVXk1RbZnQ1AABENIKMv1ksJ1y5RMMvAACBRJAJBBp+AQAICoJMINDwCwBAUBBkAoFJ8QAACAqCTCBwagkAgKAgyARC86mlI/ukxhpjawEAIIIRZAIhrr2U0Nm9XFNqaCkAAEQygkyg0CcDAEDAEWQChT4ZAAACjiATKFyCDQBAwBFkAoXZfQEACDiCTKDYermfa/ZILpextQAAEKEIMoHS7jxJFqnRIdV/bXQ1AABEJIJMoETHS8nd3Mv0yQAAEBAEmUCi4RcAgIAiyAQSc8kAABBQBJlA8jT8cuUSAACBQJAJJEZkAAAIKIJMIDX3yNR+KjmbjK0FAIAIRJAJpKRMKSpecjZKh78wuhoAACIOQSaQLFGSNce9zOklAAD8jiATaM0Nv9yqAAAAvzM0yPzhD3+QxWLxevTu3duzvb6+XhMnTlTHjh3Vrl07FRQUqKqqysCKW4G7YAMAEDCGj8hccMEF2r9/v+exadMmz7bCwkK99tprWr16tTZs2KB9+/Zp7NixBlbbCkyKBwBAwMQYXkBMjNLS0k5ZX11draVLl2rlypUaMWKEJGnZsmXq06ePtmzZoosvvjjYpbYOl2ADABAwho/IlJaWKj09Xd27d9f48eNVXl4uSdq+fbsaGxuVn5/veW/v3r2VlZWlzZs3n3F/DQ0NcjgcXg9DNQeZw+XSsSPG1gIAQIQxNMgMHTpUy5cv19q1a7Vo0SKVlZVp2LBhqqmpUWVlpeLi4tS+fXuvz6SmpqqysvKM+5w1a5bsdrvnkZmZGeCj+B7xHaW4Du7lmlJjawEAIMIYempp1KhRnuX+/ftr6NCh6tq1q55//nklJia2ap/Tp09XUVGR57XD4TA2zFgskrWX9N0Wd59Mh/7G1QIAQIQx/NTSidq3b6+ePXvq008/VVpamo4ePapDhw55vaeqquq0PTXN4uPjZbPZvB6Go+EXAICACKkgU1tbq71796pLly4aPHiwYmNjVVxc7NleUlKi8vJy5ebmGlhlK9DwCwBAQBh6aunuu+/W6NGj1bVrV+3bt0+///3vFR0drRtvvFF2u10TJkxQUVGRUlJSZLPZNHnyZOXm5obPFUvNGJEBACAgDA0yX375pW688UZ999136tSpky699FJt2bJFnTp1kiTNnTtXUVFRKigoUENDg0aOHKknn3zSyJJbxzMiw+y+AAD4k8XlcrmMLiKQHA6H7Ha7qqurjeuXOXZYej7ZvVzwrftKJgAAcEYt/f0OqR6ZiBWT5L4TtkSfDAAAfkSQCRbuuQQAgN8RZIKFhl8AAPyOIBMsNPwCAOB3BJlgsfVyPzMiAwCA3xBkgsXTI1MquZzG1gIAQIQgyARLclcpKlZqqpcOf2l0NQAARASCTLBExUjtznMvc3oJAAC/IMgEEw2/AAD4FUEmmJobfpkUDwAAvyDIBBOT4gEA4FcEmWBiUjwAAPyKIBNMzSMydZ9LTQ2GlgIAQCQgyARTQqoUY3XPI1P7mdHVAAAQ9ggywWSxnNDwy5VLAAC0FUEm2Gj4BQDAbwgywUbDLwAAfkOQCTbPpHgEGQAA2oogE2yMyAAA4DcEmWCz5rif66uko9XG1gIAQJgjyARbrE1K7OJeZlQGAIA2IcgYgT4ZAAD8giBjBC7BBgDALwgyRqDhFwAAvyDIGMFzaonZfQEAaAuCjBGab1NQs0dyuYytBQCAMEaQMUJytmSJlo7VSUf2G10NAABhiyBjhOg4d5iR6JMBAKANCDJGoeEXAIA2I8gYhYZfAADajCBjlOaGXybFAwCg1QgyRmFSPAAA2owgY5TmHpnazyRno7G1AAAQpggyRklMl6KTJNcxqfZzo6sBACAsEWSMYomSrDnu5RoafgEAaA2CjJFo+AUAoE1ifP1AfX29nnjiCa1fv15ff/21nE6n1/YdO3b4rbiIR8MvAABt4nOQmTBhgt566y399Kc/1Q9+8ANZLJZA1GUOTIoHAECb+BxkXn/9df3f//2f8vLyAlGPuXgmxSPIAADQGj73yJx77rmyWq2BqMV8mpt9j3wlNdYaWwsAAGHI5yAzZ84c3Xffffriiy8CUY+5xKdI8ee4l2tKja0FAIAw5POppSFDhqi+vl7du3dXUlKSYmNjvbYfOHDAb8WZgq2X9M237j6ZlAuNrgYAgLDic5C58cYb9dVXX+mRRx5Ramqq35p9Z8+erenTp2vKlCmaN2+eJPcVUtOmTdOqVavU0NCgkSNH6sknn1RqaqpfvjMkWHtK37xLnwwAAK3gc5B57733tHnzZg0YMMBvRWzbtk1//etf1b9/f6/1hYWF+vvf/67Vq1fLbrdr0qRJGjt2rN59912/fbfhuAQbAIBW87lHpnfv3jpy5IjfCqitrdX48eO1ZMkSdejQwbO+urpaS5cu1eOPP64RI0Zo8ODBWrZsmd577z1t2bLFb99vuOZLsB3M7gsAgK98DjKzZ8/WtGnT9M9//lPfffedHA6H18NXEydO1NVXX638/Hyv9du3b1djY6PX+t69eysrK0ubN28+4/4aGhraXFNQnTgi43IZWwsAAGHG51NLV155pSTp8ssv91rvcrlksVjU1NTU4n2tWrVKO3bs0LZt207ZVllZqbi4OLVv395rfWpqqiorK8+4z1mzZmnmzJktrsFw1h6SLFJjtdTwjZTQ2eiKAAAIGz4HmfXr1/vliysqKjRlyhStW7dOCQkJftmnJE2fPl1FRUWe1w6HQ5mZmX7bv99FJ0jJXaW6z90NvwQZAABazOcgc9lll/nli7dv366vv/5agwYN8qxramrSxo0btWDBAr355ps6evSoDh065DUqU1VVpbS0tDPuNz4+XvHx8X6pMWisPd1BpmaP1PlSo6sBACBs+BxkNm7ceNbtP/zhD1u0n8svv1y7d+/2Wnf77berd+/euu+++5SZmanY2FgVFxeroKBAklRSUqLy8nLl5ub6WnZos/WUKt/iyiUAAHzkc5AZPnz4KetOnEumpT0yVqtVffv29VqXnJysjh07etZPmDBBRUVFSklJkc1m0+TJk5Wbm6uLL77Y17JDm5UrlwAAaA2fg8zBgwe9Xjc2Nmrnzp2aMWOGHn74Yb8VJklz585VVFSUCgoKvCbEizi2Xu5nRmQAAPCJxeXyzzW/GzZsUFFRkbZv3+6P3fmNw+GQ3W5XdXW1bDab0eWcXu3n0qvZUlScNO6wFBVtdEUAABiqpb/fPs8jcyapqakqKeHUSKskZUpR8ZLzqHS43OhqAAAIGz6fWvrwww+9XrtcLu3fv1+zZ8/WwIED/VWXuURFu+eTqf63+xLsdtlGVwQAQFjwOcgMHDhQFotFJ5+Ruvjii/U///M/fivMdKw93UGmpkTSSKOrAQAgLPgcZMrKyrxeR0VFqVOnTn6d1M6Umht+uQs2AAAt5nOQ6dq16ynrDh06RJBpK+6CDQCAz3xu9n300Uf13HPPeV6PGzdOKSkpOvfcc/XBBx/4tThTsRFkAADwlc9BZvHixZ57F61bt07r1q3T2rVrNWrUKN1zzz1+L9A0mkdk6sqlY0eMrQUAgDDh86mlyspKT5B5/fXXNW7cOF1xxRXq1q2bhg4d6vcCTSP+HCm2vdR4SKr9VGrfz+iKAAAIeT6PyHTo0EEVFRWSpLVr1yo/P1+S+zLslt6eAKdhsdDwCwCAj3wekRk7dqxuuukm5eTk6LvvvtOoUaMkSTt37lSPHj38XqCpWHtK322lTwYAgBbyOcjMnTtX3bp1U0VFhR577DG1a9dOkrR//379+te/9nuBpkLDLwAAPvE5yMTGxuruu+8+ZX1hYaFfCjI1z12wCTIAALSE3+61BD/wjMhwzyoAAFqCIBNKrDnu54bv3A8AAHBWBJlQEpMsJWW4l2tKja0FAIAwQJAJNfTJAADQYq0KMocOHdLTTz+t6dOn68CBA5KkHTt26KuvvvJrcabEPZcAAGgxn69a+vDDD5Wfny+73a7PP/9cv/zlL5WSkqKXXnpJ5eXlWrFiRSDqNI/mhl8HDb8AAHwfn0dkioqKdNttt6m0tNTrjtdXXXWVNm7c6NfiTIkRGQAAWsznILNt2zb96le/OmX9ueeeq8rKSr8UZWrNtymoKZVcTmNrAQAgxPkcZOLj4+VwOE5Zv2fPHnXq1MkvRZlacjfJEiM1HZEO03MEAMDZ+Bxkrr32Wj344INqbGyUJFksFpWXl+u+++5TQUGB3ws0nagYyXqee5nTSwAAnJXPQWbOnDmqra1V586ddeTIEV122WXq0aOHrFarHn744UDUaD5WGn4BAGgJn69astvtWrdund5991198MEHqq2t1aBBg5Sfnx+I+syJhl8AAFrEpyDT2NioxMRE7dq1S3l5ecrLywtUXebW3PDLpHgAAJyVT6eWYmNjlZWVpaampkDVA4kRGQAAWsjnHpnf/va3+s1vfuOZ0RcB0DwpXl2Z1HTU2FoAAAhhPvfILFiwQJ9++qnS09PVtWtXJScne23fsWOH34ozrYQ0KaaddKxWqv1Msvc2uiIAAEKSz0FmzJgxASgDXiwW9+mlgzukmhKCDAAAZ+BzkPn9738fiDpwMlsvd5Ch4RcAgDNq1d2vEQQ0/AIA8L18HpFpamrS3Llz9fzzz6u8vFxHj3o3o9IE7Cc2ggwAAN/H5xGZmTNn6vHHH9f111+v6upqFRUVaezYsYqKitIf/vCHAJRoUp7ZfQkyAACcic9B5tlnn9WSJUs0bdo0xcTE6MYbb9TTTz+tBx54QFu2bAlEjeZkzXE/11dKjafepBMAALQiyFRWVqpfv36SpHbt2qm6ulqSdM011+jvf/+7f6szszi7lJDqXmZUBgCA0/I5yGRkZGj//v2SpPPOO09vvfWWJGnbtm2Kj4/3b3Vm13yrAvpkAAA4LZ+DzE9+8hMVFxdLkiZPnqwZM2YoJydHt9xyi37xi1/4vUBTo08GAICz8vmqpdmzZ3uWr7/+emVlZWnz5s3KycnR6NGj/Vqc6XEJNgAAZ+VzkDlZbm6ucnNz/VELTtZ8CbajxNg6AAAIUT4HmRUrVpx1+y233NLqYnCSE0dkXC73rQsAAICHxeVyuXz5QIcOHbxeNzY26vDhw4qLi1NSUlLITYjncDhkt9tVXV0tm81mdDm+aToqPZ8ouZzST/ZJiV2MrggAgKBo6e+3z82+Bw8e9HrU1taqpKREl156qf72t7+1qWicJDpOSs52L9PwCwDAKfxyr6WcnBzNnj1bU6ZM8elzixYtUv/+/WWz2WSz2ZSbm6s33njDs72+vl4TJ05Ux44d1a5dOxUUFKiqqsofJYcPGn4BADgjv900MiYmRvv27fPpMxkZGZo9e7a2b9+u999/XyNGjNB1112nf//735KkwsJCvfbaa1q9erU2bNigffv2aezYsf4qOTzQ8AsAwBn53Oz76quver12uVzav3+/FixYoLy8PJ/2dfLl2g8//LAWLVqkLVu2KCMjQ0uXLtXKlSs1YsQISdKyZcvUp08fbdmyRRdffLGvpYcnRmQAADgjn4PMmDFjvF5bLBZ16tRJI0aM0Jw5c1pdSFNTk1avXq26ujrl5uZq+/btamxsVH5+vuc9vXv39sxbc6Yg09DQoIaGBs9rhyPM71PE7L4AAJyRz0HG6XT6tYDdu3crNzdX9fX1ateundasWaPzzz9fu3btUlxcnNq3b+/1/tTUVFVWVp5xf7NmzdLMmTP9WqOhPCMyeyXnMSmqzVP/AAAQMfzWI9NavXr10q5du7R161bdeeeduvXWW/Xxxx+3en/Tp09XdXW151FRUeHHag2QdK4UnSi5jkl1nxtdDQAAIcXn/3tfVFTU4vc+/vjj3/ueuLg49ejRQ5I0ePBgbdu2TX/5y190/fXX6+jRozp06JDXqExVVZXS0tLOuL/4+PjIunmlJUqy5kiHPnQ3/Fp7GF0RAAAhw+cgs3PnTu3cuVONjY3q1cvdv7Fnzx5FR0dr0KBBnvdZWjkLrdPpVENDgwYPHqzY2FgVFxeroKBAklRSUqLy8nLz3RLB2tMdZGr2SLra6GoAAAgZPgeZ0aNHy2q16plnnvHM8nvw4EHdfvvtGjZsmKZNm9bifU2fPl2jRo1SVlaWampqtHLlSv3zn//Um2++KbvdrgkTJqioqEgpKSmy2WyaPHmycnNzzXPFUrPmhl8mxQMAwIvPQWbOnDl66623vG5V0KFDBz300EO64oorfAoyX3/9tW655Rbt379fdrtd/fv315tvvqkf//jHkqS5c+cqKipKBQUFamho0MiRI/Xkk0/6WnL44xJsAABOy+cg43A49M0335yy/ptvvlFNTY1P+1q6dOlZtyckJGjhwoVauHChT/uNODaCDAAAp+PzVUs/+clPdPvtt+ull17Sl19+qS+//FIvvviiJkyYYL5Zd4OleUTm8JfSsTpjawEAIIT4PCKzePFi3X333brpppvU2Njo3klMjCZMmKA//elPfi8QkuJTpPiOUsN3Uk2p1GGg0RUBABASLC6Xy9WaD9bV1Wnv3r2SpPPOO0/Jycl+LcxfWnob8JD3Vp707XtS3nNS13FGVwMAQEC19Pe71RPiJScnq3///urWrZvWrVunTz75pLW7QkvQJwMAwCl8DjLjxo3TggULJElHjhzRkCFDNG7cOPXr108vvvii3wvEcc19MlyCDQCAh89BZuPGjRo2bJgkac2aNXK5XDp06JDmz5+vhx56yO8F4jguwQYA4BQ+B5nq6mqlpKRIktauXauCggIlJSXp6quvVmlpqd8LxHHNp5YcJVLr2poAAIg4PgeZzMxMbd68WXV1dVq7dq2uuOIKSe7ZfRMSEvxeII5r10OSRWo8JDV8a3Q1AACEBJ+DzNSpUzV+/HhlZGQoPT1dw4cPl+Q+5dSvXz9/14dmMYlScpZ7mdNLAABIasU8Mr/+9a81dOhQlZeX68c//rGiotxZqHv37vTIBJq1p1T3hbvht1Oe0dUAAGA4n4OMJA0ePFiDBw/2Wnf11dyVOeCsPaXKdYzIAABwXKvnkYEBTmz4BQAABJmwwiXYAAB4IciEE1sv93PNp5KzydhaAAAIAQSZcJKUJUXFSc4G6XCF0dUAAGA4n4PM2rVrtWnTJs/rhQsXauDAgbrpppt08OBBvxaHk0RFS9Ye7mVOLwEA4HuQueeee+RwOCRJu3fv1rRp03TVVVeprKxMRUVFfi8QJ7HS8AsAQDOfL78uKyvT+eefL0l68cUXdc011+iRRx7Rjh07dNVVV/m9QJyEhl8AADx8HpGJi4vT4cOHJUn/+Mc/PLcoSElJ8YzUIICaG365CzYAAL6PyFx66aUqKipSXl6e/vWvf+m5556TJO3Zs0cZGRl+LxAnYUQGAAAPn0dkFixYoJiYGL3wwgtatGiRzj33XEnSG2+8oSuvvNLvBeIkzZPi1X0hNdUbWwsAAAazuFwul9FFBJLD4ZDdbld1dbVsNpvR5bSdyyW90EFqrJau2i2172t0RQAA+F1Lf7/bNI9MfX29HA6H1wMBZrFwegkAgON8DjJ1dXWaNGmSOnfurOTkZHXo0MHrgSCg4RcAAEmtCDL33nuv3n77bS1atEjx8fF6+umnNXPmTKWnp2vFihWBqBEnY0QGAABJrbhq6bXXXtOKFSs0fPhw3X777Ro2bJh69Oihrl276tlnn9X48eMDUSdOZCPIAAAgtWJE5sCBA+revbskyWaz6cCBA5Lcl2Vv3LjRv9Xh9JjdFwAASa0IMt27d1dZWZkkqXfv3nr++ecluUdq2rdv79ficAbWHPdzw7dSwwFjawEAwEA+B5nbb79dH3zwgSTp/vvv18KFC5WQkKDCwkLdc889fi8QpxHbTkpMdy/XlBpbCwAABvK5R6awsNCznJ+fr08++UTbt29Xjx491L9/f78Wh7Ow9ZKO7HP3yZwz1OhqAAAwhM9B5mRdu3ZV165d/VELfGHtKVWt5xJsAICptSjIzJ8/X3fccYcSEhI0f/78s773rrvu8kth+B6eS7Bp+AUAmFeLblGQnZ2t999/Xx07dlR2dvaZd2ax6LPPPvNrgW0VcbcoaPbV69KG0VL7AdJVu4yuBgAAv2rp73eLRmSar1I6eRkG8ozIlEoup2Rp090mAAAIS2369XO5XIrwe06GrnbZkiVGajrsbvoFAMCEWhVkli5dqr59+yohIUEJCQnq27evnn76aX/XhrOJipXauScmpOEXAGBWPl+19MADD+jxxx/X5MmTlZubK0navHmzCgsLVV5ergcffNDvReIMrD3dl1/X7JHSRhhdDQAAQedzkFm0aJGWLFmiG2+80bPu2muvVf/+/TV58mSCTDDZekr7xK0KAACm5fOppcbGRg0ZMuSU9YMHD9axY8f8UhRaiLtgAwBMzucgc/PNN2vRokWnrH/qqae483Ww2Xq5n+mRAQCYVItOLRUVFXmWLRaLnn76ab311lu6+OKLJUlbt25VeXm5brnllsBUidNrHpGpK5OajkrRccbWAwBAkLUoyOzcudPr9eDBgyVJe/fulSSdc845Ouecc/Tvf//bz+XhrBK7SDHJ0rE6d5hpHqEBAMAkWhRk1q9fH5AvnzVrll566SV98sknSkxM1CWXXKJHH31UvXr99we5vr5e06ZN06pVq9TQ0KCRI0fqySefVGpqakBqCisWi3tU5uBOd8MvQQYAYDKGTge7YcMGTZw4UVu2bNG6devU2NioK664QnV1dZ73FBYW6rXXXtPq1au1YcMG7du3T2PHjjWw6hBDwy8AwMTafPfrtli7dq3X6+XLl6tz587avn27fvjDH6q6ulpLly7VypUrNWKEe56UZcuWqU+fPtqyZYunR8fUaPgFAJhYSN2gp7q6WpKUkpIiSdq+fbsaGxuVn5/veU/v3r2VlZWlzZs3n3YfDQ0NcjgcXo+IxogMAMDEQibIOJ1OTZ06VXl5eerbt68kqbKyUnFxcWrfvr3Xe1NTU1VZWXna/cyaNUt2u93zyMzMDHTpxrIRZAAA5hUyQWbixIn66KOPtGrVqjbtZ/r06aqurvY8Kioq/FRhiLLmuJ+P7JcaI3z0CQCAkxjaI9Ns0qRJev3117Vx40ZlZGR41qelpeno0aM6dOiQ16hMVVWV0tLSTruv+Ph4xcfHB7rk0BHXXkroLNV/LdWUSimDja4IAICgMXRExuVyadKkSVqzZo3efvttZWdne20fPHiwYmNjVVxc7FlXUlKi8vJyzw0rIclKwy8AwJwMHZGZOHGiVq5cqVdeeUVWq9XT92K325WYmCi73a4JEyaoqKhIKSkpstlsnrtuc8XSCWw9pW/eoU8GAGA6hgaZ5ns2DR8+3Gv9smXLdNttt0mS5s6dq6ioKBUUFHhNiIcTNF+5xIgMAMBkDA0yLpfre9+TkJCghQsXauHChUGoKEx5LsEuMbYOAACCLGSuWkIb2E4YkWlBOAQAIFIQZCJBu/MkS5R0rEaqrzK6GgAAgoYgEwmi46Xkbu5lGn4BACZCkIkUNPwCAEyIIBMpaPgFAJgQQSZS2BiRAQCYD0EmUnAXbACACRFkIoXt+G0KavdKzmPG1gIAQJAQZCJFUoYUnSA5G6W6L4yuBgCAoCDIRApLlGTNcS87aPgFAJgDQSaS0CcDADAZgkwkIcgAAEyGIBNJmht+uQQbAGASBJlIwogMAMBkCDKRpHlSvMMV0rHDxtYCAEAQEGQiSXxHKS7FvVxTamwtAAAEAUEm0nB6CQBgIgSZSEPDLwDARAgykcbGiAwAwDwIMpHGyl2wAQDmQZCJNJ4emRLJ5TK2FgAAAowgE2msPdzPRw9KDd8ZWwsAAAFGkIk0MUlSUqZ7mT4ZAECEI8hEouYrlwgyAIAIR5CJRDT8AgBMgiATiU5s+AUAIIIRZCKRjREZAIA5EGQikWdEplRyOY2tBQCAACLIRKLkblJUrORscN8JGwCACEWQiURR0VK74/PJcHoJABDBCDKRytMnQ8MvACByEWQilZWbRwIAIh9BJlIRZAAAJkCQiVTNs/vSIwMAiGAEmUjVPCJT97nU1GBoKQAABApBJlIldJZibZJcUs2nRlcDAEBAEGQilcVCnwwAIOIRZCIZQQYAEOEIMpGMhl8AQIQjyEQyRmQAABGOIBPJmN0XABDhCDKRzJrjfm74Rjp60NhaAAAIAEODzMaNGzV69Gilp6fLYrHo5Zdf9trucrn0wAMPqEuXLkpMTFR+fr5KS0uNKTYcxVqlxC7uZQf/3AAAkcfQIFNXV6cBAwZo4cKFp93+2GOPaf78+Vq8eLG2bt2q5ORkjRw5UvX19UGuNIxZjzf80icDAIhAMUZ++ahRozRq1KjTbnO5XJo3b55+97vf6brrrpMkrVixQqmpqXr55Zd1ww03BLPU8GXrKX39T4IMACAihWyPTFlZmSorK5Wfn+9ZZ7fbNXToUG3evPmMn2toaJDD4fB6mFrzlUtcgg0AiEAhG2QqKyslSampqV7rU1NTPdtOZ9asWbLb7Z5HZmZmQOsMeZ5LsLlyCQAQeUI2yLTW9OnTVV1d7XlUVFQYXZKxbCeMyLhcxtYCAICfhWyQSUtLkyRVVVV5ra+qqvJsO534+HjZbDavh6klZ0uWaKnpsHRkn9HVAADgVyEbZLKzs5WWlqbi4mLPOofDoa1btyo3N9fAysJMdJzUrrt7mYZfAECEMfSqpdraWn366aee12VlZdq1a5dSUlKUlZWlqVOn6qGHHlJOTo6ys7M1Y8YMpaena8yYMcYVHY6sPaWaUvfppdQfGV0NAAB+Y2iQef/99/WjH/33h7WoqEiSdOutt2r58uW69957VVdXpzvuuEOHDh3SpZdeqrVr1yohIcGoksOTtaekv3OrAgBAxLG4XJHdAepwOGS321VdXW3efpnSxdK2O6X0q6XhrxtdDQAA36ulv98h2yMDP+Iu2ACACEWQMQPb8dsU1H4mORuNrQUAAD8iyJhBYroUnSS5mqTaMqOrAQDAbwgyZmCxnDAxHg2/AIDIQZAxC/pkAAARiCBjFgQZAEAEIsiYRXPDL3fBBgBEEIKMWTAiAwCIQAQZs7DluJ+P7JMaa4ytBQAAPyHImEVcBym+k3u5ptTYWgAA8BOCjJl4LsHm9BIAIDIQZMzEerzhlz4ZAECEIMiYiY2GXwBAZCHImImV2X0BAJGFIGMmJ16C7XIZWwsAAH5AkDET63mSLFKjQ6r/2uhqAABoM4KMmUQnSMld3cv0yQAAIgBBxmxsXLkEAIgcBBmzoeEXABBBCDJmwz2XAAARhCBjNszuCwCIIAQZs2kekan9VHI2GVsLAABtRJAxm+QsKSpecjZKh78wuhoAANqEIGM2lijJmuNe5vQSACDMEWTMyMaVSwCAyECQMSOuXAIARAiCjBkRZAAAEYIgY0bNs/vSIwMACHMEGTNqHpE5XC4dO2JsLQAAtAFBxoziO0pxHdzLNaXG1gIAQBsQZMzIYqFPBgAQEQgyZkWQAQBEAIKMWdHwCwCIAAQZs7IxIgMACH8EGbOyMrsvACD8EWTMytrD/Xz0gNTwnbG1AADQSgQZs4pJlpIy3Mv0yQAAwhRBxsysxxt+6ZMBAIQpgoyZ0fALAAhzBBkzo+EXABDmCDJmxqR4AIAwF2N0ATBQ86klxydS2bNS0rlSp2FSVLSxdbWGs0n65h3pyH4psQvHYTSOI7RwHKGF4/CrsAgyCxcu1J/+9CdVVlZqwIABeuKJJ/SDH/zA6LLC34Gd7mdno7T55+7lpAxp8F+kzLHG1eWripek7VOkw1/+dx3HYRyOI7RwHKGF4/A7i8vlcgX1G3303HPP6ZZbbtHixYs1dOhQzZs3T6tXr1ZJSYk6d+78vZ93OByy2+2qrq6WzWYLQsVhouIl6Z2fSjr5z29xPw17ITz+R8VxhBaOI7RwHKGF4/BJS3+/Qz7IDB06VBdddJEWLFggSXI6ncrMzNTkyZN1//33f+/nCTKn4WySXu3mnaRPFt9Juni5FBUjyeK+Y7YsJy3rLNuObz/dtrZ+vnnZ1ST9Y7hUX3mGg7BIiWnSjzdJlhAetnU1SW9dKtXvP8MbOI6g4jhCC8cRWlpyHEkZ0rVlbT7NFBFB5ujRo0pKStILL7ygMWPGeNbfeuutOnTokF555ZVTPtPQ0KCGhgbPa4fDoczMTILMiar+KRX/yOgqAACR6vL1UurwNu2ipUEmpHtkvv32WzU1NSk1NdVrfWpqqj755JPTfmbWrFmaOXNmMMoLX0fOlKRPktxVim0v9/ChS3K5TlrWqdtOec8J7zvttha+78R1kvu186jk/G9oPSNLbGg30jmbJFfj97+P4wgOjiO0cByhpaXH0dLfGT8I6SDTGtOnT1dRUZHndfOIDE6Q2KVl77t4eZsTdUC1dGRpxFscRzBwHKGF4wgtZjuOlv7O+EFIzyNzzjnnKDo6WlVVVV7rq6qqlJaWdtrPxMfHy2azeT1wkk7Djt9nyXKGN1ikpEz3+0IZxxFaOI7QwnGEFo4jYEI6yMTFxWnw4MEqLi72rHM6nSouLlZubq6BlYW5qGj3JXKSTv2X8fjrwfNCe3hT4jhCDccRWjiO0MJxBK6koH1TKxUVFWnJkiV65pln9J///Ed33nmn6urqdPvttxtdWnjLHOu+RC7pXO/1SRnhcwmgxHGEGo4jtHAcoYXjCIiQvmqp2YIFCzwT4g0cOFDz58/X0KFDW/RZLr/+HiEyM2ObcRyhheMILRxHaOE4WiQiLr/2B4IMAADhp6W/3yF/agkAAOBMCDIAACBsEWQAAEDYIsgAAICwRZABAABhiyADAADCFkEGAACELYIMAAAIWwQZAAAQtmKMLiDQmicudjgcBlcCAABaqvl3+/tuQBDxQaampkaSlJmZaXAlAADAVzU1NbLb7WfcHvH3WnI6ndq3b5+sVqsslpNvOQ7JnXozMzNVUVHB/ahCAH+P0MLfI7Tw9wgtgfx7uFwu1dTUKD09XVFRZ+6EifgRmaioKGVkZBhdRliw2Wz8hyGE8PcILfw9Qgt/j9ASqL/H2UZimtHsCwAAwhZBBgAAhC2CDBQfH6/f//73io+PN7oUiL9HqOHvEVr4e4SWUPh7RHyzLwAAiFyMyAAAgLBFkAEAAGGLIAMAAMIWQQYAAIQtgoxJzZo1SxdddJGsVqs6d+6sMWPGqKSkxOiycNzs2bNlsVg0depUo0sxta+++ko///nP1bFjRyUmJqpfv356//33jS7LlJqamjRjxgxlZ2crMTFR5513nv74xz9+73144B8bN27U6NGjlZ6eLovFopdfftlru8vl0gMPPKAuXbooMTFR+fn5Ki0tDUptBBmT2rBhgyZOnKgtW7Zo3bp1amxs1BVXXKG6ujqjSzO9bdu26a9//av69+9vdCmmdvDgQeXl5Sk2NlZvvPGGPv74Y82ZM0cdOnQwujRTevTRR7Vo0SItWLBA//nPf/Too4/qscce0xNPPGF0aaZQV1enAQMGaOHChafd/thjj2n+/PlavHixtm7dquTkZI0cOVL19fUBr43LryFJ+uabb9S5c2dt2LBBP/zhD40ux7Rqa2s1aNAgPfnkk3rooYc0cOBAzZs3z+iyTOn+++/Xu+++q3feecfoUiDpmmuuUWpqqpYuXepZV1BQoMTERP3v//6vgZWZj8Vi0Zo1azRmzBhJ7tGY9PR0TZs2TXfffbckqbq6WqmpqVq+fLluuOGGgNbDiAwkuf+lk6SUlBSDKzG3iRMn6uqrr1Z+fr7RpZjeq6++qiFDhuhnP/uZOnfurAsvvFBLliwxuizTuuSSS1RcXKw9e/ZIkj744ANt2rRJo0aNMrgylJWVqbKy0uu/W3a7XUOHDtXmzZsD/v0Rf9NIfD+n06mpU6cqLy9Pffv2Nboc01q1apV27Nihbdu2GV0KJH322WdatGiRioqK9Jvf/Ebbtm3TXXfdpbi4ON16661Gl2c6999/vxwOh3r37q3o6Gg1NTXp4Ycf1vjx440uzfQqKyslSampqV7rU1NTPdsCiSADTZw4UR999JE2bdpkdCmmVVFRoSlTpmjdunVKSEgwuhzIHfCHDBmiRx55RJJ04YUX6qOPPtLixYsJMgZ4/vnn9eyzz2rlypW64IILtGvXLk2dOlXp6en8PUyOU0smN2nSJL3++utav369MjIyjC7HtLZv366vv/5agwYNUkxMjGJiYrRhwwbNnz9fMTExampqMrpE0+nSpYvOP/98r3V9+vRReXm5QRWZ2z333KP7779fN9xwg/r166ebb75ZhYWFmjVrltGlmV5aWpokqaqqymt9VVWVZ1sgEWRMyuVyadKkSVqzZo3efvttZWdnG12SqV1++eXavXu3du3a5XkMGTJE48eP165duxQdHW10iaaTl5d3ypQEe/bsUdeuXQ2qyNwOHz6sqCjvn6zo6Gg5nU6DKkKz7OxspaWlqbi42LPO4XBo69atys3NDfj3c2rJpCZOnKiVK1fqlVdekdVq9ZzHtNvtSkxMNLg687Faraf0JyUnJ6tjx470LRmksLBQl1xyiR555BGNGzdO//rXv/TUU0/pqaeeMro0Uxo9erQefvhhZWVl6YILLtDOnTv1+OOP6xe/+IXRpZlCbW2tPv30U8/rsrIy7dq1SykpKcrKytLUqVP10EMPKScnR9nZ2ZoxY4bS09M9VzYFlAumJOm0j2XLlhldGo677LLLXFOmTDG6DFN77bXXXH379nXFx8e7evfu7XrqqaeMLsm0HA6Ha8qUKa6srCxXQkKCq3v37q7f/va3roaGBqNLM4X169ef9jfj1ltvdblcLpfT6XTNmDHDlZqa6oqPj3ddfvnlrpKSkqDUxjwyAAAgbNEjAwAAwhZBBgAAhC2CDAAACFsEGQAAELYIMgAAIGwRZAAAQNgiyAAAgLBFkAEAAGGLIAPAr4YPH66pU6caXYYXi8Wil19+2egyAAQAM/sC8KsDBw4oNjZWVqtV3bp109SpU4MWbP7whz/o5Zdf1q5du7zWV1ZWqkOHDoqPjw9KHQCCh5tGAvCrlJQUv+/z6NGjiouLa/Xn09LS/FgNgFDCqSUAftV8amn48OH64osvVFhYKIvFIovF4nnPpk2bNGzYMCUmJiozM1N33XWX6urqPNu7deumP/7xj7rllltks9l0xx13SJLuu+8+9ezZU0lJSerevbtmzJihxsZGSdLy5cs1c+ZMffDBB57vW758uaRTTy3t3r1bI0aMUGJiojp27Kg77rhDtbW1nu233XabxowZoz//+c/q0qWLOnbsqIkTJ3q+C0DoIMgACIiXXnpJGRkZevDBB7V//37t379fkrR3715deeWVKigo0IcffqjnnntOmzZt0qRJk7w+/+c//1kDBgzQzp07NWPGDEmS1WrV8uXL9fHHH+svf/mLlixZorlz50qSrr/+ek2bNk0XXHCB5/uuv/76U+qqq6vTyJEj1aFDB23btk2rV6/WP/7xj1O+f/369dq7d6/Wr1+vZ555RsuXL/cEIwChg1NLAAIiJSVF0dHRslqtXqd2Zs2apfHjx3v6ZnJycjR//nxddtllWrRokRISEiRJI0aM0LRp07z2+bvf/c6z3K1bN919991atWqV7r33XiUmJqpdu3aKiYk566mklStXqr6+XitWrFBycrIkacGCBRo9erQeffRRpaamSpI6dOigBQsWKDo6Wr1799bVV1+t4uJi/fKXv/TLPx8A/kGQARBUH3zwgT788EM9++yznnUul0tOp1NlZWXq06ePJGnIkCGnfPa5557T/PnztXfvXtXW1urYsWOy2Ww+ff9//vMfDRgwwBNiJCkvL09Op1MlJSWeIHPBBRcoOjra854uXbpo9+7dPn0XgMAjyAAIqtraWv3qV7/SXXfddcq2rKwsz/KJQUOSNm/erPHjx2vmzJkaOXKk7Ha7Vq1apTlz5gSkztjYWK/XFotFTqczIN8FoPUIMgACJi4uTk1NTV7rBg0apI8//lg9evTwaV/vvfeeunbtqt/+9reedV988cX3ft/J+vTpo+XLl6uurs4Tlt59911FRUWpV69ePtUEwHg0+wIImG7dumnjxo366quv9O2330pyX3n03nvvadKkSdq1a5dKS0v1yiuvnNJse7KcnByVl5dr1apV2rt3r+bPn681a9ac8n1lZWXatWuXvv32WzU0NJyyn/HjxyshIUG33nqrPvroI61fv16TJ0/WzTff7DmtBCB8EGQABMyDDz6ozz//XOedd546deokSerfv782bNigPXv2aNiwYbrwwgv1wAMPKD09/az7uvbaa1VYWKhJkyZp4MCBeu+99zxXMzUrKCjQlVdeqR/96Efq1KmT/va3v52yn6SkJL355ps6cOCALrroIv30pz/V5ZdfrgULFvjvwAEEDTP7AgCAsMWIDAAACFsEGQAAELYIMgAAIGwRZAAAQNgiyAAAgLBFkAEAAGGLIAMAAMIWQQYAAIQtggwAAAhbBBkAABC2CDIAACBs/T+ldHOx907ojQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "costs = fcm.costs\n",
    "iters = np.arange(1, len(costs) + 1)\n",
    "plt.plot(iters, costs, color='orange', marker='o')\n",
    "plt.xlabel(\"iteration\")\n",
    "plt.ylabel(\"bias square sum\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "很明显其实两次迭代就收敛了，说明算法收敛还是很快的。接下来看一下聚类的中心验证一下书上说的性质。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGdCAYAAAAvwBgXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA4bUlEQVR4nO3df3RU9Z3/8ddkkBA0EwEjBmYgqK3i2qoLFvq12SaVU+nx7Bc7Rlpw/QrHdasHt0ndc7pie0p7tm1qsTWpp0et7SLfbonSOK5rXdl6lNBsV05dqRyliGUlX8IQMEXMUIsDTub7x+0NmWQyP5J75/6Y5+OcnDQ3N5kP4+zOK5/P+/P+BNLpdFoAAAAOqHB6AAAAoHwRRAAAgGMIIgAAwDEEEQAA4BiCCAAAcAxBBAAAOIYgAgAAHEMQAQAAjpni9AByGRoa0uHDh1VdXa1AIOD0cAAAQAHS6bROnDihOXPmqKIi95yHq4PI4cOHFYlEnB4GAACYgL6+PoXD4Zz3uDqIVFdXSzL+IaFQyOHRAACAQiQSCUUikeH38VxcHUTM5ZhQKEQQAQDAYwopq6BYFQAAOIYgAgAAHEMQAQAAjiGIAAAAxxBEAACAYwgiAADAMQQRAADgGIIIAABwjKsbmgEAkFMqJfX0SP39Ul2d1NAgBYNOjwpFIIgAALwpFpNaWqRDh85cC4eljg4pGnVuXCgKSzMAAO+JxaTm5swQIknxuHE9FnNmXCgaQQQA4C2plDETkk6P/Z55rbXVuA+uRxABAHhLT8/YmZCR0mmpr8+4D65HEAEAeEt/v7X3wVEEEQCAt9TVWXsfHEUQAQB4S0ODsTsmEMj+/UBAikSM++B6BBEAgLcEg8YWXWlsGDG/bm+nn4hHEEQAAN4TjUpdXdLcuZnXw2HjOn1EPIOGZgAAb4pGpRUr6KzqcQQRAIB3BYNSY6PTo8Ak2L40E4/H9Td/8zeaNWuWqqqq9JGPfET//d//bffDAgAAD7B1RuT48eO65ppr1NTUpOeee061tbX6/e9/rxkzZtj5sAAAwCNsDSL33XefIpGINm3aNHxtwYIFdj4kAADwEFuXZv7t3/5Nixcv1k033aTzzz9fV111lR599NFx708mk0okEhkfAADAv2wNIm+99ZYeeughfehDH9J//Md/6M4779QXv/hFbd68Oev9bW1tqqmpGf6IRCJ2Dg8AADgskE5nO77QGlOnTtXixYv1X//1X8PXvvjFL+rll1/WSy+9NOb+ZDKpZDI5/HUikVAkEtHg4KBCoZBdwwQAABZKJBKqqakp6P3b1hmRuro6XXbZZRnXFi5cqIMHD2a9v7KyUqFQKOMDAAD4l61B5JprrtG+ffsyrr355puaP3++nQ8LAAA8wtYg8qUvfUk7d+7Ut7/9be3fv19btmzRj370I61bt87OhwUAAB5haxC5+uqr9dRTT6mzs1OXX365/umf/knt7e26+eab7XxYAADgEbYWq05WMcUuAADAHVxTrAoAAJALQQQAADiGIAIAABxDEAEAAI4hiAAAAMcQRAAAgGMIIgAAwDEEEQAA4BiCCAAAcAxBBAAAOIYgAgAAHEMQAQAAjiGIAAAAxxBEAACAYwgiAADAMQQRAADgmClODwAAANdLpaSeHqm/X6qrkxoapGDQ6VH5AkEEAIBcYjGppUU6dOjMtXBY6uiQolHnxuUTLM0AADCeWExqbs4MIZIUjxvXYzFnxuUjBBEAALJJpYyZkHR67PfMa62txn2YMIIIAADZ9PSMnQkZKZ2W+vqM+zBhBBEAALLp77f2PmRFEAEAIJu6OmvvQ1YEEQAAsmloMHbHBALZvx8ISJGIcR8mjCACAEA2waCxRVcaG0bMr9vb6ScySQQRAADGE41KXV3S3LmZ18Nh4zp9RCaNhmYAAOQSjUorVtBZ1SYEEQAA8gkGpcZGp0fhSyzNAAAAxzAjAvhcOp3Su+/26NSpfk2dWqdzz21QIMCUMgB3IIgAPjYwENP+/S1KJs90h6ysDOviiztUW0uRHQDnsTQD+NTAQEx79jRnhBBJSibj2rOnWQMDHNYFwHkEEcCH0umU9u9vkZTlsK4/X9u/v1XpNId1AXAWQQTwoXff7RkzE5IprWSyT+++y2FdQEmlUlJ3t9TZaXzm5F5qRAA/OnWqsEO4Cr0PgAViMamlJfNE33DY6N5axo3RmBEBfGjq1MIO4Sr0PgCTFItJzc2ZIUSS4nHjeqx8a7YIIhZIp1M6frxbR4926vjxbtbd4bhzz21QZWVY0jiHdSmgysqIzj2Xw7oA26VSxkxIOkvNlnmttbVsl2lYmpkktkfCjQKBoC6+uEN79jTLCCMj/x+gEU4uvridfiJAKfT0jJ0JGSmdlvr6jPvKsHsrMyKTwPZIuFltbVR/8RddqqzMPKyrsjKsv/iLLoIyUCr9BdZiFXqfzzAjMkH5t0cGtH9/q847bwV/dcIxtbVRnXfeCjqrAk6qK7AWq9D7fIYgMkHFbI+cMaOxRKMCxgoEgrwGASc1NBi7Y+Lx7HUigYDx/YbyrNliaWaC2B4JAChIMGhs0ZWM0DGS+XV7u3FfGSKITFC5bo9khxAATEA0KnV1SXMza7YUDhvXy7iPCEszE2Ruj0wm48peJxJQZWXYV9sj2SEEAJMQjUorVhi7Y/r7jZqQhoaynQkxlWxG5Dvf+Y4CgYBaW1tL9ZC2MrdH/vmr0d+V5K/tkewQAgALBIPGFt1Vq4zPZR5CpBIFkZdfflmPPPKIPvrRj5bi4UqmXLZHcoAaAMAuti/N/PGPf9TNN9+sRx99VN/85jftfriSK4ftkewQKk46nfL16wEArGR7EFm3bp2uv/56LVu2LG8QSSaTSiaTw18nEgm7h2cJv2+PZIdQ4aijAWCZVKos6klsXZp5/PHHtWvXLrW1tRV0f1tbm2pqaoY/IpGIncNDgcp1h1CxqKMBYJlYTKqvl5qapNWrjc/19b48HM+2INLX16eWlhb97Gc/07Rp0wr6mfXr12twcHD4o6+vz67hoQgcoJYfdTQALFNmJ/XaFkReeeUVvf322/rLv/xLTZkyRVOmTNGOHTv0gx/8QFOmTFEqyymDlZWVCoVCGR9wXrntEJqIYupoAGBcZXhSr21B5Nprr9Vrr72mV199dfhj8eLFuvnmm/Xqq68q6MN1Lj8rlx1CE0UdDQBLFHNSr0/YVqxaXV2tyy+/POPa2WefrVmzZo25Dm8ohx1CE0UdDQBLlOFJvXRWRVH8vkNoosqx0y4AG0zkpF6P764paRDp7u4u5cMBJWPW0ezZ0yyjbmZkGKGOBkCBij2pNxYzakpGLueEw8Yhex45v4ZD7wCLUEcDYNKKOanXJ7trAul0tsjlDolEQjU1NRocHGQHDTyDzqoAJi3bTEckYoSQaNRYjqmvH7+w1Zw5OXDAkWWaYt6/CSIAALhRrtqP7m6jyVk+27cbh+uVWDHv3xSrlhn+WgcAjzBP6s3m6acL+x0e2F1DECkjnIMCwPdOn5TOqrLvfjeIxYwlmkIUugvHQRSrlgnOQQHge688Jj30v6TBXF2ORxg8ZNz/ymN2jspaZufVfAIBo6akwf0tAwgiZYBzUAD43umT0q87pHfekh67Pn8YGTxk3PfOW8bPnT5ZmnFOVr7Oq6Z0+szuGpcjiJQBzkEB4HtnVUm3PiPNqJeO954JI6mUUdjZ2Wl8TqXOhJDjvcb9tz7jneWZQms+Wls900eEGpEywDkoADyp2I6hNWFpzbNnQsYPG6X/e1Lad/jMPZfMkf5PlXRqwAgha541fs4rCq35WLHC3nFYiBmRMsA5KAA8JxYz+mQ0NUmrVxuf6+vzN+kyw8jUWiNsfDohhf7cCCwUML4+NWB8v9AQkm1WxSlm59XRzc5MHqoNMRFEyoB5DorZanysgCorI5yDAsAdJtsx9Jw6YybknSFpZoV069lSOGh8nllhXP/p+8Z9hYxlIoHILsV0XvUIgkgZMM9B+fNXo78ryd3noKTTKR0/3q2jRzt1/Hg3RbWAn5m7QrL12jSvtbbmnpXo6TGWYza/dyaM3DYihGx+T3ojbtyXi1tbqEejUleXNDfzOAmFw8Z1j9SGmAgiZcKr56AMDMS0c2e9du9u0t69q7V7d5N27qxnuzHgV/l2haTTUl9f7hBhFnQm0tJTo3bDPHXSuD7yvmysCER2ikal3l6jc+qWLcbnAwc8F0IkilUdV8pOp7W1UZ133grPdFY1e5+M3nZs9j5xc4ACMEGF7grJdZ9Z0BkKSJ8dtRvms1XGjEginbvws5hA5EALdUm5O696CEHEQU50Og0Egpoxo9GW322l/L1PAtq/v1XnnbfCtUEKwAQUuisk130NDcbumE8nzizHPHXSCCFmzcjzNbkLOq0IRCgISzMOodNpbvQ+AcqUFbtC/thvbNEdWRNyKJVZM3LLNOO+8VgRiFAQgogD6HSaH71PgDI12V0hZrMyc4vuL0NnakISaWMmxNzam6sDqw+3yboVQcRm2XZ88Nd+fvQ+AcrYRHeFjO6Yuq5b2nMws6Dz9f9nXB/dgXU0H26TdatAOp2tJNgdEomEampqNDg4qFAo5PRwijZeDUhtbbMOHWrP+/MLF27R7NmrbByhe6XTKe3cWa9kMq7sM0cBVVaGtXTpAWpEAL8qprPq6ZPGAXbvvFVYx9SRoWXmhdKd/5W9zXssZuyeGVm4GokYIcSDO1RKpZj3b4KITcbb8WH07SjsKb/iiu2eKCy1y5nnUMp8zoy/Rtg1AyDDK48ZB9jd+kxhHVMHD0mb/1q6pkVatGb8+4ptNQ+CiNPO/DWfa/klKGlI/LWfW/ZZpYguvridEAJgrNMnizvArtj7J6OMAk0x799s37VB/hoQSTILUUfPkLi/02kpea33iR+VstcNMGnFhopShZBsSzzhsFGHUuZLPAQRGxS6k2Pu3Fb94Q9dWfqI8Nf+SF7pfeJHTvS6AXzHbBU/egHCbBXvwbbsVmJpxgbHj3dr9+6mvPddccV2nXtuA39twpVy1zlRowMUJJUyDskbr0trIGDMjBw44KtlmmLev9m+a4NiTrs1/9qfPXuVZsxoJITAFeh1A1jEirNzfI4gYgOvn3YL0OsGZSeVkrq7pc5O47NVh9nRKj4vgohNvHraLeyTrbmdW9HZFmUlFjOWT5qapNWrjc/19cb1yaJVfF4Uq9qIHR/lJdfuEq8VfdLZFmXD7kJSs1V8PD72MaQzNSJl3CqeYlXAArmChiTPFX3S2RZloVSFpGbYkTLDiNkq3oe7ZihWRUl5acnBDrlPUr5R+/b9nbxW9EmdE8pCqQpJJ3p2TplgaQaT4rUlB6sVsrvkgw+O5foNw0WfbuuVYtY5Zf/vS68b5OGFLqKlLCSNRqUVK9z/nDiAIIIJG6/PhDET0OzKJQerFdZFNz+3Fn1S54QJ8UoX0VIXkgaDUmOjNb/LR1iawYTQZ8JgVYBwc9EnvW5QFLMeYvSSh1n8acVOFKuYhaSBcXo+BQLGSbtlXEhaCgQRTAh9JgyTDxBnmtsBnpdKGTMh2fZAmNdaW63r0TFZwaAxSyONDSPm1+3tLJ/YjCCCCaHPhKGQLrpTpsz68/cp+oTPebGLKIWkjqNGBBNS6EzAWWedb/NInGXuLjFqZbKfpHzJJT+SJIo+4X9e7SJaboWkLiskJohgQsyZgPH7TBjeeGONPvQhf++gKXR3CUWf8D0vdxEtl0JSFxYS09AME3Zm14w0fhhxb9Muq+XqrAqUBbNBWL4uoj47adYzxusia0NjNRqaoSTMmYCpU+fkuKt8dtCwuwRlj+JP93JxITFBBJNSWxvVwoWb89xVHjtoAIjiT7dycSExNSKYtFOn3i7wPpcVqAGwR7kVf3qBiwuJCSKYNE5qBTBGuRR/eoWLC4kJIpi0/DtojJNaadoFIK9sW0slZlcmy+wim6+Q2IEusgQRTFohvTRo2gUgr2xbS2fNMj4fG3F4pBvPrXE7s5C4udkIHSPDiMOFxLYWq7a1tenqq69WdXW1zj//fN1www3at2+fnQ8Jh5g7aCorMwvUKivDZbF1F8AkjXdGzbFjmSFEcue5NV7g0kJiW/uILF++XJ///Od19dVX64MPPtC9996r119/Xb/73e909tln5/15+oh4D700ABTN7D+Sa1fHaMX0JHFZJ1HHleD5KOb9u6QNzQYGBnT++edrx44d+qu/+qu89xNEAKAMdHdLTU0T+9nt23MXxbqwk2g5cG1Ds8HBQUnSzJkzs34/mUwqkUhkfAAAfG4yW0Zz/ex4yz0s7bhKyYLI0NCQWltbdc011+jyyy/Pek9bW5tqamqGPyKRSKmGBw9Lp1M6frxbR4926vjx7pJ0cHXiMQHfmsyW0fF+1sWdRJGpZEszd955p5577jn953/+p8LhcNZ7ksmkksnk8NeJREKRSISlGYxrYCA2zmFz9h2058RjAr6W74yabPLViLzwgrRsWf7fk29pBxPiuqWZu+66S7/4xS+0ffv2cUOIJFVWVioUCmV8AOMxD90bGQgkKZmMa8+eZg0MWD/t6sRjAr6X64yabPJtN43FpJUrC3tsBzqJIpOtQSSdTuuuu+7SU089pRdffFELFiyw8+FQRtLplPbvb1H2BmppSWm9+eYdGho6VcLHLI/D/QBbjLe1dNasM71ETLm2m5p1Ie+8U9jjOtBJFJlsbWi2bt06bdmyRU8//bSqq6t15MgRSVJNTY2qqqrsfGj43Lvv9oyZlRjt9OkBvfRSWB/+8MOWLJnkf8wzh/vNmNE46ccDys54Z9RIhW03zVUXMpqDnUSRydYg8tBDD0mSGketv23atElr1qyx86Hhc4UeoHf69ID27Gm2pKlaoY/J4X7AJIx3Rk0hdRz5TpgdzaFOoshkaxApYYsSlJliD9Dbv79V5523YlLN1TjcD3Chkc25fve7wn5m1izpRz+ij4hLlOVZM3T/9L78B+2NZM2SyZnHzL8klA+vQcAC2ZqVFeKJJ6Rrr7VnTCha2QURtl76Q+ZBe4WZ7JJJIBDURRc9oN/97qac9/3P//yDamuj4wYLXoOABcyi1GJm3s26ELbrukpJO6s6ja2X/mIetHfWWecVdL8VSyaFPJY5+5INr0HAAsUUpZocPmE2p1TKaHPf2Wl8LrMma2UTRNh66U+1tVF9/ONxnXVWbY67AqqsjOjccydfHT+ZglVeg4BFii1KlRw/YXZcsZjRzK2pSVq92vhcX19W7efLJogUs/US3lJRMVUf/vDDkgJ//hjJ+Prii9stqcGYTMEqr0HAIoU2IfvqV6UtW4zuqQcOuDOEcBZO+QQRtl76m7lMU1mZ2QypsjJsydZdk1mwOjbwmMaffeE1CFik0CZk114rrVpl1IS4cTmGs3AklVGxKlsv/a+2Nqrzzlth626UzCLZgDKXWXLPvvAaBCzS0GAstYx3No0XmpXlW15Kp6W+PuM+nxfXls2MyGT+koWzijnpNhAIasaMRs2evUozZjTasiV2orMvvAYBi+Q6m8bNRakjFbq8VAZn4ZTNjMhk/pKFc9y61XUisy+8BgELmWfTjO4jEg4bIcRt9SCjFbq8VAZn4QTSLm5/WswxwoXK/sYW0cUXt9PDwWXMra5jd5kYb9pW1n6UEq9B+NLIDqe5zoPxy+NOVipl7I7Jt7x04IA3/j2jFPP+XXZBRKKrpRek0ynt3FmfY5dJQJWVYS1desCT/+14DcJXsnU4DYeN5RO3z0w4ydw1I2WGEXN5yY3bjQtEEIHnHT/erd27m/Led8UV2znpFnDSeB1OffBmWhLZQlwk4o3lpRyKef8umxoReAtbXQEPyLcFNRAwtqCuWOHJ5YWSiEaN58eLy0sWIYjAldjqCngAW1CtEQyW9fNTNtt34S1sdQU8gC2osABBBK5kbnX981ejvyuJra6A49iCCgsQROBapWrbDmCCzA6no5uKmQIBo/DSzR1O4ThqROBqpWjbDmCCzA6nzc1G6Mi2BdXtHU7hOGZE4HqlaNsOYILMDqdzM2cuFQ6zdRcFYUYEADA5bEHFJBBEAACTV+ZbUDFxLM0AAADHEEQAAIBjCCIAAMAxBBEAAOAYgggAAHAMQQQAADiGIAIAABxDEAEAAI4hiAAAAMcQRAAAgGMIIgAAwDEEEQAA4BiCCAAAcAxBBAAAOIYgAgAAHEMQAQAAjiGIAAAAxxBEAACAYwgiAADAMQQRAADgGIIIAABwDEEEAAA4hiACAAAcU5Ig8sMf/lD19fWaNm2alixZot/85jeleFgAAOBytgeRJ554Qnfffbc2bNigXbt26YorrtB1112nt99+2+6HBgAALmd7EPn+97+v22+/XWvXrtVll12mhx9+WNOnT9c///M/2/3QAADA5WwNIqdOndIrr7yiZcuWnXnAigotW7ZML7300pj7k8mkEolExgcAAPAvW4PIH/7wB6VSKc2ePTvj+uzZs3XkyJEx97e1tammpmb4IxKJ2Dk8AADgMFftmlm/fr0GBweHP/r6+pweEiz0/umUrfcDALzH1iBy3nnnKRgM6ujRoxnXjx49qgsuuGDM/ZWVlQqFQhkf8IfO3xzU8vZf6fC7Jwu6//C7J7W8/Vfq/M1Bm0cGAHCSrUFk6tSpWrRokV544YXha0NDQ3rhhRf08Y9/3M6Hhou8fzqlR3b8j3qP/Umf/9HOvGHk8Lsn9fkf7VTvsT/pkR3/w8wIAPiY7Uszd999tx599FFt3rxZe/fu1Z133qn33ntPa9eutfuh4RLTzgpqy+1LNW/mdB1850wYSaWk7m6ps9P4nEqdCSEH3/mT5s2cri23L9W0s4JO/xMAADaZYvcDfO5zn9PAwIC+9rWv6ciRI7ryyiu1bdu2MQWs8Lc551bp8b9bOhwy/vqBnXr78aU69GbV8D3hD5/U+Z/fqWNJI4Q8/ndLNefcqhy/FQDgdYF0Op12ehDjSSQSqqmp0eDgIPUiPnH43ZP66weMsHH6+HQd7Vyq1IkqBatPavaqnTprxp80q3K6nvkSIQQAvKqY929X7ZqB/82urtLbjy/V6ePTddaMP2n2qp2qnPvOcAg5fXy6Bp5YqtnVhBAAKAcEEZRUT4906M0qHe08E0Yu+JuXhkPI0c6l6ttXpZ4ep0cKACgFgghKqr/f+Jw6UaVjz16R8b1jz16h1ImqjPsAAP5GEEFJ1dUZn4PVJzXr+t0Z35t1/W4Fq09m3AfAPqmhlLp7u9X5Wqe6e7uVGmKrPErP9l0zwEgNDcbumKHGMzUhx569QrOu3z1cMxLcsVQNDdSIAHaK7Y2pZVuLDiUODV8Lh8LqWN6h6MKogyNDuWFGBCV19ISxRXdkTUgyPjOjZqT2czt19ERhHVgBFC+2N6bmrc0ZIUSS4om4mrc2K7Y35tDIUI4IIigZs1nZsaSxRbeie+lwTUjqRJWCO5ZqVuV0HUsW1oEVQPFSQym1bGtRWmM7N5jXWre1skyDkiGIoCRGd0x95ktL1fu7Km3fLm3ZIm3fLh3YU6VnvjS2AysA6/Qc7BkzEzJSWmn1JfrUc5CtaygNgghs9/7plFY/eiaEmB1Tg0GpsVFatcr4HAye6cBqhpHVj+7krBnAQv0nCtuSVuh9wGQRRGC7aWcF9YVPXqT6WYW1bTfDSP2s6frCJy/irBnAQnXVhW1JK/Q+YLJo8Y6Sef90qqhQUez9APJLDaVU31GveCKetU4koIDCobAOtBxQsIL/+8PE0OIdrlRsqCCEANYLVgTVsbxDkhE6RjK/bl/eTghByRBEAKDMRBdG1bWyS3NDczOuh0Nhda3soo8ISoqlGQAoU6mhlHoO9qj/RL/qquvUMK+BmRBYopj3bzqrAkCZClYE1Vjf6PQwUOZYmgEAAI4hiAAAAMcQRAAAgGMIIgAAwDEEEQAA4BiCCAAAcAxBBAAAOIYgAgAAHEMQAQAAjiGIAAAAxxBEAACAYwgiAADAMQQRAADgGE7fBVwslZJ6eqT+fqmuTmpokIKc0g7ARwgigEvFYlJLi3To0Jlr4bDU0SFFo86NCwCsxNIMXCuVkrq7pc5O43Mq5fSISicWk5qbM0OIJMXjxvVYzJlxAYDVCCJwpVhMqq+Xmpqk1auNz/X15fEGnEoZMyHp9NjvmddaW8srmFkhNZRSd2+3Ol/rVHdvt1JDPIGAG7A0A9cxZwNGvxGbswFdXf5emujpGTsTMlI6LfX1Gfc1NpZsWJ4W2xtTy7YWHUqceWLDobA6lncoutDHLyYHpYZS6jnYo/4T/aqrrlPDvAYFKyhwwljMiMBVmA0wClOtvK/cxfbG1Ly1OSOESFI8EVfz1mbF9pbBNFuJxfbGVN9Rr6bNTVodW62mzU2q76jnuUZWBBG4SjGzAX5VV2ftfeUsNZRSy7YWpTU22ZrXWre1skxjIYIfikUQgaswG2Bs0Q2HpUAg+/cDASkSMe5Dbj0He8a8IY6UVlp9iT71HPRxsi0hgh8mgiDiAeW0e4TZAKNPSEeH8b9HhxHz6/Z2+okUov9EYYm10PuQm9+CHwXOpUGxqsuVWy8JczYgHs9eJxIIGN/3+2xANGoU5Wb7b9/e7s//9naoqy4ssRZ6H3LzU/CjwLl0mBFxsXLsJcFswBnRqNTbK23fLm3ZYnw+cIAQUoyGeQ0Kh8IKKPs6V0ABRUIRNczzebItEb8Ev1LVuTDjYgik09n+7nSHRCKhmpoaDQ4OKhQKOT2ckkqljL4Z4xVumjMDBw74800520xQJMJsAIpnvqlIyqhdMMNJ18ou/sK1SGoopfqOesUT8ax1IgEFFA6FdaDlgGu38pr/hvGWmKz6N/h9xqWY929mRFyq3HePMBsAq0QXRtW1sktzQ3MzrodDYUKIxYIVQXUsN6Y0R89CmV+3L293bQiRSlPnws6iTNSIuBS7R4yZHhp2wQrRhVGtuGQFDbZKwAx+2f7ab1/e7vrgZ3edS76dRQEF1LqtVSsuWVE2r0+CiEuxewSwVrAiqMb6RqeHURa8HPzsrnMpZsalXF6vBBGXYvcIAC/zavAzC5zz1blMtMDZTzuLrGJbjUhvb69uu+02LViwQFVVVbrooou0YcMGnTp1yq6H9BV2jwBA6dld5+KXnUVWsi2IvPHGGxoaGtIjjzyiPXv26IEHHtDDDz+se++9166H9B2zl8TczBo7hcP+P/gNAJxiZ4EzW8rHKun23Y0bN+qhhx7SW2+9VdD95bx9d6RUytgd099v1IQ0NDATAgB2s+sE4XLYUl7M+3dJa0QGBwc1c+bMcb+fTCaVTCaHv04kEqUYluuxewQASs+uOhev7yyyWsmCyP79+/Xggw/q/vvvH/eetrY2feMb3yjVkAAAcISXdxZZreilmXvuuUf33Xdfznv27t2rSy+9dPjreDyuT37yk2psbNSPf/zjcX8u24xIJBLxxNIMyycoJV5vgPPsWrrxg2KWZooOIgMDAzp27FjOey688EJNnTpVknT48GE1NjZq6dKleuyxx1RRUXh9rFdqRMrtYDo4i9cb4Dy/t2ifLFuDSDHi8biampq0aNEi/cu//IuCRf7J5oUgYh5MN/pZNLfYsrsFVuL1BjjPLDYd3WfET8Wmk+WKIBKPx9XY2Kj58+dr8+bNGSHkggsuKOh3uD2IlPvBdOXG6eUQXm+A80p1KJ7XueLQu+eff1779+/XCy+8oHA4rLq6uuEPvyj3g+nKSSxmhICmJmn1auNzfb1xvVR4vQHOK8WheOXGtiCyZs0apdPprB9+wcF05cFcDhkdAuJx43qpwgivN8B5tGi3nm1BpBxwMJ3/pVJGYWi2/Gxea2017rMbrzfAebRotx5BZBLMg+lGnwVjCgSkSISD6bzMTcshvN4A59Gi3XoEkUngYLqxUimpu1vq7DQ+l2KmwE5uWg7h9QY4L9eheKbJHIpXjggik8TBdGe4oaDTam5bDuH1BjjPbNE+s2rskSXZriG3kh56Vyy3b98dyemtnU4br7+F6ec/N76fTb7nzsnn1twyG49n/7c5tWW23F9vgNNie2O6ceuNY67TS8Tgij4iVvBSECln+fpbSMabZGendNNNmdfzdQl1QxdRM2RJmWGEJmJAeaKXSH6u6COC8pGvoFMywsrKlZnLNPm2xX75yxPbNmt1nQrLIQBGopeItQgimLRiCjXNra6FbIv9/veL3zZrV51KNCr19krbt0tbthifDxwghADliF4i1pri9ADgfcUUao7c6ppvW2yumYyR22YbG41r49WpmDMok529CAbPPBaA8kUvEWsxI4JJM/tbFKq/37rtrubvcVPjMQD+ZlUvkdRQSt293ep8rVPdvd1KDZXn/4MiiGDSRva3KERdnXXbXc3f46bGYwD8LVcvEfPrfL1EYntjqu+oV9PmJq2OrVbT5ibVd9QrttfD/Q4miCACS0SjxhbdXFtIR3b+LKRLaKG/S3JX4zEA/mf2EpkbyqxiD4fCebfuxvbG1Ly1eUzBazwRV/PW5gmFES/PrlAjAss0Nxs7VVauHPu9bJ0/OzqMnwkEsm+Lvftu6f77jf+d7fsjf5fbGo8B8L/owqhWXLJCPQd71H+iX3XVdWqY15BzJiQ1lFLLthalNXYdOa20AgqodVurVlyyouCtv7G9MbVsa8kINuFQWB3LOzzRy4QZEVjqppukJ58cWzOSbatrvm2x3/1u4dtmnTqHxW8t7QEUJ1gRVGN9o1Z9ZJUa6xvzhgert/7aMbtSajQ0gy2K6fxpVWfVUjces6rZGl1S4XapoVRRf/Xb6f0P3te0KdNsu99una91anVsdd77tkS3aNVHVuW8x82N1Yp5/2ZpBrYoZqtrvnsL/V3mDEu2cNDebn0IsWKrsBs6xwK5uGnav+vNLm16fZN+ct1PdMHZF+S9/8h7R3Tbf9ymtZevVfOHx54x4UTAsnLrbzGzK431jYUOseSYEYHv2D3DkK+lfaHnz4wXZrLN4DBrAieY0/6j6xmcOE/l/Q/e143/dqMOnjio8DlhbVq+KWcYOfLeEa3dtlaH/nhI86rn6cn//WTGzIhTAcucxYgn4lnrRIqZxbBydsVqtHhHWTNnUFatMj5b/YZtxVbhYvqe+PFUY7hfvqJKSWrd1lqy3RnTpkzTT677icLnhHXoj4e0dttaHXnvSNZ7R4aQ8Dlh/eS6n4wJIU7VVVix9dfkl8ZqBBGgSFZsFS40zHzrWxM7bweYLDeep3LB2Rdo0/JNOcPI6BAyeubEDQFrMlt/R7KqsZrTCCJAkazYKlxomOnooFssnOHW81RyhZF8IURyT8CKLoyqt6VX22/dri3RLdp+63YdaDlQ1LKQlbMrTiKIAEWyYqtwoWHmnXfG/x7dYmEnN0/7Zwsjr779at4QIrkrYBW79Tcbq2ZXnMSuGaBIZkv7XM3YRjZby8YMM/F49hmPQECaMSN3EDHRLRZ2MKf98xVVOjXtb4YRM3zc8twtkpS3kNXNAWuiJtJYzU2YEQEmIF8ztnxbb0eezzN6ZsX8uqWlsLHQLRZ28MK0/wVnX6C2hraMa20NbTl30/ilrmI0K2ZXnEIQASYoGpV6e6Xt26UtW4zPBw4U3v8jX5j5ylec6RYLmNw+7X/kvSNa37M+49r6nvXj7qaRvBGwyg19RACH5eoRUupusUA2buqsahpdmNrW0Kb1Pevz1oiYsvURiYQial/e7njA8oNi3r8JIh5BQ6vyla37aiRifbdYwCvG2x1TyK6ZkdwYsPyCIOIztAEHQRTF8PMbbL6wUWwYgT0IIj5STBtwuBthAqVgZetytwWaQkMGYcR5BBGfsOpMEziPWS2UgpVnw7jpsDvJ+rNmJspt4cytCCI+0d1tnCuSz/bthZ90azX+ys+PWS2UgpVHwrvpsLuRrD59t1huC2duRhDxic5O45CzfLZsMQ54KzX+ys+PWS2USndvt5o25//LZfut23MeCW9loLHD+x+8X9TMRrH3j8et4cytOH3XJ6w408Qu5l/5HMaWmxUn9QKFsKp1uV1nsaSGUuru7Vbna53q7u2e8KFyxYYKq5ZjnD4oz89o8e5ihbQBD4dL39Aq3xH2gYBxGNuKFfyVb8VJvUAhrGpdXkygKbRewutLGsWEs1yzTciOGREXK6QNeL4zTezAX/mFc/OsFvzFqtblhQaa37/ze9V31Ktpc5NWx1araXOT6jvqFdubOR1qLmmMfiOPJ+Jq3to85n43ctNBeX5EEHG5yZ5pYgf+yi+cFSf1AoWwqnV5IYFmVtUsfb3763nDRTFLGlYt3djBjwfluQlBxAMme6aJ1fgrv3BundWCP1lxNky+QGMGiELCRaFLGt/q+VZBsytO8etBeW7BrhkUzdwJkq92hZ0gZ9CmHaVkRa+L8c5i+du//Ftt6N6Q9+e337pd/Sf6tTpWwNa/LNy2G8VcYpIyQ5jbxukWbN+F7TiMrXj0XIHXZAs0W/dsLShcbIluUV11XUFbisfj9Fbh0Tgor3DFvH+zawYTYtauZOsjwl/52QWDzjWeAyYiWBEcswukmHoJc0kjnohnXcrJx227UaILo1pxyQo6q1qMIIIJi0aNLbr8lQ+Uj3zhwpzFMN+gO5Z3qHlrc0Z9iXlfoeHETbtRsoUzr3JLu3qKVTEp5l/5q1YZnwkhgL8VuzsnVwHtNxq/UdBjshvFerG9MdcUCFMjAgAoWrH1Etn++pak+o76vLMrbqkR8YtStKunWBUAYDurduewG6V0SnWWEGfNAABsZ9ZLrPrIKjXWN07ojcuK3iconF1nCU1GSYpVk8mklixZot27d+u3v/2trrzyylI8LADAA9iNUjpubFdfkiDy5S9/WXPmzNHu3btL8XAAAI/x024UN3Nju3rbl2aee+45/fKXv9T9999v90MBAIAc3Niu3tYgcvToUd1+++366U9/qunTp+e9P5lMKpFIZHwAAABrWHU4opVsCyLpdFpr1qzRHXfcocWLFxf0M21tbaqpqRn+iEQidg0PAICy5LYC4aK3795zzz267777ct6zd+9e/fKXv9TWrVu1Y8cOBYNB9fb2asGCBTmLVZPJpJLJ5PDXiURCkUiE7bsAAFjMzs6qtvYRGRgY0LFjx3Lec+GFF2rlypV65plnFBhx9nkqlVIwGNTNN9+szZs3530s+ogAAOA9rmhodvDgwYwaj8OHD+u6665TV1eXlixZonA4nPd3EEQAAPAeV5y+O2/evIyvzznnHEnSRRddVFAIAQAA/kdnVQAA4JiSNDSTpPr6ern4WBsAAOCAkgURoNRSKamnR+rvl+rqpIYGKUjHaABwFYIIfCkWk1papEMjznYKh6WODinKGVoA4BrUiMB3YjGpuTkzhEhSPG5cj8WcGRcAYCyCCHwllTJmQrKVI5nXWluN+wAAzmNpBr7S0zN2JmSkdFrq65MefFCaPZvaEQBwGkEEvtLfX9h9X/rSmf9N7QgAOIelGfhKXV3xP0PtCAA4hyACX2loMGY4AoH895qoHQEA5xBE4CvBoLHMIhUfRvr6jBoTAEDpEETgO9Go1NUlzZ1b/M8WWmMCALAGQQS+FI1Kvb3S9u3Sli3SAw8U9nMTqTEBAEwcu2bgW8Gg1Nho/O9USvre94zC1Gw9RgIBo7akoaGkQwSAsseMCMpCrtoR8+v2dvqJAECpEURQNsarHQmHjev0EQGA0mNpBmUlGpVWrOBUXgBwC4IIys7I2hEAgLNYmgEAAI5hRgQAgDKSGkqp52CP+k/0q666Tg3zGhSscG59miACAECZiO2NqWVbiw4lzhxTHg6F1bG8Q9GFzlTsszQDAEAZiO2NqXlrc0YIkaR4Iq7mrc2K7XXm5E+CCAAAPpcaSqllW4vSGtvR0bzWuq1VqaHSn/xJEAEAwOd6DvaMmQkZKa20+hJ96jlY+pM/CSIAAPhc/4nCTvQs9D4rEUQAAPC5uurCTvQs9D4rEUQAAPC5hnkNCofCCiiQ9fsBBRQJRdQwr/QnfxJEAADwuWBFUB3LjZM/R4cR8+v25e2O9BMhiAAAUAaiC6PqWtmluaHMkz/DobC6VnY51kckkE6nx+7lcYlEIqGamhoNDg4qFAo5PRwAADyvFJ1Vi3n/prMqAABlJFgRVGN9o9PDGMbSDAAAcAxBBAAAOIYgAgAAHEMQAQAAjiGIAAAAxxBEAACAYwgiAADAMQQRAADgGIIIAABwDJ1VAQCeUooW5SgdgggAwDNie2Nq2daiQ4lDw9fCobA6lnc4dmgbJoelGQCAJ8T2xtS8tTkjhEhSPBFX89ZmxfbGHBoZJoMgAgBwvdRQSi3bWpTW2APjzWut21qVGkqVemiYJIIIAMD1eg72jJkJGSmttPoSfeo52FPCUcEKBBEAgOv1n+i39D64h61B5Nlnn9WSJUtUVVWlGTNm6IYbbrDz4QAAPlVXXWfpfXAP23bNPPnkk7r99tv17W9/W5/61Kf0wQcf6PXXX7fr4QAAPtYwr0HhUFjxRDxrnUhAAYVDYTXMa3BgdJgMW4LIBx98oJaWFm3cuFG33Xbb8PXLLrvMjocDAPhcsCKojuUdat7arIACGWEkoIAkqX15O/1EPMiWpZldu3YpHo+roqJCV111lerq6vSZz3wm74xIMplUIpHI+AAAQJKiC6PqWtmluaG5GdfDobC6VnbRR8SjbJkReeuttyRJX//61/X9739f9fX1+t73vqfGxka9+eabmjlzZtafa2tr0ze+8Q07hgQA8IHowqhWXLKCzqo+UtSMyD333KNAIJDz44033tDQ0JAk6Stf+YpuvPFGLVq0SJs2bVIgENDPf/7zcX//+vXrNTg4OPzR19c3uX8dAMB3ghVBNdY3atVHVqmxvpEQ4nFFzYj8wz/8g9asWZPzngsvvFD9/cb2qZE1IZWVlbrwwgt18ODBcX+2srJSlZWVxQwJAAB4WFFBpLa2VrW1tXnvW7RokSorK7Vv3z594hOfkCSdPn1avb29mj9//sRGCgAAfMeWGpFQKKQ77rhDGzZsUCQS0fz587Vx40ZJ0k033WTHQwIAAA+yrY/Ixo0bNWXKFN1yyy06efKklixZohdffFEzZsyw6yEBAIDHBNLp9NjOMC6RSCRUU1OjwcFBhUIhp4cDAAAKUMz7N2fNAAAAxxBEAACAYwgiAADAMQQRAADgGNt2zVjBrKPlzBkAALzDfN8uZD+Mq4PIiRMnJEmRSMThkQAAgGKdOHFCNTU1Oe9x9fbdoaEhHT58WNXV1QoEjGOeE4mEIpGI+vr62NJrA55fe/H82ovn1148v/by0/ObTqd14sQJzZkzRxUVuatAXD0jUlFRoXA4nPV7oVDI8/+h3Izn1148v/bi+bUXz6+9/PL85psJMVGsCgAAHEMQAQAAjvFcEKmsrNSGDRtUWVnp9FB8iefXXjy/9uL5tRfPr73K9fl1dbEqAADwN8/NiAAAAP8giAAAAMcQRAAAgGMIIgAAwDG+CCLPPvuslixZoqqqKs2YMUM33HCD00PynWQyqSuvvFKBQECvvvqq08Pxhd7eXt12221asGCBqqqqdNFFF2nDhg06deqU00PzrB/+8Ieqr6/XtGnTtGTJEv3mN79xeki+0NbWpquvvlrV1dU6//zzdcMNN2jfvn1OD8uXvvOd7ygQCKi1tdXpoZSM54PIk08+qVtuuUVr167V7t279etf/1qrV692eli+8+Uvf1lz5sxxehi+8sYbb2hoaEiPPPKI9uzZowceeEAPP/yw7r33XqeH5klPPPGE7r77bm3YsEG7du3SFVdcoeuuu05vv/2200PzvB07dmjdunXauXOnnn/+eZ0+fVqf/vSn9d577zk9NF95+eWX9cgjj+ijH/2o00MprbSHnT59Oj137tz0j3/8Y6eH4mv//u//nr700kvTe/bsSUtK//a3v3V6SL713e9+N71gwQKnh+FJH/vYx9Lr1q0b/jqVSqXnzJmTbmtrc3BU/vT222+nJaV37Njh9FB848SJE+kPfehD6eeffz79yU9+Mt3S0uL0kErG0zMiu3btUjweV0VFha666irV1dXpM5/5jF5//XWnh+YbR48e1e23366f/vSnmj59utPD8b3BwUHNnDnT6WF4zqlTp/TKK69o2bJlw9cqKiq0bNkyvfTSSw6OzJ8GBwclideqhdatW6frr78+4zVcLjwdRN566y1J0te//nV99atf1S9+8QvNmDFDjY2Neueddxwenfel02mtWbNGd9xxhxYvXuz0cHxv//79evDBB/WFL3zB6aF4zh/+8AelUinNnj074/rs2bN15MgRh0blT0NDQ2ptbdU111yjyy+/3Onh+MLjjz+uXbt2qa2tzemhOMKVQeSee+5RIBDI+WGur0vSV77yFd14441atGiRNm3apEAgoJ///OcO/yvcq9Dn98EHH9SJEye0fv16p4fsKYU+vyPF43EtX75cN910k26//XaHRg7kt27dOr3++ut6/PHHnR6KL/T19amlpUU/+9nPNG3aNKeH4whXtngfGBjQsWPHct5z4YUX6te//rU+9alPqaenR5/4xCeGv7dkyRItW7ZM3/rWt+weqicV+vyuXLlSzzzzjAKBwPD1VCqlYDCom2++WZs3b7Z7qJ5U6PM7depUSdLhw4fV2NiopUuX6rHHHlNFhSv/PnC1U6dOafr06erq6srYNXfrrbfq3Xff1dNPP+3c4Hzkrrvu0tNPP61f/epXWrBggdPD8YV//dd/1Wc/+1kFg8Hha6lUSoFAQBUVFUomkxnf86MpTg8gm9raWtXW1ua9b9GiRaqsrNS+ffuGg8jp06fV29ur+fPn2z1Mzyr0+f3BD36gb37zm8NfHz58WNddd52eeOIJLVmyxM4helqhz69kzIQ0NTUNz+YRQiZm6tSpWrRokV544YXhIDI0NKQXXnhBd911l7OD84F0Oq2///u/11NPPaXu7m5CiIWuvfZavfbaaxnX1q5dq0svvVT/+I//6PsQIrk0iBQqFArpjjvu0IYNGxSJRDR//nxt3LhRknTTTTc5PDrvmzdvXsbX55xzjiTpoosuUjgcdmJIvhKPx9XY2Kj58+fr/vvv18DAwPD3LrjgAgdH5k133323br31Vi1evFgf+9jH1N7ervfee09r1651emiet27dOm3ZskVPP/20qqurh+tuampqVFVV5fDovK26unpMrc3ZZ5+tWbNmlU0NjqeDiCRt3LhRU6ZM0S233KKTJ09qyZIlevHFFzVjxgynhwbk9Pzzz2v//v3av3//mGDnwhVT1/vc5z6ngYEBfe1rX9ORI0d05ZVXatu2bWMKWFG8hx56SJLU2NiYcX3Tpk1as2ZN6QcEX3FljQgAACgPLEgDAADHEEQAAIBjCCIAAMAxBBEAAOAYgggAAHAMQQQAADiGIAIAABxDEAEAAI4hiAAAAMcQRAAAgGMIIgAAwDEEEQAA4Jj/D3+nhsbgfDNDAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i in range(80):\n",
    "    if y[i] == 0:\n",
    "        plt.scatter(X[i][0], X[i][1], c='b')\n",
    "    elif y[i] == 1:\n",
    "        plt.scatter(X[i][0], X[i][1], c='r')\n",
    "    elif y[i] == 2:\n",
    "        plt.scatter(X[i][0], X[i][1], c='y')\n",
    "    else:\n",
    "        plt.scatter(X[i][0], X[i][1], c='g')\n",
    "\n",
    "for _m in fcm.m:\n",
    "    plt.scatter(_m[0], _m[1], marker='x', s = 100)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "图像中的叉表示一个聚类中心，可以发现其实这个算法并没有正确给出三个中心，而是和C均值一样把两个类混在一起了。\n",
    "\n",
    "显然书上202和203页写的结论在这个例子上不成立，我猜测原因是FCM算法的结果受初始中心值影响极大，因为初始的中心设为了普通C均值算法的中心（可以去看代码中的`initialize`函数），最终导致了这样的结果。如果把初始的中心设置为正确的中心，估计结果就会符合预期。但是这也从侧面说明了这个算法的弊端，如果我们对数据的内在规律不了解，随机生成初始的中心，算法结果可能会很差。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 五、分级聚类方法\n",
    "\n",
    "这个方法的思想很简单，书上也只是简单的介绍了一下思想，就是通过自底向上的思想逐级合并子类为父类。分级聚类得到结果后还需要人为的划分类别，还有好多细节和指标要考虑，我也就不实现了。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 六、自组织映射神经网络（SOM网络）\n",
    "\n",
    "SOM网络是最典型的竞争型神经网络，思想是：轮流通过每个样本训练网络，每个样本会找到一个最相似的神经元作为获胜节点，在这个获胜节点的某个邻域更新神经元的权值。反复迭代直到收敛，这时每个样本都只有一个获胜节点，称为样本的像。\n",
    "\n",
    "这个算法在书上介绍得很清楚，但是有两个实现算法的关键问题书上没有提及，即更新步长$\\alpha(t)$和邻域$N_c(t)$如何计算。\n",
    "\n",
    "更新步长在实现时暂且可以用固定值来处理，但是邻域是逐渐缩小的，并且到迭代的终点会缩小成只包含一个获胜节点。邻域如何缩小书上并没有介绍，这是实现SOM的最大问题。我在网上查阅资料的时候发现SOM的资料远远少于其他算法，并且也没找到太详细的实现邻域缩小算法的细节，这里先不实现了，有时间的话回来补上这个坑。"
   ]
  }
 ],
 "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
}