fastai_NLP/basics_data_processing.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%reload_ext autoreload\n",
    "%autoreload 2\n",
    "# 多行输出\n",
    "from IPython.core.interactiveshell import InteractiveShell\n",
    "InteractiveShell.ast_node_interactivity = \"all\" "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 博客\n",
    "Great NLP-related blogs:\n",
    "- [Sebastian Ruder](http://ruder.io/)\n",
    "- [Jay Alammar](https://jalammar.github.io/)\n",
    "- [Abigail See](http://www.abigailsee.com/)\n",
    "- [Joyce Xu](https://medium.com/@joycex99)\n",
    "- [Stephen Merity](https://smerity.com/articles/articles.html)\n",
    "- [Rachael Tatman](https://towardsdatascience.com/evaluating-text-output-in-nlp-bleu-at-your-own-risk-e8609665a213)\n",
    "\n",
    "Other great technical blog posts:\n",
    "- [Peter Norvig](http://nbviewer.jupyter.org/url/norvig.com/ipython/ProbabilityParadox.ipynb) (more [here](http://norvig.com/ipython/))\n",
    "- [Julia Evans](https://codewords.recurse.com/issues/five/why-do-neural-networks-think-a-panda-is-a-vulture) (more [here](https://jvns.ca/blog/2014/08/12/what-happens-if-you-write-a-tcp-stack-in-python/))\n",
    "- [Julia Ferraioli](http://blog.juliaferraioli.com/2016/02/exploring-world-using-vision-twilio.html)\n",
    "- [Slav Ivanov](https://blog.slavv.com/picking-an-optimizer-for-style-transfer-86e7b8cba84b)\n",
    "- find [more on twitter](https://twitter.com/math_rachel)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 处理流程\n",
    "\n",
    "[自然语言处理时，通常的文本清理流程是什么？ - 知乎](https://www.zhihu.com/question/268849350)\n",
    "\n",
    "1. 获取语料库\n",
    "2. 去除非文本信息\n",
    "3. 拼写检查\n",
    "4. 转化为小写，中文的编码 ，去掉标点符号，\n",
    "5. Rare word replacement: 将词频小于5的替换成一个特殊的Token <UNK>\n",
    "6. Add <BOS>, <EOS>：在句前，句后加入<BOS>和<EOS>。这个是个依赖具体模型（Model-Based）的处理策略。因为现在都用RNN，RNN的结构对于结尾不敏感（Sensitive），所以要加。\n",
    "7. Long Sentence Cut-Off：将句长超过L的变成L。这个也是Model-Based，因为RNN需要Unrolling。\n",
    "5. 词性标注与Named Entity Recognition\n",
    "6. 词性还原与词干提取\n",
    "7. 引入停用词\n",
    "8. TF-IDF的特征处理\n",
    "9. 建立分析模型"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## NLP\n",
    "\n",
    "NLP 包含以下的任务:\n",
    "\n",
    "- 词性标注\n",
    "- 命名实体识别（人名、地名）\n",
    "- 问答系统 Question answering\n",
    "- 语音辨识 Speech recognition\n",
    "- 语音合成和语音识别 Text-to-speech and Speech-to-text\n",
    "- 主题模型 Topic modeling\n",
    "- 语义分类 Sentiment classification\n",
    "- 语言模型 Language modeling\n",
    "- 翻译 Translation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# !pip install -U nltk spacy gensim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# !sudo apt-get -y install protobuf-compiler libprotoc-dev"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# !pip install pytext-nlp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```bash\n",
    "$ git clone https://github.com/facebookresearch/fastText.git\n",
    "$ cd fastText\n",
    "$ make\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 基于矩阵分解的主题模型"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最佳分解是将文档聚类成两组，组与组之间具有尽可能彼此不同的单词分布，但在组中的文档中尽可能相似。 我们将这两个组称为“主题”。 我们会根据每个主题中最常出现的词汇将这些词汇分为两组。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sklearn.datasets import fetch_20newsgroups\n",
    "from sklearn import decomposition\n",
    "from scipy import linalg\n",
    "import matplotlib.pyplot as plt\n",
    "from pathlib import Path\n",
    "from dataclasses import dataclass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "newsgroups 数据集包括18,000个新闻组帖子，其中包含20个主题"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 获取数据集"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "root = Path('/home/lab/Datasets/scikit_learn_data')\n",
    "root_nltk = Path('/home/lab/Datasets/nltk_data')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "categories = ['alt.atheism', 'talk.religion.misc', 'comp.graphics', 'sci.space']\n",
    "remove = ('headers', 'footers', 'quotes')\n",
    "newsgroups_train = fetch_20newsgroups(data_home=root, subset='train', categories=categories, remove=remove)\n",
    "newsgroups_test = fetch_20newsgroups(data_home=root, subset='test', categories=categories, remove=remove)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2034,)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "(2034,)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "2034"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newsgroups_train.filenames.shape # 样本名\n",
    "newsgroups_train.target.shape\n",
    "len(newsgroups_train.data)\n",
    "# 2034 篇文章也即2034个样本"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['data', 'filenames', 'target_names', 'target', 'DESCR'])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newsgroups_train.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 样本真实内容"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "'/home/lab/Datasets/scikit_learn_data/20news_home/20news-bydate-train/comp.graphics/38816'"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hi,\n",
      "\n",
      "I've noticed that if you only save a model (with all your mapping planes\n",
      "positioned carefully) to a .3DS file that when you reload it after restarting\n",
      "3DS, they are given a default position and orientation.  But if you save\n",
      "to a .PRJ file their positions/orientation are preserved.  Does anyone\n",
      "know why this information is not stored in the .3DS file?  Nothing is\n",
      "explicitly said in the manual about saving texture rules in the .PRJ file. \n",
      "I'd like to be able to read the texture rule information, does anyone have \n",
      "the format for the .PRJ file?\n",
      "\n",
      "Is the .CEL file format available from somewhere?\n",
      "\n",
      "Rych\n"
     ]
    }
   ],
   "source": [
    "newsgroups_train.target[0]\n",
    "newsgroups_train.filenames[0]\n",
    "print(newsgroups_train.data[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 标签与标签名"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newsgroups_train.target[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'comp.graphics'"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newsgroups_train.target_names[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 话题数与话题单词数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_topics, num_top_words = 6, 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 停用词"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 停用词 stop words 对文本的含义没有太大帮助的词，比如 the, 通常会被过滤掉[停用词 - 维基百科，自由的百科全书](https://zh.wikipedia.org/wiki/%E5%81%9C%E7%94%A8%E8%AF%8D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.feature_extraction import stop_words"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['ENGLISH_STOP_WORDS',\n",
       " '__builtins__',\n",
       " '__cached__',\n",
       " '__doc__',\n",
       " '__file__',\n",
       " '__loader__',\n",
       " '__name__',\n",
       " '__package__',\n",
       " '__spec__']"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dir(stop_words)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['a',\n",
       " 'about',\n",
       " 'above',\n",
       " 'across',\n",
       " 'after',\n",
       " 'afterwards',\n",
       " 'again',\n",
       " 'against',\n",
       " 'all',\n",
       " 'almost',\n",
       " 'alone',\n",
       " 'along',\n",
       " 'already',\n",
       " 'also',\n",
       " 'although',\n",
       " 'always',\n",
       " 'am',\n",
       " 'among',\n",
       " 'amongst',\n",
       " 'amoungst']"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sorted(list(stop_words.ENGLISH_STOP_WORDS))[:20]  # 字典本身无序，转换成列表后排序"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- stop words 通常没有太多具体的含义\n",
    "- stop words 通常没有一个明确的停用词表能够适用于所有的工具"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 词干和词性还原 Stemming and Lemmatization\n",
    "- organize, organizes, and organizing 他们的词干是删除截尾的部分字母，最后的可能不是真正的单词，词性还原更加准确 cars -> car 名词还原\n",
    "- 词形还原（lemmatization），是把一个词汇还原为一般形式（能表达完整语义），方法较为复杂；而词干提取（stemming）是抽取词的词干或词根形式（不一定能够表达完整语义），方法较为简单\n",
    "- 不同语言不同"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "import nltk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[nltk_data] Downloading package wordnet to /home/lab/nltk_data...\n",
      "[nltk_data]   Package wordnet is already up-to-date!\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# nltk.download('wordnet', download_dir=root_nltk)\n",
    "nltk.download('wordnet')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 词干提取"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nltk import stem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 基于Porter词干提取算法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "wnl = stem.WordNetLemmatizer()\n",
    "porter = stem.porter.PorterStemmer()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "word_list = ['feet', 'foot', 'foots', 'footing']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['foot', 'foot', 'foot', 'footing']"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[wnl.lemmatize(word) for word in word_list]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['feet', 'foot', 'foot', 'foot']"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[porter.stem(word) for word in word_list]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "wl = ['organize', 'organizes', 'organizing']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['organize', 'organizes', 'organizing']"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[wnl.lemmatize(word) for word in wl]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['organ', 'organ', 'organ']"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[porter.stem(word) for word in wl]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- nltk 具有更多的方法，但是方法没有经过优化，Spacy 的方法较少，但是都是经过优化的"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "import spacy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "from spacy.lemmatizer import Lemmatizer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "# lemmatizer = Lemmatizer()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 下载英语模型"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "# !pip install https://github.com/explosion/spacy-models/releases/download/en_core_web_md-1.2.0/en_core_web_md-1.2.0.tar.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 词性还原"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [lemmatizer.lookup(word) for word in word_list]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- spacy 不提供词干提取只提供词性还原，这也是一种优化，因为 词性提取是更准确的"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 停用词 stop words 在各个库中也是不一样的"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "# nlp = spacy.load(\"en_core_web_sm\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sorted(list(nlp.Defaults.stop_words))[:20]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上述方法长期以来被认为是标准技术，但如果使用深度学习，它们通常会损害您的表现。 \n",
    "\n",
    "词干化，词形还原和删除停用词都涉及丢弃信息。\n",
    "\n",
    "但是，在处理更简单的模型时，它们仍然有用"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 词频"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.feature_extraction.text import CountVectorizer, TfidfVectorizer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[nltk_data] Downloading package punkt to /home/lab/nltk_data...\n",
      "[nltk_data]   Package punkt is already up-to-date!\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nltk.download('punkt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nltk import word_tokenize\n",
    "\n",
    "@dataclass\n",
    "class LemmaTokenizer:\n",
    "    wnl = stem.WordNetLemmatizer()\n",
    "    \n",
    "    def __call__(self, doc):\n",
    "        return [self.wnl.lemmatize(t) for t in word_tokenize(doc)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "vectorizer = CountVectorizer(stop_words='english')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 向量化 - 统计词频"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "vectors = vectorizer.fit_transform(newsgroups_train.data).todense() # (文档数, 词汇量)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2034, 26576)"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vectors.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "vocab = np.array(vectorizer.get_feature_names())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['00', '000', '0000', ..., 'zware', 'zwarte', 'zyxel'], dtype='<U80')"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vocab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(26576,)"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vocab.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['applaud', 'apple', 'applelink', ..., 'binoculars', 'binomial',\n",
       "       'bio'], dtype='<U80')"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vocab[3456:4567]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SVD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. 我们显然希望在一个主题中出现频率最高的单词在另一个主题中出现频率较低，否则这个单词就不能很好地分离出这两个主题。因此，我们**希望主题是正交**\n",
    "2. SVD算法将一个矩阵分解为一个正交列矩阵和一个正交行矩阵(以及一个对角矩阵，其中包含每个因子的相对重要性)\n",
    "3. SVD 是精确分解\n",
    "- SVD的应用\n",
    "    - 语义分析\n",
    "    - 协同过滤/推荐(获得Netflix奖)\n",
    "    - 数据压缩\n",
    "    - 主成分分析"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 4min 15s, sys: 34.9 s, total: 4min 50s\n",
      "Wall time: 25.4 s\n"
     ]
    }
   ],
   "source": [
    "%time U, s, Vh = linalg.svd(vectors, full_matrices=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[4.33926985e+02, 0.00000000e+00, 0.00000000e+00, ...,\n",
       "        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],\n",
       "       [0.00000000e+00, 2.91510127e+02, 0.00000000e+00, ...,\n",
       "        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],\n",
       "       [0.00000000e+00, 0.00000000e+00, 2.40711377e+02, ...,\n",
       "        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],\n",
       "       ...,\n",
       "       [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,\n",
       "        1.95236382e-15, 0.00000000e+00, 0.00000000e+00],\n",
       "       [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,\n",
       "        0.00000000e+00, 1.83950100e-15, 0.00000000e+00],\n",
       "       [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,\n",
       "        0.00000000e+00, 0.00000000e+00, 1.10014025e-15]])"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.diag(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2034, 2034)"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "(2034,)"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "(2034, 26576)"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U.shape\n",
    "s.shape\n",
    "Vh.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f076a271518>]"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAGXVJREFUeJzt3XtwXGd5x/Hvc/YmybIdO5Yd13awAwZiSiHUTeikUCDQuLTFAQYwvXnazGQ6EzowbadNYNrSP9zS0lKghbZpYXApJbjTSzwMFFIngTJTCA4kIY5jrFxIFCu24iS+yLKk3X36x3lXWsl7k629nOPfZ0Zzjt49u/v4rPw777777jnm7oiISHpF3S5ARETaS0EvIpJyCnoRkZRT0IuIpJyCXkQk5RT0IiIpp6AXEUk5Bb2ISMop6EVEUi7b7QIAVq1a5Rs3bux2GSIiiXLfffc96+5DzbbriaDfuHEj+/fv73YZIiKJYmY/amU7Dd2IiKScgl5EJOUU9CIiKaegFxFJOQW9iEjKKehFRFJOQS8iknKJDvrRExN87OuHeGzsdLdLERHpWYkO+qMnJ/nkXcM8cXy826WIiPSsRAe9haWuby4iUl+ygz4kvYJeRKS+ZAf9TJ9eRETqSXTQV6hDLyJSX6KDfnboRlEvIlJPooO+QjEvIlJfooNeH8aKiDSX7KCfnWDZ1TpERHpZsoNePXoRkabSEfTdLUNEpKclO+g1j15EpKmWg97MMmb2fTP7cvh9pZndaWaHw3JF1ba3mtmwmR0ys+vbUXg1Dd2IiNS3kB79+4GDVb/fAuxz983AvvA7ZrYF2AG8AtgGfNrMMotT7lyzQzdKehGReloKejNbD/wC8E9VzduB3WF9N3BDVfvt7j7p7o8Dw8DVi1PuvLrCUj16EZH6Wu3Rfxz4faBc1bbG3UcBwnJ1aF8HPFW13UhoW3T6MFZEpLmmQW9mvwgcc/f7WnzMWp+QnpPFZnaTme03s/1jY2MtPnTtp9IpEERE6mulR38t8DYzewK4HXiTmf0LcNTM1gKE5bGw/Qiwoer+64Ej8x/U3W9z963uvnVoaOi8ijdNuhERaapp0Lv7re6+3t03En/Iepe7/yqwF9gZNtsJ3BHW9wI7zKxgZpuAzcC9i145GqMXEWlF9gLu+xFgj5ndCDwJvAvA3Q+Y2R7gYaAI3OzupQuutAYLXXrNuhERqW9BQe/u9wD3hPXjwHV1ttsF7LrA2prSyI2ISHOJ/mZshYZuRETqS3TQ66RmIiLNJTvoK9Mru1yHiEgvS3bQ61KCIiJNJTroKxTzIiL1JTroTReYEhFpKuFBr3n0IiLNJDvou12AiEgCJDroK/RZrIhIfYkOep2mWESkuWQH/cxpirtciIhID0t20OtSgiIiTSU76MNSPXoRkfoSHfRojF5EpKlEB72hs5qJiDST7KDXRHoRkaYSHfQV6s+LiNSX6KDXh7EiIs0lO+gr57pR0ouI1JXsoA9LxbyISH3JDnpNuhERaSrZQa9LCYqINJXooEeXEhQRaSrRQa959CIizSU66EVEpLlEB73m0YuINJfsoNc1Y0VEmkp20IelevQiIvUlO+h1mmIRkaaSHfS6lKCISFPJDnpdSlBEpKlEB32FevQiIvUlOuj1hSkRkeYSHfQiItJcooN+9sNYjd2IiNST7KDXaYpFRJpKdtCHpXJeRKS+pkFvZn1mdq+ZPWBmB8zsT0L7SjO708wOh+WKqvvcambDZnbIzK5vV/GzlxJs1zOIiCRfKz36SeBN7v4q4NXANjN7LXALsM/dNwP7wu+Y2RZgB/AKYBvwaTPLtKP42R69kl5EpJ6mQe+x0+HXXPhxYDuwO7TvBm4I69uB29190t0fB4aBqxe16kBj9CIizbU0Rm9mGTO7HzgG3Onu3wHWuPsoQFiuDpuvA56quvtIaJv/mDeZ2X4z2z82NnZexZsm0ouINNVS0Lt7yd1fDawHrjazH2+wea30PafP7e63uftWd986NDTUWrX16ruge4uIpNuCZt24+wvAPcRj70fNbC1AWB4Lm40AG6ruth44csGVNi6srQ8vIpJkrcy6GTKzS8J6P/Bm4BFgL7AzbLYTuCOs7wV2mFnBzDYBm4F7F7vw2frUoxcRaSTbwjZrgd1h5kwE7HH3L5vZ/wF7zOxG4EngXQDufsDM9gAPA0XgZncvtaf8eJxIHXoRkfqaBr27PwhcVaP9OHBdnfvsAnZdcHUtMDNNrxQRaSDR34wF9ehFRJpJftBrjF5EpKHkB33N2ZwiIlKR+KAHDd2IiDSS/KA3netGRKSRxAe9gQbpRUQaSH7Q68NYEZGGkh/0mC4lKCLSQPKD3vRhrIhII8kPejR0IyLSSPKDXuekFxFpKPFBDxq6ERFpJPFBHw/dKOlFROpJfNCjD2NFRBpKfNBrhF5EpLHkB71pHr2ISCMpCHpNrxQRaST5QY/G6EVEGkl+0OtSgiIiDSU/6LtdgIhIj0t80IOGbkREGkl80JsZZQW9iEhdiQ/6bGSUyuVulyEi0rMSH/SZyCiqSy8iUlfigz6bMUoKehGRupIf9OrRi4g0lIKgjyiVFPQiIvUkPug1Ri8i0ljigz4eo9esGxGRehIf9OrRi4g0lvygN6OoMXoRkbqSH/SRUdI5EERE6kpF0Jc1dCMiUlcqgl49ehGR+hIf9JGpRy8i0kjig17fjBURaaxp0JvZBjO728wOmtkBM3t/aF9pZnea2eGwXFF1n1vNbNjMDpnZ9W39B0Q6142ISCOt9OiLwO+6+5XAa4GbzWwLcAuwz903A/vC74TbdgCvALYBnzazTDuKh3h6ZVlj9CIidTUNencfdffvhfVTwEFgHbAd2B022w3cENa3A7e7+6S7Pw4MA1cvduEVGZ29UkSkoQWN0ZvZRuAq4DvAGncfhfhgAKwOm60Dnqq620hoa4uMKehFRBppOejNbBD4d+AD7n6y0aY12s5JYjO7ycz2m9n+sbGxVss4h6ZXiog01lLQm1mOOOS/4O7/EZqPmtnacPta4FhoHwE2VN19PXBk/mO6+23uvtXdtw4NDZ1v/WF65XnfXUQk9VqZdWPAZ4CD7v6xqpv2AjvD+k7gjqr2HWZWMLNNwGbg3sUrea6sZt2IiDSUbWGba4FfA35gZveHtg8CHwH2mNmNwJPAuwDc/YCZ7QEeJp6xc7O7lxa98iDSPHoRkYaaBr27f4va4+4A19W5zy5g1wXU1bJMhKZXiog0kPhvxmrWjYhIY8kP+ijSuW5ERBpIQdCjMXoRkQYSH/SR5tGLiDSU+KDP6DTFIiINJT7os+rRi4g0lPigjyLDHfXqRUTqSHzQZyye4q9evYhIbYkP+igKQa8evYhITYkP+qyCXkSkocQHfSYEvebSi4jUlvigz2Xif4J69CIitSU+6Gd79DopvYhILYkPeo3Ri4g0lvign+nRlxT0IiK1JD7oK2P0+jBWRKS2xAd9ZmboRmP0IiK1JD7os5peKSLSUOKDXmP0IiKNJT7osxn16EVEGkl+0EeVL0xpjF5EpJYUBL2GbkREGkl80OtcNyIijSU+6DVGLyLSWPKDXmP0IiINJT7oNb1SRKSxxAe9hm5ERBpLfNAvyWcBODkx3eVKRER6U+KDfu3yPgBGT5ztciUiIr0p8UGfzUQM5DOMTxa7XYqISE9KfNADLClkGZ9S0IuI1JKKoB8sZDk9Wep2GSIiPSkVQb+koKEbEZF60hH0+SynFfQiIjWlIugvW97H48+Od7sMEZGelIqgX7+in+fHp7pdhohIT0pF0A/ksxTLzlRR57sREZmvadCb2WfN7JiZPVTVttLM7jSzw2G5ouq2W81s2MwOmdn17Sq82kA+A8AZTbEUETlHKz36zwHb5rXdAuxz983AvvA7ZrYF2AG8Itzn02aWWbRq66icBmF8SlMsRUTmaxr07v5N4Ll5zduB3WF9N3BDVfvt7j7p7o8Dw8DVi1RrXf2hRz+hHr2IyDnOd4x+jbuPAoTl6tC+DniqaruR0HYOM7vJzPab2f6xsbHzLCO2pBAH/bi+NCUico7F/jDWarTVPH+wu9/m7lvdfevQ0NAFPenAzNCNevQiIvOdb9AfNbO1AGF5LLSPABuqtlsPHDn/8loz82GsevQiIuc436DfC+wM6zuBO6rad5hZwcw2AZuBey+sxOYqPfoz0wp6EZH5ss02MLMvAm8AVpnZCPDHwEeAPWZ2I/Ak8C4Adz9gZnuAh4EicLO7tz19K2P0Z3QaBBGRczQNend/b52brquz/S5g14UUtVADOU2vFBGpJxXfjNX0ShGR+lIR9PlsRD4TcUpDNyIi50hF0AOsvaSPIy/ourEiIvOlJugvXznAk8d1qmIRkflSFfQ/eu5Mt8sQEek5qQn61Uv7eOHMNKVyzS/iiohctFIT9DpVsYhIbekJ+sqXpjSXXkRkjtQE/arBAgBP6NqxIiJzpCbor33JKvKZiLseOdZ8YxGRi0hqgn6wkOWKoSU8Ona626WIiPSU1AQ9wPoVAzypKZYiInOkKugvXznAU89N4K4pliIiFakK+ivXLmViusSDIye6XYqISM9IVdBf+5JVANx9SB/IiohUpCrof+ySfjavHuSBp17odikiIj0jVUEP8LrNQ/zv4WcZPTHR7VJERHpC6oL+N67diAOf3Dfc7VJERHpC6oJ+w8oBfvWay/nSd5/kwBF9KCsikrqgB/jAm1/KYCHLb3/x+0wWde4bEbm4pTLoVyzJ89fveTWPjY3zsTt/2O1yRES6KpVBD3DdlWt499b1/MM3HuOjX3tEX6ISkYtWttsFtNOut78Sw/jU3Y8SmfH+6zaTzaT22CYiUlOqgz6Xifizd7yS6VKZv7lrmEPPnOKT772Kvlym26WJiHRM6ru3UWT81btfxQff+nK+/vBRbvr8fRx5QXPsReTikfqgBzAzbnr9i/nTt7+Sbz92nDf+5T189GuPcOrsdLdLExFpu4si6Ct++ZrLuet3f5ZtP34Zn7r7Ud7w0Xv4/Ld/xHSp3O3SRETa5qIKeojPWf+JHVex933X8uLVg/zhfz3E9R//JnsfOML4pC4sLiLpY70w7XDr1q2+f//+jj+vu/M/B4/xZ189yGNj4xSyEW/Zsob3/NQGXnvFpeQ0Q0dEepiZ3efuW5ttl+pZN82YGW/ZsoY3vmyIe594jv9+6BnuuP8IX35wlP5chje9fDXXXLGSazZdyktWD5KJrNsli4gs2EXdo6/l7HSJux85xreGn+VrB47y7OlJAJbkM7x87TJeftlSrly7jCvXLuNlly1lsHBRHytFpIta7dEr6Btwd0aen+Dex5/jwZEXODh6ioPPnOTU2dmx/BddOsBL1yxlw4oB1q/oDz8DrFvRz/L+XBerF5G009DNIjAzNqwcYMPKAd75k+uBOPyffmGCR0ZPcXD0JAefOcnho6f51uFnmZieewK1pYUsly3v47LlfaxZ1seaZQXWLOtj9dI+Vi8rsHppgVWDBX2BS0TaSkG/QGbG+hUDrF8xwJu3rJlpd3eeG5/i6RcmGHl+gpHnz/D08xM8c/Isz5w4y+Gjpxk7PUmpfO47qIF8hpVL8ly6JM+lg4WZ9UsG8qwYyLG8P8fyyrI/x7L+HIP5LJE+MxCRFijoF4mZcelggUsHC/zE+ktqblMqxweDoyfPcuzUWY6dnOT4+BTPjU9x/HS8fvTkWQ6OnuT4+BRTxfrz+yODJYUsy/pyDBayDPZlWVLIsrSQZbAQry8pZOjPZxjIZRjIZ+P1fGjLZ+P13Gxbfy6DmQ4eImmjoO+gTGQMLS0wtLQALG+4rbszMV3i+TPTnDgzzYmJaU5MTHFyosjJs/Hvp84WOT1Z5PTZIqcm47annz/D+GSJ05NFxqeKLPQjmP5cHPp92Yi+XIZ8NiKXichnI/KZiFxYFrIRuYzF7dmIfCZDLmsUwrYz9wnrhXC/+bflazx2Plt5/EgznUQWQduC3sy2AZ8AMsA/uftH2vVcaWRmodedZd0l/ef1GO7OZLHMmakSZ6aKTEyVwnqJieni7PrMMm47WyxxdrrM2ekSU8Uy06UyU6UyU8UyZybitqliiemSx+vhtspyMWUiiw8CGSOfzdQ4wMweOAoNDkqVA0w2Y2QjIxvF65nIyEXxASW+LV7Phdsq22UjI5uJyJjNbBtZ3J6JKtvOrlfuG0XES0PvlqRr2hL0ZpYBPgW8BRgBvmtme9394XY8n9RmZvTlMvTl4s8AOsHdmS55fHCocRCYaS+WmSyVma5x22SxXHUQqRxs4oPWnMcoxeuTxTKnzhY5Pu+gVLmtsm23J5jNOSCYkQkHkJkDRsZmDyRRRFTz4DF3PTILB6yIjEEmiuL2GveNzMhEhOeIyETxSf8qzxlVljNtzLTNuX2mrep2C/ervt2MKGLu7fOeI5qpx2rXogPkomhXj/5qYNjdHwMws9uB7YCCPuXMjHw27nEvKXS7mrmKpTLFslMsO6WSUyzP/l65rVSOD1LxMv69WC5TDNuXylAKy2K5TNmdYtiu5GH7ksft4fFKleesun+x7JTLc7eZ2c4r9VU/TnmmtonpBvetWo9rK8/UVp6puduvxMJENnuQrD6gzD94zB6QmHOgnDnozduu+oCZzUTkKu/qZtbjd3O5aPbdYC6868uF7bKR0Z/PsCQffy42kM8wWMgyUMiwciDfM9e/aFfQrwOeqvp9BLimTc8l0pL4P263q+g+d6fs8eSAsnvVgcDnHBTObZs9mFVur36MObe7h4Na1e0+eyCaub3qOc6thXOetzRv29nHDAdgJxwE44Nq5UA52+ZMFkuUPBxwS7MHxunS7AF9uhQfJKdLznT5/N4NLu/PsXpp897OG142xId+Yct5vJKta1fQ13qvNWdXmdlNwE0Al19+eZvKEJH5zCwM82hIpFWVd1OVd3/Tpdl3ehPTJcYn48+3Tk8WOTMVT5B4YOQEZ6aanyhxzbK+ttffrqAfATZU/b4eOFK9gbvfBtwG8Tdj21SHiMgFi4d5FvZ28NfaVMv5aNcA0neBzWa2yczywA5gb5ueS0REGmhLj97di2b2PuBrxNMrP+vuB9rxXCIi0ljb5tG7+1eAr7Tr8UVEpDW9MfdHRETaRkEvIpJyCnoRkZRT0IuIpJyCXkQk5XriUoJmNgb86AIeYhXw7CKVs1hUU+t6sa5erAl6s65erAl6s67FrulF7j7UbKOeCPoLZWb7W7luYiepptb1Yl29WBP0Zl29WBP0Zl3dqklDNyIiKaegFxFJubQE/W3dLqAG1dS6XqyrF2uC3qyrF2uC3qyrKzWlYoxeRETqS0uPXkRE6kh00JvZNjM7ZGbDZnZLB593g5ndbWYHzeyAmb0/tH/YzJ42s/vDz1ur7nNrqPOQmV3fxtqeMLMfhOffH9pWmtmdZnY4LFd0qi4ze1nV/rjfzE6a2Qe6sa/M7LNmdszMHqpqW/C+MbOfDPt42Mw+aRdwUdM6NX3UzB4xswfN7D/N7JLQvtHMJqr22d93sKYFv16LWVODur5UVdMTZnZ/aO/UvqqXBV39uzqHuyfyh/j0x48CVwB54AFgS4eeey3wmrC+FPghsAX4MPB7NbbfEuorAJtC3Zk21fYEsGpe218At4T1W4A/73RdVa/ZM8CLurGvgNcDrwEeupB9A9wL/DTxldS+Cvz8Itf0c0A2rP95VU0bq7eb9zjtrmnBr9di1lSvrnm3/xXwRx3eV/WyoKt/V/N/ktyjn7kAubtPAZULkLedu4+6+/fC+ingIPF1cuvZDtzu7pPu/jgwTFx/p2wHdof13cANXarrOuBRd2/05bi21eTu3wSeq/F8Le8bM1sLLHP3//P4f+c/V91nUWpy96+7e+UadN8mvkJbXZ2oqYGO7KdmdYXe77uBLzZ6jDbsq3pZ0NW/q/mSHPS1LkDeKGzbwsw2AlcB3wlN7wtvuT9b9Xatk7U68HUzu8/i6/ICrHH3UYj/MIHVXagL4iuNVf9H7Pa+goXvm3VhvVP1/SZx765ik5l938y+YWavq6q1EzUt5PXq9H56HXDU3Q9XtXV0X83Lgp76u0py0De9AHnbCzAbBP4d+IC7nwT+Dngx8GpglPitJHS21mvd/TXAzwM3m9nrG2zbsbosvqTk24B/C029sK8aqVdHJ/fZh4Ai8IXQNApc7u5XAb8D/KuZLetQTQt9vTr9Or6XuZ2Iju6rGllQd9M6z9/W/ZXkoG96AfJ2MrMc8Qv7BXf/DwB3P+ruJXcvA//I7JBDx2p19yNheQz4z1DD0fDWsPLW9Vin6yI+8HzP3Y+G+rq+r4KF7psR5g6ltKU+M9sJ/CLwK+GtPOHt/vGwfh/x+O5LO1HTebxeHdlPAGaWBd4BfKmq3o7tq1pZQI/9XSU56Lt2AfIwHvgZ4KC7f6yqfW3VZm8HKrMD9gI7zKxgZpuAzcQfvCx2XUvMbGllnfhDvYfC8+8Mm+0E7uhkXcGcHle391WVBe2b8Db8lJm9Nvwd/HrVfRaFmW0D/gB4m7ufqWofMrNMWL8i1PRYh2pa0OvViZqqvBl4xN1nhj46ta/qZQG99ne1WJ/qduMHeCvxp9yPAh/q4PP+DPHbqgeB+8PPW4HPAz8I7XuBtVX3+VCo8xCL+Gn6vLquIP5E/wHgQGWfAJcC+4DDYbmyw3UNAMeB5VVtHd9XxAeaUWCauAd14/nsG2ArcdA9Cvwt4YuHi1jTMPE4buVv6+/Dtu8Mr+sDwPeAX+pgTQt+vRazpnp1hfbPAb81b9tO7at6WdDVv6v5P/pmrIhIyiV56EZERFqgoBcRSTkFvYhIyinoRURSTkEvIpJyCnoRkZRT0IuIpJyCXkQk5f4fh5RrFslyjHwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f076a0952e8>]"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xt81PWd7/HXJxcSIECAhFsSBCGoXORiRLqibVGqa61IWyteWrbdLbtnbb203d16+nhs29Njjw+3Wve0XoqWalcr0tZdWU9rRRSvXAyKKCASbiYQSbglQO7J5/wxv+iAIRly4TeZeT8fjzxm5ju/mXzCg7znl8/8fp8xd0dERBJXStgFiIhIz1LQi4gkOAW9iEiCU9CLiCQ4Bb2ISIJT0IuIJLiYg97MUs3sLTN7Jrj9IzPbY2Ybgq8rora93cxKzGyrmV3WE4WLiEhs0k5h21uALcDAqLWfu/vPojcys4nAAmASMAp43swmuHtzV4sVEZFTF9MevZnlA58HHo5h83nAUnevd/edQAkws/MliohIV8S6R38v8M/AgBPWv2VmXwOKge+6+yEgD1gTtU1ZsHZSOTk5PmbMmBhLERERgPXr1+9399yOtusw6M3sSqDC3deb2Wei7noA+AngweXdwDcAa+NpPjFnwcwWAYsARo8eTXFxcUeliIhIFDPbHct2sbRuLgSuMrNdwFJgjpk95u773L3Z3VuAh/i4PVMGFEQ9Ph/Ye+KTuvtidy9y96Lc3A5fkEREpJM6DHp3v93d8919DJE3WV9w9xvNbGTUZvOBd4Pry4EFZpZhZmOBQmBdN9ctIiIxOpWjbk50l5lNI9KW2QX8PYC7bzKzZcBmoAm4SUfciIiEx+JhTHFRUZGrRy8icmrMbL27F3W0nc6MFRFJcAp6EZEEp6AXEUlwvTro9xyu5a5n32Pv4dqwSxERiVu9OuiP1Tdx/6rtvLKtMuxSRETiVq8O+sJhWQwbkMGrJQfCLkVEJG716qA3M2aPz+G1kv20tIR/mKiISDzq1UEPcOH4HA4ea2DLh9VhlyIiEpcSIugBXivZH3IlIiLxqdcH/YhBmRQOy1KfXkTkJHp90ENkr37dzgPUNWqkjojIiRIi6GePz6GusYU3PzgUdikiInEnIYJ+1rihpKaY+vQiIm1IiKDPykhjekE2r25T0IuInCghgh4iffqNe6qoqmkMuxQRkbiSMEE/uzAHd1i9Q3v1IiLREibopxVk079PKq+qTy8icpyECfr01BRmnTlUfXoRkRMkTNBDpE+/60ANpQdrwi5FRCRuJFTQzy6MjEN4fbv26kVEWsUc9GaWamZvmdkzwe0hZrbCzLYFl4Ojtr3dzErMbKuZXdYThbeldWzxK2rfiIh85FT26G8BtkTd/j6w0t0LgZXBbcxsIrAAmARcDtxvZqndU277WscWv779gMYWi4gEYgp6M8sHPg88HLU8D3g0uP4ocHXU+lJ3r3f3nUAJMLN7yu2YxhaLiBwv1j36e4F/Blqi1oa7ezlAcDksWM8DSqO2KwvWTovWPr3GIYiIRHQY9GZ2JVDh7utjfE5rY+0TfRQzW2RmxWZWXFnZfZ/5OnxgZGyx+vQiIhGx7NFfCFxlZruApcAcM3sM2GdmIwGCy4pg+zKgIOrx+cDeE5/U3Re7e5G7F+Xm5nbhR2ij4PE5vLHroMYWi4gQQ9C7++3unu/uY4i8yfqCu98ILAcWBpstBJ4Ori8HFphZhpmNBQqBdd1eeTs0tlhE5GNdOY7+TmCumW0D5ga3cfdNwDJgM/AscJO7n9Zd69axxTpLVkQE0k5lY3dfBawKrh8ALjnJdncAd3Sxtk5rHVusN2RFRBLszNhoGlssIhKRsEF/UTC2WOMQRCTZJWzQT9XYYhERIIGDvnVssfr0IpLsEjboQWOLRUQgwYP+Io1DEBFJ7KAfH4wtVp9eRJJZQge9xhaLiCR40ENkmuXBYw1sLtfYYhFJTgkf9BeOV59eRJJbwgd969hi9elFJFklfNCDxhaLSHJLiqC/qDAYW7xbY4tFJPkkRdBfcGYwtljtGxFJQkkR9BpbLCLJLCmCHiKHWW7cU8XhmoawSxEROa2SJ+jHR8YWr95+IOxSREROq6QJeo0tFpFklTRBr7HFIpKskiboIdKn19hiEUk2HQa9mWWa2Toze9vMNpnZj4P1H5nZHjPbEHxdEfWY282sxMy2mtllPfkDnIrZGocgIkkoLYZt6oE57n7UzNKBV83sz8F9P3f3n0VvbGYTgQXAJGAU8LyZTXD30E9LHT8si+EDI2OLF8wcHXY5IiKnRYd79B5xNLiZHny1N/N3HrDU3evdfSdQAszscqXdwMy4UGOLRSTJxNSjN7NUM9sAVAAr3H1tcNe3zGyjmS0xs8HBWh5QGvXwsmDtxOdcZGbFZlZcWVnZhR/h1Mwer7HFIpJcYgp6d29292lAPjDTzCYDDwDjgGlAOXB3sLm19RRtPOdidy9y96Lc3NxOFd8ZGlssIsnmlI66cffDwCrgcnffF7wAtAAP8XF7pgwoiHpYPrC3G2rtFhpbLCLJJpajbnLNLDu43he4FHjPzEZGbTYfeDe4vhxYYGYZZjYWKATWdW/ZXTO7MId1OzW2WESSQyx79COBF81sI/AGkR79M8BdZvZOsP5Z4DYAd98ELAM2A88CN8XDETfRZo/Pob5JY4tFJDl0eHilu28Eprex/tV2HnMHcEfXSus5F5w5lLRgbPFfBT17EZFElVRnxrbKykhj+uhs9elFJCkkZdBD5OibdzS2WESSQNIGvcYWi0iySNqgn1qQTVZGGq+ofSMiCS5pgz4ytniITpwSkYSXtEEPkT79bo0tFpEEl9RBr7HFIpIMkjroW8cWq08vIoksqYP+o7HFJfs1tlhEElZSBz1E2jeHaho1tlhEElbSB33r2GKdJSsiiSrpg374wEwmDM/SG7IikrCSPughslevscUikqgU9GhssYgkNgU9H48t1mGWIpKIFPR8PLZYfXoRSUQK+oDGFotIolLQB1rHFr+uscUikmAU9IHWscU6nl5EEk2HQW9mmWa2zszeNrNNZvbjYH2Ima0ws23B5eCox9xuZiVmttXMLuvJH6C7aGyxiCSqWPbo64E57j4VmAZcbmazgO8DK929EFgZ3MbMJgILgEnA5cD9ZpbaE8V3N40tFpFE1GHQe8TR4GZ68OXAPODRYP1R4Org+jxgqbvXu/tOoASY2a1V95CLCjUOQUQST0w9ejNLNbMNQAWwwt3XAsPdvRwguBwWbJ4HlEY9vCxYi3vjciNjixX0IpJIYgp6d29292lAPjDTzCa3s7m19RSf2MhskZkVm1lxZWVlbNX2MI0tFpFEdEpH3bj7YWAVkd77PjMbCRBcVgSblQEFUQ/LB/a28VyL3b3I3Ytyc3M7UXrP0NhiEUk0sRx1k2tm2cH1vsClwHvAcmBhsNlC4Ong+nJggZllmNlYoBBY192F95TZGlssIgkmlj36kcCLZrYReINIj/4Z4E5grpltA+YGt3H3TcAyYDPwLHCTu/easZDDNLZYRBJMWkcbuPtGYHob6weAS07ymDuAO7pcXUguHJ/D79Z+wN7DtYzK7ht2OSIiXaIzY9uw8FNjSEsxvv3EWzQ2t4RdjohIlyjo2zAmpz//50vnsn73Ie5+7v2wyxER6RIF/UlcNXUU180czYMvbefF9yo6foCISJxS0Lfjh1+YyNkjBvCdZRsor6oNuxwRkU5R0LcjMz2V+26YQX1TCzc/8RZN6teLSC+koO/AuNwsfjp/Cm/sOsQ9K9SvF5HeR0Efg6un57Hg/ALuX7Wdl96Pj3ENIiKxUtDH6IdfmMRZwwdw25Mb+LCqLuxyRERipqCPUd8+kX59XWMzNy9Vv15Eeg8F/SkYPyyL/331ZNbtPMi9z28LuxwRkZgo6E/RF2fk85WifO5bVcIr29SvF5H4p6DvhB9fNZnCYVncunQD+6rVrxeR+Kag74S+fVK57/oZ1DQ0c8vSt2jWh5SISBxT0HdS4fAB/OTqyazZcZB/X6l+vYjELwV9F3z5vHy+NCOfX7ywjVe3aX69iMQnBX0X/eTqSYzLzeLWJzdQcUT9ehGJPwr6LurXJ437b5jB0fpGbnlig/r1IhJ3FPTdYMLwAfyvqyazescBfvGC+vUiEl8U9N3kmqJ8vjg9j39fuY3X9XmzIhJHFPTdxMz4ydWTOTOnP7c8uYHKI/VhlyQiAsQQ9GZWYGYvmtkWM9tkZrcE6z8ysz1mtiH4uiLqMbebWYmZbTWzy3ryB4gn/TPSuO+GGVTXNnLbk+rXi0h8iGWPvgn4rrufA8wCbjKzicF9P3f3acHXnwCC+xYAk4DLgfvNLLUHao9LZ48YyI+vmsSrJfu578WSsMsREek46N293N3fDK4fAbYAee08ZB6w1N3r3X0nUALM7I5ie4trzy9g3rRR3Pv8+6zefiDsckQkyZ1Sj97MxgDTgbXB0rfMbKOZLTGzwcFaHlAa9bAy2n9hSDhmxh3zpzBmaH9uWfoW+4+qXy8i4Yk56M0sC/gjcKu7VwMPAOOAaUA5cHfrpm08/BPNajNbZGbFZlZcWZl4UyCzMtL45fUzOBz061vUrxeRkMQU9GaWTiTkH3f3pwDcfZ+7N7t7C/AQH7dnyoCCqIfnA3tPfE53X+zuRe5elJub25WfIW5NHDWQH31hEq9s28/9q9SvF5FwxHLUjQG/Bra4+z1R6yOjNpsPvBtcXw4sMLMMMxsLFALruq/k3uW6mQV8Yeoo7lnxPmt3qF8vIqdfLHv0FwJfBeaccCjlXWb2jpltBD4L3Abg7puAZcBm4FngJndv7pny45+Z8dP5kzljaH9uXvoWB9SvF5HTzNzD7x0XFRV5cXFx2GX0qE17q5h//+uMz83ix/Mmcf6YIWGXJCK9nJmtd/eijrbTmbGnyaRRg/jFddM5cKyeax5czTd/W0xJxZGwyxKRJKCgP40umzSCVd/7LP902Vms3n6Az/38ZW5/6h0q9HGEItKD1LoJyYGj9fzihRIeW7Ob9NQUvnnRWBZ9ehxZGWlhlyYivUSsrRsFfch27T/Gvz23lf+3sZyh/ftwy6WFXDdzNOmp+mNLRNqnHn0vMSanP/ddP4P/uulCxg/L4l+f3sTnfv4yf3qnnHh4ERaR3k9BHyemFWSzdNEslvxNEempxj8+/iZffOB11u08GHZpItLLKejjiJkx5+zh/PmWi7nrS+ey93AtX/nVav7uUR2hIyKdpx59HKttaGbJazt5YNV2ahqauPb80dx2aSHDBmaGXZqIxAG9GZtAWo/QeXztbtJSdISOiEQo6BPQ7gPH+Le/bOUZHaEjIuiom4R0xtD+/PL6GTwddYTO3Hte0hE6ItIuBX0vNDXqCJ0+aSn84+NvcuuTG2hqbgm7NBGJQ2ry9lKtR+h8esIw7nuxhHtWvE9Ts3Pvgmlq5YjIcRT0vVxqinHzJYX0TU/ljj9toamlhV9cN4M+aQp7EYlQGiSIb158Jv965UT+smkf//j4euqbkvYjAETkBAr6BPKN2WP5ybxJPL+lgn/4j/XUNSrsRURBn3C++qkx/HT+FF7cWsk3f1ussBcRBX0iuv6C0dz15XN5tWQ/33jkDWoamsIuSURCpKBPUF8pKuCer0xlzY4D/M1v3uBYvcJeJFl1GPRmVmBmL5rZFjPbZGa3BOtDzGyFmW0LLgdHPeZ2Mysxs61mdllP/gBycvOn53Pvgums332IhUvWcaSuMeySRCQEsezRNwHfdfdzgFnATWY2Efg+sNLdC4GVwW2C+xYAk4DLgfvNLLUnipeOXTV1FL+4bjobSg/ztSXrqFbYiySdDoPe3cvd/c3g+hFgC5AHzAMeDTZ7FLg6uD4PWOru9e6+EygBZnZ34RK7K6aM5L4bZvDunipufHgtVTUKe5Fkcko9ejMbA0wH1gLD3b0cIi8GwLBgszygNOphZcGahOiySSN48MbzeK/8CNc/vIZDxxrCLklETpOYg97MsoA/Are6e3V7m7ax9omJW2a2yMyKzay4srIy1jKkCy45ZziLv3Ye2yqOct1DazhwtD7skkTkNIgp6M0snUjIP+7uTwXL+8xsZHD/SKAiWC8DCqIeng/sPfE53X2xuxe5e1Fubm5n65dT9JmzhvHrhUXs3H+M6x5aQ+URhb1IoovlqBsDfg1scfd7ou5aDiwMri8Eno5aX2BmGWY2FigE1nVfydJVFxXm8puvn0/pwVoWLF5NRXVd2CWJSA+KZY/+QuCrwBwz2xB8XQHcCcw1s23A3OA27r4JWAZsBp4FbnJ3nZ4ZZ/5qXA6PfP18yqvquHbxGsqrasMuSUR6iD5hKsmt332QhUveYEj/PjyxaBZ52X3DLklEYqRPmJKYnHfGEP7jb2dyqKaBa3+1mtKDNWGXJCLdTEEvTB89mMf/7gKO1DVx7a9Ws/vAsbBLEpFupKAXAM7Nz+Z337yA2sZmrv3VGnZUHg27JBHpJgp6+cikUYN4YtEsGptbWLB4DSUVCnuRRKCgl+OcPWIgSxfNosVhweLV/Oa1nTqLVqSXU9DLJxQOH8DS4AicH//3Zi746Upu+t2bvPR+Jc0t4R+lJSKnRodXSru2lFezrLiU/3prD4dqGhk1KJMvn5fPl88rYPTQfmGXJ5LUYj28UkEvMalvambllgqefKOUl7dV4g6fOnMo155fwOWTR5CZrknUIqebgl56zN7DtTz1ZhnLisv44GANAzLTuGrqKK49v4ApeYOITM0QkZ6moJce19LirN15kN8Xl/Knd8upa2zh7BEDuKaogPnT8xjSv0/YJYokNAW9nFbVdY3899t7WfZGKW+XVZGeasydOJxrigq4uDCX1BTt5Yt0NwW9hOa9D6v5fXEZ//nWHg4ea2DEwMgbuNcU5XPG0P5hlyeSMBT0ErqGphZWbtnHsuJSXnq/khaHWWcO4bqZo7ny3FHayxfpIgW9xJXyqlqeenMPy4pL2X2ghvHDsvje5yZw2aQRevNWpJMU9BKX3J1n3/2Qnz23le2Vx5iaP4h/uuxsZhfmhF2aSK+jMcUSl8yMv54ykr/cejF3fflc9h9t4MZfr+X6h9bw1geHwi5PJCFpj15CVd/UzONrPuC+F0s4cKyBuROH873PncVZIwaEXZpI3FPrRnqVY/VNLHl1J4tf3sHRhibmT8vj1ksnaMyCSDsU9NIrHTrWwIMvb+eR13bR4s6C80fz7TnjGTYwM+zSROKOgl56tX3Vdfzfldt48o1S0lKNr184ln+4eByD+qWHXZpI3Oi2N2PNbImZVZjZu1FrPzKzPWa2Ifi6Iuq+282sxMy2mtllnf8RJJkNH5jJHfOn8Px3Ps3lk0bw4EvbmX3XC9z3Ygk1DU1hlyfSq3S4R29mFwNHgd+6++Rg7UfAUXf/2QnbTgSeAGYCo4DngQnu3tze99AevXRkS3k1dz+3lee3VJCTlcG354xnwcwCMtI0NVOSV7ft0bv7y8DBGL/vPGCpu9e7+06ghEjoi3TJOSMH8vDC8/nj//gU43L788Plm7jk7pf4w/oyfRiKSAe6chz9t8xsY9DaGRys5QGlUduUBWsi3eK8M4awdNEsfvuNmQzu14fv/f5tLr/3ZZ5990Pi4f0mkXjU2aB/ABgHTAPKgbuD9bbOZW/zt8/MFplZsZkVV1ZWdrIMSUZmxsUTcln+rQu5/4YZNLvzD4+tZ959r/Hi1goFvsgJOhX07r7P3ZvdvQV4iI/bM2VAQdSm+cDekzzHYncvcvei3NzczpQhSc7MuGLKSJ4LzrI9eKyBr//mDb70wOu8um2/Al8k0KmgN7ORUTfnA61H5CwHFphZhpmNBQqBdV0rUaR9aakpfKWogBe++xnumD+Z8qo6bvz1Wq5dvIa1Ow6EXZ5I6NI62sDMngA+A+SYWRnwQ+AzZjaNSFtmF/D3AO6+ycyWAZuBJuCmjo64EekufdJSuOGCM/jyefksXVfKfS+WcO3iNVw4fijfmXsW550xuOMnEUlAOmFKElZdYzOPrdnNgy9tZ//RBj49IZfvzJ3A1ILssEsT6RY6M1YkUNPQxG9X7+ZXL23nUE0jl54zjNvmTmDSqEFhlybSJQp6kRMcrW/ikdcig9Oq65q4fNIIbps7QZMypddS0IucRFVtI0te3cmSV3dytKGJz08Zya2XTmD8sKywSxM5JQp6kQ4crmngoVd28JvXdlHX2My8aXncckkhY3L0AebSOyjoRWJ04Gg9i1/ewaOrd9HY7Hxxeh43X1JIwRDNwpf4pqAXOUUVR+p4cNUOHlu7m5YW55qiAr49ZzyjsvuGXZpImxT0Ip30YVUd968q4Yl1H2AYo4f2Y2BmGgMy0xnYN50BmWkMzGy9TDthLZ2BfSPb9u+TillbU0FEuoeCXqSL9hyu5dHXd1F2qIbq2iaO1DVSXRdc1jbR0NzS7uNTU4ysjLRI8Gd8/AIwMDOdwf3SmThqINMKshmb018vCNIpsQZ9h2fGiiSrvOy+/M8rzjnp/XWNzVTXNXKkronq2uCyg9ulB2s4UtfE/qP11DdFXigG9U1nakE20wqymV6QzdSCbIb073O6fkxJAgp6kU7KTE8lMz2VYZ04DL+5xdlWcYQNHxxmQ2nk65cvbKN1tP7oIf2YFoT/tNHZTBw5kMx0fciKdI5aNyJx4lh9ExvLqni77PBHLwAfVtcBkJ5qTBwZafW07v2r5SPq0YskgA+r6thQeoi3SiPh/86eKmoaInMC1fIRBb1IAmqr5fP+viMftXzysvtSMKQvBYP7UTCk33HXc7MySEnRXwCJRG/GiiSg1BTj7BEDOXvEQBbMHA1EWj7v7KliQ+lh3iuvpvRQLS+9X0nFkfrjHtsnLYX8wa3BH/ViENwe1DddraAEpaAX6eX6Z6Qx68yhzDpz6HHrdY3NlB2qpfRQDWUHayg9VEvpwRpKD9WwofQwVbWNx20/ICON/CH9KBjcN3gBiFyOHtKP7H59SE810lJTSE810lNS9NdBL6KgF0lQmempjB+WddJhbdV1jZHgP1hL2aGa4EWglp37j/HytkrqGts/TyDFID01hfTUFNJSjbSUFPoELwZpqUaf49aD66kppKcY6akp9ElLYWxOf6YWDGJKXja5AzJ64p9BUNCLJK2BmelMGjWozbn87s7+ow2UBi8A1bWNNDY7TS0tNDY7jc0tNAWXH68H15ujtmlp3aaFhqYWjjU009jUQlNLC3WNLfz3xr20vk04alAmU/IHcW5+NufmD2JK3iCy++nN5e6goBeRTzAzcgdkkDsggxmje+4jGI/VN7FpbzUbyw6zsayKjWWH+cumfR/df8bQfkzJG8TU/Gym5A9ict4gsjIUW6dK/2IiEpr+GWnMHDuEmWOHfLRWVdPIO3uq2LjnMBtLq3jrg8M8s7EcADMYl5vFufmDODdvEOcW6GSyWCjoRSSuDOqXzuzCHGYX5ny0Vnmknnf3RE4me6esipff389Tb+4BIkciTRg+gKn5g5gStHxGD+mno4iidHgcvZktAa4EKtx9crA2BHgSGAPsAr7i7oeC+24H/hZoBm529790VISOoxeRU+HufFhdx9ulVbyzp7XtU3XckURZGWnkZfclb3Bf8gf3jbrej7zsvuRk9en1LwTddsKUmV0MHAV+GxX0dwEH3f1OM/s+MNjd/8XMJgJPADOBUcDzwAR3b27veyjoRaSr3J0PDtawpbyaskO1H33tORw5quhIXdNx22ekpZAXvADkD+730YtB/uDIC8KwAZmkxvkhpN12wpS7v2xmY05Yngd8Jrj+KLAK+Jdgfam71wM7zayESOivjrVwEZHOMDPOGNqfM4a2/VGQ1XWN7GkN/0M1wQtA5IVg894POXCs4bjt01ONkYM+/ktgwvAs5pw9vFd+tnBne/TD3b0cwN3LzWxYsJ4HrInarixY+wQzWwQsAhg9enQnyxARic3AzHQGjkznnJED27y/pqGJvYeP/0sg8sJQwyvbKvnD+jJ++qf3ODOnP3MnDmfuxOFMHz047vf6ofvfjG3rJ26zN+Tui4HFEGnddHMdIiKnpF+fNMYPG8D4k8yd3nu4lue37GPF5n0seW0nv3p5B0P792HO2cOYO3E4FxXm0rdPfB7909mg32dmI4O9+ZFARbBeBhREbZcP7O1KgSIi8WBUdl++9qkxfO1TY6iua+SlrZWs2LyPZzd9yO/Xl5GZnsLs8bl8buJw5pwzjJys+DnTt7NBvxxYCNwZXD4dtf47M7uHyJuxhcC6rhYpIhJPBmam84Wpo/jC1FE0NLWwbudBVmz+kBWb9/H8ln2YwYzRgz9q8YzLDbevH8tRN08QeeM1B9gH/BD4L2AZMBr4ALjG3Q8G2/8A+AbQBNzq7n/uqAgddSMiicDd2bS3+qMWz6a91QCcmRvp639u4nCmFXRfX1/z6EVEQrbncC3Pb46E/podB2hqcXKyWvv6I7ioMKdLZ/Uq6EVE4khVbSOrtlbw/JYKVr1XwZH6JjLTU/jqrDP4wecnduo59cEjIiJxZFDfdOZNy2PetDwamlpYu/MAKzbvY+Sgvj3+vRX0IiKnWZ+0FC4qzOWiwtzT8v1STst3ERGR0CjoRUQSnIJeRCTBKehFRBKcgl5EJMEp6EVEEpyCXkQkwSnoRUQSXFyMQDCzSmB3F54iB9jfTeX0tN5UK/SuelVrz+lN9famWqFr9Z7h7h2edRUXQd9VZlYcy7yHeNCbaoXeVa9q7Tm9qd7eVCucnnrVuhERSXAKehGRBJcoQb847AJOQW+qFXpXvaq15/SmentTrXAa6k2IHr2IiJxcouzRi4jISfTqoDezy81sq5mVmNn3w66nPWZWYGYvmtkWM9tkZreEXVNHzCzVzN4ys2fCrqUjZpZtZn8ws/eCf+NPhV3TyZjZbcH/gXfN7Akzywy7pmhmtsTMKszs3ai1IWa2wsy2BZeDw6yx1Ulq/bfg/8FGM/tPM8sOs8ZobdUbdd/3zMzNLKe7v2+vDXozSwXuA/4amAhcZ2ad+zyu06MJ+K67nwPMAm6K83oBbgG2hF1EjP4deNbdzwamEqd1m1koJzsHAAADC0lEQVQecDNQ5O6TgVRgQbhVfcIjwOUnrH0fWOnuhcDK4HY8eIRP1roCmOzu5wLvA7ef7qLa8QifrBczKwDmAh/0xDfttUEPzARK3H2HuzcAS4F5Idd0Uu5e7u5vBtePEAmivHCrOjkzywc+Dzwcdi0dMbOBwMXArwHcvcHdD4dbVbvSgL5mlgb0A/aGXM9x3P1l4OAJy/OAR4PrjwJXn9aiTqKtWt39OXdvCm6uAfJPe2EncZJ/W4CfA/8M9Mibpr056POA0qjbZcRxcEYzszHAdGBtuJW0614i//Fawi4kBmcClcBvglbTw2bWP+yi2uLue4CfEdlzKweq3P25cKuKyXB3L4fITgswLOR6YvUN4M9hF9EeM7sK2OPub/fU9+jNQW9trMX9IURmlgX8EbjV3avDrqctZnYlUOHu68OuJUZpwAzgAXefDhwjfloLxwl62/OAscAooL+Z3RhuVYnJzH5ApGX6eNi1nIyZ9QN+APxrT36f3hz0ZUBB1O184uxP4BOZWTqRkH/c3Z8Ku552XAhcZWa7iLTE5pjZY+GW1K4yoMzdW/9C+gOR4I9HlwI73b3S3RuBp4C/CrmmWOwzs5EAwWVFyPW0y8wWAlcCN3h8H0M+jsiL/tvB71s+8KaZjejOb9Kbg/4NoNDMxppZHyJvaC0PuaaTMjMj0kPe4u73hF1Pe9z9dnfPd/cxRP5dX3D3uN3rdPcPgVIzOytYugTYHGJJ7fkAmGVm/YL/E5cQp28cn2A5sDC4vhB4OsRa2mVmlwP/Alzl7jVh19Med3/H3Ye5+5jg960MmBH8n+42vTbogzdbvgX8hcgvyjJ33xRuVe26EPgqkb3jDcHXFWEXlUC+DTxuZhuBacBPQ66nTcFfHX8A3gTeIfI7GFdncprZE8Bq4CwzKzOzvwXuBOaa2TYiR4fcGWaNrU5S6y+BAcCK4PfswVCLjHKSenv++8b3XzUiItJVvXaPXkREYqOgFxFJcAp6EZEEp6AXEUlwCnoRkQSnoBcRSXAKehGRBKegFxFJcP8f9jLTNsGbNOoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(s[:15])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([433.92698542, 291.51012741, 240.71137677, 220.00048043,\n",
       "       182.74475229, 168.15235371, 148.15722056, 128.78642   ,\n",
       "       126.49055396, 118.34718466])"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_topics, num_top_words = 6, 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 主题 topics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [],
   "source": [
    "def top_words(t):\n",
    "    # 按 Vh 的数值对单词进行排序\n",
    "    vs = [vocab[i] for i in np.argsort(t)]\n",
    "    # 取出关键词，排名前几的，数值最大的\n",
    "    # 按照从大到小的顺序\n",
    "    kw = vs[:-num_top_words-1:-1]\n",
    "    return kw\n",
    "    \n",
    "def show_topics(a):\n",
    "    # 对每个文档计算排名靠前的词\n",
    "    topic_words = [top_words(t) for t in a]\n",
    "    return [' '.join(t) for t in topic_words]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- topic文章的关键词"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['ditto critus propagandist surname galacticentric kindergarten surreal imaginative salvadorans autotheism',\n",
       " 'jpeg gif file color quality image jfif format bit version',\n",
       " 'graphics edu pub mail 128 3d ray ftp send image',\n",
       " 'jesus god matthew people atheists atheism does graphics religious said',\n",
       " 'image data processing analysis software available tools display tool user',\n",
       " 'god atheists atheism religious believe religion argument true atheist example']"
      ]
     },
     "execution_count": 113,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "show_topics(Vh[:num_topics])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "与其限制因子正交，另一个想法是限制它们是非负的，**更容易解释，非精确分解**\n",
    "\n",
    "$$V = WH$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_topics = 6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf = decomposition.NMF(n_components=num_topics, random_state=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2034, 26576)"
      ]
     },
     "execution_count": 103,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "matrix([[0, 0, 0, ..., 0, 0, 0],\n",
       "        [0, 0, 0, ..., 0, 0, 0],\n",
       "        [0, 0, 0, ..., 0, 0, 0],\n",
       "        ...,\n",
       "        [0, 0, 0, ..., 0, 0, 0],\n",
       "        [0, 0, 0, ..., 0, 0, 0],\n",
       "        [0, 0, 0, ..., 0, 0, 0]], dtype=int64)"
      ]
     },
     "execution_count": 103,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vectors.shape\n",
    "vectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [],
   "source": [
    "W1 = clf.fit_transform(vectors)\n",
    "H1 = clf.components_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2034, 6)"
      ]
     },
     "execution_count": 106,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "array([[8.84627696e-02, 2.97289582e-02, 0.00000000e+00, 2.37507804e-02,\n",
       "        0.00000000e+00, 3.72347995e-02],\n",
       "       [0.00000000e+00, 7.58983467e-04, 3.52199796e-03, 1.36001870e-02,\n",
       "        9.86395268e-04, 1.61400988e-02],\n",
       "       [0.00000000e+00, 1.64550346e-02, 2.04416652e-04, 9.72735106e-03,\n",
       "        2.82118112e-03, 2.79858305e-02],\n",
       "       ...,\n",
       "       [9.07954483e-03, 3.00117203e-02, 3.28471200e-03, 1.43810275e-02,\n",
       "        3.39512238e-03, 4.71436750e-03],\n",
       "       [1.66879110e-02, 0.00000000e+00, 5.28028127e-05, 3.60103660e-02,\n",
       "        0.00000000e+00, 2.95061909e-01],\n",
       "       [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "        0.00000000e+00, 0.00000000e+00]])"
      ]
     },
     "execution_count": 106,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "W1.shape\n",
    "W1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(6, 26576)"
      ]
     },
     "execution_count": 107,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "array([[1.21314189e-01, 0.00000000e+00, 0.00000000e+00, ...,\n",
       "        3.18193652e-05, 6.36387303e-05, 2.94728045e-04],\n",
       "       [1.19659887e-01, 1.21895163e-01, 1.66521367e-04, ...,\n",
       "        1.16507442e-04, 2.33014883e-04, 5.20256298e-02],\n",
       "       [6.26357335e-02, 5.25135886e-01, 8.69181875e-04, ...,\n",
       "        3.11413247e-04, 6.22826494e-04, 0.00000000e+00],\n",
       "       [2.38710449e-03, 2.23453097e-01, 0.00000000e+00, ...,\n",
       "        1.55491031e-04, 3.10982062e-04, 0.00000000e+00],\n",
       "       [1.27275233e-01, 1.93687336e-01, 5.04255146e-04, ...,\n",
       "        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],\n",
       "       [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,\n",
       "        0.00000000e+00, 0.00000000e+00, 0.00000000e+00]])"
      ]
     },
     "execution_count": 107,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H1.shape\n",
    "H1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 两种方法的主题关键词对比\n",
    "- **SVD和NMF都可以获得主题与主题对应的关键词**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['jpeg image gif file color images format quality version files',\n",
       " 'edu graphics pub mail 128 ray ftp send 3d com',\n",
       " 'space launch satellite nasa commercial satellites year market data earth',\n",
       " 'jesus matthew prophecy people said messiah david isaiah psalm lord',\n",
       " 'image data available software processing ftp edu analysis images display',\n",
       " 'god atheists atheism religious believe people religion does argument true']"
      ]
     },
     "execution_count": 114,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "['ditto critus propagandist surname galacticentric kindergarten surreal imaginative salvadorans autotheism',\n",
       " 'jpeg gif file color quality image jfif format bit version',\n",
       " 'graphics edu pub mail 128 3d ray ftp send image',\n",
       " 'jesus god matthew people atheists atheism does graphics religious said',\n",
       " 'image data processing analysis software available tools display tool user',\n",
       " 'god atheists atheism religious believe religion argument true atheist example']"
      ]
     },
     "execution_count": 114,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "show_topics(H1)\n",
    "show_topics(Vh[:num_topics])\n",
    "# 有一个不重合"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 基于 TF-IDF 的主题模型"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [],
   "source": [
    "tf = TfidfVectorizer(stop_words='english')\n",
    "tfidf = tf.fit_transform(newsgroups_train.data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 148,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0., 0., 0., ..., 0., 0., 0.],\n",
       "       [0., 0., 0., ..., 0., 0., 0.],\n",
       "       [0., 0., 0., ..., 0., 0., 0.],\n",
       "       ...,\n",
       "       [0., 0., 0., ..., 0., 0., 0.],\n",
       "       [0., 0., 0., ..., 0., 0., 0.],\n",
       "       [0., 0., 0., ..., 0., 0., 0.]])"
      ]
     },
     "execution_count": 148,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tfidf.toarray()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_topics = 6\n",
    "clf = decomposition.NMF(n_components=num_topics, random_state=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "metadata": {},
   "outputs": [],
   "source": [
    "W1 = clf.fit_transform(tfidf)\n",
    "H1 = clf.components_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['don people just think like know say religion ve time',\n",
       " 'thanks graphics files image file program windows format know ftp',\n",
       " 'space nasa launch shuttle orbit lunar moon earth station data',\n",
       " 'ico bobbe tek beauchaine bronx manhattan sank queens vice blew',\n",
       " 'god jesus bible believe atheism christian does belief exist faith',\n",
       " 'objective morality values moral subjective science absolute claim freedom goal']"
      ]
     },
     "execution_count": 151,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "show_topics(H1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 156,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.00000000e+00, 2.84560108e-02, 0.00000000e+00, ...,\n",
       "        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],\n",
       "       [2.62253968e-02, 5.37739098e-03, 1.29560860e-03, ...,\n",
       "        1.34500220e-06, 2.69000440e-06, 2.70361504e-03],\n",
       "       [3.42515702e-02, 5.92040562e-02, 3.51409881e-03, ...,\n",
       "        6.36500137e-05, 1.27300027e-04, 3.91269990e-04],\n",
       "       [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,\n",
       "        0.00000000e+00, 0.00000000e+00, 1.71932455e-04],\n",
       "       [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,\n",
       "        8.37529941e-06, 1.67505988e-05, 0.00000000e+00],\n",
       "       [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,\n",
       "        0.00000000e+00, 0.00000000e+00, 0.00000000e+00]])"
      ]
     },
     "execution_count": 156,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f075fc74f60>]"
      ]
     },
     "execution_count": 158,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3X2YHFWdL/Dvj2BcX9b1hfHqDWCiZkWWCy7MRrxXxH2QNaC7kZXFyLoquxrj3eyyj+vVuKCCiAaiiCuBEJB3SAgQIZBJQkgmxCSTecn720wyeZnJZJKZyUySmcxk3n/3j66eqemp6jrVXd1VXfP9PE+e9FSfrj6nq+pXp06dc0pUFUREFC/nhJ0BIiIKHoM7EVEMMbgTEcUQgzsRUQwxuBMRxRCDOxFRDDG4ExHFEIM7EVEMMbgTEcXQuWF98XnnnacTJ04M6+uJiArS5s2bT6hqkVe60IL7xIkTUVVVFdbXExEVJBGpM0nHZhkiohhicCciiiEGdyKiGGJwJyKKIQZ3IqIYYnAnIoohBnciohhicKeMvbmvBUfausLOBhE5YHCnjH3jsQpc8+s3w84GETlgcKes9A4Mhp0FInLA4E5EFEMM7kREMcTgTkQUQwzuREQxxOBORBRDDO5ERDHE4E5EFENGwV1EpopIjYjUishslzSfFZFtIrJbRDiyhYgoRJ6P2RORcQDmAbgWQAOAShFZqqp7bGneDeBBAFNVtV5E3p+rDBMRkTeTmvsUALWqelBVewEsAjAtJc3NAJaoaj0AqGpzsNkkomyt2tOExzccCjsblCcmwX0CgCO2vxusZXZ/DuA9IrJWRDaLyNeDyiARBePbT1Xhzlf3eCekWPBslgEgDsvUYT1XALgGwNsAlInIJlXdN2JFIjMAzACACy+80H9uiYjIiEnNvQHABba/zwfQ6JBmhap2quoJAOsAXJa6IlVdoKrFqlpcVFSUaZ6JiMiDSXCvBDBZRCaJyHgA0wEsTUnzCoCrRORcEXk7gE8C2BtsVomIyJRns4yq9ovILAArAYwD8Jiq7haRmdb781V1r4isALADwCCAR1V1Vy4zTkRE7kza3KGqJQBKUpbNT/l7LoC5wWWNKHoqDrXhkgnvwtvHGx06RKHhCFUiQ8dOn8VND5fh/72wI+ysEHlicCcy1NnTDwCoPt4eck6IvDG4ExHFEIM7EY05g4OKikNtYWcjpxjciWjMeWzDIdz0cBlKa+I7UwqDOxGNOQdaOgEAjafOhpyT3GFwJyKKIQZ3IqIYYnAnIoohBnciohhicCciiiEGdyKiGGJwJyKKIQZ3IqIYYnAnojFLUx8YGiMM7kQ+xTgejBni9GTomGFwJzI2BiICxQaDOxFRDDG4ExHFEIM7EVEMMbgTEcUQgzsRUQwZBXcRmSoiNSJSKyKzHd7/rIicFpFt1r+fBJ9VIiIyda5XAhEZB2AegGsBNACoFJGlqronJekfVfWLOcgjERH5ZFJznwKgVlUPqmovgEUApuU2W0RElA2T4D4BwBHb3w3WslSfEpHtIrJcRP4ikNwREVFGPJtl4DwsL3UE9hYAH1LVMyJyPYCXAUwetSKRGQBmAMCFF17oM6tERGTKpObeAOAC29/nA2i0J1DVdlU9Y70uAfAWETkvdUWqukBVi1W1uKioKItsExFROibBvRLAZBGZJCLjAUwHsNSeQEQ+IJKYikdEpljrbQ06s0REZMazWUZV+0VkFoCVAMYBeExVd4vITOv9+QBuBPBdEekHcBbAdNU4T6ZJRBRtJm3uyaaWkpRl822vHwDwQLBZI4ooVluoAHCEKpGhsTAHOMUHgzsRUQwxuBPRmBXnFjYGdyIac8ZCCxuDOxFRDDG4ExHFEIM7EVEMMbgTEcUQgzsRUQwxuBMRxRCDOxFRDDG4ExHFEIM7EVEMMbgTEcUQgzuRT3Gej4Tig8GdyNBYmI+E4oPBnYgohhjciYhiiMGdiCiGGNyJiGKIwZ2IKIYY3ImIYsgouIvIVBGpEZFaEZmdJt1ficiAiNwYXBaJiHJE4ztqwTO4i8g4APMAXAfgYgBfFZGLXdLdA2Bl0JkkIgqSjIFBCyY19ykAalX1oKr2AlgEYJpDun8D8BKA5gDzR0REGTAJ7hMAHLH93WAtGyIiEwDcAGB+uhWJyAwRqRKRqpaWFr95JSIiQybB3ekCJrWh6n4AP1TVgXQrUtUFqlqsqsVFRUWmeSQiIp/ONUjTAOAC29/nA2hMSVMMYJEkGrLOA3C9iPSr6suB5JKIiHwxCe6VACaLyCQARwFMB3CzPYGqTkq+FpEnALzGwE5xE99+FRRHns0yqtoPYBYSvWD2AlisqrtFZKaIzMx1BokKzWs7GrGnsT3sbNAYZ1Jzh6qWAChJWeZ481RVv5l9toiix7T33KzntgIADs/5Qu4yQ+SBI1SJiGKIwZ2IKIYY3ImIYojBnYgohhjciYhiiMGdiCiGGNyJiGKIwZ2IKIYY3ImIYojBnYgohhjciWjMivNkcAzuRDTmiPFMQYWLwZ2IKIYY3Il8Uo3zxTzFBYM7kSHrSWNEBYHBnYgohhjciYhiiMGdiCiGGNyJiGKIwZ2IKIYY3ImIYojBnYgohoyCu4hMFZEaEakVkdkO708TkR0isk1EqkTk08FnlYiITJ3rlUBExgGYB+BaAA0AKkVkqarusSVbDWCpqqqIXApgMYCLcpFhIiLyZlJznwKgVlUPqmovgEUAptkTqOoZHR6T/Q7Ee7I1IqLIMwnuEwAcsf3dYC0bQURuEJFqAMsA/LPTikRkhtVsU9XS0pJJfomIyIBJcHeaUGNUzVxV/6CqFwH4EoC7nFakqgtUtVhVi4uKivzllIiIjJkE9wYAF9j+Ph9Ao1tiVV0H4CMicl6WeSMiogyZBPdKAJNFZJKIjAcwHcBSewIR+ahYU+aJyOUAxgNoDTqzRFGQvGxt7ugONR9E6XgGd1XtBzALwEoAewEsVtXdIjJTRGZayb4MYJeIbEOiZ81XlJNeU8zY2yc31p7AlLtXY8WuY6Hlhygdz66QAKCqJQBKUpbNt72+B8A9wWaNKLp2Hj0NANhSfwpTL/lgyLmhTMW5CsoRqkQUClVFU3s4TVtj4bkrDO5EFIpnyuvxyV+sxp7G9rCzEksM7kQUirIDJwAAh050hpyTeGJwJyKKIQZ3IqIYYnAnosCs338Cz5XXh50NgmFXSCIiE1/7fTkA4OZPXmj8GeU8gznBmjsRkU9dvf14aXMDojxWkzV3IgqVOM5NGG0/e3UPFlUewYT3vA1Xfvh9YWfHEWvuREQ+JQdfdfX2h5wTdwzuRDRmbTrYipW7j4edjZxgcCeiUIV5Q3X5ruP4ztObQ/v+XGJwJ/IpwvfQCkohtrUXEgZ3IkNjYbKpfGIXyNxicCcj/QODYWdhhNNdfTjS1hV2Nsa8f36iEp++Z01W62ANPjcY3MnTun0t+Ohty7H9yKmwszLkmvvW4qp7S8POxpi3proZDSfPhp0NcsDgTp7W1rQAACoPt4Wck2EnzvSGnQVfmtq7ccVdq1Db3BF2ViKHzTO5weBOlAfLdx5Da2cvni6rCzsrkRHH5piO7j5sqT8ZdjYAMLgTEQXmu89swd8/uBGdPeEPbmJwJ8ojNkAMC7M5JlfXDDsaEvel+gfC39IM7kR5IOxH6SqOzTNRwOBOFKC+iHUZLQRxvKEahTIZBXcRmSoiNSJSKyKzHd7/RxHZYf3bKCKXBZ9VouhbVHkk7fsc3TosjjX2KF2heQZ3ERkHYB6A6wBcDOCrInJxSrJDAK5W1UsB3AVgQdAZJSoE3b0DjssjdMzn3eCgor41ngPO3E7WUTiJm9TcpwCoVdWDqtoLYBGAafYEqrpRVZP9fzYBOD/YbBIVtigc7GH53ZpafGZuKQ60nBmxPApNF5lyq6FH6SRuEtwnALBfazZYy9z8C4DlTm+IyAwRqRKRqpaWFvNcEhWIQg5YuVJ+qBUAcPx0t+P7cWyeiQKT4O70yzvuwSLy10gE9x86va+qC1S1WFWLi4qKzHNJFAHZ1L6jVKOLmjieEKNQIpPH7DUAuMD29/kAGlMTicilAB4FcJ2qtgaTPSKKi9STYxxr7FEqkUnNvRLAZBGZJCLjAUwHsNSeQEQuBLAEwD+p6r7gs0lx0Ns/iJ5+5xuOhSCI2ncca6leCu2qZV5pLb73/Laws5E1z+Cuqv0AZgFYCWAvgMWqultEZorITCvZTwC8D8CDIrJNRKpylmMqWFf+cjU+dvuKsLORU25NNwUW33KiUE5sc1fWYMnWo1mtQyNwB92kWQaqWgKgJGXZfNvrbwH4VrBZo7hp6yysmRwpGG7NL4US7P0oqH7uRGTO69iOQIUuErr7BtAXgflX4ozBnShA9uDd0tGDMxGYHTAq7L/NRT9egVV7msLLTI4li9rR3YdfluxFb3/+p6VgcCfKkb+6+w1ce9+bI5ZF6Ko9b8Iq84kzPag4lL8HzGyuOzmq6fHXr+/Dw+sO4g9bG/KWjyQGd6IcOpYycIfNMvnz5Yc24qaHy3Kybqcbpi9UjZ5XqMeqsfcP5n/DM7gTuejuG8DzlfXB9HyIcJW9rbMXze3Oo0eDlO/wVpeH+WzcNmsUTuJGvWWIxqLfrEpcUv/Z28Zj6iUfMPpMBI5p3y6/axUA4PCcL4ScEwoSa+450N03gG89WYmDKRMlUWFJtp+2n+0LbJ2FGPztVBW/X38IpzP4TaLQ9ztfotDNk8E9B8oPteGNvc346dLdYWeFsnDuuMQ1dxDtpdFtlPGn7GAr7nptD25/eRdOGo5biFLf71yKWjEZ3IlcjDsncbQODJp3Y4vY8Z2V5yvrR/U2Sd4gbD/bh6OnzoaRrUAEdcLxczHS1N6NV7ZlN/LVDwb3DAwOKl7b0YhBlxqd6eXniTM9uPyuVdjT2B5k9igg556TODz81NzDvxgPzg9f2hlob5M4/TaeJweHwn7t0XLcumhb3sY+MLj7cKQtcff92fI6zHpuK5536Ppk57UDrK1pQVtnLx5dfzCwPFJwhmvucQpL+edVRx4LTfGLqxqG5rMfzFOBGdwNrdh1HFfdW4o11U1oau8BAJzo6Alm5WNg5y5EVmwfFdyDuFm291h7Qd9w5y7rRBxfAsD2I6fymxWwK6SxnUcTG2dPY3tgd8Lj1D4bR8krr+TWNpl/3KlSZn8CUfL96377RwCF1/0wq33WbcbMiB4I3X0D+JO3jMvswxE4+7HmngXXAQw+1xOB/cC3sdCtLaiYU3m4LbIBLB8Kteypz3wtNAzuOeS1TxfqTg8ACyvS32+Ik3TnsaMnC7fHCIUjX/UiBvcMBL1xTGrBHd19+Oh/laC0pjnYL89Qyc5jYWcBAHCqK4dzxFsn33TNcMsMfoeRJ/HoXvHUNnegucNsGoKO7gAHdkX0J/H7GED7dnZ+yHRW2fGNwT0L9t4wJ8704MaHNvqao8NPzX1f0xn0Dyr+e/V+P1mMtMZTZ9GeZZA40pa7mrPfg7u5vRv3rKjOUW5y73P3rcOUu1enTZPc57fW5/8GYb4FfWWd715XDO4ZcNpEiyrqUVV3Ek+WHTaqnHX29OOHL+50XZ/rdwewf9Qc70BpdfhXAP97zhpcd/8fw86GJ9Pf/LaXdzkuH3mSKOC2uCxFYUh+mLp68/v8YAZ3n7wOdPv7IsCSLQ2YOHvZqMvdZ8vr0DvgY+RjgDHh8/evwy1PVAa3wixEeZSj39/crGY29gJc2PeeMr357zdf9uRRaGpicDfkdYlub6Kx11AWWTceD7V0uq7Pz44Qxj6THHSxqDKaN1Grj7fnZNRfkDHHbxNPHEUh4PmRs23GG6rRlW4nTX3L7VI0m1pBOqXVzXhs/SF/K/fQZ11h1DYPdw1bX3si0O/IlEIx9f4/4puPVeTuOwyjklu6dNu64WTu5xxPddmdr3umWVhR77g8k3DnNVI710E/0/WLJK7Gbn9559Do9LTfk9nX5IxRcBeRqSJSIyK1IjLb4f2LRKRMRHpE5PvBZzOa/ATodDu4r53CY0+95YlK/Oy1PX7W6OjY6eHmkqjttE6q6k5iYFBx4kxAo4YxvH1Ng4NRo0xKou8t3u4rT0Ewma73R0t25iEn/lUcakNda6fr+9uOnBqx72Zra/1JPLOpHv/x/DZfn4vC/QXP4C4i4wDMA3AdgIsBfFVELk5J1gbg3wH8KvAcRlC6DTeizT3lvR0NpzBx9jJsrmsbEezfrGnGb1btS/ud+Zw29ZVtR/GpX67BpoOtefm+Ay1nsLPhdNbr+dXrNSj++RtYm6a76NqaZpQd8F+unQ2n0dWXedNPIY9psMsmZJmcJE+f7cPiNHM23fRwGa6eu9b1/S/N24Cr7il1fX9fUwdqmzu8MwL/VylR28QmNfcpAGpV9aCq9gJYBGCaPYGqNqtqJYDgOr8WAHub3FAtDzpiJ7a/XrevBQCwprp5aN4SAGjv7sdvDbs42o+Prt5+bDwQfPPIlrqTAIDqY/mZrfKaX7+Jv31gfUaftddCV+9tAgB883H3m8XffLwSX31kk6/vaO/uw98+sB63Lkxfe3MLXlE76PPNT/l/8OJ2/ODFHdh1NPOTfeosnva//uY36/C5+9YZrWdEv3Vr457p6Te6vxOF+wsmwX0CAPuptMFaNnY5bDi3my9D85KkvJ1NraCnP9Gl6ocv7cTNj5QbtQfG1V0BNEG5SW7T7r7EPYeapvQ1PrfjOdNa+1X3rsEND27w9ZkjbV05bcfPxYnK/vs0W5PxJffxcMmobXfJT1fikp+uHPrbK4j39IVXDpPg7rQ9MzovicgMEakSkaqWlpZMVhEae1OMW8BOvul0xk/tJuXUzNI/MOhZK9jfdAYfu30FXtl2FPuOJ4JNrueHjkItxM3ZEA6e/oHgB5OpKrpTynKk7eyowUJ1rZ1pn4B01b2l+HSaZokwue1GTk2ZQe5zuegKaXqSW7I1fw/nSGUS3BsAXGD7+3wAjZl8maouUNViVS0uKirKZBWRMmJoyoiA7pBWRgZ0px3n1ue3jagVpH4eGA5mb+xtNq4Rnj7bl9EJYPTpLBir9jThyY2HA1nXyMAQbL3S7fc9drobnQ4DUrKZTO3+N/bjoh+vGBqxW9/qXPu+eu5afPZXazP+njAE1fEg30ZUyEzSe00/kGcmwb0SwGQRmSQi4wFMB7A0t9mKnpH90tPcULV/RsT9Ut1h2bId5vO12PPgFVMuu/N1XPLTldhSf9Jo3bk+wL79VFVBPF/Wz6+QfhqFNGMarL9f2tIAADjd1YfWMz34zNzh2vfE2cvwzKa6ob8zeTh1FPg5+SVTTpy9DN9/IdGjyO2EZ7quzETnZOOXZ3BX1X4AswCsBLAXwGJV3S0iM0VkJgCIyAdEpAHA9wDcLiINIvKuXGY8CuwxcPhyUl2CrX2Qk/8Amlor1RGvzXbfv39wo6/vHFp/FKohNlWH2xyXu/2klYfbcjoS9vXdx3HpHa+j6rD7ydN0a6s6z9vyhMeVzqQfLcMtj+eur39SZud98w85Ncu8uDlx4rOf8JIGBxVzV1YH2gV2KC8RuorIhNHDOlS1BEBJyrL5ttfHkWiuGRNMg92a6mZMfv87Hd/Ldr850dGDaqvN/XRXMDW5I21dOGmbZTFZzigF9+6+Adw4f/i5nunyNjCo2NfUgX+w0vt+MIbhRtpoda10a/83WY29t1Um+4YqUFozfB+r9UwP3vfOt/pfURpLtjTg0An3PuaZEgH2N3Xg2t+Y9WKxW197AvNKD2Bf0xk88vVixzQZD2JKsw6vVaa/uudj9iLPpI13vzWqM7VN3m/78K2Lto74u9z2VPqbHy3H4QAOuqvuLcXfPbBhVC0oCgMyktI9rDr12ZS/fWPf0BOPMmG6hfwE47N9A5hXWuvwXR7TW5h/BQDg7pK9nmm6evt9jWb+3uLt+N2a0Xk3le6GalXdyZRlZvtcsomqtz/9PE3NHd2+H2uYuE/mnWbodcSacPiYPQMTZy/zTOOndmaa3u6gR/A+3NqJiee9w99KXXSm3HyNUs09nTPdI/O9PYCBUYDBAW5wUCef6rN0eyOQZlCqqvv3mUwbYF+PlznLq/FUWZ13wgwcOtGJSdb+mNGViGG6O19N3LtJN95DoZ5TGccRa+4ZSLfjqQLr9o/u5pm6fze3B9tGGGT8TW1rjFJsT/0d7TW8xhHPKtURA8Uy+i7D6Qe8gz/Q0Z2+t9Jws4zzyUIk+Bup2a4v3e+y3uEYSGf0djX7XJvVLbRvIPi9NJuaeBQqRAzuPo3sDWN7be0ICuAl6wYQRqQduaOkG2IdtmyDYi6c7R1w7JFiD+h2c1ZU45yU3/yVbcN9jo+eOoubHi5Le7/C9OA2SeV6c05GrkMTbXYZ58WPbNeYdhoOp2UhBTyn731th3dv7gK/n8rgngn7zrKwot6aK2Z4WRh32YOdCTKl5h6BWshnf1WKS+943bMGnPRcef2o7XDrouHpAx5YU4uKQ2141TrI9zS2489vW57RpFMmm9v1hKnJddjGQGT4HfmWdnbUNHMsRcGJMyMHggXxlKRsml1zgcE9Sz9ashNffqhs6Gbe712CrH1bd3T3Bb7x/7g/d1PwRuGGapPVjHXlL83bTtNfgaiVJpHo6U2Jh6dMXzA874y9B0s6XidzERl1FZEuV/mqHGT7PembJ7PbZ1o7g+/amM6s57aMWpbNoCTTE18uMbgHxN7m53bIJB8q/Wx5vfHBHoZRQdFjZ/zqgk3o6s3tFAi+KdL+xoNW54rUstbZBsqYbiGTydu8mrrE5bWb6uPDk7r5udGaL9nGr1nPbfVOhOBOhMt3HXdefyBrDweDu0+qBjVZhz2io7sfu44OH5BRi+32yceSQTFZSq8DtexgKzbUBjc98IJ1B7Bil/loXScKYMVu5wMWGO42abIddnj0urFvVzfnuET3isNtWLWnyfPzqUEsiBGrWbe5B1gFddoOfmeGvM9j2my/Pn1PKVbu9t42UcXg7iF17uc+23NP3WoNTkv7B0f2w03XxhfkQWOqx9ZP+JyUvcKen4mzl+HHDg+CTs1zc0d3xgNeflFSjZnPbEH/wGDGswN6zaXzgnXTO7kNUzflI+sO4tdWsPAK7l76BwbTXkV8+6mqodem0wY/s6keFYfaRjwdK9+8eo25fcL0Ks9+0vrO01VpUia4TeaWzeE0/80DaVfS3TeIikPOI6bDxuDuIXXu5+W7jjluZ6+AnHoiaDjpfuPOvqrBQTXqZ2+q33ZymjZvw9BMhOkGY6SW7GlbrdEtzZS7V+Ovs5zg6uZHy/Gx21dktQ4vbkHXa7i/H3e8uhvnGrfLOI9QdVp208Nl+Nx9b2abvYx5Bc3uvgH85+LtON4+skfTxT9xnhwvlb3MbjXoMCpCdre9vBM3PVwWyWm3OYjJpwO2B13bj7c9x+xNLsG1ufQOpB9559fn7x8+WW0/cgqb606i8dRZXPGh9wwnSsm+yfGTi2MsHzWiZMy1Fzno/uRN7T2uzTJJI9vcI9Zm5yp9V8iLfmx2YlYA8988OGr5uICOoyA6BBw73T1qSmYAOGV1pU3dZxpPnR2aDC4sDO4ZSA6csO97JTuH23cda14+1j9iUrAM98uVu4/j83/xgVHL7ScnAPjBiztw9NRZPHDzXw4ta7Qm2UrWikyykK4G1T8wiF+UVBusxVkua0VONffL7nwd73xrsIeG6dgBtxGqdRnOiAgAp7p6nZsBXfKUfobLYel7hIx+0y39OQLHJryvLHB/Ytadr+7G28eP88xjUJo7enCL7QlfToOu7MvS5T1fxlyzzK6jpx3PwH4kjxO34zWoG1WHT3Ti4z/JrFniO09vNkrXZF0yJ582BIyeldDk0nd97Qk0dzgPKNpwoBWPbci8H/4Mw7Jkwq1yGPQDULxqoWK7ie2UMpv8fOJnq3DFz98Y/Z0ue+p3n8nN7+3WFJlJz7HHNxzGvNIDMOmeblpB2lp/EoNpVlhme6ZwS8r8S6nzGkXBmAruJ8704Iu/W48fvLgjq/X4bV9PLPP/PQsr6v1/yEZVPW9eJUuSLnsmu+2z5fW48aEy74RpuN1byMV0rkn56lO+xWEaX7sRN0bzkCVVxX6XB0XvbzK7STvyCnPkXuJ0iNxdshd7Gkf3LIrKiOgbHtyIl7eZPTkp9SZ79EL7GAvuZ62n55g+tMKNV23BcYShjyNWAexoOIXOLPuOT1+wyfPmVfJy/Vev14x6LznZlemeW9/WNdy7wCYqB6+TfOWtw7DmrYq8RIrnKupdewGlfv2/L9w6YuqGoXS2hEu2jHzfrZ27tKZ51LJcnGDtXUX9/JymI6BTRbHmPibb3LPdDskN2eZjHnWnndrNnOXVriNd01m6feR8GeU+bkgec5ijZWHFESys8DcHzpzl1Zh59UeG/m442ZX1gK1cHjfnuHSFDJrp6rt6+3Hzo+U5zQuQvm9+6u+9dHvjqH0rkW444X++sB1fvuJ823vO6567cnQlIhc//e227rpdPpq0Mj3Zuz0X2TFtZl/hW8HV3Du6+/C71fszmgsi2WMh2+5TyY+79attdXiAsZ9gnUlgB5wnLAtDo+2pR9Me2OBr3ptXHYLIqS73B0Jn6/nKI5g4exm6esJ7Sr3dDRk+LStVfVsXBgc1zb6etpe60XcEddMw16O1/+sPo8dluPHq1eQm7C6ZTgqu5j5neTWeLa/HpKJ34IuX/k9fn01utmznCIriJRgQnXzZJ99q7ezF6mrzq5Z/Wzh62Hm6B3Rk6819ialpm1xuBgdl25H0be5B21x3Et99drPr4wXT7SpB/Ny/XG7eO+pbT3kPUMpGfZv5YLpMTzSDah7g83UiKLiae7LdvK61y/fovOFh9YrOnn78rztWYq2P5pKkHo+nvoQlIrEd/TmYWzvXgpw+ISpW7m5ybX5Jt6+0OVx5mvjRkp0ZfS7XDmfRjdTUoKpxM+gVP38jkFkovRRccE+2a81dWeM4Oq+0utmxq+OOhlNDN/sGNXGzsKO73/FGIgDUHO/APz4afl9VP6JSc49CH19KLxczfWbbuytXvB7BZzfHxxWH3aDq0DONTayvzd0srkkFF9zHpcnxrqOnccsYOhZ/AAAGCUlEQVQTlbjz1T2j3vu7BzYMDSlX1aFLT7fLsK8sKCu42lzyQc1EXhZXReP+TNRkOjr5kXWjR9im83gW4z5MFVxwTw3G9omlkjcy61rTt7HZ28ecuoOd7uobGlZMROSltMbfYwW35+EejNENVRGZCuC3AMYBeFRV56S8L9b71wPoAvBNVR09+30AUrsbfez2FZgy8b24+mNFWG89sGLjgVbUtXbife98KxaW1+PjH3zXiM+0dfbiuK3r34Nra/GJC96NoyfP4h+KL8BlP4ve/NhEFB95aHL3Du4iMg7APADXAmgAUCkiS1XV3vZxHYDJ1r9PAnjI+j9wTj2VKg63oeLwyJsZV89dm3Y93312+Nxz74rhdve7S/ZmlT8iIi/5uD9m0iwzBUCtqh5U1V4AiwBMS0kzDcBTmrAJwLtF5IMB5xVA7vvEsjmGiHIt05GwfpgE9wkA7MMUG6xlftMEYlyUx7ITEUWESXB3iqap1xQmaSAiM0SkSkSqWlr83YBIuubj78/oc0REUfG1Ky/M+XeY3FBtAHCB7e/zAaSOETdJA1VdAGABABQXF2fU6HTV5CIcnvOFTD5KRDRmmNTcKwFMFpFJIjIewHQAS1PSLAXwdUm4EsBpVc3uCcdERJQxz5q7qvaLyCwAK5HoCvmYqu4WkZnW+/MBlCDRDbIWia6Qt+Quy0RE5MWon7uqliARwO3L5tteK4B/DTZrRESUqYIboUpERN4Y3ImIYojBnYgohhjciYhiiMGdiCiGJKxn/4lIC4A6z4TOzgOQ+9nuwzcWyskyxsNYKCMQjXJ+SFWLvBKFFtyzISJVqlocdj5ybSyUk2WMh7FQRqCwyslmGSKiGGJwJyKKoUIN7gvCzkCejIVysozxMBbKCBRQOQuyzZ2IiNIr1Jo7ERGlUXDBXUSmikiNiNSKyOyw8+OXiBwWkZ0isk1Eqqxl7xWRVSKy3/r/Pbb0P7LKWiMin7ctv8JaT62I/LekPjk8j0TkMRFpFpFdtmWBlUlE3ioiz1vLy0VkYj7LZ+XBqYx3iMhRa1tuE5Hrbe8VYhkvEJFSEdkrIrtF5FZredy2pVs5Y7U9oaoF8w+JKYcPAPgwgPEAtgO4OOx8+SzDYQDnpSy7F8Bs6/VsAPdYry+2yvhWAJOsso+z3qsA8CkknoK1HMB1IZbpMwAuB7ArF2UC8H8BzLdeTwfwfETKeAeA7zukLdQyfhDA5dbrPwWwzypL3LalWzljtT0LreZu8rDuQjQNwJPW6ycBfMm2fJGq9qjqISTmy58iiYePv0tVyzSx9zxl+0zeqeo6AG0pi4Msk31dLwK4Jt9XKi5ldFOoZTymqlus1x0A9iLxLOS4bUu3cropyHIWWnDP24O4c0gBvC4im0VkhrXsf6j15Crr/+SDYt3KO8F6nbo8SoIs09BnVLUfwGkA78tZzv2ZJSI7rGabZHNFwZfRakb4SwDliPG2TCknEKPtWWjB3ehB3BH3f1T1cgDXAfhXEflMmrRu5S3k3yGTMkW1vA8B+AiATwA4BuDX1vKCLqOIvBPASwD+Q1Xb0yV1WFbI5YzV9iy04G70IO4oU9VG6/9mAH9AoqmpybrEg/V/s5XcrbwN1uvU5VESZJmGPiMi5wL4M5g3keSMqjap6oCqDgJ4BIltCRRwGUXkLUgEvGdVdYm1OHbb0qmccduehRbcTR7WHVki8g4R+dPkawB/A2AXEmX4hpXsGwBesV4vBTDduvM+CcBkABXWpXGHiFxpteN93faZqAiyTPZ13QhgjdXGGapkwLPcgMS2BAq0jFaefg9gr6reZ3srVtvSrZxx2555vXsbxD8kHsS9D4k71reFnR+fef8wEnfdtwPYncw/Em1xqwHst/5/r+0zt1llrYGtRwyAYiR2vgMAHoA1IC2kci1E4jK2D4kay78EWSYAfwLgBSRuZFUA+HBEyvg0gJ0AdiBxMH+wwMv4aSSaDnYA2Gb9uz6G29KtnLHanhyhSkQUQ4XWLENERAYY3ImIYojBnYgohhjciYhiiMGdiCiGGNyJiGKIwZ2IKIYY3ImIYuj/A5khlSr/vM41AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(H1[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f075fcb77f0>]"
      ]
     },
     "execution_count": 157,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFRZJREFUeJzt3WuMXOd93/Hvf2eXy4tISRRXkcqLSKW0DNV1EmnNxMilNgK1lBCAduOmTIoqVxBqo7YpaiAsjAYB/KJxgxZpYDks6whNiiBM0tgJ4dJRm0S13fpGypBlUTKllXzRirR4cSSKEsm9zL8v5iw5Ws1yh8vZszuPvh9gMec855lznodn+dszz5xLZCaSpLIMLHUDJEm9Z7hLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCjS4VBvesGFDbt26dak2L0l96bHHHjudmSPz1VuycN+6dStHjhxZqs1LUl+KiG91U89hGUkqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KNfusvn+Fzz55a9O0Y7pJUo489+hyff+7Mom/HcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkk1SrKW7RjuklSzqGEbhrskFaircI+InRFxLCLGImLvHHXeExGPR8TRiPhMb5spSboa8z5mLyIawEPAPcA4cDgiDmbmU211bgA+BuzMzG9HxM2L1WBJ0vy6OXLfAYxl5vOZOQEcAHbNqvMzwCcy89sAmXmyt82UJF2NbsJ9I/BC2/x4VdbubcCNEfF/IuKxiLi/Vw2UJF29eYdl6PzF7uxzeQaBu4EfB1YBX4iIL2bmM29YUcQeYA/Ali1brr61kqSudHPkPg5sbpvfBBzvUOcvMvO1zDwNfBb4vtkrysz9mTmamaMjIyMLbbMk9a2s5zT3rsL9MLA9IrZFxApgN3BwVp0/B340IgYjYjXwg8DTvW2qJJUhajjRfd5hmcyciogHgUeABvBwZh6NiAeq5fsy8+mI+AvgCaAJfDwzn1zMhkuS5tbNmDuZeQg4NKts36z53wR+s3dNkyQtlFeoSlKBDHdJKpDhLkkFMtwlqUY1nQlpuEtS3aKGm/4a7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalGWdNtIQ13SapZHXeFNNwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SaqRt/yVpELVcJq74S5JJTLcJalAXYV7ROyMiGMRMRYRezssf09EvBIRj1c/v9b7pkqSujU4X4WIaAAPAfcA48DhiDiYmU/Nqvq5zPyJRWijJOkqdXPkvgMYy8znM3MCOADsWtxmSZKuRTfhvhF4oW1+vCqb7d0R8dWI+HRE/J2etE6SClPTHX/nH5ah81k7s5v3FeC2zDwXEfcBfwZsf9OKIvYAewC2bNlylU2VpELUcM/fbo7cx4HNbfObgOPtFTLzbGaeq6YPAUMRsWH2ijJzf2aOZuboyMjINTRbknQl3YT7YWB7RGyLiBXAbuBge4WIuCWi9acoInZU6z3T68ZKkroz77BMZk5FxIPAI0ADeDgzj0bEA9XyfcAHgH8WEVPAeWB31vUsKUnSm3Qz5j4z1HJoVtm+tumPAh/tbdMkSQvlFaqSVCDDXZIKZLhLUs28K6QkaUEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklSTOu/KYrhLUs1quOOv4S5JJTLcJalAhrskFchwl6QCGe6SVCDDXZJqUufz6Qx3SapZ1HDTX8NdkgpkuEtSgQx3SSqQ4S5JBeoq3CNiZ0Qci4ixiNh7hXrviojpiPhA75ooSbpa84Z7RDSAh4B7gTuBn46IO+eo9xHgkV43UpJKUOOZkF0due8AxjLz+cycAA4AuzrU+xfAnwIne9g+SSrOcrkr5Ebghbb58arskojYCLwf2HelFUXEnog4EhFHTp06dbVtlSR1qZtw7/Q3Zvani98CfjUzp6+0oszcn5mjmTk6MjLSbRslSVdpsIs648DmtvlNwPFZdUaBA9H6rLEBuC8ipjLzz3rSSknSVekm3A8D2yNiG/AisBv4mfYKmbltZjoi/hvwKYNdkpbOvOGemVMR8SCts2AawMOZeTQiHqiWX3GcXZJUv26O3MnMQ8ChWWUdQz0zf+7amyVJuhZeoSpJNcka7/lruEtSzWo4zd1wl6QSGe6SVCDDXZIKZLhLUoEMd0kqkOEuSTVZbrf8lST10HK55a8kqc8Y7pJUIMNdkgpkuEtSgQx3SSqQ4S5JNanxppCGuyTVLWo4F9Jwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQXqKtwjYmdEHIuIsYjY22H5roh4IiIej4gjEfEjvW+qJPW3rPGmv4PzVYiIBvAQcA8wDhyOiIOZ+VRbtb8CDmZmRsQ7gT8G3r4YDZYkza+bI/cdwFhmPp+ZE8ABYFd7hcw8l3np2qs11HtPeknSLN2E+0bghbb58arsDSLi/RHxdeB/Ar/QaUURsacatjly6tSphbRXktSFbsK903Wybzoyz8xPZubbgfcBH+60oszcn5mjmTk6MjJydS2VJHWtm3AfBza3zW8Cjs9VOTM/C3xvRGy4xrZJkhaom3A/DGyPiG0RsQLYDRxsrxARfzuqO+FExF3ACuBMrxsrSerOvGfLZOZURDwIPAI0gIcz82hEPFAt3wf8JHB/REwC54F/3PYFqySJem/5O2+4A2TmIeDQrLJ9bdMfAT7S26ZJUpkGvOWvJJWjWR26Dyx+thvuklSXZjUsU8OBu+EuSXXJS0fuDstIUjEuH7kb7pJUjJkj9xpGZQx3SarLzKmQfqEqSQW5dLZMDeluuEtSTRxzl6QCOeYuSQWaufuAp0JKUkG8QlWSCuQVqpJUoEtj7g7LSFI5Lp/nbrhLUjGani0jSeW5dOReQ/Ia7pJUk6Z3hZSk8niFqiQVyCtUJalAXqEqSQXyClVJKlCz2XpdNleoRsTOiDgWEWMRsbfD8n8SEU9UP5+PiO/rfVMlqb8ly+gK1YhoAA8B9wJ3Aj8dEXfOqvYN4O9l5juBDwP7e91QSep3y+0K1R3AWGY+n5kTwAFgV3uFzPx8Zv5NNftFYFNvmylJ/W+5XaG6EXihbX68KpvLLwKf7rQgIvZExJGIOHLq1KnuWylJBVhuV6h2+iOTHcqIiPfSCvdf7bQ8M/dn5mhmjo6MjHTfSkkqQLPGu0IOdlFnHNjcNr8JOD67UkS8E/g4cG9mnulN8ySpHM1lNuZ+GNgeEdsiYgWwGzjYXiEitgCfAP5pZj7T+2ZKUv+r8wrVeY/cM3MqIh4EHgEawMOZeTQiHqiW7wN+DbgJ+Fj1cWMqM0cXr9mS1H/qvEK1m2EZMvMQcGhW2b626V8Cfqm3TZOksjSbXqEqScWZGXOvY1zGcJekmsxcobpcvlCVJPXAcrtCVZLUA5fPc1/8bRnuklSTy0fui78tw12SajJd4xWqhrsk1WRyqnVD96Eabi5juEtSTSanW0fuKwYNd0kqxsT0NABDDYdlJKkYk1MeuUtScS5Ot8bcDXdJKshE9YXqcKOx6Nsy3CWpJjPh7pG7JBVkshqW8QtVSSrI+clpBgeCwYZH7pJUjPMT06xasfjj7WC4S1JtXp+YYrXhLklleX1imtUrunoA3jUz3CWpJucnplk15JG7JBWldeRuuEtSUV6fXGZfqEbEzog4FhFjEbG3w/K3R8QXIuJiRHyw982UpP53vsYvVOcd2Y+IBvAQcA8wDhyOiIOZ+VRbte8C/xJ436K0UpIKsNy+UN0BjGXm85k5ARwAdrVXyMyTmXkYmFyENkpSEc6cm2D9mhW1bKubcN8IvNA2P16VSZK6NDHV5PzkNDesGqple92Ee6ebIORCNhYReyLiSEQcOXXq1EJWIUl96dzFKQDWrlw+wzLjwOa2+U3A8YVsLDP3Z+ZoZo6OjIwsZBWS1JfOXWiF+3Url8+R+2Fge0Rsi4gVwG7g4OI2S5LKcurcRQDWr6kn3Of9fJCZUxHxIPAI0AAezsyjEfFAtXxfRNwCHAHWAc2I+BXgzsw8u4htl6S+cfzl8wBsvGF1LdvravAnMw8Bh2aV7Wub/g6t4RpJUgcz4X7rDStr2Z5XqEpSDU68coG1w4OsW0Zj7pKka3T85fPccn09R+1guEtSLV46e8Fwl6SSZCYvvnyeWw13SSrHd1+b4PS5Ce64ZV1t2zTcJWmRffPMawBs21DPaZBguEvSonvuZCvcb99wXW3bNNwlaZGNnTrHUCPYvN4jd0kqxqNfP8ndt91IY6DTfRgXh+EuSYto7OQ5nj15jnvfcWut2zXcJWkRff650wC8946ba92u4S5Ji+iRo99hy/rVbF6/qtbtGu6StEi+efo1/t/YGX7yrk1E1DfeDoa7JC2a3/7rZ1k5NMBPvav+m+Ya7pK0CF46e4FPPXGCf3T3Zm69vt4hGTDcJannppvJB//kqwwE/PwPb12SNtTzpFZJeovITP79oaf53LOn+Y1/+He5faS+q1LbGe6S1CMvnb3Ahz75Nf7y6ZPc/+7b2L1jy5K1xXCXpGt09sIk/+Uzz/G7//cbTDeTf/cTd/ILSzQcM8Nwl6QFyEyeOnGWQ187wR986du8/Poku77/b/Fv7rmDLTfVdw+ZuRjuktSlsxcm+cJzZ3j06yd59NhJXjp7kYFoXX36r+95G+/YeP1SN/ESw12SZpmcbnL85fN84/RrPHXiLEdfPMuTx1/hW2deB2Dt8CA/+rYNvPeOm3nPHTczsnZ4iVv8Zl2Fe0TsBP4z0AA+npm/MWt5VMvvA14Hfi4zv9LjtkrSgjWbyWsTU5y9MMWZcxcvPR3pzLmLnD53kZOvXuTEKxd48W/O852zF5hu5qX3blm/mndsXMdPjW7mB7bcwLu2rmeosbzPJJ833COiATwE3AOMA4cj4mBmPtVW7V5ge/Xzg8DvVK+SNKfMZKqZTEw1mZxuMjHV5OLMdDU/Od0quzjZ5PzkNBenprk4WZVNTXNhssmFyWnOT07z+sVpzk1Mce7CFOcuTvHqhUlevVDNT0yR2bkdw4MD3LxumFvWrWTHtvVsvGEVW25azW3rV/P2W9Zx/eqhev9heqCbI/cdwFhmPg8QEQeAXUB7uO8Cfj8zE/hiRNwQEbdm5omet1haAplJM6GZSTOTTMi2+WYCs+aTpNlslU03L5fPTF8qa8L0zPuarfLpahuXp5PpJm3r6fD+qm4zubSey3W5vP7qtZm0TV+u8+ZtcOl9081O7UqmqzrNKqynpptMNZPJ6fbpVlBPTbfC/GI1P1fgXo3BgWDVUIM1w4OsHm6wdniQtSuHGLlumLUrB7luZWt+7fAg61YNsn7NMOvXDHHTmmE2rB1mzYpG7fd+WWzdhPtG4IW2+XHefFTeqc5GoOfh/plnTvHhTz31hrLs8Nsx5+9LhwWd6na7zk6/mNmh5ly/wN3+Yve6PXPX7a5ep5rdr2+O9nS57a7391W0J7kc1jkTzG3zJWoMBAMBAxHVdGu+MdCajwgaMTNdlUcw0Ol9A0GjKhtsBKtXDDLYCAYHBhhqBIONAYYGWstWDA4w1BhgxeAAw43L028on5luDDA02HpdOTTAqhUNVg42GB4aYHiwwYrBAVYODjC4zIdIlkI34d7pz9nsX/du6hARe4A9AFu2LOzk/uuGB7nje9a+eUGHFsz1d7jTX+hOdTv9Ie9cr7v1zdWg6LCg+20vfH1zly9hezqus7sjqmvZdgADVYgFrfCKgEa0JmbCrFXeqndpnjfODwy0thoRVQAyT1C26gxcqt/+ShWcMyF6+b1R1WlE2zpn1j/z/pl1Ve+9PF3WUarerJtwHwc2t81vAo4voA6ZuR/YDzA6Orqg46G7b7uRu2+7cSFvlaS3jG4+yxwGtkfEtohYAewGDs6qcxC4P1p+CHjF8XZJWjrzHrln5lREPAg8QutUyIcz82hEPFAt3wcconUa5BitUyF/fvGaLEmaT1fnuWfmIVoB3l62r206gV/ubdMkSQvlV8ySVCDDXZIKZLhLUoEMd0kqkOEuSQWKuS4FX/QNR5wCvrXAt28ATvewOcvVW6Gf9rEMb4U+wvLo522ZOTJfpSUL92sREUcyc3Sp27HY3gr9tI9leCv0Efqrnw7LSFKBDHdJKlC/hvv+pW5ATd4K/bSPZXgr9BH6qJ99OeYuSbqyfj1ylyRdQd+Fe0TsjIhjETEWEXuXuj1XKyK+GRFfi4jHI+JIVbY+Iv53RDxbvd7YVv/fVn09FhH/oK387mo9YxHx27GEzwiLiIcj4mREPNlW1rM+RcRwRPxRVf6liNhaZ/+qNnTq469HxIvVvnw8Iu5rW9aPfdwcEY9GxNMRcTQi/lVVXtq+nKufRe3P1iPG+uSH1i2HnwNuB1YAXwXuXOp2XWUfvglsmFX2H4C91fRe4CPV9J1VH4eBbVXfG9WyLwPvpvUQoU8D9y5hn34MuAt4cjH6BPxzYF81vRv4o2XSx18HPtihbr/28Vbgrmp6LfBM1ZfS9uVc/Sxqf/bbkfulh3Vn5gQw87DufrcL+L1q+veA97WVH8jMi5n5DVr3y98REbcC6zLzC9n67fn9tvfULjM/C3x3VnEv+9S+rv8B/Hjdn1Tm6ONc+rWPJzLzK9X0q8DTtJ6FXNq+nKufc+nLfvZbuM/1IO5+ksD/iojHovVMWYDvyerJVdXrzVX5XP3dWE3PLl9OetmnS+/JzCngFeCmRWv51XkwIp6ohm1mhiv6vo/VMMIPAF+i4H05q59Q0P7st3Dv6kHcy9wPZ+ZdwL3AL0fEj12h7lz97ed/h4X0abn293eA7wW+HzgB/MeqvK/7GBHXAX8K/Epmnr1S1Q5l/dzPovZnv4V7Vw/iXs4y83j1ehL4JK2hppeqj3hUryer6nP1d7yanl2+nPSyT5feExGDwPV0P0SyaDLzpcyczswm8F9p7Uvo4z5GxBCtwPuDzPxEVVzcvuzUz9L2Z7+FezcP6162ImJNRKydmQb+PvAkrT78bFXtZ4E/r6YPArurb963AduBL1cfjV+NiB+qxvHub3vPctHLPrWv6wPAX1djnEtqJvAq76e1L6FP+1i16XeBpzPzP7UtKmpfztXP0vZnrd/e9uKH1oO4n6H1jfWHlro9V9n222l96/5V4OhM+2mNxf0V8Gz1ur7tPR+q+nqMtjNigFFav3zPAR+luiBtifr1h7Q+xk7SOmL5xV72CVgJ/AmtL7K+DNy+TPr434GvAU/Q+s98a5/38UdoDR08ATxe/dxX4L6cq59F7U+vUJWkAvXbsIwkqQuGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBfr/uGFwEeyK3LoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(sorted(H1[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f075fa2be10>]"
      ]
     },
     "execution_count": 168,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl0Ved57/Hvo3lEEmgCIWYxHIzxIOMpxtgxBCe2sROncdLbKbmXOolz07RdjXvblaRJM7W9t5k8xNdNp7Slvqnt2A42ECe2SbBjhBnMjMAMEiAJMUlMmp77xznIB1k2R3DEPsPvs5aWzn7Pu6VnM/z2e949mbsjIiLpIyPoAkRE5NJS8IuIpBkFv4hImlHwi4ikGQW/iEiaUfCLiKQZBb+ISJpR8IuIpBkFv4hImskKuoDBlJeX+4QJE4IuQ0QkaaxZs+aQu1fE0jchg3/ChAk0NDQEXYaISNIwsz2x9o1pqsfMFprZNjNrNLMH36XPPDNbZ2abzOzlqPbdZvZm5D2luYhIwM474jezTOAhYD7QBKw2s2fcfXNUn1LgYWChu+81s8oBP+YWdz8Ux7pFROQCxTLinwM0uvsud+8ClgCLBvT5BPCku+8FcPfW+JYpIiLxEkvw1wD7opabIm3RpgJlZvaSma0xs9+Nes+B5ZH2xRdXroiIXKxYDu7aIG0Db+KfBVwNvB/IB141s9fcfTtwo7vvj0z/rDCzre7+yjt+SXinsBhg3LhxQ9kGEREZglhG/E1AbdTyWGD/IH1ecPcTkbn8V4DZAO6+P/K9FXiK8NTRO7j7Y+5e7+71FRUxnZEkIiIXIJbgXw3UmdlEM8sB7gOeGdDnp8BNZpZlZgXAtcAWMys0s2IAMysEFgAb41e+iIgM1XmD3917gAeAZcAW4Al332Rm95vZ/ZE+W4AXgA3A68Dj7r4RqAJ+ZWbrI+0/c/cXhmdT4Hsv7mD9vqPD9eNFRFKCJeIzd+vr632oF3AdO9nNB77zCm2dZ/jDuZP4/G115GZlDlOFIiKJxczWuHt9LH1T5l49JQXZLPvCXD5yVQ0Pv7STO773K43+RUQGkTLBD1CSn83f3Dubf/yDa+g43cM9D/+ab7+wldPdvUGXJiKSMFIq+M+6ZVoly/94LvdePZZHXtrJnd//Fes0+hcRAVI0+AFG5IVH///0B9fQeaaHDz/8a771vEb/IiIpG/xnzZtWybIvzOWjV9fy6Ms7uUOjfxFJcykf/BAe/X/73sv550/O4YRG/yKS5tIi+M+6eWoFy74wl9+qD4/+P/S9lazdeyToskRELqm0Cn4Ij/6/9ZHL+ZdPzuFUVy8feWQV33x+i0b/IpI20i74z5o7tYIXvjCXj11Tyw9f3sWHvreSNzT6F5E0kLbBD+HR/zc//Pbo/95HVvHNpRr9i0hqS+vgP2tuZO7/Y9eM44evaPQvIqlNwR9RnJfNNz88i3/91BxOd/dx7yOr+IZG/yKSghT8A9xUV8ELf3QT980Zx2Ov7OKD31vJmj0a/YtI6lDwD6I4L5tv3DOLH3/qWs5093Hvo6v4+s82a/QvIilBwf8e3ldXzrIvzOUTc8bxf1e+xQe/u5I1ew4HXZaIyEWJKfjNbKGZbTOzRjN78F36zDOzdWa2ycxeHsq6iawoN4uv3zOLf/vv13Kmp497H31Vo38RSWrnDX4zywQeAm4HQsDHzSw0oE8p8DBwl7vPBD4a67rJ4sYp4dH/b1/79ui/YbdG/yKSfGIZ8c8BGt19l7t3AUuARQP6fAJ40t33Qv+D1WNdN2kU5Wbx13fP4t+jRv/3PrKKH7+2hyMnuoIuT0QkJrEEfw2wL2q5KdIWbSpQZmYvmdkaM/vdIaybdG6IjP7/bOE0jp3q5i+f3sicb/yc//EvDSx984CmgUQkoWXF0McGaRv4oN4s4Grg/UA+8KqZvRbjuuFfYrYYWAwwbty4GMoKVlFuFp+ZN4VP3zyZzQeO8/TaZn66bj8rNrdQnJfFBy8bzd1X1nDtxJFkZAz2xyAiEoxYgr8JqI1aHgvsH6TPIXc/AZwws1eA2TGuC4C7PwY8BuGHrcdUfQIwM2aOKWHmmBIevH0Gr+5s56m1zTy3YT//2bCP0SV5LLqihnuurGFadXHQ5YqIYO7vnbFmlgVsJzyabwZWA59w901RfWYAPwA+AOQArwP3AVvPt+5g6uvrvaGh4QI3KTGc6uplxZYWnl7bzMvb2+jtc2aMHsE9V47hrtk1VJfkBV2iiKQQM1vj7vWx9D3viN/de8zsAWAZkAn8yN03mdn9kfcfdfctZvYCsAHoAx53942RYt6x7gVtVZLJz8nkrtljuGv2GNo7z/DchgM8tbaZbyzdyjef38oNk0dx9xU1LLysmuK87KDLFZE0ct4RfxBSYcT/bt46dIKn1zbz9Lpm9rSfJDcrg/mhKu65soa5UyvIztQ1dSIydEMZ8Sv4A+LurN13lKfXNvPs+v0cOdnNyMIc7rg8fFD4ytpSzHRQWERio+BPMt29fbyyvY2n1jazYnMLZ3r6GD+qgLuvqOHuK2uYWF4YdIkikuAU/Ems43Q3L2w8yNPrmlm1sx13uKK2lHuurOH2WdVUFuugsIi8k4I/RRw8dppn1jfz1Nr9bDlwHIAJowq4anwZ9eNHUj+hjCkVRbpOQEQU/Klo68HjvLytjYY9R3hjzxHaI7eIGJGXxVXjy7h6XBlXTyjjitpSCnJiuTxDRFJJXE/nlMQwvXoE06tH8IeEDwzvbj/Jmj1HWLPnMGv2HOGlbW0AZGYYodEjuHp8WeSTQRljSvODLV5EEopG/Cni2Mlu3th3hDW7j7BmzxHW7TvKqcg9g0aX5HH1+DKujkwRzRhdTJZOGxVJKRrxp6GSgmxumVbJLdMqgfCZQlsPdNAQ+USwZs8RnttwAID87EyuqC2lfkL4U8FV48ooyddFZCLpQiP+NLL/6Kn+YwQNew6z5UAHvX3hv/+pVUWRTwUjuXp8GRNGFeg6ApEkooO7EpMTZ3pYv+8oa/YcCe8Q9h6h43QPAOVFOfzNvZdz6/SqgKsUkVhoqkdiUpibxQ1TyrlhSjkAfX3OjtZO1uw5wkO/bOThX+5U8IukIAW/9MvIMKZVFzOtupjWjtN898UdHOo8Q3lRbtCliUgc6dQOGdSCUDXu8OKWlqBLEZE4U/DLoGaMLqamNJ/lmxT8IqlGwS+DMjMWzKxiZeMhTpzpCbocEYmjmILfzBaa2TYzazSzBwd5f56ZHTOzdZGvL0W9t9vM3oy061SdJLIgVE1XTx8rd7QFXYqIxNF5D+6aWSbwEDCf8DN0V5vZM+6+eUDXle5+x7v8mFvc/dDFlSqX2jUTyigtyGb5phYWXjY66HJEJE5iGfHPARrdfZe7dwFLgEXDW5YkgqzMDG6dXsmLW1vp6e0LuhwRiZNYgr8G2Be13BRpG+h6M1tvZs+b2cyodgeWm9kaM1t8EbVKABaEqjl2qpvXdx8OuhQRiZNYgn+w6/YHXu77BjDe3WcD3weejnrvRne/Crgd+KyZzR30l5gtNrMGM2toa9OccqKYO7Wc3KwMnd0jkkJiCf4moDZqeSywP7qDux93987I66VAtpmVR5b3R763Ak8Rnjp6B3d/zN3r3b2+oqJiyBsiw6MgJ4ub6spZsbmFRLy9h4gMXSzBvxqoM7OJZpYD3Ac8E93BzKotckcvM5sT+bntZlZoZsWR9kJgAbAxnhsgw29BqJrmo6fYHHkKmIgkt/Oe1ePuPWb2ALAMyAR+5O6bzOz+yPuPAvcCnzazHuAUcJ+7u5lVAU9F9glZwL+7+wvDtC0yTG6dUYkZLN/UwswxJUGXIyIXSXfnlJh89NFVnDjTy9LP3xR0KSIyiKHcnVNX7kpM5oeq2HzgOPsOnwy6FBG5SAp+icn8UDUAKzbr7B6RZKfgl5hMLC9kalWRgl8kBSj4JWbzQ1W8vvswR092BV2KiFwEBb/EbEGomt4+5xdbW4MuRUQugoJfYjarpoTqEXm6ilckySn4JWYZGcZtoUpe3t7G6e7eoMsRkQuk4JchWRCq5lR3L79u1F22RZKVgl+G5LpJoyjOzdJ0j0gSU/DLkORkZTBveiU/39JCb1/iXfUtIuen4JchWxCqov1EF2v3Hgm6FBG5AAp+GbJ50yrIzjSW62IukaSk4JchK87L5vrJ5SzfdFD36BdJQgp+uSALQlXsbj9JY2tn0KWIyBAp+OWCzA9VAWi6RyQJxRT8ZrbQzLaZWaOZPTjI+/PM7JiZrYt8fSnWdSU5VY3IY3ZtKcs3HQy6FBEZovMGv5llAg8Rflh6CPi4mYUG6brS3a+IfH11iOtKEloQqmJ90zEOHjsddCkiMgSxjPjnAI3uvsvdu4AlwKIYf/7FrCsJ7gMzw9M9K7ZoukckmcQS/DXAvqjlpkjbQNeb2Xoze97MZg5xXUlCkyuKmFheqHv0iySZWILfBmkbeA7fG8B4d58NfB94egjrhjuaLTazBjNraGtri6EsCZqZsSBUxas7D3H8dHfQ5YhIjGIJ/iagNmp5LLA/uoO7H3f3zsjrpUC2mZXHsm7Uz3jM3evdvb6iomIImyBBmh+qorvXeWmbdtYiySKW4F8N1JnZRDPLAe4DnonuYGbVZmaR13MiP7c9lnUluV05rozyohxN94gkkazzdXD3HjN7AFgGZAI/cvdNZnZ/5P1HgXuBT5tZD3AKuM/Dl3QOuu4wbYsEIDPDuG1GFc9tOMCZnl5yszKDLklEzuO8wQ/90zdLB7Q9GvX6B8APYl1XUsv8UBVLVu/jtV2HuXmqpulEEp2u3JWLduOUcgpyMlmxWRdziSQDBb9ctLzsTG6eWsGKzS306R79IglPwS9xMT9URcvxM2xoPhZ0KSJyHgp+iYtbp1eSmWGa7hFJAgp+iYvSghyunThSz+IVSQIKfomb+aEqdrR28tahE0GXIiLvQcEvcXP2Hv2a7hFJbAp+iZuxZQXMHDNC0z0iCU7BL3E1P1TFmr1HaOs4E3QpIvIuFPwSVwtC1bjDL7Zq1C+SqBT8ElczRhcztixf0z0iCUzBL3FlZswPVbGy8RAnzvQEXY6IDELBL3G3IFRNV08fK3foHv0iiUjBL3F3zYQySguyNd0jkqAU/BJ3WZkZ3Dq9khe3ttLT2xd0OSIygIJfhsWCUDXHTnXz+u7DQZciIgPEFPxmttDMtplZo5k9+B79rjGzXjO7N6ptt5m9aWbrzKwhHkVL4ps7tZzcrAxN94gkoPMGv5llAg8BtwMh4ONmFnqXft8m/JjFgW5x9yvcvf4i65UkUZCTxU115azY3EL4KZwikihiGfHPARrdfZe7dwFLgEWD9Psc8F9AaxzrkyS2IFRN89FTbD5wPOhSRCRKLMFfA+yLWm6KtPUzsxrgHuBR3smB5Wa2xswWv9svMbPFZtZgZg1tbToNMBXcOqMSMzTdI5JgYgl+G6Rt4Gf37wBfdPfeQfre6O5XEZ4q+qyZzR3sl7j7Y+5e7+71FRV6YHcqKC/KpX58GSs2K/hFEkkswd8E1EYtjwX2D+hTDywxs93AvcDDZnY3gLvvj3xvBZ4iPHUkaWJBqJrNB46z7/DJoEsRkYhYgn81UGdmE80sB7gPeCa6g7tPdPcJ7j4B+AnwGXd/2swKzawYwMwKgQXAxrhugSS0t+/Rr1G/SKI4b/C7ew/wAOGzdbYAT7j7JjO738zuP8/qVcCvzGw98DrwM3d/4WKLluQxobyQqVVFCn6RBJIVSyd3XwosHdA22IFc3P33o17vAmZfRH2SAhaEqnnk5Z0cPdlFaUFO0OWIpD1duSvDbn6oit4+5xdbdaavSCJQ8Muwm1VTQvWIPJ3WKZIgFPwy7DIywvfof3l7G6e7BzvjV0QuJQW/XBLzQ1Wc6u7l142Hgi5FJO0p+OWSuG7SKIpzszTdI5IAFPxySeRkZXDL9Ep+vqWF3j7dtE0kSAp+uWTmh6poP9HF2r1Hgi5FJK0p+OWSmTetguxMY7ku5hIJlIJfLpnivGxumFzO8k0HdY9+kQAp+OWSmh+qYnf7SRpbO4MuRSRtKfjlkjp70zZN94gER8Evl1TViDxm15ayfNPBoEsRSVsKfrnkFoSqWN90jIPHTgddikhaUvDLJfeBmZF79G/RdI9IEBT8cslNrihiYnmhpntEAhJT8JvZQjPbZmaNZvbge/S7xsx6zezeoa4r6cPMWBCq4rVd7Rw/3R10OSJp57zBb2aZwEOEH5YeAj5uZqF36fdtwk/qGtK6kn4WzKyiu9d5aVtb0KWIpJ1YRvxzgEZ33+XuXcASYNEg/T4H/BfQegHrSpq5oraM8qIcPZJRJACxBH8NsC9quSnS1s/MaoB7gIGPYzzvulE/Y7GZNZhZQ1ubRoGpLjPDuG1GFb/c2sqZHt2jX+RSiiX4bZC2gdfbfwf4orsP/B8cy7rhRvfH3L3e3esrKipiKEuS3YKZVXSe6eG1XYeDLkUkrcTysPUmoDZqeSywf0CfemCJmQGUAx80s54Y15U0dcPkcgpyMlmx+SA3T9XOXuRSiWXEvxqoM7OJZpYD3Ac8E93B3Se6+wR3nwD8BPiMuz8dy7qSvvKyM7l5agUrNrfQp3v0i1wy5w1+d+8BHiB8ts4W4Al332Rm95vZ/Rey7sWXLaliwcwqWo6fYUPzsaBLEUkbsUz14O5LgaUD2gYeyD3b/vvnW1fkrFumVZKZYazYfJArakuDLkckLejKXQlUaUEO104cqWfxilxCCn4J3IJQFTtaO9l2sCPoUkTSgoJfAnfn7DEU52Xxzee3BF2KSFpQ8EvgRhXl8vn31/HStjZ+sVVTPiLDTcEvCeH3bpjA5IpCvvbcFrp6+oIuRySlKfglIWRnZvClO2fy1qET/OOv3wq6HJGUpuCXhHHz1AreP72S7/+ikdYOPZ1LZLgo+CWh/OUdIc709PK3L2wLuhSRlKXgl4QysbyQT75vIv9vTRPr9x0NuhyRlKTgl4TzuVvrqCjO5SvPbtI9fESGgYJfEk5RbhZfXDidtXuP8vS65qDLEUk5Cn5JSB++sobZtaV86/mtdJ7pCbockZSi4JeElJFhfOXOEK0dZ3jol41BlyOSUhT8krCuHFfGh6+q4R9WvsXuQyeCLkckZSj4JaE9uHA62ZnG15fqPj4i8RJT8JvZQjPbZmaNZvbgIO8vMrMNZrYu8sD090W9t9vM3jz7XjyLl9RXOSKPB26tY8XmFlbuaAu6HJGUcN7gN7NM4CHgdiAEfNzMQgO6vQjMdvcrgE8Cjw94/xZ3v8Ld6+NQs6SZT75vAhNGFfBXz26mu1f38RG5WLGM+OcAje6+y927gCXAougO7t7p7mdPuC4EdPK1xE1uViZ/+aEQja2d/Oure4IuRyTpxRL8NcC+qOWmSNs5zOweM9sK/IzwqP8sB5ab2RozW3wxxUr6ev+MSm6qK+fvf76d9s4zQZcjktRiCX4bpO0dI3p3f8rdpwN3A1+LeutGd7+K8FTRZ81s7qC/xGxx5PhAQ1ub5nLlXGbGl+8Mcaqrl79bvj3ockSSWizB3wTURi2PBfa/W2d3fwWYbGblkeX9ke+twFOEp44GW+8xd6939/qKiooYy5d0MqWymN+9fgJLVu9lY/OxoMsRSVqxBP9qoM7MJppZDnAf8Ex0BzObYmYWeX0VkAO0m1mhmRVH2guBBcDGeG6ApJfP31bHyIIcvvrsZt4+rCQiQ3He4Hf3HuABYBmwBXjC3TeZ2f1mdn+k20eAjWa2jvAZQB+LHOytAn5lZuuB14GfufsLw7Ehkh5K8rP50w9M4/Xdh3luw4GgyxFJSpaIo6b6+npvaNAp/zK43j7nrh/8iiMnunjxT+aRn5MZdEkigTOzNbGeMq8rdyXpZGYYX75zJvuPneaRl3cGXY5I0lHwS1KaM3Ekd84eww9f3knTkZNBlyOSVBT8krT+/PbpmME3dB8fkSFR8EvSGlOaz2fmTWHpmwdZtfNQ0OWIJA0FvyS1xXMnUVOaz1ef3UyP7uMjEhMFvyS1vOxM/uJDM9h6sIP/WL3v/CuIiIJfkt/tl1Vz3aSR/O/l2zh6sivockQSnoJfkl74Pj4zOX6qm79fofv4iJyPgl9SwozRI/jta8fz49/sZdvBjqDLEUloCn5JGX88fypFuVn81bObdB8fkfeg4JeUUVaYw58smMqqne0s29QSdDkiCUvBLynlE3PGMa2qmK8v3czp7t6gyxFJSAp+SSlZmRl8+c4Q+w6f4vGVu4IuRyQhKfgl5dwwpZyFM6t56Jc7OXDsVNDliCQcBb+kpL/40Ax63fnW81uDLkUk4cQU/Ga20My2mVmjmT04yPuLzGyDma2LPDf3fbGuKzIcakcW8IdzJ/HTdftp2H046HJEEsp5g9/MMgk/Vet2IAR83MxCA7q9CMx29yuATwKPD2FdkWHx6XmTqR6Rx1ee3URvn07vFDkrlhH/HKDR3Xe5exewBFgU3cHdO/3tE6cLAY91XZHhUpCTxZ9/cDobm4/zkzW6j4/IWbEEfw0Q/b+mKdJ2DjO7x8y2Aj8jPOqPeV2R4XLX7DHUjy/jb5dt4/jp7qDLEUkIsQS/DdL2js/N7v6Uu08H7ga+NpR1AcxsceT4QENbW1sMZYmcn5nxlbtm0n6ii+/9fEfQ5YgkhFiCvwmojVoeC+x/t87u/gow2czKh7Kuuz/m7vXuXl9RURFDWSKxuaymhI/V1/JPq3bT2NoZdDkigYsl+FcDdWY20cxygPuAZ6I7mNkUM7PI66uAHKA9lnVFLoU//cA08rMz+dpzm3UfH0l75w1+d+8BHgCWAVuAJ9x9k5ndb2b3R7p9BNhoZusIn8XzMQ8bdN3h2BCR91JelMvnb6vj5e1t/GJra9DliATKEnH0U19f7w0NDUGXISmmq6ePhd99hb4+Z9kX5pKblRl0SSJxY2Zr3L0+lr66clfSRk5WBl+6I8Tu9pN89+c7ONR5JuiSRAKRFXQBIpfSvGmVzA9V8fBLO3n4pZ2MLMyhrrKIadXF1FUVM7WyiKlVxZQV5gRdqsiwUfBL2nnoE1fxm7fa2d7SyY6WDra1dPDkG810nunp71NRnMvUqvBOIPxVRF1VMSPysgOsXCQ+FPySdnKyMriproKb6t4+bdjdOXDsNNtaOtjR0sH2lk62t3Sw5PV9nIq6r//okjzqqoqZFtkRTK0qpq6yiMJc/VeS5KF/rSKEL/QaU5rPmNJ8bplW2d/e1+c0Hz3FtoMdbG/tYEdLJ9sOdvDarna6evr6+40tyw/vBKqKmBbZIUypLCIvWweQJfEo+EXeQ0aGUTuygNqRBdwWqupv7+1z9h4+ybaDkU8IreFpo5U72ujuDZ8pZwbjRxZQV1XMFbWl3DB5FLNqSsjK1DkVEiwFv8gFyMwwJpYXMrG8kIWXVfe3d/f2saf9BNsjnwx2tHaw9WAHKzaHnwFclJvFnIkjuWHyKK6bNIrQ6BFkZAx2ZxOR4aPgF4mj7MwMplQWM6WymA/OGt3ffqjzDK/taufVneGvsxeRlRZkc+3EkdwwuZzrJ4+irrKIyEXwIsNGwS9yCZQX5XLH5WO44/IxABw8dppXdx3i1Z3trNrZzrJNLf39rpv09o5gwqgC7Qgk7nTlrkgC2Hf4ZGQncIhXd7XTcjx8cdnokjyunzSK6yeHv8aWFQRcqSSqoVy5qxG/SAI4ewD5t66pxd3ZdehE/7TQS9vbeHJtMwDjRhZw/aRR3DBlFNdPGkXliLyAK5dkpBG/SILr63O2t3b0Twv9Zlc7x0+HLzabXFHYPy103aRRjNQVx2lrKCN+Bb9IkuntczbvP86ruw6xamc7q986zImu8EVm06uLuW7SKCZVFFJTms/YsgJqyvIp0gVmKU/BL5JGunv72NB0jNd2hY8RvLHn6DlXGwOUFWSHdwKl+YwtO/tVwNiR+dSU5lOsW1EkPQW/SBpzdw51dtF05CRNR07RdOQUzUffft105CSnu/vOWackP/ucHcLbO4jwzkH3KEp8cT+4a2YLge8CmcDj7v6tAe//NvDFyGIn8Gl3Xx95bzfQAfQCPbEWJiIXxsyoKM6lojiXK8eVveN9d6f9RBfNUTuCs993tZ3gle2H3vGJYUReFjVlBed+WigLf1oYN6pAO4Ykc97gN7NMwk/Vmk/4GbqrzewZd98c1e0t4GZ3P2JmtwOPAddGvX+Lux+KY90icoHMjPKiXMqLcpldW/qO992dIye7z9khnN1J7G0/yarGQ/3HFM6qGpFLXWX4/kR1VUX9r3WwOTHFMuKfAzS6+y4AM1sCLAL6g9/dV0X1f43wQ9VFJAmZGSMLcxhZmMPlYwffMRw92d0/hfTWoZM0tnbS2NrBEw37OBm1UxhVmNO/M5hSEb6jaV1lERXFubowLUCxBH8NsC9quYlzR/MDfQp4PmrZgeVm5sAP3f2xIVcpIgnDzCgrzKGsMIdZY0vOea+vzzlw/DQ7WjoiO4NOdrR28sy6/f2noEJ46qiuqjiyMyiK7ByKGVOSpx3CJRBL8A/2tzDoEWEzu4Vw8L8vqvlGd99vZpXACjPb6u6vDLLuYmAxwLhx42IoS0QSTUaGUVManvufF3V7a3enrfMMjS3hHcGO1vCO4cWtLfxnw9vjyoKcTKZURnYEleFPB1Mqi6gdWUCmbmYXN7EEfxNQG7U8Ftg/sJOZXQ48Dtzu7u1n2919f+R7q5k9RXjq6B3BH/kk8BiEz+oZwjaISIIzMyqL86gszuOGKeXnvHf4RFfUp4PwDmFVYztPvtHc3ycnK4PJFeGdQGj0CC4fW8JlNSWU5Oug8oWIJfhXA3VmNhFoBu4DPhHdwczGAU8Cv+Pu26PaC4EMd++IvF4AfDVexYtI8htZmMOciSOZM3HkOe3HT3ezMzJVdHbHsHbvEZ5d//a4c8KoAmaNLeXymhIuH1vCzJoSXawWg/P+Cbl7j5k9ACwjfDrnj9x9k5ndH3n/UeBLwCjg4cj83NnTNquApyJtWcC/u/sLw7IlIpJSRuRlc+W4snecknr0ZBcbm4+zvukobzYd4409b+8MzGByRRGX15Qwa2x4ZxAaXUJ+jp6EFk0XcInuhJPqAAAFAElEQVRI0jvUeYY3m4/xZtMxNjQdY0PTUVo7wnc4zTCYWlXMrMingsvHljJ9dDG5Wam1M9CVuyKS9lqOn2ZD0zHebDrKhubwDuHwiS4AsjONadXFzKop5fKxJcyqKWFadTHZSfxYTAW/iMgA7k7z0VPhTwX9nw6O9p9mmpOVwYzRI86ZJppSUZQ0z0hW8IuIxMDd2Xv4ZPiTQfMx1u87ysbmY/1XJudkZlCUl0V+dib5OZnnfs/OpCAnk7ycTAoi7XmRtoH9C/rfy+pfNz8nk+xMi9t1C3oQi4hIDMyM8aMKGT+qkDtnhx+L2dcXfhDOm81H2XqwgxNnejjZ1cvp7l5OdvVyqquXoye72N/Vy6nuqPbuXoY6js7MMAqywzuP/OxMqkfk8cT91w/Dlp5LwS8iEiUjw/ovIhsKd+dMTx+nIjuBU929b7/u6u3feZzqjnrd9fZO43R3L3nZl2ZaScEvIhIHZkZednhK5533RE0syXHUQkRE4kbBLyKSZhT8IiJpRsEvIpJmFPwiImlGwS8ikmYU/CIiaUbBLyKSZhLyXj1m1gbsucDVy4FDcSwnkWjbklcqb5+2LTGMd/eKWDomZPBfDDNriPVGRclG25a8Unn7tG3JR1M9IiJpRsEvIpJmUjH4Hwu6gGGkbUteqbx92rYkk3Jz/CIi8t5SccQvIiLvIWWC38wWmtk2M2s0sweDrieezKzWzH5pZlvMbJOZfT7omuLNzDLNbK2ZPRd0LfFkZqVm9hMz2xr5+xv+xytdQmb2hci/yY1m9h9mlhd0TRfKzH5kZq1mtjGqbaSZrTCzHZHviX6r/ZikRPCbWSbwEHA7EAI+bmahYKuKqx7gT9x9BnAd8NkU2z6AzwNbgi5iGHwXeMHdpwOzSaFtNLMa4H8C9e5+GZAJ3BdsVRfln4CFA9oeBF509zrgxchy0kuJ4AfmAI3uvsvdu4AlwKKAa4obdz/g7m9EXncQDo+aYKuKHzMbC3wIeDzoWuLJzEYAc4F/AHD3Lnc/GmxVcZcF5JtZFlAA7A+4ngvm7q8Ahwc0LwL+OfL6n4G7L2lRwyRVgr8G2Be13EQKBWM0M5sAXAn8JthK4uo7wJ8BfUEXEmeTgDbgHyPTWI+bWWHQRcWLuzcDfwfsBQ4Ax9x9ebBVxV2Vux+A8AAMqAy4nrhIleC3QdpS7nQlMysC/gv4I3c/HnQ98WBmdwCt7r4m6FqGQRZwFfCIu18JnCBFpgoAIvPdi4CJwBig0Mz+W7BVSSxSJfibgNqo5bEk8UfOwZhZNuHQ/zd3fzLoeuLoRuAuM9tNeIruVjP7cbAlxU0T0OTuZz+d/YTwjiBV3Aa85e5t7t4NPAncEHBN8dZiZqMBIt9bA64nLlIl+FcDdWY20cxyCB9geibgmuLGzIzwPPEWd/8/QdcTT+7+5+4+1t0nEP57+4W7p8So0d0PAvvMbFqk6f3A5gBLire9wHVmVhD5N/p+UujgdcQzwO9FXv8e8NMAa4mbrKALiAd37zGzB4BlhM8s+JG7bwq4rHi6Efgd4E0zWxdp+1/uvjTAmiQ2nwP+LTIg2QX8QcD1xI27/8bMfgK8QfjMs7Uk8ZWuZvYfwDyg3MyagC8D3wKeMLNPEd7RfTS4CuNHV+6KiKSZVJnqERGRGCn4RUTSjIJfRCTNKPhFRNKMgl9EJM0o+EVE0oyCX0QkzSj4RUTSzP8HjQiQdIakaSQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(sorted(H1[0], reverse=True)[:12])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 基本上也就数值最大的10个词有很强的区分性"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 169,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "43.645567658711265"
      ]
     },
     "execution_count": 169,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clf.reconstruction_err_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 可以基于 **词频** 进行主题建模也可以基于 **TF_IDF** 进行主题建模"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Truncated SVD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 2min 48s, sys: 396 ms, total: 2min 49s\n",
      "Wall time: 2min 49s\n"
     ]
    }
   ],
   "source": [
    "%time u, s, v = np.linalg.svd(vectors, full_matrices=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {},
   "outputs": [],
   "source": [
    "# !pip install fbpca"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[facebook/fbpca: Fast Randomized PCA/SVD](https://github.com/facebook/fbpca)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.utils.extmath import randomized_svd\n",
    "import fbpca"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 308 ms, sys: 603 ms, total: 912 ms\n",
      "Wall time: 81.9 ms\n"
     ]
    }
   ],
   "source": [
    "# Computes a truncated randomized SVD\n",
    "%time u, s, v = randomized_svd(tfidf, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 81 ms, sys: 146 ms, total: 227 ms\n",
      "Wall time: 21.3 ms\n"
     ]
    }
   ],
   "source": [
    "%time u, s, v = fbpca.pca(tfidf, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- SVD和NMF都可以进行稀疏矩阵的分解"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}