LNGTutorial.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:fcdf56a20c61123aa904ee5f9279fc094ee19466a8778538e017a5481e66eccc"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "A tutorial on Local n-grams for Authorship Attribution"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "By Dr. Robert Layton"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In this tutorial, I will show you how to use local n-grams (LNG) methods for authorship attribution.\n",
      "LNG based methods are highly effective for authorship attribution, and represent one of the significant differences in techniques between authorship studies and other document based machine learning, such as topic modelling.\n",
      "If you have any requests, or spot any bugs, please send me a message via [Github](https://github.com/robertlayton/authorship_tutorials).\n",
      "I've setup a project just for tutorials like this."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "By the way, this tool is called the ipython notebook.\n",
      "It's free, it's easy to use and allows for using python in-line.\n",
      "I won't introduce it here, but for more information, have a look [here](http://ipython.org/notebook.html).\n",
      "You can download this notebook, and run it from your own computer, allowing you to change things in-line.\n",
      "Alternatively, you can copy and paste the code segments into a new file."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "You can run this code by either downloading this workbook, using the python interactive interpreter, or by copying the lines into a file and running them with python.\n",
      "I'll be assuming knowledge of how to run python scripts, and therefore won't be covering beginner knowledge here.\n",
      "\n",
      "Please also note that this code is **slow**.\n",
      "I have really focused on making readable code, rather than fast code, and left optimisations out that could confused."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If you find any mistakes or corrections, please contact me at r.layton@icsl.com.au"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Just a configuration, ignore me\n",
      "%pylab inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].\n",
        "For more information, type 'help(pylab)'.\n"
       ]
      }
     ],
     "prompt_number": 135
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Installation Dependencies"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "I'll be using the python programming language, using numpy, scipy and scikit-learn as the 'stack'.\n",
      "You can install scikit-learn by following [this page](http://scikit-learn.org/stable/install.html), which will in turn install numpy and scipy.\n",
      "If you can run the following code, you are ready to run this tutorial.\n",
      "I've used version scikit-learn version 0.14 to do this tutorial, but it *should* work with version 0.12 and above."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import sklearn\n",
      "print(sklearn.__version__)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.15-git\n"
       ]
      }
     ],
     "prompt_number": 136
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "You will also need the [matplotlib](http://matplotlib.org/users/installing.html) library, which has the pyplot interface.\n",
      "This generates the graphs in this tutorial.\n",
      "You can test it's installed by running the next line:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import matplotlib.pyplot as plt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 137
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "0) Utility Functions"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The first utility function we need is a function for counting the frequency of character n-grams in a string.\n",
      "A character n-gram is a sequence of *n* characters in a row, present in a string.\n",
      "For these examples, we assume `n=3`, with these n-grams also called *tri-grams*.\n",
      "For example:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 3  # Length of n-gram\n",
      "suess = \"Sam! If you will let me be, I will try them. You will see.\"\n",
      "for i in range(5):  # Print the first five n-grams\n",
      "    print(\"n-gram #{}: {}\".format(i, suess[i:i+n]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "n-gram #0: Sam\n",
        "n-gram #1: am!\n",
        "n-gram #2: m! \n",
        "n-gram #3: ! I\n",
        "n-gram #4:  If\n"
       ]
      }
     ],
     "prompt_number": 138
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We won't do any preprocessing in this tutorial, so capital letters, numbers and punctuation will all remain as they are in the original documents.\n",
      "\n",
      "Now, we will use the following function to extract and count all n-grams from a given string, which I will refer to as a document.\n",
      "The counts will be optionally normalised so that the figures are the frequencies, summing to 1.0."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# A defaultdict works the same as a normal dict, except\n",
      "# that if the key is not found, it is given a default value.\n",
      "from collections import defaultdict\n",
      "\n",
      "\n",
      "def count_ngrams(document, n, normalise=False):    \n",
      "    counts = defaultdict(float)  # Default to 0.0 if the key not found\n",
      "    # Iterate through all n-grams in the document\n",
      "    for i in range(len(document) - n + 1):\n",
      "        ngram = document[i:i+n]\n",
      "        # Update the count of this n-gram.\n",
      "        counts[ngram] = counts[ngram] + 1\n",
      "    if normalise:\n",
      "        # Normalise so that sums equal 1\n",
      "        normalise_factor = float(len(document) - n + 1)\n",
      "        for ngram in counts:\n",
      "            counts[ngram] /= normalise_factor\n",
      "    return counts"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 139
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "counts = count_ngrams(suess, 3)\n",
      "# Print some examples\n",
      "for ngram in [\" I \", \"Sam\", \"ill\", \"abc\"]:\n",
      "    print(\"{0} occurs {1:.0f} time(s)\".format(ngram, counts[ngram]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " I  occurs 1 time(s)\n",
        "Sam occurs 1 time(s)\n",
        "ill occurs 3 time(s)\n",
        "abc occurs 0 time(s)\n"
       ]
      }
     ],
     "prompt_number": 140
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Let's graph the most frequent\n",
      "from operator import itemgetter\n",
      "top_ngrams = sorted(counts.items(), key=itemgetter(1), reverse=True)[:10]\n",
      "ngrams, ng_counts = zip(*top_ngrams)\n",
      "y_pos = np.arange(len(ngrams))\n",
      "\n",
      "plt.figure()\n",
      "plt.barh(y_pos, ng_counts, align='center', alpha=0.4)\n",
      "plt.yticks(y_pos, ngrams)\n",
      "plt.xlabel('Count')\n",
      "plt.title('Most frequent n-grams')\n",
      "plt.xticks(np.arange(0, max(ng_counts)+1, 1.0))\n",
      "\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAEVCAYAAAAfCXWSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGxJJREFUeJzt3XtUVXX+//HXAQEvHAVTyCtoFolcBBMTuZmm0xTlmFNU\namWTOk0Xp7TJ5YwyrdTVyOikrhTNTKtpHG2sbKamGiUuXUyWpV3EwA4iauk0mQipwOf3Rz/PVxSC\nAwcObJ6PtVrr7LP3+Xw+W/LF9rP3+bxtxhgjAIAleHl6AAAA9yHUAcBCCHUAsBBCHQAshFAHAAsh\n1AHAQgh1tAq///3v1bNnT/Xu3dvTQwHaNEIdCg0NlZ+fn/773//WeD8mJkZeXl46ePBgk9r38vLS\ngQMH6tx/8OBBLV26VPv27dPhw4eb1FdLS0lJ0bp16zw9DMCJUIdsNpsGDhyol156yfne3r17VVFR\nIZvN5pY+fuo7bgcPHtQll1yiSy65pNb9lZWVbhlDc3DXn09jVVVVebR/tD6EOiRJkydP1saNG53b\nGzZs0NSpU2uE8YkTJzR16lQFBQUpNDRUCxcudO4vLCxUcnKyAgIC1LNnT912222SpKSkJElSdHS0\n7Ha7Nm/eXKPfd955R+PGjdPhw4dlt9s1bdo0FRcXy8vLS88++6xCQkI0duxYSdKzzz6r8PBwde/e\nXT/72c9q/Avi7bff1pVXXqmAgAA98MADSk5Odl5Bp6ena8qUKc5jHQ6HvLy8VF1d7Tyve+65R717\n91bfvn31hz/8wbnvueeeU0JCgubMmaPu3btr4MCBevPNNyVJ8+bNU05Oju6//37Z7XY9+OCDF/25\nnutr48aNCgkJUc+ePbVo0aKf/FmcO7ZHjx564oknFBoaqu3btzvPZdKkSZoyZYq6deumDRs26KOP\nPtLIkSMVGBio3r1764EHHtDZs2ed7Xl5eWnVqlW6/PLL1bVrV82fP19FRUUaOXKkAgIClJaW5jz+\n+PHjuuGGGxQYGKhLLrlESUlJP/kLGa2QQbsXGhpq3nnnHRMWFma++OILU1lZafr27WuKi4uNzWYz\nxcXFxhhjpkyZYiZMmGDKysqMw+EwV1xxhVm3bp0xxpi0tDSzaNEiY4wxp0+fNnl5ec72bTabKSoq\nqrP/rKws07dvX+f2V199ZWw2m7nzzjtNeXm5qaioMK+88ooZNGiQ2bdvn6mqqjJPPPGEiY+PN8YY\nc+zYMWO3283LL79sKisrzbJly0yHDh2cY0tPTzeTJ0++qP2qqipjjDETJkwwM2fONOXl5eabb74x\ncXFxJjMz0xhjzPr1642Pj4955plnTHV1tVm1apXp3bu3s62UlBRnP7U519f06dPNDz/8YD755BPj\n5+dnvvjii1qP/+yzz4y/v7/Jy8szZ86cMbNnzzY+Pj7mP//5jzHGmAULFhgfHx/z6quvGmOMqaio\nMPn5+ebDDz80VVVVxuFwmMGDB5u//OUvNf78J0yYYE6ePGk+++wz4+vra0aPHm2++uorc+LECRMe\nHm42bNhgjDHmscceMzNnzjSVlZWmsrLS5Obm1nluaJ24UofTlClTtHHjRr399tsKDw9Xnz59nPuq\nqqq0adMmLV68WF26dFFISIgeeeQRPf/885IkX19fORwOlZaWytfXV/Hx8Q3u19RxJZienq5OnTqp\nY8eOWr16tebOnauwsDB5eXlp7ty5+vjjj3Xw4EH961//UkREhCZOnChvb2/NmjVLl156ab3tS9LX\nX3+tN954Q8uWLVOnTp3Us2dPzZo1S3/729+cx4SEhOiee+6RzWbT1KlTdeTIEX3zzTcNav+cBQsW\nyM/PT1FRUYqOjtYnn3xS63FbtmzRjTfeqPj4ePn4+Ojxxx+/aIonPj5eN954oySpY8eOio2NVVxc\nnLy8vBQSEqLp06fr3XffrfGZRx99VP7+/goPD1dkZKSuu+46hYaGqmvXrrruuuu0e/duST/+HI8c\nOSKHwyFvb2+NGjWq3nND60KoQ9KPc8NTpkzRiy++WOvUy/Hjx3X27FmFhIQ43+vfv79KS0slSX/6\n059kjFFcXJwiIiK0fv36Jo+pX79+ztfFxcV66KGHFBgY6JwakKTS0lIdOXJEffv2rfOzP6W4uFhn\nz55Vr169nG3PnDlTx44dcx5z/i+Izp07S5LKysqc7zVkXv3CNk6dOiVJ8vf3l91uV9euXVVSUnLR\nuXTq1Omiew0Xnuv+/ft1ww03qFevXurWrZvmzZt30U3v4ODgGm1euH3ufObMmaNBgwZp3Lhxuuyy\ny/Tkk0/We25oXQh1OPXv318DBw7UG2+8oYkTJ9bY16NHD/n4+MjhcDjfO3jwoDNggoODtWbNGpWW\nliozM1P33XffTz7x0hDnh2X//v21Zs0a/e9//3P+d+rUKY0cOVK9evVSSUmJ81hjTI1tf39/lZeX\nO7ePHj3qfN2vXz/nkz/n2j1x4oT27t3r8hgbo6ysTCdPntT333+vfv36qVevXjp06JBzf0VFxUUB\nfWGfv/71rxUeHq7CwkKdOHFCCxcudN4TcJW/v78yMjJUVFSk1157TUuXLnXO56NtINRRw7p167R9\n+3Z16tSpxvve3t665ZZbNG/ePJWVlam4uFjLli3T5MmTJUmbN292hlFAQIBsNpu8vH783ys4OFhF\nRUVNGtfMmTO1aNEiff7555J+vLl57qbrz3/+c3322WfaunWrKisrtXz58hrBPXToUGVnZ6ukpEQn\nTpzQ4sWLnft69eqlcePG6eGHH9bJkydVXV2toqIiZWdnN2hcjT23uqZsJk2apG3btun999/XmTNn\nlJ6eXu/0TllZmex2uzp37qx9+/Zp1apVLvV//uvXX39dhYWFMsaoa9eu8vb2lre3dwPPCq0BoY4a\nBg4cqNjYWOf2+VeFK1asUJcuXTRw4EAlJibqjjvu0LRp0yRJu3bt0tVXXy273a6bbrpJy5cvV2ho\nqKQf58bvvPNOBQYGasuWLbX2e+HV54XbEyZM0O9+9zulpaWpW7duioyM1L///W9JP/4rYvPmzXrs\nscfUo0cPFRYWatSoUc6wGjt2rG699VZFRUVp+PDhSk1NrdH+xo0bdebMGeeTNb/85S+dvxRsNttP\nju2hhx7Sli1b1L17d82aNatB51bXe5IUHh6uFStWKC0tTb1795bdbldQUJD8/PzqHE9GRob++te/\nqmvXrpo+fbrS0tJqHFNf/+e3WVhYqGuvvVZ2u13x8fH6zW9+o+Tk5FrHitbJZhpylwdoY0aPHq0p\nU6Y4f+m0VWVlZQoMDFRhYWGN+xlAXbhSh2W11euVbdu2qby8XKdOndLs2bMVFRVFoKPBCHVYlqe/\n7dlYr732mvr06aM+ffqoqKioxuOVQH2YfgEAC+FKHQAspIOnOm6r/zQGAE+qb3LFY6Eutd0bWe1d\nenq60tPTPT0MNBI/v7arIRfDTL8AgIUQ6gBgIYQ6XJaSkuLpIaAJ+PlZm8ceabTZbMypA4ALGpKb\nXKkDgIUQ6gBgIYQ6AFgIoQ4AFkKoA4CFePQbpfPmrfFk92ikoCA/PfTQnZ4eBoBaeDTUQ0Kme7J7\nNFJxMb+MgdaK6RcAsBBCHQAspEmhPnnyZK1evdq5/eGHHyo6OlpVVVVNHhgAwHVNCvWlS5dqyZIl\nOn78uKqrq/XAAw9o1apV8vb2dtf4AAAuaFKoBwUFafbs2Xr00Ue1evVqRUdHa9iwYbr77rsVFRWl\n2NhYZWVluWmoAID6NPnpl5kzZ2rDhg169913tWvXLq1cuVLe3t7as2ePCgoKNG7cOH355Zfy9fW9\n6LPbtqU7X19xRYrCwlKaOhwAsIysrCyXL4ybHOo2m00zZsxQfn6+AgMDlZeXpwcffFCSFBYWppCQ\nEBUUFCgyMvKiz6ampje1ewCwrJSUlBpLJf/xj3+s9zNuefrFZrPVKLN04dKQ1CMFgJbh9kcaExMT\n9eKLL0qS9u/fr4MHDyosLMzd3QAAauH2K/X77rtP1dXVioqKUlpamjZs2CAfHx93dAMAqIdHKx9l\nZlL5qC0qLl6jhQtZ4gFoaVQ+AoB2hlAHAAvx6CqNrPbXNgUF+Xl6CADq4NE5dQ91DQBtEnPqANDO\nUPkILqPyEdB6UfkILuNeCNB6Mf0CABZCqAOAhRDqAGAhDQ51h8OhK6+8UnfffbfCwsJ0xx136K23\n3tKoUaN0xRVX6KOPPpIknTp1StOmTdOIESMUGxur1157rdkGDwCoyaUr9aKiIs2ePVv79u1TQUGB\nNm3apLy8PGVkZGjRokWSpIULF2rMmDH68MMPtX37ds2ZM0fl5eXNMngAQE0uPf0yYMAADRkyRJI0\nZMgQjR07VpIUEREhh8MhSXrrrbe0bds2ZWRkSJJOnz6tkpISlt8FgBbgUqj7+f3f18O9vLycJeq8\nvLxUWVnp3PePf/xDl19+eb3tUc4OAOrmkXJ2Fxo/fryWL1+uFStWSJJ2796tmJiYWo+lnB0A1K3Z\ny9ldWJbu/O1zr//whz/o7NmzioqKUkREhBYsWOBKFwCAJqBIBlxGkQzAM1jQCwDaGUIdACyEUAcA\nC6HyEVxG5SOg9aLyEQC0EdwoBYB2hlAHAAuhnB1cRjk7oPWinB1cxg1uoPVi+gUALIRQBwALcUuo\np6en689//nOt+0JDQ/Xtt9+6oxsAQD3cEuoXrt7Y0H0AAPdi+gUALMSjT79Q+QgA6ubRykeNmWah\n8hEA1K0xlY/cEupUNwKA1oE5dQCwELeEemZmpp5//vla9/H0CwC0HLdMv8yYMaPOfQcOHHBHFwCA\nBmD6BQAshFAHAAuhnB1cRjk7oPWinB0AtBGUswOAdoZQBwALoZwdXEY5O6D1opwdXMYNbqD1YvoF\nACyEUAcACyHUAcBCXAp1f39/Pfroo4qIiNC1116rDz74QMnJybrsssu0bds2SVJVVZXmzJmjuLg4\nRUdHa80a5l8BoKW4FOrl5eUaM2aMPv30U9ntds2fP1/bt2/X1q1bNX/+fEnSunXrFBAQoJ07d2rn\nzp1au3atHA5Hc4wdAHABl55+8fX11fjx4yVJkZGR6tixo7y9vRUREeEM7rfeekt79+7Vli1bJEnf\nf/+9CgsLFRoaelF7lLMDgLo1ezk7Hx8f52svLy/5+vo6X1dWVjr3rVy5Utdee2297VHODgDq1phy\ndm6/UTp+/Hg9/fTTzpDfv3+/ysvL3d0NAKAWLl2pX1jF6Pztc69/9atfyeFwKDY2VsYYBQUFaevW\nrW4YKgCgPh5dpTEzk1Ua26Li4jVauJBvAwMtjVUaAaCdIdQBwEKofASXUfkIaL2ofAQAbQRz6gDQ\nzhDqAGAhVD4C2hkqV1kblY+AdoYHFKyN6RcAsBBCHQAshFAHAAtpUqgvXbpUkZGRioyM1FNPPSVJ\ncjgcioyMdB6TkZHRoOUiAQBN1+gbpfn5+Xruuee0c+dOVVdXa8SIEUpOTlZAQECN4y5c2REA0Hwa\nfaWem5uriRMnqlOnTurSpYsmTpyonJwcQhwAPKjRV+oXfl3VGCObzaYOHTqourra+X5FRUWdbVDO\nDgDq1phydo1e+2X37t2666679MEHH6i6ulpXX321XnjhBYWHh6t3794qKChQly5dlJycrOuuu04L\nFiyo2THrqQMewXr4bVdD1n5p9JV6TEyM7rrrLsXFxUmS7r33XkVHR0uS5s+fr7i4OPXp00fh4eFM\nyQBAC6HyEdDOcKXedrFKIwC0M4Q6AFgIoQ4AFkI5O6CdoRyhtVHODgDaCG6UAkA7Q6gDgIVQzg5o\nZyhnZ22UswPaGR5QsDamXwDAQgh1ALAQQh0ALKTeUG9IebqqqioNHDhQkvTdd9/J29tbubm5kqSk\npCQVFRW5c8wAgDq4fKVe2zK63t7eCgsL0+eff67c3FwNGzZM2dnZOn36tA4dOqTLLrvMLYMFAPw0\ntz39kpiYqOzsbH311VeaO3eu1q5dq+TkZA0fPrzOz1D5CADq1pjKR/WGekPL0yUlJenpp5/WkSNH\n9Pjjj2vJkiXKyspSUlJSnW2npqa7NFgAaE9SUlKUkpLi3L5w6rs29U6/BAcH65tvvtG3336r06dP\n6/XXX6/1uLi4OL333nvy9vaWn5+foqOjlZmZ+ZOhDgBwr3pD3cfHx1mebty4cXWWp/P19VX//v11\n9dVXS/rxyr2srKzGTVYAQPOinB3QzlDOru1ilUYAaGcIdQCwEEIdACyEcnZAO0M5O2ujnB0AtBHc\nKAWAdobKRwBgIVQ+AoA2Y0a9RzD9AgAWQqgDgIU0e6jfe++9+uKLL5q7GwCAWmBOfe3atc3dBQDg\n/2vSlfrmzZv1yCOPSJKeeuopZ4WjAwcOKCEhQdKP6wHn5+c3cZgAgIZoUqgnJSUpJydHkpSTk6Me\nPXro8OHDysnJUXJysqQfH5avbaleAID7NWn6JTg4WGVlZSorK9OhQ4d0++23Kzs7W7m5ubr55pvr\n/Tzl7ACgbgUFWdq/P8ulzzR5Tj0+Pl7r169XWFiYEhIStG7dOr3//vtaunRpvZ+lnB0A1C0srObF\n7uuvu6GcXX0SExO1ZMkSJScnKyYmRjt27FDHjh1lt9ub2jQAwEVNvlJPSEhQaWmpkpKS5OXlpf79\n+2vw4MHuGBsAwEWUswOANmLGDFZpBIB2hVAHAAuh8hEAWAiVjwCgjaDyEQC0M4Q6AFgI5ewAwEIo\nZwcAbQbl7ACgXSHUAcBC3Brqo0aNkiQ5HA5FRkZKkrKyspSamurObgAAdXBrqOfl5bmzOQCAi9wa\n6v7+/u5sDgDgIrc+/eJq2ToqHwFA3TxS+agpqHwEAHXzSOUjAEDrQagDgIW4NdTPn1Ov6zUAoPlQ\nzg4A2gjK2QFAO0OoA4CFUM4OACyEcnYA0EZQzg4A2hlCHQAshHJ2AGAhlLMDgDaDcnYA0K4Q6gBg\nIc0W6ucKZpxf2g4A0LyaLdRZxAsAWh7TLwBgIR59+oVydgBQN8rZAYCFUM4OANo5Qh0ALKRFnn7h\nSRgAaBnNFurff/+9JCk0NFR79uxprm4AAOdh+gUALITKRwBgIVQ+AoA2gspHANDOEOoAYCFUPgIA\nC6HyEQC0GVQ+AoB2hVAHAAtpllC//vrrnd8oPVcBCQDQ/JplTv2f//yn8zXrvgBAy2nUlfqSJUu0\nYsUKSdJvf/tbjRkzRpK0fft23XHHHRowYIC+/fZb940SANAgjQr1pKQk5eTkSJJ27dqlU6dOqbKy\nUrm5uUpOTnbrAAEADdeo6ZfY2Fjl5+fr5MmT6tixo6666irt2rVLOTk5Wr58uRYvXtygdihnBwB1\na0w5u0av/TJ27FjddNNNOn78uKKiolRQUKBnnnlGBw4c0IABA5Sfn6/u3bvLbrfr5MmTF3dssykz\nk7VfAKChZsxoxrVfEhMTlZGRoeTkZCUmJmr16tWKiYlpbHMAADdoUqgfPXpUI0eOVFBQkDp16qTE\nxMSLjuPpFwBoOR5depfpFwBouGadfgEAtD6EOgBYCOXsAMBCKGcHAG0E5ewAoJ0h1AHAQgh1ALAQ\nQh0ALIRQBwALIdQBwEIIdQCwEEIdACyEUAcACyHUAcBCCHUAsBBCHQAshFAHAAsh1AHAQgh1ALAQ\nQh0ALIRQBwALIdQBwEIIdQCwEEIdACyEUAcACyHUAcBCCHUAsBBCHQAshFCHy7Kysjw9BDQBPz9r\nI9ThMkKhbePnZ22EOgBYCKEOABZiM8YYj3Rss3miWwBo0+qL7A4tNI6LeOh3CQBYGtMvAGAhhDoA\nWEiLh/qbb76pK6+8UpdffrmefPLJlu4eTTBt2jQFBwcrMjLS00OBi0pKSjR69GgNGTJEERERWr58\nuaeHBBf88MMPGjFihIYOHarw8HDNnTu3zmNb9EZpVVWVwsLC9M4776hPnz4aPny4XnrpJQ0ePLil\nhoAmyMnJkb+/v6ZOnaq9e/d6ejhwwdGjR3X06FENHTpUZWVlGjZsmF555RX+7rUh5eXl6ty5syor\nK5WQkKCMjAwlJCRcdFyLXqnv3LlTgwYNUmhoqHx8fJSWlqZXX321JYeAJkhMTFRgYKCnh4FGuPTS\nSzV06FBJkr+/vwYPHqzDhw97eFRwRefOnSVJZ86cUVVVlbp3717rcS0a6qWlperXr59zu2/fviot\nLW3JIQDtnsPh0O7duzVixAhPDwUuqK6u1tChQxUcHKzRo0crPDy81uNaNNR5Nh3wrLKyMk2aNElP\nPfWU/P39PT0cuMDLy0sff/yxDh06pOzs7DqXe2jRUO/Tp49KSkqc2yUlJerbt29LDgFot86ePaub\nb75ZkydP1oQJEzw9HDRSt27ddP3112vXrl217m/RUL/qqqv05ZdfyuFw6MyZM9q0aZNuvPHGlhwC\n0C4ZY3TPPfcoPDxcs2bN8vRw4KLjx4/ru+++kyRVVFTo7bffVkxMTK3Htmiod+jQQStXrtT48eMV\nHh6uW2+9lbvvbchtt92m+Ph47d+/X/369dP69es9PSQ0UF5enl544QXt2LFDMTExiomJ0Ztvvunp\nYaGBjhw5omuuuUZDhw7ViBEjlJqaqjFjxtR6rMfWfgEAuB/fKAUACyHUAcBCCHUAsBBCHQAshFCH\npR09elRpaWkaNGiQrrrqKl1//fX68ssv3db+u+++q/fff99t7QFNRajDsowx+sUvfqFrrrlGhYWF\n2rVrlxYvXqyvv/7abX3s2LFD7733ntvaA5qKUIdl7dixQ76+vpo+fbrzvaioKCUkJGjOnDmKjIxU\nVFSU/v73v0uSsrKylJqa6jz2/vvv14YNGyRJoaGhSk9P17BhwxQVFaWCggI5HA5lZmZq2bJliomJ\nUW5ubsueIFALj5WzA5rbp59+qmHDhl30/ssvv6xPPvlEe/bs0bFjxzR8+HAlJSVddJzNZnOuV2Sz\n2dSzZ0/l5+dr1apVysjI0Nq1azVz5kzZ7XY9/PDDzX4+QENwpQ7LqmsBuby8PN1+++2y2WwKCgpS\ncnKyPvroo3oXnJs4caIkKTY2Vg6Hw/k+399Da0Kow7KGDBmi/Pz8WvddGMQ2m00dOnRQdXW1872K\niooax/j5+UmSvL29VVlZ6ebRAu5BqMOyrrnmGp0+fVpr1651vrdnzx4FBARo06ZNqq6u1rFjx5Sd\nna24uDj1799fn3/+uc6cOaPvvvtO27dvr7cPu92ukydPNudpAC5hTh2WtnXrVs2aNUtPPvmkOnbs\nqAEDBmjZsmUqKytTdHS0bDablixZoqCgIEnSLbfcooiICA0YMECxsbG1tnn+XHtqaqomTZqkV199\nVStXrtSoUaNa7NyA2rCgFwBYCNMvAGAhhDoAWAihDgAWQqgDgIUQ6gBgIYQ6AFjI/wMQ0EMQjnqR\nIgAAAABJRU5ErkJggg==\n"
      }
     ],
     "prompt_number": 141
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "frequencies = count_ngrams(suess, 3, normalise=True)\n",
      "# Print some examples\n",
      "for ngram in [\" I \", \"Sam\", \"ill\", \"abc\"]:\n",
      "    print(\"{0} has frequency {1:.3f}\".format(ngram, frequencies[ngram]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " I  has frequency 0.018\n",
        "Sam has frequency 0.018\n",
        "ill has frequency 0.054\n",
        "abc has frequency 0.000\n"
       ]
      }
     ],
     "prompt_number": 142
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We now need a corpus to do the classification on.\n",
      "The format I'll use for this tutorial is to use one list of strings for the \"documents\" and one array of integers as the \"classes\", or \"targets\".\n",
      "Each document has a corresponding class value, such that all documents with the same class were written by the same person.\n",
      "\n",
      "\n",
      "First, grab [this zip file](https://github.com/robertlayton/authorship_tutorials/raw/master/books.zip) containing the dataset and unzip it into a directory somewhere.\n",
      "Then update `default_folder` value below, changing the value from `None` to the path of the data."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import os\n",
      "from os.path import expanduser\n",
      "import codecs\n",
      "import numpy as np\n",
      "\n",
      "default_folder = None  # Change this to the folder name containing the dataset\n",
      "\n",
      "def get_single_corpus(folder=None):\n",
      "    if folder is None:\n",
      "        folder = os.path.join(expanduser(\"~\"), \"Data\", \"books\")\n",
      "    documents = []\n",
      "    authors = []\n",
      "    training_mask = []\n",
      "    authornum = 0\n",
      "    i = 0\n",
      "    subfolders = [name for name in os.listdir(folder)\n",
      "                  if os.path.isdir(os.path.join(folder, name))]\n",
      "    for subfolder in subfolders:\n",
      "        sf = os.path.join(folder, subfolder)\n",
      "        #print \"Author %d is %s\" % (authornum, subfolder)\n",
      "        for document_name in os.listdir(sf):\n",
      "            i += 1\n",
      "            with codecs.open(os.path.join(sf, document_name), encoding='utf=8') as input_f:\n",
      "                documents.append(cleanFile(input_f.read()))\n",
      "                authors.append(authornum)\n",
      "                training_mask.append(True)\n",
      "        authornum += 1\n",
      "    # There is an author that appears only once. To simplify this tutorial, we remove it now.\n",
      "    min_docs = 10\n",
      "    c = np.bincount(authors)\n",
      "    validauthors = [authornum for authornum, count in enumerate(c) if count >= min_docs]\n",
      "    documents = [d for d, c in zip(documents, authors) if c in validauthors]\n",
      "    authors = [c for c in authors if c in validauthors]\n",
      "    assert len(documents) == len(authors)\n",
      "    return documents, np.array(authors, dtype='int')\n",
      "\n",
      "\n",
      "def cleanFile(document):\n",
      "    # Removes the stuff pg put in the ebooks\n",
      "    lines = document.split(\"\\n\")\n",
      "    start = 0\n",
      "    end = len(lines)\n",
      "    for i in range(len(lines)):\n",
      "        line = lines[i]\n",
      "        if line.startswith(\"*** START OF THIS PROJECT GUTENBERG\"):\n",
      "            start = i + 1\n",
      "        elif line.startswith(\"*** END OF THIS PROJECT GUTENBERG\"):\n",
      "            end = i - 1\n",
      "    return \"\\n\".join(lines[start:end])\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "documents, classes = get_single_corpus()\n",
      "documents = np.array(documents, dtype=object)  # For fancy indexing using numpy"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 143
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Example document, first 20 lines\n",
      "print(documents[0].split(\"\\n\")[:20])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[u'\\r', u'\\r', u'\\r', u'\\r', u'Produced by David Brannan.  HTML version by Al Haines.\\r', u'\\r', u'\\r', u'\\r', u'\\r', u'\\r', u'\\r', u'\\r', u'\\r', u'\\r', u'The Adventure of the Cardboard Box\\r', u'\\r', u'\\r', u'By\\r', u'\\r', u'Sir Arthur Conan Doyle\\r']\n"
       ]
      }
     ],
     "prompt_number": 144
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Let's graph the most frequent n-grams of the first document\n",
      "counts = count_ngrams(documents[0], 3, normalise=True)\n",
      "top_ngrams = sorted(counts.items(), key=itemgetter(1), reverse=True)[:30]\n",
      "ngrams, ng_counts = zip(*top_ngrams)\n",
      "\n",
      "plt.figure()\n",
      "plt.bar(x_pos, ng_counts, alpha=0.4)\n",
      "plt.xlabel('n-gram (by rank)')\n",
      "plt.ylabel('Normalised Frequency')\n",
      "plt.title('Frequency of the most common n-grams')\n",
      "plt.xlim(0, len(ngrams))\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEXCAYAAABsyHmSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtYlGX+P/D3IHhCdBQFlUFYHVBQxAOGsZpjhiwUZJaK\nmiJaohtaW5a1fbdwy1Yst4NkQat42L5G21ZAEauuomUhlqT9xBTUQQ6eFRUVgeHz+4PL5+vIYYbJ\n4TC+X9fFdc0zz33f87l5cN4+pxmViAiIiIgsYNfSBRARUdvFECEiIosxRIiIyGIMESIishhDhIiI\nLMYQISIiizFEqM05fPgwhg0bhq5duyIhIcGsPnZ2djh27JiVKyO6+zBEWhlPT0907twZTk5OcHJy\nQteuXXHq1KmWLqtVWblyJSZMmIDLly8jNja2znqdToe1a9e2QGV3ji3Mge4ODJFWRqVS4auvvsKV\nK1dw5coVXL58Gb179zZqU11d3ULVtQ6FhYXw9fVtcL1KpWrGaqzDFubQmhkMhpYuwXYItSqenp7y\n3//+t87zKpVK3n//fdFqtdK/f38REUlPTxd/f39Rq9USFBQkBw4cUNrv27dPhg8fLk5OTjJt2jSZ\nNm2a/M///I+IiCQnJ8uYMWPqjH/06FEREamoqJDnnntO+vXrJ66urrJgwQK5fv26iIjs2LFD3Nzc\nZNWqVeLi4iJ9+vSR5ORkZZxr167Js88+Kx4eHtKtWzcZO3asXL9+XcLCwmT16tVGr+nn5ydffvll\nvb+H1NRU8fX1FbVaLTqdTg4dOiQiIuPHj5d27dpJx44dxcnJSfLz8436/fnPf1bWd+nSRRYtWqTM\n78MPPxQvLy9Rq9Xy1FNPGfVbu3at+Pj4SPfu3SUkJEQKCwvrrev48eOiUqkkOTlZ3N3dpUePHvLB\nBx9ITk6O+Pn5iVqtltjYWKV9TU2NvPbaa+Lh4SEuLi4ye/ZsuXTpkoiIXL9+XWbOnCnOzs6iVqtl\n1KhRcvr06QbncLtvv/1W7r33XlGr1eLu7i7r168XEZGysjKZNWuW9OrVSzw8POT111+XmpoaEand\n9kFBQfKnP/1J1Gq1DBgwQHbv3i3r1q0Td3d3cXFxkQ0bNiivERUVJQsXLpTQ0FDp0qWLjBkzRk6e\nPCmLFy8WtVotgwYNktzcXKV9Xl6ejBs3TtRqtQwePFjS0tKMxvrjH/8oDz74oDg5OUlgYKDyN9fQ\n73nDhg3Sr18/6dmzpyxfvrzetjfdbOvs7Kz8zm/+W3r11Vfl0Ucflccff1y6du0qa9eulZycHBk9\nerSo1Wrp06ePxMbGSmVlpTKeSqWSNWvWiFarFScnJ/nLX/4iBQUFMnr0aOnWrZtMmzZNaX/27Fl5\n8MEHRa1WS48ePWTs2LHK79zWMURaGU9PT9m2bVud51UqlUycOFEuXrwoFRUVsm/fPnFxcZGcnByp\nqamRDRs2iKenp1RWVsqNGzekX79+8s4770h1dbV89tln4uDgIH/5y19ExHSIPPPMM/Lwww/LxYsX\n5cqVKxIeHi4vvfSSiNSGiL29vbz66qtSXV0tGRkZ0rlzZykrKxMRkT/+8Y8yfvx4KS0tFYPBID/8\n8IPcuHFDPv30UwkMDFRe7+effxZnZ2epqqqqM9fDhw+Lo6OjbNu2Taqrq2XlypWi1WqVtjqdTtau\nXdvg77C+9SqVSsLDw+XSpUty4sQJ6dWrl2RmZoqIyJdffilarVZ+/fVXMRgM8vrrr0tQUFC9Y998\nc1u4cKHcuHFDtmzZIu3bt5dJkybJ2bNnpaSkRFxcXGTnzp0iUhtOWq1Wjh8/LuXl5TJ58mSZNWuW\niIh8+OGHEh4eLtevX5eamhrZt2+fXL582aw56vV6cXJykk8++USqq6vl/Pnz8vPPP4uIyKxZs2TS\npElSXl4uer1evL29lbGSk5PF3t5e1q9fLzU1NfI///M/4ubmpryBbtmyRZycnOTq1asiUvvG37Nn\nT9m3b59UVFTI/fffLx4eHrJp0yal//jx40VEpLKyUgYMGCB/+9vfpKqqSrZv3y5OTk5y+PBhZSxn\nZ2fZu3evVFdXy8yZMyUyMrLR3/P8+fOloqJC9u/fLx06dFD+M3G7gwcPSpcuXWT37t1SWVkpS5Ys\nEQcHB6MQcXBwkNTUVBGpDfCffvpJ9uzZIwaDQfR6vfj4+Mg777xj9DczadIkuXLlihw8eFDat28v\n48ePl+PHj8ulS5fE19dXCdwXX3xRFixYINXV1VJdXS3fffddg9vO1jBEWhkPDw/p0qWLqNVqUavV\n8sgjj4hI7R/0jh07lHYLFixQQuGmgQMHys6dO2Xnzp3St29fo3VBQUFmhUhNTY04Ojoa/Q/x+++/\nl9/97nciUhsinTp1EoPBoKx3cXFR/jF26tTJaI/opuvXr0v37t2loKBARESee+65OnsDN/31r3+V\nadOmKcs1NTXi5uamvDHrdDr5xz/+UW/fhtarVCrZvXu3sjx16lSJj48XEZE//OEPRm/YBoNBOnfu\nLCdOnKgz9s03t9LSUuU5Z2dn+fTTT5XlRx99VN59910REbn//vvlgw8+UNYdPnxYHBwcpLq6Wtat\nW1dnD7KxOdzqjTfekMmTJ9d5vrq6Wtq3b2/0ZpuYmCg6nU5Eare9l5eXsu7AgQOiUqnkzJkzRvPZ\nv3+/iIjMmTNH5s+fr6xbvXq1+Pr6GvVXq9UiIrJr1y7p3bu3UT3Tp0+XuLg4EakNkSeffFJZl5GR\nIYMGDap3fjd/zyUlJcpz99xzj3zyySf1tl+2bJnMmDFDWb527Zq0b9/eKETGjRtXb9+b3n77beXf\nm0jt38z333+vLI8cOVJWrlypLD/33HPyzDPPiIjIK6+8Ig8//LDy93034TmRVkalUiE1NRUXL17E\nxYsX8fnnnyvr3N3dlceFhYVYtWoVunfvrvwUFxfj5MmTKC0thZubm9G4Hh4eZr3+2bNnce3aNYwc\nOVIZNzQ0FOfOnVPaODs7w87u//50OnfujPLycpw7dw4VFRUYMGBAnXE7duyIqVOnYtOmTRARfPLJ\nJ5g1a1a9NZw8eRL9+vUz+p24u7ujpKTE6LnG1Lf+1nNLN2sGan+XTz/9tDJfZ2dnADB6vdu5uroq\njzt16lRn+ebYJ0+eNPrd9+vXD9XV1Thz5gxmzZqFkJAQREZGws3NDUuXLjU639XYHIuLi9G/f/86\nz587dw5VVVV1XvPWudxeKwD06tWr3voBwMXFRXncsWNHo+Vb25aWlhr9jQK1f3elpaXKfBr6PTXk\n9m129epVAECXLl2UC0+Kiopw8uRJaDQao7Fvbsebbl0PAEeOHMFDDz2EPn36oFu3bnj55Zdx/vx5\nozbmbufnn38eWq0WEydOxIABAxAfH9/ovGwJQ6QNufVNpV+/fnj55ZeVsLl48SLKy8sxbdo09OnT\np84bYGFhofLY0dER165dU5ZvvfqrZ8+e6NSpE/Ly8pRxy8rKcPnyZZP19ezZEx07dkRBQUG966Oi\novDxxx9j27Zt6Ny5MwIDA+tt17dvX6N6RQRFRUV1grEhTT0p3a9fPyQlJRn9Lq9evYrRo0c3aZz6\n9O3bF3q9Xlk+ceIE7O3t4erqCnt7e7zyyis4ePAgvv/+e3z11VfYuHGjWXNwd3fH0aNH6zzfs2dP\nODg41HnN299AraFv374oKiqC3PLB4IWFhWZvt6YoLy9XLjxxd3dHnz59UFxcrKy/fv16nUC4/Xe6\ncOFC+Pr6oqCgAJcuXcLy5ctRU1NjUT1dunTBW2+9haNHjyItLQ1///vfsX37dovGamsYIm3Uk08+\niQ8//BA5OTkQEVy9ehVff/01ysvLERQUBHt7e7z33nuoqqrC559/jr179yp9/f39cfDgQezfvx8V\nFRWIi4tT1tnZ2eHJJ5/EM888g7NnzwKo/R/5li1bTNZkZ2eHuXPn4tlnn8XJkydhMBjwww8/oLKy\nEgBw7733QqVSYcmSJZg9e3aD40ydOhVff/01tm/fjqqqKqxatQodO3ZEUFCQ0kYa+QYDV1fXet9g\nbyW1h3IBAAsWLMAbb7yBvLw8AMClS5fwr3/9y+R8TY0PANOnT8fbb78NvV6P8vJy/PnPf0ZkZCTs\n7OyQlZWFX375BQaDAU5OTnBwcEC7du3MmsPMmTOxbds2/Otf/0J1dTXOnz+P/fv3o127dpg6dSpe\nfvlllJeXo7CwEG+//TYef/zx3zQPcwQGBqJz585YuXIlqqqqkJWVha+++gqRkZFNHqup9Tz22GNI\nT09X/t7i4uJMvl55eTmcnJzQuXNn/Prrr/jggw+a9Pq3Pv7qq69QUFAAEUHXrl3Rrl07ZVvaOoZI\nG3H7/6JGjhyJjz76CLGxsejRowe8vLyU/8U6ODjg888/x/r16+Hs7IxPP/0UkydPVv7ovb298cor\nr+CBBx7AwIEDMXbsWKPx4+PjodVqMXr0aHTr1g3BwcE4cuRIg7Xc6q233oKfnx9GjRoFZ2dnvPTS\nS0b/u5s9ezZ++eWXRt/UvL298c9//hOLFi1Cr1698PXXXyM9PR329vZm1fD000/js88+Q48ePfDM\nM8/U20alUiljTJo0CUuXLkVkZCS6desGPz8//Oc//2lwfHP2dG62mTt3LmbNmoX77rsP/fv3R+fO\nnbF69WoAtXuAU6ZMQbdu3eDr6wudTqcc4jM1B3d3d2RkZGDVqlVwdnbG8OHDceDAAQDA6tWr4ejo\niP79+2Ps2LGYOXMmoqOj68zbnPnc3r6x/u3bt0d6ejq++eYb9OrVC7Gxsdi0aRO8vb0tfm1zngMA\nX19frF69GpGRkejbty+cnJzg4uKCDh06NPjab731Fv73f/8XXbt2xfz58xEZGVlnro29/q1jFhQU\nIDg4GE5OTggKCsJTTz2FcePGNTg3W6KSO/HfA2r1oqOjodFo8Nprr7VoHZs2bcJHH32EXbt2tWgd\nZNvKy8vRvXt3FBQUmH0+kCxj1T2RzMxMDBo0CF5eXg2eaFq8eDG8vLzg7++P3Nxc5fm5c+fC1dUV\nfn5+dfqsXr0aPj4+GDJkCJYuXWq1+m1Ja/i/wrVr1/D+++9j/vz5LV0K2aD09HRcu3YNV69exZIl\nSzB06FAGSDOwWogYDAbExsYiMzMTeXl52Lx5Mw4dOmTUJiMjAwUFBcjPz0dSUhIWLlyorIuOjkZm\nZmadcXfs2IG0tDQcOHAA/+///T8sWbLEWlOwKfXtzjen//znP3BxcUGfPn0wY8aMFquDbFdaWhrc\n3Nzg5uaGo0eP4pNPPmnpku4KVjuc9cMPP2DZsmVKEKxYsQIA8OKLLyptFixYgPHjx2PatGkAgEGD\nBiErK0u5rE+v1yM8PBy//PKL0mfq1KlYsGAB7r//fmuUTURETWC1PZGSkhKja8Y1Gk2dy07NaXO7\n/Px87Nq1C6NHj4ZOp8OPP/54ZwsnIiKz2ZtuYhlzD53cviNkql91dTUuXryI7Oxs7N27F1OnTq33\nI775AXZERJZpygEqq+2JuLm5oaioSFkuKiqqc8PT7W2Ki4tN3pik0WgwefJkAMCoUaNgZ2dX56ai\nm27eC2CLP6+++mqL18C5cX6cn+39NJXVQiQgIAD5+fnQ6/WorKxESkoKIiIijNpEREQo9zZkZ2dD\nrVYbfaxAfSZNmqTcCXrkyBFUVlbW+XgDIiJqHlYLEXt7eyQkJCAkJAS+vr6YNm0afHx8kJiYiMTE\nRABAWFgY+vfvD61Wi5iYGKxZs0bpP336dAQFBeHIkSNwd3dHcnIygNpLf48dOwY/Pz9Mnz5dCSEi\nImp+NnuzoUqlsmjXrK3IysqCTqdr6TKswpbnBnB+bZ2tz6+p750MESIiUjT1vZOfnUVERBZjiBAR\nkcUYIkREZDGGCBERWYwhQkREFmOIEBGRxRgiRERkMYYIERFZjCFCREQWY4gQEZHFrPZ9Iq3Byy8n\nNbrexaUDnn46qpmqISKyPTYdIh4e8xtdX1jYeMgQEVHjeDiLiIgsxhAhIiKLMUSIiMhiDBEiIrIY\nQ4SIiCzGECEiIotZNUQyMzMxaNAgeHl5IT4+vt42ixcvhpeXF/z9/ZGbm6s8P3fuXLi6usLPz6/e\nfqtWrYKdnR0uXLhgldqJiMg0q4WIwWBAbGwsMjMzkZeXh82bN+PQoUNGbTIyMlBQUID8/HwkJSVh\n4cKFyrro6GhkZmbWO3ZRURG2bt0KDw8Pa5VPRERmsFqI5OTkQKvVwtPTEw4ODoiMjERqaqpRm7S0\nNERF1d4xHhgYiLKyMpw6dQoAMHbsWHTv3r3esZ999lmsXLnSWqUTEZGZrBYiJSUlcHd3V5Y1Gg1K\nSkqa3OZ2qamp0Gg0GDp06J0tmIiImsxqH3uiUqnMaiciZve7du0a3njjDWzdurXB/rdKT49THnt7\n6zBwoM6smoiI7hZZWVnIysqyuL/VQsTNzQ1FRUXKclFRETQaTaNtiouL4ebm1uCYR48ehV6vh7+/\nv9J+5MiRyMnJgYuLS5324eFxv3EWRES2TafTQafTKcvLli1rUn+rHc4KCAhAfn4+9Ho9KisrkZKS\ngoiICKM2ERER2LhxIwAgOzsbarUarq6uDY7p5+eH06dP4/jx4zh+/Dg0Gg327dtXb4AQEZH1WS1E\n7O3tkZCQgJCQEPj6+mLatGnw8fFBYmIiEhMTAQBhYWHo378/tFotYmJisGbNGqX/9OnTERQUhCNH\njsDd3R3Jycl1XsPcQ2ZERGQdKmnspEIbplKpkJjY+NQKC5OwfHnjHxdPRHQ3UalUjZ5rvh3vWCci\nIosxRIiIyGIMESIishhDhIiILMYQISIiizFEiIjIYgwRIiKyGEOEiIgsxhAhIiKLMUSIiMhiDBEi\nIrIYQ4SIiCzGECEiIosxRIiIyGIMESIishhDhIiILMYQISIiizFEiIjIYgwRIiKymNVDJDMzE4MG\nDYKXlxfi4+PrbbN48WJ4eXnB398fubm5yvNz586Fq6sr/Pz8jNo///zz8PHxgb+/PyZPnoxLly5Z\ndQ5ERFQ/q4aIwWBAbGwsMjMzkZeXh82bN+PQoUNGbTIyMlBQUID8/HwkJSVh4cKFyrro6GhkZmbW\nGXfixIk4ePAg9u/fD29vb/ztb3+z5jSIiKgBVg2RnJwcaLVaeHp6wsHBAZGRkUhNTTVqk5aWhqio\nKABAYGAgysrKcOrUKQDA2LFj0b179zrjBgcHw87OTulTXFxszWkQEVED7K05eElJCdzd3ZVljUaD\nPXv2mGxTUlKC3r17m/Ua69atw/Tp0+tdl54epzz29tZh4ECd+cUTEd0FsrKykJWVZXF/kyFy/vx5\nODs7WzS4SqUyq52IWNRv+fLlaN++PWbMmFHv+vDwOLPGISK6W+l0Ouh0OmV52bJlTepv8nDW6NGj\nMWXKFGRkZNR5szfFzc0NRUVFynJRURE0Gk2jbYqLi+Hm5mZy7PXr1yMjIwMff/xxk2oiIqI7x2SI\nHD58GE8++SQ2btwIrVaLl156CUeOHDFr8ICAAOTn50Ov16OyshIpKSmIiIgwahMREYGNGzcCALKz\ns6FWq+Hq6trouJmZmXjzzTeRmpqKjh07mlULERHdeSZDxM7ODhMnTsQnn3yCjz76CBs2bMCoUaMw\nbtw4fP/99432tbe3R0JCAkJCQuDr64tp06bBx8cHiYmJSExMBACEhYWhf//+0Gq1iImJwZo1a5T+\n06dPR1BQEI4cOQJ3d3ckJycDABYtWoTy8nIEBwdj+PDh+OMf//hbfgdERGQhlZg4RnXu3Dl8/PHH\n2LhxI1xdXfHEE08gPDwc+/fvx2OPPQa9Xt9MpTaNSqVCYmLjh98KC5OwfPn8ZqqIiKj1U6lUTTp1\nYfLEelBQEB5//HGkpqYanc8ICAjAggULLKuSiIhsgskQOXz4cINXS7344ot3vCAiImo7TJ4TmThx\nIsrKypTlCxcuICQkxKpFERFR22AyRM6ePQu1Wq0s9+jRA6dPn7ZqUURE1DaYDJF27dqhsLBQWdbr\n9cpHjhAR0d3N5DmR5cuXY+zYsbjvvvsAALt27UJSUpLVCyMiotbPZIj84Q9/wE8//YTs7GyoVCq8\n88476NmzZ3PURkRErZxZH8BYWVmJHj16oLq6Gnl5eQCg7JkQEdHdy2SILF26FCkpKfD19UW7du2U\n5xkiRERkMkS++OILHD58GB06dGiOeoiIqA0xGSIDBgxAZWWlTYfIu+9uwJkzNxpt4+LSAU8/HdVM\nFRERtQ0mQ6RTp04YNmwYJkyYoASJSqXCe++9Z/XimsuZMzfg4dH4Z2gVFvKKNCKi25kMkYiICERE\nRCgffSIiZn9pFBER2TaTITJnzhxcu3YNJ06cwKBBg5qjJiIiaiNM3nqelpaG4cOH4w9/+AMAIDc3\nt84XSxER0d3J5J5IXFwc9uzZg/HjxwMAhg8fjmPHjlm9sNaKJ+GJiP6PyRBxcHAw+gBGAHf1Z2fx\nJDwR0f8xmQaDBw/Gxx9/jOrqauTn52PRokUICgpqjtqIiKiVMxkiq1evxsGDB9GhQwdMnz4dXbt2\nxTvvvNMctRERUStnMkQcHR3xxhtv4Mcff8SPP/6I5cuXo2PHjmYNnpmZiUGDBsHLywvx8fH1tlm8\neDG8vLzg7++P3Nxc5fm5c+fC1dUVfn5+Ru0vXLiA4OBgeHt71/nCLCIial4mQ2T8+PF1fu6//36T\nAxsMBsTGxiIzMxN5eXnYvHkzDh06ZNQmIyMDBQUFyM/PR1JSEhYuXKisi46ORmZmZp1xV6xYgeDg\nYBw5cgQTJkzAihUrzJknERFZgckT62+++abyuKKiAv/+979hb2/6w39zcnKg1Wrh6ekJAIiMjERq\naip8fHyUNmlpaYiKqr2KKTAwEGVlZTh16hR69+6NsWPHQq/X1xk3LS0NO3fuBABERUVBp9MxSIiI\nWojJNAgICDBaHjNmDEaNGmVy4JKSEri7uyvLGo0Ge/bsMdmmpKQEvXv3bnDc06dPw9XVFQDg6ura\n6Ff1pqfHKY+9vXUYOFBnsm4iortJVlYWsrKyLO5vMkQuXLigPK6pqcGPP/6Iy5cvmxzY3I9GERGL\n+t1s21j78PA4s8ciIrob6XQ66HQ6ZXnZsmVN6m8yREaMGKG8Udvb28PT0xNr1641ObCbmxuKioqU\n5aKiImg0mkbbFBcXw83NrdFxXV1dlUNeJ0+ehIuLi8laiIjIOkyeWNfr9Th+/DiOHz+O/Px8bN26\nFWPGjDE5cEBAAPLz86HX61FZWYmUlJQ6H5cSERGBjRs3AgCys7OhVquVQ1UNiYiIwIYNGwAAGzZs\nwKRJk0zWQkRE1mFyT+Tf//53nUNGNw9BqVQqTJ48uf6B7e2RkJCAkJAQGAwGzJs3Dz4+PkhMTAQA\nxMTEICwsDBkZGdBqtXB0dERycrLSf/r06di5cyfOnz8Pd3d3/PWvf0V0dDRefPFFTJ06FWvXroWn\npyc+/fRTiydPRES/jckQWbduHb7//nvcf//9EBFkZWXh3nvvVQ4jNRQiABAaGorQ0FCj52JiYoyW\nExIS6u27efPmep/v0aMHtm3bZqpsIiJqBiZDpLKyEnl5eejTpw8A4OTJk4iKijLaayAioruTyXMi\nRUVFRpfcurq64sSJE1YtioiI2gaTeyIPPPAAQkJCMGPGDIgIUlJSEBwc3By1ERFRK2cyRFavXo0v\nvvgC3377LYDacxqPPPKI1QsjIqLWz2SIqFQqjBgxAk5OTggODsa1a9dw5coVODk5NUd9RETUipk8\nJ5KUlIQpU6ZgwYIFAGpvCOS9GUREBJgRIu+//z6+++47dO3aFQDg7e2NM2fOWL0wIiJq/UyGSIcO\nHdChQwdlubq6ukmfb0VERLbLZIiMGzcOy5cvx7Vr17B161ZMmTIF4eHhzVEbERG1ciZDJD4+Hr16\n9YKfnx8SExMRFhaG119/vTlqIyKiVq7Rq7Oqq6sxZMgQ/Prrr5g/f35z1URERG1Eo3si9vb2GDhw\nIAoLC5urHiIiakPM+lKqwYMH45577oGjoyOA2ntH0tLSrF4cERG1bg2GSHV1Nezt7fHaa69BpVIZ\nfQMhr84iIiKgkRC55557sG/fPuh0OixatAirV69uzrqIiKgNaPCcyK17Ht99912zFENERG2LyUt8\niYiIGtLg4axff/0Vfn5+AICjR48qj4HacyIHDhywfnVERNSqNRgihw4das46iIioDWrwcJanp2ej\nP+bIzMzEoEGD4OXlhfj4+HrbLF68GF5eXvD390dubq7Jvjk5ObjnnnswfPhwjBo1Cnv37jVzqkRE\ndKdZ7ZyIwWBAbGwsMjMzkZeXh82bN9fZu8nIyEBBQQHy8/ORlJSEhQsXmuz7wgsv4LXXXkNubi7+\n+te/4oUXXrDWFIiIyASrhUhOTg60Wi08PT3h4OCAyMhIpKamGrVJS0tDVFQUACAwMBBlZWU4depU\no3379OmDS5cuAQDKysrg5uZmrSkQEZEJJu9Yt1RJSQnc3d2VZY1Ggz179phsU1JSgtLS0gb7rlix\nAmPGjMGSJUtQU1ODH374ocEa0tPjlMfe3joMHKj7jbMiIrItWVlZyMrKsrh/gyFy69VYtzPn6ixz\n72q/9X4Uc8ybNw/vvfceHnnkEfzrX//C3LlzsXXr1nrbhofHNWlsIqK7jU6ng06nU5aXLVvWpP4N\nhkh6ejoAYM2aNQCAWbNmQUTw8ccfmzWwm5sbioqKlOWioiJoNJpG2xQXF0Oj0aCqqqrBvjk5Odi2\nbRsA4LHHHsMTTzxhVj1ERHTnmbw6a8uWLVi5ciX8/PwwdOhQxMfHY8uWLSYHDggIQH5+PvR6PSor\nK5GSkoKIiAijNhEREdi4cSMAIDs7G2q1Gq6uro321Wq12LlzJwBg+/bt8Pb2tnjyRET025g8JyIi\n+O677zBmzBgAwO7du806BGVvb4+EhASEhITAYDBg3rx58PHxQWJiIgAgJiYGYWFhyMjIgFarhaOj\nI5KTkxvtCwBJSUl46qmncOPGDXTq1AlJSUkWT56IiH4bkyGybt06REdHK1dEqdVq5c3elNDQUISG\nhho9FxN20ztHAAAYI0lEQVQTY7SckJBgdl+gdg/n9hP0RETUMkyGyMiRI3HgwAFcunQJIgK1Wt0c\ndRERURtg8j6RU6dOYd68eZg2bRrUajXy8vKwdu3a5qiNiIhaOZMhMmfOHEycOBGlpaUAAC8vL7z9\n9ttWL4yIiFo/kyFy7tw5TJs2De3atQMAODg4wN7eavcoEhFRG2IyRLp06YLz588ry9nZ2ejWrZtV\niyIiorbB5C7FqlWrEB4ejmPHjiEoKAhnz57FZ5991hy1ERFRK2fW1Vk7d+7E4cOHAQADBw6Eg4OD\n1QsjIqLWz+ThrE8//RTXr1/HkCFD8MUXX2DatGnYt29fc9RGREStnMkQee2119C1a1d89913+O9/\n/4u5c+diwYIFzVEbERG1ciZD5OZVWV999RWefPJJPPTQQ6iqqrJ6YURE1PqZDBE3NzfMnz8fKSkp\nePDBB1FRUYGamprmqI2IiFo5s86JhISEYMuWLVCr1bh48SLefPPN5qiNiIhauQavzrp8+TK6du2K\nGzduYPz48QCACxcuoEOHDggICGi2AomIqPVqMESmT5+Or7/+GiNGjKjzLYUqlQrHjh2zenFERNS6\nNRgiX3/9NQBAr9c3Vy1ERNTGNBgipu4FGTFixB0vhoiI2pYGQ+TZZ5+tcxjrVjt27LBKQbbk3Xc3\n4MyZG422cXHpgKefjmqmioiI7qwGQyQrK6sZy7BNZ87cgIfH/EbbFBby632JqO0y6zPdf/nlFxw6\ndAgVFRXKc7Nnz7ZaUURE1DaYvE8kLi4OixcvRmxsLHbs2IEXXngBaWlpZg2emZmJQYMGwcvLC/Hx\n8fW2Wbx4Mby8vODv74/c3Fyz+q5evRo+Pj4YMmQIli5dalYtRER055ncE/nss8+wf/9+jBgxAsnJ\nyTh9+jRmzpxpcmCDwYDY2Fhs27YNbm5uGDVqFCIiIuDj46O0ycjIQEFBAfLz87Fnzx4sXLgQ2dnZ\njfbdsWMH0tLScODAATg4OODs2bO/7TdAREQWM7kn0qlTJ7Rr1w729va4dOkSXFxcUFRUZHLgnJwc\naLVaeHp6wsHBAZGRkUhNTTVqk5aWhqio2pPKgYGBKCsrw6lTpxrt+8EHH+Cll15SPo6+V69eTZ40\nERHdGSb3REaNGoWLFy/iySefREBAABwdHREUFGRy4JKSEri7uyvLGo0Ge/bsMdmmpKQEpaWlDfbN\nz8/Hrl278Oc//xkdO3bEW2+91eAd9Onpccpjb28dBg7UmaybiOhukpWV9ZsupDIZImvWrAEALFiw\nACEhIbhy5QqGDh1qcuDGLg++lYiY1e6m6upqXLx4EdnZ2di7dy+mTp3a4N3z4eFxTRqbiOhuo9Pp\noNPplOVly5Y1qb9ZV2ft378fer0eBoMBIoKCggJMnjy50T5ubm5Gh72Kioqg0WgabVNcXAyNRoOq\nqqoG+2o0GuW1R40aBTs7O5w/fx7Ozs7mTKXVauo9Jaba8/4TImoOJkMkOjoav/zyCwYPHgw7u/87\nhWIqRAICApCfnw+9Xo++ffsiJSUFmzdvNmoTERGBhIQEREZGIjs7G2q1Gq6urnB2dm6w76RJk7B9\n+3aMGzcOR44cQWVlZZsPEKDp95SYan9rW970SETWYjJE9uzZg4MHD5p9eEoZ2N4eCQkJCAkJgcFg\nwLx58+Dj44PExEQAQExMDMLCwpCRkQGtVgtHR0ckJyc32hcA5s6di7lz58LPzw/t27fHxo0bmzrn\nu05TAqopgcNwIiKzTqzn5eVh8ODBTR48NDQUoaGhRs/FxMQYLSckJJjdFwAcHBywadOmJtdC5mlK\n4PCOfCIy63DWvffei969e6NDhw4Aak+aHzhwwOrFke3gXguRbTIZIvPmzcM///lPDBkyxOicCFFT\ncK+FyDaZDBEXFxdEREQ0Ry1ERNTGmAyRYcOGYcaMGQgPD0f79u0B1B7OMnV1FpGleOiLqO0wGSIV\nFRXo0KEDtmzZYvQ8Q4SshYe+iNqORkPEYDCgR48eWLVqVXPVQ0REbUijIdKuXTvs3r0bItLk+0SI\nmgMPfRG1LLPOiTz88MOYMmUKOnfuDIDnRKj14KEvopZl1jmRHj16YPv27UbPM0SIiMhkiKxfv74Z\nyiAiorbI5N2DRUVFeOSRR9CrVy/06tULjz76KIqLi5ujNiIiauVMhkh0dDQiIiJQWlqK0tJShIeH\nIzo6ujlqIyKiVs5kiJw9exbR0dFwcHCAg4MD5syZgzNnzjRHbURE1MqZPCfi7OyMTZs2YcaMGRAR\nfPLJJ+jZs2dz1EZ0x/GSYKI7y2SIrFu3DosWLcKzzz4LAAgKClK+94OoreF3qxDdWSZDxNPTE+np\n6c1RC1Grwu9WITKtwRBp6Mvab965/sorr1inIiIiajMaDBFHR8c6H3Vy9epVrF27FufOnWOIEFmI\nh77IljQYIkuWLFEeX758Ge+99x6Sk5MRGRmJ5557rlmKI7JFPPRFtqTRcyLnz5/H22+/jY8//hiz\nZ8/Gvn370L17d7MHz8zMxDPPPAODwYAnnngCS5curdNm8eLF+Oabb9C5c2esX78ew4cPN6vvqlWr\n8Pzzz+PcuXPo0aOH2TURtSXca6HWrtE9kS+++ALz58/HgQMH4OTk1KSBDQYDYmNjsW3bNri5uWHU\nqFGIiIiAj4+P0iYjIwMFBQXIz8/Hnj17sHDhQmRnZ5vsW1RUhK1bt8LDw8PCaRO1Dda6msyS9kT1\naTBE/v73v6N9+/Z4/fXX8frrrxutU6lUuHz5cqMD5+TkQKvVwtPTEwAQGRmJ1NRUoxBJS0tDVFTt\nH2hgYCDKyspw6tQpHD9+vNG+zz77LFauXImHH364yRMmslVNPUzGy53pTmgwRGpqan7TwCUlJXB3\nd1eWNRoN9uzZY7JNSUkJSktLG+ybmpoKjUaDoUOHmqwhPT1OeeztrcPAgToLZ0N0d+PlzrYrKysL\nWVlZFvc3eZ+Ipcz9EisRMXvM69ev44033sDWrVvN6h8eHmf22ETU/LjX0vJ0Oh10Op2y3NDtHQ2x\nWoi4ubmhqKhIWS4qKoJGo2m0TXFxMTQaDaqqqurte/ToUej1evj7+yvtR44ciZycHLi4uFhrKkRk\nJdxrafusFiIBAQHIz8+HXq9H3759kZKSgs2bNxu1iYiIQEJCAiIjI5GdnQ21Wg1XV1c4OzvX29fH\nxwenT59W+v/ud7/DTz/9xKuziO4C3GtpnawWIvb29khISEBISAgMBgPmzZsHHx8fJCYmAgBiYmIQ\nFhaGjIwMaLVaODo6Kp/J1VDf2/F734nuHtxraZ2sFiIAEBoaitDQUKPnYmJijJYTEhLM7nu7Y8eO\n/bYCicgm3enLnbmH0zCrhggRUUu405c7cw+nYQwRIqIm4D0zxhgiRERNYK17Ztpq4DBEiIhagbZ6\n4QBDhIiojWlNey0MESKiNqY1HSZjiBAR2TBrHyazs7gnERHd9RgiRERkMYYIERFZjCFCREQWY4gQ\nEZHFGCJERGQxhggREVmMIUJERBZjiBARkcUYIkREZDGGCBERWYwhQkREFrN6iGRmZmLQoEHw8vJC\nfHx8vW0WL14MLy8v+Pv7Izc312Tf559/Hj4+PvD398fkyZNx6dIla0+DiIjqYdUQMRgMiI2NRWZm\nJvLy8rB582YcOnTIqE1GRgYKCgqQn5+PpKQkLFy40GTfiRMn4uDBg9i/fz+8vb3xt7/9zZrTICKi\nBlg1RHJycqDVauHp6QkHBwdERkYiNTXVqE1aWhqiomo/xz4wMBBlZWU4depUo32Dg4NhZ2en9Cku\nLrbmNIiIqAFW/T6RkpISuLu7K8sajQZ79uwx2aakpASlpaUm+wLAunXrMH369HpfPz09Tnns7a3D\nwIE6C2dCRGSbDh/Owt696YiLK7Wov1VDRKVSmdVORCwaf/ny5Wjfvj1mzJhR7/rw8DiLxiUiulsM\nHKhDx45HEBdX+8VVy5Yta1J/q4aIm5sbioqKlOWioiJoNJpG2xQXF0Oj0aCqqqrRvuvXr0dGRgb+\n+9//WnEGRETUGKueEwkICEB+fj70ej0qKyuRkpKCiIgIozYRERHYuHEjACA7OxtqtRqurq6N9s3M\nzMSbb76J1NRUdOzY0ZpTICKiRlh1T8Te3h4JCQkICQmBwWDAvHnz4OPjg8TERABATEwMwsLCkJGR\nAa1WC0dHRyQnJzfaFwAWLVqEyspKBAcHAwDuvfderFmzxppTISKielg1RAAgNDQUoaGhRs/FxMQY\nLSckJJjdFwDy8/PvXIFERGQx3rFOREQWY4gQEZHFGCJERGQxhggREVmMIUJERBZjiBARkcUYIkRE\nZDGGCBERWYwhQkREFmOIEBGRxRgiRERkMYYIERFZjCFCREQWY4gQEZHFGCJERGQxhggREVmMIUJE\nRBZjiBARkcWsGiKZmZkYNGgQvLy8EB8fX2+bxYsXw8vLC/7+/sjNzTXZ98KFCwgODoa3tzcmTpyI\nsrIya06h1SosPNzSJViNLc8N4PzaOlufX1NZLUQMBgNiY2ORmZmJvLw8bN68GYcOHTJqk5GRgYKC\nAuTn5yMpKQkLFy402XfFihUIDg7GkSNHMGHCBKxYscJaU2jVTpw40tIlWI0tzw3g/No6W59fU1kt\nRHJycqDVauHp6QkHBwdERkYiNTXVqE1aWhqioqIAAIGBgSgrK8OpU6ca7Xtrn6ioKHz55ZfWmgIR\nEZlgtRApKSmBu7u7sqzRaFBSUmJWm9LS0gb7nj59Gq6urgAAV1dXnD592lpTICIiU8RKPvvsM3ni\niSeU5U2bNklsbKxRm4ceeki+++47ZXnChAny448/1um7ceNGWbRokYiIqNVqozG6d+9e7+sD4A9/\n+MMf/ljw0xT2sBI3NzcUFRUpy0VFRdBoNI22KS4uhkajQVVVVZ3n3dzcANTufZw6dQq9e/fGyZMn\n4eLiUu/r1+YIERFZk9UOZwUEBCA/Px96vR6VlZVISUlBRESEUZuIiAhs3LgRAJCdnQ21Wg1XV9dG\n+0ZERGDDhg0AgA0bNmDSpEnWmgIREZlgtT0Re3t7JCQkICQkBAaDAfPmzYOPjw8SExMBADExMQgL\nC0NGRga0Wi0cHR2RnJzcaF8AePHFFzF16lSsXbsWnp6e+PTTT601BSIiMqVJB7/agG+++UYGDhwo\nWq1WVqxY0dLl3HEeHh7i5+cnw4YNk1GjRrV0Ob9ZdHS0uLi4yJAhQ5Tnzp8/Lw888IB4eXlJcHCw\nXLx4sQUr/G3qm9+rr74qbm5uMmzYMBk2bJh88803LVih5U6cOCE6nU58fX1l8ODB8u6774qI7Wy/\nhuZnK9vv+vXrcs8994i/v7/4+PjIiy++KCJN3342FSLV1dUyYMAAOX78uFRWVoq/v7/k5eW1dFl3\nlKenp5w/f76ly7hjdu3aJfv27TN6k33++eclPj5eRERWrFghS5cubanyfrP65hcXFyerVq1qwaru\njJMnT0pubq6IiFy5ckW8vb0lLy/PZrZfQ/Ozle0nInL16lUREamqqpLAwED59ttvm7z9bOpjT8y5\nN8UWiA1dNDB27Fh0797d6DlbuheovvkBtrENe/fujWHDhgEAunTpAh8fH5SUlNjM9mtofoBtbD8A\n6Ny5MwCgsrISBoMB3bt3b/L2s6kQMefelLZOpVLhgQceQEBAAD766KOWLscq7oZ7gVavXg1/f3/M\nmzfPJj66R6/XIzc3F4GBgTa5/W7Ob/To0QBsZ/vV1NRg2LBhcHV1xfjx4zF48OAmbz+bChGVStXS\nJVjd7t27kZubi2+++Qbvv/8+vv3225YuyapUKpXNbdeFCxfi+PHj+Pnnn9GnTx8899xzLV3Sb1Je\nXo5HH30U7777LpycnIzW2cL2Ky8vx2OPPYZ3330XXbp0santZ2dnh59//hnFxcXYtWsXduzYYbTe\nnO1nUyFizr0pbV2fPn0AAL169cIjjzyCnJycFq7ozrt5LxCARu8FaqtcXFyUf5xPPPFEm96GVVVV\nePTRRzFr1izlcntb2n435/f4448r87Ol7XdTt27d8OCDD+Knn35q8vazqRAx596UtuzatWu4cuUK\nAODq1avYsmUL/Pz8WriqO8/W7wU6efKk8viLL75os9tQRDBv3jz4+vrimWeeUZ63le3X0PxsZfud\nO3dOORR3/fp1bN26FcOHD2/69rPmmf+WkJGRId7e3jJgwAB54403WrqcO+rYsWPi7+8v/v7+Mnjw\nYJuYX2RkpPTp00ccHBxEo9HIunXr5Pz58zJhwoQ2f4moSN35rV27VmbNmiV+fn4ydOhQefjhh+XU\nqVMtXaZFvv32W1GpVOLv7290uautbL/65peRkWEz2+/AgQMyfPhw8ff3Fz8/P1m5cqWISJO3n0rE\nRi4zICKiZmdTh7OIiKh5MUSIiMhiDBEiIrIYQ4SIiCzGECFqRkuWLMHOnTsBAJ6enrhw4UKL1dKl\nS5c6z50+fRphYWEtUA21VQwRogYYDIY7Ot6VK1ewa9cujBs3DkDt3cB34uJIS+us705kV1dXdO/e\nHfv27futZdFdgiFCNkOv18PHxwfz58/HkCFDEBISgoqKinrbvvbaaxg0aBDGjh2LGTNmYNWqVQAA\nnU6HP/3pTxg1ahTeffddfPXVVxg9ejRGjBiB4OBgnDlzBgAQFxeHqKgo3HffffD09MTnn3+OJUuW\nYOjQoQgNDUV1dXWd10xNTcUDDzxg9NzKlSsxdOhQBAYG4ujRo7hy5Qr69++v9L98+TL69+9fJyjm\nzJmDBQsWYPTo0Vi6dCn27t2LoKAgjBgxAr///e9x5MgRAMD69esxefJkhIaGwtvbG0uXLq1T17lz\n5xAUFIRvvvkGQO3Ngps3b27Kr57uZla/o4WomRw/flzs7e1l//79IiIydepU+ec//1mnXU5Ojgwb\nNkxu3LghV65cES8vL+WjvXU6nTz11FNK21tvtProo4/kueeeE5Ha75QYO3asVFdXy/79+6VTp06S\nmZkpIiKPPPKIfPnll3Ved8GCBfLvf/9bWfb09FRuGN24caM89NBDIlL7HSQ3+ycmJsqSJUvqjDVn\nzhwJDw+XmpoaERG5fPmyVFdXi4jI1q1b5dFHHxURkeTkZOnfv79cvnxZKioqxMPDQ4qLi0VEpEuX\nLnL69GkJDAyUbdu2KWMfO3ZM7rnnnoZ+zURGrPbNhkQt4Xe/+x2GDh0KABg5ciT0en2dNrt378ak\nSZPQvn17tG/fHuHh4Ubrp02bpjwuKirC1KlTcerUKVRWVqJ///4Aag8FhYaGol27dhgyZAhqamoQ\nEhICAPDz86v3dQsLC5XPPrtp+vTpAIDIyEj86U9/AgA88cQTWLlyJR5++GGsX78e//jHP+qd65Qp\nU5RDUmVlZZg9ezYKCgqgUqmM9oQmTJigfDCir68vCgsL4ebmhsrKSkyYMAFr1qzB2LFjlfZ9+vSp\nt36i+vBwFtmUDh06KI/btWuH6upqFBcXY9iwYRg+fDgSExPrnIuQ285LODo6Ko8XLVqExYsX48CB\nA0hMTMT169eVde3btwdQ+0moDg4OyvN2dnb1Hs4Caj96uyE3AyEoKAh6vR5ZWVkwGAzw9fWtt/3N\n74IAgL/85S+YMGECfvnlF6SnpxvVWd/vBAAcHBwQEBCAzMxMo3FFpM1/8i41H4YI2TyNRoOff/4Z\nubm5iImJwe9//3ukp6fjxo0bKC8vx9dff23U/tZQuXz5Mvr27Qug9vxCfW3M5eHhoXw66s0xUlJS\nAAApKSkICgpS1s2ePRszZ87E3LlzzRr71jqTk5PN6qNSqbBu3Tr8+uuvWLlypfL8yZMn4eHhYdYY\nRAwRsim3/w+6vv9RBwQEICIiAkOHDkVYWBj8/PzQrVu3evvExcVhypQpCAgIQK9evZR1t3/Pgjmv\nO2bMGPz4449GbS5evAh/f3+sXr0ab7/9trJuxowZuHjxonK4y9RcX3jhBbz00ksYMWIEDAZDg3Xe\n3l+lUmHz5s3Yvn07PvzwQwC13xB63333Nfi6RLfiBzDSXenq1atwdHTEtWvXMG7cOHz00UfKV6Fa\nS3l5OcaPH4+9e/eabPvZZ58hPT1d+Uju5jRz5kwsWbIEw4cPb/bXpraHJ9bprjR//nzk5eWhoqIC\nc+bMsXqAALU3940fPx47duzA+PHjG2y3aNEi/Oc//0FGRobVa7rdmTNnUFZWxgAhs3FPhIiILMZz\nIkREZDGGCBERWYwhQkREFmOIEBGRxRgiRERkMYYIERFZ7P8Den7Z1nYLasAAAAAASUVORK5CYII=\n"
      }
     ],
     "prompt_number": 145
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "# Print stats about the authors\n",
      "print(\"\\n\".join(\"Author {} appears {} times\".format(*r) for r in enumerate(np.bincount(classes))))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Author 0 appears 51 times\n",
        "Author 1 appears 44 times\n",
        "Author 2 appears 22 times\n",
        "Author 3 appears 10 times\n",
        "Author 4 appears 11 times\n",
        "Author 5 appears 29 times\n",
        "Author 6 appears 10 times\n"
       ]
      }
     ],
     "prompt_number": 146
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "That's all the setup that is needed for this tutorial."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "1) Source Code Author Profiles (SCAP)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "SCAP was not the first LNG method to be created, but it is the easiest to explain.\n",
      "In SCAP, documents are *profiled* by taking the *L* most frequent n-grams that occur in the document.\n",
      "An author is profiled in a similar way, except the *L* most frequent n-grams from all known documents from that author are used.\n",
      "Two profiles are compared by finding the number of n-grams each profile has in common, which can be normalised by dividing by L.\n",
      "The more n-grams they have in common, the more similar the profiles are.\n",
      "When trying to determine which author wrote a particular document, the most similar profile.\n",
      "To make this calculation consistent with other methods that will be introduced later, we use a distance instead of similarity.\n",
      "The distance is just `1 - s`, where s is the similarity.\n",
      "In this format, the lower the distance, the more likely the author, and the author with the lowest distance to a document is the most likely author of that document.\n",
      "\n",
      "For two profiles $P_1$ and $P_2$ (which can be either author or document profiles) the distance is given below.\n",
      "An n-gram $x$ is in $P_i$ if it is in the profile, and the intersection of two profiles is the set of n-grams that are in both.\n",
      "\n",
      "*As a similarity:*\n",
      "$$s_{SCAP}(P_1, P_2) = \\frac{1}{L}|P_1 \\cap P_2|$$\n",
      "\n",
      "*As a distance:*\n",
      "$$d_{SCAP}(P_1, P_2) = min(0, 1 - s_{SCAP})$$ \n",
      "\n",
      "\n",
      "The class below implements SCAP using the base classes from scikit-learn (imported as sklearn).\n",
      "This will be important later, but for now, the important thing to acknoweldge is the `fit` and `predict` functions.\n",
      "The `fit` function takes a set of training data and creates profiles for each author in that data.\n",
      "The `predict` function predicts which of those authors is most likely to have written each document.\n",
      "The class also has utility functions for creating profiles and comparing profiles.\n",
      "\n",
      "The reference for this paper is below:"
     ]
    },
    {
     "cell_type": "raw",
     "metadata": {},
     "source": [
      "@inproceedings{frantzeskou2006scap,\n",
      "\ttitle = {Effective identification of source code authors using byte-level information},\n",
      "\tbooktitle = {Proceedings of the 28th international conference on Software engineering},\n",
      "\tauthor = {Frantzeskou, Georgia and Stamatatos, Efstathios and Gritzalis, Stefanos and Katsikas, Sokratis},\n",
      "\tyear = {2006},\n",
      "\tpages = {893\u2013896},\n",
      "}"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from operator import itemgetter\n",
      "from sklearn.base import BaseEstimator, ClassifierMixin  # For using scikit-learn's test framework. Ignore for now\n",
      "\n",
      "class SCAP(BaseEstimator, ClassifierMixin):\n",
      "    def __init__(self, n, L):\n",
      "        self.n = n\n",
      "        self.L = L\n",
      "        self.author_profiles = None  # Will be trained by fit()\n",
      "\n",
      "    def create_profile(self, documents):\n",
      "        # Creates a profile of a document or list of documents.\n",
      "        if isinstance(documents, basestring):\n",
      "            # documents can be either a list of documents, or a single document.\n",
      "            # if it's a single document, convert to a list\n",
      "            documents = [documents,]\n",
      "        # profile each document independently\n",
      "        profiles = (count_ngrams(document, self.n, normalise=False)\n",
      "                    for document in documents)\n",
      "        # Merge the profiles\n",
      "        main_profile = defaultdict(float)\n",
      "        for profile in profiles:\n",
      "            for ngram in profile:\n",
      "                main_profile[ngram] += profile[ngram]\n",
      "        # Normalise the profile\n",
      "        num_ngrams = sum(main_profile.values())\n",
      "        for ngram in main_profile:\n",
      "            main_profile[ngram] /= num_ngrams\n",
      "        # Return the profile with only the top L n-grams\n",
      "        return self.top_L(main_profile)\n",
      "\n",
      "    def top_L(self, profile):\n",
      "        # Returns the profile with only the top L most frequent n-grams\n",
      "        if self.L >= len(profile):\n",
      "            return profile\n",
      "        threshold = sorted(profile.values())[-self.L]\n",
      "        return {ngram: profile[ngram]\n",
      "                for ngram in profile\n",
      "                if profile[ngram] >= threshold}\n",
      "\n",
      "    def compare_profiles(self, profile1, profile2):\n",
      "        # Number of n-grams in both profiles, divided by L\n",
      "        similarity =  len(set(profile1.keys()) & set(profile2.keys())) / float(self.L)\n",
      "        # Slight edge case here, similarity could be higher than 1.\n",
      "        # Just make it equal to 1 if it is over.\n",
      "        similarity = min(similarity, 1.0)\n",
      "        distance = 1 - similarity\n",
      "        return distance\n",
      "\n",
      "    def fit(self, documents, classes):\n",
      "        # Fits the current model to the training data provided\n",
      "        # Separate documents into the sets written by each author\n",
      "        author_documents = ((author, [documents[i] \n",
      "                                      for i in range(len(documents))\n",
      "                                      if classes[i] == author])\n",
      "                            for author in set(classes))\n",
      "        # Profile each of the authors independently\n",
      "        self.author_profiles = {author:self.create_profile(cur_docs)\n",
      "                                for author, cur_docs in author_documents}\n",
      "    \n",
      "    def predict(self, documents):\n",
      "        # Predict which of the authors wrote each of the documents\n",
      "        predictions = np.array([self.predict_single(document) for document in documents])\n",
      "        return predictions\n",
      "\n",
      "    def predict_single(self, document):\n",
      "        # Predicts the author of a single document\n",
      "        # Profile of current document\n",
      "        profile = self.create_profile(document)\n",
      "        # Distance from document to each author\n",
      "        distances = [(author, self.compare_profiles(profile, self.author_profiles[author]))\n",
      "                      for author in self.author_profiles]\n",
      "        # Get the nearest pair, and the author from that pair\n",
      "        prediction = sorted(distances, key=itemgetter(1))[0][0]\n",
      "        return prediction"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 147
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "OK, let's see it in action. I'll point out here that this method of testing should not be used in the real world. The fundamental mistake I'll be making is to use training data to test the method, which should never be done in practice. Anyway, let's do it now:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "model = SCAP(n=4, L=2000)\n",
      "# Training the model\n",
      "model.fit(documents, classes)\n",
      "# Predict the author of each of the documents\n",
      "y_pred = model.predict(documents)  # Notice how I am predicting the data I used for training? Don't do that.\n",
      "# Compute the accuracy of these predictions.\n",
      "print(\"Accuracy is {:.1f}%\".format(100. * np.mean(classes == y_pred)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Accuracy is 80.8%\n"
       ]
      }
     ],
     "prompt_number": 148
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Not a bad result, but as I said, don't use training data for testing.\n",
      "Here is a function that uses scikit-learn to perform the testing, which automatically splits the data into training and testing data before computing the accuracy, through a process called cross fold validation.\n",
      "More information on this is available in the [scikit-learn documentation](http://scikit-learn.org/stable/modules/cross_validation.html), but I won't go into the details here.\n",
      "In a nutshell, what cross validation does is split the data into two sets, training and testing.\n",
      "It does this multiple times with different splits of the data (depending on the exact method used)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# First, setup the cross fold validation object\n",
      "from sklearn import cross_validation\n",
      "cv = cross_validation.KFold(len(documents), n_folds=5, random_state=0)\n",
      "# Now, compute the score using it\n",
      "scores = cross_validation.cross_val_score(model, documents, classes, cv=cv)\n",
      "print(\"Accuracy is {:.1f}%\".format(100. * np.mean(scores)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Accuracy is 54.2%\n"
       ]
      }
     ],
     "prompt_number": 149
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "A fair bit lower, but not too bad (the chance rate is around 28%). This is a better representation of what we can expect the algorithm to do if applied to real data with these parameters."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "2) Common n-grams (CNG)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "CNG was the first LNG method proposed in the literature.\n",
      "The method is almost entirely the same as SCAP (should that be the other way around?), except for the method used to compute the distance.\n",
      "This distance computation uses the frequencies of the n-grams, where SCAP simply used the fact they were in the top *L* n-grams.\n",
      "\n",
      "$$ d_{CNG}(P_1, P_2) = \\displaystyle\\sum_{x \\in X_{P_1} \\cup X_{P_2}} \\Biggl(\\frac{2\\cdot(P_1(x) - P_2(x))}{P_1(x) + P_2(x)}\\Biggr)^2 $$\n",
      "\n",
      "Note that $P_i(x)$ is the frequency of n-gram $x$ in profile $P_i$, which is zero if the n-gram is not in the profile.\n",
      "(This is important, as it is zero if it is not in the top L, even if it does have a frequency in the original document.)\n",
      "\n",
      "To simplify the coding, we simply derive our new class from the SCAP class defined above. The reference for this method is below."
     ]
    },
    {
     "cell_type": "raw",
     "metadata": {},
     "source": [
      "@inproceedings{kevsel2003cng,\n",
      "\ttitle = {N-gram-based author profiles for authorship attribution},\n",
      "\tvolume = {3},\n",
      "\tbooktitle = {Proceedings of the Conference Pacific Association for Computational Linguistics, {PACLING}},\n",
      "\tauthor = {Ke{\\textbackslash}vselj, Vlado and Peng, Fuchun and Cercone, Nick and Thomas, Calvin},\n",
      "\tyear = {2003},\n",
      "\tpages = {255\u2013264}\n",
      "}"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "class CNG(SCAP):\n",
      "    def compare_profiles(self, profile1, profile2):\n",
      "        ngrams = set(profile1.keys() + profile2.keys())\n",
      "        d1 = np.array([profile1.get(ng, 0.) for ng in ngrams])\n",
      "        d2 = np.array([profile2.get(ng, 0.) for ng in ngrams])\n",
      "        return np.mean(4 * np.square((d1 - d2) / (d1 + d2 + 1e-16)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 150
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Thanks to the fact we inherited from the above `SCAP` class, testing this method uses almost exactly the same code.\n",
      "The only difference in the code below has been highlighted."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Setup a new model, this is the only difference\n",
      "model = CNG(n=4, L=2000)\n",
      "# First, setup the cross fold validation object\n",
      "from sklearn import cross_validation\n",
      "cv = cross_validation.KFold(len(documents), n_folds=5, random_state=0)\n",
      "# Now, compute the score using it\n",
      "scores = cross_validation.cross_val_score(model, documents, classes, cv=cv)\n",
      "print(\"Accuracy is {:.1f}%\".format(100. * np.mean(scores)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Accuracy is 54.2%\n"
       ]
      }
     ],
     "prompt_number": 151
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Slightly better than before. One problem here is that the values of the parameters (`n` and `L` in this case) make a huge difference in the accuracy of the algorithm.\n",
      "To find good algorithms, we can perform a grid search of different parameters.\n",
      "\n",
      "Thanks to the format of the classes, we can use scikit-learn's grid search parameter for this purpose.\n",
      "I've greatly restricted the parameter spaces, as this can take a long time to compute."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# This code snippet is a shortened version of the one in the scikit-learn docs, available here:\n",
      "# http://scikit-learn.org/stable/auto_examples/grid_search_digits.html#example-grid-search-digits-py\n",
      "from sklearn.grid_search import GridSearchCV\n",
      "from sklearn.metrics import classification_report\n",
      "parameters = [{'n': [3, 4], 'L':[500, 1000]},]  # These are the parameters to search through.\n",
      "\n",
      "# The n and L values here don't matter, they will be updated for training.\n",
      "clf = GridSearchCV(CNG(n=1, L=1), parameters, cv=5, scoring='accuracy')\n",
      "clf.fit(documents, classes)\n",
      "print(\"The best model found has n={} and L={}\".format(clf.best_estimator_.n, clf.best_estimator_.L))\n",
      "print(\"The scores for the different parameter sets were:\")\n",
      "for params, mean_score, scores in clf.grid_scores_:\n",
      "    print(\"{}: {:.3f} (+/-{:.3f})\".format(params, mean_score, scores.std() / 2))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The best model found has n=4 and L=500\n",
        "The scores for the different parameter sets were:\n",
        "{'L': 500, 'n': 3}: 0.751 (+/-0.054)\n",
        "{'L': 500, 'n': 4}: 0.814 (+/-0.050)\n",
        "{'L': 1000, 'n': 3}: 0.768 (+/-0.046)\n",
        "{'L': 1000, 'n': 4}: 0.808 (+/-0.037)\n"
       ]
      }
     ],
     "prompt_number": 152
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "From this, we can increase the parameter options and do larger searches.\n",
      "I recommend moving to a more distributed system, and having a look at the `n_jobs` parameter, which can run the search across multiple processes (see the [docs](http://scikit-learn.org/stable/modules/generated/sklearn.grid_search.GridSearchCV.html) for more details).\n",
      "\n",
      "The returned model is a valid scikit-learn estimator, so you can call the below code, which will use this model.\n",
      "Additionally, you can retrieve the actual best estimator instance by calling `clf.best_estimator_`.\n",
      "One final word of caution here about the next snippet: I have used the documents list, which was used in training the model, for the evaluation.\n",
      "You shouldn't do that, instead, you should split the data up into a training and evaluation set, where the evaluation set is used *only at the end* and never in any training or cross validation."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "y_true, y_pred = classes, clf.predict(documents)\n",
      "print(classification_report(y_true, y_pred))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "             precision    recall  f1-score   support\n",
        "\n",
        "          0       0.92      0.92      0.92        51\n",
        "          1       1.00      0.61      0.76        44\n",
        "          2       0.86      0.86      0.86        22\n",
        "          3       0.75      0.90      0.82        10\n",
        "          4       0.85      1.00      0.92        11\n",
        "          5       0.63      0.90      0.74        29\n",
        "          6       0.91      1.00      0.95        10\n",
        "\n",
        "avg / total       0.87      0.84      0.84       177\n",
        "\n"
       ]
      }
     ],
     "prompt_number": 153
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can have a look at the effect that the L parameter has on the accuracy of CNG by plotting it for different n-values.\n",
      "The code below only does n=3, but can be increased to do more values.\n",
      "\n",
      "This code will take a while to run, depending on your computer, maybe half an hour or more."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from sklearn.cross_validation import KFold\n",
      "\n",
      "n_values = [3, ]  # Add more values here to test different n-values\n",
      "L_values = np.arange(500, 3001, 500)  # From 500 to 3000 inclusive, in increments of 500.\n",
      "\n",
      "# Takes too long on the full set of documents.\n",
      "# Set sample_size = None if you want to run on the full dataset.\n",
      "sample_size = None\n",
      "\n",
      "if sample_size is None:\n",
      "    sample_documents = documents\n",
      "    sample_classes = classes\n",
      "else:\n",
      "    sample_documents = documents[:sample_size]\n",
      "    sample_classes = classes[:sample_size]\n",
      "assert len(sample_documents) == len(sample_classes)\n",
      "# The number of folds is the size of the smallest class\n",
      "class_counts = np.bincount(sample_classes)\n",
      "# class_counts could be zero if class doesn't exist, so remove those values\n",
      "class_counts = class_counts[np.nonzero(class_counts)]\n",
      "n_folds = min(class_counts)\n",
      "\n",
      "plt.figure()\n",
      "for n in n_values:\n",
      "    scores = []\n",
      "    for L in L_values:\n",
      "        cv = KFold(len(sample_documents), n_folds=n_folds, random_state=0)\n",
      "        cur_model = CNG(n=n, L=L)\n",
      "        cur_score = cross_validation.cross_val_score(cur_model, sample_documents, sample_classes, cv=cv)\n",
      "        cur_score = np.mean(cur_score)\n",
      "        scores.append(cur_score)\n",
      "        print(\"CNG(n={}, L={}) = {:.4f}\".format(n, L, cur_score))\n",
      "    plt.plot(L_values, scores, label='n={}'.format(n))\n",
      "plt.ylim(0, 1.0)\n",
      "plt.xticks(L_values)\n",
      "plt.show()\n",
      "                                                "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "CNG(n=3, L=500) = 0.6137\n",
        "CNG(n=3, L=1000) = 0.6359"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "CNG(n=3, L=1500) = 0.6193"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "CNG(n=3, L=2000) = 0.6190"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "CNG(n=3, L=2500) = 0.6252"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "CNG(n=3, L=3000) = 0.6088"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAD9CAYAAABdoNd6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAE0xJREFUeJzt3V9s3XX9x/HXmecQ95+NsTLOaazs1PZs604x7QohmINC\nOmc8KMOkJM5llqVpgnPGmBlvaE0clMSLhXpR/gqodV6o5aI76AhHJroWpbDEEVJ0DafH0aRzY92G\nnHr6+V3st9P132m7ftuufT8fyQnn9Hx2zud8+O75/fLtt8XnnHMCACxqS+Z7AgCA2UfsAcAAYg8A\nBhB7ADCA2AOAAcQeAAyYNPbf/va3VVBQoLKysgnH7Nu3T8XFxYpGo+rq6vJ0ggCAmZs09nv27FEi\nkZjw+fb2dr3//vvq7u7WU089pfr6ek8nCACYuUljf/fdd2vNmjUTPv/yyy9r9+7dkqSqqiqdO3dO\nfX193s0QADBjMz5nn06nVVhYmHscCoXU29s705cFAHjI78WLjP6NCz6fb8yY8b4GAJicF7/VZsZH\n9sFgUKlUKve4t7dXwWBw3LHOOW7O6dFHH533OVwvN9aCtWAt8t+8MuPYx+Nxvfjii5Kk48eP68Yb\nb1RBQcGMJwYA8M6kp3Eeeugh/elPf1J/f78KCwvV2NiowcFBSVJdXZ127Nih9vZ2hcNhLV++XM8/\n//ysTxoAMD2Txr61tXXSF2lubvZkMlbEYrH5nsJ1g7UYxloMYy2853NenhTK90Y+n6fnnwDAAq/a\nya9LAAADiD0AGEDsAcAAYg8ABhB7ADCA2AOAAcQeAAwg9gBgALEHAAOIPQAYQOwBwABiDwAGEHsA\nMIDYA4ABxB4ADCD2AGAAsQcAA4g9ABhA7AHAAGIPAAYQewAwgNgDgAHEHgAMIPYAYACxBwADiD0A\nGEDsAcAAYg8ABhB7ADCA2AOAAcQeAAwg9gBgALEHAAOIPQAYQOwBwIBJY59IJFRaWqri4mI1NTWN\neb6/v1/bt29XeXm5tmzZop///OezMU8AwAz4nHNuoiez2axKSkp09OhRBYNBVVZWqrW1VZFIJDem\noaFBn3zyiR577DH19/erpKREfX198vv9I9/I51OetwIAjMOrduY9su/s7FQ4HFZRUZECgYBqamrU\n1tY2YsyGDRt0/vx5SdL58+d10003jQk9AGB+5a1yOp1WYWFh7nEoFFJHR8eIMXv37tUXv/hF3Xrr\nrRoYGNBvfvObCV+voaEhdz8WiykWi13brAFgkUomk0omk56/bt7Y+3y+SV/g4MGDKi8vVzKZ1D//\n+U/dd999euedd7Ry5coxY6+OPQBgrNEHwo2NjZ68bt7TOMFgUKlUKvc4lUopFAqNGPOXv/xF3/jG\nNyRJGzdu1Gc/+1m99957nkwOAOCNvLGvqKhQd3e3enp6lMlkdPjwYcXj8RFjSktLdfToUUlSX1+f\n3nvvPd12222zN2MAwLTlPY3j9/vV3Nys6upqZbNZ1dbWKhKJqKWlRZJUV1enH/3oR9qzZ4+i0aiG\nhob0xBNPaO3atXMyeQDA1OS99NLTN+LSSwCYtjm59BIAsDgQewAwgNgDgAHEHgAMIPYAYACxBwAD\niD0AGEDsAcAAYg8ABhB7ADCA2AOAAcQeAAwg9gBgALEHAAOIPQAYQOwBwABiDwAGEHsAMIDYA4AB\nxB4ADCD2AGAAsQcAA4g9ABhA7AHAAGIPAAYQewAwgNgDgAHEHgAMIPYAYACxBwAD/PM9gcUmk5HO\nnMl/k6TVq6VVq4ZvEz1esUL61Kfm9zMBWPh8zjk3J2/k82mO3soTzkkffTR+rP/zn4lD/skn0tq1\n0k03Df9z9E2Szp+/fPvoo+H74z2+cEFavjz/DmEqj9lpAAuTV+00EfupHG2Pvp09Ky1dOn6s891W\nrpR8Pu/mPjQkDQzk3yFMZadx8eLwTuNadxirVl3+fEs4+YdpGhqSBgcv/10cHBy+Xf346vuf/vTI\n7W7pUm//Xi0kJmOf72g731H31Ufb+W5Xj1m7VrrhBo8+/HVgNnca091hsNOYHuekbHZqkZxqTOd6\n3NCQFAgM3264Yfz7V27//e/l7e3KNjs4OHYbGu+2cmX+5xfif+Eu+Nhf69H2smXjx3kuj7Yt83qn\nMXqHOvrfU77HszV2vt4nm504oNJwFEfHMV84r5fn/P6Z/R3MZMZud1O5jf4zV58WnclOY9WquTsY\nnLPYJxIJ7d+/X9lsVg8//LAOHDgwZkwymdT3vvc9DQ4Oat26dUomk+NO+POfd9M+2h595B0IzPgz\n4zpwZadxJWTS5SPYq+V7PFtj5/N9/P6JI7rQjkavV0NDl4M/kx3GlYMXv3/mO4xVqy4fwObbEc5J\n7LPZrEpKSnT06FEFg0FVVlaqtbVVkUgkN+bcuXO666679MorrygUCqm/v1/r1q0bd8Jvvuk42gaw\n4Dk3fKrpWnYWV98ymfw7hmee8Sb2eS+97OzsVDgcVlFRkSSppqZGbW1tI2L/q1/9Sjt37lQoFJKk\ncUN/RUXFjOcLAPPO57v8TeOlS6WCgpm91uBg/h3CM894M+e8sU+n0yosLMw9DoVC6ujoGDGmu7tb\ng4ODuueeezQwMKDvfve72rVrlzezA4BFLhC4fIp67drxn6+r8+Z98sbeN4XzLIODg3rrrbf06quv\n6tKlS7rzzjt1xx13qLi4eMzYhoaG3P1YLKZYLDbtCQPAYpZMJsf9vudM5Y19MBhUKpXKPU6lUrnT\nNVcUFhZq3bp1Wrp0qZYuXaovfOELeueddyaNPQBgrNEHwo2NjZ68bt4rnSsqKtTd3a2enh5lMhkd\nPnxY8Xh8xJj7779ff/7zn5XNZnXp0iV1dHRo06ZNnkwOAOCNvEf2fr9fzc3Nqq6uVjabVW1trSKR\niFpaWiRJdXV1Ki0t1fbt27V161YtWbJEe/fuJfYAcJ1ZUD9BCwDWeNVOfmAdAAwg9gBgALEHAAOI\nPQAYQOwBwABiDwAGEHsAMIDYA4ABxB4ADCD2AGAAsQcAA4g9ABhA7AHAAGIPAAYQewAwgNgDgAHE\nHgAMIPYAYACxBwADiD0AGEDsAcAAYg8ABhB7ADCA2AOAAcQeAAwg9gBgALEHAAOIPQAYQOwBwABi\nDwAGEHsAMIDYA4ABxB4ADCD2AGAAsQcAA4g9ABgwaewTiYRKS0tVXFyspqamCce9+eab8vv9+u1v\nf+vpBAEAM5c39tlsVo888ogSiYROnjyp1tZWvfvuu+OOO3DggLZv3y7n3KxNFgBwbfLGvrOzU+Fw\nWEVFRQoEAqqpqVFbW9uYcU8++aQefPBB3XzzzbM2UQDAtcsb+3Q6rcLCwtzjUCikdDo9ZkxbW5vq\n6+slST6fbxamCQCYCX++J6cS7v379+vxxx+Xz+eTcy7vaZyGhobc/VgsplgsNuWJAoAFyWRSyWTS\n89f1uTx1Pn78uBoaGpRIJCRJjz32mJYsWaIDBw7kxtx22225wPf392vZsmV6+umnFY/HR77R/+8M\nAABT51U788b+f//7n0pKSvTqq6/q1ltv1bZt29Ta2qpIJDLu+D179uirX/2qHnjggVmbMABY4lU7\n857G8fv9am5uVnV1tbLZrGpraxWJRNTS0iJJqqurm/EEAACzL++RvadvxJE9AEybV+3kJ2gBwABi\nDwAGEHsAMIDYA4ABxB4ADCD2AGAAsQcAA4g9ABhA7AHAAGIPAAYQewAwgNgDgAHEHgAMIPYAYACx\nBwADiD0AGEDsAcAAYg8ABhB7ADCA2AOAAcQeAAwg9gBgALEHAAOIPQAYQOwBwABiDwAGEHsAMIDY\nA4ABxB4ADCD2AGAAsQcAA4g9ABhA7AHAAGIPAAYQewAwgNgDgAFTin0ikVBpaamKi4vV1NQ05vlf\n/vKXikaj2rp1q+666y6dOHHC84kCAK6dzznn8g3IZrMqKSnR0aNHFQwGVVlZqdbWVkUikdyYv/71\nr9q0aZNWr16tRCKhhoYGHT9+fOQb+Xya5K0AAKN41c5Jj+w7OzsVDodVVFSkQCCgmpoatbW1jRhz\n5513avXq1ZKkqqoq9fb2znhiAADv+CcbkE6nVVhYmHscCoXU0dEx4fhnn31WO3bsGPe5hoaG3P1Y\nLKZYLDb1mQKAAclkUslk0vPXnTT2Pp9vyi/22muv6bnnntMbb7wx7vNXxx4AMNboA+HGxkZPXnfS\n2AeDQaVSqdzjVCqlUCg0ZtyJEye0d+9eJRIJrVmzxpPJAQC8Mek5+4qKCnV3d6unp0eZTEaHDx9W\nPB4fMeaDDz7QAw88oF/84hcKh8OzNlkAwLWZ9Mje7/erublZ1dXVymazqq2tVSQSUUtLiySprq5O\nP/7xj3X27FnV19dLkgKBgDo7O2d35gCAKZv00kvP3ohLLwFg2ubs0ksAwMJH7AHAAGIPAAYQewAw\ngNgDgAHEHgAMIPYAYACxBwADiD0AGEDsAcAAYg8ABhB7ADCA2AOAAcQeAAwg9gBgALEHAAOIPQAY\nQOwBwABiDwAGEHsAMIDYA4ABxB4ADCD2AGAAsQcAA4g9ABhA7AHAAGIPAAYQewAwgNgDgAHEHgAM\nIPYAYACxBwADiD0AGEDsAcAAYg8ABkwa+0QiodLSUhUXF6upqWncMfv27VNxcbGi0ai6uro8n+Ri\nk0wm53sK1w3WYhhrMYy18F7e2GezWT3yyCNKJBI6efKkWltb9e67744Y097ervfff1/d3d166qmn\nVF9fP6sTXgzYkIexFsNYi2Gshffyxr6zs1PhcFhFRUUKBAKqqalRW1vbiDEvv/yydu/eLUmqqqrS\nuXPn1NfXN3szBgBMW97Yp9NpFRYW5h6HQiGl0+lJx/T29no8TQDATPjzPenz+ab0Is65Kf25qb6e\nBY2NjfM9hesGazGMtRjGWngrb+yDwaBSqVTucSqVUigUyjumt7dXwWBwzGuN3iEAAOZO3tM4FRUV\n6u7uVk9PjzKZjA4fPqx4PD5iTDwe14svvihJOn78uG688UYVFBTM3owBANOW98je7/erublZ1dXV\nymazqq2tVSQSUUtLiySprq5OO3bsUHt7u8LhsJYvX67nn39+TiYOAJgG55HPfOYzrqyszJWXl7vK\nykrnnHNnzpxx9957rysuLnb33XefO3v2bG78wYMHXTgcdiUlJe6VV17xahrzYs+ePW79+vVuy5Yt\nua9dy2f/29/+5rZs2eLC4bDbt2/fnH4Gr4y3Fo8++qgLBoOuvLzclZeXu/b29txzi3ktPvjgAxeL\nxdymTZvc5s2b3aFDh5xzNreNidbC4rbx8ccfu23btrloNOoikYj74Q9/6Jyb/e3Cs9gXFRW5M2fO\njPjaD37wA9fU1OScc+7xxx93Bw4ccM45949//MNFo1GXyWTcqVOn3MaNG102m/VqKnPu9ddfd2+9\n9daIwE3nsw8NDTnnnKusrHQdHR3OOee+/OUvuyNHjszxJ5m58daioaHB/fSnPx0zdrGvxenTp11X\nV5dzzrmBgQH3uc99zp08edLktjHRWljdNi5evOicc25wcNBVVVW5Y8eOzfp24emvS3Cjvgl79TX4\nu3fv1u9//3tJUltbmx566CEFAgEVFRUpHA6rs7PTy6nMqbvvvltr1qwZ8bXpfPaOjg6dPn1aAwMD\n2rZtmyTpW9/6Vu7PLCTjrYU0/jfoF/ta3HLLLSovL5ckrVixQpFIROl02uS2MdFaSDa3jWXLlkmS\nMpmMstms1qxZM+vbhWex9/l8uvfee1VRUaGnn35aktTX15f7Zm1BQUHuh63+/e9/j7iqZ7zr9xe6\n6X720V8PBoOLak2efPJJRaNR1dbW6ty5c5JsrUVPT4+6urpUVVVlftu4shZ33HGHJJvbxtDQkMrL\ny1VQUKB77rlHmzdvnvXtwrPYv/HGG+rq6tKRI0f0s5/9TMeOHRvxvM/ny3ud/WK+Bn+yz77Y1dfX\n69SpU3r77be1YcMGff/735/vKc2pCxcuaOfOnTp06JBWrlw54jlr28aFCxf04IMP6tChQ1qxYoXZ\nbWPJkiV6++231dvbq9dff12vvfbaiOdnY7vwLPYbNmyQJN188836+te/rs7OThUUFOjDDz+UJJ0+\nfVrr16+XNPVr8xey6Xz2UCikYDA44iePF9OarF+/PrfxPvzww7lTdhbWYnBwUDt37tSuXbv0ta99\nTZLdbePKWnzzm9/MrYXlbUOSVq9era985Sv6+9//PuvbhSexv3TpkgYGBiRJFy9e1B/+8AeVlZUp\nHo/rhRdekCS98MILuX/B8Xhcv/71r5XJZHTq1Cl1d3fnzjstFtP97LfccotWrVqljo4OOef00ksv\n5f7MQnf69Onc/d/97ncqKyuTtPjXwjmn2tpabdq0Sfv378993eK2MdFaWNw2+vv7c6erPv74Y/3x\nj3/U7bffPvvbhRffWf7Xv/7lotGoi0ajbvPmze7gwYPOucuXEn3pS18a91Kin/zkJ27jxo2upKTE\nJRIJL6Yxb2pqatyGDRtcIBBwoVDIPffcc9f02a9cRrVx40b3ne98Zz4+yoyNXotnn33W7dq1y5WV\nlbmtW7e6+++/33344Ye58Yt5LY4dO+Z8Pp+LRqO5SwuPHDlictsYby3a29tNbhsnTpxwt99+u4tG\no66srMw98cQTzrlr6+V01sLnHL/HAAAWO/5PVQBgALEHAAOIPQAYQOwBwABiDwAGEHsAMOD/AKfU\nur5eqdZXAAAAAElFTkSuQmCC\n"
      }
     ],
     "prompt_number": 162
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "3) Recentred Local Profiles (RLP)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "RLP was created by myself, and has a number of differences to the first two algorithms presented.\n",
      "At its core though, it is based on the principles as CNG (and therefore SCAP).\n",
      "\n",
      "In the original derivation of CNG, the work of Bennet (1976), which used a similar model, with a language default value.\n",
      "This language default is effectively a profile of all documents in that language.\n",
      "Then, rather than compare two profiles, we first normalise the profiles against what is the normal value for that n-gram.\n",
      "Intuitively, while the `the` n-gram is very common for most authors, it is expected to be, and therefore the presence of this n-gram is hardly surprising.\n",
      "By comparing to the language default, we see that this n-gram is expected and therefore are less surprised.\n",
      "A profile is then computed by first computing the CNG-esque profile of the author, then subtracting (recentering) the values from the corpus default profile.\n",
      "These profiles are then compared in RLP.\n",
      "\n",
      "RLP does a number of other things differently as well.\n",
      "Rather then setting n-gram frequencies to zero if they are not present in the profile (as in CNG), the actual frequency value is used for comparisons, even if that n-gram is not present.\n",
      "Profiles are compared using the following equation, where $||P_i||_2$ is the Euclidean distance of the profile (i.e. the sum of the squares of each frequency value, and $P_i \\cdot P_j$ is the dot product of the n-gram values in each profile, for all n-grams in the profiles' union.\n",
      "Another major change is that the `top L` n-grams are ranked by absolute, rather than actual highest.\n",
      "This means an n-gram with a value of -0.15 occurs higher than an n-gram with value 0.14.\n",
      "\n",
      "                                                                                                            \n",
      "Two profiles are then compared using the following equation, which is also known as the cosine distance.                                                                                                             \n",
      "\n",
      "$$ d_{RLP}(P_1, P_2) = 1 - \\frac{P_1 \\cdot P_2}{||P_1||_2 ||P_2||_2}  $$\n",
      "\n",
      "References:\n",
      "    "
     ]
    },
    {
     "cell_type": "raw",
     "metadata": {},
     "source": [
      "@book{bennett1976scientific,\n",
      "  title={Scientific and engineering problem-solving with the computer},\n",
      "  author={Bennett, William Ralph},\n",
      "  year={1976},\n",
      "  publisher={Prentice Hall PTR}\n",
      "}\n",
      "\n",
      "@article{layton2012recentred,\n",
      "  title={Recentred local profiles for authorship attribution.},\n",
      "  author={Layton, Robert and Watters, Paul Andrew and Dazeley, Richard},\n",
      "  journal={Natural Language Engineering},\n",
      "  volume={18},\n",
      "  number={3},\n",
      "  pages={293--312},\n",
      "  year={2012},\n",
      "  publisher={Cambridge Univ Press}\n",
      "}"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.spatial import distance\n",
      "\n",
      "\n",
      "class RLP(SCAP):\n",
      "    \n",
      "    def __init__(self, n, L):\n",
      "        super(RLP, self).__init__(n, L)\n",
      "        self.topLonly = False  # Don't limit profiles on creation\n",
      "        self.language_profile = None\n",
      "    \n",
      "    def compare_profiles(self, profile1, profile2):\n",
      "        # All n-grams in the top L of either profile.\n",
      "        ngrams = set(self.top_L(profile1).keys() + self.top_L(profile1).keys())\n",
      "        # Profile vector for profile 1\n",
      "        d1 = np.array([profile1.get(ng, 0.) for ng in ngrams])\n",
      "        # Profile vector for profile 2\n",
      "        d2 = np.array([profile2.get(ng, 0.) for ng in ngrams])\n",
      "        return distance.cosine(d1, d2)\n",
      "    \n",
      "    def top_L(self, profile):\n",
      "        threshold = sorted(map(abs, profile.values()))[-self.L]\n",
      "        copy = defaultdict(float)\n",
      "        for key in profile:\n",
      "            if abs(profile[key]) >= threshold:\n",
      "                copy[key] = profile[key]\n",
      "        return copy\n",
      "    \n",
      "    def fit(self, documents, classes):\n",
      "        self.language_profile = self.create_profile(documents)\n",
      "        super(RLP, self).fit(documents, classes)\n",
      "        \n",
      "    def create_profile(self, documents):\n",
      "        # Creates a profile of a document or list of documents.\n",
      "        if isinstance(documents, basestring):\n",
      "            # documents can be either a list of documents, or a single document.\n",
      "            # if it's a single document, convert to a list\n",
      "            documents = [documents,]\n",
      "        # profile each document independently\n",
      "        profiles = (count_ngrams(document, self.n, normalise=False)\n",
      "                    for document in documents)\n",
      "        # Merge the profiles\n",
      "        main_profile = defaultdict(float)\n",
      "        for profile in profiles:\n",
      "            for ngram in profile:\n",
      "                main_profile[ngram] += profile[ngram]\n",
      "        # Normalise the profile\n",
      "        num_ngrams = float(sum(main_profile.values()))\n",
      "        for ngram in main_profile:\n",
      "            main_profile[ngram] /= num_ngrams\n",
      "        if self.language_profile is not None:\n",
      "            # Recentre profile.\n",
      "            for key in main_profile:\n",
      "                main_profile[key] = main_profile.get(key, 0) - self.language_profile.get(key, 0)\n",
      "        # Note that the profile is returned in full, as exact frequencies are used\n",
      "        # in comparing profiles (rather than chopped off)\n",
      "        return main_profile\n",
      "        "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 155
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we test with out overfitting code from earlier."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "model = RLP(n=4, L=2000)\n",
      "# Training the model\n",
      "model.fit(documents, classes)\n",
      "# Predict the author of each of the documents\n",
      "y_pred = model.predict(documents)  # Notice how I am predicting the data I used for training? Don't do that.\n",
      "# Compute the accuracy of these predictions.\n",
      "print(\"Accuracy is {:.1f}%\".format(100. * np.mean(classes == y_pred)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Accuracy is 75.1%\n"
       ]
      }
     ],
     "prompt_number": 156
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Have a look at some of the values in the language profile\n",
      "for ngram in ['the ', \"'ell\", \" and\", \"est \", \"ever\"]:\n",
      "    print(\"{} = {:.6f}\".format(repr(ngram), model.language_profile[ngram]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "'the ' = 0.007964\n",
        "\"'ell\" = 0.000011\n",
        "' and' = 0.005158\n",
        "'est ' = 0.000345\n",
        "'ever' = 0.000779\n"
       ]
      }
     ],
     "prompt_number": 157
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Have a look at one of the author profiles\n",
      "ap = model.author_profiles[1]\n",
      "print(\"Length of profile is {}\".format(len(ap))) # Notice how this is much higher than L (which is currently 2000)\n",
      "for ngram in ['the ', \"'ell\", \" and\", \"est \", \"ever\"]:\n",
      "    print(\"{} = {:.6f}\".format(repr(ngram), ap[ngram]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Length of profile is 158363\n",
        "'the ' = -0.001139\n",
        "\"'ell\" = 0.000021\n",
        "' and' = -0.000118\n",
        "'est ' = 0.000039\n",
        "'ever' = -0.000027\n"
       ]
      }
     ],
     "prompt_number": 158
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Notice how the n-gram `est ` has a negative frequency.\n",
      "This means the frequency for this author was 0.000089 lower than the frequency for the language default.\n",
      "RLP uses this information, where it may be lost in other LNG methods.\n",
      "The cost, as you may be discovered if you are running the code as we go, is that RLP is *significantly* slower than CNG or SCAP.\n",
      "Experiments I've conducted show that RLP is typically higher than CNG (which is in turn typically higher than SCAP) for most corpora.\n",
      "\n",
      "Next, we do a more appropriate test using cross validation."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Setup a new model\n",
      "model = RLP(n=4, L=2000)\n",
      "# First, setup the cross fold validation object\n",
      "from sklearn import cross_validation\n",
      "cv = cross_validation.KFold(len(documents), n_folds=5, random_state=0)\n",
      "# Now, compute the score using it\n",
      "scores = cross_validation.cross_val_score(model, documents, classes, cv=cv)\n",
      "print(\"Accuracy is {:.1f}%\".format(100. * np.mean(scores)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Accuracy is 50.7%\n"
       ]
      }
     ],
     "prompt_number": 159
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "That's lower than before, but with different models comes different parameters for improving performance.\n",
      "A grid search through some parameter values will give a better score (but will take a while to run).\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from sklearn.grid_search import GridSearchCV\n",
      "from sklearn.metrics import classification_report\n",
      "parameters = [{'n': [3, 4], 'L':[500, 1000]},]  # These are the parameters to search through.\n",
      "\n",
      "# The n and L values here don't matter, they will be updated for training.\n",
      "clf = GridSearchCV(RLP(n=1, L=1), parameters, cv=5, scoring='accuracy')\n",
      "clf.fit(documents, classes)\n",
      "print(\"The best model found has n={} and L={}\".format(clf.best_estimator_.n, clf.best_estimator_.L))\n",
      "print(\"The scores for the different parameter sets were:\")\n",
      "for params, mean_score, scores in clf.grid_scores_:\n",
      "    print(\"{}: {:.3f} (+/-{:.3f})\".format(params, mean_score, scores.std() / 2))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The best model found has n=3 and L=1000\n",
        "The scores for the different parameter sets were:\n",
        "{'L': 500, 'n': 3}: 0.723 (+/-0.038)\n",
        "{'L': 500, 'n': 4}: 0.695 (+/-0.048)\n",
        "{'L': 1000, 'n': 3}: 0.734 (+/-0.037)\n",
        "{'L': 1000, 'n': 4}: 0.701 (+/-0.047)\n"
       ]
      }
     ],
     "prompt_number": 160
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 160
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 160
    }
   ],
   "metadata": {}
  }
 ]
}