{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pivoted Document Length Normalization\n",
    "\n",
    "## Background\n",
    "\n",
    "In many cases, normalizing the tfidf weights for each term favors weight of terms of the documents with shorter length. The _pivoted document length normalization_ scheme counters the effect of this bias for short documents by making tfidf independent of the document length.\n",
    "\n",
    "This is achieved by *tilting* the normalization curve along the pivot point defined by user with some slope.\n",
    "Roughly following the equation:\n",
    "\n",
    "`pivoted_norm = (1 - slope) * pivot + slope * old_norm`\n",
    "\n",
    "This scheme is proposed in the paper [Pivoted Document Length Normalization](http://singhal.info/pivoted-dln.pdf) by Singhal, Buckley and Mitra.\n",
    "\n",
    "Overall this approach can in many cases help increase the accuracy of the model where the document lengths are hugely varying in the entire corpus.\n",
    "\n",
    "## Introduction\n",
    "\n",
    "This guide demonstrates how to perform pivoted document length normalization.\n",
    "We will train a logistic regression to distinguish between text from two different newsgroups.\n",
    "Our results will show that using pivoted document length normalization yields a better model (higher classification accuracy)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# Download our dataset\n",
    "#\n",
    "import gensim.downloader as api\n",
    "nws = api.load(\"20-newsgroups\")\n",
    "\n",
    "#\n",
    "# Pick texts from relevant newsgroups, split into training and test set.\n",
    "#\n",
    "cat1, cat2 = ('sci.electronics', 'sci.space')\n",
    "\n",
    "#\n",
    "# X_* contain the actual texts as strings.\n",
    "# Y_* contain labels, 0 for cat1 (sci.electronics) and 1 for cat2 (sci.space)\n",
    "#\n",
    "X_train = []\n",
    "X_test = []\n",
    "y_train = []\n",
    "y_test = []\n",
    "\n",
    "for i in nws:\n",
    "    if i[\"set\"] == \"train\" and i[\"topic\"] == cat1:\n",
    "        X_train.append(i[\"data\"])\n",
    "        y_train.append(0)\n",
    "    elif i[\"set\"] == \"train\" and i[\"topic\"] == cat2:\n",
    "        X_train.append(i[\"data\"])\n",
    "        y_train.append(1)\n",
    "    elif i[\"set\"] == \"test\" and i[\"topic\"] == cat1:\n",
    "        X_test.append(i[\"data\"])\n",
    "        y_test.append(0)\n",
    "    elif i[\"set\"] == \"test\" and i[\"topic\"] == cat2:\n",
    "        X_test.append(i[\"data\"])\n",
    "        y_test.append(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gensim.parsing.preprocessing import preprocess_string\n",
    "from gensim.corpora import Dictionary\n",
    "\n",
    "id2word = Dictionary([preprocess_string(doc) for doc in X_train])\n",
    "train_corpus = [id2word.doc2bow(preprocess_string(doc)) for doc in X_train]\n",
    "test_corpus = [id2word.doc2bow(preprocess_string(doc)) for doc in X_test]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1184 787\n"
     ]
    }
   ],
   "source": [
    "print(len(X_train), len(X_test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We perform our analysis on top k documents which is almost top 10% most scored documents\n",
    "k = len(X_test) // 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gensim.sklearn_api.tfidf import TfIdfTransformer\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from gensim.matutils import corpus2csc\n",
    "\n",
    "# This function returns the model accuracy and indivitual document prob values using\n",
    "# gensim's TfIdfTransformer and sklearn's LogisticRegression\n",
    "def get_tfidf_scores(kwargs):\n",
    "    tfidf_transformer = TfIdfTransformer(**kwargs).fit(train_corpus)\n",
    "\n",
    "    X_train_tfidf = corpus2csc(tfidf_transformer.transform(train_corpus), num_terms=len(id2word)).T\n",
    "    X_test_tfidf = corpus2csc(tfidf_transformer.transform(test_corpus), num_terms=len(id2word)).T\n",
    "\n",
    "    clf = LogisticRegression().fit(X_train_tfidf, y_train)\n",
    "\n",
    "    model_accuracy = clf.score(X_test_tfidf, y_test)\n",
    "    doc_scores = clf.decision_function(X_test_tfidf)\n",
    "\n",
    "    return model_accuracy, doc_scores"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get TFIDF scores for corpus without pivoted document length normalisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9682337992376112\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/misha/envs/gensim/lib/python3.7/site-packages/sklearn/linear_model/logistic.py:433: FutureWarning: Default solver will be changed to 'lbfgs' in 0.22. Specify a solver to silence this warning.\n",
      "  FutureWarning)\n"
     ]
    }
   ],
   "source": [
    "params = {}\n",
    "model_accuracy, doc_scores = get_tfidf_scores(params)\n",
    "print(model_accuracy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Normal cosine normalisation favors short documents as our top 78 docs have a smaller mean doc length of 1668.179 compared to the corpus mean doc length of 1577.799\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Sort the document scores by their scores and return a sorted list\n",
    "# of document score and corresponding document lengths.\n",
    "def sort_length_by_score(doc_scores, X_test):\n",
    "    doc_scores = sorted(enumerate(doc_scores), key=lambda x: x[1])\n",
    "    doc_leng = np.empty(len(doc_scores))\n",
    "\n",
    "    ds = np.empty(len(doc_scores))\n",
    "\n",
    "    for i, _ in enumerate(doc_scores):\n",
    "        doc_leng[i] = len(X_test[_[0]])\n",
    "        ds[i] = _[1]\n",
    "\n",
    "    return ds, doc_leng\n",
    "\n",
    "\n",
    "print(\n",
    "   \"Normal cosine normalisation favors short documents as our top {} \"\n",
    "   \"docs have a smaller mean doc length of {:.3f} compared to the corpus mean doc length of {:.3f}\"\n",
    "   .format(\n",
    "       k, sort_length_by_score(doc_scores, X_test)[1][:k].mean(), \n",
    "       sort_length_by_score(doc_scores, X_test)[1].mean()\n",
    "   )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get TFIDF scores for corpus with pivoted document length normalisation testing on various values of alpha."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Score for slope 0.0 is 0.9720457433290979\n",
      "Score for slope 0.1 is 0.9758576874205845\n",
      "Score for slope 0.2 is 0.97712833545108\n",
      "Score for slope 0.30000000000000004 is 0.9783989834815756\n",
      "Score for slope 0.4 is 0.97712833545108\n",
      "Score for slope 0.5 is 0.9758576874205845\n",
      "Score for slope 0.6000000000000001 is 0.9733163913595934\n",
      "Score for slope 0.7000000000000001 is 0.9733163913595934\n",
      "Score for slope 0.8 is 0.9733163913595934\n",
      "Score for slope 0.9 is 0.9733163913595934\n",
      "Score for slope 1.0 is 0.9682337992376112\n",
      "We get best score of 0.9783989834815756 at slope 0.30000000000000004\n"
     ]
    }
   ],
   "source": [
    "best_model_accuracy = 0\n",
    "optimum_slope = 0\n",
    "for slope in np.arange(0, 1.1, 0.1):\n",
    "    params = {\"pivot\": 10, \"slope\": slope}\n",
    "\n",
    "    model_accuracy, doc_scores = get_tfidf_scores(params)\n",
    "\n",
    "    if model_accuracy > best_model_accuracy:\n",
    "        best_model_accuracy = model_accuracy\n",
    "        optimum_slope = slope\n",
    "\n",
    "    print(\"Score for slope {} is {}\".format(slope, model_accuracy))\n",
    "\n",
    "print(\"We get best score of {} at slope {}\".format(best_model_accuracy, optimum_slope))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9783989834815756\n"
     ]
    }
   ],
   "source": [
    "params = {\"pivot\": 10, \"slope\": optimum_slope}\n",
    "model_accuracy, doc_scores = get_tfidf_scores(params)\n",
    "print(model_accuracy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "With pivoted normalisation top 78 docs have mean length of 2077.346 which is much closer to the corpus mean doc length of 1577.799\n"
     ]
    }
   ],
   "source": [
    "print(\n",
    "   \"With pivoted normalisation top {} docs have mean length of {:.3f} \"\n",
    "   \"which is much closer to the corpus mean doc length of {:.3f}\"\n",
    "   .format(\n",
    "       k, sort_length_by_score(doc_scores, X_test)[1][:k].mean(), \n",
    "       sort_length_by_score(doc_scores, X_test)[1].mean()\n",
    "   )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing the pivoted normalization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since cosine normalization favors retrieval of short documents from the plot we can see that when slope was 1 (when pivoted normalisation was not applied) short documents with length of around 500 had very good score hence the bias for short documents can be seen. As we varied the value of slope from 1 to 0 we introdcued a new bias for long documents to counter the bias caused by cosine normalisation. Therefore at a certain point we got an optimum value of slope which is 0.5 where the overall accuracy of the model is increased.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/misha/envs/gensim/lib/python3.7/site-packages/sklearn/linear_model/logistic.py:433: FutureWarning: Default solver will be changed to 'lbfgs' in 0.22. Specify a solver to silence this warning.\n",
      "  FutureWarning)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABDAAAAHwCAYAAABQRJ8FAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xm8LFV5L/zf4wEBFSHgiTIfE6eor+MJxqiRaJyCgZs4ROMARiXx6lUzmIAxTlfvi75JTOKME6iJExpFwfhq4hjjcDA4EiImGECUI8gkiqLr/lG1oWn20Huf3WfX7v39fj712d21VlWvrq7d9fRTq1ZVay0AAAAAQ3aDtW4AAAAAwFIkMAAAAIDBk8AAAAAABk8CAwAAABg8CQwAAABg8CQwAAAAgMGTwGBZqur5VfXWtW7HelJV51TVr01Qb0tVtaraZWe0CwCGQnwxXVV1YlW9aMK6E8UtAGtBAoOZVlUnVNVZVfXTqjp6ibon9gmEI8fmv6yfv+jyrK2q2qeq/qGqvl9V36yq31mk7t5VdVJVXdhPz5+nzjOq6r/69Z1ZVbcZKftffdllVbWtqu49UvasqvpKVV3e13nW2Ho/WlXb+2W/OLq/VdWvVtWXq+qSqrqofz8HjJR/taquGJmurqr3j5T/Rv/aV1TVp6vq9iNlu/X78req6ntV9aqq2nWk7A39dru8qs6oqoeMLHv7/n1+r58+MrbuD46160dV9eVJtmdV7VdVp/TtalW1ZWy5R/bv5cqq+thY2W2q6n399ry4qj5UVbed7zMHWE1VdZeqOr3/bjq9qu6yQL1Fv1/nqX90/134srH5R/bzT1zlt8Iqqs5L+mP4Rf3jWqTun1XVf/cxwdur6qYj5Qse86vqPmNlV/T7x8P68jv2x8TvVlUbe90l98n+2HtmX/61qvofI2WPqi62vrS6GOqksXa/taou6N/Tf1TVk0bK5k7Wjbb7z+fZNvv0x/ZPTbrsYvFCX77g/2y/TV5TVd/p44n310j8xbBIYDDrvpjkfyb5woT1/yPJ4+eeVNcb4pFJvrH6TVs/+oPs0L8vXpnkR0lunuQxSV5dVXdYoO7LktwoyZYkhyZ5XFU9Ya6wP9g+McnhSW6S5KFJvtuX3SPJ8UkenmSvJG9I8g9VtWlu8XT70M8keXCSp1XVo0Ze+xlJ9mut3TTJMUneWlX79WVfS/Kg1treSfZP8vUkr55bsLV2h9baTVprN0myZ5Jzk7yrb9etk/xdkt9PsneS9yc5pa7t0XNskq1J7pjkNknuluQ5fdku/bru27+n5yR5Z12bTPhW/373SXKzJKckeftIux4y166+bZ+ea9dS2zPJT5P8Y5KHZX4XJ/nrdNt83N59W26b7nP/XJL3LbAegFVRVTdM913z1nTf9ScleV8/f9xS36/z+UaSR9Z1e2QelS5G2dBq+L1Uj0nyP5LcOcmdkvxGkt9boO7jkzwuyb3SHfP3SPLyucLFjvmttU+OHXcfmuSKdMfTJPlxknemO/aOW3Sf7H+4vzXJHya5aZJnJfn7qvrZfvl/SXKv1tpeSX6uX99o757/N8mWPs45IsmLquruY23Ye6T9/3ueNr4kyZnzbrWFl10wXpjgf/YZSe6Z7jPbP8n3MvJZMDCtNZPpelOSP01yfpLLk5yV5P79/OcneetIvSOSfDXJJUk+luQXRsrOSXJcuh9l30vypiS7j5Q/NMkZ/bKfTnKnKb6fTyU5eok6Jyb5iyTfSfIzI2384Ojy6RJ/z0nyzSQXJnlzkr1G1vO4vuyiJH/Wb4dfG1n22HTByUXpDi779GVbkrQkuyzQvrnlLu+36W+OlT853Zf9XPnd+vkHJXlPku39a75igc/yOq/ff54vTneg+kGSWyV5wshr/GeS3xtrw5H9Z3pZ39YHJ3lEktPH6v1hkvet4ud743TJi9uMzHtLkuMXqP/dJL848vzZST458hmdm36fn2fZ307yubHXbumSEvPV/9skL1+g7NAkP0xy6Dxlu6ULAr62wLL37T+HG/fPn5bk1JHyG/Sf29z/7rYkjxgp/50k5y6yTb+U5GHzzN8lyVOTXLnAcluS/CRd8LLk9hxbb5tbbp7yJyX52BLr2Kdfx76rtW+ZTKbVnTID8UWSB/bvoUbm/XeSB0+4/Lzfr33Z0elijn9Mcng/b58k307y/yU5ccJtdNd0J28uT/KOdEnnF02yjTISt8zTvsOT/Fu64/y5SZ4/Vn7vfn2X9OVH9/P3SPKX6eKjS/v3uEeSw5KcN7aOa16/3y9OTvfD87L+WHBokn/tX+OCJK9IcsOR5e+Q5MPpftB+J90x/hZJrhw9PqRL5G9Psusq7hufTnLMyPMnJvnMAnVPTvKskee/nC4muNE8da9zzJ+n/E1J3jTP/FslacvZJ5PcI8mFY+Xbk9xznuVuki4OPm2B9d62/4we2T/fkkVi3ZHt8K/pYs5Pjcxfctm+3vXihSzxP5vuZNFLx/bzs1ZrvzCt7jT0M6qsgeq6YD8t3Q+8PZM8KN3BZLzebZK8Lckzk2xOclqS94+dgXhMv/zPpzvr+5x+2bsmeWO6rPS+SV6b7mzxbgu06UvVdaufb3rVarzv3g/TZWjnzpg/Pt0X86ij++lX02Web5Lu4JnqutW/Ol0SY//+vR04suz/SpeZv2+uzfC+csK2fSPJfdJly1+QkTP3VfWIdAf5x6fLlh+R5KK+V8AH0gUMW5IckJEz5xN4XLqzCXvm2oTNQ/vXeEKSl1XV3fo2HJpuWz0r3ZnxX0m335yS5JZV9Qtj6x3frunX86pFPusvLdDO2yS5urU2enbqi+mCmIXU2OM79o8P7Kc7VtW51V328IKRHigfTLKpqu7Rb9/fTRcEfnue91LpPrOvjs3/QFX9MMln0wWd20bKDq6qS9IlH/44yUsXaP9RSd7dWvv+Iu9p9H3NV35gVe01T7tvnm6bjrf7knT/Iy9P8n8WaNfj0yWDzumfL7U9V9OvJPl2a+2iKawb2EEzFF/cIcmXWv9Lp/elLH7MmXu9eb9f5/HmXNsj9FHpYpOrRtaz4Dbqt9N70yXy90l31v5hI8suaxuN+X7frr3T/ch7SvWXF1TVIemOkS/v23SXdMfHpDtBdPd0P073SfIn6XrgTeLIdD/2907X0/AnSf4gXY/Aeya5f7retqmqPZN8JF0CaP90P+D/qbX27XTH20eOrPdxSd7eWvvx+AtW1e8ssl9cUlUHL9DWO6SLP+YsNxbZLcmt56k33zF/rq03TtdL8qRFXmfhBlx/n9yW5MyqOqKqNvWf71Xp9vG5Ze5dVZemS6o8LF3Ph9F1vqqqrkzy7+kSGKeNvew3q+q8qnpTVd1sZLlN6WLqp6VLVsxn3mWXsNT/7BuS3Kuq9q+qG6X7fvnghOtmZ1vrDIppeFO6L/sLk/xaxrLSGTlDkuTPk7xzpOwG6bKbh/XPz0ny+yPlv57kG/3jVyf532PrPivJfaf0nibtgfGidGcP/jXdgfI76c4QjPbA+Kck/3Nkudum66q3S5LnpjsYzpXN9QyYO5NwZkbORCfZb2TZLZkgszyy7BlJjuwffyjJM+apc890WfPrrTOT9cB44RJteO/c66YLgF62QL1XJ3lx//gO6RI3u63i53ufdD9cR+c9OQucsU93Juc96RIzt0qXHLqqL/vlfjuc2u8DW9J1231yX17pzub8OMnVGevNMfY6L0gXvFzvvSbZNclDkvzhAsvuk+5M5S/NU3ajdGeiDhuZd7t0geVhSW6Y7v/zp0mO68tflK43zeZ0Z6I+m3l6jvTt+kiS1y7QrhunCxQPX6D87Iz8ry21PUfq7VAPjHRJkvOTPHq19iuTybS6U2Ykvujb9/axeX+Xsd4I8yy36PdrX+foXNs74TvpTlp8Jt1lBi9K3wNjsW2ULpn7rVz3bPOn0/fAWGobZZEeGPO096/TH/vT9Yr5h3nqzPUIvPM8ZYdl6R4Yn1iiDc+ce90kj07ybwvU++0k/9I/3pTuxMP1ekDu4L7xkyS3G3l+6/7YVvPUfVK64+GW/nM+pa97z7F61zvmj5U/Lsl/LfAai/bAWGifTNdz5Ip0cc6VWfiYf0D/Gd1mnrJN6WLq56T/f0930m9rumP+zdMlpj40sswfJHn16P/CSNmiy45t14+NzVv0f7bf/m/vt//V6XoZ7bOa+4Zp9SY9MLie1trZ6Q4Gz09yYXWDCu0/T9X9052Vn1vup+m6C44OenPuyONv9sskySFJ/mg0m53uUof5Xmenaq19Kt2PvD9L8oHW2g/GqlznffeP575M98/Ie25dpnz0bPAh6cZLmHvPZ6Y72N18qXZV1eOrG2hpbtk7pjv7kHTbbr5xOg5K8s3W2tVLrX8Bo59fquohVfWZ6gY4uiRd0LhUG5LurMDv9D0SHpcu6LpqgborcUW6XiGjbpruzMB8np4umPp6urNab0tyXl8293m/tLV2Set6Erw23XtNuoP6E9IlYm6Y5LFJPjD+P1JVT0t3lurw+d5ra+3HrbUPJnlgVR0xT/nFufYazfFrfn8rXdfYj4/U//d0Z2heke5sx83Sda+ee18vTndAPiNdIPvedEmY74y0+Qbpztj9KN3Zj+vp9+nXJHnzyPWwc8vfO11y5OSR2Uttzx1WVZuT/P9JXtVae9tqrRdYXTMUXyz3mDPR9+uoPvY4Nd2Pv31ba/8yVmWxbbR/kvNb/8usNxq3rHgb9b0P5wajvjTduEtLxQE3S7L7AmWTGI9FbtP3ZPx2VV2WrkfgJLHI+5LcvqpumeQBSS5trX1uhW1ayPi+cdMkV4x9FnPemC7++Fi63g8f7eefN1bvesf8MUclefMCr7GghfbJ6u5A89Jce0LkvkleX/MMVNtaOz9db5fr9e5trf2kj6kPTPKUft4VrbVtrbWrW2vf6V/3gVW1Z/9d8PR08ff1LLbsBG93qf/ZV6br/bJvuhM174keGIMlgcG8Wmt/31q7d7qDXEs3mM64b/XlSa7pLn9QujMAcw4aeXxwv0zSHYxe3Frbe2S60UI/Pur6IzGPTq9Z+Ttd0FuT/FHmv8zhOu873fu6Ot0PwQsy8p77bmj7jtQ9N8lDxt737v0BYEF9t8zXpfuy3rd1gzx+Jdd2PTw3XTfacecmOXieH8BJd7b+RiPPbzFPnWsOhn3X0nen6wZ6874Np03QhrTWPpPuAHmfdGMvvGW+ev3rvGaRz3qhLrf/kWSX6gaynHPnLNBFt7V2cWvtMa21W7TW7pDuu3AuiDmrb+toIDD6+C7pElv/0Vr7aWvtH9N97r888h5+N92YJfdvrY0HIuN2yQLbrS/72Vz/oDtvsNJaO7m1dsfW2r5JnpfurM7n+7IftNae1lo7oLX2c+kSa6f3Qe/c/+8b0iXTHtbm6VI74gbp9p3xEbqPSvKe1toVI/OW2p47pKp+Jl3y4pTW2otXa73AdMxIfPHVJHfq2zXnTlngmLPM79dRb04Xi8x3e9nFttEFSQ4Ya9/oJQ/L2kZj/j5dT4GDWjeI42uydBzw3XSXH85Xdp1YpL+EYPNYnfFjxqvTXZpw69YNFPnssTb83HwNb639MN3YY49NdzJlsVjkMYvsF1cscgnJV9PFH3MWi0V+2lp7XmttS2vtwL7e+bnufp4skqCoqoPSJRrmvSx3kfe32D55l3S9Xrb1bfx8ul6bC91ad7E4Zqnyufd0g3Rjm+yX5GtV9e0kf5Pk0D5RtWmJZZey1P/sXdL1brq4P+n08v61J71EhZ1pWl07TOt3SndJxP3SZSJvmC5DfFJf9vxc28XztukOPPdP1wXtj9MN7HjDvvycJF9Ol3ndJ12XyP/Tl21Nd5C5R7qDzo3TXUu55yq/lxumy/r/S7pLCnZPcoMF6p6Ya7tX7tO/r+qfj15C8qR0Z+5vma4728kj2+QO6bK89+5f+y/SJTfmukL+QbpM+yH988259jKQLVngEpIkt0938L9tui55T+jX+6S+/BH99rx7vz1vlS6w2ZTuEoa/6Lfx7ulGjk66sw/fTRfU7JXuzMQ1r9+380kjbdgzXW+R+/av8ZB03Qrnttmh6QbUun+6g8kBuW43yj9Ld73hN6a037493ZmMG6frantpkjssUPfn0yWWNvXv47ujddMFAh/o3/OB6QKlJ/ZlR6VLmPxcvx0e0G+H2/Xlj0nXLfUX5nnd2/Wvt0e6/5nHpvtxPzfg6m/1n/EN+n3jnUm+MLaOA/vP/ufnWf/d+/c0t+zfj5TNnZWrJL/U7y8PHCl/TbpuyjeZZ70PSDcg3KZ0yZS/TRc8jw6at0e/ze83z/ILbs++fPdcOxjqbcfWu6kv//0kn+gfz3VFvWm6xNMrprFPmUym1Z0yI/FF3/ZvprtzwW7pTi58MyMDSY7VX/D7dZ66R6fvNt+3//65drDv0UtIFtxG/fTffft2TXds+XGuPV4vuo2y+CCeFyY5qn98aP987nM7ON0Z7Uem+9G6b5K79GWvTHcJ7v799/o9+223V/pLFPq2Pi/XjZuu2S9G2vC5dJfsVrrj6lkj22zPdAmcZ/br3zPJPUaWvVeuHRD9kCns47+frnft3DH3qxm53Gms7j7p4pFKF+d9JSMDgPZ1Fjzm9+XPzjyX2PTr3L1fb+sf7zZSvtgx/77p4qK5z+6u6U56PLB//pgkB/ePD0nXM+Q9/fOfTTdmy036z/lB6fbTI/rye+TaOGffdAPMfrQv2y3dybS56RnpEie3WGrZvnyxeGHR/9l0g6C+u98fd+236/mrvX+YVun/bK0bYBrelC4j+bn+y/3idD889u/Lnp/rjpvwm+m6qV/af4GN/gg8J9eOEn5Juu7wNxopf3C6s8Nzo0i/K6ufwPhY/8U9Oh22QN0TMzJC91jZ+F1Inpvu4L89/S2ZRuoelS5wWOguJH+Y7mB7ebqD6FzQtSWL34Xkxf3n8d0kf9Vv79EEw+/3670i3UHwrv38g9NdLnBRv+zfjizzyn77n50uwbNgAqOf99R0PU0uSXfmYnxU899Ml6S4vF/ng0bKDk43JsMLprTf7tO/z+/32/93Rsruk64L59zzR6b7AX5luksqHjS2rpv27+3y/nN+bq5NZlWSF/avcXm6QOVxI8v+V7pA8YqR6TV92S+kOxhf3m/Dz2fkbjLpBnn9r/49fLtvwyFjbTsu/R1TFthP5/5vX5uR0cpz7aCqV/b7yWNGyg7pP/sfjrX7MX35I9IlHa5It8+fmrFR/dNdd/zNue006fbsy8f/R9tI2dHzlJ848r/W+u012u6Dp7GPmUymHZsyW/HFXZOcnu4yuS+kP+b2Zc9O8sH+8aLfr/Os9+iMXPc/VnZNAmOCbbQ13WWDc3cheUeue7xecBtl8QTGw/vv+sv7z+8VY5/bfdId5+buUnJUP3+PdONlnN+39xNJ9hh5zxekS4b8ca4/BsZ4AuNXcu0x6ZPpjsmjYyXcMV2y5HvpjqXHji3/9SQfn9I+Xukuv7i4n16a6x7vrkhyn/7xbdIdj6/st+n1xsTKIsf8vvw6JwRG5m/J9Y+d50y6T6b7gX92rr3r3B+NlL043WUu3+//npD+7i7pTqB8vN+vLkuXaHzyyLKPzrVxzgXpTnDcYpL/haWWzSLxwgT/s/umGxPjwr7tn8oqj49iWr1pLiCHVVdV56T7AfyRtW4La6+q9kh3YLhba+3ra90eANYn8QU7oqr+OV0PxdevdVuA5ZvvuniAaXhKks9LXgAAa6GqfjHJ3dLdmhVYh6Y+iGd19w/+t6r6wDxlR/ejGJ/RT0+adnuAna8/W/aMdIORAayYuAJYiao6Kd0tQ5/ZWlvwjjHAsO2MHhjPSHeN+Pgo+nPe0Vpb8nZSrD+ttS1r3QaGwb4ArCJxxQbnmMJKtNaOWus2ADtuqj0wqurAdKMKu8YMANgh4goA2Nim3QPjr5P8SbpbGC3kYVX1K+luS/gHrbVzxytU1TFJjkmSG9/4xne/3e1uN422AgATOv3007/bWtu8k192VeKKRGwBs+zL5196zeP/54C91rAlDJn9ZMes9vabNK6YWgKjqh6a5MLW2ulVddgC1d6f5G2ttauq6vfS3QbrfuOVWmsnpLtFT7Zu3dq2bds2pVYDAJOoqm/u5NdbtbgiEVvAnC3HnnrN43OOP3wNW7J6Rt/Tthl5T6w++8mOWe3tN2lcMc1LSO6V5Ih+8L63J7lfVb11tEJr7aLW2lX909cnufsU2wMArF/iCgDY4KaWwGitHddaO7AfaOlRSf65tfbY0TpVtd/I0yPSDcoFAHAd4goAYGfcheQ6quqFSba11k5J8vSqOiLJ1UkuTnL0zm4PALB+iSuAaZrFS2xgPdspCYzW2seSfKx//NyR+cclOW5ntAEAmA3iCoCNRzKJZMq3UQUAAABYDTv9EhIAAGC2jJ4dT5whB6ZDDwwAAABg8CQwAAAAgMGTwAAAAAAGTwIDAAAAGDyDeAIAAMAUuP3r6tIDAwAAABg8CQwAAABg8CQwAAAAgMEzBgYAAMAaMD4CLI8eGAAAAMDgSWAAAAAAgyeBAQAAAAyeBAYAAAAweAbxBAAAWCEDccLOowcGAAAAMHgSGAAAAMDgSWAAAAAAgyeBAQAAAAyeBAYAAAAweBIYAAAAwOBJYAAAAACDJ4EBAAAADJ4EBgAAADB4EhgAAADA4ElgAAAAAIMngQEAAAAM3i5r3QAAAGB92HLsqdc8Puf4w9ewJcBGJIEBAAAA69xGSDC6hAQAAAAYPAkMAAAAYPAkMAAAAIDBMwYGAAAMzEa4lh1gufTAAAAAAAZPAgMAAAAYPJeQAAAAwE42eqlY4nKxSeiBAQAAAAyeBAYAAAAweBIYAAAAwOBNPYFRVZuq6t+q6gPzlO1WVe+oqrOr6rNVtWXa7QEA1i9xBQBsXDujB8Yzkpy5QNkTk3yvtXarJC9L8pKd0B4AYP0SVwDABjXVBEZVHZjk8CSvX6DKkUlO6h+fnOT+VVXTbBMAsD6JKwBgY5t2D4y/TvInSX66QPkBSc5Nktba1UkuTbLveKWqOqaqtlXVtu3bt0+rrQDAsK1KXJGILQBgPZpaAqOqHprkwtba6Tu6rtbaCa21ra21rZs3b16F1gEA68lqxhWJ2AIA1qNdprjueyU5oqp+PcnuSW5aVW9trT12pM75SQ5Kcl5V7ZJkryQXTbFNAMD6JK4AsuXYU695fM7xh69hS4C1MLUeGK2141prB7bWtiR5VJJ/HgsykuSUJEf1jx/e12nTahMAsD6JKwCAafbAmFdVvTDJttbaKUnekOQtVXV2kovTBSQAABMRVwDAxrFTEhittY8l+Vj/+Lkj83+Y5BE7ow0AwGwQVwCsfy4HYiV2eg8MAADYaPxYA9hx076NKgAAAMAOk8AAAAAABk8CAwAAABg8CQwAAABg8CQwAAAAgMGTwAAAAAAGz21UAQAABsatd+H69MAAAAAABk8CAwAAABg8CQwAAABg8CQwAAAAgMGTwAAAAAAGTwIDAAAAGDwJDAAAAGDwJDAAAACAwZPAAAAAAAZPAgMAAAAYPAkMAAAAYPAkMAAAAIDB22WtGwAAADCrthx76jWPzzn+8DVsCax/emAAAAAAgyeBAQAAAAyeS0gAAIAVGb08AmDa9MAAAAAABk8CAwAAABg8CQwAAABg8CQwAAAAgMGTwAAAAAAGTwIDAAAAGDwJDAAAAGDwJDAAAACAwZPAAAAAAAZPAgMAAAAYPAkMAAAAYPAkMAAAAIDBk8AAAAAABk8CAwAAABg8CQwAAABg8KaWwKiq3avqc1X1xar6alW9YJ46R1fV9qo6o5+eNK32AADrm9gCADa2Xaa47quS3K+1dkVV7ZrkU1X1wdbaZ8bqvaO19rQptgMAmA1iCwDYwKaWwGittSRX9E937ac2rdcDAGab2AIANrapjoFRVZuq6owkFyb5cGvts/NUe1hVfamqTq6qgxZYzzFVta2qtm3fvn2aTQYABkxsAQAb11QTGK21n7TW7pLkwCSHVtUdx6q8P8mW1tqdknw4yUkLrOeE1trW1trWzZs3T7PJAMCAiS0AYOOa5hgY12itXVJVH03y4CRfGZl/0Ui11yd56c5oDwCwvoktmLPl2FOveXzO8YevYUuA+fgfZTVN8y4km6tq7/7xHkkekOTfx+rsN/L0iCRnTqs9AMD6JrYAgI1tmj0w9ktyUlVtSpcoeWdr7QNV9cIk21prpyR5elUdkeTqJBcnOXqK7QEA1jexBQBsYNO8C8mXktx1nvnPHXl8XJLjptUGAGB2iC0AYGPbKWNgAAAAHWMCAKyMBAYAAIPjRz6zbHT/BiY31duoAgAA0JG4gB0jgQEAAAAMngQGAAAAMHgSGAAAAMDgGcQTAABYlwz2ChuLHhgAAADA4ElgAAAAAIMngQEAAAAMnjEwAADYYcYiAGDa9MAAAAAABk8CAwAAABg8l5AAAADAALgcb3F6YAAAAACDpwcGAAAbkjOdAOuLHhgAAADA4ElgAAAAAIMngQEAAAAMnjEwAACYCXNjWhjPAlafMWMYAgkMAAAAZoJEy2xzCQkAAAAweBIYAAAAwOC5hAQAAKLrOcDQSWAAAMAKjSY9kukkPrYce6qECkBcQgIAAACsAxIYAAAAwOC5hAQAAIB1ydg1G4seGAAAAMDgSWAAAAAzYXxQVWC2SGAAAAAAgyeBAQAAAAyeQTwBAGAdMFghsNHpgQEAALCTGKcDVk4CAwCAQfEDD4D5SGAAAADAlG059lQJ2h0kgQEAAAAMngQGAAAAMHgSGAAAAMDgTS2BUVW7V9XnquqLVfXVqnrBPHV2q6p3VNXZVfXZqtoyrfYAAOub2AIANrZdprjuq5Lcr7V2RVXtmuRTVfXB1tpnRuo8Mcn3Wmu3qqpHJXlJkt+eYpsAgPVLbAGsqdEBGM85/vA1bAlsTFPrgdE6V/RPd+2nNlbtyCQn9Y9PTnL/qqpptQkAWL/EFgCwsU11DIyq2lRVZyS5MMmHW2ufHatyQJJzk6S1dnWSS5PsO802AQDrl9gCADauaV5CktbaT5Lcpar2TvIPVXXH1tpXlrueqjomyTFJcvDBB69yKwGA9UJswY4H754SAAAgAElEQVQavQQAgIUN8ZKpnXIXktbaJUk+muTBY0XnJzkoSapqlyR7JblonuVPaK1tba1t3bx587SbCwAMnNgCmMSWY0+9ZgLWv2nehWRzf3YkVbVHkgck+fexaqckOap//PAk/9xaG7+WFQBAbMGa82MYYG1N8xKS/ZKcVFWb0iVK3tla+0BVvTDJttbaKUnekOQtVXV2kouTPGqK7QEA1jexBQzIrCdyZv39wXo0tQRGa+1LSe46z/znjjz+YZJHTKsNAMDsEFvA2vFjHhiCnTIGBgAAAMCOmOpdSAAAYCMa4uj9rF/2J2bRSnp26YEBAAAADJ4EBgAAADB4LiEBAGDmbcRBKF12sHwbcT+B9UQCAwAA1jGJCmCjcAkJAADAOrTl2FOvmdaz9d5+dh49MAAAYAD8iANYnB4YAAAAwOBJYAAAAACD5xISAAAgiQFBgWHTAwMAAAAYPAkMAAAAYPAkMAAAAIDBk8AAAAAABk8CAwAAABg8CQwAAABg8CQwAAAAgMGTwAAAAAAGTwIDAAAAGDwJDAAAAGDwJDAAAACAwZPAAAAAAAZPAgMAAAAYPAkMAAAAYPAkMAAAAIDBmyiBUVX3rqon9I83V9Utp9ssAGCWiS0AgOVaMoFRVc9L8qdJjutn7ZrkrdNsFAAwu8QWAMBKTNID4zeTHJHk+0nSWvtWkj2n2SgAYKaJLQCAZZskgfGj1lpL0pKkqm483SYBADNObAEALNskCYx3VtVrk+xdVU9O8pEkr5tuswCAGSa2AACWbZelKrTW/qKqHpDksiS3TfLc1tqHp94yAGAmiS0AgJVYNIFRVZuSfKS19qtJBBYAwA4RWwAAK7XoJSSttZ8k+WlV7bWT2gMAzDCxBQCwUkteQpLkiiRfrqoPpx8tPElaa0+fWqsAgFkmtgAAlm2SBMZ7+gkAYDWILQCAZZtkEM+TquqGSW7Tzzqrtfbj6TYLAJhVYgsAYCWWTGBU1WFJTkpyTpJKclBVHdVa+8R0mwYAzCKxBQCwEpNcQvKXSR7YWjsrSarqNkneluTu02wYADCzxBYAwLIteheS3q5zAUaStNb+I8mu02sSADDjxBYAwLJNksDYVlWvr6rD+ul1SbYttVBVHVRVH62qr1XVV6vqGfPUOayqLq2qM/rpuSt5EwDAuiK2AACWbZJLSJ6S5KlJ5m5t9skkr5pguauT/FFr7QtVtWeS06vqw621r43V+2Rr7aETtxgAWO/EFgDAsk2SwNglyd+01v4qSapqU5LdllqotXZBkgv6x5dX1ZlJDkgyHmQAABuL2AIAWLZJLiH5pyR7jDzfI8lHlvMiVbUlyV2TfHae4ntW1Rer6oNVdYcFlj+mqrZV1bbt27cv56UBgOERWwAAyzZJAmP31toVc0/6xzea9AWq6iZJ3p3kma21y8aKv5DkkNbanZO8PMl751tHa+2E1trW1trWzZs3T/rSAMAwiS0AgGWbJIHx/aq629yTqrp7kh9MsvKq2jVdgPF3rbX3jJe31i6bC2Baa6cl2bWqbjZRywGA9UpsAQAs2yRjYDwzybuq6ltJKsktkvz2UgtVVSV5Q5Iz565xnafOLZJ8p7XWqurQdAmViyZtPACwLoktAIBlWzKB0Vr7fFXdLslt+1lntdZ+PMG675XkcUm+XFVn9POeneTgfr2vSfLwJE+pqqvTnXl5VGutLfM9AADriNgCAFiJJRMYVfWIJP/YWvtKVT0nyd2q6kWttS8stlxr7VPpzqosVucVSV6xnAYDAOub2AIAWIlJxsD48/5WZfdOcv90XTdfPd1mAQAzTGwBACzbJAmMn/R/D0/yutbaqUluOL0mAQAzTmwBACzbJAmM86vqtekG1zqtqnabcDkAgPmILQCAZZskWHhkkg8leVBr7ZIk+yR51lRbBQDMMrEFALBsk9yF5Mok7xl5fkGSC6bZKABgdoktAICV0F0TAAAAGDwJDAAAAGDwJkpgVNUhVfVr/eM9qmrP6TYLAJhlYgsAYLmWTGBU1ZOTnJzktf2sA5O8d5qNAgBml9gCAFiJSXpgPDXJvZJcliStta8n+dlpNgoAmGliCwBg2SZJYFzVWvvR3JOq2iVJm16TAIAZJ7YAAJZtkgTGx6vq2Un2qKoHJHlXkvdPt1kAwAwTWwAAyzZJAuPYJNuTfDnJ7yU5LclzptkoAGCmiS0AgGXbZYI6eyR5Y2vtdUlSVZv6eVdOs2EAwMwSWwAAyzZJD4x/ShdUzNkjyUem0xwAYAMQWwAAyzZJAmP31toVc0/6xzeaXpMAgBkntgAAlm2SBMb3q+puc0+q6u5JfjC9JgEAM05sAQAs2yRjYDwzybuq6ltJKsktkvz2VFsFAMwysQUAsGxLJjBaa5+vqtsluW0/66zW2o+n2ywAYFaJLQCAlZikB0aS/GKSLX39u1VVWmtvnlqrAIBZJ7YAAJZlyQRGVb0lyc8nOSPJT/rZLYkgA9hhW4499ZrH5xx/+Bq2BNhZxBYAwEpM0gNja5Lbt9batBsDAGwIYgsAYNkmuQvJV9INrgUAsBrEFgDAsk3SA+NmSb5WVZ9LctXczNbaEVNrFQAwy8QWAMCyTZLAeP60GwEAbCjPX+sGAADrzyS3Uf14VR2S5NattY9U1Y2SbJp+0wCAWSS2AABWYskxMKrqyUlOTvLaftYBSd47zUYBALNLbAEArMQkg3g+Ncm9klyWJK21ryf52Wk2CgCYaWILAGDZJklgXNVa+9Hck6raJd292gEAVkJsAQAs2yQJjI9X1bOT7FFVD0jyriTvn26zAIAZJrYAAJZtkgTGsUm2J/lykt9LclqS50yzUQDATBNbAADLNsldSH6a5HX9BACwQ8QWAMBKLJjAqKovZ5HrUVtrd5pKiwCAmSS2AAB2xGI9MB7a/31q//ct/d/HxkBbAMDyiS0AgBVbMIHRWvtmklTVA1prdx0p+tOq+kK661cBACYitgAAdsQkg3hWVd1r5MkvT7gcAMB8xBYAwLItOYhnkicmeWNV7ZWkknwvye9OtVUAwCwTWwAAyzbJXUhOT3LnPshIa+3SqbcKAJhZYgsAYCUWuwvJY1trb62qPxybnyRprf3VlNsGAMwQsQUAsCMW64Fxo/7vnitZcVUdlOTNSW6ebmTxE1prfzNWp5L8TZJfT3JlkqNba19YyesBy7Pl2FOveXzO8YevYUuADURsAQCs2GIJjJ/v/36ttfauFaz76iR/1Fr7QlXtmeT0qvpwa+1rI3UekuTW/XSPJK/u/wIAs0dsAQCs2GIjfv96fxbjuJWsuLV2wdwZj9ba5UnOTHLAWLUjk7y5dT6TZO+q2m8lrwcADJ7YAgBYscV6YPxjulHBb1JVl43MrySttXbTSV+kqrYkuWuSz44VHZDk3JHn5/XzLhhb/pgkxyTJwQcfPOnLAgDDIrYAAFZswR4YrbVntdb2TnJqa+2mI9OeywwwbpLk3Ume2Vq7bKn6C7TlhNba1tba1s2bN69kFQDAGhNbAAA7YrFLSJIkrbUjV7ryqto1XYDxd62198xT5fwkB408P7CfBwDMKLEFALASSyYwquq3qurrVXVpVV1WVZePdftcaLlK8oYkZy5yW7RTkjy+Or+U5NLW2gUL1AUAZoDYAgBYicXGwJjz0iS/0Vo7c5nrvleSxyX5clWd0c97dpKDk6S19pokp6W7zdnZ6W519oRlvgYAsP6ILQCAZZskgfGdFQQYaa19Kt2gXIvVaUmeutx1AwDrmtgCAFi2SRIY26rqHUnem+SquZkLXHcKALAUsQUAsGyTJDBumq4L5gNH5rUkggwAYCXEFgDAsi2ZwGituXYUAFg1YgsAYCUmuQvJgVX1D1V1YT+9u6oO3BmNA5iGLceees0E7HxiCwBgJZZMYCR5U7pbku3fT+/v5wHATiPxNFPEFgDAsk2SwNjcWntTa+3qfjoxyeYptwsAmF1iCwBg2SZJYFxUVY+tqk399NgkF027YQDAzBJbAADLNsldSH43ycuTvCzdCOGfTmLwLWBwRi8tOOf4w9ewJcASxBYAwLJNcheSbyY5Yie0BYAZIZnEYsQWAMBKTHIXkpOqau+R5z9TVW+cbrMAgFkltgAAVmKSMTDu1Fq7ZO5Ja+17Se46vSYBADNObAEALNskCYwbVNXPzD2pqn0y2dgZAADzEVsAAMs2SbDwl0n+tare1T9/RJIXT69JsD7tjGv+jSsAzAixBQCwbJMM4vnmqtqW5H79rN9qrX1tus0CAGaV2AIAWImJumv2QYXAAtjw9IKB1SG2AACWa5IxMAAAAADWlAGzYAHOtAMAAAyHBAYAzDgJWQBgFriEBAAAABg8CQwAAABg8FxCwrqiGzQAAMDGJIEBACMkSgEAhsklJAAAAMDgSWAAAAAAgyeBAQAAAAyeMTBYE64xBwAAYDkkMNipRhMXwPoi8QgAwFqSwACANSAhBACwPBIYAGP8sGQh9g0AgLUjgTFwgmUAAABwFxIAAABgHdADA2acXjwAAMAs0AMDAAAAGDwJDAAAAGDwXEICsApcqgMAANMlgQFMjR/1AADAanEJCQAAADB4emCwKtbTmfaVtHU9vT8AAIBZpAcGAAAAMHhT64FRVW9M8tAkF7bW7jhP+WFJ3pfkv/pZ72mtvXBa7eG6ZqEXwtDas16MbjeA9URsAQAb2zQvITkxySuSvHmROp9srT10im0AAGbHiRFbAMCGNbVLSFprn0hy8bTWDwBsLGILANjY1noMjHtW1Rer6oNVdYc1bgsAsP6JLQBgRq1lAuMLSQ5prd05ycuTvHehilV1TFVtq6pt27dv32kNBICdwdg0q0ZsAQAzbM0SGK21y1prV/SPT0uya1XdbIG6J7TWtrbWtm7evHmnthMAWB/EFgAw29YsgVFVt6iq6h8f2rflorVqDwCwvoktAGC2TfM2qm9LcliSm1XVeUmel2TXJGmtvSbJw5M8paquTvKDJI9qrbVptYeObsoArFdiCwDY2KaWwGitPXqJ8lekuxUaAMCSxBYAsLGt9V1IAAAAAJYkgQEAAAAMngQGAAAAMHhTGwMDVsPooKPnHH/4GraEjc4AuAAAsLYkMBistfrB6IcqAADA8EhgADuNHjUAAMBKGQMDAAAAGDw9MNaRjXD22uUbAAAAzEcCA2AGbYSEJwAAG4tLSAAAAIDBk8AAAAAABk8CAwAAABg8CQzYQAySCgAArFcG8QTgOgwACgDAEOmBAcykLceeqscJAADMEAkMdpgfiQAAAEybS0gAWDMuVwEAYFJ6YAAzQU8gAACYbRIYAAAAwOC5hASmZDld48d7D+hKDwAAcF0SGMC65bIRAADYOCQwmDl6MwA7k4FIAQB2DgkMVsSZb2Al/NgHAGClJDCYGZIqwzKUz2OuHX4sAwDA+uYuJKy5LceeOpgfuzvL3HveaO+b2WHfBQBgZ9MDA3qr+YNsrX/crfXrAwAArDY9MAAAAIDBk8AAAAAABk8CAwAAABg8CQwAAABg8CQwAAAAgMGTwAAAAAAGTwIDAAAAGDwJDAAAAGDwJDAAAACAwZPAAAAAAAZPAgMAAAAYPAkMAAAAYPAkMAAAAIDBk8AAAAAABm9qCYyqemNVXVhVX1mgvKrqb6vq7Kr6UlXdbVptAQDWP7EFAGxs0+yBcWKSBy9S/pAkt+6nY5K8eoptAQDWvxMjtgCADWtqCYzW2ieSXLxIlSOTvLl1PpNk76rab1rtAQDWN7EFAGxsazkGxgFJzh15fl4/73qq6piq2lZV27Zv375TGgcArDtiCwCYYetiEM/W2gmtta2tta2bN29e6+YAAOuc2AIA1p+1TGCcn+SgkecH9vMAAFZCbAEAM2wtExinJHl8P2L4LyW5tLV2wRq2BwBY38QWADDDdpnWiqvqbUkOS3KzqjovyfOS7JokrbXXJDktya8nOTvJlUmeMK22AADrn9gCADa2qSUwWmuPXqK8JXnqtF4fAJgtYgsA2NjWxSCeAAAAwMYmgQEAAAAMngQGAAAAMHgSGAAAAMDgSWAAAAAAgyeBAQAAAAyeBAYAAAAweBIYAAAAwOBJYAAAAACDJ4EBAAAADJ4EBgAAADB4EhgAAADA4ElgAAAAAIMngQEAAAAMngQGAAAAMHgSGAAAAMDgSWAAAAAAgyeBAQAAAAyeBAYAAAAweBIYAAAAwOBJYAAAAACDJ4EBAAAADJ4EBgAAADB4EhgAAADA4ElgAAAAAIMngQEAAAAMngQGAAAAMHgSGAAAAMDgSWAAAAAAgyeBAQAAAAyeBAYAAAAweBIYAAAAwOBJYAAAAACDJ4EBAAAADJ4EBgAAADB4EhgAAADA4ElgAAAAAIMngQEAAAAMngQGAAAAMHhTTWBU1YOr6qyqOruqjp2n/Oiq2l5VZ/TTk6bZHgBgfRNbAMDGtcu0VlxVm5K8MskDkpyX5PNVdUpr7WtjVd/RWnvatNoBAMwGsQUAbGzT7IFxaJKzW2v/2Vr7UZK3Jzlyiq8HAMw2sQUAbGDTTGAckOTckefn9fPGPayqvlRVJ1fVQfOtqKqOqaptVbVt+/bt02grADB8YgsA2MDWehDP9yfZ0lq7U5IPJzlpvkqttRNaa1tba1s3b968UxsIAKwrYgsAmFHTTGCcn2T0rMeB/bxrtNYuaq1d1T99fZK7T7E9AMD6JrYAgA1smgmMzye5dVXdsqpumORRSU4ZrVBV+408PSLJmVNsDwCwvoktAGADm9pdSFprV1fV05J8KMmmJG9srX21ql6YZFtr7ZQkT6+qI5JcneTiJEdPqz0AwPomtgCAjW1qCYwkaa2dluS0sXnPHXl8XJLjptkGAGB2iC0AYONa60E8AQAAAJYkgQEAAAAMngQGAAAAMHgSGAAAAMDgSWAAAAAAgyeBAQAAAAyeBAYAAAAweBIYAAAAwOBJYAAAAACDJ4EBAAAADJ4EBgAAADB4EhgAAADA4ElgAAAAAIMngQEAAAAMngQGAAAAMHgSGAAAAMDgSWAAAAAAgyeBAQAAAAyeBAYAAAAweBIYAAAAwOBJYAAAAACDJ4EBAAAADJ4EBgAAADB4EhgAAADA4ElgAAAAAIMngQEAAAAMngQGAAAAMHgSGAAAAMDgSWAAAAAAgyeBAQAAAAyeBAYAAAAweBIYAAAAwOBJYAAAAACDJ4EBAAAADJ4EBgAAADB4EhgAAADA4ElgAAAAAIMngQEAAAAMngQGAAAAMHgSGAAAAMDgTTWBUVUPrqqzqursqjp2nvLdquodfflnq2rLNNsDAKxvYgsA2LimlsCoqk1JXpnkIUlun+TRVXX7sWpPTPK91tqtkrwsyUum1R4AYH0TWwDAxjbNHhiHJjm7tfafrbUfJXl7kiPH6hyZ5KT+8clJ7l9VNcU2AQDrl9gCADawaq1NZ8VVD0/y4Nbak/rnj0tyj9ba00bqfKWvc17//Bt9ne+OreuYJMf0T++Y5CtTafRsulmS7y5Zi8S2Wi7ba3K21fLYXpNby211SGtt8858QbHFIPj/XB7ba3K21fLYXpOzrZZnrbbXRHHFLjujJTuqtXZCkhOSpKq2tda2rnGT1g3ba3K21fLYXpOzrZbH9pqcbbVyYouVsa2Wx/aanG21PLbX5Gyr5Rn69prmJSTnJzlo5PmB/bx561TVLkn2SnLRFNsEAKxfYgsA2MCmmcD4fJJbV9Utq+qGSR6V5JSxOqckOap//PAk/9ymdU0LALDeiS0AYAOb2iUkrbWrq+ppST6UZFOSN7bWvlpVL0yyrbV2SpI3JHlLVZ2d5OJ0gchSTphWm2eU7TU522p5bK/J2VbLY3tNbkNtK7HFINhWy2N7Tc62Wh7ba3K21fIMentNbRBPAAAAgNUyzUtIAAAAAP5vO/cf61Vdx3H8+RogYDCBMkfohhqzsIyAnDR1LRcqW9MWNTYz0LZW1pZ/2CJZpW3NbKuVZWkuJzqWP0gWa2uARVotUUF+XDDkgmQykq2EcCmivvvjvC8eL9/vvfdcLpzv/X5fj+3s+/l+zvme8znv+/me+97n+zlnSHgAw8zMzMzMzMxa3rAawJB0maTtkrolLa67Pa1A0m5JWyRtlPRU1k2StEbSjnydmPWSdFvGb7OkmfW2/viTdLekfZK6SnWV4yNpYW6/Q9LCRsca7prE6iZJe7J/bZQ0r7Tumxmr7ZIuLdW3/fdU0hmS1kraJmmrpK9lvftWA33Ey/2rAUljJD0haVPG6+asP1PSujz3B/Ihlkgane+7c/3U0r4axtEKndCfBsO5RXPOK6pxbjFwzi2qcW4xcG2XV0TEsFgoHta1EzgLOAnYBEyvu111L8Bu4F296n4ALM7yYuDWLM8Dfg8IuABYV3f7T0B8LgZmAl2DjQ8wCdiVrxOzPLHucztBsboJuKHBttPzOzgaODO/myM65XsKTAZmZnk88GzGxH2rWrzcvxrHS8C4LI8C1mW/eRBYkPV3AF/O8nXAHVleADzQVxzrPr9WWTqlPw0yNrtxbtEsNs4rjj1evvY3jpVzi6GJl/vX0efeVnnFcJqBcT7QHRG7IuI14H7giprb1KquAJZmeSlwZan+3ig8DkyQNLmOBp4oEfEYxVPoy6rG51JgTUT8JyJeAtYAlx3/1p9YTWLVzBXA/RFxKCKeA7opvqMd8T2NiL0RsSHLB4FngCm4bzXUR7ya6fT+FRHxcr4dlUsAHweWZ33v/tXT75YDl0gSzeNohY7oT0PIuQXOK6pybjFwzi2qcW4xcO2WVwynAYwpwD9L71+g707aKQJYLWm9pC9m3WkRsTfL/wJOy7JjWKgan06P21dzauLdPdMWcayOyGl1H6YYzXbf6keveIH7V0OSRkjaCOyjSD53Avsj4vXcpHzuR+KS6w8A76SD4jVIjk9zzi2q8bW/Ol/7++DcohrnFv1rp7xiOA1gWGMXRsRM4HLgK5IuLq+MiKBIRKwBx6dfvwDOBmYAe4Ef1tuc1iJpHPAb4PqI+G95nfvW0RrEy/2riYh4IyJmAKdT/LrxvpqbZJ3FucUgOTYD4mt/H5xbVOPcYmDaKa8YTgMYe4AzSu9Pz7qOFhF78nUfsIKiQ77YM30zX/fl5o5hoWp8OjZuEfFiXvDeBO7irWliHR8rSaMo/mEui4iHs9p9q4lG8XL/6l9E7AfWAnMopgePzFXlcz8Sl1x/CvBvOjBeFTk+TTi3qMzX/gp87W/OuUU1zi2qa4e8YjgNYDwJTMunpZ5E8UCRlTW3qVaS3iFpfE8ZmAt0UcSl54nDC4HfZnkl8Pl8avEFwIHSlLROUjU+q4C5kibmNLS5Wdf2et3H/CmK/gVFrBbkU4rPBKYBT9Ah39O8D/BXwDMR8aPSKvetBprFy/2rMUmnSpqQ5bHAJyju7V0LzM/Nevevnn43H/hj/krXLI5W6Ij+VJVzi0Hxtb8CX/sbc25RjXOLgWu7vCJa4MmoA10onrb7LMU9O0vqbk/dC8XTcjflsrUnJhT3KP0B2AE8AkzKegG3Z/y2ALPrPocTEKNfU0wfO0xxn9YXBhMf4FqKB9V0A9fUfV4nMFb3ZSw2U1y0Jpe2X5Kx2g5cXqpv++8pcCHFFM7NwMZc5rlvVY6X+1fjeJ0HPJ1x6QK+nfVnUSQK3cBDwOisH5Pvu3P9Wf3F0Uvn9KdBxMS5Rd/xcV5x7PHytb9xrJxbDE283L+OjlVb5RXKhpiZmZmZmZmZtazhdAuJmZmZmZmZmXUoD2CYmZmZmZmZWcvzAIaZmZmZmZmZtTwPYJiZmZmZmZlZy/MAhpmZmZmZmZm1PA9gmLUpSTdJuqHudgyEpEWS3tNk3T2S5jdad4zHvLFUniqpq6/tzczMOp1zi36P6dzC7DjzAIaZtYJFQMMk4zi6sf9NzMzMbJhahHMLs7bjAQyzNiJpiaRnJf0FOKdUP0PS45I2S1ohaWLWv1fSI5I2Sdog6WxJH5P0u9JnfyZpUZZ3S7pF0kZJT0maKWmVpJ2SvlT6zNclPZnHuznrpkp6RtJdkrZKWi1pbP4CMhtYlvsd28f5zZL0qKT1edzJWf8nSbdKeiLP/6KsP1nSg5K25XmvkzRb0veBsXm8Zbn7Eb3bNjR/FTMzs+HLuYVzC7NW4gEMszYhaRawAJgBzAM+Ulp9L/CNiDgP2AJ8J+uXAbdHxIeAjwJ7B3Co5yNiBvBn4B5gPnAB0JNMzAWmAednW2ZJujg/Oy2Pdy6wH/h0RCwHngKuiogZEfFKk/MbBfwUmB8Rs4C7ge+VNhkZEecD15fO7zrgpYiYDnwLmAUQEYuBV/J4VzVr2wBiYWZm1racWzi3MGs1I+tugJkNmYuAFRHxPwBJK/P1FGBCRDya2y0FHpI0HpgSESsAIuLV3L6/46zM1y3AuIg4CByUdEjSBGBuLk/nduMo/oE/DzwXERuzfj0wtcL5nQN8AFiTbRzB25Oihxvs90LgJ3l+XZI297H/Y2mbmZlZO3JucfR+nVuY1cgDGGbW2+u8fXbWmF7rD+Xrm6Vyz/uRgIBbIuLO8ockTe21/RtAlamUArZGxJwm63v2/QaDu7YdS9vMzMysOecWzi3MhoRvITFrH48BV+a9n+OBTwJExAHgpZ57N4GrgUfz140XJF0JIGm0pJOBfwDT8/0E4JKK7VgFXCtpXO53iqR39/OZg8D4frbZDpwqaU7ud5Skc/v5zF+Bz+b204EPltYdzqmjZmZm1phzi6M5tzCrkWdgmLWJiNgg6QFgE7APeLK0eiFwRyYRu4Brsv5q4E5J3wUOA5+JiF2SHgS6gOd4a7rmQNuxWtL7gb/ldMyXgc9R/PLQzD3ZvleAOY3uVY2I1/KhXLfl1NWRwI+BrX3s9+fAUknbgL/ntgdy3S+BzZI2AEsqnKKZmVlHcG7RkHMLsxopIupug5nZcSFpBDAqIl6VdDbwCHBORLxWc9PMzMxsGHJuYS1Jgk8AAABnSURBVFYvz8Aws3Z2MrA2p3MKuM4JhpmZmR0D5xZmNfIMDDMzMzMzMzNreX6Ip5mZmZmZmZm1PA9gmJmZmZmZmVnL8wCGmZmZmZmZmbU8D2CYmZmZmZmZWcvzAIaZmZmZmZmZtbz/Aw14DwPQxm80AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x504 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as py\n",
    "\n",
    "best_model_accuracy = 0\n",
    "optimum_slope = 0\n",
    "\n",
    "w = 2\n",
    "h = 2\n",
    "f, axarr = py.subplots(h, w, figsize=(15, 7))\n",
    "\n",
    "it = 0\n",
    "for slope in [1, 0.2]:\n",
    "    params = {\"pivot\": 10, \"slope\": slope}\n",
    "\n",
    "    model_accuracy, doc_scores = get_tfidf_scores(params)\n",
    "\n",
    "    if model_accuracy > best_model_accuracy:\n",
    "        best_model_accuracy = model_accuracy\n",
    "        optimum_slope = slope\n",
    "\n",
    "    doc_scores, doc_leng = sort_length_by_score(doc_scores, X_test)\n",
    "\n",
    "    y = abs(doc_scores[:k, np.newaxis])\n",
    "    x = doc_leng[:k, np.newaxis]\n",
    "\n",
    "    py.subplot(1, 2, it+1).bar(x, y, width=20, linewidth=0)\n",
    "    py.title(\"slope = \" + str(slope) + \" Model accuracy = \" + str(model_accuracy))\n",
    "    py.ylim([0, 4.5])\n",
    "    py.xlim([0, 3200])\n",
    "    py.xlabel(\"document length\")\n",
    "    py.ylabel(\"confidence score\")\n",
    "    \n",
    "    it += 1\n",
    "\n",
    "py.tight_layout()\n",
    "py.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above histogram plot helps us visualize the effect of `slope`. For top k documents we have document length on the x axis and their respective scores of belonging to a specific class on y axis.  \n",
    "As we decrease the slope the density of bins is shifted from low document length (around ~250-500) to over ~500 document length. This suggests that the positive biasness which was seen at `slope=1` (or when regular tfidf was used) for short documents is now reduced. We get the optimum slope or the max model accuracy when slope is 0.2.\n",
    "\n",
    "# Conclusion\n",
    "\n",
    "Using pivoted document normalization improved the classification accuracy significantly:\n",
    "\n",
    "- Before (slope=1, identical to default cosine normalization): 0.9682\n",
    "- After (slope=0.2): 0.9771"
   ]
  }
 ],
 "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}