Corex_w2v_v1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/muruwu/opt/anaconda3/lib/python3.8/site-packages/gensim/corpora/dictionary.py:11: DeprecationWarning: Using or importing the ABCs from 'collections' instead of from 'collections.abc' is deprecated since Python 3.3, and in 3.9 it will stop working\n",
      "  from collections import Mapping, defaultdict\n",
      "/Users/muruwu/opt/anaconda3/lib/python3.8/site-packages/scipy/sparse/sparsetools.py:21: DeprecationWarning: `scipy.sparse.sparsetools` is deprecated!\n",
      "scipy.sparse.sparsetools is a private module for scipy.sparse, and should not be used.\n",
      "  _deprecated()\n",
      "/Users/muruwu/opt/anaconda3/lib/python3.8/site-packages/gensim/models/doc2vec.py:73: DeprecationWarning: Using or importing the ABCs from 'collections' instead of from 'collections.abc' is deprecated since Python 3.3, and in 3.9 it will stop working\n",
      "  from collections import namedtuple, defaultdict, Iterable\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy.sparse as ss\n",
    "import random\n",
    "import pandas as pd\n",
    "\n",
    "from corextopic import corextopic as ct\n",
    "from corextopic import vis_topic as vt # jupyter notebooks will complain matplotlib is being loaded twice\n",
    "\n",
    "from sklearn.datasets import fetch_20newsgroups\n",
    "from sklearn.feature_extraction.text import CountVectorizer\n",
    "\n",
    "import pyLDAvis\n",
    "import pyLDAvis.gensim_models\n",
    "pyLDAvis.enable_notebook()\n",
    "\n",
    "import gensim\n",
    "from sklearn.datasets import fetch_20newsgroups\n",
    "from gensim.utils import simple_preprocess\n",
    "from gensim.parsing.preprocessing import STOPWORDS\n",
    "from gensim.corpora import Dictionary\n",
    "\n",
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pickle\n",
    "import pyLDAvis\n",
    "import pyLDAvis.gensim_models\n",
    "\n",
    "\n",
    "from nltk.stem import WordNetLemmatizer, SnowballStemmer\n",
    "from nltk.stem.porter import *\n",
    "import numpy as np\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import warnings\n",
    "warnings.simplefilter('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 10372 documents.\n",
      "95% of the documents are below: 605 words.\n",
      "Solid green line indicates median, dotted red line indicates 95 percentile. Outliers may be cropped.\n",
      "19039\n",
      "In the dataset there are 19475 textual documents\n",
      "And this is the first one:\n",
      " skip to main content search upload sign up log in books video audio software imagesabout blog projects help donate contact jobs volunteer people search metadata search text contents search\n",
      "There are 19475 documents.\n",
      "95% of the documents are below: 53 words.\n",
      "Solid green line indicates median, dotted red line indicates 95 percentile. Outliers may be cropped.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAY40lEQVR4nO3db4wc933f8fenpP7QsimR0ZE43tEhnR7UUkIs8Q4MDRWGa1rROeEfAa2As5GIKBQeK9CFjQZIyAYtyQdE1T4wErUVIVJ2RCG2CEaJqxMDySEYq0ULQdQdTZuiJJYnixVPdyYvSl0rDqJU9LcP9kd7cty72xne7exyPi9gMbO/ndn9LMHvd3d/M7uniMDMzKrhH5QdwMzMmsdN38ysQtz0zcwqxE3fzKxC3PTNzCpkYdkBZnP77bfHqlWrmvZ4I+Mj9K7onX3D8XFYsaI5j2Vmc2MO6rZdjIyM/GVEdEwdV6ufstnX1xfDw8NNezztFbG7gX+TkRHovbaG3fBjmdncmIO6bReSRiKib+q4p3fMzCrETb+ovqteQM2s1blu3fTNzKrETd/MrELc9IvavbvsBGaWl+vWTb+wPXvKTmBmeblu3fQLq8i5vmbXFdetm35hExNlJzCzvFy3bvpmZlXipp+159bGt127dv5ymNn8cN266Rc2MlJ2AjPLy3U7e9OXdIekU5nLjyV9RdJSSccknUvLJZl9dkkalXRW0v2Z8V5Jp9Ntj0nSfD2xeTc4WHYCM8vLdTt704+IsxFxd0TcDfQCfwN8C9gJHI+IHuB4uo6kNcAAcCfQDzwuaUG6u/3AINCTLv1z+3Sa6ODBshOYWV6u29zTOxuAtyLifwNbgENp/BDwQFrfAhyOiA8i4m1gFFgnqRNYHBEvR+2nPZ/O7GNmZk2Qt+kPAM+k9eURMQGQlsvSeBdwIbPPWBrrSutTx68iaVDSsKThycnJnBHNzGw6DTd9STcCm4E/nm3TOmMxw/jVgxEHIqIvIvo6Oq76GwCt4d13y05gZnm5bnO90/88cDIiLqbrF9OUDWl5KY2PASsz+3UD42m8u854e/JZAGbtx3Wbq+l/gZ9P7QAMAVvT+lbgucz4gKSbJK2mdsD2RJoCel/S+nTWzkOZfdrP5s1lJzCzvFy3jf2NXEkfAe4DtmeGHwWOSHoYeAd4ECAizkg6ArwOfAjsiIjLaZ9HgKeARcAL6WJmZk3SUNOPiL8BfmHK2HvUzuapt/0+YF+d8WHgrvwxzcxsLvgbuUU98UTZCcwsL9etm35h/mafWftx3brpF9bGvyBhVlmuWzd9M7MqcdM3M6sQN/2iNm4sO4GZ5eW6ddMv7Pnny05gZnm5bt30C9u0qewEZpaX69ZNv7CjR8tOYGZ5uW7d9M3MqsRN38ysQtz0i4q6fwrAzFqZ69ZNv7ADB8pOYGZ5uW7d9Avbvn32bcystbhu3fTNzKrETd/MrELc9IsaGio7gZnl5bp10y+st7fsBGaWl+vWTb+wrq6yE5hZXq5bN30zsyppqOlLuk3Ss5LelPSGpE9JWirpmKRzabkks/0uSaOSzkq6PzPeK+l0uu0xyX/GxsysmRp9p/8HwIsR8Y+ATwJvADuB4xHRAxxP15G0BhgA7gT6gcclLUj3sx8YBHrSpX+OnkfzbdtWdgIzy8t1O3vTl7QY+DTwNYCI+LuI+BGwBTiUNjsEPJDWtwCHI+KDiHgbGAXWSeoEFkfEyxERwNOZfdqPv9ln1n5ctw290/8EMAn8oaTvSnpS0i3A8oiYAEjLZWn7LuBCZv+xNNaV1qeOX0XSoKRhScOTk5O5nlDT+CwAs/bjum2o6S8E1gL7I+Ie4CekqZxp1JunjxnGrx6MOBARfRHR19HR0UDEEpw8WXYCM8vLddtQ0x8DxiLilXT9WWovAhfTlA1peSmz/crM/t3AeBrvrjNuZmZNMmvTj4gfAhck3ZGGNgCvA0PA1jS2FXgurQ8BA5JukrSa2gHbE2kK6H1J69NZOw9l9mk/nZ1lJzCzvFy3LGxwu38FfEPSjcAPgH9B7QXjiKSHgXeABwEi4oykI9ReGD4EdkTE5XQ/jwBPAYuAF9KlPY37Q4pZ23HdNnbKZkScSnPsvxwRD0TE/4mI9yJiQ0T0pOVfZbbfFxG/FBF3RMQLmfHhiLgr3faldBZPe9qzp+wEZpaX69bfyC1s796yE5hZXq5bN30zsypx0zczqxA3/aKGh8tOYGZ5uW7d9M3MqsRNv6i+vrITmFlerls3fTOzKnHTNzOrEDf9onbvLjuBmeXlunXTL8zf7DNrP65bN/3CVqwoO4GZ5eW6ddMvbGKi7ARmlpfr1k3fzKxK3PSLWru27ARmlpfr1k2/sJGRshOYWV6uWzf9wgYHy05gZnm5bt30Czt4sOwEZpaX69ZN38ysStz0zcwqpKGmL+m8pNOSTkkaTmNLJR2TdC4tl2S23yVpVNJZSfdnxnvT/YxKekyS5v4pNcm775adwMzyct3meqf/TyPi7oi48tukO4HjEdEDHE/XkbQGGADuBPqBxyUtSPvsBwaBnnTpv/anUBKfBWDWfly31zS9swU4lNYPAQ9kxg9HxAcR8TYwCqyT1AksjoiXIyKApzP7tJ/Nm8tOYGZ5uW4bbvoB/LmkEUlXznlaHhETAGm5LI13ARcy+46lsa60PnXczMyaZGGD290bEeOSlgHHJL05w7b15uljhvGr76D2wjII8PGPf7zBiGZmNpuG3ulHxHhaXgK+BawDLqYpG9LyUtp8DFiZ2b0bGE/j3XXG6z3egYjoi4i+jo6Oxp9NMz3xRNkJzCwv1+3sTV/SLZI+dmUd+FXgNWAI2Jo22wo8l9aHgAFJN0laTe2A7Yk0BfS+pPXprJ2HMvu0H3+zz6z9uG4bmt5ZDnwrnV25EPhmRLwo6VXgiKSHgXeABwEi4oykI8DrwIfAjoi4nO7rEeApYBHwQrq0Jwmi7uyUmbUq1+3sTT8ifgB8ss74e8CGafbZB+yrMz4M3JU/ppmZzQV/I9fMrELc9IvauLHsBGaWl+vWTb+w558vO4GZ5eW6ddMvbNOmshOYWV6uWzf9wo4eLTuBmeXlunXTNzOrEjd9M7MKcdMvquJf8DBrS65bN/3CDhwoO4GZ5eW6ddMvbPv2shOYWV6uWzd9M7MqcdM3M6sQN/2ihobKTmBmeblu3fQL6+0tO4GZ5eW6ddMvrMt/3tes7bhu3fTNzKrETd/MrELc9Ivatq3sBGaWl+vWTb8wf7PPrP24bt30C/NZAGbtx3XbeNOXtEDSdyUdTdeXSjom6VxaLslsu0vSqKSzku7PjPdKOp1ue0yS5vbpNNHJk2UnMLO8XLe53ul/GXgjc30ncDwieoDj6TqS1gADwJ1AP/C4pAVpn/3AINCTLv3XlN7MzHJpqOlL6gZ+HXgyM7wFOJTWDwEPZMYPR8QHEfE2MAqsk9QJLI6IlyMigKcz+7Sfzs6yE5hZXq7bht/p/z7wO8BPM2PLI2ICIC2XpfEu4EJmu7E01pXWp45fRdKgpGFJw5OTkw1GbLLx8bITmFlertvZm76kjcCliBhp8D7rzdPHDONXD0YciIi+iOjr6Oho8GGbbM+eshOYWV6u24be6d8LbJZ0HjgMfFbSHwEX05QNaXkpbT8GrMzs3w2Mp/HuOuPtae/eshOYWV6u29mbfkTsiojuiFhF7QDtX0TEbwBDwNa02VbgubQ+BAxIuknSamoHbE+kKaD3Ja1PZ+08lNnHzMyaYOE17PsocETSw8A7wIMAEXFG0hHgdeBDYEdEXE77PAI8BSwCXkgXMzNrklxNPyJeAl5K6+8BG6bZbh+wr874MHBX3pAtaXi47ARmlpfr1t/INTOrEjf9ovr6yk5gZnm5bt30zcyqxE3fzKxC3PSL2r277ARmlpfr1k2/MH+zz6z9uG7d9AtbsaLsBGaWl+vWTb+wiYmyE5hZXq5bN30zsypx0y9q7dqyE5hZXq5bN/3CRhr9pWkzaxmuWzf9wgYHy05gZnm5bt30Czt4sOwEZpaX69ZN38ysStz0zcwqxE2/qHffLTuBmeXlunXTL8xnAZi1H9etm35hmzeXncDM8nLduumbmVWJm76ZWYXM2vQl3SzphKTvSTojaW8aXyrpmKRzabkks88uSaOSzkq6PzPeK+l0uu0xSZqfp9UETzxRdgIzy8t129A7/Q+Az0bEJ4G7gX5J64GdwPGI6AGOp+tIWgMMAHcC/cDjkhak+9oPDAI96dI/h8+lufzNPrP247qdvelHzV+nqzekSwBbgENp/BDwQFrfAhyOiA8i4m1gFFgnqRNYHBEvR0QAT2f2aT9t/CHFrLJct43N6UtaIOkUcAk4FhGvAMsjYgIgLZelzbuAC5ndx9JYV1qfOl7v8QYlDUsanpyczPN8zMxsBg01/Yi4HBF3A93U3rXfNcPm9V5KY4bxeo93ICL6IqKvo6OjkYhmZtaAXGfvRMSPgJeozcVfTFM2pOWltNkYsDKzWzcwnsa764y3p40by05gZnm5bhs6e6dD0m1pfRHwOeBNYAjYmjbbCjyX1oeAAUk3SVpN7YDtiTQF9L6k9emsnYcy+7Sf558vO4GZ5eW6beidfifwHUnfB16lNqd/FHgUuE/SOeC+dJ2IOAMcAV4HXgR2RMTldF+PAE9SO7j7FvDCHD6X5tq0qewEZpaX65aFs20QEd8H7qkz/h6wYZp99gH76owPAzMdD2gfR4+WncDM8nLd+hu5ZmZV4qZvZlYhbvpFRd2zTc2slblu3fQLO3Cg7ARmlpfr1k2/sO3by05gZnm5bt30zcyqxE3fzKxC3PSLGhoqO4GZ5eW6ddMvrLe37ARmlpfr1k2/sK66vwptZq3Mdeumb2ZWJW76ZmYV4qZf1LZtZScws7xct276hfmbfWbtx3Xrpl+YzwIwaz+uWzf9wk6eLDuBmeXlunXTNzOrEjf9ojo7y05gZnm5bt30CxsfLzuBmeXlunXTL2zPnrITmFlertvZm76klZK+I+kNSWckfTmNL5V0TNK5tFyS2WeXpFFJZyXdnxnvlXQ63faYJM3P02qCvXvLTmBmebluG3qn/yHw2xHxj4H1wA5Ja4CdwPGI6AGOp+uk2waAO4F+4HFJC9J97QcGgZ506Z/D52JmZrNYONsGETEBTKT19yW9AXQBW4DPpM0OAS8Bv5vGD0fEB8DbkkaBdZLOA4sj4mUASU8DDwAvzOHzKWbPrWUnMDNrilxz+pJWAfcArwDL0wvClReGZWmzLuBCZrexNNaV1qeO13ucQUnDkoYnJyfzRGye4eGyE5hZXq7bxpu+pI8CfwJ8JSJ+PNOmdcZihvGrByMORERfRPR1dHQ0GtHMzGbRUNOXdAO1hv+NiPjTNHxRUme6vRO4lMbHgJWZ3buB8TTeXWe8PfX1lZ3AzPJy3TZ09o6ArwFvRMRXMzcNAVvT+lbgucz4gKSbJK2mdsD2RJoCel/S+nSfD2X2MTOzJpj1QC5wL/CbwGlJp9LYvwEeBY5Iehh4B3gQICLOSDoCvE7tzJ8dEXE57fcI8BSwiNoB3PIP4pqZVUgjZ+/8D+rPxwNsmGaffcC+OuPDwF15Aras3bvLTmBmeblu/Y3cwvzNPrP247p10y9sxYqyE5hZXq5bN/3CJibKTmBmeblu3fTNzKrETb+otWvLTmBmeblu3fQLGxkpO4GZ5eW6ddMvbHCw7ARmlpfr1k2/sIMHy05gZnm5bt30zcyqxE3fzKxC3PSLevfdshOYWV6uWzf9wnwWgFn7cd266Re2eXPZCcwsL9etm76ZWZW46ZuZVYibfj17bq1dZvLEE83JYmZzx3Xrpl+Yv9ln1n5ct276hWm6PyZmZi3Ldeumb2ZWJW76M1i188+untufba7fzKyFzdr0JX1d0iVJr2XGlko6JulcWi7J3LZL0qiks5Luz4z3SjqdbntMav3PWedv/uL0N27c2LwgZjY3XLcNvdN/CuifMrYTOB4RPcDxdB1Ja4AB4M60z+OSFqR99gODQE+6TL3PlrVq55/97PKz63f+y5JTmVluzz9fdoLSzdr0I+K/A381ZXgLcCitHwIeyIwfjogPIuJtYBRYJ6kTWBwRL0dEAE9n9mlLTz67t+wIZpbXpk1lJyhd0Tn95RExAZCWy9J4F3Ahs91YGutK61PH65I0KGlY0vDk5GTBiHPn/M1fvGqq53NvvVpSGjMr7OjRshOUbq4P5Nabp48ZxuuKiAMR0RcRfR0dHXMWzsys6hYW3O+ipM6ImEhTN5fS+BiwMrNdNzCexrvrjLe1K3P8U51/9NebnMTMrDFFm/4QsBV4NC2fy4x/U9JXgRXUDtieiIjLkt6XtB54BXgI+E/XlPwaTG3W529ubL/sFM+q3/XHRLO2E9NOMFRGI6dsPgO8DNwhaUzSw9Sa/X2SzgH3petExBngCPA68CKwIyIup7t6BHiS2sHdt4AX5vi5NNUXTr1YdgQzy+vAgbITlG7Wd/oR8YVpbtowzfb7gH11xoeBu3Kla2H//tv/mWfubpuzTs0MYPv2yv/+jr+Ra2ZWIUXn9G0GPsBrZq3K7/QLevif/duyI5hZXkNDZSconZt+QaeX/8OyI5hZXr29ZSconad3Cjrx+Nbcp21eNe2z6Odjnvoxa4Kursqftul3+gXM+OubZmYtzE3fzKxCPL1T1Nob5vTufMaPWRNs21Z2gtK56Re1aRH87fw/jF8MzOaQv5Hr6Z3CDvx12QnMLC+fveOmX9jET8tOYGZ5nTxZdoLSeXqnTU037TMdTweZGfidfnEflU/dNGs3nZ1lJyidm35Rv/2xshOYWV7jbf+3m66Zm35RLzXh1B0zm1t79pSdoHSe0y/qv/0dfKbBP7nVAnzqpxmwd2/lG3/lm37V5+V9QNisWirf9C0ff2Iwa2/Vavp7boU9/3du7mvbLXNzP9cJf2KwtjA8XHaC0lWr6VvL8CcGs3I0velL6gf+AFgAPBkRjzblgffc+rPl+Zth1d9+89ru7+BPYPfia89lf0/eTwwz8QuIXaWvr/K/p9/Upi9pAfBfgPuAMeBVSUMR8fq8PeiVZj/F+Zu/eO2Nfw7vx+beXL6AzIXpXoT8qceaqdnv9NcBoxHxAwBJh4EtwPw0/Wka/hVVP3PHmivvi1CrvWhdD86T/0V2ro5XtcpxL0UTP+pI+udAf0T8Vrr+m8CvRMSXpmw3CAymq3cAZws+5O3AXxbct9mcdX60U1Zor7zOOj/mKusvRkTH1MFmv9NXnbGrXnUi4gBwzT98LWk4Ivqu9X6awVnnRztlhfbK66zzY76zNvtnGMaAlZnr3YB/DMPMrEma3fRfBXokrZZ0IzAADDU5g5lZZTV1eiciPpT0JeDb1E7Z/HpEnJnHh2ynv43mrPOjnbJCe+V11vkxr1mbeiDXzMzK5Z9WNjOrEDd9M7MKuS6bvqR+SWcljUraWXYeAElfl3RJ0muZsaWSjkk6l5ZLMrftSvnPSrq/iTlXSvqOpDcknZH05VbNmh77ZkknJH0v5d3b4nkXSPqupKOtnDM9/nlJpyWdkjTcynkl3SbpWUlvpv+7n2rFrJLuSP+eVy4/lvSVpmaNiOvqQu0A8VvAJ4Abge8Ba1og16eBtcBrmbH/COxM6zuB/5DW16TcNwGr0/NZ0KScncDatP4x4H+lPC2XNT2+gI+m9RuAV4D1LZz3XwPfBI626v+BTNbzwO1TxloyL3AI+K20fiNwW6tmzWReAPwQ+MVmZm3qk2zSP+SngG9nru8CdpWdK2VZxd9v+meBzrTeCZytl5na2U6fKinzc9R+K6kdsn4EOAn8Sivmpfa9lOPAZzNNv+VyZh6zXtNvubzAYuBt0okprZx1Sr5fBf5ns7Nej9M7XcCFzPWxNNaKlkfEBEBaLkvjLfEcJK0C7qH27rlls6Ypk1PAJeBYRLRq3t8Hfgf4aWasFXNeEcCfSxpJP40CrZn3E8Ak8Idp6uxJSbe0aNasAeCZtN60rNdj02/opx5aXOnPQdJHgT8BvhIRP55p0zpjTc0aEZcj4m5q76TXSbprhs1LyStpI3ApIkYa3aXOWLP/H98bEWuBzwM7JH16hm3LzLuQ2tTp/oi4B/gJtSmS6ZT+b5u+nLoZ+OPZNq0zdk1Zr8em304/9XBRUidAWl5K46U+B0k3UGv434iIP23lrFkR8SPgJaCf1st7L7BZ0nngMPBZSX/Ugjl/JiLG0/IS8C1qv5LbinnHgLH0CQ/gWWovAq2Y9YrPAycj4mK63rSs12PTb6efehgCtqb1rdTmz6+MD0i6SdJqoAc40YxAkgR8DXgjIr7ayllT3g5Jt6X1RcDngDdbLW9E7IqI7ohYRe3/5F9ExG+0Ws4rJN0i6WNX1qnNP7/Winkj4ofABUl3pKEN1H6uveWyZnyBn0/tXMnUnKzNPnjRpAMkv0btrJO3gN8rO0/K9AwwAfw/aq/eDwO/QO3A3rm0XJrZ/vdS/rPA55uY859Q+/j4feBUuvxaK2ZNj/3LwHdT3teAf5fGWzJvevzP8PMDuS2Zk9o8+ffS5cyVOmrhvHcDw+n/wX8FlrRw1o8A7wG3ZsaaltU/w2BmViHX4/SOmZlNw03fzKxC3PTNzCrETd/MrELc9M3MKsRN38ysQtz0zcwq5P8D+YcMrJ6HhWwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#just for convenience of outputting all files to this specific project folder for this model\n",
    "os.chdir('Corex_w2v/Amanda_model')\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "def describe_training_documents(list_of_docs):\n",
    "    print('There are',len(list_of_docs),'documents.')\n",
    "    document_lengths = list(map(lambda x: len(x.split()),list_of_docs))\n",
    "    percentile_50 = int(np.percentile(document_lengths,50))\n",
    "    percentile_95 = int(np.percentile(document_lengths,95))\n",
    "    print('95% of the documents are below:',percentile_95,'words.')\n",
    "    plt.axvline(percentile_50, lw=1, color='g')\n",
    "    plt.axvline(percentile_95, lw=1, color='r', linestyle='--')\n",
    "    _ = plt.hist(document_lengths, bins=50, range=(0,percentile_95+100))\n",
    "    print('Solid green line indicates median, dotted red line indicates 95 percentile. Outliers may be cropped.')\n",
    "\n",
    "def flatten_list(l):\n",
    "    return [item for sublist in l for item in sublist]\n",
    "\n",
    "def tokenize(text):\n",
    "    return [token for token in simple_preprocess(text) if token not in STOPWORDS]\n",
    "\n",
    "def chop(text):\n",
    "    temp=text.split()\n",
    "    num=len(temp)//30\n",
    "    new=[]\n",
    "    for i in range(num):\n",
    "        if i !=num-1:\n",
    "            a=temp[30*i:30*(i+1)]\n",
    "        else:\n",
    "            a=temp[30*(num-1):]\n",
    "        new.append(\" \".join(a))\n",
    "    return new\n",
    "\n",
    "# Get 20 newsgroups data\n",
    "newsgroups = fetch_20newsgroups(subset='train', remove=('headers', 'footers', 'quotes'))\n",
    "documents_train = list(np.load('../../training/train.npy')) # historical materials 4451 documents\n",
    "documents_train2 = list(np.load('../../training/train2.npy'))   # census bureau 4226 documents\n",
    "df_occsc = pd.read_csv('../../OCC_pairs.csv').rename(columns={'OCC_DES':'Full Occupation'})\n",
    "assert(df_occsc['Full Occupation'].nunique() == len(df_occsc))\n",
    "occ_list = list(set(list(df_occsc['Full Occupation'])))\n",
    "\n",
    "occ_list2=[]\n",
    "#eliminate confusing words in the occ_list\n",
    "confusing_words=['occupation','occupations','Occupations','men']\n",
    "for occ in occ_list:\n",
    "    a=occ.split(' ')\n",
    "    for word in confusing_words:\n",
    "        if word in a:\n",
    "            a.remove(word)\n",
    "    occ_list2.append(' '.join(a))\n",
    "occ_list=occ_list2\n",
    "\n",
    "with open(\"../../training/nyt_text_modified.txt\",'r') as f:\n",
    "    nyt_text2=[]\n",
    "    lines=f.readlines()\n",
    "    for line in lines:\n",
    "        nyt_text2.append(line)\n",
    "\n",
    "# Q: why select only 4000 \n",
    "\n",
    "twenty_news=[]\n",
    "for news in newsgroups.data:\n",
    "    if len(news.split())>15:\n",
    "        twenty_news.append(news)\n",
    "describe_training_documents(twenty_news)\n",
    "random.seed(10)\n",
    "twenty_news=random.sample(twenty_news,4000)\n",
    "\n",
    "# but fit the model with the training dataset\n",
    "# corex model requires them to be the same shape \n",
    "documents = []\n",
    "documents.extend(documents_train)\n",
    "documents.extend(documents_train2)\n",
    "documents.extend(nyt_text2)\n",
    "\n",
    "documents_mod=[]\n",
    "for doc in documents:\n",
    "    chopped=chop(doc)\n",
    "    for a in chopped:\n",
    "        tokens=simple_preprocess(a)\n",
    "        new_tokens=[]\n",
    "        for word in tokens:\n",
    "            if not word.isdigit():\n",
    "                new_tokens.append(word)\n",
    "        new=\" \".join(new_tokens)\n",
    "        documents_mod.append(new)\n",
    "\n",
    "document_total=documents_mod[:]\n",
    "document_total.extend(occ_list)\n",
    "print(len(documents_mod))\n",
    "print(\"In the dataset there are\", len(document_total), \"textual documents\")\n",
    "print(\"And this is the first one:\\n\", documents_mod[0])\n",
    "\n",
    "describe_training_documents(document_total)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "below are merely test functions\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Q: but this isn't successfully parsed for some reason"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'second thoughts on voluntary code making by american industry do not reveal very clear cut situation the wage hour price structure under the nra codes is it is presumed that wages and working hours can be embodied in voluntary codes but where is the line of demarcation when it comes to voluntary stabilization of prices when the anti trust laws are considered'"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len_arr=np.asarray([len(item.split(' ')) for item in document_total])\n",
    "document_total[np.argmax(len_arr)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Q: And there're documents with non-English chars \n",
    "we can clean with basic regrex maybe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "穏raphic division the proportion was higher for the foreignborn white mafos than for the males of either of the nath積 white classes and in ach of the geographic dhisions except\n"
     ]
    }
   ],
   "source": [
    "print(max(document_total))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'doc_length')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEHCAYAAAC0pdErAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd5wU9f348df7CnDA0Y9ejiZIUAROLEhRkYBEURMTNTG2r8REjSamEDWWWEKSX2JMYqLEhsaCJVZAQaKADTialEPqUQ846h3t6uf3x87ebZmtt7szd/d+Ph732N3Z2Zn3zc7Oez5lPiPGGJRSSqlAaU4HoJRSyp00QSillLKlCUIppZQtTRBKKaVsaYJQSillK8PpAKLVoUMHk5ub63QYSilVryxbtmy/MSYnns/WmwSRm5tLfn6+02EopVS9IiLb4v2sVjEppZSypQlCKaWULU0QSimlbGmCUEopZUsThFJKKVuaIJRSStnSBKGUUsqWJgilVEocL6/k8837nQ5DxUAThFIqJX795mqu+fdidhw87nQoKkqaIJRSKbFhTykAx8orHY5ERUsThFIKgNKTFTS2O0yeKK+ioqo6IcuqqjYcLYs9+RljKDlZkZAYEk0ThFKK/UfLOO2Bufz9f5ucDiWlTr3vA77/9OKELOvB99Yy+P4PKa+MLeG8+OU2Tn9gLtsOHEtIHImU9AQhIm1E5A0RWS8iBSJyjoi0E5F5IrLRemyb7DiUUqHtKykDYPbqIocjSb0lWw8mZDlvLtsJQHmMJZJ56/YCUHjAfW0zqShBPA58YIwZCAwBCoCpwHxjTH9gvvVaKaWUiyQ1QYhIK2A08AyAMabcGHMYmAzMsGabAVyWzDjiMfHxRfzrk80JX+5r+TsY8chHMdX1bthbSu7UWfzkpWUR512z6wh9755N0ZETdQnTlabNWc9V079wOgyVZHtLTtLv7tms3nkk5s8ePl7OKffOIXfqLJ75dGsSorN35EQFx8qrUra+VEl2CaIPUAw8JyIrRORpEWkBdDLGFAFYjx3tPiwiU0QkX0Tyi4uLkxyqv4KiEv7wwfqEL3fqm1+xr7SM6hjaAl9bugOA2av3RJz3P19uo6ra8MnXqd1eqfDkgs18uSUx1QEq9QzR7fQLNhRTWW2Y8UVhzOtYueNwTRvAQ++vi/nz8VqzK/ZkVh8kO0FkAMOAfxljhgLHiKE6yRgz3RiTZ4zJy8mJ64ZISimXEcTpEFSUkp0gdgI7jTHebgJv4EkYe0WkC4D1uC/JcYRUVW14dHYB+0pP2r7/xrKdLNhQezb+/GdbufQfn3LwWHnEZc9bt5d3V+22fW/34RP84YP1NVVNm4uPMvHxRcxeXcTRskoefG8tJyuiK7L+65PNrNtdAsCrVmnDzpbio0x+4jPeWrEzquUGKq+s5nfvreOFLwqZX7A35HxFR07w+zkFVFvFpP+t38vExxdx28vLKSgqsf3MyYoqHnxvbVTdBI+XV/KjF/N54N21NeuI5J2Vu8LGnAgVVZ7tcyiKfSOc97/azYdrI5cW7by8eDtfbD4QNP3JBbX7CEBlVTW/fXsNP31lRUxdM/eVnOTR2QVURdjuuw6f4OevreS+d9bwwLtrORZH9883lu3keIzXTARG5f19zVu3l0H3fcCJGKqByio9++SR4xVB//eLX25jaWHk0mxVteH3swvYW3KSpxZsDippvJa/g083uffq8qTectQYs0dEdojIAGPM18CFwDrr7zpgmvX4TjLjCOfTTfuZvnALW4qP8vR1Zwa9/4vXVwFQOG0SAA+85ym2PvT+Oh773hlhl33zC55bpF46pGvQe7e9vJxVO48w6bQuDO7Wmqunf8m+0jJ+8tJyfjSmD899Vkj3ts256bzeSIQTrj98sJ4/frierb+fFHa+Hzy9mN1HTvKzmYe5fGj38Au1MWv1bp79rLZe17tNAv1s5kq+3HKQ8YM6MbxXO2583rMdCopK+GDNHjY9enHQZ/7z5Tae+6yQphnpTJ04MGwc0xdu4cO1noP9FcO6cXr3NhFjv+PVlWFjToTZq4t49rOtlJys4P9dOSTu5dz28gogvljvfmu17WenzVnPtDnra6Yv2FDMi1967kTZo10Wk04L3kft/PKNr1iwoZixp+Rwbr8OYf6H5azYfrjmdatm8R1qnlm0ldsv7B/XZ8Fz4tWvY3bNb3HGF4XcMqZvVJ99a/kunvuskMoqw7aDx1m4oZixA3I4t28Hfvv2GiDyd7R46wGeWriFgj2lLLRONH0/86s3vorjv0qdVNyT+nbgJRFpAmwBbsBTcnlNRG4CtgNXpiAOW9XWGUZFVWwXCFXG0ohgoyygr7Rv17gqK5Zoz44BomnzjrX7XaBoPx6uH3io7eY9M6uO4h/xPXt103Vd3rginV27gW+MsewW3ovKIv2HgRefVcS5Teq6KQP3j1i+G9/fS4V3nw7x8VDncN71VyboYrxUS3qCMMasBPJs3row2ev2tefISXYcOo4x0Kt9c7KapLNi+2HWF3ku/9+4t7Rm3livJp23bi+DurZi24FjnNu3A59t2k/fnJa288ayw3t30PkFwTVwR45X8I+PN3LbBfZnV9PmrKdZZhpHT1Yyond70gT2Hw1d9XG8vJIvNh/ggoEd+XDtHrq2yeKdlbs5f0BHzutvf6ZYVW2Yu3YPEwZ3RmyKOTOX7mDWV9FVlWzcdxSAA0fLufft1Xx7WHeym2WyfPshurbOChlDQVEJQ3p4ShBrdx+hWWZ6yG0PsLTwID3bNWddUQl5vdqS3SzT7/2SkxU8s2gr3z2zB93aZEWM+9CxctbuLgkZ3+7DJ5i5dAedWzdjcNfWtMrKoFf7FkDk7Rfo4LFyCopKGNmvA8YYXlmyg+5tsxh9irPtc6t2HKZt8yb0bN885DxzVheRZv2PizYWk7/tINeM6ImIsOvwCfYcOcnwXp7LoRYHdEQ4WVHF/IJ9ZKQLYwfk8PLi7aSneVoyJg/tRivf7zDg92Wsz0fy2ab9tM7KpLyqmmE9PXFUWidq6Wm1383uIydZtq02vpKTFSzfdsivGho8x5AP1uyhedPQh9jFW4KrAt0mFSUIV7joLwsotepB2zbPZHC31izaWFv3t/tIbRvEB2uir/9dueNwTfEVYOMjE4OuzFy4oTimH/HnVh3ynz78mgsGdmTL/uArLK986nM27D0aslfPkRMV/GzmqqjXefd/V/P2yt38ZuJAfj+ntvfWM59uDVmMfvbTrTwyu4DHrzqDyWd0C3r/tfzo2zresC4yenO55/E/X273e7/gdxNsPzf1v6u5akRPACb97VMgfLH/yidru8mOO7UTT1/nf+5y+8srWLChmMfnb4yqiuf655awaucR1j80wbY0M/ZPnwSV3LzLnfF5Ib97fx1/vnII3x4eucrvB08vZl1RCZsemcis1UU11UkbHp5IkwznBkWY/MRnQPjt7nsR2MOzCgBo36IpEwZ35rw//A9jaj/v3Qe8Hnp/HS8t9uwPg7u1Ys2u2raUeQX7eOHGESHXaww8Yq0vHN/frDeOymrP95bhkyC8Vc5ed7yygo8DegxuP3CcwgPH+MlLyxnZr33IdX5v+pcR43Jaoxlqo9SnkezQ8Qo2WWesdvbH0Mh45IT/GCp2VSSHjsfWaHnYZ/7Sk/aNdBv2euJfv8e+0TdWW60f8LYYRtosspJqcWmZ3/RkVLBUJaEuaev+4H1gc3Ho/cKO93vwrbrwLQuEq9bbU+LZfvuPloWcx39dnlKuwX+bR9t9NJWi6al04Jjnfwj31YrAdp990jc5gKfjha/AbWEwfp+Phbc6ND1dQm7jzcXBJ28nKqpqvtOiI/adX+qLRpEgoulxBOHrzktPVvj1gNh16DgFRSV+B/NQApOIZ3meA39FVTX7j5b5LXt3mJ3qaFllnXvJADW9Q4pLy9hz5GTN/3EgyoMV+B/YjpZVsrfkJDsPHbf9f30dPFbO0bLKmmGfYx27JtqeKCfKqygNMwja/qPlMVcnFpeW1Sy3uLSspgPBsbJKdh32XJxYcrKSnYfCH5QOHSuvqcLw3T+PHPePNzD5ekVTJeWNy7cnUFmlZ9sdPh56u5RVVkUVx4nyKgoDSrcHjpZFXc9/orwq4iB10fR+Ki4tY9m2gxwvrwxKNoGvj5dXUlxaRkVVddjfUXW1YYt18M9Mi+0wWXKyouazJdZvodKnjbOiqjrsCaqbNPgqphPlVQx7aF5U8978Qj4zQhRXT3tgrt/r5dsPM/HxRUHzzV0b3JXyvnfW8sNzcv2meQ8ml//z86hi8xp8/4cxzR/KiEfm8/Blg7lz5kq/6R/axP/2il1cNrRbUJdD3y68scTl+33k3zuOX74euSpskU8d79NRXiE75k8fs6+0LGTVx5ETFby6dAdXW1VUkcxeXcRPXlpu+96IR+fXPP+oYC8fFewNWS1WUVXNUJ9t8NTCLdx6QT9aNctkyO9q97OFG4r54bNLePqHeYwb1CmqGAN94/4P8akh4fpnl/LKlLP51Zv+vWeWbz8EwPo9pVz79BKWFB6s2W7vrtrNT19Zweu3nFMzvzGege587Th4nFF//Jgpo/tEFdvDswpqqptCeWrhFkaFaN/xOvORj2qePxNQZRiYIJ74eDNPfBx5hIR/frKpptozLS226zZueG5pzXNvu98Sny6xl/7jM9vu3vPW7WGMw+1JgRp8CeJElNcSAEENTfHw/tDc7mhZJYujHKTMO9+xssQPJXDgaHlQHa6dVXEMu7AvxNm3r88DrhkId2KeXxjbdxuqcdRueGm7M+Wvdnq6ia7YEbxe3zAjFYJ8T+i/sGkYFcHvgLUkoH//kq2ez6wvKgm7fXYe8pz0zF27J2LX7FSJt/rNd5+MMT9EFOpaIN9uwW7R4BNEKMnqHummbpdKRauh7LaRqpgStdzGotEmCG8DoZ26nDA8/3mh7fSXF2+3nR7Jd58KPzid7/UbuVNnxbTsV5ZEF9O7K3eRO3VWUsamila4M9LcqbN4z6e6K3fqLP4y9+ua1+HaA95btZvcqbN4etEW2+XOScLw14Pus6+O8+0+CfD/5m7we+17jArcHicrqqL+/i//52dB0+z2z7/N3+hXffPbd9by2aboumZ+FUeJL3fqLNsxjXx7G0byfz49CuP1+zkFLNsWXWkx3gbw+qLBJwi33CHL2yWxPoo0SmW0DabJdPsrK/xe/83nxjcLN0Q+wHjrwgN73zw+f2PcMcW65wV27Q1FCD6JCdWYbSfaqoy/zNsQ03IT4e0VuxK6PGPCn1zYeWpB8MlCY9XgE4RyNzd20XSCDmCXHLp/1U2D78UUq1tfXs6sr5JzV61bQ/SAqe/qMqzy/82IrkqgLvfmmDYn8oVSAE8v2hJUZbB+TynHyyv53Xvrwg6EaOfZGO5H8NePNoTsJPHEx5tplpFe03203z1zGHeq/wj5r+X7x3bXa9FfJOm9o1kodoXwm2YsDZo2c6mnBFSXO6MlepSS3YdPJGTo+1SkGZdUdvjRBBEgWckBYFYjvJ1jJN6eL8lUEuJiw0Chuly+8MW2mJMDwD8+jv7+zpGW/+d5/u0RHwUMvxJ4L+nAq5HDiadPfuBYYgBvr7QfuTgWq3cltifPL15P0GB4Ljx4p4JWMSkVgRvP7Hy5PT4n1YeBE92swScI3T1UXUVztbyTYh0epCGJVAJNVBvE1gPBQ2o0Bg0+QShVV08tdHevlm/9/VOnQ0iYRJeGErE8QVLSm8uNJ7OaIJRSYbnxwBWtRMSeqp5QbumS70sThFLKNRJ+SY37jrn1SoNPEC5MykrVK8fjuJ90vBJdlVPXuyg2dg0+QSil6uZ4lMOrJ0K0Q/OnUmM+ydQEoZRyDTcM2+IUN/7vmiCUUq4R6WZTKrU0QSilwmrENSwppb2YHKCDdSlVNy48bqVUY/73G3yCUErVzUcF4QfzUw1X0gfrE5FCoBSoAiqNMXki0g6YCeQChcB3jTH1416dSinVSKSqBHG+MeYMY4z3juJTgfnGmP7AfOu1Ukq5TyOuY3OqimkyMMN6PgO4LGlrarzfrVJK1UkqEoQB5orIMhGZYk3rZIwpArAeO9p9UESmiEi+iOQXF9f9ph9KKeVWbiyopOKGQSONMbtFpCMwT0TWR/tBY8x0YDpAXl5eXJvPhdtcKaWCuLHHZdJLEMaY3dbjPuAtYASwV0S6AFiP+0IvoW5K9MIbpVQdpOqw7cYSRFIThIi0EJFs73NgPLAGeBe4zprtOuCdZMVw0WMLk7VopZRq0JJdxdQJeMsaYyQDeNkY84GILAVeE5GbgO3AlUmOQymlVIySmiCMMVuAITbTDwAXJnPdSiml6kavpFZKKWVLE4RSSilbmiCUUsoFXNiJSROEUkqFU1BU6nQIjtEEoZRSYTTm0Ww1QSillAvoDYOUUkrVG5oglFJK2dIEoZRSypYmCKWUcgH3tUBoglBKKRWCJgillHIDFxYhNEEopZSypQlCKaVcYMv+Y06HEEQThFJKKVuaIJRSStnSBKGUUsqWJgillFK2NEEopZSypQlCKaWULU0QSimlbGmCUEopZUsThFJKKVspSRAiki4iK0Tkfet1OxGZJyIbrce2qYhDKaVU9FJVgrgDKPB5PRWYb4zpD8y3XiullHKRpCcIEekOTAKe9pk8GZhhPZ8BXJbsOJRSSsUmFSWIvwK/Aqp9pnUyxhQBWI8dUxCHUkqpGCQ1QYjIt4B9xphlcX5+iojki0h+cXFxgqNTSikVTrJLECOBS0WkEHgVuEBE/gPsFZEuANbjPrsPG2OmG2PyjDF5OTk5SQ5VKaWUr6QmCGPMb4wx3Y0xucBVwP+MMT8A3gWus2a7DngnmXEopZSKXdQJQkSusLqlHhGREhEpFZGSONc7DbhIRDYCF1mvlVJKuUhGDPP+EbjEGFMQcU4bxphPgE+s5weAC+NZjlJKqdSIpYppb7zJQSmlVP0TsQQhIldYT/NFZCbwNlDmfd8Y898kxaaUUspB0VQxXeLz/Dgw3ue1ATRBKKVUAxQxQRhjbgAQkZHGmM983xORkckKTCmllLNiaYP4e5TTlFJKNQDRtEGcA5wL5IjIz33eagWkJyswpZRSzoqmDaIJ0NKaN9tnegnwnWQEpZRSynnRtEEsABaIyPPGmG0piEkppZQLxHKh3D9ExARMOwLkA08ZY04mLiyllGpcspvGcjhOjVgaqbcAR4F/W38lwF7gFOu1UkqpeInTAQSLJWUNNcaM9nn9nogsNMaMFpG1iQ5MKaWUs2IpQeSISE/vC+t5B+tleUKjUkop5bhYShB3AZ+KyGY8haHewE9EpAW1tw9VSikVj8AWXheIOkEYY2aLSH9gIJ4Esd6nYfqvyQhOKaWUc2JtNh8O5FqfO11EMMa8kPColFJKOS7qBCEiLwJ9gZVAlTXZAJoglFKqAYqlBJEHDDLGuLCmTCml6jc3Hlhj6cW0BuicrECUUkq5SywliA7AOhFZgv8Ngy5NeFRKKaUcF0uCeCBZQSilVGPnxtr7WLq5LhCRXkB/Y8xHItIcHe5bKaUSwn3pIYY2CBG5GXgDeMqa1A3P/amVUko1QLE0Ut8KjMQzSB/GmI1Ax2QEpZRSynmxJIgyY0zNmEsikoE7S0VKKVXvuLAJIqYEsUBE7gayROQi4HXgvXAfEJFmIrJERFaJyFoRedCa3k5E5onIRuuxbfz/glJK1X/GhefbsSSIqUAxsBr4ETAbuDfCZ8qAC4wxQ4AzgAkicra1rPnGmP7AfOu1UkopF4mlF1M1tTcLivYzBs9NhgAyrT8DTAbGWtNnAJ8Av452uUop1dC4sYopYoIQkdWEaWswxpwe4fPpwDKgH/CEMWaxiHQyxhRZny8SEdvGbhGZAkwB6Nmzp90sSinVINx+QT+nQwgSTQniW3VZgTGmCjhDRNoAb4nI4Bg+Ox2YDpCXl+fC/KqUUonROivT6RCCREwQxpht0SxIRL4wxpwTZjmHReQTYAKwV0S6WKWHLsC+aANWSimVGrE0UkfSLHCCiORYJQdEJAsYB6wH3gWus2a7DngngXEopVS9k5YmTocQJNYbBoVjVwXUBZhhtUOkAa8ZY94XkS+A10TkJmA7cGUC41BKqXonTRp2gghijPkKGGoz/QBwYTLXrZRSqm4SWcXkvvSnlFIqbrEM1tdbRJr5vM4SkVyfWa5NYFwJc8eF/Z0OQSmlInLjdRCxlCBeB6p9XldZ0wAwxqxJVFCJ1KKpjkiulPJo3kSPB7GIJUFk+A7WZz1vkviQEqtNc9eHqJRKkZ7tmjsdQr0SS4IoFpGa24uKyGRgf+JDSqwrh3eveZ6sKxXbt3BfErp30qmc27e902Eo5YhubbJsp0+/No+rR/Tkvm8Nsp1nSPfWIT8L0LW1f2/+P30n7EASManvg/XdAtwtIjtEZAeesZOmJCesxBGfrmM3juxd87xfx5YJW8f/7hqbsGXF68dj+/q9vmJYd16++Wx+NKZP0Lz3XzKI0afk1LwunDaJwmmTkh6jUgCnd2+d9HX0aGd/kM9ulsHvrziNG8/rzUv/d1bN9CeuGQZA1zZZfDb1AtvPzrljFFeN8B/y58q8HgmK2J1iGaxvM3C2iLQExBhTmryw6hc3ZP7ALmTe+9vaNXwJUFlVHfyGUg1ENA2+sf5qXXiZQtLF0ouptYj8Bc/Iqx+LyJ9FJPmnAgkwdeJAXrxphN+0p64dXvP8zNz4b0dx86jeQWOovH3ryLiXZ+fJHwwL+/4fv+1fzL1gYEfaWdVedjdCT0sT/vl9zzJ/+61BCYrS39UjGvaZVUMwdkBO5JnqIFwPwqvO9Owf15+bm5B1Pfa9IX6vbx4VXHKG6A/yf7866PItAG48rzcds5syqEsr7r/E89t58gfDeeiywdx98cCg+ft1bMldF50CwG3n21dxe5fjRrFUMT0LlALftf5KgOeSEVSi3TKmL6P6+/8Y+ua0ZGDnbAAevDTq8QP9FE6bxD2TBvlVY43q34EzerSpef23EDtaLCYM7hL2/e+e2cNvx3/2+jNrYqoOUYJo07wJhdMmcdN5vYNnCBBr9VOfDi34/RWJq5tViVc4bRLP3zAi8ox18DPrwBjIGJj27dMpnDaJwd0Sc455+dDufq9bNotcOWJ38uR1yZCujOrfIWh6y6YZLLlnHLPvGMUNVpX1hMGdufbsXlw53P+kaMI3OvPRz8dw+4X9KZw2iV98c0DNe0OsarZe7Zsz6bQuVjwRQ065WK6k7muM+bbP6wdFZGWiA1KxiXRGZLvTNcaysnKlarszmASIdQ+PpppYEnktsO9v0MU/x1hKECdE5DzvCxEZCZxIfEjJE81ZRbwmn9EV8G8Ih7pVX0Xj3knhi6eXDfXElZPdtGbamP7JrVpIlDvH9ddc1sCdY/W0G9g5m17t7bug3jwqfCn3p1bvxGvP7lUzLfBw773+wfcg37l10PiiCd3frj2nV9C0zHTPCn42zlP99tML3H0hbyxHzFuAF3zaHQ5ROyJrvZCZnhayuqRPTgu2FB+Le9mPXzWUx68Krk7q3Cp4JwRompFGWWXdG4q9VUShzm5O796GwmmT+O5TX1BcWsarU86mZ4gfYiy82/GxeRt4fP7GhCzvx/9Zxpw1e/yWf+e4U1i/p4QJf10U13LfuOUcvvPkF3WOL1pXDOvGf5fvCpo+uFsr1uwqSVkcdVU4bRK5U2clfT092jX3+0161/ns9Xnc+Hw+HbObcs+kQfx70VbbGH09dNlgXvzS/u4EGTYjpTZvEtsJY7TJo3VWJqvuH2/73sZHLq557o1/X+nJmOJIpWjuKPdzn5cvAC2s58fwDN/9VRLiavDcWN8Yj0SecTWEbZLQagjVqLhx94+miinb+svDU4poBbQGfgS4t/k9ChcM9NzptEPL2C50i/YCtIsGdQr5nl2dpzeeeAyxGsbTQ4wpP8a67qFr69AXAdkJXFxWpv9QBXYHdbv/I/BzsajLQTfUj87bMJhoeWGqFJtmJHJszPDGh9n3UqFJeu3/enafdn7vjTkl9irOEb3bRZ7J0s7mwtWJVkePppn230GfDp7ros7uU/vb9nb08HZmaRvlqAyxnjS5+aQimjvKPQggInOBYd7rH0TkAXzGYqqP7ho/gOvPzaVjq2ZRf0U52U157oYzOXqyMux8+feOo1WzTL8eTr7s2uamXzucA8fKufH5pazd7V8dkX/vOPIe/ijk+i4a1In3bz+P/p3sLwD88Zi+fGd4dzqFqPICWHz3hZz16HxEoH/HlmzYe5Snrs0LiqMyQsPiby4+NWjaryYM4MH31oX9nNdfvjsk8kx4ukjuP1rG+18V+U2/ZUxfZnxeyImKqpClkse+dwazVhcFTZ/7s9HsP1rGwWPlDOnehsz0NMY/toASm+/7x2P78uOxfUkTYfD9HwKetqjf/Hd10LyCcNv5/fjzvA1R/W9ei351PqP++HHQ9Dl3jKJTq2Zs2FvKVdO/DHr/H9cMo+RkBdnNMnhkVgEvfLGN/h1bsnHfUdv1PH/DmVz/3FK/adnNMiiNsJ+D56r9h2cV+E3zVrEs/+1FNG+SzsmKKkpPViICXaI4SfEeNL1f3ws3jmDgbz+I+LlV940nI11YveuI3/SHLx/MXd88hWYhTlQGdW3Fl7+5kE6tmga998tvDuC0bq392vEai1gq4XoC5T6vy4HchEaTYulpQkfrgBlt8a5rmyyaZqTTtGX4M+IOLcPvTNU2R66M9DQ6tWpmO6BYpOUBYbsMpqVJ2OQA0KqZ53oOAevgc5QmAWe9LZr67zJ2+c+uFGNXBxxKtKWN3PbNKTlRETS9aUYap3VrzZLCgyG7Mgb+X16ndMrmlE7ZUa0/Mz2tZptlZaZzoqIq7PzebZXdNIPSssgHXvDU0dvp1jaLVs0y/c54fTXJSKvZZ07t0gqAYT3bhkwQdiML+H5jrbMyOWKzrcFzcA2UYTXGes/mm2WmxzYuWsDuEurAHqh1c/v7Omemp9ExO/z+H9ho7Q0hzec40djEUuZ9EVgiIg+IyP3AYmBGcsJyr/QElQYbQn07kNCzqq7WGDjR3ry9bYsmtj1RoPbAGpjQEtQb16UAABSLSURBVMl3V+jb0dM0F+6uYN7vvH2MVZp20sOsJ7CKxTtnuK6cTTOCD8B9cmqTxoDO0SVNL6fvjtYykd97qn6rLjwoxDLUxiMiMgcYZU26wRizIjlhuVeoOv54XXt2r6CeF9795N5Jp9apXSJeseymV5/Zk+xmmfx85koqqw0f/2Js3Ov99cQB5OW25dx+wRcoea347UV8vvkA1cbwrdO7UDHY0K1tFl3bZPHf5bt4b9VuAB667BtcNKiTX6nqwztHs2nfUbq3ja0dJtCQHm1YteOw33Z64cazWL3riN+Zbrc2Wew67OkJ7nu8nHR6FwZ2bsXry3aycENxzfRpV5zGVJ/qqVAXMT5xzTC/xPfhnaPZdfg4zZtkUG1MTX26VzTH6pzspjU9h7yeu/5MVu48zNGTlUwc3JkXvtjGpuKjvLx4e9hl9enQIuG/k1gN7taarq2bsfvIybBXdL9968iQJySpynFu7sodU5o1xiwHlicplnoh0WdGQ3u2CU4Q1uOQHm38zuKSLZ5/LS1NuHRIV+56zXPNZNc28RfFm2akc7FN47E3rv4dW9K2RRMmnV47T5MM4ftnefqbr9h+uGb+5k0ymDC4s99yBnTOjvlM2M75A3JYteOw37R2LZoENb52b1ubIPz+H4RLhnTl7RX+3WFPCYgtr5d9g/eYgCEyov2/Ip2gXjDQv2G7bYsmnD+g9gTlRith+SaIdi2aBJ1RfDNguztl/Dc68/znhbQKUyL1HfUgFDeMteaU1HWraCBC1V3Hy+5MyzvNyROLjBhj8B58nKxa8Fa7JDqGzHT/77x2PeE/5xtGmghpaf6fC/zuo9/WsR2w0tPSamJIhsBoErEW7zJiabsKWkYdA8mwtlu0vYy864s3ZjemoeRV0NZT38vrwcz8HTWv/3GN5+K3x+ZtYHPxsaBhteM1deJAig6fYOLgLtyB/4glf/3eGTzz6VaG9UzuVdjh/OHbp/OvBZuj7tLrbXQP+dNIQeL4v1G92X+0LKrxpQDuvnggj85eH3G+V6eczburdvP3/20CPGfS+0rLQg4IF2hwt1b8/eqhtG/ZhKIjJ5gyxrMPPXz5YLq2yeI7w7vzUcFehnRvw82jegddFPafm85i1c7DLNpYzBk92sZcv37JkC58tfMwPxt3ClfmdWelT+ln5pSzg3rMRfLcDWeybncJf/rwayC4ZPKNrnUfX6ljdjNuPb+v3xhL799+Hp98vY//Nze2nmDxevTywXRvm2U7JpOd1lmZ3H5BPy4d0jWm9bi4hkkTRKCbR/fxSxDfOt3zZf/rk81AbU+fumqTlcktY+yTTdc2WUkbZTVaHVs14/5LvhH1/N5jRKhuvanQomkGD10W/cCLU0b3jSpB9O+UzV3jB9QkiFjXc/fFp9Y0mj982Wk10ztmN+OBSz3b2NtWcs+kQew4eIIP1u6pme+8/h04r38Hbg0xGmgkTTPS+d1kT7x5LdqRl1t7TcFZfdpzVoieUKGcP6Ajp3VrXZsgAs59E9H8kJ4m/PKb/qOjDu7WmsHdWqcsQXRsVfv9RENEuGv8gMgz1iNaxRQlb9f/RB3/3Fic9IqnM4X3M24+G3KMm79sl3K4jVtZkpogRKSHiHwsIgUislZE7rCmtxOReSKy0Xp0ri4lSj+/6BQy0oTc9i0izxwF34Nwx+ymPODiMeGj8ejlp9E6KzOuBNonJ/w27dG2OU3S0+I+O/vB2T2DruaNx4/H9mVwt+A+/4G+0bUVPx7b19VXyIbyw3N6cVYUVy37/menBVx/41tCiZeTJVGnuLCXa9KrmCqBu4wxy0UkG1gmIvOA64H5xphpIjIVmIrnFqbOifDlXDSoE5sevTj8TDGtrnaFS+4Zl7DlOuWas3pyzVk9Q74f7ud+98TgK699ZTVJZ8MjE+OMzL9apy5+PWEgv54QfFOYQLN+6ukJfrXNFc6xcOJ44a2KioXvBXCJunVtYsf4cuGR14ebk2FSSxDGmCKrayzWEB0FQDdgMrUX2c0ALktmHLFI1Xfl5n02GdvAxb8B19FtlZgeV/WxBOc2KWuDEJFcYCieK7A7GWOKwJNEANurwURkiojki0h+cXGx3Sz1lovzg61ubeK/uMxuYDzfgQwb6gHxwlM9u3X3tnUfXt1tvENlX3J6cgY91DYId0hJLyYRaQm8CdxpjCmJtkhljJkOTAfIy8urb8fU8NxchLAx/64xVMV496+vHhhPVZUhu1mGX88wgH9+fxjX/PtLlhYeSmSYrnLTeb35zvDusY1BVE9kNUln1f3jEzukhQ+nh+pwghurwpKeIEQkE09yeMkY819r8l4R6WKMKRKRLsC+ZMfhNu7bFWrZ7afRDpbmK1yX4Mz0NLIT1GXYrUSkQSYHr2jHzIpHY8oPbv5Xk5ogxFNUeAYoMMb8xeetd/HcjW6a9fhOMuNwo5woRmdNtVT/KN14xqTcIZ4TkkDeYV/q20isuVNncdWZPWiWmc7znxcmrOE/HskuQYwErgVWi4j3cuG78SSG10TkJmA7cGWS44hL4P2lEylwnKCGLlyDYWM6W1ThfXDnKBZt2B/V8PaR3DiyN73at2Dcqakf8LKuXl26I/JMKZDUBGGM+ZTQJagLk7nuRBjQOXkD5bm5a5tSThnYuRUDO0e+1iQaaWkS9q6ObuPG8rReSW1pahVpfQ/biewml8rbTdYX3nsQaLK0pzVwjYN391+0cT+5U2cFvT/+sQVsCnGjp2TTsZgs068dzpvLd9K7g89VvQk6bt1z8amMGZDD+McWJmaBDcQjlw+md04LRveP/R7FDZnmy8bpf+vt++ps2HuU8srqFEfjoQnC0qNdc+4cd0pSln3z6OhG/XRaMi8ssjvotW/ZNKork5VSztB6D6WUcpCbr/jWEkQY7RpwH3al6qubR/Vm56HgO/U1ZJuKjzKoa2Ia72OhCSKMjHT3Zvb6RrekSpR7JtXvkY/jUVXtTBuEVjEppZTLOVUNpQkijD4dkncdRGOjPXOid+FAT9/9ZF6Ho9yjIorSgVO/H61iCqNn++SMwvn7KxJzf4JE04O4O3x7eHcmnta5ZsRU1bBVxzgIZippCcIBOpSxikSTg/Ll1MWkmiAc0KmeDR6mlGqc9DTFxpM/GE5OduK7uK687yJeWrydsQPq3+BhDdUnvxhLmgij//Sx06GoxiqKwoFTlQ6aIGwka6TVNs2bcOv5/ZKybLdz68VAub5DqyjlUk61D2oVk6rhzkO4Ukq7uSqllHIVTRBKKeVyWsWkGjatv1IqbtpIrZSqt8YP6sR5/Ts4HUa9FE37gl5JrRyXzItxtADRsE3/YZ7TIdRb0f3stJFaKaWUDW2DUEop5SqaIFRKODWWjFINgVO/Hk0QSinlcjpYn1JKKVsNsgQhIs+KyD4RWeMzrZ2IzBORjdZj22TGoKKnlUBKuVNDbaR+HpgQMG0qMN8Y0x+Yb71WSqlGKZpjf4NMEMaYhcDBgMmTgRnW8xnAZcmMQbmDlk6UsufmDhxOtEF0MsYUAViPIW+OICJTRCRfRPKLi4tTFqBSSrmJNlLbMMZMN8bkGWPycnJynA5HKaUckdaIEsReEekCYD3ucyAGZcPFJV2lGrXGNFjfu8B1wDTr8R0HYlApNm5QJ4b0aAPAjSNznQ1GqXqmQQ7WJyKvAGOBDiKyE7gfT2J4TURuArYDVyYzBuUOrbMyeefWkU6HoVS95NQd5ZKaIIwxV4d468JkrlcppRqSBtnNVSmlVHhRXQeR9CjsaYJQSim30xKEcpqbL9hRqqFy889OE4RSSjmoSYZ7D8PujUwppRqBrMx0p0MISROEUkq5nFPdXDVBKKWUsqUJQimllC1NEEoppWxpglBKKWVLE4RSSilbmiCUUspB0VygqmMxKaWUchVNEEop5XL9OrZ0ZL1O3DBIKaWUj8JpkzhyooIhD84lu2kGqx/8ptMhAVqCUEopdzDWo4sG79MEoZRSLtA003M4Hjugo8OR1NIqJqWUcoFmmeks+tX5dGzV1OlQamiCUEopl+jRrrnTIfjRKiallFK2NEEopZSypQlCKaWULU0QSimlbDnWSC0iE4DHgXTgaWPMNKdiqS9+PLYvRYdPMLhba3p3aMHnmw8kfB0/Gt2H8d/olPDlut1/bjqLpYUHnQ5DKVcRY0zkuRK9UpF0YANwEbATWApcbYxZF+ozeXl5Jj8/P0URKqVUwyAiy4wxefF81qkqphHAJmPMFmNMOfAqMNmhWJRSStlwKkF0A3b4vN5pTfMjIlNEJF9E8ouLi1MWnFJKKecShN1oI0F1XcaY6caYPGNMXk5OTgrCUkop5eVUgtgJ9PB53R3Y7VAsSimlbDiVIJYC/UWkt4g0Aa4C3nUoFqWUUjYc6eZqjKkUkduAD/F0c33WGLPWiViUUkrZc+w6CGPMbGC2U+tXSikVnl5JrZRSypYjF8rFQ0SKgW1xfrwDsD+B4SSSxhYfjS0+bo4N3B1ffY2tlzEmrm6g9SZB1IWI5Md7JWGyaWzx0dji4+bYwN3xNcbYtIpJKaWULU0QSimlbDWWBDHd6QDC0Njio7HFx82xgbvja3SxNYo2CKWUUrFrLCUIpZRSMdIEoZRSylaDTxAiMkFEvhaRTSIyNQXr6yEiH4tIgYisFZE7rOkPiMguEVlp/V3s85nfWPF9LSLf9Jk+XERWW+/9TUTsRsGNNb5Ca5krRSTfmtZOROaJyEbrsW2qYxORAT7bZqWIlIjInU5uNxF5VkT2icgan2kJ21Yi0lREZlrTF4tIbh1j+5OIrBeRr0TkLRFpY03PFZETPtvwSQdiS9j3mITYZvrEVSgiK1O93ST0ccPZ/c0Y02D/8IzztBnoAzQBVgGDkrzOLsAw63k2njvnDQIeAH5hM/8gK66mQG8r3nTrvSXAOXiGR58DTExAfIVAh4BpfwSmWs+nAn9wIraA720P0MvJ7QaMBoYBa5KxrYCfAE9az68CZtYxtvFAhvX8Dz6x5frOF7CcVMWWsO8x0bEFvP9n4L5UbzdCHzcc3d8aegki5XeuM8YUGWOWW89LgQJsbobkYzLwqjGmzBizFdgEjBCRLkArY8wXxvONvgBclqSwJwMzrOczfNbjVGwXApuNMeGunE96bMaYhUDgjaoTua18l/UGcGG0pR272Iwxc40xldbLL/EMox9SKmMLw/Ht5mUt47vAK+GWkYzYwhw3HN3fGnqCiOrOdcliFeGGAoutSbdZxf9nfYqKoWLsZj0PnF5XBpgrIstEZIo1rZMxpgg8OyrQ0aHYvK7C/0fqhu3mlchtVfMZ68B+BGifoDhvxHP26NVbRFaIyAIRGeWz/lTGlqjvMVnbbRSw1xiz0WdayrdbwHHD0f2toSeIqO5cl5QVi7QE3gTuNMaUAP8C+gJnAEV4irLhYkxW7CONMcOAicCtIjI6zLypjg3x3B/kUuB1a5Jbtlsk8cSTlFhF5B6gEnjJmlQE9DTGDAV+DrwsIq1SHFsiv8dkfcdX439ikvLtZnPcCDlriPUkNLaGniAcuXOdiGTi+ZJfMsb8F8AYs9cYU2WMqQb+jaf6K1yMO/GvIkhI7MaY3dbjPuAtK469VtHUW3ze50RslonAcmPMXitOV2w3H4ncVjWfEZEMoDXRV83YEpHrgG8B37eqGLCqIQ5Yz5fhqa8+JZWxJfh7TMZ2ywCuAGb6xJzS7WZ33MDh/a2hJ4iU37nOqtN7BigwxvzFZ3oXn9kuB7y9KN4FrrJ6GPQG+gNLrOJkqYicbS3zh8A7dYythYhke5/jadRcY8VwnTXbdT7rSVlsPvzO4tyw3QIkclv5Lus7wP+8B/V4iMgE4NfApcaY4z7Tc0Qk3Xrex4ptS4pjS+T3mNDYLOOA9caYmuqZVG63UMcNnN7fIrVi1/c/4GI8PQI2A/ekYH3n4Sm2fQWstP4uBl4EVlvT3wW6+HzmHiu+r/HpcQPk4fkhbQb+gXXlex1i64On58MqYK13e+Cph5wPbLQe26U6NmuZzYEDQGufaY5tNzyJqgiowHP2dVMitxXQDE9V2iY8PU/61DG2TXjqmL37nbfHyret73sVsBy4xIHYEvY9Jjo2a/rzwC0B86ZsuxH6uOHo/qZDbSillLLV0KuYlFJKxUkThFJKKVuaIJRSStnSBKGUUsqWJgillFK2NEEopZSypQlCNVjiGWL6F07HEY54hpfu4HQcStnRBKFUPWUNl6BU0miCUA2KiNwjnhuofAQMsKadISJfSu2NdNpa0/uJyEciskpElotI3xDLHCsin4jIG+K5Ic9L3mGSfUsAIpInIp9Yzx8QkRkiMtea5woR+aN4buTygTXujtcvRWSJ9dfP+nyOiLwpIkutv5E+y50uInPxDOWsVNJoglANhogMxzPe1lA8A6+dab31AvBrY8zpeIZ7uN+a/hLwhDFmCHAuniEYQhkK3InnRi19gJFRhNQXmIRnHP7/AB8bY04DTljTvUqMMSPwDIvwV2va48Bjxpgz8Qz58LTP/MOBycaYa6KIQam4aRFVNSSjgLeMNVCdiLwLtADaGGMWWPPMAF63Bi3sZox5C8AYczLCspcYayA38dySMhf4NMJn5hhjKkRkNZ675H1gTV9tfd7rFZ/Hx6zn44BBUns/l1begRaBd40xJyKsW6k60wShGppoBxeL9T7VZT7Pq6j97VRSWxJvZvcZY0y1iFSY2oHPqvH/7Rmb52nAOYGJwEoYx2KMXam4aBWTakgWApeLSJZ1tn0JnoPpIam9G9i1wALjuRnLThG5DGpu6N48jnUW4qnyAU9VUDy+5/P4hfV8LnCbdwYROSPOZSsVN00QqsEwnnv6zsQzVPKbwCLrreuAP4nIV3juaPY7a/q1wE+t6Z8DneNY7YPA4yKyCE/JIh5NRWQxcAfwM2vaT4E8q2F9HXBLnMtWKm463LdSSilbWoJQSillSxuplbKIyGl47nzmq8wYc5YT8SjlNK1iUkopZUurmJRSStnSBKGUUsqWJgillFK2NEEopZSy9f8B3cbRkFZ7EjcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(0)\n",
    "plt.plot(len_arr)\n",
    "plt.xlabel('doc_number')\n",
    "plt.ylabel('doc_length')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "test codes end\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Transform data into a sparse matrix¶"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# count vectorizer \n",
    "(filtering out: useless phrases like 'aa','ab' with token_pattern, & min_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "vectorizer = CountVectorizer(analyzer='word',       \n",
    "                             min_df=2,\n",
    "# minimum reqd occurences of a word \n",
    "                             stop_words='english',             \n",
    "# remove stop words\n",
    "                             lowercase=True,                   \n",
    "# convert all words to lowercase\n",
    "                             token_pattern='[a-zA-Z]{3,}',  \n",
    "# num chars > 3\n",
    "                             max_features=20000,             \n",
    "# max number of uniq words    \n",
    "                            )\n",
    "data_vectorized = vectorizer.fit_transform(document_total)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "doc_word = vectorizer.fit_transform(document_total)\n",
    "doc_word = ss.csr_matrix(doc_word)\n",
    "fixed_vocabulary=vectorizer.vocabulary_\n",
    "words = list(np.asarray(vectorizer.get_feature_names()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Q: try find frequency of fixed voc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING: Some words never appear (or always appear)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<corextopic.corextopic.Corex at 0x7f8852c31550>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "topic_num=60\n",
    "topic_model = ct.Corex(n_hidden=topic_num, words=words, eps=1e-2,max_iter=100, verbose=False, seed=1)\n",
    "topic_model.fit(doc_word[:len(documents_mod)], words=words)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('american', 0.03570231998520268, 1.0),\n",
       " ('war', 0.02542871871758147, 1.0),\n",
       " ('world', 0.018345561285183373, 1.0),\n",
       " ('country', 0.012495536049215904, 1.0),\n",
       " ('society', 0.01144207092559819, 1.0),\n",
       " ('says', 0.010938772561319636, 1.0),\n",
       " ('movement', 0.010286919155662554, 1.0),\n",
       " ('action', 0.008177295178333346, 1.0),\n",
       " ('street', 0.00756682912906009, 1.0),\n",
       " ('federation', 0.007255233062166837, 1.0)]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Print a single topic from CorEx topic model\n",
    "topic_model.get_topics(topic=20, n_words=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I aded '~' here since the original function also prints out words whose absence is significant to a topic (they have s<0), so here, I added ~ before them to distinguish these absent words from the other topic words\n",
    "\n",
    "For details, you can check here: https://github.com/gregversteeg/corex_topic#running-the-corex-topic-model\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Print all topics from the CorEx topic model\n",
    "\n",
    "topics = topic_model.get_topics()\n",
    "for n,topic in enumerate(topics):\n",
    "    # w: word, mi: mutual information, s: sign\n",
    "    topic_words,_,_ = zip(*[(w,mi,s) if s > 0 else ('~'+w,mi,s) for w,mi,s in topic])\n",
    "    topic_words,_,_ = zip(*topic)\n",
    "    print('{}: '.format(n) + ','.join(topic_words))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(436, 60)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results=topic_model.predict(doc_word[len(documents_mod):])\n",
    "results.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Total Corelation and Model Selection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.365385775407672"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "topic_model.tc "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "tc_lst=[[i,v] for i,v in enumerate(topic_model.tcs)]\n",
    "tc_df=pd.DataFrame(tc_lst,columns=['topic','tc'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "<div id=\"altair-viz-b349fd614c6145558eefd68515d8f6f7\"></div>\n",
       "<script type=\"text/javascript\">\n",
       "  (function(spec, embedOpt){\n",
       "    let outputDiv = document.currentScript.previousElementSibling;\n",
       "    if (outputDiv.id !== \"altair-viz-b349fd614c6145558eefd68515d8f6f7\") {\n",
       "      outputDiv = document.getElementById(\"altair-viz-b349fd614c6145558eefd68515d8f6f7\");\n",
       "    }\n",
       "    const paths = {\n",
       "      \"vega\": \"https://cdn.jsdelivr.net/npm//vega@5?noext\",\n",
       "      \"vega-lib\": \"https://cdn.jsdelivr.net/npm//vega-lib?noext\",\n",
       "      \"vega-lite\": \"https://cdn.jsdelivr.net/npm//vega-lite@4.8.1?noext\",\n",
       "      \"vega-embed\": \"https://cdn.jsdelivr.net/npm//vega-embed@6?noext\",\n",
       "    };\n",
       "\n",
       "    function loadScript(lib) {\n",
       "      return new Promise(function(resolve, reject) {\n",
       "        var s = document.createElement('script');\n",
       "        s.src = paths[lib];\n",
       "        s.async = true;\n",
       "        s.onload = () => resolve(paths[lib]);\n",
       "        s.onerror = () => reject(`Error loading script: ${paths[lib]}`);\n",
       "        document.getElementsByTagName(\"head\")[0].appendChild(s);\n",
       "      });\n",
       "    }\n",
       "\n",
       "    function showError(err) {\n",
       "      outputDiv.innerHTML = `<div class=\"error\" style=\"color:red;\">${err}</div>`;\n",
       "      throw err;\n",
       "    }\n",
       "\n",
       "    function displayChart(vegaEmbed) {\n",
       "      vegaEmbed(outputDiv, spec, embedOpt)\n",
       "        .catch(err => showError(`Javascript Error: ${err.message}<br>This usually means there's a typo in your chart specification. See the javascript console for the full traceback.`));\n",
       "    }\n",
       "\n",
       "    if(typeof define === \"function\" && define.amd) {\n",
       "      requirejs.config({paths});\n",
       "      require([\"vega-embed\"], displayChart, err => showError(`Error loading script: ${err.message}`));\n",
       "    } else if (typeof vegaEmbed === \"function\") {\n",
       "      displayChart(vegaEmbed);\n",
       "    } else {\n",
       "      loadScript(\"vega\")\n",
       "        .then(() => loadScript(\"vega-lite\"))\n",
       "        .then(() => loadScript(\"vega-embed\"))\n",
       "        .catch(showError)\n",
       "        .then(() => displayChart(vegaEmbed));\n",
       "    }\n",
       "  })({\"config\": {\"view\": {\"continuousWidth\": 400, \"continuousHeight\": 300}}, \"data\": {\"name\": \"data-3ce031d01c3fb7d01b896ed45ee6d5c1\"}, \"mark\": \"bar\", \"encoding\": {\"tooltip\": [{\"type\": \"quantitative\", \"field\": \"topic\"}, {\"type\": \"quantitative\", \"field\": \"tc\"}], \"x\": {\"type\": \"quantitative\", \"field\": \"topic\"}, \"y\": {\"type\": \"quantitative\", \"field\": \"tc\"}}, \"width\": 800, \"$schema\": \"https://vega.github.io/schema/vega-lite/v4.8.1.json\", \"datasets\": {\"data-3ce031d01c3fb7d01b896ed45ee6d5c1\": [{\"topic\": 0, \"tc\": 0.6643436353991138}, {\"topic\": 1, \"tc\": 0.5597549334731112}, {\"topic\": 2, \"tc\": 0.46578686047479273}, {\"topic\": 3, \"tc\": 0.3211614267037636}, {\"topic\": 4, \"tc\": 0.3005400429960456}, {\"topic\": 5, \"tc\": 0.27632966493786293}, {\"topic\": 6, \"tc\": 0.2529251881364536}, {\"topic\": 7, \"tc\": 0.2439695242694419}, {\"topic\": 8, \"tc\": 0.17464486878524782}, {\"topic\": 9, \"tc\": 0.17271963614301}, {\"topic\": 10, \"tc\": 0.14952968930776953}, {\"topic\": 11, \"tc\": 0.13937257058001468}, {\"topic\": 12, \"tc\": 0.13701273125283084}, {\"topic\": 13, \"tc\": 0.12396385099195875}, {\"topic\": 14, \"tc\": 0.0808804970735721}, {\"topic\": 15, \"tc\": 0.08039454497782571}, {\"topic\": 16, \"tc\": 0.07430695787940803}, {\"topic\": 17, \"tc\": 0.07212931929812084}, {\"topic\": 18, \"tc\": 0.0698651660453549}, {\"topic\": 19, \"tc\": 0.06328750994323272}, {\"topic\": 20, \"tc\": 0.05909269558492548}, {\"topic\": 21, \"tc\": 0.056953508415820184}, {\"topic\": 22, \"tc\": 0.05181966921937532}, {\"topic\": 23, \"tc\": 0.05154655425956186}, {\"topic\": 24, \"tc\": 0.05010662938744367}, {\"topic\": 25, \"tc\": 0.0490120014417236}, {\"topic\": 26, \"tc\": 0.04184237778255194}, {\"topic\": 27, \"tc\": 0.041232705949794383}, {\"topic\": 28, \"tc\": 0.04106111134722826}, {\"topic\": 29, \"tc\": 0.03852416202969196}, {\"topic\": 30, \"tc\": 0.036258147849562564}, {\"topic\": 31, \"tc\": 0.03323562348736094}, {\"topic\": 32, \"tc\": 0.03219578734952477}, {\"topic\": 33, \"tc\": 0.02880976073600325}, {\"topic\": 34, \"tc\": 0.028701209569136073}, {\"topic\": 35, \"tc\": 0.02359333449678278}, {\"topic\": 36, \"tc\": 0.023480448841801893}, {\"topic\": 37, \"tc\": 0.016988122711807344}, {\"topic\": 38, \"tc\": 0.016784623105894937}, {\"topic\": 39, \"tc\": 0.016486945245159045}, {\"topic\": 40, \"tc\": 0.014271148827089531}, {\"topic\": 41, \"tc\": 0.01267054207373514}, {\"topic\": 42, \"tc\": 0.012555205462381643}, {\"topic\": 43, \"tc\": 0.01164167897001826}, {\"topic\": 44, \"tc\": 0.011140055556811227}, {\"topic\": 45, \"tc\": 0.011121512876298171}, {\"topic\": 46, \"tc\": 0.01057246435925399}, {\"topic\": 47, \"tc\": 0.010528404753580162}, {\"topic\": 48, \"tc\": 0.010469222859488656}, {\"topic\": 49, \"tc\": 0.010311291008189522}, {\"topic\": 50, \"tc\": 0.010045147426976778}, {\"topic\": 51, \"tc\": 0.009870102956236406}, {\"topic\": 52, \"tc\": 0.009735354262542269}, {\"topic\": 53, \"tc\": 0.009567982445959019}, {\"topic\": 54, \"tc\": 0.009233123926694727}, {\"topic\": 55, \"tc\": 0.009162551805967381}, {\"topic\": 56, \"tc\": 0.008265533877775797}, {\"topic\": 57, \"tc\": 0.008039798013457928}, {\"topic\": 58, \"tc\": 0.007966975830765029}, {\"topic\": 59, \"tc\": 0.007573640634369094}]}}, {\"mode\": \"vega-lite\"});\n",
       "</script>"
      ],
      "text/plain": [
       "alt.Chart(...)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import altair as alt\n",
    "chart=alt.Chart(tc_df).mark_bar().encode(x='topic',y='tc',tooltip=tc_df.columns.tolist()).properties(\n",
    "    width=800)\n",
    "chart"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you hover over, you can see that topic 6 is the cut-off boundary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "chart.save('Tc_distr_v1.html')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pointwise Document TC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can decompose total correlation further. The topic correlation is the average of the pointwise total correlations for each individual document. The pointwise total correlations can be accessed through log_z."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 7.09779218e-01,  4.20599071e-01,  3.71959114e-01, ...,\n",
       "         4.46974806e-03,  3.36120146e-03,  3.41164400e-03],\n",
       "       [ 7.09779218e-01, -5.95204711e-01,  3.71959114e-01, ...,\n",
       "        -4.23488360e-03,  3.68688331e-03,  3.91125372e-03],\n",
       "       [ 7.09779218e-01,  1.99244834e+00,  3.71959114e-01, ...,\n",
       "        -4.31722060e-02,  2.43334735e-03,  3.86500588e-03],\n",
       "       ...,\n",
       "       [ 7.09779218e-01,  4.20599071e-01, -2.55576444e-01, ...,\n",
       "        -1.06394289e-03,  2.55240972e-03,  6.63173132e-04],\n",
       "       [ 7.09779218e-01, -5.95204711e-01,  3.71959114e-01, ...,\n",
       "        -7.11688698e-04,  3.52173435e-03,  2.08494015e-03],\n",
       "       [-1.84844036e+00,  4.20599071e-01,  3.71959114e-01, ...,\n",
       "        -9.75686014e-04,  2.93997932e-03,  1.96996661e-03]])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "topic_model.log_z # n_docs x k_topics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.66434364 0.55975493 0.46578686 0.32116143 0.30054004 0.27632966\n",
      " 0.25292519 0.24396952 0.17464487 0.17271964 0.14952969 0.13937257\n",
      " 0.13701273 0.12396385 0.0808805  0.08039454 0.07430696 0.07212932\n",
      " 0.06986517 0.06328751 0.0590927  0.05695351 0.05181967 0.05154655\n",
      " 0.05010663 0.049012   0.04184238 0.04123271 0.04106111 0.03852416\n",
      " 0.03625815 0.03323562 0.03219579 0.02880976 0.02870121 0.02359333\n",
      " 0.02348045 0.01698812 0.01678462 0.01648695 0.01427115 0.01267054\n",
      " 0.01255521 0.01164168 0.01114006 0.01112151 0.01057246 0.0105284\n",
      " 0.01046922 0.01031129 0.01004515 0.0098701  0.00973535 0.00956798\n",
      " 0.00923312 0.00916255 0.00826553 0.0080398  0.00796698 0.00757364]\n",
      "[0.66434364 0.55975493 0.46578686 0.32116143 0.30054004 0.27632966\n",
      " 0.25292519 0.24396952 0.17464487 0.17271964 0.14952969 0.13937257\n",
      " 0.13701273 0.12396385 0.0808805  0.08039454 0.07430696 0.07212932\n",
      " 0.06986517 0.06328751 0.0590927  0.05695351 0.05181967 0.05154655\n",
      " 0.05010663 0.049012   0.04184238 0.04123271 0.04106111 0.03852416\n",
      " 0.03625815 0.03323562 0.03219579 0.02880976 0.02870121 0.02359333\n",
      " 0.02348045 0.01698812 0.01678462 0.01648695 0.01427115 0.01267054\n",
      " 0.01255521 0.01164168 0.01114006 0.01112151 0.01057246 0.0105284\n",
      " 0.01046922 0.01031129 0.01004515 0.0098701  0.00973535 0.00956798\n",
      " 0.00923312 0.00916255 0.00826553 0.0080398  0.00796698 0.00757364]\n"
     ]
    }
   ],
   "source": [
    "print(np.mean(topic_model.log_z, axis=0)) #The pointwise total correlations in log_z represent the correlations within an individual document explained by a particular topic. These correlations have been used to measure how \"surprising\" documents are with respect to given topics\n",
    "print(topic_model.tcs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introducing Anchoring in the semi-supervised topic mode"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "CorEx is a discriminative model, whereas LDA is a generative model. This means that while LDA outputs a probability distribution over each document, CorEx instead estimates the probability a document belongs to a topic given that document's words. As a result, the probabilities across topics for a given document do not have to add up to 1. The estimated probabilities of topics for each document can be accessed through log_p_y_given_x or p_y_given_x.\n",
    "\n",
    "Hierarchical Topic Models The labels attribute gives the binary topic expressions for each document and each topic. We can use this output as input to another CorEx topic model to get latent representations of the topics themselves. This yields a hierarchical CorEx topic model. Like the first layer of the topic model, one can determine the number of latent variables to add in higher layers through examination of the topic TCs.\n",
    "\n",
    "Anchored CorEx is an extension of CorEx that allows the \"anchoring\" of words to topics. When anchoring a word to a topic, CorEx is trying to maximize the mutual information between that word and the anchored topic. So, anchoring provides a way to guide the topic model towards specific subsets of words that the user would like to explore.\n",
    "\n",
    "Anchoring a single set of words to a single topic. This can help promote a topic that did not naturally emerge when running an unsupervised instance of the CorEx topic model. For example, one might anchor words like \"snow,\" \"cold,\" and \"avalanche\" to a topic if one suspects there should be a snow avalanche topic within a set of disaster relief articles.\n",
    "\n",
    "Anchoring single sets of words to multiple topics. This can help find different aspects of a topic that may be discussed in several different contexts. For example, one might anchor \"protest\" to three topics and \"riot\" to three other topics to understand different framings that arise from tweets about political protests.\n",
    "\n",
    "Anchoring different sets of words to multiple topics. This can help enforce topic separability if there appear to be chimera topics. For example, one might anchor \"mountain,\" \"Bernese,\" and \"dog\" to one topic and \"mountain,\" \"rocky,\" and \"colorado\" to another topic to help separate topics that merge discussion of Bernese Mountain Dogs and the Rocky Mountains."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0: products,value,establishments,industry,total,reported,manufacture,statistics,increase,industries\n",
      "1: fig,piece,inch,cut,mold,pieces,surface,inches,cement,dry\n",
      "2: wage,earners,cent,age,employed,proportion,females,engaged,gainful,males\n",
      "3: precinct,transfers,museum,assignments,patrolmen,exhibition,annum,cable,sterling,bills\n",
      "4: new,york,north,south,jersey,west,central,massachusetts,ohio,east\n",
      "5: quantity,used,sugar,steel,purchased,iron,milk,fuel,electric,rolling\n",
      "6: boulle,marquetry,dye,bronze,stock,ebony,packages,louis,amsterdam,pint\n",
      "7: water,add,hot,boil,boiling,liquid,soda,pan,white,oven\n",
      "8: today,president,washington,government,federal,department,jan,mrs,nov,commission\n",
      "9: avenue,police,supreme,oct,act,address,tomorrow,thomas,tonight,vote\n",
      "10: union,local,members,national,employers,organized,trades,secretary,membership,shops\n",
      "11: company,market,bank,old,night,july,sales,trading,sept,reserve\n",
      "12: introduction,borne,mind,estimate,animals,ivory,ornament,compare,decorated,architecture\n",
      "13: silver,wife,currency,intelligent,worth,tells,conf,magazine,reptd,universal\n",
      "14: number,average,figures,census,shows,individual,gives,comparison,earlier,summary\n",
      "15: furniture,style,cabinet,decoration,examples,seen,rich,beauty,maker,artistic\n",
      "16: states,united,pennsylvania,leading,percentage,rank,corporations,corporate,controlled,greatest\n",
      "17: week,labor,hours,prevailing,month,months,service,april,december,representative\n",
      "18: method,described,simple,prepared,girls,food,teacher,cook,cooking,vegetables\n",
      "19: sides,gold,applied,fine,cold,till,lay,care,concrete,set\n",
      "20: american,war,world,country,society,says,movement,action,street,federation\n",
      "21: church,fashion,examination,roman,religion,popular,instances,novelty,letters,tuberculosis\n",
      "22: international,committee,executive,institute,evening,duty,boston,aug,charles,afternoon\n",
      "23: small,hand,machine,clothing,owing,doing,felt,closely,newspaper,gains\n",
      "24: trade,conditions,industrial,manufacturers,worker,home,living,garment,health,occupation\n",
      "25: soft,clean,tin,plate,chair,carefully,shape,heat,turn,remove\n",
      "26: kindred,annual,university,held,france,century,late,editor,feeling,came\n",
      "27: man,young,woman,ago,death,play,artists,justice,title,servant\n",
      "28: yesterday,organization,association,meeting,convention,john,public,appointed,announced,operative\n",
      "29: locals,court,dec,unions,education,congress,party,office,council,army\n",
      "30: table,cost,different,classes,establishment,expenses,importance,branches,considerable,separate\n",
      "31: story,german,italian,press,town,england,eighth,judge,henry,san\n",
      "32: division,prices,foreign,domestic,occupied,agricultural,exchange,mechanical,middle,varied\n",
      "33: wood,design,paper,articles,feet,designs,brass,woods,solid,light\n",
      "34: school,law,people,study,brought,problem,present,need,practical,facts\n",
      "35: question,strike,demand,said,money,commercial,financial,recent,cents,saturday\n",
      "36: operations,presented,important,principal,available,data,include,bureau,chief,arranged\n",
      "37: city,business,building,factory,chicago,price,published,philadelphia,statement,closed\n",
      "38: period,class,tables,say,left,actual,away,mass,assigned,periods\n",
      "39: work,good,way,quite,french,shall,taste,little,experience,things\n",
      "40: materials,added,covered,board,printing,avoid,animal,lumber,electrical,apparatus\n",
      "41: means,better,point,right,line,strong,thing,prevent,force,contributed\n",
      "42: time,day,life,working,common,adopted,makers,particular,subject,course\n",
      "43: make,end,kept,length,inside,fixed,frame,longer,trouble,takes\n",
      "44: house,opening,hands,film,words,written,word,sign,speak,raised\n",
      "45: greater,larger,smaller,employees,returns,somewhat,practically,contract,greatly,affected\n",
      "46: shown,times,days,known,just,modern,open,free,opened,aid\n",
      "47: coal,fourth,usual,steady,iowa,naturally,range,altogether,ranks,anthracite\n",
      "48: tho,art,effect,shop,interesting,view,cause,picture,clock,coming\n",
      "49: best,position,lines,thoroughly,sufficient,easily,clear,taken,foundation,filled\n",
      "50: kind,basis,returned,single,natural,officials,complete,type,mainly,built\n",
      "51: air,wages,square,run,factories,cross,upper,fact,rail,deep\n",
      "52: use,material,various,using,processes,finishing,machinery,machines,treated,portion\n",
      "53: workers,women,change,farm,agriculture,finished,summer,majority,parts,result\n",
      "54: form,placed,joint,purpose,year,body,true,ordinary,spirit,consists\n",
      "55: nature,serve,status,thought,unless,numerous,survey,brings,source,failure\n",
      "56: season,early,beginning,started,successful,continue,developed,limited,homes,additional\n",
      "57: according,state,formed,possible,does,hour,taking,help,account,short\n",
      "58: includes,silk,manufactures,valued,vessels,rosin,textile,exports,exclusive,wore\n",
      "59: lower,book,places,named,figure,receipts,led,newspapers,broad,rapidly\n"
     ]
    }
   ],
   "source": [
    "#to automatically generate anchor words: for each label in a data set, \n",
    "#we find the words that have the highest mutual information with the label.\n",
    "# we took a very simple to automatically generate the anchor words to create a semi-supervised model\n",
    "anchor_words=[]\n",
    "for n,topic in enumerate(topics):\n",
    "    topic_words,_,_ = zip(*topic)\n",
    "    anchor_words.append(list(topic_words[:3]))\n",
    "    print('{}: '.format(n) + ','.join(topic_words))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "#manually adjust the anchor words based on the observation\n",
    "\n",
    "anchor_words[0]=['manufacturing','mechanical','industry']\n",
    "anchor_words[1]=['water','wood','boil']\n",
    "anchor_words[2]=['workers','organizations','committee']\n",
    "anchor_words[3]=['furniture','wood','carve','carpenters','workmen'] \n",
    "anchor_words[7]=['forgemen','harmmermen']\n",
    "anchor_words[8]=['market','trade','stock','sales']\n",
    "anchor_words[11]=['university','institute','education','professional'] \n",
    "anchor_words[13]=['forestry','agriculture','husbandry'] \n",
    "anchor_words[14]=['electric','power','mechanic','engines'] \n",
    "anchor_words[18]=['iron','blacksmiths','steel']\n",
    "anchor_words[20]=['art','gallery','artists']\n",
    "anchor_words[26]=['machinery','machinists','apparatus']\n",
    "anchor_words[31]=['meat','butchers','slaughtering']\n",
    "anchor_words[32]=['textile','fabrics','clothing','tailors']\n",
    "anchor_words[37]=['family','house','personal']\n",
    "anchor_words[38]=['coal','miners']\n",
    "anchor_words[39]=['court','justice','lawyers','judge']\n",
    "anchor_words[47]=['household','domestic','house']\n",
    "anchor_words[49]=['public','school','teachers']\n",
    "anchor_words[51]=['clerks','shop']\n",
    "anchor_words[56]=['book','publishing','authors']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Predict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_occsc = pd.read_csv('../../OCC_pairs.csv').rename(columns={'OCC_DES':'Full Occupation'})\n",
    "assert(df_occsc['Full Occupation'].nunique() == len(df_occsc))\n",
    "occ_original = list(set(list(df_occsc['Full Occupation'])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_predict_result(results):\n",
    "    pairs=[]\n",
    "    for i in range(len(results)):\n",
    "        for j in range(len(results[i])):\n",
    "            if results[i][j]==True:\n",
    "                pairs.append([i,j])\n",
    "    return pairs\n",
    "\n",
    "#might change the threshold\n",
    "#get log_p_y_given_x for each doc\n",
    "def get_predict_proba(results_proba):\n",
    "    pairs_proba={i:[] for i in range(len(results_proba))}\n",
    "    #i: document num\n",
    "    for i in range(len(results_proba)):\n",
    "        #j: topic num for each doc\n",
    "        for j in range(len(results_proba[i])):\n",
    "            #threshold\n",
    "            if results_proba[i][j] >= 1e-3:\n",
    "                pairs_proba[i].append((j,results_proba[i][j]))\n",
    "    return pairs_proba\n",
    "\n",
    "def count_topics(pairs):\n",
    "    available={}\n",
    "    for a,b in pairs:\n",
    "        if b not in available.keys():\n",
    "            available[b]=1\n",
    "        else:\n",
    "            available[b]+=1\n",
    "    return available \n",
    "\n",
    "def get_topic_content(pairs,topic):\n",
    "    result=[]\n",
    "    for a,b in pairs:\n",
    "        if b == topic:\n",
    "            result.append(occ_list[a])\n",
    "    return result \n",
    "\n",
    "def get_topic_list(pairs):\n",
    "    result={}\n",
    "    for a,b in pairs:\n",
    "        if b not in result.keys():\n",
    "            result[b]=[a]\n",
    "        else:\n",
    "            result[b].append(a)\n",
    "    return result\n",
    "\n",
    "def save_topic(pairs):\n",
    "    result_dic=get_topic_list(pairs)\n",
    "    with open(\"train_result_30.txt\",'w') as f:\n",
    "        for m,n in result_dic.items():\n",
    "            topic_words,_,_ = zip(*anchored_topic_model.get_topics(topic=m))\n",
    "            title=str(m)+\":\"+(','.join(topic_words))\n",
    "            f.write(title+'\\n')\n",
    "            for file in n:\n",
    "                occ_title=occ_list[file]\n",
    "                f.write(occ_title+'\\n')\n",
    "            f.write(\"\\n \\n \\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING: Some words never appear (or always appear)\n"
     ]
    }
   ],
   "source": [
    "m1 = ct.Corex(n_hidden=10, words=words)\n",
    "m1.fit(doc_word[:len(documents_mod)], words=words);\n",
    "results=m1.transform(doc_word[len(documents_mod):])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_proba,mutual=m1.transform(doc_word[len(documents_mod):],details=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "CorEx is a *discriminative* model, whereas LDA is a *generative* model. This means that while LDA outputs a probability distribution over each document, CorEx instead estimates the probability a document belongs to a topic given that document's words. \n",
    "\n",
    "As a result, the probabilities across topics for a given document do not have to add up to 1. \n",
    "\n",
    "The estimated probabilities of topics for each document can be accessed through **`log_p_y_given_x`** or **`p_y_given_x`**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The pointwise total correlations in **`log_z`** represent the correlations within an individual document explained by a particular topic. \n",
    "\n",
    "These correlations have been used to measure how \"surprising\" documents are with respect to given topics (see references below)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "results=m1.transform(doc_word[len(documents_mod):])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'clear_output' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-88-8eab51df8be8>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mpairs_proba\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mget_predict_proba\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mresults_proba\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0mpairs_proba\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 7\u001b[0;31m \u001b[0mclear_output\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'clear_output' is not defined"
     ]
    }
   ],
   "source": [
    "#just an idea of how pairs_proba looks like\n",
    "results=m1.transform(doc_word[len(documents_mod):])\n",
    "results_proba,mutual=m1.transform(doc_word[len(documents_mod):],details=True)\n",
    "pairs=get_predict_result(results)\n",
    "pairs_proba=get_predict_proba(results_proba)\n",
    "pairs_proba\n",
    "clear_output()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train_model(topic_num,anchor_w=anchor_words):\n",
    "    anchored_topic_model= ct.Corex(n_hidden=topic_num, seed=2)\n",
    "    anchored_topic_model.fit(doc_word[:len(documents_mod)], words=words, anchors=anchor_w, anchor_strength=10)\n",
    "    print(\"-\"*12,\"topic words\",\"-\"*12)\n",
    "    for n in range(topic_num):\n",
    "        topic_words,_,_ = zip(*anchored_topic_model.get_topics(topic=n))\n",
    "        print('{}: '.format(n) + ','.join(topic_words))\n",
    "    print(\"-\"*12,\"print total correlation\",\"-\"*12)\n",
    "    print(anchored_topic_model.tc)\n",
    "    #predict\n",
    "    results=anchored_topic_model.predict(doc_word[len(documents_mod):])\n",
    "    results_proba,mutual_info=anchored_topic_model.predict_proba(doc_word[len(documents_mod):])\n",
    "    \n",
    "    pairs=get_predict_result(results)\n",
    "    pairs_proba=get_predict_proba(results_proba)\n",
    "    \n",
    "    has_categorized=[]\n",
    "    \n",
    "    #a is the document index, d is the topic number it belongs to\n",
    "    for a,b in pairs:\n",
    "        if a not in has_categorized:\n",
    "            has_categorized.append(a)\n",
    "    print(\"-\"*12,\"occ that have been categorized into certain types\",\"-\"*12)\n",
    "    print(len(has_categorized))\n",
    "\n",
    "    has_not_categorized=[]\n",
    "    for i in range(len(occ_original)):\n",
    "        if i not in has_categorized:\n",
    "            has_not_categorized.append(occ_original[i])\n",
    "          \n",
    "    print(\"-\"*12,\"has_not_categorized\",\"-\"*12)\n",
    "    print(has_not_categorized)\n",
    "\n",
    "    print(\"-\"*12,\"corresponding counts of topic\",\"-\"*12)\n",
    "    print(count_topics(pairs))\n",
    "    \n",
    "    return anchored_topic_model,pairs,pairs_proba,has_not_categorized"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_pairs_proba(anchored_model):\n",
    "    #predict\n",
    "    results=anchored_topic_model.predict(doc_word[len(documents_mod):])\n",
    "    results_proba,mutual_info=anchored_topic_model.predict_proba(doc_word[len(documents_mod):])\n",
    "    \n",
    "    pairs=get_predict_result(results)\n",
    "    pairs_proba=get_predict_proba(results_proba)\n",
    "    \n",
    "    has_categorized=[]\n",
    "    \n",
    "    #a is the document index, d is the topic number it belongs to\n",
    "    for a,b in pairs:\n",
    "        if a not in has_categorized:\n",
    "            has_categorized.append(a)\n",
    "    print(\"-\"*12,\"occ that have been categorized into certain types\",\"-\"*12)\n",
    "    print(len(has_categorized))\n",
    "\n",
    "    has_not_categorized=[]\n",
    "    for i in range(len(occ_original)):\n",
    "        if i not in has_categorized:\n",
    "            has_not_categorized.append(occ_original[i])\n",
    "          \n",
    "    print(\"-\"*12,\"has_not_categorized\",\"-\"*12)\n",
    "    print(has_not_categorized)\n",
    "\n",
    "    print(\"-\"*12,\"corresponding counts of topic\",\"-\"*12)\n",
    "    print(count_topics(pairs))\n",
    "    \n",
    "    return pairs,pairs_proba,has_not_categorized"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING: Some words never appear (or always appear)\n",
      "WARNING: Anchor word not in word column labels provided to CorEx: harmmermen\n",
      "------------ topic words ------------\n",
      "0: industry,manufacturing,mechanical,branches,branch,summary,presents,principal,chemical,growth\n",
      "1: wood,water,boil,add,dry,stain,hot,till,boiling,logwood\n",
      "2: workers,committee,organizations,garment,organize,amalgamated,jewish,lisa,examined,studied\n",
      "3: wood,furniture,workmen,carpenters,chairs,carved,decoration,ivory,decorated,carve\n",
      "4: new,york,north,jersey,carolina,massachusetts,editor,england,mexico,hampshire\n",
      "5: used,quantity,sugar,fuel,beet,cane,consumed,coal,kinds,molasses\n",
      "6: boulle,marquetry,dye,style,century,taste,renaissance,designs,louis,design\n",
      "7: san,newspapers,francisco,periodicals,stage,screen,lady,hollywood,conviction,sort\n",
      "8: trade,market,stock,sales,exchange,prices,foreign,posted,trading,sterling\n",
      "9: avenue,police,supreme,today,yesterday,street,court,john,precinct,technology\n",
      "10: teacher,lesson,lessons,girl,seen,competent,language,discuss,rome,looking\n",
      "11: acid,sulphuric,skins,hides,nitric,acids,mixed,leather,treated,pipe\n",
      "12: west,east,old,grade,eighth,village,born,eastern,miles,referring\n",
      "13: age,engaged,gainful,years,sex,persons,status,distribution,rico,porto\n",
      "14: city,department,report,business,metropolitan,issued,industrial,cities,chicago,bureau\n",
      "15: men,women,conditions,factory,jobs,experience,shop,employ,efficiency,strength\n",
      "16: states,united,leading,gives,rank,statestable,exports,higher,imports,variation\n",
      "17: book,makes,details,nature,provide,useful,property,acting,requirements,brief\n",
      "18: females,occupations,males,south,divisions,occupied,gainfully,central,geographic,atlantic\n",
      "19: ohio,pennsylvania,illinois,district,california,indiana,michigan,state,wisconsin,virginia\n",
      "20: proportion,larger,population,female,male,varied,formed,silks,considerably,large\n",
      "21: said,makers,things,modern,people,know,heart,office,improved,famous\n",
      "22: agriculture,proportions,animal,specified,forestry,husbandry,vary,mines,averages,include\n",
      "23: power,electric,according,classification,steam,motors,engines,occupational,introduction,rented\n",
      "24: party,religion,senate,house,political,books,letter,article,schools,art\n",
      "25: parts,bodies,frames,fancy,boxes,various,variety,rubber,cloths,mirrors\n",
      "26: university,students,school,college,education,universities,cambridge,oxford,professor,examinations\n",
      "27: special,stated,account,regard,schedule,object,certain,peculiar,merely,favor\n",
      "28: quite,come,showed,free,contract,near,easy,provided,numbers,merchant\n",
      "29: long,say,need,lines,clear,require,watch,heavy,seasons,trying\n",
      "30: steel,iron,furnaces,turpentine,avoid,works,rosin,blast,disclosure,pig\n",
      "31: employment,smaller,gave,period,nearly,mainly,proper,required,dropped,lowest\n",
      "32: figures,census,expenses,relative,comparison,ownership,corporate,importance,cost,earlier\n",
      "33: establishments,value,table,reported,statistics,manufacture,shows,industries,separately,censuses\n",
      "34: earners,wage,number,cent,total,increase,average,employed,decade,increased\n",
      "35: shown,indicated,shall,matter,clay,built,operators,change,course,employer\n",
      "36: month,december,nearest,maximum,representative,minimum,day,year,months,july\n",
      "37: union,local,members,locals,organization,unions,executive,strike,employers,federal\n",
      "38: labor,week,wages,paid,working,earnings,ended,overtime,majority,fatigue\n",
      "39: different,represented,scope,oils,covers,farms,mining,constitutes,attributed,exception\n",
      "40: great,times,fact,did,especially,outside,natural,order,timber,tendency\n",
      "41: mind,borne,described,method,obtained,tables,spirits,appropriate,previously,chapter\n",
      "42: president,american,war,international,national,world,convention,meeting,secretary,congress\n",
      "43: miss,picture,story,father,son,young,husband,title,radio,engineer\n",
      "44: milk,sand,butter,vegetables,mixture,cream,heat,fruits,condensed,salt\n",
      "45: life,play,attention,theatre,question,prove,feeling,success,original,army\n",
      "46: left,close,hold,receive,continue,raise,holding,danger,joined,foot\n",
      "47: finished,finishing,assigned,tin,plate,polishing,hands,cleaning,measure,turning\n",
      "48: like,interesting,film,called,think,appears,addition,eyes,skill,started\n",
      "49: trades,movement,big,record,season,evidence,manufacturers,action,difficulty,supply\n",
      "50: law,health,examination,physical,legal,scientific,human,novelty,study,seek\n",
      "51: make,necessary,better,prevent,foundation,readily,kept,applying,ones,carefully\n",
      "52: goods,cotton,woolen,silk,wool,tho,worsted,knit,hosiery,shoes\n",
      "53: facts,concerning,brought,shops,worker,girls,homes,condition,inquiry,garments\n",
      "54: products,materials,mills,included,manufactured,rolling,duplication,primarily,wire,product\n",
      "55: joint,board,adopted,projects,letters,point,contact,log,embossed,mitre\n",
      "56: fig,inch,mold,sides,cut,piece,inches,ends,pieces,surface\n",
      "57: division,middle,professional,columbia,column,twelfth,fourth,sixth,incorporated,buffalo\n",
      "58: good,care,hard,best,strong,bad,pretty,leave,taking,results\n",
      "59: cents,serve,price,tax,cooking,foods,summer,engineers,worth,meal\n",
      "------------ print total correlation ------------\n",
      "22.217659015596137\n",
      "------------ occ that have been categorized into certain types ------------\n",
      "436\n",
      "------------ has_not_categorized ------------\n",
      "[]\n",
      "------------ corresponding counts of topic ------------\n",
      "{0: 185, 22: 37, 28: 436, 1: 11, 57: 61, 8: 87, 23: 4, 54: 28, 26: 15, 50: 24, 9: 30, 41: 14, 6: 6, 12: 6, 34: 11, 3: 26, 10: 5, 33: 6, 42: 13, 24: 13, 37: 27, 56: 16, 14: 21, 18: 50, 59: 30, 19: 3, 36: 1, 49: 14, 44: 11, 11: 6, 35: 8, 53: 6, 29: 2, 43: 5, 31: 6, 38: 7, 45: 4, 52: 17, 25: 11, 13: 9, 30: 19, 47: 9, 16: 5, 32: 2, 55: 6, 27: 2, 40: 1, 21: 2, 15: 4, 48: 5, 7: 3, 17: 4, 51: 3, 2: 5, 39: 1}\n"
     ]
    }
   ],
   "source": [
    "anchored_topic_model,pairs,pairs_proba,has_not_categorized =train_model(60,anchor_words[:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save model:\n",
    "#must have hyperparameter set to \"False\" to preserve word labels!!\n",
    "model_name='m1.dat'\n",
    "anchored_topic_model.save(model_name,False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "loaded from old\n",
      "------------ occ that have been categorized into certain types ------------\n",
      "436\n",
      "------------ has_not_categorized ------------\n",
      "[]\n",
      "------------ corresponding counts of topic ------------\n",
      "{0: 185, 22: 37, 28: 436, 1: 11, 57: 61, 8: 87, 23: 4, 54: 28, 26: 15, 50: 24, 9: 30, 41: 14, 6: 6, 12: 6, 34: 11, 3: 26, 10: 5, 33: 6, 42: 13, 24: 13, 37: 27, 56: 16, 14: 21, 18: 50, 59: 30, 19: 3, 36: 1, 49: 14, 44: 11, 11: 6, 35: 8, 53: 6, 29: 2, 43: 5, 31: 6, 38: 7, 45: 4, 52: 17, 25: 11, 13: 9, 30: 19, 47: 9, 16: 5, 32: 2, 55: 6, 27: 2, 40: 1, 21: 2, 15: 4, 48: 5, 7: 3, 17: 4, 51: 3, 2: 5, 39: 1}\n"
     ]
    }
   ],
   "source": [
    "import pickle as pkl\n",
    "\n",
    "#load model:\n",
    "model_path='m1.dat'\n",
    "if os.path.exists(model_path):\n",
    "    print('loaded from old')\n",
    "    with open(model_path, 'rb') as f:\n",
    "        anchored_topic_model = pkl.load(f) \n",
    "        pairs,pairs_proba,has_not_categorized=get_pairs_proba(anchored_topic_model)\n",
    "else:\n",
    "    print('created from old')\n",
    "    anchored_topic_model,pairs,pairs_proba,has_not_categorized =train_model(60,anchor_words)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/Users/muruwu/Desktop/nyu/sophomore/2nd_sem/humanities_research/topic_modeling/Coding_files/Corex_w2v/Amanda_model'"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "os.getcwd()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#just an idea of what a topic_helper dictionary looks like\n",
    "topic_helper"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Score Imputation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dictionary formation: {Topic: Score}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#new\n",
    "topic_score_dict = {} # key: topic, value: occscore\n",
    "topic_helper={i:[] for i in range(len(anchored_topic_model.get_topics()))}\n",
    "#topic_helper: {topic_num:[(score,proba)...]}\n",
    "for job_num in pairs_proba.keys():\n",
    "    #job_num: index of each occ despcriction\n",
    "    topic_prob=pairs_proba[job_num]\n",
    "    #topic_prob: the list of topcis & prob that an occ belongs to\n",
    "    #s: sum of weights (prob)\n",
    "    occ=occ_original[job_num]\n",
    "    #occ: each original occ description (without being cleaned of confusing words)\n",
    "    score = list(df_occsc[df_occsc['Full Occupation']==occ]['OCCSCORE'])[0]\n",
    "    #debug\n",
    "    #score: occ score\n",
    "    for pair in topic_prob:\n",
    "        #pair: (topic,prob) \n",
    "        weight = pair[1]\n",
    "        #weight: proba\n",
    "        topic_helper[pair[0]].append((score,weight))\n",
    "\n",
    "for topic in topic_helper.keys():\n",
    "    lst = topic_helper[topic]\n",
    "    #lst: [(score,proba)...] of a topic_num\n",
    "    s=0\n",
    "    for pair in lst:\n",
    "        occ_score = 0\n",
    "        s += pair[1]\n",
    "    occ_score = 0\n",
    "    for pair in lst:\n",
    "        occ_score += pair[0] * (pair[1]/s)\n",
    "    topic_score_dict[topic] = occ_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "topic_score_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "topic_score_dict.dumps('topic_score_dict.pkl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load modified occ description\n",
    "\n",
    "### Load Ancestry census data and calculate occ score¶"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file = open('../../common_occupation_correction.pkl','rb')\n",
    "mapping = pickle.load(file)\n",
    "def mod_occ(occ):\n",
    "    if occ in mapping.keys():\n",
    "        return mapping[occ]\n",
    "    else:\n",
    "        return occ\n",
    "\n",
    "df_ward = pd.read_csv('../../1910_New York_Manhattan Ward 9.csv')\n",
    "df_ward['Occupation'] = df_ward['Occupation'].apply(lambda x : mod_occ(x))\n",
    "df_ward['Full Occupation'] = df_ward['Industry'] + ' ' + df_ward['Occupation']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#drop all nans from full occupation\n",
    "df_ward.drop(df_ward[df_ward['Full Occupation'].isna()].index, inplace = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 32771 full occupations are nan values\n",
    "print(df_ward['Full Occupation'].isna().sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(df_ward['Full Occupation'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute occupation score of the given occupations¶"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#new \n",
    "def plural(word):\n",
    "    if word.endswith('y'):\n",
    "        return word[:-1] + 'ies'\n",
    "    elif word[-1] in 'sx' or word[-2:] in ['sh', 'ch']:\n",
    "        return word + 'es'\n",
    "    elif word.endswith('an'):\n",
    "        return word[:-2] + 'en'\n",
    "    else:\n",
    "        return word + 's'\n",
    "    \n",
    "#debug\n",
    "nan_count=topics_emp=0\n",
    "not_in_dict_topicId=[]\n",
    "empty_score=0\n",
    "empty_score_occ_list=[]\n",
    "not_in_dict_occ_with_dup=[]\n",
    "#-------\n",
    "\n",
    "def match_occ(occ):\n",
    "     #debug\n",
    "    global nan_count\n",
    "    global topics_emp\n",
    "    global not_in_dict_topicId\n",
    "    global empty_score\n",
    "    global not_in_dict_occ\n",
    "    global empty_score_occ_list\n",
    "    global not_in_dict_occ_with_dup\n",
    "    #-------\n",
    "    \n",
    "    if type(occ) == float:\n",
    "        \n",
    "        #debug\n",
    "        nan_count+=1\n",
    "        #-------\n",
    "        \n",
    "        return 0\n",
    "    else:\n",
    "        \n",
    "        vectorizer = CountVectorizer(       \n",
    "                             stop_words='english',             \n",
    "# remove stop words\n",
    "                             lowercase=True,                   \n",
    "# convert all words to lowercase\n",
    "                             token_pattern='[a-zA-Z]',  \n",
    "# num chars > 3\n",
    "                             max_features=20000,\n",
    "                             vocabulary=fixed_vocabulary,\n",
    "                             binary=True         \n",
    "# max number of uniq words    \n",
    "                            )\n",
    "        #vectorizer = CountVectorizer(stop_words='english', max_features=20000, vocabulary=fixed_vocabulary,binary=True)\n",
    "        #debug\n",
    "#         print(occ)\n",
    "#         print()\n",
    "        vectorizer = CountVectorizer(stop_words='english', max_features=20000, vocabulary=fixed_vocabulary,binary=True)\n",
    "        doc_word1 = vectorizer.fit_transform([occ])\n",
    "        doc_word1 = ss.csr_matrix(doc_word1)\n",
    "        proba, mutual_info =anchored_topic_model.predict_proba(doc_word1)\n",
    "        pairs=get_predict_proba(proba)\n",
    "        #debug\n",
    "#         print(pairs)\n",
    "        #pairs: {0:[(topic_num,proba)...]}\n",
    "        \n",
    "        if pairs == {}:\n",
    "            \n",
    "            #debug\n",
    "            topics_emp+=1\n",
    "            #-------\n",
    "            return 0\n",
    "        else:\n",
    "            score_list = []\n",
    "            s = 0\n",
    "            for pair in pairs[0]:\n",
    "                #print(pairs[0])\n",
    "                if pair[0] in topic_score_dict.keys():\n",
    "                    score = topic_score_dict[pair[0]]\n",
    "                    weight = np.e**pair[1]\n",
    "                    s += weight\n",
    "                    score_list.append((score,weight))       \n",
    "            #debug\n",
    "                else:\n",
    "                    if pair[0] not in not_in_dict_topicId:\n",
    "                        not_in_dict_topicId.append(pair[0])\n",
    "                    not_in_dict_occ_with_dup.append(occ) \n",
    "            if score_list==[]: #when occ's score_list is completely empty (i.e. occ_score must be 0)\n",
    "                empty_score+=1\n",
    "                empty_score_occ_list.append(occ)\n",
    "                \n",
    "            #-------\n",
    "            occ_score = 0\n",
    "            for pair in score_list:\n",
    "                occ_score += pair[0] * (pair[1]/s)\n",
    "                \n",
    "            return occ_score\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_ward['Occupation Score'] = df_ward['Full Occupation'].apply(lambda x : match_occ(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(nan_count)\n",
    "print(topics_emp)\n",
    "print(not_in_dict_topicId)\n",
    "print(empty_score)\n",
    "print(empty_score_occ_list)\n",
    "print(not_in_dict_occ_with_dup)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(df_ward['Occupation Score'])\n",
    "plt.title(\"OCC scores of test data\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_ward.to_csv('with_occ.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train Word2Vec model "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_labels(label):\n",
    "    labels=[]\n",
    "    for i in range(len(label)):\n",
    "        no_value=True\n",
    "        for j in range(topic_num):\n",
    "            if label[i][j]==True:\n",
    "                labels.append(j)\n",
    "                no_value=False\n",
    "                break\n",
    "        if no_value:\n",
    "            labels.append('nothing')\n",
    "    return labels\n",
    "\n",
    "labels=generate_labels(anchored_topic_model.labels)\n",
    "label_for_use=[]\n",
    "nothing_index=[]\n",
    "for i in range(len(labels)):\n",
    "    if labels[i]=='nothing':\n",
    "        nothing_index.append(i)\n",
    "    else:\n",
    "        label_for_use.append(labels[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_words=[]\n",
    "for s in range(len(document_total)):\n",
    "    if s not in nothing_index:\n",
    "        temp=tokenize(document_total[s])\n",
    "        all_words.append(temp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "from gensim.models import Word2Vec\n",
    "startTime = time.time()\n",
    "word2vec_model = Word2Vec(all_words[:len(label_for_use)], size=200, iter=10, min_count=20)\n",
    "usedTime = time.time() - startTime\n",
    "print('形成word2vec模型共花费%.2f秒' %usedTime)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "word2vec_model.wv.most_similar('clerks')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "word2vec_model.similarity('clerical','clerks')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def get_contentVector(cutWords, word2vec_model):\n",
    "    vector_list = [word2vec_model.wv[k] for k in cutWords if k in word2vec_model]\n",
    "    contentVector = np.array(vector_list).mean(axis=0)\n",
    "    return contentVector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "\n",
    "startTime = time.time()\n",
    "contentVector_list = []\n",
    "for i in range(len(all_words[:len(label_for_use)])):\n",
    "    cutWords = all_words[i]\n",
    "    content_vector=get_contentVector(cutWords, word2vec_model)\n",
    "    contentVector_list.append(content_vector)\n",
    "contentVector_list.pop(1875)\n",
    "X = np.array(contentVector_list)\n",
    "X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import LabelEncoder\n",
    "\n",
    "labelEncoder = LabelEncoder()\n",
    "y=label_for_use[:]\n",
    "y.pop(1875)\n",
    "y = labelEncoder.fit_transform(y)\n",
    "print(X.shape, y.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "logisticRegression_model = LogisticRegression()\n",
    "logisticRegression_model.fit(X,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "contentVector_predict = []\n",
    "for i in range(len(all_words[len(label_for_use):])):\n",
    "    cutWords = all_words[len(label_for_use):][i]\n",
    "    content_vector=get_contentVector(cutWords, word2vec_model)\n",
    "    contentVector_predict.append(content_vector)\n",
    "X_predict = np.array(contentVector_predict)\n",
    "results=logisticRegression_model.predict(X_predict)\n",
    "word2vec_result=[]\n",
    "for i in range(len(results)):\n",
    "    word2vec_result.append([i,results[i]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "predict_same=0\n",
    "for i in range(len(word2vec_result)):\n",
    "    word2vec=word2vec_result[i]\n",
    "    index=word2vec[0]\n",
    "    for a,b in pairs:\n",
    "        if a==index:\n",
    "            if b==word2vec[1]:\n",
    "                predict_same+=1\n",
    "            print(occ_list[i],word2vec[1],b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "topic_score_dict_word2vec = {} # key: topic, value: occscore\n",
    "topic_helper_word2vec={i:[] for i in range(topic_num)}\n",
    "for a,b in word2vec_result:\n",
    "    occ=occ_original[a]\n",
    "    score = list(df_occsc[df_occsc['Full Occupation']==occ]['OCCSCORE'])[0]\n",
    "    topic_helper_word2vec[b].append(score)\n",
    "\n",
    "for topic in topic_helper_word2vec.keys():\n",
    "    lst = topic_helper_word2vec[topic]\n",
    "    if len(lst)!=0:\n",
    "        occ_score=sum(lst)/len(lst)\n",
    "        topic_score_dict_word2vec[topic] = occ_score\n",
    "topic_score_dict_word2vec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "def word2vec_predict(occ):\n",
    "    if type(occ) == float:\n",
    "        return 0\n",
    "    temp=tokenize(occ)\n",
    "    contentVector_list = []\n",
    "    content_vector=get_contentVector(temp, word2vec_model)\n",
    "    if type(content_vector)==np.float64:\n",
    "        return 0\n",
    "    contentVector_list.append(content_vector)\n",
    "    X_predict = np.array(contentVector_list)\n",
    "    result=logisticRegression_model.predict(X_predict)\n",
    "    if result[0] in topic_score_dict_word2vec.keys():\n",
    "        occ_score=topic_score_dict_word2vec[result[0]]\n",
    "        return occ_score\n",
    "    else:\n",
    "        return 0\n",
    "\n",
    "df_ward['Occupation Score_Word2Vec'] = df_ward['Full Occupation'].apply(lambda x : word2vec_predict(x))\n",
    "df_ward['Occupation Score_Word2Vec'] =  df_ward['Occupation Score_Word2Vec'].apply(lambda x:  round(x, 5))\n",
    "df_ward"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot histogram without the 0 values:\n",
    "selection=(df_ward['Occupation Score_Word2Vec']!=0)\n",
    "df_ward[selection]['Occupation Score_Word2Vec'].hist()\n",
    "plt.title(\"old: subset=test,without 0 scores\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "word2_score=list(df_ward['Occupation Score_Word2Vec'])\n",
    "corex_score=list(df_ward['Occupation Score'])\n",
    "difference=[]\n",
    "for i in range(len(word2_score)):\n",
    "    if (word2_score[i]!=0 and corex_score[i]!=0):\n",
    "        difference.append(abs(word2_score[i]-corex_score[i]))\n",
    "sum(difference)/len(difference)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#merge Corex with Word2Vec\n",
    "for i in range(len(corex_score)):\n",
    "    if corex_score[i]==0:\n",
    "        if word2_score[i]!=0:\n",
    "            corex_score[i]=word2_score[i]\n",
    "df_ward['Occupation Score Merged']=corex_score\n",
    "df_ward"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot histogram without the 0 values:\n",
    "selection=(df_ward['Occupation Score Merged']!=0)\n",
    "df_ward[selection]['Occupation Score Merged'].hist()\n",
    "plt.title(\"old: subset=test,without 0 scores\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Deal with Categories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import copy\n",
    "occ_layer1=[\"manufacturing industries\",\n",
    "            \"professional services\",\n",
    "            \"trade & retail & business\",\n",
    "            \"agriculture farmers\",\n",
    "            \"government clerical officials\",\n",
    "            'domestic household',\n",
    "            'transportation']\n",
    "\n",
    "occ_sort=copy.deepcopy(document_total[len(documents_mod):])\n",
    "occ_sort.sort()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "occ_category={3:[range(11)],4:[range(11,23),range(287,304)],5:[range(23,60)],0:[range(60,249)],1:[range(249,287)],\n",
    "              2:[range(304,384)],6:[range(384,436)]}\n",
    "dic_occ={}# mapping occ to its bigger category in occ_layer1\n",
    "for i in range(len(occ_sort)):\n",
    "    for key,value in occ_category.items():\n",
    "        for range_ in value:\n",
    "            if i in range_:\n",
    "                dic_occ[occ_sort[i]]=key"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#try with word2vec\n",
    "topic_category_word2vec={}\n",
    "topic_category_helper_word2vec={i:[] for i in range(topic_num)}\n",
    "for a,b in word2vec_result: # b is topic, a is occ index\n",
    "    job_title=document_total[len(documents_mod):][a]\n",
    "    job_category=dic_occ[job_title]\n",
    "    topic_category_helper_word2vec[b].append(job_category)# without probability\n",
    "for topic,lst in topic_category_helper_word2vec.items():\n",
    "    if len(lst)!=0:\n",
    "        total=len(lst)\n",
    "        temp={}\n",
    "        for a in range(len(occ_layer1)):\n",
    "            temp[occ_layer1[a]]=lst.count(a)/total\n",
    "        topic_category_word2vec[topic]=temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#try with corex\n",
    "topic_category_corex={}\n",
    "topic_category_helper_corex={i:[] for i in range(topic_num)}\n",
    "for job_num in pairs_proba.keys():\n",
    "    topic_prob=pairs_proba[job_num]\n",
    "    for pair in topic_prob:\n",
    "        job_title=document_total[len(documents_mod):][job_num]\n",
    "        job_category=dic_occ[job_title]\n",
    "        weight = pair[1]\n",
    "        topic_category_helper_corex[pair[0]].append((job_category,weight))\n",
    "\n",
    "for topic in topic_category_helper_corex.keys():\n",
    "    lst = topic_category_helper_corex[topic]\n",
    "    if len(lst)!=0:\n",
    "        total=0\n",
    "        temp={}\n",
    "        for a,b in lst: # a is category, b is proba\n",
    "            total += b\n",
    "        for i in range(len(occ_layer1)):\n",
    "            count=0\n",
    "            for a,b in lst:\n",
    "                if i==a:\n",
    "                    count+=b\n",
    "            temp[occ_layer1[i]]=count/total\n",
    "        topic_category_corex[topic]=temp\n",
    "\n",
    "topic_category_corex           "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def word2vec_category(occ):\n",
    "    if type(occ) == float:\n",
    "        return 0\n",
    "    temp=tokenize(occ)\n",
    "    contentVector_list = []\n",
    "    content_vector=get_contentVector(temp, word2vec_model)\n",
    "    if type(content_vector)==np.float64:\n",
    "        return 0\n",
    "    contentVector_list.append(content_vector)\n",
    "    X_predict = np.array(contentVector_list)\n",
    "    result=logisticRegression_model.predict(X_predict)\n",
    "    if result[0] in topic_category_word2vec.keys():\n",
    "        category_proba=topic_category_word2vec[result[0]]\n",
    "        sort_dict={k: v for k, v in sorted(category_proba.items(), key=lambda item: item[1],reverse=True)}\n",
    "        category=list(sort_dict.keys())[0]\n",
    "        return category\n",
    "    else:\n",
    "        return 0\n",
    "df_ward['Occupation Category Word2Vec'] = df_ward['Full Occupation'].apply(lambda x : word2vec_category(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def corex_category(occ):\n",
    "    if type(occ)!=float:\n",
    "        vectorizer = CountVectorizer(stop_words='english', max_features=20000, vocabulary=fixed_vocabulary,binary=True)\n",
    "        temp=[i.lower() for i in occ.split()]\n",
    "        occ=\" \".join(temp)\n",
    "        doc_word1 = vectorizer.fit_transform([occ])\n",
    "        doc_word1 = ss.csr_matrix(doc_word1)\n",
    "        proba, mutual_info =anchored_topic_model.predict_proba(doc_word1)\n",
    "        pairs=get_predict_proba(proba)\n",
    "        if pairs == {}:\n",
    "            return 0\n",
    "        else:\n",
    "            temp={i:0 for i in occ_layer1}\n",
    "            #p=0\n",
    "            for pair in pairs[0]:\n",
    "                if pair[0] in topic_category_corex.keys():\n",
    "                    proba_dict=topic_category_corex[pair[0]]\n",
    "                    for a,b in proba_dict.items():\n",
    "                        temp[a]+=pair[1]*b  \n",
    "            sort_dict={k: v for k, v in sorted(temp.items(), key=lambda item: item[1],reverse=True)}\n",
    "            category=list(sort_dict.keys())[0]\n",
    "            return category\n",
    "    else:\n",
    "        return 0 \n",
    "    \n",
    "df_ward['Occupation Category Corex'] = df_ward['Full Occupation'].apply(lambda x : corex_category(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_ward"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(df_ward['Occupation Category Corex']).count(0)"
   ]
  },
  {
   "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.8.3"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}