graph_model.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import jieba\n",
    "import re\n",
    "import jieba.posseg as pseg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Demonstration of inference in HMM. \n",
    "`./data/initial_vector.txt` contains initial probability  $\\pi$   \n",
    "`./data/emit_probability.txt` contains the transition probabilities $A$  \n",
    "`./data/emit_probability.txt` contains the emission probability $p(y_{t+1}|y_t)$  \n",
    "This is a program for recognizing chinese organization, the calculation of above probability is based on `People's daily`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['始##始', '人民日报', '出版社', '头版头条', '末##末']\n",
      "['S', 'I', 'D', 'A', 'A']\n",
      "===识别的实体===\n",
      "人民日报出版社\n",
      "--------------------------------------------------------------------------------\n",
      "['始##始', '新', '中国', '成立', '啦', '末##末']\n",
      "['S', 'I', 'D', 'A', 'A']\n",
      "===识别的实体===\n",
      "新中国\n",
      "--------------------------------------------------------------------------------\n",
      "['始##始', '中国移动', '研究所', '在', '京', '召开大会', '末##末']\n",
      "['S', 'I', 'D', 'A', 'A']\n",
      "===识别的实体===\n",
      "中国移动研究所\n",
      "--------------------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "import jieba\n",
    "class OrgRecognize:\n",
    "    def __init__(self):\n",
    "        self.hidden_states = [\"A\", \"B\", \"C\", \"D\",\"F\",\"G\",\"I\",\"J\",\"K\",\"L\",\"M\",\"P\",\"S\",\"W\",\"X\",\"Z\"]\n",
    "        self.initial_vector = self.load_initial_vector()\n",
    "        self.transision_matrix = self.load_transition_matrix(hidden_states=self.hidden_states)\n",
    "        self.emission_matrix = self.load_emission_matrix(hidden_states=self.hidden_states)\n",
    "\n",
    "    def load_patterns(self):\n",
    "        \"\"\"\n",
    "        organization pattern\n",
    "        :return: list all patterns\n",
    "        \"\"\"\n",
    "        result = []\n",
    "        with open(\"./data/nt.pattern.txt\", \"r\") as file:\n",
    "            datas = file.readlines()\n",
    "            for line in datas:\n",
    "                result.append(line.strip())\n",
    "        return result\n",
    "\n",
    "    def load_transition_matrix(self,hidden_states):\n",
    "        \"\"\"\n",
    "        Load transition matrix\n",
    "        :return: dict：key first state，value dict--key is next state，value is the corresponding probability\n",
    "        \"\"\"\n",
    "        result = {x: {} for x in hidden_states}\n",
    "        with open(\"./data/transition_probability.txt\",\"r\") as file:\n",
    "            datas = file.readlines()\n",
    "            for line in datas:\n",
    "                split_line = line.strip().split(\",\")\n",
    "                result[split_line[0]][split_line[1]] =  split_line[2]\n",
    "        return result\n",
    "    \n",
    "    def load_initial_vector(self):\n",
    "        \"\"\"\n",
    "        Load initial probabilities\n",
    "        :return: dict：key is the hidden state，value is the probability\n",
    "        \"\"\"\n",
    "        result = {}\n",
    "        with open(\"./data/initial_vector.txt\",\"r\") as file:\n",
    "            datas = file.readlines()\n",
    "            for line in datas:\n",
    "                split_line = line.strip().split(\",\")\n",
    "                result[split_line[0]] =  split_line[2]\n",
    "        return result\n",
    "\n",
    "    def load_emission_matrix(self,hidden_states):\n",
    "        \"\"\"\n",
    "        Load emission matrix\n",
    "        :param hidden_states: list of hidden state\n",
    "        :return: dict：key is the hidden state，value is a dict key is observable var，value is the probability\n",
    "        \"\"\"\n",
    "        result = {x:{} for x in hidden_states}\n",
    "        with open(\"./data/emit_probability.txt\",\"r\") as file:\n",
    "            datas = file.readlines()\n",
    "            for line in datas:\n",
    "                split_line = line.strip().split(\",\")\n",
    "                result[split_line[0]][split_line[1]] = split_line[2]\n",
    "        return  result\n",
    "    def get_observed_states(self,sentence):\n",
    "        return sentence\n",
    "\n",
    "    def viterbi(self, input_sentence):\n",
    "        \"\"\"\n",
    "        Inference the best hidden state sequence\n",
    "        input_setence: observable tokens\n",
    "        :return: best state sequence\n",
    "        \"\"\"\n",
    "        hidden_states=self.hidden_states\n",
    "        initial_probability=self.initial_vector\n",
    "        transition_probability=self.transision_matrix\n",
    "        emit_probability=self.emission_matrix\n",
    "        observed_states = self.get_observed_states(sentence=input_sentence)\n",
    "        self.observed_states = observed_states\n",
    "        result = []\n",
    "        compute_recode = [] #记录每一次的计算结果\n",
    "        #初始化\n",
    "        tmp_result = {}\n",
    "        for state in hidden_states:\n",
    "            if observation[0] in emit_probability[state] :\n",
    "                tmp_result[state] = eval(initial_probability[state])*eval(emit_probability[state][observation[0]])\n",
    "            else:\n",
    "                tmp_result[state] = 0\n",
    "        compute_recode.append(tmp_result)\n",
    "\n",
    "        #对于之后的词语，继续计算\n",
    "        for index,word in enumerate(observation[1:]):\n",
    "            tmp_result = {}\n",
    "            for current_state in hidden_states:\n",
    "                #取最大值：上一次的所有状态(x)*转移到当前状态（current_state）*发射概率\n",
    "                if word in emit_probability[current_state]:\n",
    "                    tmp_result[current_state] = max([compute_recode[index][x]*eval(transition_probability[x][current_state])*\n",
    "                                                              eval(emit_probability[current_state][word]) for x in hidden_states])\n",
    "                else:\n",
    "                    tmp_result[current_state] = 0\n",
    "            compute_recode.append(tmp_result)\n",
    "\n",
    "        #返回概率最大的标签序列\n",
    "        tag_sequence = []\n",
    "        for recode in compute_recode:\n",
    "            tag_sequence.append(max(recode, key=recode.get))\n",
    "        return tag_sequence\n",
    "    def get_organization(self, observation, sequence, patterns):\n",
    "        \"\"\"\n",
    "        得到识别的机构名\n",
    "        :param observation: 单词序列\n",
    "        :param sequence: 标注序列\n",
    "        :param patterns: 模式串\n",
    "        :return: list，机构名\n",
    "        \"\"\"\n",
    "        org_indices = []  # 存放机构名的索引\n",
    "        orgs = [] # 存放机构名字符串\n",
    "        tag_sequence_str = \"\".join(tag_sequence)  # 转为字符串\n",
    "        for pattern in patterns:\n",
    "            if pattern in tag_sequence_str:\n",
    "                start_index = (tag_sequence_str.index(pattern))\n",
    "                end_index = start_index + len(pattern)\n",
    "                org_indices.append([start_index,end_index])\n",
    "        if len(org_indices)!=0:\n",
    "            for start,end in org_indices:\n",
    "                orgs.append(\"\".join(observation[start:end]))\n",
    "        return orgs\n",
    "\n",
    "    \n",
    "def ner(orgrecog, setence):\n",
    "    sentence = [\"始##始\"]\n",
    "    sentence.extend(list(jieba.cut(sentence_str)))\n",
    "    sentence.append(\"末##末\")\n",
    "    print(sentence)\n",
    "    tag_sequence = orgrecog.viterbi(sentence)\n",
    "    print( tag_sequence )\n",
    "    patterns = orgrecog.load_patterns()\n",
    "    results = orgrecog.get_organization(sentence,tag_sequence,patterns)\n",
    "    \n",
    "    if len(results) == 0:\n",
    "        print (\"未识别到机构名\")\n",
    "        print (tag_sequence)\n",
    "    else:\n",
    "        print('===识别的实体===')\n",
    "        for result in results:\n",
    "            print (result)\n",
    "    print('-'*80)\n",
    "            \n",
    "if __name__ == '__main__':\n",
    "    orgrecog = OrgRecognize()  \n",
    "    # Note the corpus is too old, so some new organization cannot be recognized\n",
    "    sentence_str = \"人民日报出版社头版头条\"\n",
    "    ner(orgrecog, sentence_str)\n",
    "    sentence_str = \"新中国成立啦\"\n",
    "    ner(orgrecog, sentence_str)\n",
    "    sentence_str = \"中国移动研究所在京召开大会\"\n",
    "    ner(orgrecog, sentence_str)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Demonstration of CRF for sequence labeling\n",
    "reference: [link](http://www.albertauyeung.com/post/python-sequence-labelling-with-crf/)  \n",
    "The task is to identify country names, the code include training and test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "blocker (N) || energy (N) || corp (N) || said (I) || an (I) || offering (I) || of (I) || 20 (I) || mln (I) || common (I) || shares (I) || is (I) || underway (I) || at (I) || 2.625 (I) || dlrs (I) || per (I) || share (I) || through (I) || underwriters (I) || led (I) || by (I) || drexel (N) || burnham (N) || lambert (N) || inc (N) || and (I) || alex. (N) || brown (N) || and (N) || sons (N) || inc (N) || absb. (I) || the (I) || company (I) || is (I) || offering (I) || 19.7 (I) || mln (I) || shares (I) || and (I) || shareholders (I) || the (I) || rest. (I) || before (I) || the (I) || offering (I) || it (I) || had (I) || about (I) || 33.6 (I) || mln (I) || shares (I) || outstanding. (I)\n",
      "              precision    recall  f1-score   support\n",
      "\n",
      "           I       0.98      0.98      0.98      2805\n",
      "           N       0.85      0.89      0.87       437\n",
      "\n",
      "   micro avg       0.96      0.96      0.96      3242\n",
      "   macro avg       0.92      0.93      0.92      3242\n",
      "weighted avg       0.96      0.96      0.96      3242\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import codecs\n",
    "import numpy as np\n",
    "import nltk\n",
    "import pycrfsuite\n",
    "from bs4 import BeautifulSoup as bs\n",
    "from bs4.element import Tag\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import classification_report\n",
    "\n",
    "# Read data file and parse the XML\n",
    "with codecs.open(\"n3-collection/reuters.xml\", \"r\", \"utf-8\") as infile:\n",
    "    soup = bs(infile, \"html5lib\")\n",
    "\n",
    "docs = []\n",
    "for elem in soup.find_all(\"document\"):\n",
    "    texts = []\n",
    "\n",
    "    # Loop through each child of the element under \"textwithnamedentities\"\n",
    "    for c in elem.find(\"textwithnamedentities\").children:\n",
    "        if type(c) == Tag:\n",
    "            if c.name == \"namedentityintext\":\n",
    "                label = \"N\"  # part of a named entity\n",
    "            else:\n",
    "                label = \"I\"  # irrelevant word\n",
    "            for w in c.text.split(\" \"):\n",
    "                if len(w) > 0:\n",
    "                    texts.append((w, label))\n",
    "    docs.append(texts)\n",
    "\n",
    "\n",
    "data = []\n",
    "for i, doc in enumerate(docs):\n",
    "\n",
    "    # Obtain the list of tokens in the document\n",
    "    tokens = [t for t, label in doc]\n",
    "\n",
    "    # Perform POS tagging\n",
    "    tagged = nltk.pos_tag(tokens)\n",
    "\n",
    "    # Take the word, POS tag, and its label\n",
    "    data.append([(w, pos, label) for (w, label), (word, pos) in zip(doc, tagged)])\n",
    "\n",
    "\n",
    "def word2features(doc, i):\n",
    "    word = doc[i][0]\n",
    "    postag = doc[i][1]\n",
    "\n",
    "    # Common features for all words\n",
    "    features = [\n",
    "        'bias',\n",
    "        'word.lower=' + word.lower(),\n",
    "        'word[-3:]=' + word[-3:],\n",
    "        'word[-2:]=' + word[-2:],\n",
    "        'word.isupper=%s' % word.isupper(),\n",
    "        'word.istitle=%s' % word.istitle(),\n",
    "        'word.isdigit=%s' % word.isdigit(),\n",
    "        'postag=' + postag\n",
    "    ]\n",
    "\n",
    "    # Features for words that are not\n",
    "    # at the beginning of a document\n",
    "    if i > 0:\n",
    "        word1 = doc[i-1][0]\n",
    "        postag1 = doc[i-1][1]\n",
    "        features.extend([\n",
    "            '-1:word.lower=' + word1.lower(),\n",
    "            '-1:word.istitle=%s' % word1.istitle(),\n",
    "            '-1:word.isupper=%s' % word1.isupper(),\n",
    "            '-1:word.isdigit=%s' % word1.isdigit(),\n",
    "            '-1:postag=' + postag1\n",
    "        ])\n",
    "    else:\n",
    "        # Indicate that it is the 'beginning of a document'\n",
    "        features.append('BOS')\n",
    "\n",
    "    # Features for words that are not \n",
    "    # at the end of a document\n",
    "    if i < len(doc)-1:\n",
    "        word1 = doc[i+1][0]\n",
    "        postag1 = doc[i+1][1]\n",
    "        features.extend([\n",
    "            '+1:word.lower=' + word1.lower(),\n",
    "            '+1:word.istitle=%s' % word1.istitle(),\n",
    "            '+1:word.isupper=%s' % word1.isupper(),\n",
    "            '+1:word.isdigit=%s' % word1.isdigit(),\n",
    "            '+1:postag=' + postag1\n",
    "        ])\n",
    "    else:\n",
    "        # Indicate that it is the 'end of a document'\n",
    "        features.append('EOS')\n",
    "\n",
    "    return features\n",
    "\n",
    "# A function for extracting features in documents\n",
    "def extract_features(doc):\n",
    "    return [word2features(doc, i) for i in range(len(doc))]\n",
    "\n",
    "# A function fo generating the list of labels for each document\n",
    "def get_labels(doc):\n",
    "    return [label for (token, postag, label) in doc]\n",
    "\n",
    "\n",
    "X = [extract_features(doc) for doc in data]\n",
    "y = [get_labels(doc) for doc in data]\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)\n",
    "\n",
    "trainer = pycrfsuite.Trainer(verbose=False)\n",
    "\n",
    "# Submit training data to the trainer\n",
    "for xseq, yseq in zip(X_train, y_train):\n",
    "    trainer.append(xseq, yseq)\n",
    "\n",
    "# Set the parameters of the model\n",
    "trainer.set_params({\n",
    "    # coefficient for L1 penalty\n",
    "    'c1': 0.1,\n",
    "\n",
    "    # coefficient for L2 penalty\n",
    "    'c2': 0.01,  \n",
    "\n",
    "    # maximum number of iterations\n",
    "    'max_iterations': 200,\n",
    "\n",
    "    # whether to include transitions that\n",
    "    # are possible, but not observed\n",
    "    'feature.possible_transitions': True\n",
    "})\n",
    "\n",
    "# Provide a file name as a parameter to the train function, such that\n",
    "# the model will be saved to the file when training is finished\n",
    "trainer.train('crf.model')\n",
    "\n",
    "# Generate predictions\n",
    "tagger = pycrfsuite.Tagger()\n",
    "tagger.open('crf.model')\n",
    "y_pred = [tagger.tag(xseq) for xseq in X_test]\n",
    "\n",
    "# Let's take a look at a random sample in the testing set\n",
    "i = 12\n",
    "outs = []\n",
    "for x, y in zip(y_pred[i], [x[1].split(\"=\")[1] for x in X_test[i]]):\n",
    "    outs.append(\"%s (%s)\" % (y, x))\n",
    "print(' || '.join(outs))\n",
    "\n",
    "# Create a mapping of labels to indices\n",
    "labels = {\"N\": 1, \"I\": 0}\n",
    "\n",
    "# Convert the sequences of tags into a 1-dimensional array\n",
    "predictions = np.array([labels[tag] for row in y_pred for tag in row])\n",
    "truths = np.array([labels[tag] for row in y_test for tag in row])\n",
    "\n",
    "# Print out the classification report\n",
    "print(classification_report(\n",
    "    truths, predictions,\n",
    "    target_names=[\"I\", \"N\"]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Demonstration of PLSA\n",
    "Here we provide a implementation of PLSA from [link](https://github.com/isnowfy/plsa.git)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 141,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 iter\n",
      "likelihood 2.123584 \n",
      "1 iter\n",
      "likelihood -12.995615 \n",
      "2 iter\n",
      "likelihood -11.750363 \n",
      "3 iter\n",
      "likelihood -10.234484 \n",
      "4 iter\n",
      "likelihood -8.967368 \n",
      "5 iter\n",
      "likelihood -8.107269 \n",
      "6 iter\n",
      "likelihood -7.478987 \n",
      "7 iter\n",
      "likelihood -7.054530 \n",
      "8 iter\n",
      "likelihood -6.857159 \n",
      "9 iter\n",
      "likelihood -6.811286 \n",
      "10 iter\n",
      "likelihood -6.807646 \n",
      "11 iter\n",
      "likelihood -6.807588 \n",
      "12 iter\n",
      "likelihood -6.807588 \n",
      "13 iter\n",
      "likelihood -6.807588 \n",
      "0 iter\n",
      "likelihood 4.392996 \n",
      "1 iter\n",
      "likelihood -14.747961 \n",
      "2 iter\n",
      "likelihood -14.542116 \n",
      "3 iter\n",
      "likelihood -14.353453 \n",
      "4 iter\n",
      "likelihood -14.086266 \n",
      "5 iter\n",
      "likelihood -13.601759 \n",
      "6 iter\n",
      "likelihood -12.697529 \n",
      "7 iter\n",
      "likelihood -11.290563 \n",
      "8 iter\n",
      "likelihood -9.747169 \n",
      "9 iter\n",
      "likelihood -8.513129 \n",
      "10 iter\n",
      "likelihood -7.593646 \n",
      "11 iter\n",
      "likelihood -7.028783 \n",
      "12 iter\n",
      "likelihood -6.834483 \n",
      "13 iter\n",
      "likelihood -6.808457 \n",
      "14 iter\n",
      "likelihood -6.807592 \n",
      "15 iter\n",
      "likelihood -6.807588 \n",
      "16 iter\n",
      "likelihood -6.807588 \n"
     ]
    }
   ],
   "source": [
    "\n",
    "import math\n",
    "import operator\n",
    "import random\n",
    "import gzip\n",
    "import sys\n",
    "import marshal\n",
    "from functools import reduce\n",
    "\n",
    "def cos_sim(p, q):\n",
    "    sum0 = sum(map(lambda x:x*x, p))\n",
    "    sum1 = sum(map(lambda x:x*x, q))\n",
    "    sum2 = sum(map(lambda x:x[0]*x[1], zip(p, q)))\n",
    "    return sum2/(sum0**0.5)/(sum1**0.5)\n",
    "\n",
    "def _rand_mat(sizex, sizey):\n",
    "    ret = []\n",
    "    for i in range(sizex):\n",
    "        ret.append([])\n",
    "        for _ in range(sizey):\n",
    "            ret[-1].append(random.random())\n",
    "        norm = sum(ret[-1])\n",
    "        for j in range(sizey):\n",
    "            ret[-1][j] /= norm\n",
    "    return ret\n",
    "\n",
    "class Plsa:\n",
    "\n",
    "    def __init__(self, corpus, topics=2):\n",
    "        self.topics = topics\n",
    "        self.corpus = corpus\n",
    "        self.docs = len(corpus)\n",
    "        self.each = list(map(sum, map(lambda x:x.values(), corpus)))\n",
    "        self.words = max(reduce(operator.add, map( lambda x:list(x.keys()), corpus)))+1\n",
    "        self.likelihood = 0\n",
    "        self.zw = _rand_mat(self.topics, self.words)\n",
    "        self.dz = _rand_mat(self.docs, self.topics)\n",
    "        self.dw_z = None\n",
    "        self.p_dw = []\n",
    "        self.beta = 0.8\n",
    "\n",
    "    def save(self, fname, iszip=True):\n",
    "        d = {}\n",
    "        for k, v in self.__dict__.items():\n",
    "            if hasattr(v, '__dict__'):\n",
    "                d[k] = v.__dict__\n",
    "            else:\n",
    "                d[k] = v\n",
    "        if sys.version_info[0] == 3:\n",
    "            fname = fname + '.3'\n",
    "        if not iszip:\n",
    "            marshal.dump(d, open(fname, 'wb'))\n",
    "        else:\n",
    "            f = gzip.open(fname, 'wb')\n",
    "            f.write(marshal.dumps(d))\n",
    "            f.close()\n",
    "\n",
    "    def load(self, fname, iszip=True):\n",
    "        if sys.version_info[0] == 3:\n",
    "            fname = fname + '.3'\n",
    "        if not iszip:\n",
    "            d = marshal.load(open(fname, 'rb'))\n",
    "        else:\n",
    "            try:\n",
    "                f = gzip.open(fname, 'rb')\n",
    "                d = marshal.loads(f.read())\n",
    "            except IOError:\n",
    "                f = open(fname, 'rb')\n",
    "                d = marshal.loads(f.read())\n",
    "            f.close()\n",
    "        for k, v in d.items():\n",
    "            if hasattr(self.__dict__[k], '__dict__'):\n",
    "                self.__dict__[k].__dict__ = v\n",
    "            else:\n",
    "                self.__dict__[k] = v\n",
    "\n",
    "    def _cal_p_dw(self):\n",
    "        self.p_dw = []\n",
    "        for d in range(self.docs):\n",
    "            self.p_dw.append({})\n",
    "            for w in self.corpus[d]:\n",
    "                tmp = 0\n",
    "                for _ in range(self.corpus[d][w]):\n",
    "                    for z in range(self.topics):\n",
    "                        tmp += (self.zw[z][w]*self.dz[d][z])**self.beta\n",
    "                self.p_dw[-1][w] = tmp\n",
    "\n",
    "    def _e_step(self):\n",
    "        self._cal_p_dw()\n",
    "        self.dw_z = []\n",
    "        for d in range(self.docs):\n",
    "            self.dw_z.append({})\n",
    "            for w in self.corpus[d]:\n",
    "                self.dw_z[-1][w] = []\n",
    "                for z in range(self.topics):\n",
    "                    self.dw_z[-1][w].append(((self.zw[z][w]*self.dz[d][z])**self.beta)/self.p_dw[d][w])\n",
    "\n",
    "    def _m_step(self):\n",
    "        for z in range(self.topics):\n",
    "            self.zw[z] = [0]*self.words\n",
    "            for d in range(self.docs):\n",
    "                for w in self.corpus[d]:\n",
    "                    self.zw[z][w] += self.corpus[d][w]*self.dw_z[d][w][z]\n",
    "            norm = sum(self.zw[z])\n",
    "            for w in range(self.words):\n",
    "                self.zw[z][w] /= norm\n",
    "        for d in range(self.docs):\n",
    "            self.dz[d] = [0]*self.topics\n",
    "            for z in range(self.topics):\n",
    "                for w in self.corpus[d]:\n",
    "                    self.dz[d][z] += self.corpus[d][w]*self.dw_z[d][w][z]\n",
    "            for z in range(self.topics):\n",
    "                self.dz[d][z] /= self.each[d]\n",
    "\n",
    "    def _cal_likelihood(self):\n",
    "        self.likelihood = 0\n",
    "        for d in range(self.docs):\n",
    "            for w in self.corpus[d]:\n",
    "                self.likelihood += self.corpus[d][w]*math.log(self.p_dw[d][w])\n",
    "\n",
    "    def train(self, max_iter=100):\n",
    "        cur = 0\n",
    "        for i in range(max_iter):\n",
    "            print ('%d iter' % i)\n",
    "            self._e_step()\n",
    "            self._m_step()\n",
    "            self._cal_likelihood()\n",
    "            print ('likelihood %f ' % self.likelihood)\n",
    "            if cur != 0 and abs((self.likelihood-cur)/cur) < 1e-8:\n",
    "                break\n",
    "            cur = self.likelihood\n",
    "\n",
    "    def inference(self, doc, max_iter=100):\n",
    "        doc = dict(filter(lambda x:x[0]<self.words, doc.items()))\n",
    "        words = sum(doc.values())\n",
    "        ret = []\n",
    "        for i in range(self.topics):\n",
    "            ret.append(random.random())\n",
    "        norm = sum(ret)\n",
    "        for i in range(self.topics):\n",
    "            ret[i] /= norm\n",
    "        tmp = 0\n",
    "        for _ in range(max_iter):\n",
    "            p_dw = {}\n",
    "            for w in doc:\n",
    "                p_dw[w] = 0\n",
    "                for _ in range(doc[w]):\n",
    "                    for z in range(self.topics):\n",
    "                        p_dw[w] += (ret[z]*self.zw[z][w])**self.beta\n",
    "            # e setp\n",
    "            dw_z = {}\n",
    "            for w in doc:\n",
    "                dw_z[w] = []\n",
    "                for z in range(self.topics):\n",
    "                    dw_z[w].append(((self.zw[z][w]*ret[z])**self.beta)/p_dw[w])\n",
    "            # m step\n",
    "            ret = [0]*self.topics\n",
    "            for z in range(self.topics):\n",
    "                for w in doc:\n",
    "                    ret[z] += doc[w]*dw_z[w][z]\n",
    "            for z in range(self.topics):\n",
    "                ret[z] /= words\n",
    "            # cal likelihood\n",
    "            likelihood = 0\n",
    "            for w in doc:\n",
    "                likelihood += doc[w]*math.log(p_dw[w])\n",
    "            if tmp != 0 and abs((likelihood-tmp)/tmp) < 1e-8:\n",
    "                break\n",
    "            tmp = likelihood\n",
    "        return ret\n",
    "\n",
    "    def post_prob_sim(self, docd, q):\n",
    "        sim = 0\n",
    "        for w in docd:\n",
    "            tmp = 0\n",
    "            for z in range(self.topics):\n",
    "                tmp += self.zw[z][w]*q[z]\n",
    "            sim += docd[w]*math.log(tmp)\n",
    "        return sim\n",
    "\n",
    "######### unittest #################################\n",
    "\n",
    "\n",
    "\n",
    "def test_train():\n",
    "    corpus = [{0:2,3:5},{0:5,2:1},{1:2,4:5}]\n",
    "    p = Plsa(corpus)\n",
    "    p.train()\n",
    "    assert cos_sim(p.dz[0], p.dz[1])>cos_sim(p.dz[0], p.dz[2])\n",
    "    assert p.post_prob_sim(p.corpus[0], p.dz[1])>p.post_prob_sim(p.corpus[0], p.dz[2])\n",
    "\n",
    "def test_inference():\n",
    "    corpus = [{0:2,3:5},{0:5,2:1},{1:2,4:5}]\n",
    "    p = Plsa(corpus)\n",
    "    p.train()\n",
    "    z = p.inference({0:4, 6:7})\n",
    "    assert abs(cos_sim(p.dz[0], p.dz[1])-cos_sim(p.dz[0], z))<1e-8\n",
    "test_train()\n",
    "test_inference()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Demonstration of LDA  for text analysis\n",
    "\n",
    "Analyze the topic in chapters of 《Journey to the West》  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_weicheng(filename='story.txt'):\n",
    "    lines = open(filename).readlines()\n",
    "    books = []\n",
    "    for line in lines:\n",
    "        line = line.strip()\n",
    "        if not line:\n",
    "            continue\n",
    "        if line.startswith('第') and line.endswith('章'):\n",
    "            books.append(line)\n",
    "        else:\n",
    "            books[-1] += line\n",
    "    for i in range(len(books)):\n",
    "        books[i] = ' '.join(jieba.cut(books[i]))\n",
    "    return books\n",
    "\n",
    "def load_honglou(filename='honglou.txt'):\n",
    "    lines = open(filename).readlines()\n",
    "    book = [[]]\n",
    "    for line in lines:\n",
    "        line = line.strip()\n",
    "        if not line:\n",
    "            continue\n",
    "        if '书香屋' in line:\n",
    "            book[-1].pop()\n",
    "            book.append([])\n",
    "        else:\n",
    "            book[-1].append(line)\n",
    "    return [' '.join(jieba.cut(''.join(chp))) for chp in book]\n",
    "\n",
    "def load_xiyou(filename='xiyouji_wuchengen.txt'):\n",
    "    \n",
    "    def _is_sep(line):\n",
    "        pat = re.compile(u'第.{1,3}回')\n",
    "        if pat.search(line):\n",
    "            return True\n",
    "        return False\n",
    "    \n",
    "    lines = open(filename).readlines()\n",
    "    book = [[]]\n",
    "    for line in lines:\n",
    "        line = line.strip()\n",
    "        if not line:\n",
    "            continue\n",
    "        if _is_sep(line):\n",
    "            book.append([])\n",
    "        else:\n",
    "            book[-1].append(line)\n",
    "    print(len(book))\n",
    "    book = [''.join(chp) for chp in book if len(chp)]\n",
    "    book_noun = []\n",
    "    for chp in book:\n",
    "        words = pseg.cut(chp)\n",
    "        book_noun.append( ' '.join([word for word, flag in words if flag[0] == 'n']))\n",
    "        # book_noun.append( ' '.join(jieba.cut(chp)))\n",
    "    return book_noun"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "122\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/users/sulixin/anaconda2/envs/py36/lib/python3.6/site-packages/sklearn/decomposition/online_lda.py:314: DeprecationWarning: n_topics has been renamed to n_components in version 0.19 and will be removed in 0.21\n",
      "  DeprecationWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "LDA\n",
      "Topic 0:三藏 师父 行者 形容 大圣 这里 一个 八戒 意思 后来 不知 唐僧 那里 不是 所以 古代 小龙 土地 故事 甚么\n",
      "Topic 1:行者 八戒 师父 三藏 一个 大圣 沙僧 妖精 那里 怎么 唐僧 我们 不知 不是 菩萨 和尚 呆子 两个 只见 三个\n",
      "Topic 2:意思 这里 佛教 道教 形容 后来 古代 一种 传说 称为 故事 比喻 神仙 所以 一个 这是 就是 地方 又称 认为\n",
      "Topic 3:行者 八戒 师父 三藏 一个 唐僧 沙僧 怎么 那里 我们 大圣 和尚 妖精 不知 两个 菩萨 甚么 长老 不是 国王\n",
      "Topic 4:菩萨 太宗 御弟 南无 女王 袈裟 玄奘 三藏 长老 唐王 太师 圣僧 法师 行者 真经 锡杖 唐僧 取经 师父 徒弟\n",
      "Topic 5:大圣 菩萨 行者 天王 悟空 玉帝 一个 猴王 如来 不知 太宗 只见 那里 如何 两个 太子 大王 陛下 哪吒 龙王\n",
      "Topic 6:行者 菩萨 一个 这里 八戒 甚么 师父 唐僧 称为 就是 不知 袈裟 乃是 后来 猴王 两个 原来 那个 一种 今日\n",
      "Topic 7:八戒 国王 公主 三藏 师父 长老 唐僧 行者 菩萨 沙僧 两个 徒弟 和尚 陛下 如何 意思 不敢 不知 驸马 一个\n",
      "Topic 8:光蕊 玄奘 丞相 婆婆 母亲 龙王 我儿 和尚 打死 一个 师父 今日 夫人 只见 父母 长老 唐王 孩儿 夜叉 报仇\n",
      "Topic 9:八戒 行者 菩萨 师父 一个 佛教 大圣 意思 这里 三藏 后来 悟空 怎么 取经 甚么 不知 不是 故事 形容 沙僧\n"
     ]
    }
   ],
   "source": [
    "from sklearn.feature_extraction.text import TfidfVectorizer, CountVectorizer\n",
    "from sklearn.datasets import fetch_20newsgroups\n",
    "from sklearn.decomposition import NMF, LatentDirichletAllocation\n",
    "\n",
    "def display_topics(model, feature_names, no_top_words):\n",
    "    for topic_idx, topic in enumerate(model.components_):\n",
    "        print (\"Topic %d:\" % (topic_idx) + \" \".join([feature_names[i]\n",
    "                        for i in topic.argsort()[:-no_top_words - 1:-1]]) )\n",
    "\n",
    "\n",
    "documents = load_xiyou()\n",
    "print()\n",
    "no_features = 1000\n",
    "\n",
    "# NMF is able to use tf-idf\n",
    "tfidf_vectorizer = TfidfVectorizer(max_df=0.95, min_df=2, max_features=no_features, stop_words='english')\n",
    "tfidf = tfidf_vectorizer.fit_transform(documents)\n",
    "tfidf_feature_names = tfidf_vectorizer.get_feature_names()\n",
    "\n",
    "# LDA can only use raw term counts for LDA because it is a probabilistic graphical model\n",
    "tf_vectorizer = CountVectorizer(max_df=0.95, min_df=2, max_features=no_features, stop_words='english')\n",
    "tf = tf_vectorizer.fit_transform(documents)\n",
    "tf_feature_names = tf_vectorizer.get_feature_names()\n",
    "\n",
    "no_topics = 10\n",
    "\n",
    "\n",
    "# Run LDA\n",
    "lda = LatentDirichletAllocation(n_topics=no_topics, max_iter=20, learning_method='online', learning_offset=50.,random_state=0).fit(tf)\n",
    "\n",
    "no_top_words = 20\n",
    "print('\\nLDA')\n",
    "display_topics(lda, tf_feature_names, no_top_words)"
   ]
  }
 ],
 "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}