{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "98ed234f",
   "metadata": {},
   "source": [
    "# 机器学习与社会科学应用"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61b2ffaa",
   "metadata": {},
   "source": [
    "# 第八章 特征工程入门与实践"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d36a17af",
   "metadata": {},
   "source": [
    "<font face=\"宋体\" >郭峰    \n",
    "    教授、博士生导师  \n",
    "上海财经大学公共经济与管理学院  \n",
    "上海财经大学数实融合与智能治理实验室  \n",
    "    邮箱：guofengsfi@163.com</font> "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "759d31a2",
   "metadata": {},
   "source": [
    "<font face=\"宋体\" >本章目录：  \n",
    "第一节 特征工程简介  \n",
    "第二节 特征理解：探索性分析  \n",
    "第三节 特征增强：清洗数据  \n",
    "第四节 特征构造：生成新数据  \n",
    "第五节 特征选择：筛选属性  \n",
    "第六节 特征转换：数据降维</font> "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3801049d-5f60-453a-b174-ca0c001f7871",
   "metadata": {},
   "source": [
    "## 第一节 特征工程简介"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90154359-a67f-4d45-8d12-10e5249b0da5",
   "metadata": {},
   "source": [
    "维基百科对特征工程的定义\n",
    ">Feature engineering or feature extraction or feature discovery is the process of using domain knowledge to extract features (characteristics, properties, attributes) from raw data. The motivation is to use these extra features to improve the quality of results from a machine learning process, compared with supplying only the raw data to the machine learning process.\n",
    "\n",
    "特征工程（Feature Engineering）是从原始数据中提取特征并将其转换为适合机器学习的格式，从而改善机器学习性能，即提高预测准确度，同时减少运行时间。业界流传着一句话：数据和特征决定了机器学习的上限，而模型和算法只是逼近这个上限而已。那特征工程到底是什么呢？顾名思义，其本质是一项工程活动，目的是最大限度地从原始数据中提取有价值特征以供算法和模型使用。\n",
    "\n",
    "特征工程的性能包含模型的预测性能与模型的元指标两方面。\n",
    "其中预测性能根据任务的不同而不同，分类任务可以使用如下指标：\n",
    "- 真阳性率和假阳性率\n",
    "- 灵敏度（真阳性率）和特异性；\n",
    "- 假阴性率和假阳性率。\n",
    "回归任务则可以使用：\n",
    "- 平均绝对误差\n",
    "- R方\n",
    "\n",
    "元指标，即不直接与模型预测性能相关的指标，包括：\n",
    "- 模型拟合、训练所需的时间\n",
    "- 拟合后的模型预测实例的时间\n",
    "- 需要持久化（永久保存）的数据大小\n",
    "\n",
    "为了提取知识和做出预测，机器学习使用数学模型来拟合数据。这些模型将特征作为输入。\n",
    "\n",
    "特征就是原始数据某个方面的数值表示。在机器学习流程中，特征是数据和模型之间的纽带。\n",
    "  \n",
    "经典特征工程包括特征理解、特征增强、特征构建、特征选择和特征转换等五个步骤，从而为进一步提高机器学习性能作准备。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c562b6d-7524-48f2-bc51-e4b0646e11a9",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 第二节 特征理解：探索性分析"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2364b40-364e-439e-82cb-28c208d7c8dc",
   "metadata": {},
   "source": [
    "特征理解就是理解我们的数据特征，对数据有初步、总体的认识。具体而言，就是考察数据规模（行数、列数）、数据类型（结构化还是非结构化，定量还是定性）、数据等级（定类等级、定序等级、定距等级、定比等级）、数据的统计信息、数据的分布、数据之间的关系等。下面就通过一个经典的鸢尾花案例演示如何进行探索性分析。探索性分析（EDA，Exploratory Data Analysis）本质上就是刚拿到数据，对数据的整体规模、结构、数据之间的关系等还不清楚，需要进行探索，从而为进一步的数据分析打下基础。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a920567-ec02-41f2-b519-e782c2660c03",
   "metadata": {},
   "source": [
    "### 2.1 查看数据概况"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b7d17f3-59e3-430a-94f9-5aee40425beb",
   "metadata": {},
   "source": [
    "这部分主要是查看数据的规模、列名、非空行数、数据类型（数值型还是字符型），从而对数据有一个总体的认识。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48e0d3df-ab99-450c-a60d-2d96601061e9",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "import pandas as pd \n",
    "import numpy as np \n",
    "import matplotlib.pyplot as plt \n",
    "import seaborn as sns\n",
    "\n",
    "path = \"D:/python/机器学习与社会科学应用/演示数据/08特征工程入门与实践/\"\n",
    "iris = pd.read_csv(path+\"iris.csv\")\n",
    "print(iris.shape)\n",
    "\n",
    "# 查看数据基本信息:\n",
    "iris.info()\n",
    "\n",
    "# 显示缺失值数量\n",
    "iris.isnull().sum()\n",
    "\n",
    "# 查看前5行，自己可以设置\n",
    "iris.head()\n",
    "\n",
    "# 随机查看5行\n",
    "iris.sample(5)\n",
    "\n",
    "# 查看末尾五行\n",
    "iris.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "812c9db3-8e2d-4c97-89f5-c5d52d6c71b7",
   "metadata": {},
   "source": [
    "### 2.2 查看数据统计信息"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22fbc19b-2e35-4fca-870e-0a840549bf1f",
   "metadata": {},
   "source": [
    "####  描述性统计"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fd153d6-8979-44d1-ab19-6f1c83471577",
   "metadata": {},
   "source": [
    "数据等级\n",
    "\n",
    "等级 | 属性 | 例子 | 描述性统计 | 图表 \n",
    "- | - | - | - | -\n",
    "定类 | 离散, 无序 | 种类 | 频率/占比, 众数 | 条形图, 饼图\n",
    "定序 | 有序类别, 比较 | 考试等级 |　频率, 众数, 中位数, 百分数 | 条形图，饼图，茎叶图\n",
    "定距 | 数字差别有意义 | 摄氏度,年份 | 频率,众数,中位数,均值,标准差 | 条形图, 饼图, 茎叶图, 箱线图, 直方图\n",
    "定比 | 连续, 可以做加减乘除运算 | 金钱, 重量 | 均值, 标准差 | 直方图, 箱线图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "148f1ec5-85d7-4732-8867-c5453c82dd3a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 查看连续变量数据的统计特征：频数、均值、标准差、最大、最小、分位数\n",
    "print(iris.describe().T)\n",
    "\n",
    "# 查看离散变量数据的分布，判断属于什么数据等级\n",
    "print(iris.species.describe())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1f63792-2d13-41f1-87e1-1354a499fc10",
   "metadata": {},
   "source": [
    "通过这个部分的分析可知，前四列数据为数值便来给你，属于定比等级，最后1列为定类等级，即类别变量。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acaea6e9-9d77-43cf-ac05-83c1e354dc64",
   "metadata": {},
   "source": [
    "#### 变量分布"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ad3d987-2abb-4ae6-a380-6f660f69cb66",
   "metadata": {},
   "source": [
    "通过图形，我们可以直观地看到每列数据的分布"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02349eb2-a927-42d7-aeb1-a35378f4eda2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 直方图\n",
    "iris.hist()\n",
    "# 直方图设置，将数据分为20组，使用蓝色画图，并添加标签与标题\n",
    "# iris.sepal_length.hist(bins=20, color='b')\n",
    "# plt.xlabel('sepal_length')\n",
    "# plt.title('Iris Data')\n",
    "\n",
    "# 核密度图\n",
    "iris.plot.density()\n",
    "# iris.plot.density(subplots=True)\n",
    "# iris.sepal_length.plot.density()\n",
    "\n",
    "# 同时画sepal_length的直方图与核密度图\n",
    "sns.distplot(iris.sepal_length, rug=True)\n",
    "\n",
    "# 箱型图\n",
    "iris.plot.box()\n",
    "# iris.sepal_width.plot.box()\n",
    "# 画不同品种（species）的鸢尾花的花萼宽度（sepal_width）\n",
    "# sns.boxplot(x='species', y='sepal_width', data=iris)\n",
    "\n",
    "# 饼图（适用于类别变量）\n",
    "# 单独展示，否则与上图有重叠，如果分类较多，可以设置前五项\n",
    "iris.species.value_counts()\n",
    "iris.species.value_counts().sort_values(ascending=False).head(5).plot(kind='pie')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33b5ef30-9f87-44c3-8a00-84dc3cb4565a",
   "metadata": {},
   "source": [
    "#### 分组统计"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c1270fc-759a-4939-8f09-5bf989268c84",
   "metadata": {},
   "outputs": [],
   "source": [
    "iris_grouped = iris.groupby('species')\n",
    "\n",
    "# 查看相关系数\n",
    "iris_grouped.corr()\n",
    "\n",
    "# 查看所有分组的统计指标\n",
    "iris_grouped.describe()\n",
    "\n",
    "# 查看petal_length变量的统计指标\n",
    "iris_grouped.petal_length.describe()\n",
    "\n",
    "# 查看所有变量的分组平均值\n",
    "iris_grouped.mean()\n",
    "\n",
    "# 计算所有变量的分组平均值与分组标准差，使用agg()方差（表示aggregate，即合计）\n",
    "iris_grouped.agg(['mean', 'median', 'std'])\n",
    "\n",
    "# 计算不同变量的统计值\n",
    "iris_grouped.agg({'sepal_length': 'mean', 'sepal_width': 'mean', 'petal_length': 'std', 'petal_width': 'std'})\n",
    "\n",
    "# 自定义四分位距指标\n",
    "def iqr(x):\n",
    "    return x.quantile(0.75) - x.quantile(0.25)\n",
    "iris_grouped.agg(iqr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21ce4cc7-acec-410c-8985-aa0859a1bac7",
   "metadata": {},
   "source": [
    "## 2.3 变量相关性探索"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acd4e794-7775-4e86-ae2e-840a69edb788",
   "metadata": {},
   "source": [
    "###  整体相关系数矩阵"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4671bf81-735c-4746-83b2-3087e0709af0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 考察变量之间的相关系数矩阵\n",
    "corr = iris.corr()\n",
    "print(corr)\n",
    "\n",
    "sns.heatmap(corr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "842f2e7f-fba0-4c58-a652-2fe62c0a4833",
   "metadata": {},
   "source": [
    "#### 分组相关系数矩阵"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d706ac3-f91e-47ab-98f3-7cebe31a6c32",
   "metadata": {},
   "outputs": [],
   "source": [
    "iris_grouped = iris.groupby('species')\n",
    "# 查看分组相关系数\n",
    "iris_grouped.corr()\n",
    "\n",
    "sns.heatmap(iris_grouped.corr())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "391402a3-ca6c-4f6b-93bc-0d2d4014d5f6",
   "metadata": {},
   "source": [
    "#### 散点图 \n",
    "如果我们更关心变量之间的关系，而两个变量之间的散点图（scatter plot）是常用的画图工具。可用seaborn的scatterplot()函数画散点图。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "627376a1-9571-484a-8f33-93f3427d465d",
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.scatterplot(x='petal_length', y='petal_width', data=iris)\n",
    "\n",
    "# 用不同的颜色表示不同的鸢尾花品种\n",
    "sns.scatterplot(x='petal_length', y='petal_width', data=iris, hue='species')\n",
    "\n",
    "# 以不同颜色和不同图标区分鸢尾花品种\n",
    "sns.scatterplot(x='petal_length', y='petal_width', data=iris, style='species', hue='species')\n",
    "\n",
    "# 画两两变量的散点图\n",
    "sns.pairplot(data=iris, height=2)\n",
    "\n",
    "# 画部分变量的散点图，并在主对角线上画核密度图\n",
    "sns.pairplot(iris, diag_kind='kde', vars=['sepal_length', 'sepal_width', 'petal_length'], hue='species')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "399e655d-0278-45a3-8ac2-81370a32e603",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 第三节 特征增强：清洗数据"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd3086c4-3f0f-40c8-80ff-e220d8756678",
   "metadata": {},
   "source": [
    "该阶段的目的是清洗和增强数据，将数据转化为机器学习模型能够直接使用的格式。清洗数据是指调整已有的列和行，增强数据是指在数据集中删除和添加新的列。特征增强的意义是，识别有问题的区域，并确定哪种修复方法最有效，。\n",
    "\n",
    "通过特征提取，我们能得到未经处理的特征，这时的特征可能存在以下问题：\n",
    "- **不属于同一量纲**：即特征的规格不一样，不能够放在一起比较。无量纲化可以解决这一问题。\n",
    "- **信息冗余**：对于某些定量特征，其包含的有效信息为区间划分，例如学习成绩，假如只关心“及格”或不“及格”，那么需要将定量的考分，转换成“1”和“0”表示及格和未及格。二值化可以解决这一问题。\n",
    "- **定性特征不能直接使用**：某些机器学习算法和模型只能接收定量特征的输入，那么需要将定性特征转换为定量特征。最简单的方式是为每一种定性值指定一个定量值，但是这种方式过于灵活，增加了调参的工作。通常使用哑编码的方式将定性特征转换为定量特征：假设有N种定性值，则将这一个特征扩展为N种特征，当原始特征值为第i种定性值时，第i个扩展特征赋值为1，其他扩展特征赋值为0。哑编码的方式相比直接指定的方式，不用增加调参的工作，对于线性模型来说，使用哑编码后的特征可达到非线性的效果。\n",
    "- **存在缺失值**：缺失值需要补充。\n",
    "- **异常值需要处理**：异常值可能会影响模型的预测结果，因而需要处理。\n",
    "- **信息利用率低**：不同的机器学习算法和模型对数据中信息的利用是不同的，之前提到在线性模型中，使用对定性特征哑编码可以达到非线性的效果。类似地，对定量变量多项式化，或者进行其他的转换，都能达到非线性的效果。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39500c3a-74b5-4fb9-b5b8-6bb35662eaad",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 3.1 缺失值处理"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c09c02d1-4ff1-4947-bcfc-848045201ba4",
   "metadata": {},
   "source": [
    "有些特征可能因为无法采样或者没有观测值而缺失。例如位置特征，用户可能禁止获取地理位置或者获取地理位置失败，形成缺失值，缺失值的存在会对模型的预测结果产生影响，因而需要进行处理。  \n",
    "详细参考：https://zhuanlan.zhihu.com/p/537865122"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7f01722-dfaa-4840-89fb-e46074246753",
   "metadata": {},
   "source": [
    "#### 识别缺失值"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5d3a1db-0fef-46dc-9432-a195ac1de7bb",
   "metadata": {},
   "source": [
    "在原始数据中，缺失值可能直接以空值显示，这种情况下可以直接进行统计；另一种情况是空值可能以某些特殊的值表示，比如0、99、？等，这时就需要数据使用者手动进行转换。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b07a5fd3-6362-4a10-9a2f-24cacd7e53ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "df = pd.DataFrame([[1, np.nan, pd.NaT], [None, 'b', np.nan]])\n",
    "print(df)\n",
    "\n",
    "# 统计每列的空值\n",
    "df.isnull().sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33caa279-d9b3-4997-9944-23a8d7bbced2",
   "metadata": {},
   "source": [
    "自定义缺失值识别"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d0db52d-e9fe-446f-abbe-55c878e685d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "df = pd.DataFrame([[1, '?', 3], ['?', 4, '?'], [7, 8, '?']], index=list(\"ABC\"))\n",
    "print(df)\n",
    "\n",
    "# 统计自定义缺失值的数量\n",
    "df.isin(['?']).sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "157a7174-e245-42fb-87f8-bfb541d22ff0",
   "metadata": {},
   "source": [
    "#### 缺失值处理"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa3c80ac-86c8-4f84-aead-13f99c695152",
   "metadata": {},
   "source": [
    "对缺失值处理主要有两种方式：  \n",
    "- 删除缺失值，这种一般适用于数据集中缺失值较少的情况，如果缺失值比例超过10%，那么删除数据会对数据分析结果有很大影响，不合理；\n",
    "- 填充值，填充的值可以是均值、中位值、前值、后值，对于分类变量而言，可以是众数，此外还可以进行插值，如二次曲线插值，某些情况下也可以填充0，可以根据实际情况选择合适的填充值。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0f9a22c-18b7-40cf-abf5-6c73d6fbd351",
   "metadata": {},
   "source": [
    "####  删除缺失值\n",
    "\n",
    "删除缺失值，必然会导致数据量的减少，如果缺失值占数据的比例较大，比例超过了数据的10%（具体标准根据项目来定），删除数据对数据分析的结果会有很大影响，可以对其进行合理的填充。\n",
    "\n",
    "语法：dropna(axis=0, how='any', thresh=None, subset=None, inplace=False)\n",
    "\n",
    "- axis: 默认为0（'index')，按行删除，即删除有空值的行。将axis参数修改为1或'columns'，则按列删除，即删除有空值的列；\n",
    "- how: 参数默认为any，只有有一行（或列）数据中有空值就会删除该行（或列）。将how参数修改为all，则只有一行（或列）数据中全部都是空值才会删除该行（或列）。\n",
    "- thresh: 表示删除空值的阈值，传入一个整数。如果一行（或列）数据中存在少于thresh个非空值（non-NA values），则删除该行（或列）。\n",
    "- subset：删除空值是，只判断subset指定的列（或行）的子集，其他列（或行）中的空值忽略，不处理。\n",
    "- inplace: 默认为False，返回元数据的一个副本。将inplace参数修改为True，则会修改数据本身。 \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c029c6e4-8c9f-4200-8dc9-789b2998b58c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd \n",
    "\n",
    "df = pd.DataFrame([[1, 2, 3], [np.nan, 5, pd.NaT], [7, 8, 9],\n",
    "                  [np.nan, np.nan, np.nan], [None, 14, 15]], \n",
    "                 index=['one', 'two', 'three', 'four', 'five'],\n",
    "                 columns=list('ABC'))\n",
    "print(df)\n",
    "\n",
    "# 将超过2个空值的行删除\n",
    "df.dropna(thresh=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b93edcf-3eff-49d8-8537-a476f3ae79c5",
   "metadata": {},
   "source": [
    "####  填充缺失值"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f222abee-c682-4a7a-9795-74ae55dd965e",
   "metadata": {},
   "source": [
    "语法： fillna(value=None, method=None, axis=None, inplace=False, limit=None)\n",
    "\n",
    "- value: 表示填充的值，可以是一个指定值，也可以是字典，series或DataFrame\n",
    "- method: 填充的方式，默认为None。有ffill、pad、bfill、backfill四种方式。ffill和pad表示用缺失值的前一个值填充，如果axis=0，则用空值上一行的值填充，如果axis=1，则用空值左边的值填充。加入空值在第一行或第一列，以及空值前面的值全都是空值，则无法获取到可用的填充之，填充后依然为空值。bfill和backfill表示用缺失值的后一个值填充，axis的用法以及找不到填充值的情况同ffill和pad。\n",
    "- axis: 通常配合method参数使用，axis=0表示行，axis=1表示按列\n",
    "- limit: 表示填充执行的次数。如果是按行填充，则填充一次表示执行一次，按列同理。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41493423-b984-4c0c-b5a0-180ce7c6bad4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd \n",
    "\n",
    "df = pd.DataFrame([[1, 2, 3], [np.nan, 5, pd.NaT], [7, 8, 9],\n",
    "                  [np.nan, np.nan, np.nan], [None, 14, 15]], \n",
    "                 index=['one', 'two', 'three', 'four', 'five'],\n",
    "                 columns=list('ABC'))\n",
    "print(df)\n",
    "\n",
    "# 前一个值填充\n",
    "df.fillna(method='ffill')\n",
    "\n",
    "# 不同列填充不同值\n",
    "values = {'A':100, 'B': 200, 'C': 300}\n",
    "df.fillna(value=values, limit=2)\n",
    "\n",
    "# 该列均值进行填充\n",
    "df.fillna(df.mean())\n",
    "\n",
    "# 该列的中位值进行填充\n",
    "df.fillna(df.median())\n",
    "\n",
    "# 该列的众数进行填充\n",
    "df['C'].fillna(df['C'].mode()[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f933a94c-8727-4ba4-9fe6-52a2c3a8fa8f",
   "metadata": {},
   "source": [
    "####  插值"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5762078-4059-44f1-b86a-b3a690c345c0",
   "metadata": {},
   "source": [
    "除了直接使用统计量或特定值进行填充，还可以使用具体函数进行插值填充。Series与DataFrame对象都可以用interpolate()函数进行插值，默认对缺失值进行线性插值。此外可以通过参数method进行调整：\n",
    "\n",
    "- 如果数据为一个时间序列且有增长的趋势，那么使用二次曲线，即method='quadratic'可能比较合适\n",
    "- 如果数据的值接近一个累积分布函数，则method='pchip'应该比较合适 - 如果为了平滑曲线，考虑method='akima'\n",
    "- 此外，还可以插入样条曲线，即method=\"spline\", 或多项式,method=\"polynomial\"，但这两个还需要指定阶数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79fde7b1-bab0-43cf-80f7-92b02b102bfb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 插值，默认使用线性插值\n",
    "import pandas as pd \n",
    "\n",
    "df = pd.DataFrame(\n",
    "    {\n",
    "        \"A\": [1, 2.1, np.nan, 4.7, 5.6, 6.8],\n",
    "        \"B\": [0.25, np.nan, np.nan, 4, 12.2, 14.4],\n",
    "    }\n",
    ")\n",
    "print(df)\n",
    "\n",
    "# 使用线性插值\n",
    "df.interpolate()\n",
    "\n",
    "# 二次样条曲线\n",
    "df.interpolate(method='spline', order=2)\n",
    "\n",
    "\n",
    "# 二次多项式\n",
    "df.interpolate(method='polynomial', order=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9a2bc69-7802-4771-a8ed-88398266c643",
   "metadata": {},
   "source": [
    "#### 不处理"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d15cd081-a72d-489d-8c69-af541db164a9",
   "metadata": {},
   "source": [
    "LightGBM和XGBoost都能将NaN作为数据的一部分进行学习，所以不需要处理缺失值。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7334652-e86b-4696-84bf-aa2e90b5a7d4",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 3.2 异常值处理"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9a5549d-ee5e-434d-9d32-042a23ab0fb0",
   "metadata": {},
   "source": [
    "异常值，也称为噪声数据，是指那些在数据集中存在的不合理的值，需要注意的是，不合理的值是偏离正常范围的值，不是错误值。比如人的身高为2.4m，体重为1吨等，虽然这些值确实可能存在，但是过于偏离正常范围，因而属于异常值。异常值会对结果造成偏差，所以在数据分析过程中需要重视。从本质上说，根据中心极限定理，总体均值近似服从正态分布，大部分数据其实集中在一个较小的范围内，而我们要研究的目标是适合大样本的规律，因而，对于这些异常值，其实不是我们的研究对象，但是在技术上，如果我们不对异常值处理，可能使结果产生偏差，因而，我们需要对异常值进行处理。  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a6866d3-e0a5-4eeb-9c33-367fb14a0936",
   "metadata": {},
   "source": [
    "#### 识别异常值"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7dd7b5f0-e4db-44ef-ac43-e0c55e2b0c31",
   "metadata": {},
   "source": [
    "#### 简单统计分析\n",
    "\n",
    "简单的统计分析。最常用的统计量是最大值与最小值，可以用来判断这个变量的取值是否超出合理的范围。例如，使用df.describe()查看数据的统计信息是否合理，该数据为美国波士顿的住房数据为例，展示数据是否合理，相关变量的定义如下：\n",
    "CRIM：每个城镇的人均犯罪率\n",
    "ZN：超过25,000平方英尺的住宅用地所占比例\n",
    "INDUS：每个城镇的非零售商业用地比例\n",
    "CHAS：查尔斯河虚拟变量（如果地块与河流接壤，则为1；否则为0）\n",
    "NOX：一氧化氮浓度（每千万分之一）\n",
    "RM：每个住宅的平均房间数\n",
    "AGE：建于1940年之前的自住单位所占比例\n",
    "DIS：到波士顿五个就业中心的加权距离\n",
    "RAD：到径向高速公路的可达性指数\n",
    "TAX：每1万美元的全额财产税率\n",
    "PTRATIO：每个城镇的师生比例\n",
    "B：1000(Bk - 0.63)^2，其中Bk是城镇中黑人的比例\n",
    "LSTAT：人口中地位较低的人所占的百分比\n",
    "MEDV：自住房屋的中位数价值（以千美元为单位）\n",
    "具体示例代码为：："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2e20aeb-945c-4918-aa54-4d647ead7554",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 加载数据\n",
    "import pandas as pd\n",
    "from sklearn.datasets import fetch_openml\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "path = \"D:/python/机器学习与社会科学应用/演示数据/08特征工程入门与实践/\"\n",
    "\n",
    "df = pd.read_csv(path+'boston.csv')\n",
    "print(df.head())\n",
    "df.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a4af0fe-c08f-46bf-b2c5-d3b890bb16ec",
   "metadata": {},
   "source": [
    "####  3$\\sigma$原则"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bb5d9f2-fe84-41fc-b307-ff492c74ae7a",
   "metadata": {},
   "source": [
    "如果样本近似服从正态分布，根据正态分布的特点，如果样本值位于均值加上或减去3倍标准差之外的概率为0.003，属于小概率事件，因而我们可以将其视为异常值。在实践中，首先通过数据的直方图和密度图观察样本分布，然后计算z-score进行定量判断: \n",
    "$z-score=\\frac{x-mean(X)}{std(X)}$，其中,X某列数据,x为该列的某行的值."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fa860bc-ec17-4e08-bf4a-1a0c17499454",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 将样本值位于大于或小于均值加上或减去3倍标准差的值定义为异常值，并进行标记为空值\n",
    "def detect_outliers(filename, attribute, threshold=3, fill_value=np.nan):\n",
    "    \"\"\"功能:探测异常值:\n",
    "    Args:\n",
    "        filename: 为数据集,\n",
    "        attribute: 为要识别的特征,\n",
    "        threshold: 为要设定的门槛值,默认为3倍标准差,\n",
    "        fill_value: 为对异常值进行标记值\n",
    "    \"\"\" \n",
    "    outliers_ll = filename[attribute].mean() - threshold*filename[attribute].std()\n",
    "    outliers_ul = filename[attribute].mean() + threshold*filename[attribute].std()\n",
    "    # 将异常值标记为\n",
    "    filename[attribute].where((filename[attribute]<outliers_ul) | (filename[attribute]>outliers_ll), fill_value)\n",
    "    return filename\n",
    "\n",
    "# 对数据集df的'TNDUS'的异常值标记为True\n",
    "df = detect_outliers(df, 'INDUS', True)\n",
    "# 显示异常值的行\n",
    "df[df['INDUS']==True]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09ed8370-91f4-40ed-9715-1786feab9406",
   "metadata": {},
   "source": [
    "#### 箱型图分析\n",
    "\n",
    "箱型图非常适合做异常值观察的图形，箱型图的五根线从上到下依次表示为最大值、最小值、上四分位、下四分位和中位数，最大值区间：上四分位+1.5IQR，最小值区间：下四分位-1.5IQR，其中IQR=上四分位-下四分位，高于最大值或小于最低值被认为是异常值。 因此，异常值通常被定义为小于QL-1.5IQR或者大于QU+1.5IQR的值，QL为下四分位数，QU为上四分位数。"
   ]
  },
  {
   "attachments": {
    "a1156988-5cc8-4673-a839-ca766b28c527.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUMAAAFaCAIAAAALiKgDAAAgAElEQVR4nOydd3xVVbb4196n3pvbkpseQu9dUMYCShgEbKg4Y/cNFhh1sDwd/Y3giB3HEccCjMJH8WFHHyMoKoNIcRDkAYMMhKLBJKSQfvs9de/fHzu5HBLABG4gxPP93A+fw8np56y99lp77bUQpRRsbGzOcPDpvgAbG5skYEuyjU1nwJZkG5vOgC3JNjadAVuSbWw6A7Yk/6IhhJzuS7BJDrYk/0JRVRUATNO0rqSUHmdU0hb7jowtyb9QmAwnhJMt6LpeV1dnXZNYIIRYxZ5SSgjRdZ3tFYvFWBMQCoXi8Xg4HA4Gg4qiJNoF2sSpubtfILYk/0JxOp0AYJpmLBarqanBGAMAh7n09PR4PK5pGlujKAolFAAooYIgMKk2DRMADh06hFHj97Nv375IJBKPx1VV3b9/P0LI6/VijCmliqIwgbdpV/jTfQE2pw1N0yRRikQjGRkZX3311d8X/H3p0qUAsGHDhvfff3/gwIGEkJSUFMMwAoHAOeecM3HixJqaGo7j0v3p11133axZs+LxeG5uriRJ+/bu+/TTT2fNnAUAjzzyyOeff67reigUSk9Pl2XZVsWnAFuSf6GYhhkOh/1+v8fjAQCn08kLvGEagWBAluVDhw5Nnz69vr5+wIAB1dXVX3755aZNmyZOnJiRkUEpVTUVY6zreq9evSZMmMBxnKIowWDw66+/9vl8JSUlN9xwAyFk0KBBjzzyCMdxCKHTfbudH1uSf6EQSnxeHwAEAoFt27Zt27atoaHh+++/r6mp6datm8vlGj58+Nq1a3fs2DF58uRvv/02LS2ttLS0V69eACCKoqqqpmlSSr/4/Iufin96//33fT7fHXfc4XA4MjMzFyxY4PP5EEKKogAAx3EAwOSZ6WdbtpOObSf/EqGURqNRjucURVmzZs3fF/x937595eXlTz311Nq1aysqKurq6hwOR0lJSWFhoSAIlNL9+/f36tXLNMy5c+f26dNn8+bNN9xwQ3Z29s7/7Pzyyy979uwpy/Jzzz23ZMmSxx9//K233lq3bp2u65Ik2R7vU4Mtyb9EKKUOh6OhoUEUxfHjxy9duvTWW2/t1q3b8uXLH3nkkeHDh3u9Xl3Xg8FgLBYDAIxxnz59AIDjuenTp+/fv1+SpMWLFxcXF2/YsOGtt94aNGiQw+G47bbbsrKyBgwY4PF4vvnmm1AohBBKqF/bd92uJEeS7Td0pmB9U6qqUkI9Ho9JzNTU1Pz8/JqaGr/fv3PnTgAQRdHpdHbp0iUSiQiCcPDgQeaydrvde/bsUVX1vffeczgciqL4/f7NmzcfPHjwnXfe2bJlS2Fh4datWwVBYBa4IAin62Z/USTHTqaUGobB3pmiKHPnzvX5fNFoVFVVjuOYmWRz2uF5/oEHHlAUheO4SCTi8Xg4nps9e3ZWVtYPP/ywZs0aQRACgcCkSZM0TSOEiKIYDocdskPXdbfbbZhGVXVVbm7uhx9+OHHixGg0umDBAkEQwuFwJBLRdZ3n+bS0NEmSYrGYy+USBEFRlGg0ygxyzNkdwHYkOZKMMWajjuFwuLq6+pNPPrnhhhscDofb7WY+D5uOgMfjQQhVVlbm5eX5/f66urpQKHT77beXl5crimIYxu9///uGhoZoNKrrOsa4tLS0a9euhBJBEOLxuCRJ2dnZ8Xh8xYoV48eP/9WvfsW62bquT5kyhVJaV1cXi8UkSVIUpbq6WtO02tra3Nxc9m3YtCtJ811zPAcAbrdb0zSXy3XHHXeIoijLMosKtOkI1NbW6rreo0cPAIjFYn6/nxBiGIYoiuvWrfN4PHPnzn3rrbe+/fZbNgisqqokSYIg1NXV+f3+WCzmdDofe+yxGTNmxGKx8ePHn3POOXv37p0yZconn3xSVVX16aefTp8+vaysrEePHgUFBaqqut1uhBAFCrbXup1JjiSbhqnpmiRJGGNBEPbu3evxeJgMS5KUlFPYnDypqakY4dra2vT0dMMwAABjXF9f7/V6CwsLH3744RUrVrz66qtXXXXV1KlTEUKbN29WFMU0zauvvlpVVV3XA4FAQ0PD5MmT58yZc+DAAbfbrarq5ZdfHovF9u7d+/nnn5977rkYY1EUq6ur169ff/nllwMAoQQjTAlFGFFKEWr893Q/j05FciQZYeRwOEKhkMfjCYfDOTk5iRg9URQPb2a/vNMHpVSWZUpoqi/VNEyHw8FWZmdnHzp0qKioaOLEiQMHDpwyZcq11157zTXXLFiwIDMzs1+/fhUVFUOHDq2oqIhEIhkZGXPnzvV6vSNGjFBVVRCE6urqQ4cOlZSUBAKB0tLSJ598cvTo0bFYrEuXLixOOx6P87wdttDuJM1Ojsfjuq6zOFvm5aKU8pz9CjsiHM9xwDEXRjwef+aZZyZPnuz3+30+35o1a7Kyst58883PPvts9uzZhYWFpmnqur5s2bLvvvvu7bffFkUxGAyy1z1kyJDBgwcjhLZt2zZ+/PhXXnll8+bNM2fO1DSNw5xhGg0NDampqZqmESAAgAlG2G7N24WkuRMdDoff70cI+f1+JskYY9td2dFAGCGMCCGqqsqyLPDCG2+8UVVVddttt7E5EoIg/P3vf3/mmWf+8pe/DBo0iOf5SCQiimJNTc3w4cNZr8rr9e7bt4/1zysrK7/55pvbb7/9+eefv+qqq6ZMmTJhwoQXX3wRACRJkmU5Ho8jhDDCGNli3I4kR9LYlwEAkUiE6WSMMc/xdne6w8K6vgij4cOHv/TSS16vl+O5SCTy5z//edmyZevWrRsyZIjL5ZJlefv27X369Nm3b9+UKVMAIBQK6bqekpLidDrfeOONCy+8cO3ate++++60adMEQTjvvPNee+21UCjUp28fFrxpxxqcGlBSHjSllJgEYYQxDgaDw4YNKywslGWZzWs7fDJbsE8r1igr0zRZGGazl2INrrRGaFmpqakpKSkBgG7duqWkpDBNbj2LdcEwDA5zAMAUMjug/SUknSR5vBDieM40TEVTqqqqhgwZ4nA44vE4mwRr0wE5VugVm5Z8fNLT05nAp6enW2WSiW4z3cDzvD2efApIjiTruq7rOkLI4XBkZWWxueZOp5PF/STlFDZJoZkyPDHdiBDKysqyrjn+NAmrNrZpJ5IjZjzPl5WVMdOroqLCNM0ffvihvr6e5/kRI0Yk5RQ2HYdEn5wQwkaJT/cV2SSvd52WlnbDDTf89NNPFRUVqqqOHj1aluVrrrkmIcl2k3zaScorsJrWGGOKaJvCPOzPoJ1I5ijR3XffbRiGLMsDBgxgAfR33HEH+5P9/joNLV+l/XI7AkmTZEmSLr300kmTJkmSVF1dLUnStddem5ubC/abtrFpf5IzCsXQNO3AgQMTJ05k8yi2bNni8/nsCNvOTVu/H/tjaCeS2bsWRbF///4PPPBAVVXVn//8Z5/Px7JG2LmOOzGojZzu6+20JHmIqLq6esKECXv37r3++usBwO122wJsY3MKSKZOZtPTBwwY8Nvf/pYFHthtsI3NqeEE7eRm1i+bla4ois/ni8ViHMcFg8HMzEzTMK2TKCilrQkhsrGxaStJ8HixoGsWf08IicfjaWlp1tAuhJCu6ywBRWK9ra5tbJJI2yTZNEzDNERRbCaHmqaxggaJNawgGDs4SzbAEg/Ykmxj0x60WScHg0Gv10spraqq4jiOJetKyHYkEjFN0+v1WneJx+OUUjabInE6W5JtbJJIm33XLIlxNBr94YcfGhoasrOze/XqxbTxunXrvvvuO9M0CwoKunXrxsJC5s2bd+DAgZ49e1511VVdunRhyXQJIT8ryc1sbFvybWyOQ5v9T0ypBoNBluO6oqKCDTWFQqHVq1ffeeedDz300HvvvReNRgkhs2bNysnJefHFF3v37r169WpomuBOCWV5ghKZN0kT8Xic/Vc39GbntcelbWyORdskOeGyzsnJsU5yMgxj8+bNvXr1ikajdXV1F110UU1NTTgcRgj1798fACZNmvR///d/0FTsi+0iy3Ii8yZugh1fURSWpcCWXhub1tA2SWZhOoQQjLHL5bL+KR6P9+vXLzc31+fzybK8efPmioqKgwcPDho0iNUWMk2TECIIgq7rCCPDMILBYMts2AIv1NTUOBwOuzttY9N6TiTGS9M0RVGYwcxg1Umi0egFF1ygqirrOefn56ekpLDUuRUVFaIosiTpLO8Mz/Pr1q3bs2ePy+Wy1rwXRfHgwYNPPPEEM6eTcIs2Nr8A2iDJiRoCHMcxMXY4HCzxIgDIspyZmck2EwThggsuaGho8Hq9oihSSnNzc62pYQghHMcRQiilKSkpVs1cWlq6adMmWyHb2LSJNutkSinP841+r0Awsf6ss856/fXX+/Xrl5qaum7duqlTp+bm5GqatmHDht69e5eWll5yySXA8rM1mcqjR48eN26cJElslgVbWVxcXFVVBQC6rtul4WxsWskJzqBgOjMrOwtjzNKU9+zZs6CgYOnSpcXFxRdeeGH//v0bGhruuuuuL7744scff1QUJZF1gIExdrvdGGPTMCnQxICT0+ksLCwEAFVV7YR+NjatpA2RIS23ZC4rn8+nqqrD4cAIYw5HIhFJkgKBQGZmJqsJFgwGJUliRjJLmMpg6Th1Q8cYJ2K/wuHwJZdc8s2Gbziea5ax9ajLNjY20Cad3DIfqs/nY/+1lnFzu90AwGxmplQTIV9WyWRTKTieY0UeE7DRKTsbo41NmzilKWwTiex/VkSbpby3sbE5PidlJ9vY2HQQTnVaebsJsLFpD+x5/zY2nQFbkm1sOgO2JNvYdAZsSbax6QzYkmxj0xmwJdnGpjNgS7KNTWfAlmQbm86ALck2Np2BUx3jZWPzs1hD7u2gwFbSLpJ8/EqrLadGWDdOJOWzph846pY27YFdJfcMJZmS3JqPwJ7h9AvBmtK8rU2DrZNPgOTYyezRa5oGALW1tcQkiXxd1jnJrIJUy9y3x8+Ga1ffPZVQSt98883vvvvuWBuYhgkAW7dufeyxx9gyAIRCIUVRYrFYLBaLx+PxeJylNGcdK03TVFUlhKiqGolETMMMBoO6rpuGqTXB/sQSP1FCKTni22j/+z7jSYJOZnk/dF03TVMUxVRfKqGEGEQQBIwxSxsCAHV1dW6321r5zaZjsn37dk3TfvWrXzVbTwjRNA1jjAgKhUIlJSUcz1155ZX79+/nOO7gwYN+v1+W5VAodNddd82aNYvthTGeP3/+a6+91qNHD7/fv379+quuuurdd9+VZdnayufl5S1durR3796apiGEMMKUUIRRojmw2/HjkwS54niOUAIAsizH43FJlADAMAwAwAiz9H0IIY/Hw1byPM9aXPaeoEVmH7sy6+mltLTUmgKZUVtbm5KSwrI4xWKxvXv3RsIRAFi+fLmmafn5+cXFxampqQBw2223jRw50rrviBEjhg8fvnjxYqfTecEFF1x44YVbt2794osv2PaMm266qVu3boSQWCyWkpLS/nfZ2UiOhuR5Xtd1RVH+/e9/67o+cuRIt9ttGAbP84WFheFwOCsrq2/fvoZhsBb39CbNNA2zWcqhU0bCevxZDcPKA5ySizoCwzAkSbKmc2KkpaWx7Itr1qz517/+deDAAcM0nn32WYfDce8990qSJMtyTU0NS1c+adIkthdrxMPhMM/zTqdz1apV3bt379279+7duz1uTzAYtJ5CEAQA8Pl8mqYRIACACWZ5oGx+luRIMkLI6XS++uqr559/fm1t7aJFi26++eac7Jxt27etWLGiX79+P/zwQ3Fx8ZgxY067vk0Y8KZpmqaJEGpZRLb9Tk0JBQzW7mJZWVmXLl0URRF4wSQmqylPCElJSTktz0rTtPT09PT09Egk4pAdifWxeMztdmuaFolEvF7vgQMHXC5Xenp6165dDdOQZVnX9YyMjLq6Olb6S1EUWZbZbfbr16+urg4Ali1bdvfdd9fX16uqeu9990qStGfPHrZZfX19JBJJSUmxe9EnRtKs1u+++y49PX306NHRaPTtt9/eunXrVVdd9fnnn994443du3cvKytbt27dmDFjOI7D6LQJs2mYJjFZMn1BEJgS0HWd/bf9viHTMDVd43leEIRmJl+XLl2gyVkIAAIvvL307aVLl/7zn/9sp4s5PhjjaDS6efPmeDxubUruvPNO9tfLLrvM4XC8/fbbhJBf//rXmZmZgiAoirJjx45Ro0Z9+eWXHo/nnXfe2b1791133ZWXm6ebem1trSiKsVisoqJi2LBhGzduHDNmzCsvv8J6RoFAwOfzXXXVVZIkqaoajUZZUkebNpE0SQ4EAtnZ2cQkLpdrwIABW7ZsqaioKC8v79Gjh8PhyMjI2Lp16y233MJzfCv7S+3hseR4rry0PDs7e9++feFwODMzU9d1Vv5GEIR27fNXVlYOHjw4MzOzZXtRW1vr9XoFQVi9evUjjzxSW1s7YMCA9ruS4yOJkqIolFJJkqxVfpjbkuM40zR37drVvXv3aDQ6Y8aML774YsGCBb/61a8effTRlStXTpw4sXv37hjjRYsW3X333YQSVVVTU1Pdbnc0Gk1LS1u/fn1WVtamTZsK9xTu27cvNzc3GAw6HI7KyspYLOb1emVZTrRrNq0naZI8duzY//f//t/FF1+89butGzZsSE1NZfUWEUKhUMjlcomiaG3jmWfyWEdjvd/2cHR37dpV1/X3339/zZo19fX1giBIksRx3FEDUZKI2+2+4447rr/ueszhZsKcnp7+448//vWvf/3888/j8Xh2dnZNTU37XQnjWGO2HM9RSocMGXLPPfckBpkY1dXV6enpkiQtWLAAIXT99dcHAoEFCxb4/f6cnJycnJzdu3efe+656enpy5cvHzx4cH5+PhuGrKioCAaDGRkZjz766O9+97vZs2d7vd4FCxZs377d4/F07dq1vLw8LS1NlmVVVXVdZ/WJbNpE0kRFFMWHH374jTfeSEtL6927NwB06dLF6XRijD0eTywWo5RifNiBgTB68sknN27cOGjQIMMwmG5k7paKigrWKrfTkBWltLq62jTNrKwsVVXj8biiKDzPNztdwqJOLDudzqysrNraWk3TmDXYegKBgNvt1g1d5uVmTcZzzz33/PPPA0BOTo4oij/99FO3bt0uv/zyk7nHxGUnxmMppRMmTLjxxhsTNbqa3SmDtWher1fXdasdVF1dnZGRoapqZWVlbW1tQUFBRUXFrFmz9u/fv3z58m7duqWmps6bN2/EiBGmab722msffvghAFBCMcYZGRm5ubm6rvfp06d///41NTU8z7/88stFRUVffvnl/fffjxBauXLltGnT5s+fb5cQOjGSJiqRSCQtLe3mm28+cODA119/PXHixHg8jhCqra3Nzc0tLCzs1auX9cvAGF977bUXX3yx3+8PBoNMNwqC8P3335eWlgq8EFfiLT2oJw8zUwcMGFBSUoIQysjI8Hq9PM/H43FroTl0JKxzUVZWVlNTE4/H+/bt29bz+nw+n8/HxuGs1NbWPvjgg5deeukLL7zw2WefsYrTCKFHH330ZG4z4SG3Rr+mp6fn5uYCgKZpzEdg3ZjR6JZrgT/NjxAKBAK33nrrs88+u2nTJqfTqes6czWPGDFi0qRJS5YsWbt2bTQaHTp0qMfjCYVCAODxeIqKimpqagzDYM6I0tJSl8tFKe3SpcuOHTsQQvfcc09paenMmTMppe3UF+v0JO2RuVyuvXv3BgKBbdu2DR8+fODAgbquX3jhhR988MHw4cOLiorOPvvshEI2DMM0zf79+7OBDdM0McbMh5yamrpkyRKO59rJc8uiDmbMmHHvvfeyloWVpDqWu6tRqRFKKHnggQe2b9/u8/muvPLKm266qU3nraioyM3NZfFMVinKyMhoaGgYNmzY22+/vXbt2qeffnrnzp0DBgw499xzT+Iuj4mqqpIkUUrZGGHrd9R0DXTIzs6eMWPGr0b9asOGDTzPcxyXk5OzePHiRx55pKamZv78+VOnTnU4HOvXr7fum5WVJQjCp59++vHHHwcCgXHjxjU0NAAAQmjs2LFdu3adNm3aX//6V2YhM38kMQkA2ENQrSc5ksy+9kOHDmmaduWVVzIDieO4sWPHyrKclpbGcdyFF16YiOnhed5qCzEzNR6PszY+EokAAIswSRbWoD+EkciJrLcPTTb5sWZrMGlHGHHAOZ3OmpqajIyMnj17tqzhfhx4jmd9WiY/Vvd1PB5PTU2NRCKEkIKCgtGjR69YseKtt946+VtmHO5dE8rxnKZpkiSJokgJ1XVdkqRWFqkWBRFhpKrqFVdcUVZe5vF4mA389ddfT5gwQZIkt9v99ddf//TTT2PHji0rK/N6vR6Ph1IaDAZZ5bCBAwc+99xzfr9/9+7dTqdTkiSM8Y033rhixYqhQ4dGo1H2TNhTOo0DHGcoSRtPBoDRo0ezt8taVkVR3C73xIkTrQYnWLQf+6DZwAwlVBIlSZTCkTDTFcz1kpTLa3lSazFnjHHr9b/L5WILTK+2su5cwl7leT4RQM7OzmTJITsIJWw8bPLkyZMnTzYNE+EkhLux01FKEYcopez6EUIUaOL2ASDRE25oaOjWvRvG2NpxIIQYhmHqJkKINc2FhYUjR47EGP/1r3999dVXg8Hg3Llzt2zZsmXLllWrVt10001PPvlknz598vLyWAQYQmjw4MHsaF6vt6Kigj3Dzz77bPHixeecc84TTzxx44032lGZJ0zSWj4mk5qmBYNBURRDoZAsy7qhH0saT1dY/Ml8KJFIJBqNGobBeg0nAM/zx9I2GGGMMIc5Rvt1LBsFmMOyLIfD4YEDB0qSNHbs2LFjx0qSdM4558yfP/+5557r3bt3ZhPZ2dldu3adP3++LMsNDQ1ut7t///51dXVLly4FgI8//njcuHGxWOyjjz7q2bPnXXfdde+9995yyy0fffQRAIiiuHv37ng8DgCxWCwUCtXU1Ph8vgcffHDcuHELFixwOBxfffXVK6+8csMNNxQVFZmGaSvkEyAJOplZNYQSSqksy06HU1EUj8cTDAa9Xq91dtuxQAglmhRKaSv7eyfDickz81p5vV5rwHBbT4q4o831w803a1fYKXRddzqcu3fvRggxwzU1NTUej1NKJVFCGJmmmXgXkiTV1taybQCAxWlNnjyZ9aUnTZp01llnsQ1cLtfkyZNHjx69c+fOYDDIcZyqqgUFBaqqmqbp8XhWrVp1xx135OTk3HbbbYMGDQKAtLS0jz/++N133xVF0SQmhzmEUSv7OzaMJA+ixmIxQgjrHrMIPmYDH/P0R3Z0ASAYDP7mN79Zs2ZNYmWy3uKx7rT1x7///vvXrVvndDqfe+650aNHN9u9TckVjrNBe3zBx5ouGovFmH8bY8xxHLOlCSENDQ1+v9+6MZvb4HK5FEVhQSOGYTidzlgsJssyxri2tjY9PV3XdWYuMaemLMtsqhzHcay7HgqFPB5PIBBwOp2iKMbjcSbqoigyQyMxu8aW5DaRZHc/iwRKwF5qs22avZWfnbN2rK+8rTQbazmBo+m6jhDSNO2w86x10yGshnHL/57AlbSVZsdPnLrZ+2JdA4xxMzFmK5mNLcsyAIiiyHyWiSOkp6eDxX2QCIZtViKbzbJKVN52OBzsaImzUHTEB2PLcCtpd4OkM70JNosr4fROrDyNl3RiJLxuHZCOfG0dmVMxBN9KlcUghDTrb1tJVn/7tHwrLTsjp/4azhTsh9NWbCfhiWNrD5uOgy3JJ4gtxjYdCluSbWw6A7YknyC2QrbpUNiSbGPTGbAl2camM2BLso1NZ8CWZBubzoAtyTY2nQFbkm1sOgMdImGSXcLLxuYkOZ2SbAuwjU2yaEdJPmplI5ZCuVlZJqtIU0rtPKk2Nm0lmXYyy7euKMrBgwfZAitioKpqopaXpmuafkSFgUS+uETVXFarlWVXtLGxaQ1J08mxWKy8vHzWrFnDhg3r3r17fn6+qqqhUGjJkiWCIMiyfMUVV+Tl5bEgx2YpgVgOWrbMcRwhpL6+nhCSn5+frMuzsencJE2SDcN49913Fy5cqKrq888/P2zYsLy8vHXr1vXv379Pnz6hUGjDhg1TpkxhyehZ8g3Wi26WJ52lbs3Jzklubk0bm85N0nrXdXV1ifLZ3bt3j8VihmG8+eabkydPHjBgwIgRIwoLC1mWH9SiAhPL5sfgOE7X9fqG+rq6OmuSaptTA6XUNEzTOEoC8KTDskoc9S1bi8u1FUqpruu6rifqWiXSV3RikqaTu+Z3TU1NnTVrVrdu3TRN69u3bygUcjgclFJWdLOoqEgURGjyhBFC4vE4S3DPsmSz49TU1MiynJGRkXgNpzLflQ1CyCRH8UomF0VRWA2wRB2MZm/ZNM1oNOrxeE4g43fCggMMiPxSPpikSbKqqbW1tddcc01dXd2yZcuKiopcLhdLxZiSkuJ0OjmOYzUQWL9aUZSFCxcuW7YsFosl8rMBQF5e3vbt2wGA4zkmzHZJkVPMKSjLxNL6WVVlsxSloiiybcDiVWl9O25Nu/8LIWnvbMOGDWefffaIESNYevcvvvhi6tSpBw8erKur69Kli67rKSkpCTEGAKfDed999z3wwAPNGt1gMDht2jQAUBRF4DvW+2B2QTsVrOo4sEKZCUFKIqxwD1hScLbs9DK5Zb02VVWZOca6bBzH8Tx//OfPsqxDUwWscDgMAL+E0upJk2SO46qqqljS40OHDvXt29flcl1yySVlZWV+v//f//732WefnRBjlhKd53iMsaqqVlc2pTQSibA6j8m6tiTCsnkn0eg6+SzcSYRSqmlafX09AKSlpSVqZbJyPydQ3cY0TOa5ZIUaMcbFxcWSJKWkpMiyLPCCtegf6wskxLuysrKoqCg/Pz8jI4PVM2B12I8vz2VlZUyfs3pUbGSUlcJKbGP93jrmZ3YCJP+yzcUAACAASURBVE29XHzxxdFo9G9/+9uSJUv69u170UUXAcCFF164bt26F198sbi4eOzYsUyMWbVEnuNZkURJkmRZlppIZJNmFVUO11vuAHmzEqm5T16SGytAWoozN/P9UAsnea7Wwx7yQw89xArBHL4YQgklrMw1pfSbb77Ztm3bt99+u2nTph9++OGf//znypUrv//++6+++mrHjh3bt2/fs2dPYWGhpmmGaSRuh5XgWbhw4Ycffujz+WRZvmLyFbfccgvLm8+kKxgMEpOwtPh5eXkPPvjg3r17WVvv9XonT5584MABTdMwxvF43DTMSCQSCoUIIXPmzHnhhRcA4IorrqCU/utf//qv//ovSumW77YsWbIEIaQoSiwWAwBWG7TZcz5lT7j9SKZF9Jvf/Mbr9bIuDc/xANC7d++H/vhQKBwihPi8PvbCoKnQ6XHgOK5jFgcyTTMpkpyAqTvAltE4DHD6dEXCwmwcI2StKs9xwBFCFEWJRqM33HDDBRdcYJpmRkbGpk2bsrKy+vTp8+233wYCgVGjRjFX1oYNGxYtWjRx4sQdO3YMHz6cZbrfv3//4sWLf/zxR3aKxYsXT5o0acuWLaNGjeI5PhAI+Hy+WCwmCEI0Gq2urpYkafTo0W63W1VVVVXvvvvu11577emnnwYAh8MRiUScTidrfe6///5+/frde++9NTU1hJCHHnpo5syZPMfv2buna9euuq6zDPusAoZpmOyZQyfSycmRZNaI+tP8LMCDObo4womiyPGcz+c7geeFOdz+jSVzmLe2ycAYt9pOtgaodcQmqU2oqqrrOsZYlmUmPPPmzWNFVf/4xz/efvvtZ5111ksvvSSK4u+n/17VVIxxbm5u9+7dQ6HQgAEDKKWqqlZXVz/44IMPP/yw0+GMxWK6rmdlZU2fPv2FF1548803HbIjcbq6ujqv1/vZZ5+NGTNm6dKlFRUVqqrG4/Gzzjpr5cqVgwYNcrvdwWAwHo9PmjSpZ8+erJzAtm3bKisrHQ5HMBicM2dOVlZWOBLetm2bpmnvvPNOZWWlIAh5eXlXXnml2+UGAGYvdBpJTs5HxvO8KIgcz9XV1TET1+FwCIKgqirrjuq6zl5eUk53Mpxk+Thmrf3sVk2/5iQ6zGfWByQIgkN2CIJgGEYsFguHw36/X9M0wzB27dolSVIkElFVtaKiguM5QRA4jsMYB4NBVieZ+UdefPFFl8v14IMPYg5HIhFXigsArrzySlVVH330UVVTfT5fIBCQRMnv9xcXF7/88suPP/44q6KenZ2dm5sbCARuvfXWmpoaZmxHIhFRFFnh9Ztuuun666+/+uqrKaULFy7ct2/fq6++umrVqvLycqa9Dxw48M0336iqego886eFpNVPxhw2DTMjPQMscycaayNTyiqJMr+ltY7E8bVu0j93VpQYAHRDFXgJAFNqArTW58ScugghURQTw91H29O0iDEGAI4TAJq7WFgAgyiKmMehUMjtdn/11Vf33Xffjh07rO6ZxMZtutOj7tIaXwOzQisrKz///HO3280GCAOBwNixYzHBJjUlUQoGgxhj5hD2+Xzdu3d3Op1ut1uWZWaL8jzPSrpxHNfQ0JCamjp79uwdO3b84x//AICffvopPz+fmdC5ublvv/32tGnTZs6cOXPmTOYc4RH/3HPPlZWVUUrHjh2bkZHRu3dva6mqH3744cCBAzfeeKPD4airqxs/fvzFF19cWlr69NNPp6amFhYWnn/++Xfeeec1U64ZNGjQY489pijKihUrNm3adPXVV3fWoYfkSHJLt8FRVd9pn+SUiHYIhUI8x3u9noQLs+XGjdqTUEIblx2OFI4TDIMEA2GMj3EvrEDZ4TJlCCgGoIRQRVFkWWaypCgKx3HMp+B2u2Ox2KxZs5YvX87KF7aU5BNo1BJtZdNdEEEQ6urqCCE+n+9YI64Y45KSkr1791ZWVlZVVTH3taqqF114EQBwmAOA3Nxc5pQuKSkpKytjnjBd1xVFkUQJAKKxqNvtDofDrLzjuHHjCCFLliy54YYbMMYej6d37947duzgeT4rK2vfvn333nvviy++OHXq1Llz5/r9/j/96U/p6ekZGRm6rgeDwdLS0oULF/79739nF4Aw2rx58zfffHPWWWchhJiEf/bZZ3/605/mzp375ptvzp49+29/+9v777/vTHF+9tlnACCKYk1NTXZ2dlpqGot7gSZnTbLqE512ktnTOINCONavX19eVi5JidtPiGVjg40OwyHEYYQRQtVVtZLoUBV947ebysoqjnFsAsgAZAAQAAwUs2PyvHjOOecMHjyYbcRGaxVF4Xn+scceY18wpZSZeWxE1MrJ62SMsK7rLpeLjQYda69QKJSdnV1QUHDXXXfpuq6qKlhqKTIkSTJMo76h/sMPP7z66qtZlC7P8263mzlKmo3fnnvuuQ8//DAA3H///d26dZszZ44oijNmzMjJyampqXnvvff69ev3zTff3HPPPWzizdq1a997771Vq1a53W6nwzl8+PDHH3+cObcN0zA1c/Xq1ZdddllmZibrGW3dunX+/PkLFiwYMmTIhx9+mJOT8/rrr99xxx2PPvqopmlffvllQUHBhg0bZsyYwfGcqR3uTHUCAU7QOW2GY0IBKJi6tvKzlRs2rnemODAFRAGzyZQUAwDQxq+cNA4OATDvCOIEQXCKYiwW/2b11xub/PBHQhBSMDQOvZgIE+B1jCkCTdGfeGz24MGDmcvA6/WGw+Hvvvvuvvvuo5R27969urqa53me56+66io2ANt4ySc0FoUxZtYpG+NhyzfffPPNN98sSRLTqIefiuXgfr+f4zh2dkEQmEHEJDlhUNTV1ZmmGQgEli5dumLFCkKIrusHDhwoLy9nW/bp04f1NVwuV0NDw5///GdREE1iTpgwwTCM7777bsmSJazBqqioKCkpCYfDCKFHHnmEDWgvW7YsNTW1urqaEEIo6dOnjyzLRUVF/fv1RxgRQnbv3j1nzhxW/FkEceDAgQ8//LAoik8++aTf7y8rK9u6deu0adN69+49ZsyYVatWnX/++aWlpWPGjGE+cFY+Fo4xhf4MJWl2MrSoD9wRoQAEOE6UBJEXRRBQelpq3cEKMKlHdiGKgfKWjjbrfREAAA6AYiAAiubBHMRVsOjxwzsghTOCsVBDempOIBJ3eFMjoNVq8bSsDKrqbrfbNEyMMFO5brd7/Pjxu3fvfuqpp1555RWv16tpWigU2rJlCwtZ53me9YGtg8+N99FqwWZbslAKpo1NwySUIIKss1kSB4zH42yBlThnksn0YaLP5fF4/vWvf73wwgsvv/xyenq6YRgIoTlz5qxZs+b7778HgLy8PK/XW15ePmTIkNTUVE3TTGKysNzCwsLu3buz8eF4PJ6bm+twOJh1nZeXx46vaZrT6YzH4yyCCACuvPLKjz/+ePbs2aZh7tu3z+Vy5ebmJu7xP//5z+rVqzVN27lz53XXXff+++8vXbr0/vvvP++882bMmNGlSxfTNPv168cc4D6fzzRMTdccDodhGLYkn7FwACaYCAOAx5EyfMAA/+BhqZIzWhfiCAZqeSAUAwCFY3u2WoBAyfZLsVAw3ZcbCBsBw1Cc4hf//i6OwNR0Xdc1XeM4jrM0AuFw+KGHHrriiivmz5//3nvvMeM5Ly+PGYTH+s7aKsnQNHuU9X6bjdU3jmkDAIDT6YxGo2yo9lhtRyQSOf/885999tmzzz47sQ1CaOLEiQUFBQAgSVJ1dXVubi5GOBaLsbFc1igsWrSooKCAxfYBQCgU2rFjR7OJ6Mz07datG/OYmoY5dOjQZ555Zvbs2RzPrVy5curUqdbthwwZ0q9fvyeeeOKuu+6aMmWKYRjvvffepZde2tDQIEnSpEmTPvzww507d6qq6vP5Dh065PV6HQ4HG75u5WPs+CTVTu6wqjgBBgDQedA5AACZ0GF53Yfk5GdLKQ7KCyY+vBEA8z9TSpuPJ6FjjWMRgg2EzWgk4vfmBiKqJsmHDHXjli2BhroUR4osywmPVwKnw6lq6rBhwxYtWnT55ZevXLmyuro6PT09MUX0qLT+UR92lSNKKcXksAw308atx+12cxx39tlng2V6uW7quq6zMZ5QKJSWlpaWlsbxXENVAwvM2Llzp8/n27t37/333w8A8Xic4zhZln0+XzO/QCAQoJTW1dUxv3c4Eh49erTD4QCA/fv3f/DBBytWrLBu73Q6HQ5HUVERADQ0NGzfvv3WW2/NzMzkeT4UClVVVWVmZm7YsGHs2LEAkJqaKkkSE+N4PM4O2wn4ZelkNuKkA1BEEACJKi6CMyifEtMciHAEjhxgJ3BYJ5Nm65uDCAAhCEwMXNzgJE1UdFFwSnETwooDeE3TmY/XGjNMCDFMwzRNNhw1ceLEkSNHdunSBQBMw0w4VymhzOOaiEv5WfFLSCmLZrGGGTOFbD2ONWeLGlNTUlJcLhfrXTc/LKEIo3A4/OabbzocjsrKypSUFPYnNsqoKAoApKSkNDQ0HDx4sEePHn6/XxTFQ4cOGYYxYcKE2bNnd+3aNRaLuV1uhNH69esDgQCb/Zo4iyiKTqdT1/XU1FQAYCr96aefVlU1HA6PGDEiJyfHelWRSMTlcr3zzjsVFRXXX3/9oUOHZs2aFY/Hq6qqZs6cmZ+f/95775177rnz5s0rKCjYvXv34MGDmRh3prHlpN3JGaCQLSAKPCG8Zjp04tCJBBQ4ah7lDgilTTEeqDXCDNjEPmcKxFUnJxKEU1NSeJO6RbleVxsn1jfFCTIBY2GMbGwZY5yXl8fcXc1G7NjwD8KIuaASAnm8e8QIANj2bHdCCUbYOsRQXl7+wAMPyLJcWVnZp0+fiooKQRD27t27a9euDRs2NEuldu11106ZMqW4qDg3NxchFA6HU1JS2PGhyXvEvNaGYQwfPjwajQIASx3xyiuvfPTRR2+88ca5554rSdLmzZuLfypWNXX9+vUTJkxoFgDLvANsx2nTprFYo9TU1AULFuzcuRMApk6dKklSVVXVwIEDn332WZfL9T//8z8NDQ3vvPPO3XffPWbMmA8++GDx4sWLFi166aWXCgoKCCEfffTRnx/98xNPPLHmqzW6oauqijHG6Ihm7oym87RJrQEBYAABgKPAUeCBchQ4SjigiB75OpvkljQ5opviR34mRKzpeyQUEZ1okuxU1Th2eDHRGz3Jx4g5Z59vy1l4CCHWMWYfXGLCUGtv2bKl9R6Z1yczM/MPf/iD2+2ur6/Pycmpr693OBzMgk0429j2Xq+X2b2GYVx33XWXXXZZdnZ2YrIU29IwDNa+GKbhcDguueQSdsFbtmxxOBxvvPHG2LFjQ6EQx3GjRo364osvgsHgyJEjr7/+etM0rUmg6uvrMzIyBg0apKrqjBkzKisr09LSdF1PXBtLUVBdXT106FBZlvfs2bNmzZoxY8asXLnS7XZjjM877zw2Pt+rV6+ampq8vLyBAwfOXzB/9uzZJaUlLJ9cY0ICk7CYzTNdnpM5GSApVFdX/+53v1u+fPlxJouf+EOnAAiAwN13/37L+rU5An/f9Teek9fNTSnG6OixXvRnRNdyWYRiomNMASPCm4iP8xKXmT75tqmu/LwGRXlu7l/HjBljjfE6fBJK2TQj1qFlGiPxBFoGeLb+CTQ7S9M9UZOYhJDEjH8WtaLrOs/xbIAqMZc4saMgCKxHynxXzNOe6ktl22AOW68qFovV19e73W4mWolxL3YiSikbhaqsrGReaGvMHMdztbW1lZWVAwcMZHtpmsY65Oy/7GhWK7euro75yQoLC/v17cc2Y71u6wK7EZZGDgBYD6WtaQw6JkmL8YKm7uJR/5TgqM/rVD9EClzjGDJQBAQRoAhoS6FlneCW+x+jW0sBEFAAAwMPhAAgCg5RiuqqhIH+XMMvCAIllCVjOH7wxsmDMOLxYXPdNE1JlDRNYzHSLDz+qNeQmFTE5vGnp6dDkxCyIWhmFLDRHZZhor6+3uPxqKrqdrs1TWOtQCQSwRgbhpGbmxsOh2VZthoL4XDYMIzBgwdrmsZC0wzD8Pl8iqKEwiEAcDgcCCOHw8FGrdmVs8MOGDDAMIyG2oZ4PO71eoPBYFVVVd++fZnYI4SOF2Z7JpO03rXVAmzG4UFLQilQq53WmujrZMLOjA+LJ6KAKQagkLCHAaDxqyIACUnGh4+ALKuPPLiJgaDGNkIgIGO+uqQ8L79Lg6mTphCPZnPcE9GCrBfNPFstRahZE3BiT8zaHUicglLKcicwcx0AeJ4/VouT8IFZ+9UslwCGw3PXeJ6nhLLBc6a0U1JSCCGCILAjJ8wHSqlVWzJY3igASGzPmg/m/Ldej9vtZtfPYc40TIfsYEH+aWlp1m1Y76PxUWOEaKOX4QyKSvxZkhd3TSihJDHI0WwstPHLwEAIOSIt7umYI0oQNMvThuhhEQU4Mgr7aJdGofEIFsOTACKJwyIKmBIOYYEAoqCbJgtyICY5VqY7JudMijpsjI1VgE+MNnmY2sl87ZhT30+S9vJ4maaZyMxg7XsjhLimuQencXYobSHMAIfjNFuORbUYUj66hAMABgIUADCiAIQAMSklQE5Vp6N1dLQG4uRpk3ZFGLHhNGh6FJ3ggSRJJ7f4UjHCuqGLomjVzMxMYqn2OljHBjdGcv4sCGirLpwAIgBmU8x0x5Jkm85HkuKuMWIp1xJtWyQSSWSxxhgLvBCLx5glZp1Hdoa1hS0u9kjjFTCAiRIjWASoSdnvqDuzIxzjCZyYp/oEQE3Jyawr2zrX7yhbtujAWq30Zs6Ck6GtMf+N99uJsoUwkmMwYIxFUSwrK6uqqmKeT7YGYyyJUnV1NQCIoigIQjgcdjgcJjGJSU5b1xrafUoqpgDocNe6ow312XQ+kmYn7969e+PGjYSQWCyGMR46dOi4ceOCweDixYv37NkzZMgQlq+vlUdrjaY6Gawze5JLU4oB0uQSNxt7IZ1KAfw8p7iZbtPpOpk2ZiRNkgcNGpSRkZGZmVlbW7tw4UImtMuWLXO5XM8+++z27duXLFly5513ut1uFl3ws8VKiEl+NgXnCUCarPr2d2CSRKAYogQdZRKkjU3SSNrXrOt6Zmamrut79uwxDGPQoEG6rm/fvv3yyy93uVxnnXVWSUkJm/LO4l0VRWFpsRiJ+fTRaJRFFyTKLCMLJ3+d1uwACGGcODLCh3+HU+qRRjd1C2c1bootsa6z3IhJKQU2SIwo6wIgxLGOALVg3T+5d9p60DE4meOczF/beg3H2vf4xzxdT7udSJpOZvHAgiCUlJSce+65siyHw+H6+vrU1FRFUZwOZ3V1NSvmxoIQZVn++uuvf/rpp9zcXKtIV1dXV1RUAICqqu1R0IRSCu1nJ7MedWNZMcKCxvAvr2ttc+pJmiRHo9FoNKooSlVV1eTJk+HIwN2W/PTTT/v27SsuLg6Hw1ZJZhPZ2q/rSylBqH0rWmCLqxq3ZmSrM9KRFV1HvrYTJpmRIV6vd82aNcOGDWM5zUVR9Hg8jZPFw+HMzMzElqZpZmdn33rrrcy/TSlNTKArPVi6fPlywzRYgF7SOUInJ9uljA6bKy0nQtrYtCNJU30sLduuXbtGjx7NUqXFYrFRo0YVFxeHQqHCwkKWYqLxrAiLgijwAovOtc6DZSnRoUU+x44AbREWwqzlxO/YdMLwQJsORTJnUGzbtm3UqFGyLLPEEX6/f+LEia+99lp+fn4sFps0aVJiYxYud9TjRCKRlJQUSZISqWROARQBQk2zIhILJ8vRy1DY2LQHSRMVURQHDRrEutCJSbZ5eXnPPPMMALC0SSdZyaV9sUyfOBlBJgCAgGuK6yV2B9vmlJC0Xp8kSZmZmQ0NDQBAKU3UBAMAlh61lUWhHA5HPB5XVfU0lKU/9ryItoABgDKPHXNhI6AIiN3BtmlPkqaTmdAyzzMAIIRSUlLYeClL7MAynlp2SNaZkwOFw2J8RJ3F1itollUXm+zeCGLTMnBjJYpO6C616UAkX54Ss9VPbHfTNHmeZ+WXz1wIAgDcpJ+b1lJsmThpY5NMzmyBOQVYpzG3Uj8TBBxt3JagRv1sq2SbdsVWEe0CZTqZibH9kG3an+TrZHTkfFdrPM2ZEltj1b1HSS1yfBpNYnz4vzY27Y/9nbUf+Mh0QvajtmlHbDu5EQSHA0LstAA2Zxy2ovgZ2jAKZWNz+mgXndwae/hY2zROGG3Ke9iy/MIpAFmklyJbmG3OAGydbGPTGbDt5CZsxWtzJmPrZBubzoCtk3+GZjbzUdfb2Jx2bEluPwggcjiLvR2vadOe2JLchEXQTkrdImLJMZBYtq0Ym/YlaV9Ys7K0R61S28pMA9aaoGcspEmkLRxW0TY2SSaZAkMIqaqqAgBVVU1ihsNhtr6hoUFV1XA4fNSMP5RSVkHKNEzTMGOxmGEYp75cNW1K08XKOCZ+VhA9/GuxMwAluGVhR0RsAbY5BSStd83xXDwez8rK2rVrV11d3UUXXcTzvKqqhYWFPM/rut63b9+WdSdYJl1WW5il1BVFUdf1QDAAAH6/P1mXdwppyqrZfO6ELcw27Ugy7eS6urqXX375rLPOwgizAqvl5eWrVq3Kzc1FCAUCgXHjxjXbBWOMMbYm+hEEwef1+by+JF5YKzmsgY82F6oNkV7NNTD7r11NxqYdSZokm4b56quvPvXUU/F4XJIk0zRDodC77757xRVXZGRkVFdX//Of/0xIciIAk2W6JpQQQliykfLy8vqGesM0JElK1rW1L/ZwlE0HIGl2clV11cCBA9esWfOPf/wjFAoZhiGKYnFxcf/+/bOzsvv06VNWVpbYmKXsAwCEkGEapmnyHC/LsizL+fn5Xbt2RQix5H4dhKMoZNr0syw32s/N+9X2lEabdidpX9iqVav27t2bm5s7adKkl156ac2aNWVlZfn5+bIsR2NRjHF9fT00ua9ZjzoWizEFLssyx3O6rkcikV27djU0NIii2PoSrTY2NknrXQ8ZMsTlcvXv35/n+enTp3/66acjR44sLS0FAI/HAwD5+fmUUtM0McaEkKKioi+//PL7778XRTESiTB5FkUxLS1t3759qqpKknQGFxCnv6CcIYnXdKbkhOmUJE2SA4EAIQQhpOt6UVGRw+HIz893Op07duwYOnRocXExS2fP3joltE+fPocOHerevfvgwYMPHjxICDEMgxASj8d79OghiuKhQ4eysrKSdXmnENzo4qL48HLj+k7IGdzadi6SJsnjx4+fP3/+hg0bMjMzf/zxx5EjRwLAlClTtmzZEgwGCwsLr7vuOqaT2fa6ro8ZMwYAKKXdu3dPfBD19fWvvfYaQogNQTHf2CmboszsYWJZbmx62J9ZEcZjfrqYNG6LAQDhppKQiEOAMcKJApHsjtrxNtqfSCSiqqrb7baOO1hvytbPp5ikSXIwGLz55ptrampqampuvPFGURRLS0vHjRvncrn8fr8gCLm5udD0sgn9RQ6uUgy0M+TLjUQilFK3282q8J3prVLnIDmSTCl1pbgQRoIgpKenu91uXdd9Ph8AjBo1KhKJdOnShdm9CZ3cmgrJ1k+ELZ85LT1usdB5YNpYFMXGDottJ3cAkvOdIYQ4ngsGg06n0+fzMX+Vw+FgQ03MQa3r+qnP43MM2lO6jvB1ndliTFvA1ie0cXV1tbViLoO9YtJEYv1xInBbnsjW820laZ8apdTn87F3IEkSIYTjOEEQCCF+v59SyurLsHgPjuNYpi5oevHoSMCitK3vNQnvuFHS2qYzjxJrfTwwACCEgWUkQwhj3LjYgUXbKkVMCFs+7cRfAYAQku5PZ++RxeoRk7CCuwBACW38UQoAsVjMJCYAKIoSCoWCwSBbX1NTw6KDjnWuU3LrnYFT+mUhhJg8J8T1FEPAOjsCt5DqpIAt2jix0GJmRccj0YxWVFRUVlZijBOi2ExPGoYBACUlJZquMRkuLS2trq6OxWOsmh8xiUlMQgmhxDCMQCDA87woirW1tSz+x+v1BgIBAHC73aFQiAn5Uc9l00pOj45opSQ309JJgVjKIyOECEFHzqAkLX4WWl+irWXXuuUkxw5JbW3tqlWrVq1aVVtbe6xtXC4XAMyfP//ll19WVbWhoSE1NXXmzJmEELZXVXWV9ZXV1tauXbsWAP7zn//Mnz8fAF544QWfz1dSUnLppZc6HA5FUY46Tw5szdxqOlxvL/Hm2s+/TRCmgClgijAhBCGOnRE3DTKxPx9xVUfkIcDH+AFBWEcUCzwxKEIcADYbdySUmoBIa6K0mT2Z9M830TFupvSaTRqXZbmysrKyspKN/1thwQIAUF1dHY/HA4FAz549VVX1+Xxut3v48OF33nlneno6AGRlZiVOCgCapi1evDgWi+3Zs4cQsnbt2h07dhQXF7MRDUppijPliCdgmqZp2jLcJpImyei4WLc5/r6J99fsLSZLP1PAiQ42xcgglMMcJQhRQKTRGEaNGvUIvcqmLpuY/fARP8Szn4FBBUIFwTAIBhEcTsEpSw5RdohNCrn5xTNBYl1NXddVVSWUxGKx43/ER3URHZXme5HD6ymlwWAwFoux0xFCQqGQy+XieZ7neZfLhTDieI4SyrrQlFLTNIlJMjMzFUXZuHHj1Vdf7fF42BuZOnVqdXW1pmmKonB8o4sEAARBGDhwYM+ePXft2qVpWiwWW/j6wieffLJHjx6LFi0655xzCCEcz3E8Z5gGi+qlTeFDkBjPt0X65+hw2X+OJcnJOHSTHNFGvWsipFNCMGcC4iiyTCo+Qk0d/g8ijUUYgRxR86lJPE2KMIepSXiQgKJAfR1N85SUHeQzM8jRbofNA8MYu11uaCo6jRCqq6tzOp1JvflGMSaUgAkAIAhCKBTyer21tbVjx45VVZVSGg6HdV1PS0sDgL/85S89e/asra01DENRlDlz5vzud78TdJ1uegAAIABJREFUBbGuvs6f5v/b3/720EMPLVy48PHHH3c4HPF4HGNsmmaXLl1cLpckSZFI5Iknnvjdf/2upqbmo48+ys3NXbt27b59+0Kh0MizR77yyisFBQVffvllIBB4+eWX09LSdu/e3b9//yeeeOKSSy5p9KgBwQQnPKM2x6fDSbKV9vCKIQSIJQPhMEGgYKpJmMSxgBGiPLNmaeO0ZAIAFKi1a00AABlHFmQ8bDzzGGGKzXDcJbtNw3C63LzX7XS5BI8rEoi2HELneA5THI/HRUFUFGXVqlVPPfXU1q1bs7Ozk37jzTAN0+Px1NTUpKSkHDx4sLKy0ul0UkojkcjChQsBYPr06W63GwAikciLL75oGIYgCAihzMzMDRs2fPDBB++//37Pnj3/8Ic/UEoRQqZhGqYhiiJ7awcOHGB9b6/XO3z4cFmWP/nkky1btlx33XUzZ84EgJtuuikYDO7evds0zU2bNj322GNrvloDAKFwyNqKUUIB2yPVP0+HluSWavmk3igCsERCEw5riGgcVngMIiaYQ0AoIgBAKQEgLMdAo9FutZMPJyRo6n7TxnLnmAIvcJG4Lqc4g5GoahqHDpYjxNVUVTukI0xB1DQ9OxwOy7K8/4f9L7zwwurVqx0ORygUYnNO2g8WYqkoitfrjUajfr+fuf0QQk6Hk1nIToeT+ZxdLldaWlo8HkcIEUIKCwuXLl1KCFFVNTU1denSpYFAYPr06XElzkR969atzz///NKlSxVFUVVVFMXzzz//66+/Xr9+/aRJkw4dOnTBBRc8/PDDGOOioqKampq0tDSv15uTk2MSkxByxsxL72CctrpQrYG2iOtquaatIAQIoD4YoCIXVZS3//HR6k8/TXU6o2qscWiq0ULD0FSbCgAAo4TuJQltbPFjEwQUEZMQSk2nLIfrgl63x+1MqQ0HNR0QAJa4hPTSppSDzIG0aNGiJ5980jRN1uNlwXAt79HarrXS9LAexDTNysrKtLQ0l8tlGibCiEmsKIqsNTENk9mrFRUVfr/fJCYL7wGA2trabt26sSbm9ddfz8vLk2XZ5/MFAoFrplyTm5d7yy23mKbpcrmqqqpqa2uZkcyCwCRJ+stf/rJ06dJXXnll586diqI899xzt99+ezgcvv322z/99NOLL754y5Yt3bp143k+FouxvRqfEocxh22F3Bo6tk5mIxOWPunJv1REARBgjAkCzilTk9SGwgerK8AhmhgDxQRhdl6EOIKgKUyFI+hwVIll2XK1mCBEBYETTcPEEFXih378kRN4T4Y/M81X3dDoyCEmYfnMwuFwTU3NJZdcUl1dzXFcLBbr1q2bqqojR45UVZU5gY84ftslGQBqa2uj0ajT6czIyLjnnnsmTJjgcrmisSil1Ol0MuXsdrtNw1Q11ck7AUBRFEVRBEEQRVFVVYRQRkYG064AMGLEiIsvvvj1118/dOjQwIEDAcDn823ZsuWiiy7SNC0rK2v58uVjx45lZ2cK9rzzzps+fXpZWdlXq79aunQpx3N//OMfu3fvvmnTpv/93/+94447du/ePWjQIE3TBEFIzIHD2BbjNtDhJDnhnWZBUck9eGJ4CSOEESouLu6VleVwOZGuI0lACBPAOCEnmOcS3WBm41IMABSxlFzN1TIAICBUVRsagj6XT6d6Vl5eTFFU3QiUV8gpKQhxAEAo4Sw5vXbt2rVw4cJ58+Y1NDSUl5ebprl582aXy9UyZYp1uKhN7sCUlJRYLFZTU5OVlcXC71jvvaqqiuO49PR0TdP6D+hfUlKSm5sbjUZN0/T7/XPmzOE4jhCSmZnJ9CcLCLnppptEUXQ6nSkpKaFQSBCEm2++ee3atWPGjGFXuHnz5oceesh6AU888UTv3r23bt2qadr4i8dnZmb6/f558+YVFBTceOONc+fO/eKLL1566aVgMOj1ert27cqCfNnssdbf5i+cDifJpwAEICDAxBw+dPClBePOHzpcwnxMUQEwAQwAlJoAwAQPIUzREX0Bi818RDoBTAEDcTqk+rq67NwugUAoPSvzx6Kf3nrvndKDB3XLyBDrx7rdbkmSBEG46aabzj///FdffXXz5s2iKLJISY+7ualsHWBvvSSzid8pKSlslqimacyHRAjJysqKxWJ79+7Vdb20tBRjHAgEfD7fM888k5aWNn3adEKJIAiqqv7pT3/Kz893uVxsAwAoLy9n/mqHw3HRRRc9/fTTACDL8vfffw8A3bt3P/y0EVqyZInH47nsssu2bdsGANOmTRsyZEgkEsEYP/zww9dcc819993HrhMAAoGALMscx2GuwwU7dGTOAElOVsNMoNHphQiIiMM6CdTW98zv1rt3bxnxoBNMMGJat1FOGiWHtOX8FOEeOV0ogux0v0nRqBFnffjxh7qhYp5HlDRltG/UyRzmdF13u9xDhwxdtHDR6q9WL1q4iA0XtfyOcWM490k9DUEQrL4Gp9PZv39/lqUcAHw+HwuYj0ajhBKe5ymloijOnTsXAEzDdLvcpmHWN9S73W5N09ismB49euzcuVPTNFmWi4qKfvvb3yYaBU3TOMylpaVddtllo0eP1nV99+7dGzdunDdvniRJ8Xg8MzMzFov16NGDZZJhU+hUVRUEQdd1XdeTPhrXWTkDJLk9EBDGJhUQYCDshwjhGkOwcFvSZVr0JCZAMUVgIDARYAAKCEFjOnuaCMsgFLjDE4OsQ1OTJk3q1auX1+sNhUJsBCjpWI2XEz6Ix+NpaGjo169fampqMBj0eDx33XVXdXV1Tk7Oo48+umzZMto0nYa505wO59dff/3OO+/8+te//uGHH+bNm1dRUZGXl1dUVLRmzZrbbrtt/vz5w4YNy83NDQQClNLU1FQAqKqq6tKlS9LuvLPT4ST5CDu5PcwkSjEgDigmBoeBI4AooCPcaq2Kqfw5EhFdBIBgILQV+a6DwWCfPn0IIS0jJTsUkiQZhlFeXr5x48a+ffsihIYNG1ZcXLx3716fz1dZWRmLxbxe78aNGy+55JK01LRYPBYOhw8cOKDr+q233vrYY48tWLDg22+/femll/77v//7umuvc7lcBQUFCxcuHDNmTCQSYcnSs7KyEp15m5+lHSWZENKsvNNRB1dOPYgCjzlEKJtliBOxX0k/UdPpAKApTpIcx8SNRqOuFJemayzY67RwrJkMVuLxeH5+vqqq33///SeffFJRUdG37/9n78vDo6rO/885d5m5k9myTBISQthBEDEICMriglIUqiiIdbcqKmpxa2uLuwVp/aqIv1pFUKECVorWFZFFERSwAiJlD0sChGSyzD5z79x7zvn98SaXIQsGmIQI83l4eG5m7tz1vOd9z7t83u66rh86dKhnz54vvPCCIAgIIbvd/pvf/EYQhXvvvXft2rVDhw5dunSp3W4vKip65JFH8vLy/v73vw8cOFDX9fHjxxcWFj7yyCPDhw8HS/7QoUP5+fktfrenEZI5YqhBKaMmIwylVNd1URAh6NJGxLgOEBxuYtQ2/zoT9sQYI4QRxgJGCGNW+y3GWCAI81oeL3IMnzykdgmC0BxClZbDhRdeCBv1U9/rLpsQMnny5AEDBuTm5nbu3NncAd41aFSEUDweh88fe+yxHj16SJIUj8cJIWPGjLn44otVVW2X244alGAiS/IFF1ywdu1aVKcAsrOzTdqZFJqDpLH/YIwFURCQwBiLRqMQjQQlBHHCtiTGCUgus0dtHxkzi/NEblkUxeYoxhYC4+zss89GCVyI9XaAaoqbbrrJYrEkijFCSBCFcDhst9vj8TilFLJKotGo2+0GK0MgAuccSpRr/fh1d2pOXiD/sL6Ao7XsDZ8uSI4kQ9ptMBSE7hOMscLCQoRQKBQqKyuLRqMOh6Ndu3Yt5MVpWTTVZvFYhconW495qsoGGGOSJIHDCSEEqWAN9zEMA2MsSZLP57PZbGbuWjwet9vtoVAIXnQ0GmWM2e120/+cmOlR6w0RcGLlE0qQYVRXCJ1Cc5DMXo2bNm3atWuX2+Vul9euoKBA1/XDhw+vWLFCVdXc3NyePXv2798fdm6OpQ1rrZbFUZxbJ10OzQmCtXHdnWGOuFn2fLJHbyXwhJZdjU4oGGOr1Qr9N9PT06lBIdAtEAETHI/HIWMMIWS1WiG2BO4S88hQPgnZnYkkuwDIYAXhpwZNpXk1E8lcJ2/ZsmXSpEmSJAUCAavVquv6ggULRo8eXVRUdPDgwXfeead///6appnNGXVdFwTBfM1wEMMwIpEIQogzjoWW8F2bPLVHr5MTdWxz2h0n7mP+FjOEEOGIM4Iwoqadzc15ok3T9EERcuInDaUIPhFFEci3BFHA7EhYC+LVmGCoam50ysYYwwIYhLb2Q3MRzghcCUKIMtrG3CttF0kbVZqmlZeXr1mzpqKiQlEUXdeB+WHAgAGRSKRTp05er5ca1GKxBAIB+IkkSYQQatB4PG7UweVyKYoSCARiaixZ19YIEquLk9jw5ehD1W/0dny0fqcexyaNEEQhMV6YCEIIuDnr/aTeRqMuekwwJhi4TUAhp2gGmoOk6WSLxXLZZZcFAoG3334bPCI1NTUHDhzQdd1isQSDwVAoJIhCKBTKysoCkwzsLkmSEvOQS0tLXS6XIAiqqqalpR3jjK0DThrXz5gdS/4xQoSjJllh2zAaylvzf3Iyp0uJ60kiOZLMGKOUQntkv9//7rvvtmvXThAE6PDEGFMUJTc3V9d14FK02+3lh8rfeuutDz/8sKampkOHDuahLBbL7t27CSGtkh3RIoYurgtQA7du2zWmUziNkKQoFOOaphmGIQiCxWKprKz0er3wVUVFRUZGhtfrzcrKgloZh8PBKMvNzX18yuOPT3kcHe2qNQzj9ttvB4LopFxb02g6nny8qLdmxojwprK1G1J2nrmS3pxXnFokNxPJkWRYFO3cudPr9cqy7Ha7bTabw+Ho27fv/Pnzr7jiijVr1owYMcJisVRVVWVlZRGBNGV6lpSUANOqLMlJubZTi1q13ICsM4UUkoukKYTKysod23d4vV6Px3PjjTdCifk111xTVFRUVlY2cuTIonOLOOdZWVlVVVW1nhKoEiZHOUvat2/v8XigN3oLz8dtgn26Ua9SShGlcLxImscrNzf3+t9cjxACHnNN06A71MUXX4wQCgQCgigA70RWVhYwUYikkbPv37//GJzpSQCsYrHZFALKDJuc0TA/sUytWgC37hFe3pSEptAySJokQ2oeQggyhOrxqrlcrsQPgbmiUc3j8XiOZCa0pmrCrPFwFCdHQkeJ6+GG8WTMECccUcQRRYQhxDninANXfF1jKIwI5gl+gabiUmeaWj7T7jfpOGXulnplUg2RmDbQ4mhOKkizj2Mq3tTSOIVWw6msTz7V0zBrCzmUp/ohpHCa4MwNgaSQwumEti7JLamyjv/eoVdjw46NmB2xz2Ej8f8UUmh5tDn2n1OJpAkeQ4ggzFKu6hRaDWeyJNdbJ7OEVjMJqG0K19aNlxTOcKQGaCJSxnAKv1S0aZ3cIiHlpuqTGwF0mWjqOHUx5OND/bbMKaSQFLRpSW5BtIK1zE+QxyuFFE4AZ6okN4J662RzmyT8nyTzO1HCU8KeQjJwhpt5rXb7idNBCikkH2eyToZ1cgIN0JH/UYPt5PD1IY5T8pxCSyBpoyrF3tI4ksunnUIKTSA5g+yXKcYs4f8UUvhlI5k6OVGeG5XtNiHwUEaIU2KcwmmFZK6Tzf4jQOiFMYauIqqq2mw2QgiQ0bf56p8GOdWosahVKqc6hbaEZEqyyVQOxPTUqO32JktyJBpxOBxtnrs4taZN4ZeKpPWF4pyrqoox9vl8VVVVZ599dlyPIx3t3LnTZrPV1NT06tUrzZYGNJptQy23ifrkFFJICpKmkymlgUAgOzv7008/lSTp7LPPNgyjpKTkk08+6dChA/SqP+ussxx2B6i9poQZzHLUYOGdQgopHANJk2RgpX///fdDoVAoFEII2e32jz/++De/+U3Xrl1//PHHH374oU+fPsCkiepa/hmGgRASRdFsLAJtRFCdqP9sm6KTQL2864Rub432fEohhTaMpA1TkMm9e/eOGDEiKyvL6/Vyzg8cOJCdnY0Q6tix408//QQdOsEUB5ouq9WqKAp0VwZomkYpRdBrl/FT0Ek4ua4szFIe8hRaAUnTyaFQ6IsvvrjiiisQQlVVVU6n87///a/VanU6neFw2Ol0+v1+k8mZEAI8m9SgBjMSFa/f7w+Hw61iV7dufXItiwhHiB1NqJlaq6eQBCRNkletWqVpWm5u7vbt2zHGX3zxxeWXXz5nzhyEUCQSYYzl5OQk2sZVVVV+v7+ystLv9yd2chNFMRQKATN2qyvkYxFfHzcwzBSJpRcp5ZxCSyFpklxQUJCZmblq1apIJEIpLSoqisVilNLy8vLs7Ox9+/ZBg2xBFKBtss1me+eddzZv3lxYWAgNkwGyLNfU1NSyYSe2NE7WCvkU1yc37AvV8g3fUzgDkDRJHjJkSDwej8VimqZ99tlnTqczPT191KhR77333llnnVVaWnr99df7A/6MjIzKysqMjAybzXbbbbcJgpCeng6d3wB79uz56aefBEEIh8Mt2HW1tfxYHLOm+ramkEISkTRJxhhbLBZJkiil55xzjs1mQwgNHz68tLTU5/P17NmzW7duoFfbtWuHEKIGzczMRAhpmpbYj6Jbt25utxsOlaxrax6SWZ/McENdzmopOOur8VTfxhSSgKT1T47H4wghjDGltGfPntBiIiM9A/rIQCMoznlZWVl+fn44HJYkidQpRlEUCa7drqioqKioiMViiqJQg6K6nqyJoankoWUlh9UZ8gyzFH1ICi2KpOlk00ElCiLjDNW1gwHZs1qtIOp5eXmMMVEUMcYEEyIQdHTQWBAEMKpNkxsOmyjPiTgJ2T6u+uSGnzSej80RbaJzcgoptCCSma157H0EImCCGWWMs8T+b43+kBCi6zoodlDXifL8CwLDp6fD2nxrrZx4e1znbe0mgacUyaygMMWMsMatVs4440ySpETpPcYsADKcKL1JfTEQJWoRQeMIMXx6tng7Rr1qi4pNvfPWG0Itd95fClqb/QeyMk/lZFm/f3ILwlwnA/hpFFD+pejkZu55GiBp1jVqHpGAmePV1uZUjDGqvaSjDYom41UN4smYIYQwEhBBGENSCEYYY4IRwbXtkxGkuZ36+z0xwFur/Z9xhBBDDDU2AJL7ThPPm/hJw/OGw2GHw4EQCoVCNsXGEScCqWf3tYXxlnS0qk6uZye3vVKn5Ncn87pw1GnmBjPT72r9F0Ir3Z553qY8Jg6HAwp4HA6Hruucc8qoQATwrZ7GSOo62Zzqmnho9eZC88+mRLphMXPzlX8z0IL1yQxYNNvaTJUMcN5IWcvP2rE/axX/7IKr3nlrt49ZIStJEmOMUso4QxSZdXinJc5MltyWWrFyRBJ0L2GI47ZaFHmSrgoINKI6r+QxzsIoS5yRTZnXdT0x+cfcDSFECIHE3sRDYYw54uZ5BeFIqVyiRRAKhUzr2mKxiKIoSRIkJpzeSNo4w83Ayfw28QjNOWbTqCfDCfHkWkpq0mSf5IZI3Kd2gzCw0o/w4xJ2zOF+qgARfl3XwQo1qR1UVVVVFSFEDQr/zJ8YhsEZh+R5SqkkSUQgsBCth+rq6ng8zjk3qEEp5YwzyhhlcLrS0lLGWDQahYPH43HGGUgmCGdcj8MZYX9Yi1FGDcOQJKle/h/GuLKy0jAMVVUlSdJ1vaysDNIEOeOapiFUq5BPdMz8AtAWR1hLw0y3OiJgPEGMk4O2/mAxxrFYDCGkaVpifThKYGsRRAH+IYSi0ShCSJIkQRSqqqpsNpvVakVNz6qZmZkYY0qpLMsEE8YZ/FNVVZbl9u3bWywWm80WU2MVFRVQLWfuA4dljH377beBQABK3+Px+IEDByKRiK7r4XAYhN/cv7KyUpKkUCgkiVJNTU16erosyxjjYChYVVWFECKEnK4yDGhdj1cbeJQgxgwjhhHHkLmRZKkjvL5/i3Azh7sN6QRd16lBV61apSiK+aHT6QwGg7ABXKiAQCDQr18/p9OJEIpEIlarNS0tjVFG2RGNnXhrZWVlGRkZYOgm7vPee++tWbPmjTfeUBTl/vvvP++8826//Xav1wtZQCYsFouqqm+//Xb//v0nTpxYWVmJMX7ooYc+/vhjhJDP50tPTzd3DgaDq1atWrZs2YMPPjh48OBvv/22uLj45Zdffuqpp9avX//6668vWbIkeY+tjeJMXCdDFjQIG0cEIcQxYZgmIzrUSIJnwwhKg58wxli9odzS0DTNYrFocW3x4sWgbwE1NTWmJNeTcIfDcc4552iaNmzYMHOfxDq2RElWFGXq1Knjxo2rd95Zs2Z98803iqKoqnrDDTc88MAD/fr169u3L+T2iqIIGhhjrOv6smXLZs+eXVpa2q5dO7/ff955573++us33nhjdnY27G9e25VXXnn99df/9re/LSsro5TOmDGja9euubm569evv+aaa5L76NomWkSS247aqQeOEMNMqJM3VY8rik22KCF/wCZZEEeEM+DzIHVOGowxqgs1g9nJmuA/4JxzTDlnBkHMwAwRg7P0TGc8bsjWtJgR5wyOeURiCSKolrqMISR6vd7s7OyKigqPx9PSKzpREEVFpIy+9dZb9b668cYbO3bsOHXq1IqKCpDncDisWGulGjZKSkpg6rn//vvvvffePn36zJw502Kx/O6B32GCDcPo2bNn+/bt4/G4JEmyLEONzcUXXzxs2DCbzeb3+xVF6dOnT//+/efMmTNz5kxJklRVBQcVQqiysvKjjz664oorNm7ciBCqqakpLy8fPXr0zTffPHDgQMMwgsGgx+Pp3LmzTbFpmpaXlzf5d5MtFosgCKqq5uTkPPTQQ9SgX3/99aBBg9auXatpGqjxHj16wLrgNMOZqZMRQ4gj4g/4/aEwEsT03Nx4jR9BOgeY3wxEF5wwtW5ejAnGmDQWyeQEEyRgghFBMhaYIHCCRUkyEFf1eHn5YWdmRiNXghhiTNM0RbExZmRnZ1ODZmVmxeNxQRBatK4TVr9CA54Dn8+HMdY0LR6PZ2dnQ2wWeMtNwAMhhFBKgfIFYyxJUjgcDoVDGGPQqybAn/z222/n5ua+9NJLmqa53W5qUI1rb7zxRv/+/Z9++ukpU6ZYrdZQKOT1et1ud2Zm5u9///tdu3bNmjVr44aNwJcOS+sXXniBcy6KIuf87rvvHjZs2Ecff3Tttdd269btmWefiUQit9xyi81m+/Wvfz169OgdO3bMfWeuFtfC4fC2bduuvvrqKVOmFBQUtLIF1Ao44ySZH6k0RBlZWTo1fD6f06IIkIvFa8UUEjo44hwhjhjQl2AkYNS4Dc4oZ5wyxHVsMI51hCnigmSJM8PhdNrsdtOeh0mBIMwQN6iOOVIUWywWVZS0jRs3vvbaa7NmzdI0DapBWx8ul2v37t2bNm3auHHjuHHjJk2ahOrypcx9fD4fLIBLS0urqqocDgchRJZlq9XqdrsbHtNqtX766acvvfQSrFeBUFkQBKfTqev6O++8c88997Rv3/7OO++klHbp0qWqqurBBx8sKChwOBxTpkzZunVr7969EULg6I5EI7Isx+NxXdczMjIYY5dddhnnnBr08pGXd+rUSRTFJ554omfPnpMmTbrkkkvmzZtns9m+/PLLJUuWPP/88+FwOBqN2u320yyXs0UkuanMuLaQMScTyTjigCG7i/f6qmuYpoocowSzmeA665rUtshBCHGGOeesyXoPjBDiBHOMOBE4InFKHS5nZcCX5nREdQ0hxBDHCHPEEcIccUqZIivhSNCe5nzssT989NEn4XA4FAq1qPlnmhhNhZQzMjKee+45zvkzzzyzcuXKhx56aODAgWYgFyGUl5dXWVlJCMnPyw8EAoqiYIxFUQwEAuFwGHQ1qluMcM6nTp361ltvffLJJ/Pnz/d4PD6fTxAEWZZVVQ0EAr169brppptmzpy5bNmyl19+2e12f/nll6Wlpbt377ZardSgX3755ZtvvvnMM8+4XC5Ylv/lL38pKSmZPXs2xjgcDoPJMHTY0Jdffvnpp59+8cUXJ0yYcOedd27fvn3r1q0YY2rQH3/8sWvXrlarlRDi9/sVq9KWmigkAcmUZFVVIdyXaIyFQiFZliGyl5gxd0qkGgjyBCJQaoRCIQ0LH3+2hCAWi4Up51hExtE2F6OMCMSm2Hx+nz3NTgSiWK1VVdWKojT0YBFOBIYIxoQjhMAlTihGHMNhWdzQGa9LU8UYIWQYOpKVue/MnfnqTL/fzyhyu93gH070JJ0MKKW1eU6MCYIARrt58Q31UiAQKCgoUBRlyJAhw4cPnzlz5sMPP3z++ee/+OKLiYd1uVyhUGj+gvlXXnmlIAgCEXw+H0LIYrEYhuFwOAzDqJVjxr/++uuVK1dKkuRyuSilS5cuHTJkSIcOHVRV7dmz5yuvvPLqq68uWLDgz3/+c3V1dXp6+jPPPDNjxoydO3cihARRuPzyy++8806n0+nz+Ww2m8Ph2LNnz4033igKIufc4XAEg8HCwsL333/f5XJdfPHFHo9nzeo1Y68Z++CDD37wwQebNm264IILiouL77333lgsxjl3u92UUZGIjT6BXyiSJsmxWMxisVitVpN6HgBmGEIIEdSUn7M1oVI9HtetVmtOTk6wOiAQESGGRJljakiEHs25BS5l3SIYFtGwitFo1LCIhiJaMl0QjDXBOUcUI8qYwRDjGNf2ssMEI1EQERIlQhk7otEx4pwbBh09etSPP25WFMXtdh88UBYKhW666aaqqqrECNAJAx64JElWq1WW5RtuuKFPnz4FBQW6rkNRmqmcze309PRoNJqWlnb48GGPx/O73/3u1ltvXbVqFUgl7OP1ejVNKy4ufu1ZwVmFAAAgAElEQVS11xYuXOh0Og1qOBwOcFNZLBYwhusEmX311VcQvp44caIoip999tmYMWMGDx6MEAoGgx988IHVaj377LMXLFgQjUYVRZk3b16PHj18Pl8gEHC5XL17947FYrt37+7UqRPnPB6Pl5SU9OrVSxCFeDweiURcLtfhw4cxxjNmzKiurrZYLJRS8MkXFxd//vnn/fr1W7du3bRp0xRFqU0UwZgzjgk22xL+0oU5aZKsKApjDIZIPX0VDAZFUQRmr+NCUhlCaqGqKud88ODBvXqeRThx2uwIMVWNxomBRKGBTqZEwDZF8fn9dnva/PkLa6prFEkZMnhI+/btj7owgglHIiciwgRjxjBCiCFMOaeIIYQo1TsUFDDKJEmmrNYhRKkxffpf33333VlvzrJYLHl5eaIoTp48ORKJnMDjaghJkiwWC7TsisfjWVlZubm54JEC2paGz9Pn88HZ27VrhzEmmOTk5Fx33XXUoAatveyzzjrro48+mjFjxh/+8IecnByCiWSRzjnnnPnz5+/bty8rKwsh5HQ64fgEk3A4DCt/QRAOHz7s9XoHDx4MV+V0OjVNi0aj0WhUkiToc3D++eeDaIGzQFXVUaNGvfXWW88++yyldOvWrRkZGZ06dYrH47IsR6PRjRs3vv322zt27CgvLz/77LNvvvnmb775ZuzYsTNnzhwxYsS9997bqVOn3r17w9FEQYT4NuOMMHLaVFYk07resmWLz+fzeDxnnXUWIQQUwtatWyORSGZmZpcuXTDG8XjcYrGA3m5KLFt0jiREcCq2kSNHpjvdcVVX7FaEkKHGDEyxQOpRfIBOtgkWX8xvV+wb/7tp9epvLJJ8yUUXXXrRiMTDcsgO45xATjchnDOGBARubYwx4owZjDKEkECOPPbCjoXTp/9tzJgxTz715IYfNlmt1j59+oCBevI3a/b6EEVR1/VYLAah48RAcT0oipKXl3f48OEOHTrApFybjEkwN2pV96ZNmyZMmNCtW7fzzz8fIaTrek1lzfDhw8877zxKqc1m45zfc889uq5rmgYGsJlH+Y9//OO6665DCFVVVXk8nng8vmXLFofDIQiCxWKBlA+YNcLhsN/vt1gsFovl/PPP/+CDD0KhkCiKa9asOf/88zHGkUhElmW3292/f/+ePXtOmzatqKho9OjRxcXFO3bsmDp16tatW3v16tW5c+dHH310xYoVkPhpkS1HghG/NAqaYyBpklxRUVFWVlZSUrJ58+aNGzfefPPNoVBo//79ixcvvuCCCz788MNbbrmlZ8+eP3scWMVBop/ZLCqJsCs2hlBWRiZCSEkTEEOIMNFqgTOxo/g0ERIQR5Qh5lKcCCGkU5mI2KAyEUkjhU4EHNxIQAgRjAVyVDEjJkQSECV1xyecZLgyOePUoAMHDvpiyZfvvvvuvHnzIEiTlJvFGJs+CxCJn/0JOIRkWbbZbGCImmyK5jFhGQxiTA1KMPF4PEBjjhAC+TddCZRSzrnH42GMHThwYOnSpZ999hlCyOPxBAIByC3p0aMHzPsgxoIoBEPBzMxMp9MJ9Olnn332E088IYqiy+X68MMPFy1aBDvDr2DK2Llj5759+9xu9+zZs0eMGCHLMmiUyspKeAiQsA3KH8LXBBMstKGsu5NB0kQlJydnwIABv/rVr8Lh8P/93/8hhFwu14IFC+644w5FUXr16jV9+vQZM2ZYLBbgzURNpD0FAgGokoF0gmRdXpOol1Z5zD0EjCUicMYFhIVGYlGJ8anaLaEeC25CGLN2m8BuAkII1rFpaWkNK4FaDbIs+3w+t9ut67osy8C7pmmarutWq9VisWCMPR7PDz/8AMoW3NTwWg3DoJSmpaVxzisrKwVBsKfZKaOxWAxjHI1Gf/vb306ePBmMMsMwnE7n2rVrrVYr59xms4HrAcZGeno6pMc4HA7OedeuXV9//XWXy7Vhw4YePXrALKOqqqIosizDJPXGrDe2bt06adIkVVUnTJgQDodtNtvdd99ttVpXrlx58803f/jhh5CPLUlSMBh0Op2qqlqEn5/dfhFIZl+orKysaDT6ww8/DBo0CKRU1/WcnBxYn0B+38/C4/HYFFssFovH4+DFbYnVctM4Jp9mbXkTPy5OguavwywWS//+/RFCuq4znbU643d9GIZBMIFOmvDMVVUNBoOqqq5evToej4NVjBCSZRla8xmGkZaWxhjr1q2boijBUNAwDI/HM2fOnFdfffXKK68cM2YMeJtXr17tcrnWrFlz0UUXWeSjxAlWBD/88APGeMyYMfF43GazZWRkzJ8/f926ddXV1ffdd58kSdXV1YMHD77rrrtsNltxcXFZWdmkSZNuuumm/v37P/roow8//PC6devcbve8efMQQrNmzTr33HP/9Kc//eEPf4D1ObjTTsVzbREkje/a4XBUVVW99tprlZWVt99+O0Jo7969uq67XC4w0mRZhoA++MZgwdzwUIFAIK7HoYdj2yMVST4SA0IHDx7MycmBxwI2ausbfrqup6en+/1+WNpAKSJ8BY4iu93+8MMP33LLLZIkGYYByZtmnSNIMuc8Nze3f//+8XicUvr6669//PHHr7zyyvDhw8FW9/l8ixcvhlStcePGCaKgaRpE0SEhPBQKDRw4kDE2bdo0jHFxcXG3bt0OHjx42223KYpSXV1ttVpjsZjb7SaEFBcXT5w4sUePHv/85z8HDRqkadrChQtvueWWJ598Es4oCmL//v3XrFnz5BNP7tixo0uXLqqqulwu8Jm18hNuISSHgocaNK7H4RFnZGRMmzZt7NixGRkZr7766tSpU8Ghfccdd8yZMwcCzowxVVUXLFiwbdu2/Px8qIYFCIKwbNmyFStWoLogSsPRnLTxXRcQSkBDe/6ITp006Z7//vBfzvmLL744fNjwY+z582euS+ROSDupLaOHKgLQdXVZoi0uz+b1UEoJJvVcQVAzTBltmAUA/ze8WnAvWa1W6JJts9lggWpSMkJ1tCRK4FGDnn4IIbCZEULhcNhisRBMzFIq8+xg0oN7PxwOI4QkSRIFEVK+4amC21wQBfNcsJyhlEJGGjjPj+H8+2UhOTpZEAVFVAKBQGZmJjWoKIoFBQWyLFdVVa1bt27gwIEbN27s0qUL59xMS7DZbN26dQPXSKKREwqFqqurNU0jhCQlptrGYZLamElUpmCjOjaM1sxeMLk44HqguwAWMOdcJCLIeXOOo+s6VEo4nU5ZliGbWhRF4AYx2+4yykD2CCGQrQ1U57CyaOpctTJp1C7RTd0AXlJoaQKTUSgU4pzb7XYzG1wQBDis1+vNyclJxjNrE0jaOhkWKoWFhT/99FPnzp3BD3HllVdu2LDB7XavXr36wgsvNMci0LsUFRUVFRUpigKjB77yB/zvvfeeqqrhcDgvLy9Zl9c4GhGNYw9TIBLhCXwgJ30JdS0pE0nCT+EU1vB6QPYSdW9zYLVa/X6/2+2ORqOMMZfLJctydXW1zWYjhIhC7cAzeQXgT4IJxxwkTRCE2oTKBp14QOUCrwjGWFVVQRAgE5NzrmmawQ3ofADBPIi9gTkA7Hyapp022hiQNEnu2rXr3r179+3b16VLl8zMTEhIGDly5JYtW5YvX3755Zf36tUL1akXalDGmcPhgJxYMJYAnHPDMCwWiz3NnqxrS+EkwRk/rokLDA2n0wkWNUJIFEWfz+d2uRFCBjUMaphMyfVsDfPPRPGGDYZq0xBA5YJjzGKxyLIMS3TDMGr7lhAiy7IgCLU2fCJhGGeYY13XnU5ndXU1tBk8DZDMgO3ll19e7xNwxoI/NhGYYIi71CYeJEy6kBtMyOmTfNMozCXlEbdW07fbCnZ1M6+n4ZUkSmPit7D6hSwR2A3WUBhjEYmMM4IJwqh2wzx+3cBAoJPhsHVLDMIIQ8xsJ1ZrgdeJLmNMlmUwuRlnMIrA1QKzhum6M8Psp40Yoxatamx0/MFYqcdcnwJq+nG1/pUc49THSMtr6k98NLE55zxRXAk7yruWmAde20so8cjkqP3N2Fjtl4SgurprzLAp4UddWN1vTxuXtYlTUJ98eqTUJAtt7Wkc7/U0c39zt3ozeD0O6sTdGrLkomZT5DdKJGBaeW3tmScFp7MFm0IKZw5OGWfIMew0E618SSm0EOqlbTe1T6PbDY9zYgPj9ChdPAbOOPafFE4YPyuNKZxCtDlJ5gmAT1LjJgUTJzMYTu+BlFonn+lofhCh3gzbHEAQKBWnaAW0CUk+9tq4oZZOIYkw2Qig0AUlPHBq0Gg0CmwB9SpMgbWrsrIScubD4XBi8jxCiBq0ts7RqGXnS/wWN0BL3uIZgTYhySmcQpgJz5DSaIoc5FHYbDZo7KJpmtmBLR6P2+32srKyr776qri4mHOelpYGtLXmYQVRKC4uFohgUMNMyUih5dDm1slADJRaJ7cmFEVZtWrV888/D4V+5ufBYPDw4cMjRox4/PHHE/swQVrF999//49//KN///5bt25FCDkcjv/973+KohR2KOzYqWNpaWnnzp3//Oc/Q8WSme+ZepsthDYnyacBzFKeelWZbWcQmxmU8CeQbPfv3/8vf/lLoiTLsjxjxgzGGIjxEcprAz3y6CPFxcUffPBBenp6RUUFFBWNHz/+7rvvvuSSS0KhkGEYmZmZmqYZhiFLp1tCVRtESpKTD8MwfD5fdnZ2bYFXnS5qU/FMs5oSIeR0Ojt27Dht2rR77723srLS3Mdms+3atWvkyJFApgctl5YvXz579uxLLrlkxowZXq8XITRz5sxzzjmnuLh49erVl112WTgc3rx5c6dOnW655RbGWNu55dMbKUlOMoAHIzs7m/O6El9Wy9fVdsY0iLGpYw+WHAyFQqNHj77vvvuAcQmgquqiRYs0TTt06FB+fj7Q5YwYMQJj3LdvX4TQ4cOHCSGff/651Wq12WyvvPLKoUOHotFoMBh8//33r7vuOkrpmVBk3haQkuQkA2MMpvX27du7det2qi+nWSgsLDx48GBxcfFbb72VmLEcCARWrVo1ePDg/Px8XdcPHz6MEHI4HFdeeSUhZO3atffdd9/kyZNramqeeOKJ6upqr9frcrmGDBny3XfflZeXW61WYKIH93UKLYrk+K7BQZXorKoXOkr0ap72YIwxxtauXXuqL6Rx1Bb6CbXV/PBh+/btu3btCrQbJkBohw4dqmmaJEldunSBdjCU0v/85z8PPvjgiy++eOutt8bj8XA4nJGR8Z///GfDhg12u91utwP3C0IoceGdQsshaTq5dt2VMDNgjKGvUrJO0TYBUxWwyXLOZVl+9dVXP/jgA1he1tuzXn3fke0mejKfJIB5SxREs64IplSCCREJcPHs3r173759eXl5gwYNisfjVqvVbrf379//66+/vuiii+LxuKqqGOPy8nJPlgd6sj7zzDPbtm2bMWNGUVERQigWi9ntdoSQ1WqVJAl+YrFYgOAWDHmCyWnWvaWtITmSbK67EilsEKktFgWCa7PmOylnbDvAGANFY8gf2rFjx5QpUw4fPlxdXT1y5Mh6nBjHKOJtobi+WQYcj8fj8TghxCJbBFFQVZUyKsuyFtd27NixePHiTp06OZ3OwsLCf/7zn7169ZIk6c4773z66ae3bdt23nnnLV++fP/+/UuWLElPT//3v/991llnPf7440COBRkgmqZByxiMsSiKJmtqSlxbDcmRZIyhk+gZClmWt23b9tZbb82bNy8rKyscDnPOvV7vN998gxCCdQegSZ3cMulrplBBNKh9+/Y9evQQkCDLMhDf2Wy2gQMH+v3+cePGQTuIkpISt9tdUlJy2WWX3XDDDfF4nHPeq1evxx57zOv12u32q666CpjoZVl+6KGHRo8e3a5dO1VVMzMzgeuHEGI2tRKIQFlqkdwaaFmPF4QuFEUx+06cNkjMXqSUlpaWlpeXy7K8c+fOvLw8Xdc3b95cU1ODjmbGSTRJEo/QQpJMKdV1HWMMBJTnn3++1Wrt1q2bYRigPFVV3bFjx4svvvjtt99yziORyIEDB0aMGAHNDcePG1/hrXA4HNDrHFq3EUJg1fDuu+9+/fXXkydPrq6urqmpYYwdOnQoFotBF+XEy2gmHWcKJ4OkSbIgCuWHyiVJikQi6enpbrebMRaLxUpKSgRBCAaDRUVFzelL9EtB4sqWc37ZiMv69eu3du3aZ599dseOHTabbdSoUXfccUd6evqpXVA0nCOoQU0Hu9Vqzc3Nve666yZPngycW9OmTUtLS7PZbDk5OYv+vYgQoqrqunXrXn75ZUJIRUWF2+22WCwPPvjgt99+u2jRok6dOq1fv75Tp06VlZVffvllSUnJ8uXLu3Xr9tRTTwHRIjBs1lJzIiGVZd1CSJok79q1a+nSpYZhZGdnR6PRu+66S1XVLVu2fPPNNz169PD7/eFw+NJLLzX3P51ep0AETHBWVtboK0f369fvpZdeWrx4sSiKmZmZqqqqqmqz2Rq2hmlpj1ciP5ZJnJZIowW+Zc55u3bt0tLS4Apzc3NBkxcWFgLRLBjhmZmZfr8/JyentLR0ypQpVVVVn376aUZGRkVFRY/uPXbv3n3XXXdNmjRp8+bN8+fPj0Qif/3rXyORSEZGhtVqpZRCmhdnnIgp/dwiSNo6ubCw8Oabb4Ymg4888khpaaksy5s3bx47dmxOTs7BgwdXrFhhSvLp57qsdfhhkp+f//zzzw8ZMiQcDh84cKCgoACapDREK3i8aq+tsST2QCCQnp7OGd+6deucOXNWrlwJTO6HDh36zW9+4/P5Pvzww7Vr18KS2Ofz/frXvxYEgRr0ww8/7Nmz55QpU6hBg6Hgzp07J782WVXV8ePHT5w4ESE0d+7cuXPn3nLLLWPHjn300UedTqcZTyYCOf1efRtB0nQy0A4jhFRV9fv9giDk5uauXLly0qRJkLg7Z86c41oNtvGePY203qW1vnqr1TpixAhFUQQiQGMkWZKP3Xux9Qs2XS4XYwwTnJubO3369Isvvhj6Jy5dutTlcuXl5U2ZMuV3v/tdIBCwWCyU0mnTpsEPb731VpvN5vP50tLS0tPTv/rqq1tvvfWcc87Jz88PBoOMMafTeeutt/bt2/enn36qqKgAzvpWvrszEEmTZFPwDh8+nJGRkZ+fv3r1avByBQIB6A9k7swZZ4glLiAhzgm9Zmw2m67rQF+M6hR425/IE+1nsz+OGQdqiFNebg3Pv3v37t27d0d1bO9jxoyBbydPnowQMjs5m5IMn5jC+dRTT5kHTMz0PPfcc88991zYTvTttf33+AtF0iRZlmVd17dt2/bSSy/NmjXL5/MNHTr0448/VlXVbrczxkpLS83VICa4vLx84cKFy5Ytgx1g6AcCAYzxrl277Ha7WfjeNnHyI7I1x/SpnTVS0tsKSGYU6sCBA8uXL7/nnnui0SjUwXHODxw40LlT58MVh4cNG2buiTHOysq67bbbrrvuOohtQtwiPT29pKTkL3/5C+c8Go2eTu1tTy3aiCy1kcs4LZE0SQ4EAq+++mpBQUFRUdH27dtBknv16rVs2bK+ffv++OOPI0eOPFLgSpEkSXa7Hfp9EUI8Hg9CKB6Pd+7cWdO0tt+lOjUom4/Us2oFJE2SoVGjx+P529/+ZrFYXC5X586db7/99tWrVxuG0a9fv969e5tLZYYYM5goiFjCuq6bth/nvLy8HCHkdrtVVU15SlJIoZlImiQ7HI4HH3wQcrnC4TDwPEFePqW0XmkbiC7lFHNMCDFzgAgh2dnZUKCTkuQUUmg+khbHhBJWEMK0tLS9e/dC92NI3EcIlVeUowakigjSkhOagPv9/kAggBK8pimkkMLPImmSDAmAkKCHMW7fvj1CKBqNQvf6eDyen58PzebhX72fm9kLGRkZ0Pb2lAdpUkjhF4QkV1CYZRJgGEMqL1TtImhyixiUp5r1j9ARF7bhK0EQRFFMSXIKKTQfrZ0Ee8SwrsuRAjEGgTdb5p5+ZcwppNCiaG0er8QSP7MX7pGsYIJRAiVNCimk0Ey0OdVnNi5pIUKcFFI4LdHiOvnYvXAb7mPy+DHOjpG0nEIKKSSizbHk1iPlTOEMROKrT+WHNROnUpLNl3QmC20rj9pUt63TFW1unXwmwxQzatB6XU7PKKTsshNAm7CuE/VDo211Tz8FcoSLJ8FfgBCKx+MWiwUTzBnn+KT0Z2I2u/lh4qGqqqqysrKCwaDT6YxGo4ZhOJ1On8+Xnp5ODRqOhBN516ChhM1mq66uzszM3LJly7Rp0+bOnTt79uycnJz9+/eHQqHc3FyEUDgc7t2796hRo8zfbtu2bfHixfn5+aFQKBKJDBgwIC8vr3fv3olX4vV6O3TowBiTJAmogiCQUY9pOIWm0CYkOQWEEGdcEAVJkkBUTDE+MTDGOOPAXG+y84qiaNIhUINmZWVVVFRccMEFZWVlnTt3VhTF6/VKklRZWWm1WhVFSZTk3bt3Q/M3URQRQsFgsKamRpblvn37rlq16q9//euMGTMqKioQQl988QUQ2ZvYu3fvN998M3LkSF3XFy5cqKrqc889N2vWrOrq6lgsZrVaPR7P/Pnz33333fz8/Gg0mmK6PwGkJLlV0XBE1vbuQIhxFglGnE5nenq6pmknWT1CCEGkPmOJaa9ijBlnIX/I4/Hs2bMHIcQYg1YyDodj4sSJHTt2fPjhh0GSG15wRUXFrFmzXC5XOBweMGDAhRde+Pbbb48bN87hcMTj8TVr1iTWoiOE/H5/Wlrao48+GovFlixZMmzYsOzs7LvuugshpGkanGXDhg35+fmc83A4nJmRiep0cmICQgrHwBknyYyxk0wgg9zyE4MpSA1XgJIkSZIEtq4syyephY6E8vgRunzguGaUIYIopYFAwO12+3w+m80GTRWBcX7//v1AcF3vAgKBgD3N/txfnsvLyxMEYc+ePXfffXfv3r3/9Kc/RaNRh8NRXV1tsVj27dvXpUsX81eMsT59+nz//fcIoQULFgwcOFBRFL/fP3v2bFVVFUVRVbWysnLTpk1whR6PxzCMlE4+XrQ5j1dCNiept1o+edTmnHAejUZhIxgMIoQ0TdM0TW8A4LiF30aj0Wg0Gg6HQZKpQZv/D45ADcoZZ5RB1xXOeCQS4YyXHiiFHV5//fWioiKfzxcKhdDRjp/jBTxJMNcNwzAMA2YK0MbA0FJQUMAYA/5qq9VKCAHGJVEUV65c+cADD+zZsycej8MjggcSU2Pjxo3LysoqKSkZPHjw//t//+/mm282DAPy7QVBKC8v79y5M2PM7/fDZcC3hw4eKikpmTdv3qOPPpqWlkYIicViFotl7969qqqCMg8Gg/v27YNXj1I6+ThxZulkUxvv2bMnGo1mZGSEw2G32x2LxVBjPavMSmmEUFpamt1uz83NjcVi1dXVO3ftPK5THzhwoGvXrh0LO2KCLaIFIUQZhTGdn5+/bv266dOnr1+/nhCSlpYmimJSMs8ZY5RS0LSmj41g4nQ6TYHXdd0wDEIIsPmGQqH27ds/9thjS5YsGT9+/E033fTwww+XlZVBuzaLxdKjR4/u3bv//ve/37lz57XXXnvxxRdv2LDBarW+/vrrEydOXLhw4erVq/Pz83NzcxctWtS7d28gGJUtcjAY9Pl8hBDDMHr06DF+/Hir1WrWri5dutTpdAKnX6IMp9BMnFmSjBBijPl8vjfeeOO///0vxjgSibhcLpDVhvneIE6mJFdWVgKVwltvvQW0782HzWabPHkymJ2BQEBRFFmWqUF9Ad8TTzyxYsUKXdclScrMzPzhhx8uuOCCxGalJzymKaWMMdNHxTl32B1EIGAmgMDIsmzSG8fjcbvdXllZKYris88+e+ONN06fPv3CCy+cPn36+eefDwe0Wq1Lly4tKCjIzs5+7LHHZs+ebbVae3TvMXPmzBEjRtx777333ntvOBy+8847gQEqHAkfPnyYc96pU6cpU6bMnz+/Z8+ePp9v9erV3333Xffu3TVNs1qtoVAIqClQg1jGid37mYYWkeSmVjXNX+0wxoDlHP5M4uvEGGdmZhJCao3kOqBGG69QiupKNcPhsKIobrfbMAxokXNclR4wKQCXMPDRRqPRVatWzZo1a9OmTaqqgu/H6/XedNNNkiTpum6yi57w7UPriUgkghByOByTJ0++7777LKKl9EBpYWFh7T3WGf+yLFutVmrQgwcP2mw2wzAKCgpee+21pUuXPvfcc0uXLpVl2e/367o+Z86cYcOG5eXltW/f/p133ikqKho2fFjnLp3Xr1/ftWtXhNArr7wyZMgQEEtwTVNKFasyYcKE7t27T58+HWNcVVVFKV23bl2PHj1EUbz++uvT0tJCoZDD4Tixmz3DkUxJBs5qQRA44wY1oBWYmeEAg/hnB6XD4YAIilkphRrI2EkK9qZNm7p06RKLxUKhEK5rvPaza3JgDozH44ZhCILQqNfqGCCEQK8shJAZ3bn00ktHjRr10ksvzZ49GyG0Z88eRVF27dqVSBydiBOb2hKvE6jUOnbsmHBliFIKXgngss/LywsGg4qigL09bty4cePGGYYRjUYddsfBQwcLCwstFgvn/Iorrhg6dOiMGTMuu+yyiy666MYbb7x+wvXhSHjlypWfffYZ2DJAtyjLMiTSDx8+/LvvvhME4d5771VV9YUXXnj44YcRQr///e8nT5785ptvorqgOkrwbKfws0imJIOfg1KamZkZCATsdrskSTt37mSMVVRUXHjhhU1VTSRi586duq5Dj9/ESEyy1DLGuKioCOgBMcagV2sdbORYkgxfUYPCJEUwOUIV2gwIggATAecclo6m0Ttp0qRLLrlk7ty5S5cujUQiDoej3ghOYrZTQwoH8yHA54nr88SnEQwG3S53KBwaPnz4pk2b3njjjfLy8nA4bLVaV6xYMXLkyL59+44ePfqPj/3R5XI9/cIYCi4AACAASURBVPTT5g8ZY7t27aqsrIxEIqIohkKhLl266Loei8Xi8bjf7y8rK/vjH//YpUuXF198URRFEGPwydWLS6dwDCRNkquqqlauXLl7926n0/nAAw+kp6dLkrR9+/ZPPvkkIyMjGAxaLJYLLrjgZ4/TtWvXdu3aIYR0XW8JRj5wAlFKISQjiuJxzRGCKGBSS4RwAqVatZ5YgjHHCCGCCXTn6Nu378svv7xkyZJPP/0U+jOZV9v4EZp9zcmaAiRJEkTB5XK999576enp4Hl2Op0HDx5MS0sbNGjQ9u3bp0+ffvnll3fo0GHo0KGmzx8hBG9/8+bN77//vsPhGDRo0Pz58xVF4ZzfeeedRUVF8+bNu+qqqzRNq6ysdLvdwOiIEDI3UvhZJE2Ss7KyLr/88vT09K1btyKEoOXq/Pnzx48fn5OTU1ZW9t1335mSfIwFMzRGAy8IpEwkN2sPjGSCCXToBnV3XLEuyKNsvnWdWChSj8CMMkopNQwDFNGoUaMGDRqUnp7eBnlFHQ6H1+vVNK1fv36MMUVRPB6PruuLFi361a9+FQ6HCwoK/vjHPwYCgS1btnz00UcjR440f9u7d++cnJzhw4cPHDhQ07SamhrwgUmS1K5du6eeeuof//jHVVddRSnNzs7WNM3sQJIS4+YjafFkalCn09mzZ0/OOSyYYaNPnz45OTnnnnvu5s2bzZ0549FoFGI/AFYHWFYRQkx1xBk3WQcaDZweLzjnRCCSJJmmNToeC1YQBVMxNgfmDyFKDttASygQQRRFRVGguamu63a7nRpUEqXaWSzhCI0e82fRzItMPHI8Hq/XXxImu/T09KysLIgkIYREUWSMzZo1a9y4cV999dXQoUMVRVmzZs3777//wgsvPPLII+vXr5dEyWq1BgKBQ4cOgVcPY7x9+3ZIGp03b95rr732wAMPDBw4cOzYsevWrUMIlZeXAx8rIQQqSZK4uDiNkcxO6DDdlpWVQUBF07Rt27YRQoLBYJotDbJ2OeeEEEEUbKJt06ZNoVDorLPO8vl85nGqqqrg56qqwkK0lugrQT+fDE7tmGhKCBuSjR570d6iyM7OhmwQAAgSowwhJIpiRkbGwYMH8/Pzd+zYMWvWrAceeODJJ5+sqqp64YUXunXrFo1G27Vr99lnn7388suff/758OHDUUKhiM1mk2V5zZo1HTp0OPvssy+66KKioiLO+VNPPvXJp59MmDDhzTffvPrqqxMjcCk0E8ns1WixWNq1a9e9e3eMcXV1taIovXr1QnUt/LKzsxEIJEEIoT179nzxxRfr16+32+2J+Y9ZWVnffvstWNeJYc+WQOuLCq5LPMRCrXEOi+1EqqOmcpta7WqvvPLK3NzceisgTDBhBHI/3W53bm6uJEpjrxkLnjlI+TS9dLIs/+lPf6qsrATXnWEYEydOVBQFIsZjx47t0qXLrbfeCtHpUCiUlpZ2zjnnlJSUQMAMNTa1pXBsJE8nEwFckYcOHbJYLLm5uYyxmpqa4uLirl27VldXgx8LQhGMsfbt2z/yyCNgZUGTYThOJBI5cOAAxljTNFE4UxJXWm62Oi6A6I4bNw4C0Y3uIxAhLS3tsssuM/vsapoGSSCwraqqIAh2ux04zzVNG9B/gCAKfr8/IyOjqqpq4MCBRUVFkD3q9Xo9Hg+jLDs7G1qFRaNR8OqncFxI2iMzqIHqUikopYQQn8938803L1++/ODBg9u2bZswYQLG2GTbkyTJDHgkRlzC4bCZO9lwfJ+8Xko8QrK8aCdwAfU0Hm+shLjhRotekgnzvST6AmtjVPjIV6ZPLvH1EULqtebDGBOBMMbgc1iBm8WVYKkJogDOLWigDZY8OnqJ0fy0ojMTSZNkxpjFYiksLLziiivMxMPevXvLspyZmSkIAujk5i9TE7NKTnu0zTHa6FUd23xIbAcPgLqu4zovEYgpzCk0E8mRZMaYLMmccZvNdl6/88zPMzMzHQ5HKBTKycmpTY5v3pBtHW15CuWn3qkb1cmnFvUuA2N8MoGO5rvcwZVAhNQ6+fiQHEkmhHDMNU2LxqI2xQY5+nE9jhBSFEUQBBDjVCyhKbQR6U3hl4tkuhYwxjbFhhBilFFGgY2JMQY1dOY+qG7987OZm6nx3UaAG6v1b77dhI8zR705YyOFekiODQMhR9PVzDhr9AXwBgR0KaSQQlKQVJ1sJjAxgppwjQiCkJLkMxBNqdamPNK/dFVskky1mss9Fbg702FyiTQccI0OQfAqC6KA6mqbG3qnzMAV5As1ZYdzzoGWqJ4tnbi/KRIQyEhMmDHPfjLWeFN65djHMR9aYuqheW2MMcMwZFmurc9tLF0v6eKdHElOfAE/u1sKbQ0wIjmqPyibeqGUUcxw7QauretKBEg7EYiqqhaLhVEG3HroaFHknFNKJUmC/esZcWYCPwSuoQoVkkZMDmAEuUZH58+jZoy0xPmr3qKv4b03tVSEs6M6IxSR2hU+Z9xMrKj9NmGua6awNOcuEpHSyWc6YOQBxTQy80YJQk2MJEEUMMOQVpmWloYxbhj71Q2dc66Iis1mU1WVECIgIbGmDQ7OWW1fziPycPSIF4gAlN2JV4sQ0nUd6rEQQrFYjFIqidIxrvlnbx+RhITZZhzHDKollnObdXWMM5OFhnOOSJ2oN/v4J4Ak+65rtxrMQCm0WWCMQRsnEVDKBglChmEwxgTSeC33UbxLDUYKzBpQ4QiEM7AnWOxQ5mG326Fq6oSv9gRSZc0qd9PmBxoZkxDKZHE94as6XqR0cgo/j0TjMx6PM8aAbSsWixFCGqbHC0SQZdnr9WZlZWGMJbF+4texkch0DfWe6GhJBuqVNFsabBNCmpopmo/EwrtjXVud9Q40FYmfN4cInXFmmuK1P2TJSYNJWiYNbgLJOn4KLQcgPxREgQgEEyyIQi0dQ90bNIuEzeUlqFxZkkVB1A1dEAUtrhnUEEQhGArG9Xg8Hk9PT4fB7Q/4w5EwHNAwDF3XwSCPxmqLJaDnBuNHmgqArQuns9vt4EMyk/YZY/F4PK7HTU50QRRgddrMUWcu1xN35pxDR65gMIgxDgQCJtMbQigajWqaVk/NwhF27txZXV3NOQ8EAmD5Y4xDoZCmaeAIMKhRT1zNqnvGWSwW03XdFJlYLBYOh+FpQ4W22Q/oWC/xZ+85hdMefr//iy++WLt27Y8//vjNN998XYeNGzeuX7/+888/ByofGNOiKMqSTAiZO3fuwUMHBVGwWq0333wzVDWuW7du+fLlgiAEg0Fd1wkmkHjvdDrnz58/7tpxQC0ajUWrq6vfeuutJUuWAAkBNShwdx1hlWC1jqhQKORyuYBZBRxdgiBMnTpV0zSfzwf1G4sWLXruL8+ZWvEEyAnAMPZ6vfMXzBdFccOGDf/85z+dTudDDz1UWVnJOR86dKjX6zUPbjrkKisrFy9evGDBAoSQzWZbtWrVli1bdu7cuXPnzl27dv33v/9dvnz5J598Ul1dXdumizPzH2ecYDJhwoTvvvsOXoSqqjabzeFwRKPRQCAQiUaAvdicyJpCyro+04ExDofDzz///NVXXy1JUmJhkyAIXq+3uLj40ksvbfjDPXv27N+//9FHH3W73V988QV8+Oabbw4ZMkQUxczMTJB8v99vsVhEUfzhhx+GDR82ffp0aHkhy/JXX33Vs2dPIKzfu3dvKBTq2bNnu9x25imoQQVRGD58+OOPP37ttddSg4K3HCHUq1ev559//q9//StCKBqNPvnkk3/6059OjO7fNNqBRmbhwoXjxo2jlFZWVpaVlVVVVQHX2pYtW6AZpQmo6/R4PGlpaUDYtHbt2oceeqhz587mPi6XS1GUtWvXXnDBBU2lRWGMYVIAHv9wOAxlocd1FylJTgLAtYOS0XTqlKCgoGDfvn033XRTeno6NKwAaJq2bNmyjRs3Wq1WsBhhIM54ZYbD4cjPz3/uuedgPdy9e/fXXnttx44dGzZsaN++/YIFCyorKy+55JK+ffsCvVFxcfHcuXO3b9++aNEiSmksFotEIqqqVlVVFRcXi6K4f/9+xpggCFDnCKV1jLIdO3bs2bPn2muvDYfDJssypfS66657+umnKysrPR7PK6+8MmrUqFtuucW88hOrwKmurrbb7SNHjlyyZEnXrl09Hs+cOXOuvvpqh8Px2muvzZw5MxwOJ9ZsmgTdlFKYAW02W0FBwccff5x4WJ/P17dv37y8PBDXhsAYT5o0qbq6GniUS0pK+vbt+9gfH7vq6qtgmYMQAo/aMe6oRST5tFwer1+/fvbs2Vu2bKGUjhkz5p577nE6nZqmuVwuWMZgjMETY7VaYdAXFxcvXrz4scceg6bEp/oOakENWlVdlZOTAz5hoKQNhUJfffUVsETIslxVVZWdnW0Yxr59+xBCqqpCXTGQSfztb3+bNGmSy+WaNm1aOBx2Op1AUu/xeCZOnEgpjcfjc+fOhTardrs9FostX758zJgxOTk5kyZNgnatGOM33njjs08/C4aCQNapKArMF/AYOeehcGj27Nnjxo1DCAEJwf79+8FAyMjIUFX1o48+2rt3b8eOHb1e7wcffKCq6uDBgz/44APwmcH9+v1+t9vdlDKEmzIMwyJb7Hb77NmzA4FAPB5fvny5t8J76YhLly9fHggEysvL582bN336dErpwYMHu3btetNNNz3xxBMLFy789ttvf/rpJ0EQVn+z+neTf2exWNauXQvU/9ChGg4IzBmUUxgqsVhMlmXJKlVUVPh8vg8//BD6adlsNliiu91uWHKbSTvHThdL6eSfRyAQWLFixYwZMyiluq57PJ7333/f6/U+9dRTaWlpgUDA5XJpmhaPx9PS0jRNY5SVV5S/+uqr77333q233gp641TfRC0YY6FwKCMjw+z2ABuFhYUY45EjR0YiEUmSsrKy9u7d27179x07dnTs2PHw4cOU0q5du8qynJ2dbbFY7rzzzkAgsGjRIoQQxjgvL2/37t1paWl+vz87O/v666//+uuvg8GgJEmhUMgwjKlTp/75z39GCBFCHnjggaFDh27cuPHCCy+c8coMi8WyZcsWRVFefvllmEFEQRQsQigUstlss2fPXrt2LVw8pbR9+/Y2m624uLi6utpqtW7dunX69OmzZs2CJ7x06dLXX38dmS5ighBCLpeLGrTRRCuQDVmWBSL4A/5AIHDhhRdaLJYlS5Zs27atsLDwhhtucLvdGzZsmDp16uLFi3v06GG32zt37vzDDz8oiqJp2tixYydOnDht2rRgMPjggw8eOnSIc/79998bhnHw4ME9e/b06dNn4MCBf/zjH3VdN538oiiCGkAIpaWldezYcf/+/Wf1PIsyGgqFoHaQAxMraS4RUkqSfwbUoLIsv/3224IgwCgvLi4uKChYsWLFRRddNG7cuFqyMYw559FolFK6YsWKu+66S5KkWCzmcrmOLcatGXKsqanJzMyEfuggYwgh6JMYiUTGjx+/d+/ec889t6ysLCsrq3PnzpqmdejQASF01113jR49un379oZhWCyWxx9/3OPxfPjhh2VlZffccw/nvKysbNiwYf/73/9UVX3iiScmTJhw8cUXd+3ataKiIi0tbeLEiUVFRQcOHKiqqtI0bf369TNmzLh74t0xNRYMBt1u965du+666y54huZShVL65Zdf9uvX76yzztJ1PRKJwDKyrKxsxowZLpdr27ZtkUiEMTZjxgzOucfjAdbHRICqT8wtTQRlVNM06CaZmZlpsVg6depUVla2aNGi3//+9//+979vu+22oUOHbt++Xdd1SimsXRljkUhElmXDMDweT2lpKXAGw9J64MCBkydPVlV19+7d27Ztu+KKK2DG9Hq9NpsNMYQQMgyjS5cuoVAoFov16tWroqLi448/hogdGAilJaU61Ruu1I6xamhVSW61bPIkwuf32Ww2sK8yMjI0TYMOGwUFBd99993w4cM9Hg/nHMxFr9d7xx13rFy50u12a5rmcDg2bdo0depU1Fj7OEBrSrLdbn/ggQdMNmkYYZqmhcNhh8Nx//339+nT51//+ld2dnZVVZUsy5RSt9u9evXqb7755v777zepc2+77TZYHl980cVnn302waRfv36qqnbp3OXLZV+Wl5dbrdbx48drmpaTk7Nx48ZYLHbVVVdVV1dnZWVBMvbDDz9cXl6em5tbWFi4bt06SumAAQMQQunp6SYxq81m+/vf/37NNdfAUwIx1jQNjIJIJJKXl1daWrp27dqRI0fCrGR20uKcM8T27Nrz+eefG4ZhtVqbes4ZGRkVFRXdunUbM2aMpmkfffTRs88+u2TJknXr1l166aVjx46dM2fOwYMHX3rppeeff37BggW7du0aMmRIu3btMMbQCq9Dhw7gVV6xYoXVav3+++/BCD906NCePXu2bNkiy7LFYrHZbPfddx+lVJZkxpmu6wcOHEB1VFkQHqeMyrKcn58viALQaTUfrSTJkCv/S2RagxY5VqsVPK7xeByCouXl5d27dwc/B5QBEEw45+3bt8/NzaWUapoG9P2KojidThhqDdGakmxS88C8A8FSu91usVhWrFixa9cuTdOqKqt27NwB+m3YsGG6rnfq1GnhwoUIIbj3WCwGvBE5OTn/Xvzvhe8tdLvdy5cvLyoqysjIYIzl5OQAcxNjjBq0Z8+es2fPXrp0abt27UB6VVWdP38+Qsjv98uybLPZli1b9s4773i9XnB3IYRUVfV6vXv37r3zzjvhT6Dvc7vdkiRNmjQpGAza7fZ9+/b5fL477rgDXsSmTZuefPJJjDHcaUVFRWVlJUhyRkZGo88ENOrWrVsvvfRSRVE6d+48d+7cwg6F11577fLly51O5x/+8If+/fvn5eUVFxcrivLdd98VFRXBG5dlWVEUr9f73bff7dq9a+TIkR07duzQoUP//v0lSaqsrFy3bt0111xjrjJqS0E48/v98BwYY9nZ2ZAxFolG4MFCLK0e63hiFm0iTNXYgqJFDQrZvPB/Q4anXwRAcem67nK5QqGQJEnl5eUdOnQoLS0dOXIkfAuGNKU0Kytr9uzZ69evf+SRR2pqasLhcJ8+fe6//36EUFOzWCI/bksjGovCBoQ94OKj0SjGOC0tbejQoaNGjXrooYcefvhh6Ip+7rnn7tmzB1o0WiwWRVFCoRDElp1O53333Tdx4sRoNOp0OgcMGLB48WIQJ6/XCyFQp9NZWVmZkZ7hyaq1ez1ZHsaYJElz5szZu3dv165do9Eo7KzGVHAOgXPY5XLNnj17wIABFoslGAw6nU5VVXNycn788cf9+/dff/31sVhMkqTq6urdu3ffddddwOIYCoVyc3N9Pp/L5SKEDBo0qG/fviAeIC0Nnwkke2CMbTZbdXX1v/71r+rq6rKyMkrp448/XlJSwhibM2dORkZGZmbm+vXr5/3/9q49tqmyjb89pz09Lb1fxtptQB3RTRA3hbREyIKrs5IozkCMCVIDWZyaD6lREyNjH2i8EAJkaOYfijINf5DJFlEjGgkYYrxFEhgOgY2LXe3Fdu3anp625/L98czzHdu1DEKBb9/5/dWdnct7Ls/zvs/z/J7n6et777334GNWKpXvvvtuT0/P4sWLV61a1dnZ+cknn7S2trrd7tHRUYIgHPMcjY2NyWQSNBSQ3tLptNVqPXv2LM/zsB2agQnhD7vdzjLs1ar4Skkyy7BUhqJpOhQKMQyjVqttNhs0/q7QFSsEIOV1dXX19PSEQiGapm+//fZgMOh2u1tbW2VC/WqZDMdxpVIZCAScTueRI0d279598OBBlUqVSCTKmco38HkIcUuEUCKRsFgs2Wz24Ycfrq2ttVgsfr9/bGxs9+7dYJGCW37fvn1A87h8+fKWLVsMBgNJkizL1tTUIISi0WhdXV08Ho/FYvPnzwc7nGVZvV5/6uQpcPWBHIKkQeO+YDBI0/SiRYtOnz5dV1fHcdyZM2cKhprNZr/++utNmzaJN05MTDQ0NEC7ScDRo0f7+vr27t3L83w6ndZoNEKYAGY8jUaTSCQQQkqlckrSKEEQ4+Pj0IjLbDY///zzdrv9gQce+Pbbb3U63RdffNHb21tTUxMOh1999dVNmzZBZxyo481xXEdHh8fjGRwchOXY999/73a7Ozo6WJYdGRkZHh4GP4LD4QCFLkCn03k8nqGhodmzZ4fD4dra2kwmwzAM3MJYYGzOnDmlolZTolKSLMNkNE1Di+B8Pl89u/qOhjuamprEGWRlbOZbZwKHkSxfvpxhmL179547d46iqLVr165du1atVgMJGcMwWAilUqnZVbMpilKr1T6fz+v1fvnllxaLhaIoSNwpxo1cXUMfXPgN8yfLsi+//LLRaAyHwxqNZuvWrcuWLXO73dlstrq6ur29HSwImUzmdDoNBkM0GjWbzZgMoygKOoekUqlYLNbZ2TkwMCDH5bgcpyjK4XDkmbzJONnZL5/PQ5g6n8+TJDl37tznnnuOpunHH38cBvPjjz/29PSA4zqfz1MU5ff7w+HwQw89BHM7rKU///zz33//vaqqKhgMqlQquVx+9OjRP/74o7u722azxWIxo9EIRE6Xy+V0OmHakMlkOp2uFEFKgSvELS8cDgf0J+Q4LpVK7dixY8uWLePj4waDoaqqanh4+NFHH62urk4kEmD6sgxbV1cHxiNCKB6Pz5kz57PPPguMBcDNhhD67rvvduzYoVKpNBoNlOmPRCIDAwMul4tl2DyTB0OMoijQj4ODg+BoFPqBCplqZVCpTujQ3qm7u9vj8TQ3N/v9/g8//LC5uRnKmkOYe8pjwWy41aZuhULR1tbW1tZWsB3H/nEjGo2G53noQiiXy61Wq9frhQizOPVPnHxzI12AYtaHWqXO5/NyuXzFihWQABwIBCiKeuWVV4R9wuFwS0sLMJnA7oXpDpfj0OYaIQQOW6PRCKoBZmDgXbIcC6F1hBDEosAMuXTp0saNG0dHR7Va7eXLl5csWRIIBEKhEESV0+m0xWLZsGGD1+tFCOEYzvO8VquVyWRWqzWRSAQCAehHz7JsY2MjGM/gj6BpGrx00MSbQxz6m7Zd6plwLGe32cGtDWLz+uuvHzhw4IUXXvj555/XrFlz//33x+PxP//80+v17tmz56OPPlq6dOmyZcsQQkBWxxEejUbBu7lw4UKzySyTyYCLHolEoMVKLpeDjl8ajSaXy/X29ppMpqampt7eXq1Wu379+uPHj7/22muHDh1SKpVvvvnm+fPnX3zxRXi8ClKBGISulLNVWRdUOp1evHhxIpGYP39+OBwGlQzhcvFuYp/26OgomJT/E45ueJGl/isev7ixe6l9bgym5CQHAoFnnnmmvb1dvFFoBwPxWNgIyQxCrgKGYWNjY6lUiuO4n3766auvvlKpVAaDQVz+niAIi8Uil8shT8DlcnV1del0OqVSSVGUQqE4cuTIp59+GovFrFYrSZIQv9m5cycSFSTheb6lpWXJkiVg78CBv/7668jIyMZ/bYR9oPYAQRBXtS4VAD2u8/m80Wg8c+ZMe3v7gQMHVqxYEQwG33nnna6urpaWlubmZrfbvW3bto6ODng46XS6oaGBZdmtW7euXr1ahsmAfQ03mMvlkskkSZKRSEQmk1kslqGhoR9++GHfvn0kSTY2Nh47dgzDsOXLl5tMpmPHjt1111379+9vbW198skn9Xo9SZKgHK9Y+72CHC+apuPxOEQgx8fHwV00MTFhsVhyuRyO4QcHDp49exYq2guIRCIcx8HoKzE2CZO6A0MYh8GstWbNmhMnThw+fLjAnl+wYAHIPJWhZs2aJWynKOqRRx6B8xAEQRBEbW0tcGZqamowDNu8ebP4PBcuXLh06ZLZbD58+HB/f//Jkyd9Ph+GYelU2mgyEgRx7ty5CxcukCS5fv36BQsW7N+/3+VywQqzYNgqlapYSgUtgyMcJL+44ft0gGHYG2+8cejQIZfL1d/fb7fb33777ba2tlgstmfPnoULF2q1Wr1e/8033zzxxBN+v3/z5s1DQ0Onh0739vYCg8Dj8Zw4cSIYDL7//vs4joNbbmxsjCRJq9WazWb9fv+uXbu8Xq/JaIpGozqdLhQKJZNJpVLZ0NAwPDzsdDrtdvvKlSt37dr11ltvCd4+MR9+SlRwToYAOkmS4Ca12WyCzTNLPUuGTYYKxA4MhJDZbF60aBEos1uth/BMxfbt2202m9/vhz8FGVi1ahWsqwvcFjiO+3w+hFAmk8Fx/LbbbmtubuY4zuFwPPXUU7A/y/7X++pwOCCQQxBEZ2fn008/nU6nQ6FQfX09x3HZbBbHcaBY1NfXI4T6+/sHBwchbickEkAFdblcDl+FilSxHDt37tzH2h+DwkCTNU94nmO5PJMXkpanWUiA4zkuz61bt87j8dxzzz0IoVwul0qltm/f3tTUJIyB47g7G+8cGho6deqUXC7X6/UnT5384IMPqqurN2zYkMlkrFbrzp07V69eDQtPlmFHRkd+++036DhbVVXl8/nuvffeZDJpMpnq6+vvu+8+rVY7MTGxbt26ZDIJYfOXXnrp4sWL2WwWxFghV5TKvhBQwZq1LMN2/7v72WefNZvNf/31V19fn8/nAy4+lJuZ8iiFQgGhCNBDt/4C+38Rk1l5fwfAZJhsfHzcbDanUikVqRI+/YsXL86bNw/kRMxsEZJ7gSCN/lkQ+9pemUKh8Pv9ZrP5l19++fjjj7dt26bT6cSGPcuwkEctHgMSKnIVBfPASzQdMRbyhIH1rVBMik2BfSekagP/BLwh4JNXq9WpVEqj0cTjcWEtKR4S1AMSdCLDMDiOT4adWA46YCKEQPVk6Aywd3ieVyqVmUxGSSgF/1kpVLb69ODg4Pnz51euXHn8+PG7777b6XRCZAKCe1O64+CehZargpkkifR1RIEkczwHnCqwY8X1N1iOxTG8QB6E+UGQ2UrZ5gAAAkZJREFUAcHLiq5UdqO4NAdskcvlkF8JP8BNKKasC0UOpridEjH5ac7GcLggKpOFwRQKEGnxXSOEhIh3IpFQKBSgayiKAoPWYDDA7+JRcTwHAXCEkNlsFm4EXKHixyJIMkEQkEoB0ZDyd1Ep3zWUR126dKlGowkEAg8++KDNZuM4DshDQFItPhCyT0mSzGQyhIKQmgNVCJOTJy5DCEEKu1qtpiiKIAiBwQIvCMfwKd+CYGxfL2aLkP0Dy2ng1RV4RlEphV7iM5mm9ode1jRNgzNfiBdOqTX0ej0EqPV6PUVR5a4rGhU8JYjhwQ+BcwZVSoTf4vNkMhmDwRCPxwvaX0590UrMyRzHQZ0HIRsO4vWQjY1hGE3TsHIoOJDn+Fw+B2wKvV4vHps0J19fCM8WEj80Gk00Gi1IvSxVq7ngJJNlt0p8joWHiApEiwGecJANmH94no9EIgJ/s8yHWupf0/9m4AxQSJBhGAijFtDyZH8XshZWxUCaQEWrayEJseAqDMPA/hAMF19a/FjgdbAMy7AMy7IqlQqYzuVvpyKSPM0AUrFjvXwutYTriyn5OcXas/gLEXtfpv/9lDnb1YpcJVAq57H8lilRPEhe1Im6+F8Fxn+p3cpf+mZWGijeTRLjWwfX/BLLQ+ixcs1nuIm45qGW11bTOe0V97mZ/damVM83ZSQSrhlX+/0Uz8z/Dy+9/Ixa7JC/BkidEyVImAmQnMMSJMwESJIsQcJMgCTJEiTMBEiSLEHCTIAkyRIkzARIkixBwkzAfwAM2tA8n8X1uQAAAABJRU5ErkJggg=="
    }
   },
   "cell_type": "markdown",
   "id": "3877a75f-816a-478c-a21a-26bba19acb70",
   "metadata": {},
   "source": [
    "![image.png](attachment:a1156988-5cc8-4673-a839-ca766b28c527.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63a70e95-0233-4f35-aa8f-5ea4f7ae792f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 查看箱型图与散点图\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import warnings \n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "path = \"D:/python/机器学习与社会科学应用/演示数据/08特征工程入门与实践/\"\n",
    "\n",
    "df = pd.read_csv(path+'boston.csv')\n",
    "df.head()\n",
    "# 箱型图\n",
    "df['INDUS'].plot.box()\n",
    "# 散点图\n",
    "df.plot.scatter(x='INDUS', y='TAX')\n",
    "# 可以看到异常值分布"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1092265f-3da8-4e44-991f-c7203f745817",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 将位于上极限与下极限之外的点定义为异常点，并标记为空值\n",
    "def detect_outliers(filename, attribute, fill_value=np.nan):\n",
    "    \"\"\"功能:探测异常值:\n",
    "    Args:\n",
    "        filename: 为数据集,\n",
    "        attribute: 为要识别的特征,\n",
    "        fill_value: 为对异常值进行标记值\n",
    "    \"\"\" \n",
    "    q1 = filename[attribute].quantile(0.25)\n",
    "    q3 = filename[attribute].quantile(0.75)\n",
    "    iqr = q3 - q1\n",
    "    # 下界\n",
    "    outliers_ll = q1 - 1.5*iqr\n",
    "    # 上界\n",
    "    outliers_ul = q3 + 1.5*iqr\n",
    "    # 将异常值标记为fill_value\n",
    "    filename[attribute].where((filename[attribute]<outliers_ul) | (filename[attribute]>outliers_ll))\n",
    "    filename[filename[attribute]==fill_value]\n",
    "    return filename\n",
    "\n",
    "# 对数据集df的'INDUS'列的异常值标记为True\n",
    "df = detect_outliers(df, 'INDUS', True)\n",
    "# 显示异常值的行\n",
    "df[df['INDUS']==True]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "962d6ad3-6e95-4a59-b29b-8bc91ac33cf9",
   "metadata": {},
   "source": [
    "#### 异常值处理"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5500d745-a1fa-4afc-ba98-826d7668f9a2",
   "metadata": {},
   "source": [
    "#### 删除\n",
    "\n",
    "直接将含有异常值的记录删除，通常有两种策略：整条删除或成对删除。这种方法最简单易行，但缺点也不容忽视，一是在观测值很少的情况下，这种删除操作会造成样本量不足；二是直接删除可能会对变量的原有分布造成影响，从而导致统计模型不稳定。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a03ac16e-6d75-4bd5-9e7e-638aa52adc09",
   "metadata": {},
   "outputs": [],
   "source": [
    "def detect_outliers(filename, attribute, threshold=3, fill_value=np.nan):\n",
    "    \"\"\"功能:探测异常值:\n",
    "    Args:\n",
    "        filename: 为数据集,\n",
    "        attribute: 为要识别的特征,\n",
    "        threshold: 为要设定的门槛值,默认为3倍标准差,\n",
    "        fill_value: 为对异常值进行标记值\n",
    "    \"\"\" \n",
    "    outliers_ll = filename[attribute].mean() - threshold*filename[attribute].std()\n",
    "    outliers_ul = filename[attribute].mean() + threshold*filename[attribute].std()\n",
    "    # 将异常值标记为\n",
    "    filename[attribute].where((filename[attribute]<outliers_ul) | (filename[attribute]>outliers_ll), fill_value)\n",
    "    return filename\n",
    "\n",
    "# 将数据集df的\"INDUS\"列的异常值删除\n",
    "df = detect_outliers(df, 'INDUS')\n",
    "df.dropna(subset=['INDUS'], inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1eb9664-adde-4429-b23c-a8e0788034b9",
   "metadata": {},
   "source": [
    "####  视为缺失值\n",
    "\n",
    "利用处理缺失值的方法来处理。这一方法的好处是能够利用现有变量的信息，来填补异常值。需要注意的是，将该异常值视为缺失值处理，需要根据该异常值（缺失值）的特点来进行，针对该异常值（缺失值）是完全随机缺失的、随机确实还是非随机缺失的不同情况进行不同的处理。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01336fa5-e644-46f7-b201-d5d3c864c293",
   "metadata": {},
   "outputs": [],
   "source": [
    "def detect_outliers(filename, attribute, threshold=3, fill_value=np.nan):\n",
    "    \"\"\"功能:探测异常值:\n",
    "    Args:\n",
    "        filename: 为数据集,\n",
    "        attribute: 为要识别的特征,\n",
    "        threshold: 为要设定的门槛值,默认为3倍标准差,\n",
    "        fill_value: 为对异常值进行标记值\n",
    "    \"\"\" \n",
    "    outliers_ll = filename[attribute].mean() - threshold*filename[attribute].std()\n",
    "    outliers_ul = filename[attribute].mean() + threshold*filename[attribute].std()\n",
    "    # 将异常值标记为\n",
    "    filename[attribute].where((filename[attribute]<outliers_ul) | (filename[attribute]>outliers_ll), fill_value)\n",
    "    return filename\n",
    "\n",
    "# 将数据集df的\"INDUS\"列的异常值表示为缺失值\n",
    "df = detect_outliers(df, 'INDUS')\n",
    "df[df['INDUS']==np.nan]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "671b2da2-1aea-464d-8063-eaa573ca7e39",
   "metadata": {},
   "source": [
    "#### 缩尾或盖帽法\n",
    "\n",
    "在经济学中，这种方法被称为“缩尾”，在数据分析领域，该方法也被称为“盖帽法”。该方式是将整列替换数据框99%以上和1%以下的点，将99%以上的点值=99%的点值；小于1%的点值=1%的点值。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e294bdb-b595-423d-bd86-26658c10e5bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "def winsor(x, quantile=[0.01, 0.99]):\n",
    "    \"\"\"缩尾法处于异常值\n",
    "    Args:\n",
    "        x: pd.Series列,连续变量\n",
    "        quantile: 指定缩尾的上下分位数范围\n",
    "    \"\"\"\n",
    "    # 生成分位数\n",
    "    Q01, Q99 = x.quantile(quantile).values.tolist()\n",
    "    print(Q01, Q99)\n",
    "    # 替换异常值为指定的分位数\n",
    "    if Q01 > x.min():\n",
    "        x = x.copy()\n",
    "        x.loc[x<Q01] = Q01\n",
    "    if Q99 < x.max():\n",
    "        x = x.copy()\n",
    "        x.loc[x>Q99] = Q99\n",
    "    return x\n",
    "\n",
    "sample = pd.DataFrame({'normal': np.random.randn(1000)})\n",
    "sample.hist(bins=50)\n",
    "winsor_sample = sample.apply(winsor)\n",
    "winsor_sample.hist(bins=50)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "748f9889-a49b-4fe3-9591-50392c5e7465",
   "metadata": {},
   "source": [
    "#### 分箱法\n",
    "\n",
    "分箱法通过考察数据的“近邻”来光滑有序数据的值。有序值分布到一些桶中或箱中。包括等深分箱：每个分箱的样本量一致；等宽分箱：每个分箱中的范围一致。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef2c0c41-ddc3-4cc4-8742-7d13e62a6c7c",
   "metadata": {},
   "source": [
    "####  不处理\n",
    "\n",
    "根据该异常值的性质和特点，使用更加稳健模型来修饰，然后直接在该数据集上进行数据挖掘。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60c53f02-e9c6-4c59-a921-8334c3da9b34",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 3.3 类别变量编码"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "329a530e-aa43-4316-8790-395bb05ef14c",
   "metadata": {},
   "source": [
    "特征的数据类型可分为两大类，即连续型和离散型。其中离散型特征既有数值型的，也有类别型的。例如，性别（男、女）、成绩等级（A、B、C）等。连续型特征的原始形态就可以作为模型的输入，无论是LR、神经网络，还是SVM、GBDT、xgboost等。但是除了决策树等少数模型等直接处理字符串形式的类别型特征外，逻辑回归、支持向量机等模型的输入必须是数值型特征才能进行训练。因此，就需要将离散型中的类别型特征转换成数值型的。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d8a093d-ba8e-411a-b6eb-ed121b669a37",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### 序号编码（Ordinal Encoding）\n",
    "\n",
    "**序号编码通常用于处理类别间具有内在大小顺序关系的数据**，对于一个具有m个类别的特征，我们将其对应地映射到[0, m-1]的整数。例如对于”学历”这样的类别，”学士”、”硕士”、”博士” 可以很自然地编码成 [0,2]，因为它们内在就含有这样的逻辑顺序。但如果对于“颜色”这样的类别，“蓝色”、“绿色”、“红色”分别编码成[0,2]是不合理的，因为我们并没有理由认为“蓝色”和“绿色”的差距比“蓝色”和“红色”的差距对于特征的影响是不同的。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ee26f74-80b1-4683-8c4e-cf170a27aad2",
   "metadata": {},
   "source": [
    "sklearn.preprocessing中的OrdinalEncoder进行处理，可以通过两种方式进行编码：\n",
    "- 方法1：先将数据categorical_df载入encoder中，再将载入的数据进行OrdinalEncoder编码\n",
    "- 方法2：直接载入encoder并编码"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "859aa90e-039f-471d-8f27-a77d99966fd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd \n",
    "from sklearn.preprocessing import OrdinalEncoder \n",
    "\n",
    "categorical_df = pd.DataFrame({'my_id': ['101', '102', '103', '104'],\n",
    "                              'name': ['allen',  'chartten', 'bob','dory'],\n",
    "                              'place': ['third', 'second', 'first', 'second']})\n",
    "\n",
    "print(categorical_df)\n",
    "print('-'*30)\n",
    "\n",
    "encoder = OrdinalEncoder()  # 创建OrdinalEncoder对象\n",
    "\n",
    "# 方法1\n",
    "encoder.fit(categorical_df)  # 将数据categorical_df载入encoder中进行转换\n",
    "categorical_df = encoder.transform(categorical_df)  \n",
    "\n",
    "# 方法2\n",
    "#categorical_df = encoder.fit_transform(categorical_df)  # 代替上面的fit()和transform()\n",
    "\n",
    "print(categorical_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "931b4f19-345f-4a82-b4f9-a48be9dc86c8",
   "metadata": {},
   "source": [
    "####  独热编码（One-hot Encoding）"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba4f7499-07eb-4737-9890-86c055e6d803",
   "metadata": {},
   "source": [
    "如果类别特征本身有顺序（例：优秀、良好、合格、不合格），那么可以保留单列自然数编码，但是如果类别特征没有明显的顺序（例：红、黄、蓝），则可以使用独热编码（One-Hot Encoding），即每个类作为1个特征。输出的矩阵是稀疏的，含有大量的0.独热编码通常用于处理类别间不具有大小关系的特征。可以通过导入sklearn.preprocessing中的OneHotEncoder，创建哑变量进行处理。对于独热编码，存在类别就生成几个特征。Pandas与sklearn都提供了独热编码方式。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "464b0019-442f-4d33-bc7d-eb8a695d9de4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# 方法一\n",
    "import pandas as pd \n",
    "\n",
    "df = pd.DataFrame({'f1': ['A', 'B', 'C'],\n",
    "                  'f2': ['Male', 'Female', 'Male']})\n",
    "\n",
    "print(df)\n",
    "\n",
    "df = pd.get_dummies(df, columns=['f1', 'f2'])\n",
    "print(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0717f907-14f3-49e9-b94b-d4317597ca45",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# 方法二\n",
    "from sklearn.preprocessing import OneHotEncoder\n",
    "\n",
    "enc = OneHotEncoder(handle_unknown='ignore')  # handle_unknown\n",
    "\n",
    "df = pd.DataFrame({'f1': ['A', 'B', 'C'],\n",
    "                  'f2': ['Male', 'Female', 'Male']})\n",
    "\n",
    "array = enc.fit_transform(df).toarray()\n",
    "print(array)\n",
    "df = pd.DataFrame(array)\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1589ea31-685c-43ab-881f-5ddb498fa240",
   "metadata": {},
   "source": [
    "#### 标签编码（Label Encoding）"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a5ede8a-515b-4cd0-847d-76256ee8187f",
   "metadata": {},
   "source": [
    "Label Encoding是给某一列数据编码，而Ordinal Encoding是给所有的特征编码。因此，Label Encoding常用于给标签（label）编码，而Ordinal Encoding常用于给数据集中的特征编码。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "023ff301-f255-452d-82d3-210bd9677c50",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd \n",
    "from sklearn.preprocessing import LabelEncoder \n",
    "\n",
    "categorical_df = pd.DataFrame({'my_id': ['101', '102', '103', '104', '102'],\n",
    "                              'name': ['allen', 'chartten', 'bob', 'dory', 'bob'],\n",
    "                              'place': ['third', 'second', 'first', 'second', 'third']})\n",
    "\n",
    "print(categorical_df)\n",
    "print('-'*30)\n",
    "\n",
    "encoder = LabelEncoder()\n",
    "categorical_df['place'] = encoder.fit_transform(categorical_df['place'])\n",
    "categorical_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c1cdd4f-b99f-40b9-b3c2-6f3dff28e371",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 3.4 特征缩放"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2031d5b7-27c6-450c-9c41-526120458196",
   "metadata": {},
   "source": [
    "有些特征的值是有界限的，比如经度和维度，但有些数值型特征可以无限制地增加，比如计数值。有些模型是输入的平滑函数，比如线性回归模型、逻辑回归模型或包含矩阵的模型，它们会受到输入尺度的影响。相反，那些基于树的模型则根本不在乎输入尺度有多大。如果模型对输入特征的尺度很敏感，就需要进行特征缩放。特征缩放会改变特征的尺度，有些人将其称为特征归一化。特征缩放通常对每个特征独立进行。\n",
    "\n",
    "对于数据尺度不相同, 可以选用某种归一化操作, 在机器学习流水线上处理该问题. 特征缩放操作旨在将行和列对齐并转化为一致的规则.例如, 归一化的一种常见形式是将所有定量列转化为同一个静态范围中的值(例如,所有数都位于0-1). 我们也可以使用数学规则,例如所有列的均值和标准差必须相同, 以便在同一个直方图上显示. 标准化通过确保所有行和列在机器学习中得到平等对待,让数据保持一致. 常见的标准化方法有: - z分数标准化 - min-max标准化 - 行归一化\n",
    "\n",
    "不论使用何种缩放方法，特征缩放总是将特征除以一个常数（称为归一化常数）。因此不会改变但特征分布的形状。\n",
    "\n",
    "当一组输入特征的尺度相差很大时，就需要进行特征缩放。例如，一个人气很高的商业网站的日访问量可能是几十万次，而实际购买可能只有几千次。如果两个特征都被模型所使用，那么模型就需要在确定如何使用它们时先平衡一下尺度。如果输入特征的尺度差别非常大，就会对模型的训练算法带来数值稳定性方面的问题。在这种情况下，就应该对特征进行标准化。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "161094d2-d17e-4419-b512-d37ad704b8d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 导入数据\n",
    "from sklearn.datasets import load_iris\n",
    " \n",
    "#导入IRIS数据集\n",
    "iris = load_iris()\n",
    "\n",
    "#特征矩阵\n",
    "iris.data\n",
    "\n",
    "#目标向量\n",
    "iris.target"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eeb8ea3b-4324-4824-a461-fab5c972a5c0",
   "metadata": {},
   "source": [
    "#### Z分数标准化"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66eb4614-b563-43d2-90a2-c3acd6eaf187",
   "metadata": {},
   "source": [
    "z分数标准化是最常见的标准化技术,利用统计学里简单的z分数(标准分数)思想. z分数标准化的输出会被重新缩放,使均值为0, 标准差为1. 通过缩放特征, 统一化均值和方差(标准差的平方), 可以让KNN这种模型达到最优化,而不会倾向于较大比例的特征.\n",
    "$$z=\\frac{x-\\mu}{\\sigma}$$\n",
    "其中z为新的值(z分数), x是原数据, $\\mu$是该列的均值, $\\sigma$是列的标准差"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "997a82c2-5fd5-4c50-9137-f3f2bc66d8ab",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import StandardScaler\n",
    " \n",
    "#标准化，返回值为标准化后的数据\n",
    "standard = StandardScaler()  # 对象实例化\n",
    "zscore = standard.fit_transform(iris.data)  # 调用实例方法\n",
    "zscore"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2cbe0349-65ef-4c26-833e-4b1874e9b1ee",
   "metadata": {},
   "source": [
    "此外，也可以手动计算"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe63f96c-85b9-4d01-ae9c-d2a43eb3804d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# 创建一个示例数据集\n",
    "data = np.array([[1.0, 2.0, 3.0],\n",
    "                 [4.0, 5.0, 6.0],\n",
    "                 [7.0, 8.0, 9.0]])\n",
    "\n",
    "# 计算特征的均值和标准差\n",
    "mean = np.mean(data, axis=0)\n",
    "std = np.std(data, axis=0)\n",
    "\n",
    "# 执行标准化操作\n",
    "standardized_data = (data - mean) / std\n",
    "\n",
    "print(\"原始数据：\")\n",
    "print(data)\n",
    "print(\"标准化后的数据：\")\n",
    "print(standardized_data)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40bc091d-53af-4b68-b562-1029b2a46a61",
   "metadata": {},
   "source": [
    "#### 区间缩放（min-max标准化）"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c88a5043-4db2-4581-85f0-de8c4f2800c5",
   "metadata": {
    "tags": []
   },
   "source": [
    "min-max标准化与z-score标准化类似, 它是将原数据缩放为[0,1]的区间.\n",
    "公式为:$m = \\frac{x-x_{min}}{x_{max}-x_{min}}$, 其中m为新的值,x为原数据, $x_{min}$为该列最小值, $x_{max}$为该列最大值."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d1a5cd6-2f37-415c-af23-9f727fbd5310",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import MinMaxScaler\n",
    "\n",
    "#区间缩放，返回值为缩放到[0, 1]区间的数据\n",
    "zoom = MinMaxScaler()\n",
    "minmax=zoom.fit_transform(iris.data)\n",
    "minmax"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d543f617-c9ee-4172-8cb5-9b129d544236",
   "metadata": {
    "tags": []
   },
   "source": [
    "此外，也可以手动计算"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d50b9f7-8fa6-4f18-8212-36b253b4821d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# 创建一个示例数据集\n",
    "data = np.array([[1.0, 2.0, 3.0],\n",
    "                 [4.0, 5.0, 6.0],\n",
    "                 [7.0, 8.0, 9.0]])\n",
    "\n",
    "# 计算特征的最小值和最大值\n",
    "min_val = np.min(data, axis=0)\n",
    "max_val = np.max(data, axis=0)\n",
    "\n",
    "# 执行归一化操作（将特征缩放到[0, 1]范围内）\n",
    "normalized_data = (data - min_val) / (max_val - min_val)\n",
    "\n",
    "print(\"原始数据：\")\n",
    "print(data)\n",
    "print(\"归一化后的数据：\")\n",
    "print(normalized_data)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "897bc508-e8cd-4e8d-ab11-ec6eb18af597",
   "metadata": {},
   "source": [
    "####  行归一化"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c869e4d-5e0a-48ce-921a-4f4fda19fc93",
   "metadata": {},
   "source": [
    "该标准化方法是关于行的,而不是关于列的.行归一化不是计算每列的统计值(均值,标准差等), 而是会保证每行有单位范数(unit norm), 意味着每行的向量长度相同. 如果每行数据都在一个n为空间内,那么每行都有一个向量范数(长度). 也就是说, 每行都是空间内的一个向量:$x=(x1, x2, ..., xn)$. 通俗地理解, 行归一化就是把每一行看成一个向量,然后转化为单位向量.\n",
    "\n",
    "L2范数:$||x||=\\sqrt{(x^2_1 + x^2_2 + ...+ x^2_n)}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b49f0d9-d4ce-40e9-9928-72970894e060",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import Normalizer\n",
    "\n",
    "#归一化，返回值为归一化后的数据\n",
    "Normalizer().fit_transform(iris.data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12bd6b50",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#连续特征离散化"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb9b3ecb-706a-456c-8f8f-f150d0d168fe",
   "metadata": {},
   "source": [
    "有时连续的特征对于模型不是必须的，比如只关心及格与不及格，而不关心具体的分数，或者只关心年龄段而不关心具体的年龄，这时就需要将连续的特征根据需要进行分箱。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f5aab01-45d6-4787-8db4-408d335e4246",
   "metadata": {},
   "source": [
    "### 对定量特征二值化"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e650772-42e4-4553-958e-a78f174e844b",
   "metadata": {},
   "source": [
    "定量特征二值化的核心在于设定一个阈值，大于阈值的赋值为1，小于等于阈值的赋值为0，公式表达如下：\n",
    "$$ \n",
    "x'=\\left\\{\n",
    "\\begin{matrix}\n",
    "1, x>threshold \\\\\n",
    "0, x<=threshold \n",
    "\\end{matrix}\n",
    "\\right.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6760177d-0671-4a43-b515-ccd8f2690900",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import Binarizer\n",
    "\n",
    "#二值化，阈值设置为3，返回值为二值化后的数据\n",
    "Binarizer(threshold=3).fit_transform(iris.data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95957aae-f8de-463a-82a6-783bd000975a",
   "metadata": {},
   "source": [
    "#### 将连续特征分箱"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac9bd4b3-3af2-434c-afc1-6b0e3660bfda",
   "metadata": {},
   "source": [
    "有时，如果数据是连续的，那么将其转换为分类变量可能也是有意义的。比如将人群中的年份分段。\n",
    "Pandas的cut()函数，可以将数据分箱（binning），亦称为分桶（bucketing），意思是创建数据的范围。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a52d9fc2-4757-4502-8aa6-f5bcf16e2563",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 默认的类别名就是分箱\n",
    "df = pd.DataFrame(iris.data)\n",
    "pd.cut(df[1], bins=3, labels=False)  # 将iris.data的第1列划分为3个等距段"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20c81fec-f3c6-4e59-b513-d1103f4d0fd7",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 第四节 特征构造：生成新数据"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1fa631fb-9241-4001-acf0-e6ab84d070af",
   "metadata": {},
   "source": [
    "从原始数据中构造新特征，在机器学习或统计学习中，又称为变量选择、属性选择或变量子集选择，是在模型构建中，选择相关特征并构成特征子集的过程。根据已有特征生成新特征，增加特征的非线性。常见的数据变换有基于多项式的、指数函数的、对数函数等。特征工程中引入的新特征，需要验证它确实能提高预测精度，而不是加入一个无用的特征增加算法运算的复杂度。构造新的特征要求你在样本数据上花费大量的时间并且思考问题的本质，数据的结构，以及怎么最好的在预测模型中利用它们。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce340181-83b7-461f-a670-d76f56cd31c5",
   "metadata": {},
   "source": [
    "### 4.1 数值特征的简单变换"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f49e276b-5bc2-4a1d-9248-30697c8a6c5e",
   "metadata": {},
   "source": [
    "- 单独特征乘以一个常数（constant multiplication）或者加减一个常数：对于创造新的有用特征毫无用处，只能作为对已有特征的处理。\n",
    "- 任何针对单独特征列的单调变换（如对数）：不适用于决策树算法。对于决策树而言，$X, X^3, X^5$之间没有差异，$|X|,X^2,X^4$之间没有差异，除非发生舍入误差。\n",
    "- 线性组合（linear combination）：仅适用于决策树以及基于决策树的ensemble（如gradient boosting, random forest），因为常见的axis-aligned split function不擅长捕获不同特征之间的相关性；不适用于SVM、线性回归、神经网络等。\n",
    "- 多项式特征（polynomial feature）：\n",
    "- 比例特征（ratio feature）：X_1/X_2\n",
    "- 绝对值（absolute value)|\n",
    "- max(X1, X2), min(X1, X2), X1orX2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54604b0f-bcfd-4c4d-be3b-44c89bee4d54",
   "metadata": {},
   "source": [
    "常见的数据变换有基于多项式的、基于指数函数的、基于对数函数的。\n",
    "\n",
    "多项式变换是一种常见的技术，用于创建原始特征的高次多项式特征。例如，度为2的多项式转换公式如下：如果由a和b两个特征，则它的2次多项式为（1,a,b,a^2,ab, b^2）这可以帮助模型捕捉特征之间的非线性关系，从而提高模型的性能。多项式变换的基本思想是将输入特征的组合作为新特征添加到数据集中，以更好地拟合数据。需要注意的是，多项式特征变换可能导致特征维度的急剧增加，特别是在较高次数的情况下。因此，在应用多项式特征变换时，要谨慎选择多项式的次数，以避免维度灾难。通常，你可以使用交叉验证等技术来确定最佳的多项式次数。\n",
    "多项式特征变换通常在数据预处理阶段用于改进模型性能，特别是在涉及非线性关系的问题中，如回归和分类。以下是多项式变换的示例代码："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ce991b5-1f8d-45f6-bab5-68aad579586f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "\n",
    "# 创建一个简单的示例数据集\n",
    "X = np.array([[1, 2],\n",
    "              [3, 4],\n",
    "              [5, 6]])\n",
    "\n",
    "# 创建PolynomialFeatures对象，将特征升高到2次多项式\n",
    "poly = PolynomialFeatures(degree=2)\n",
    "\n",
    "# 进行多项式特征变换\n",
    "X_poly = poly.fit_transform(X)\n",
    "\n",
    "# 打印变换后的数据\n",
    "print(\"Original Data:\")\n",
    "print(X)\n",
    "print(\"\\nPolynomial Transformed Data:\")\n",
    "print(X_poly)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e434b988-e7c0-4144-b235-2bd1bf309f3b",
   "metadata": {},
   "source": [
    "对数变换(logarithmic transformation)是一种常用的数值特征转换方法,主要应用在数据分布偏斜的场景。对数变换的主要作用有：一是，减少数据的偏斜性，对数变换可以压缩数据中极大值和极小值的距离,使得偏斜分布的数据变得更加对称。这有助于降低少数极端数据对模型的影响。二是使数据更加符合正态分布，许多算法对输入数据有正态分布的假设，对数变换可以使偏离正态分布的数据更加贴近正态分布。三是稳定方差，对数变换可以降低数据方差不稳定的问题，对数压缩后,数据分布两端的方差不同会得到一定程度的纠正。四是，提高模型解释力和预测力，通过上述作用,对数变换可以帮助提高模型的解释力和预测性能。\n",
    "对数变换的常用形式有自然对数变换、10为底对数变换等。需要根据具体问题和数据分布情况,选择合适的对数底。，使用preproccessing库的FunctionTransformer对数据进行对数函数转换的代码如下："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adfba417-2f04-402b-9fc0-62c487f2d476",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "# 创建一个示例数据集\n",
    "data = pd.DataFrame({'A': [1, 2, 3, 4, 5],\n",
    "                     'B': [10, 20, 30, 40, 50]})\n",
    "\n",
    "# 对特征A进行对数变换\n",
    "data['A_log'] = np.log(data['A'])\n",
    "\n",
    "print(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a0d88fb-ae32-4dcc-a232-81692010f8ed",
   "metadata": {},
   "source": [
    "### 4.2 类别特征与数值特征的组合"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fced4db-3270-4f2c-9771-84033e17f5ab",
   "metadata": {},
   "source": [
    "用N1和N2表示数值特征，用C1和C2表示类别特征，利用pandas的groupby操作，可以创造出以下几种有意义的新特征：\n",
    "\n",
    "median(N1)_by(C1)  \\\\ 中位数  \n",
    "mean(N1)_by(C1)  \\\\算数平均值  \n",
    "mode(N1)_by(C1)  \\\\众数\n",
    "min(N1)_by(C1)  \\\\最小值  \n",
    "max(N1)_by(C1)  \\\\最大值  \n",
    "std(N1)_by(C1)  \\\\标准差  \n",
    "var(N1)_by(C1)  \\\\方差  \n",
    "freq(C2)_by(C1)  \\\\频数"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "878bf041-3a68-4654-b3d2-feaf070fc746",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 第五节 特征选择：筛选属性"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "770a5e50-9d5d-4d62-a0be-e6473515c488",
   "metadata": {},
   "source": [
    "当特征构造完后，需要选择有意义的特征输入机器学习的模型进行训练。通常来说，从两个方面来选择特征：\n",
    "- 特征是否发散：如果一个特征不发散，例如方差接近于0，也就是说样本在这个特征上基本没有差异，那么这个特征对于样本的区分没有什么用。\n",
    "- 特征与目标的相关性：这点显而易见，与目标相关性高的特征，应当优先选择。除方差法外，本文介绍的其他方法均从相关性考虑。\n",
    "\n",
    "根据特征选择的形式又可以将特征选择方法分为3种：\n",
    "- Filter：过滤法，按照发散性或者相关性对各个特征进行评分，设定阈值或者待选择阈值的个数，选择特征。\n",
    "- Wrapper：包装法，根据目标函数（通常是预测效果评分），每次选择若干特征，或者排除若干特征。\n",
    "- Embedded：嵌入法，先使用某些机器学习的算法和模型进行训练，得到各个特征的权值系数，根据系数从大到小选择特征。类似Filter方法，但是是通过训练来确定特征的优劣。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "899383ff-093d-4189-aba9-e0c33f34abd6",
   "metadata": {},
   "source": [
    "### 5.1 Filter"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c329d12-7400-4877-8d61-691a591cfcf6",
   "metadata": {},
   "source": [
    "#### 方差选择法"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c193ad4c-fd90-4cbc-9bd9-8d6f46dfed9f",
   "metadata": {},
   "source": [
    "　使用方差选择法，先要计算各个特征的方差，然后根据阈值，选择方差大于阈值的特征。使用feature_selection库的VarianceThreshold类来选择特征的代码如下："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c36e70bd-ed51-4d8c-9894-4c8940111d38",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sklearn.feature_selection import VarianceThreshold\n",
    "\n",
    "# 创建方差选择器对象，设置阈值\n",
    "selector = VarianceThreshold(threshold=0.1)\n",
    "\n",
    "# 应用方差选择器进行特征选择\n",
    "X_selected = selector.fit_transform(iris.data)\n",
    "\n",
    "# 打印选择后的特征\n",
    "print(X_selected)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9899c2be-353e-4d75-bbd9-9a8f6142a1b2",
   "metadata": {},
   "source": [
    "####  相关系数法"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97253784-e61a-4410-bae9-2887be105260",
   "metadata": {},
   "source": [
    "使用相关系数法，先要计算各个特征对目标值的相关系数以及相关系数的P值。有两种方法：一是pandas提供了df.corr()计算相关系数；二是用feature_selection库的SelectKBest类结合相关系数。通过对相关系数排序，然后选择相关系数高的特征。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ddc9cfb-6908-40aa-832e-faff0f4a046f",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "# 方法一\n",
    "from sklearn.feature_selection import SelectKBest\n",
    "from sklearn.feature_selection import f_classif\n",
    "\n",
    "# 创建相关系数选择器对象\n",
    "selector = SelectKBest(score_func=f_classif, k=4)\n",
    "\n",
    "# 应用相关系数选择器进行特征选择\n",
    "X_selected = selector.fit_transform(iris.data, iris.target)\n",
    "\n",
    "# 打印选择后的特征\n",
    "print(X_selected)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c19b47a-fe0e-4bb2-afea-faeaa6a804d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 方法二\n",
    "import pandas as pd \n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "\n",
    "x = pd.DataFrame(np.random.rand(100, 8))\n",
    "print(x.corr())\n",
    "sns.heatmap(x.corr())\n",
    "\n",
    "# 假定第7列为响应变量，查看与响应变量最相关的特征\n",
    "x.corr()[7].sort_values(ascending=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5ad8ff1-83ec-4153-9375-5934cbf80aa5",
   "metadata": {},
   "source": [
    "卡方检验"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c043b62-91fd-4b9c-9eac-3947c3b5dbf8",
   "metadata": {},
   "source": [
    "卡方检验适用于分类问题，计算每个属性与目标便来给你之间的卡方统计量，选择卡方值较高的属性。经典的卡方检验是检验定性自变量对定性因变量的相关性。具体的操作代码如下："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "162f3f3a-ef55-4394-be76-b6b6a4487853",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sklearn.feature_selection import SelectKBest\n",
    "from sklearn.feature_selection import chi2\n",
    "\n",
    "# 创建卡方选择器对象\n",
    "selector = SelectKBest(score_func=chi2, k=3)\n",
    "\n",
    "# 应用卡方选择器进行特征选择\n",
    "X_selected = selector.fit_transform(iris.data, iris.target)\n",
    "\n",
    "# 打印选择后的特征\n",
    "print(X_selected)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6ac6ec4-e857-4a73-9c43-5763b3ab71ac",
   "metadata": {},
   "source": [
    "互信息\n",
    " \n",
    " 该方法是计算每个属性与目标变量之间的互信息，选择互信息较高的属性。互信息是一种用于度量两个随机变量之间的相关性的方法。它基于信息论的概念，衡量了两个变量之间的共享信息量。互信息可以用来衡量一个变量对另一个变量的贡献程度，或者说它们之间的相互依赖程度。它的数值范围为0到正无穷，数值越大表示两个变量之间的相关性越高。互信息的计算方法可以根据具体的情况选择不同的算法，常用的包括经验互信息和最大信息系数等。在特征选择过程中，互信息可以作为一种过滤法的评价指标，用于选择与目标变量相关性较高的特征。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54fd151a-af5e-4945-845a-d1a5c86b3182",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sklearn.feature_selection import SelectKBest\n",
    "from sklearn.feature_selection import mutual_info_classif\n",
    "\n",
    "# 创建互信息选择器对象\n",
    "selector = SelectKBest(score_func=mutual_info_classif, k=4)\n",
    "\n",
    "# 应用互信息选择器进行特征选择\n",
    "X_selected = selector.fit_transform(iris.data, iris.target)\n",
    "\n",
    "# 打印选择后的特征\n",
    "print(X_selected)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3946643b-dc45-49ff-8a54-92402dcb56b3",
   "metadata": {},
   "source": [
    "方差分析（Analysis of Variance, ANOVA)\n",
    "适用于分类问题，计算每个属性对目标变量的方差贡献，选择方差贡献较大的属性。相关的操作代码如下："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19b0d66e-0e09-4d6d-94bd-6459d258f937",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sklearn.feature_selection import SelectKBest\n",
    "from sklearn.feature_selection import f_classif\n",
    "\n",
    "# 创建方差分析选择器对象\n",
    "selector = SelectKBest(score_func=f_classif, k=2)\n",
    "\n",
    "# 应用方差分析选择器进行特征选择\n",
    "X_selected = selector.fit_transform(iris.data, iris.target)\n",
    "\n",
    "# 打印选择后的特征\n",
    "print(X_selected)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "faee05ff-b179-44db-a8a1-ebc7fb1b6d56",
   "metadata": {},
   "source": [
    "### 5.2 Wrapper\n",
    "\n",
    "包装法是一种特征选择方法，它通过构建不同的特征子集，然后使用一个特定的机器学习模型进行评估，来确定哪些特征对模型性能的提高最为显著。包装法通常比过滤法和嵌入法更加计算密集，因为它需要多次训练模型，但它可以更准确地选择最佳特征子集。包装法的一个常见类型是递归特征消除（Recursive Feature Elimination，RFE）。以下是一些包装法中常用的筛选属性的方法。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef0ea0e3-14e0-4af2-a9ff-8de74dd7328b",
   "metadata": {},
   "source": [
    "递归消除特征法使用一个基模型来进行多轮训练，每轮训练后，消除若干权值系数的特征，再基于新的特征集进行下一轮训练。使用feature_selection库的RFE类来选择特征的代码如下："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "736526b8-d2ba-47b0-9589-a1c50bd3e70a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sklearn.feature_selection import RFE\n",
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "# 创建线性回归模型\n",
    "model = LinearRegression()\n",
    "\n",
    "# 创建RFE对象，指定要选择的特征数量\n",
    "rfe = RFE(estimator=model, n_features_to_select=2)  \n",
    "\n",
    "# 拟合RFE模型\n",
    "fit = rfe.fit(iris.data, iris.target)\n",
    "\n",
    "# 打印所选特征的排名\n",
    "print(\"Num Features: %s\" % (fit.n_features_))\n",
    "print(\"Selected Features: %s\" % (fit.support_))\n",
    "print(\"Feature Ranking: %s\" % (fit.ranking_))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d70e28f-8c1e-47d7-800a-262a6daaf2a2",
   "metadata": {},
   "source": [
    "前向选择（Forward Selection）：\n",
    "前向选择从一个空特征集开始，然后逐步添加一个特征，每次选择对模型性能提升最大的特征，直到达到预定的特征数量或性能指标。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52f5870f-d848-4f0c-a25b-72095142e106",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sklearn.feature_selection import SequentialFeatureSelector\n",
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "# 创建线性回归模型\n",
    "model = LinearRegression()\n",
    "\n",
    "# 创建SequentialFeatureSelector对象，选择前3个最重要的特征\n",
    "sfs = SequentialFeatureSelector(model, n_features_to_select=2, direction='forward')\n",
    "\n",
    "# 拟合前向选择模型\n",
    "sfs.fit(iris.data, iris.target)\n",
    "\n",
    "# 打印所选特征的索引\n",
    "print(\"Selected Features: %s\" % (sfs.get_support(indices=True)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20605f38-d370-4c2c-8fb9-c598fb335b9a",
   "metadata": {},
   "source": [
    "后向消除（Backward Elimination）：\n",
    "\n",
    "后向消除从包含所有特征的全集开始，然后逐步剔除对模型性能贡献最小的特征，直到达到预定的特征数量或性能指标。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db6dd32a-9cbc-4f3d-ad39-dfdeeb6890f7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sklearn.feature_selection import SequentialFeatureSelector\n",
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "# 创建线性回归模型\n",
    "model = LinearRegression()\n",
    "\n",
    "# 创建SequentialFeatureSelector对象，选择前2个最重要的特征\n",
    "sfs = SequentialFeatureSelector(model, n_features_to_select=2, direction='backward')\n",
    "\n",
    "# 拟合后向消除模型\n",
    "sfs.fit(iris.data, iris.target)\n",
    "\n",
    "# 打印所选特征的索引\n",
    "print(\"Selected Features: %s\" % (sfs.get_support(indices=True)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "563f2aa4-463c-4476-a76f-44b250e7c192",
   "metadata": {},
   "source": [
    "正向逐步回归（Forward Stepwise Regression）：\n",
    "\n",
    "正向逐步回归是一种用于线性回归等模型的包装法方法，它逐步添加一个特征并调整模型，然后评估性能，然后继续添加下一个特征，直到达到特征数量或性能指标的预定条件。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b8b2806-4295-498f-a124-dff6dafb0257",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from mlxtend.feature_selection import SequentialFeatureSelector\n",
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "# 创建线性回归模型\n",
    "model = LinearRegression()\n",
    "\n",
    "# 创建SequentialFeatureSelector对象，选择前2个最重要的特征\n",
    "sfs = SequentialFeatureSelector(model, k_features=2, forward=True)\n",
    "\n",
    "# 拟合正向逐步回归模型\n",
    "sfs.fit(iris.data, iris.target)\n",
    "\n",
    "# 打印所选特征的索引\n",
    "print(\"Selected Features: {}\".format(sfs.k_feature_idx_))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99b704f7-9af8-470d-8d71-92596205be16",
   "metadata": {},
   "source": [
    "### 5.3 Embedded"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b13c882-b4eb-4aaa-a1a7-0e7231f9d42d",
   "metadata": {},
   "source": [
    "特征工程中的嵌入法是一种特征选择方法，它与模型训练过程相结合，通过训练模型来评估每个特征的重要性，并选择最重要的特征。嵌入法通常更高效，因为它考虑了特征之间的交互作用，但它要求你选择一个合适的机器学习算法，以便在训练过程中估计特征的重要性。以下是一些常见的嵌入法示例："
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b08ab4d-04cd-4f49-b9e3-b9e557e52acd",
   "metadata": {},
   "source": [
    "####  基于正则化的特征选择法\n",
    "\n",
    "正则化方法包括L1正则化（Lasso）和L2正则化（Ridge，在线性模型中被广泛使用，它们可以通过特征的系数来评估特征的重要性。较小的系数通常表示较不重要的特征。使用正则化的基模型，除筛选出特征外，同时也进行降维。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44b4a231-2657-4a7f-aae2-068d2e88ce0b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.datasets import load_iris\n",
    "\n",
    "# 加载示例数据集\n",
    "iris = load_iris()\n",
    "X, y = iris.data, iris.target\n",
    "\n",
    "# 创建L1正则化逻辑回归模型\n",
    "lr = LogisticRegression(penalty='l1', solver='liblinear')\n",
    "\n",
    "# 拟合模型\n",
    "lr.fit(X, y)\n",
    "\n",
    "# 打印特征的系数（重要性）\n",
    "coefficients = lr.coef_\n",
    "print(\"Feature Coefficients:\")\n",
    "for i, coef in enumerate(coefficients[0]):\n",
    "    print(f\"Feature {i+1}: {coef}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12a3e366-f0b1-4959-b253-8fea06ea820d",
   "metadata": {},
   "source": [
    "####  基于树模型的特征选择法"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "617f3a8b-1520-4039-aad1-a0b538f0d2d3",
   "metadata": {},
   "source": [
    "决策树、随机森林和梯度提升树等基于树的模型可以提供特征重要性分数。你可以使用这些分数来选择最重要的特征。树模型种GBDT也可用来作为基模型进行特征选择，使用feature_selection库的SelectFromModel类结合GBDT模型，来选择特征的代码如下："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d31bc21a-b196-4749-89f8-e1d6380adfc2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.datasets import load_iris\n",
    "\n",
    "# 加载示例数据集\n",
    "iris = load_iris()\n",
    "X, y = iris.data, iris.target\n",
    "\n",
    "# 创建随机森林分类器\n",
    "rf = RandomForestClassifier(n_estimators=100, random_state=0)\n",
    "\n",
    "# 拟合模型\n",
    "rf.fit(X, y)\n",
    "\n",
    "# 打印特征重要性分数\n",
    "feature_importances = rf.feature_importances_\n",
    "print(\"Feature Importances:\")\n",
    "for i, importance in enumerate(feature_importances):\n",
    "    print(f\"Feature {i+1}: {importance}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acb3eb94-7d2c-4402-b401-67cd2a82fb0e",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 第六节 特征转换：数据降维"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ec9ac08-f74e-498d-b52f-b7d1fb949e92",
   "metadata": {},
   "source": [
    "当特征选择完成后，可以直接训练模型了，但是可能由于特征矩阵过大，导致计算量大，训练时间长的问题，因此降低特征矩阵维度也是必不可少的。常见的降维方法除了以上提到的基于L1惩罚项的模型以外，另外还有主成分分析法（PCA）和线性判别分析（LDA），线性判别分析本身也是一个分类模型。PCA和LDA有很多的相似点，其本质是要将原始的样本映射到维度更低的样本空间中，但是PCA和LDA的映射目标不一样：PCA是为了让映射后的样本具有最大的发散性；而LDA是为了让映射后的样本有最好的分类性能。所以说PCA是一种无监督的降维方法，而LDA是一种有监督的降维方法。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5364fb88-39b5-45be-a20c-ff713b383481",
   "metadata": {},
   "source": [
    "### 6.1 主成分分析法（PCA）\n",
    "\n",
    "主成分分析（Principal Component Analysis，PCA）是一种常用的降维技术，它可以用于数据预处理、可视化、特征选择以及去除数据中的噪音。PCA的主要目标是通过线性变换将高维数据投影到低维空间中，同时最大程度地保留原始数据的方差。这些新的投影维度被称为主成分。\n",
    "\n",
    "以下是使用Python中的Scikit-Learn库进行主成分分析（PCA）的基本示例："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70d3a126-bdab-46cb-9637-30f0e025b8f1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "from sklearn.datasets import load_iris\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# 加载示例数据集\n",
    "iris = load_iris()\n",
    "X, y = iris.data, iris.target\n",
    "\n",
    "# 创建PCA模型，降维到2维\n",
    "pca = PCA(n_components=2)\n",
    "\n",
    "# 拟合模型\n",
    "X_pca = pca.fit_transform(X)\n",
    "\n",
    "# 打印投影后的数据\n",
    "print(\"PCA Projection:\")\n",
    "#print(X_pca)\n",
    "\n",
    "# 绘制投影后的数据\n",
    "plt.scatter(X_pca[:, 0], X_pca[:, 1], c=y)\n",
    "plt.xlabel('PCA Component 1')\n",
    "plt.ylabel('PCA Component 2')\n",
    "plt.title('PCA Projection of Iris Dataset')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d327f8c9-c52e-4713-9998-a12343bbff76",
   "metadata": {},
   "source": [
    "### 6.2 线性判别分析法（LDA）\n",
    "\n",
    "线性判别分析（Linear Discriminant Analysis，简称LDA）是一种用于降维和分类的统计技术，通常用于模式识别和机器学习任务。它是一种监督学习方法，主要用于解决分类问题，但也可以用于数据压缩和特征选择。LDA的主要目标是找到一个线性变换，将多维数据投影到低维空间中，以最大程度地减少类别之间的散布（使得不同类别的数据点尽可能分开），同时最大化同一类别数据点的散布。这样做的结果是在低维空间中更容易实现数据分类。以下是使用Python中的Scikit-Learn库进行线性判别分析（LDA）的基本示例："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82c2e7fb-6c54-4096-aee6-26eb0c2085e2",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sklearn.discriminant_analysis import LinearDiscriminantAnalysis\n",
    "from sklearn.datasets import load_iris\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# 加载示例数据集\n",
    "iris = load_iris()\n",
    "X, y = iris.data, iris.target\n",
    "\n",
    "# 创建LDA模型，降维到2维\n",
    "lda = LinearDiscriminantAnalysis(n_components=2)\n",
    "\n",
    "# 拟合模型\n",
    "X_lda = lda.fit_transform(X, y)\n",
    "\n",
    "# 打印投影后的数据\n",
    "print(\"LDA Projection:\")\n",
    "print(X_lda)\n",
    "\n",
    "# 绘制投影后的数据\n",
    "plt.scatter(X_lda[:, 0], X_lda[:, 1], c=y)\n",
    "plt.xlabel('LDA Component 1')\n",
    "plt.ylabel('LDA Component 2')\n",
    "plt.title('LDA Projection of Iris Dataset')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31f1a674-0c5d-41ab-9a8c-6c87ecba330f",
   "metadata": {},
   "source": [
    "## 参考文献"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bcb07785-50ed-4afd-aa50-9be85eefb51c",
   "metadata": {},
   "source": [
    "[1] https://feature-engine.readthedocs.io/en/1.3.x/index.html  \n",
    "[2] 锡南·厄兹代米尔, 迪夫娅·苏萨拉（庄嘉盛译），《特征工程入门与实践》，北京：人民邮电出版社，2019年。  \n",
    "[3] 爱丽丝·郑，阿曼达·卡萨丽（陈光欣译），《精通特征工程》，北京：人民邮电出版社，2019年。  \n",
    "[4] https://www.zhihu.com/question/29316149/answer/110159647  \n",
    "[5] https://zhuanlan.zhihu.com/p/26444240"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3081c86f-6939-4b24-b810-65a409cb4329",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
