defactoring-pace-of-change.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Defactoring Pace of Change\n",
    "\n",
    "At a practical level, we have copied the code from Underwood and Sellers' Pace of Change git repository and restructured their code into a Jupyter Notebook. Transforming Underwood and Seller's code from a set of Python scripts into the notebook involved a series of transformations to their code. The provenance of this effort has been preserved because we did all of our work in a fork of their version repository, visible via the git history command.\n",
    "\n",
    "Underwood and Sellers's wrote their code to make it easy to replicate their results by running a single command. One of our first tasks was to trace their code following the path of execution from the initial command to the output data files. By reading the code that spanned multiple python scripts, we could obtain a rough idea of how Underwood and Sellers's analysis of poetry volumes was performed in practice. The first step of reading their code involved translating their python scripts into a Jupyter Notebook. Plain text files leave much to be desired with trying to critically interrogate code, the Jupyter Notebook format and platform enables a much richer reading experience. Translating Underwood and Sellers's scripts involved copying snippets of code from multiple Python files into a single executable notebook. For code that was not already authored in a notebook, this requires significant effort curating and restructuring the code to fit the linear structure of the notebook while preserving the functionality of Underwood and Sellers's scripts.\n",
    "\n",
    "Not all of the code in the git repository was copied into the notebook for two main reasons. First, not all of the code included in the files was actually needed to replicate the analysis. The repository includes additional code, one could call it \"cruft\", from Underwood and Seller's exploratory analysis or earlier iterations of their analytical process. For example, the file SonicScrewdiver.py is never called upon although, based upon the name, we might hypothesize it is an important catch-all module for a variety of different functions. Other sections of the code are commented out and thus never executed (such as binormal_select discussed below). As part of the defactoring process, we opted to not include unused functions or code blocks. Another reason that the code of Sellers and Underwood was not copied verbatim is because we do not represent all possible executions paths of the code. Their code included execution paths that trained multiple models of the data broken into twenty-year periods, but their final analysis trained a model on the entire dataset. In this notebook we follow the default option, \"full\", which models the entire dataset.\n",
    "\n",
    "Underwood and Seller's code is stored within four Python files, `replicate.py`, `parallel_crossvalidate.py`, `metafilter.py`, and `modelingprocess.py`. Through the process of defactoring we have copied some, but not all, code from the Python text files to a Jupyter Notebook. This process has transformed the structure and flow of the code, but not the functionality. Like *refactoring*, the inputs and outputs and overall functionality remain the same, but our process has restructured the code to be much more readable as a narrative with our own expository text intermixed with executable blocks of code."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Defactoring code from python files to the notebook](notebook_resources/defactoring.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order for the code to execute seamlessly within the notebook, we had to make minor changes and tweaks to the code. These changes fall into 4 categories:\n",
    "\n",
    "* *defactoring functions* - This is the most significant of the changes. When we defactor a function we take the function's code and move it to the global namespace. This has the effect of eliminating the function and just making it part of the main execution path. \n",
    "* *defactoring definitions* - Not all functions can be fully defactored. Functions that are called more than once or those that are short have been kept as re-usable functions. Defactoring Definition cells define the functions above the code cells that use them (preventing errors).\n",
    "* *defactoring namespace* - Because we have defactored some of the functions and their function calls some of the variables in the namespace need to be mapped to each other. This happens, for instance, when the return value of a defactored function needs to be stored in a differently named variable so the rest of the code can access it (resolving the fact we have moved a function's internal variables into a global namespace).\n",
    "* *defactoring inspections* - When we want to inspect the state of the process, we insert an inspection cell that prints the values of the variables of interest. These allow us the ability to inspect the state of the process as it is running; essentially debugging intermediate states.\n",
    "* *defactoring import* - Because the code is reliant upon external and third party functions, we need to import that code into the global namespace. This cell contains all of those imports."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### How to Read this Notebook\n",
    "\n",
    "One of the advantages to defactoring a function is that it affords us the ability to insert critical commentary into the code constituting the function itself. An unfortunate side-effect is that this makes keeping track of one's place in the code a little more difficult. We have restructured the code into an executable narrative (when read as a Jupyter Notebook) for both human and machine readers.\n",
    "\n",
    "Jupyter Notebooks are divided into a series of cells, *code cells* and *text cells*. In our defactoring, the code cells contain Python code mostly written by Underwood and Sellers, but with a few modifications by us. Code cells can be executed when reading this document using the Jupyter Platform (as opposed to static document). We have broken up the code of Underwood and Sellers and inserted text cells, written by us, which discuss and explain what is happening in the code cells and underlying computational process. In most cases the expository text is below the code cell. In the critical examination of the code we will be using a fair bit of coding language jargon. Although a very deep grasp of these coding concepts will not be needed, a certain level of coding literacy (Vee 2013) will facilitate a speedier comprehension of our narrative."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Program Structure\n",
    "\n",
    "The code expressed below has nine steps:\n",
    "* [Importing Libraries](#Importing-Libraries) - Loads the necessary Python libraries needed for the analysis.\n",
    "* [Setting Parameters](#Setting-Parameters) - Specifies parameters for the loading, cleaning, and labeling of data as well as sets conditions for the logistic regression.\n",
    "* [Preparing Metadata](#Preparing-MetaData) - Generates a list of *.tsv files from the `poems/` directory. \n",
    "    * [Cleaning Metadata](#Cleaning-Metadata) - Loads the metadata file, `poemetadata.csv` and performs some cleaning of the metadata to make labeling easier.\n",
    "    * [Sorting Training Data](#Sorting-Data) - Sort the volumes into two bins, reviewed and not reviewed using the cleaned metadata.\n",
    "* [Transforming Words into Features](#Transforming-Words-into-Features) - Identifies the 3,200 most common words in the corpus. Those most common words will be the features for the regression.\n",
    "    * [Filtering Authors](#Filtering-Authors) - Removes poems by authors who have been reviewed.\n",
    "    * [Filtering Words](#Filtering-Words) - Remove any words from the poem data that are not in the most-common feature list.\n",
    "* [Training Predictive Models](#Training-Predictive-Models) - Run a separate logistic regression for each volume, using a single volume as held-out data and measure each model's predictive power.\n",
    "* [Modeling Coefficients](#Modeling-Coefficients) - Run a single logistic regression over all the data to inspect the salient coefficients.\n",
    "* [Plotting Results](#Plotting-Results) - Generate a plot showing the accuracy of the predictive models.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Importing Libraries\n",
    "\n",
    "This section loads a series of libraries used in the Pace of Change analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "### DEFACTORING IMPORT\n",
    "import os\n",
    "import csv \n",
    "import random \n",
    "from collections import Counter \n",
    "\n",
    "import numpy as np \n",
    "import pandas as pd \n",
    "\n",
    "#from multiprocessing import Pool \n",
    "### Defactoring Import\n",
    "from multiprocess import Pool \n",
    "\n",
    "\n",
    "\n",
    "import matplotlib.pyplot as plt \n",
    "%matplotlib inline\n",
    "\n",
    "from sklearn.linear_model import LogisticRegression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We begin examination of the code by importing a series of Python libraries into working memory. This is the boundary between the layers of bespoke code and existing general purpose (os, csv, random) and scientific computing libraries (numpy, pandas, sklearn). Following from Hinsen's model of layers of scientific software, what is missing is the inclusion of libraries from the *disciplinary* layer. The most specific library in terms of use in the Hinsen model is the LogisticRegression model from Scikit-Learn, but we'd argue this lives in layer two, scientific software, because it is broadly applicable across a variety of disciplines. This begs the question, what or where are the disciplinary Python libraries for literary history or digital humanities? What functions would they perform? What domain specific tasks or methods need to encoded into a disciplinary library? Perhaps it is too early for such libraries to exist as the practices of computational and data intensive research are still new (in literary history)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting Parameters\n",
    "\n",
    "The first section of the code sets a series of parameters specifying what data to process, where data are located, and parameters for the logistic regression. While there is no complex logic or work being done in this section, many assumptions and important distinctions that shape the execution of subsequent code are defined here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "## PATHS.\n",
    "\n",
    "sourcefolder = 'poems/'\n",
    "extension = '.poe.tsv'\n",
    "classpath = 'poemeta.csv'\n",
    "outputpath = 'mainmodelpredictions.csv'\n",
    "\n",
    "## EXCLUSIONS.\n",
    "\n",
    "excludeif = dict()\n",
    "excludeif['pubname'] = 'TEM'\n",
    "# We're not using reviews from Tait's.\n",
    "\n",
    "excludeif['recept'] = 'addcanon'\n",
    "# We don't ordinarily include canonical volumes that were not in either sample.\n",
    "# These are included only if we're testing the canon specifically.\n",
    "\n",
    "excludeifnot = dict()\n",
    "excludeabove = dict()\n",
    "excludebelow = dict()\n",
    "\n",
    "excludebelow['firstpub'] = 1700\n",
    "excludeabove['firstpub'] = 1950\n",
    "sizecap = 360\n",
    "\n",
    "# For more historically-interesting kinds of questions, we can limit the part\n",
    "# of the dataset that gets TRAINED on, while permitting the whole dataset to\n",
    "# be PREDICTED. (Note that we always exclude authors from their own training\n",
    "# set; this is in addition to that.) The variables futurethreshold and\n",
    "# pastthreshold set the chronological limits of the training set, inclusive\n",
    "# of the threshold itself.\n",
    "\n",
    "## THRESHOLDS\n",
    "\n",
    "futurethreshold = 1925\n",
    "pastthreshold = 1800\n",
    "\n",
    "# CLASSIFY CONDITIONS\n",
    "\n",
    "positive_class = 'rev'\n",
    "category2sorton = 'reviewed'\n",
    "datetype = 'firstpub'\n",
    "numfeatures = 3200\n",
    "regularization = .00007\n",
    "\n",
    "\n",
    "paths = (sourcefolder, extension, classpath, outputpath)\n",
    "exclusions = (excludeif, \n",
    "              excludeifnot, \n",
    "              excludebelow, \n",
    "              excludeabove, \n",
    "              sizecap)\n",
    "thresholds = (pastthreshold, \n",
    "              futurethreshold)\n",
    "classifyconditions = (category2sorton, \n",
    "                      positive_class, \n",
    "                      datetype, \n",
    "                      numfeatures, \n",
    "                      regularization)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The parameters defined in the code cell above are a set of knobs and switches used to tweak the performance and execution of the computational modeling process. Underwood and Sellers have collected the parameters into four categories: *paths*, *exclusions*, *thresholds*, and *classifyconditions*. These categories are simultaneously distinguished discursively through the code comments (the lines beginning with a #) and technologically through four variable assignments, `exclusions`, `thresholds`, and `classifyconditions`. Because technically speaking the grouping of parameters is not strictly necessary, each of these four variables embody stylistic choices of the authors as a means of organizing and structuring the information they are encoding in Python.    \n",
    "\n",
    "The variables in `paths` specify the location of the data and metadata files as well as where to write the output files at the completion of the analysis. The variables in `exclusions` specify data and types of data to be excluded from the analysis, such as reviews from *Tait's Endinburgh Magazine* (https://en.wikipedia.org/wiki/Tait%27s_Edinburgh_Magazine), which we infer from the author's comments. Additional exclusions specify temporal boundaries from 1700 to 1950. A further set of two variables in `thresholds` also articulates a temporal boundary from 1800 to 1925. The comments indicate this distinguishes the temporal window for datasets used in *training* versus those used during *prediction.* The variables in `classifyconditions` are important parameters for the logistic regression, specifying the number of variables to train the model upon as well as setting the regularization parameter (`regularization`) for the logistic regression. What is not well documented here, is why the value .00007 was chosen over other values.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preparing Metadata\n",
    "\n",
    "With the preparation of metadata we begin to see some logical work of Pace of Change being conducted. The code in this section has two subsections, one to clean the metadata and another to sort the training data. All of the work in this section focuses on preparing the metadata, identified in the `classpath` variable and the filenames of the individual data files in the `sourcefolder`. The main task of this section is to organize the metadata of the volumes and their associated labels (positive or negative) for training the logistic regression. All of the code in this section attends to the cleanliness of the metadata; we will not start digging into the data itself until the next section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "### DEFACTORING FUNCTION\n",
    "### def create_model(paths, exclusions, thresholds, classifyconditions):\n",
    "''' This is the main function in the module.\n",
    "It can be called externally; it's also called\n",
    "if the module is run directly.\n",
    "'''\n",
    "verbose = False\n",
    "\n",
    "if not sourcefolder.endswith('/'):\n",
    "    sourcefolder = sourcefolder + '/'\n",
    "\n",
    "# This just makes things easier.\n",
    "\n",
    "# Get a list of files.\n",
    "allthefiles = os.listdir(sourcefolder)\n",
    "# random.shuffle(allthefiles)\n",
    "\n",
    "volumeIDs = list()\n",
    "volumepaths = list()\n",
    "\n",
    "for filename in allthefiles:\n",
    "\n",
    "    if filename.endswith(extension):\n",
    "        volID = filename.replace(extension, \"\")\n",
    "        # The volume ID is basically the filename minus its extension.\n",
    "        # Extensions are likely to be long enough that there is little\n",
    "        # danger of accidental occurrence inside a filename. E.g.\n",
    "        # '.fic.tsv'\n",
    "        path = sourcefolder + filename\n",
    "        volumeIDs.append(volID)\n",
    "        volumepaths.append(path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This code assembles a list of volume identifiers (`volumeIDs`) and file paths (`volumepaths`) by readings the directory listing of files in the `poems/` directory (`sourcefolder`). The filenames are in and of themselves a source of metadata, but as we will see in the code below, they need to be reconciled with the metadata stored separately from the data files. \n",
    "\n",
    "We are curious about the contents of the `volumeIDs` and `volumepaths` variables. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The first item in volumeIDs is:  uc2.ark+=13960=t06w98f3m\n",
      "The first item in volumepaths is:  poems/uc2.ark+=13960=t06w98f3m.poe.tsv\n"
     ]
    }
   ],
   "source": [
    "### DEFACTORING INSPECTION\n",
    "### Inspect the two variables defined in the codecell above.\n",
    "### We know they are lists so lets just look at the first item.\n",
    "print(\"The first item in volumeIDs is: \", volumeIDs[0])\n",
    "print(\"The first item in volumepaths is: \",volumepaths[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code has created an alignment between identifiers in the metadata records and the filename identifiers of the TSV data files themselves (located in the `poems/` folder). These identifiers, `dul1.ark+=13960=t5fb5xg2z`, are the threads that stitch together the various representations of the (meta)data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Cleaning Metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "### DEFACTORING FUNCTION DEFINITION\n",
    "### we need these helper functions for execute the next code cell\n",
    "\n",
    "def dirty_pairtree(htid):\n",
    "    period = htid.find('.')\n",
    "    prefix = htid[0:period]\n",
    "    postfix = htid[(period+1): ]\n",
    "    if '=' in postfix:\n",
    "        postfix = postfix.replace('+',':')\n",
    "        postfix = postfix.replace('=','/')\n",
    "    dirtyname = prefix + \".\" + postfix\n",
    "    return dirtyname\n",
    "\n",
    "def forceint(astring):\n",
    "    try:\n",
    "        intval = int(astring)\n",
    "    except:\n",
    "        intval = 0\n",
    "\n",
    "    return intval"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This code cell defines two functions used in the code below. The first is `dirty_pairtree()`, which cleans up the identifiers in the data. This issue arises from the fact that the HathiTrust (where Underwood and Sellers got their data) uses IDs that cannot be expressed on the filesystem and Underwood and Sellers encoded ID metadata in filenames. The `/` and `:` characters in the IDs cannot be part of a file name. So, because the volumes are stored as individual files they have a `+` and an `=` instead. However, the IDs are stored in the original format in the metadata file so the IDS have to be transformed back into the original HathiTrust format.  The second function is called `forceint()` and transforms numbers expressed as Python strings into the python integer data type with a bit of error handling in the case values that throw and error when being cast as an integer.\n",
    "\n",
    "What does the metadata look like? We can inspect the beginning of file to get a sense of the material conditions of the metadata."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The metadata files has 728 rows and 22 columns.\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>docid</th>\n",
       "      <th>actualdate</th>\n",
       "      <th>inferreddate</th>\n",
       "      <th>firstpub</th>\n",
       "      <th>recept</th>\n",
       "      <th>recordid</th>\n",
       "      <th>OCLC</th>\n",
       "      <th>author</th>\n",
       "      <th>imprint</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>loc.ark+=13960=t8sb4zz1q</td>\n",
       "      <td>1921</td>\n",
       "      <td>1921</td>\n",
       "      <td>1921</td>\n",
       "      <td>addcanon</td>\n",
       "      <td>537314.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Lawrence, D. H.</td>\n",
       "      <td>New York;T. Seltzer;1921.</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>uc1.b3342759</td>\n",
       "      <td>1919</td>\n",
       "      <td>1919</td>\n",
       "      <td>1919</td>\n",
       "      <td>random</td>\n",
       "      <td>7930862.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Wigren, Bessie C.</td>\n",
       "      <td>Boston;The Poet Lore Company;c1919</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>uc1.b4100590</td>\n",
       "      <td>1918</td>\n",
       "      <td>1918</td>\n",
       "      <td>1918</td>\n",
       "      <td>reviewed</td>\n",
       "      <td>6154122.0</td>\n",
       "      <td>2143179.0</td>\n",
       "      <td>Waugh, Alec,</td>\n",
       "      <td>London;G. Richards;1918.</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>uc1.b3340220</td>\n",
       "      <td>1918</td>\n",
       "      <td>1918</td>\n",
       "      <td>1918</td>\n",
       "      <td>reviewed</td>\n",
       "      <td>7917249.0</td>\n",
       "      <td>12688503.0</td>\n",
       "      <td>Nightingale, M.</td>\n",
       "      <td>Oxford [Oxfordshire;B.H. Blackwell;1918.</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>uc2.ark+=13960=t0ft8gj1k</td>\n",
       "      <td>1918</td>\n",
       "      <td>1918</td>\n",
       "      <td>1918</td>\n",
       "      <td>reviewed</td>\n",
       "      <td>7657411.0</td>\n",
       "      <td>2518108.0</td>\n",
       "      <td>Faber, Geoffrey,</td>\n",
       "      <td>Oxford;B. H. Blackwell;New York;Longmans, Gree...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                      docid actualdate  inferreddate  firstpub    recept  \\\n",
       "0  loc.ark+=13960=t8sb4zz1q       1921          1921      1921  addcanon   \n",
       "1              uc1.b3342759       1919          1919      1919    random   \n",
       "2              uc1.b4100590       1918          1918      1918  reviewed   \n",
       "3              uc1.b3340220       1918          1918      1918  reviewed   \n",
       "4  uc2.ark+=13960=t0ft8gj1k       1918          1918      1918  reviewed   \n",
       "\n",
       "    recordid        OCLC             author  \\\n",
       "0   537314.0         NaN    Lawrence, D. H.   \n",
       "1  7930862.0         NaN  Wigren, Bessie C.   \n",
       "2  6154122.0   2143179.0       Waugh, Alec,   \n",
       "3  7917249.0  12688503.0    Nightingale, M.   \n",
       "4  7657411.0   2518108.0   Faber, Geoffrey,   \n",
       "\n",
       "                                             imprint  \n",
       "0                          New York;T. Seltzer;1921.  \n",
       "1                 Boston;The Poet Lore Company;c1919  \n",
       "2                           London;G. Richards;1918.  \n",
       "3           Oxford [Oxfordshire;B.H. Blackwell;1918.  \n",
       "4  Oxford;B. H. Blackwell;New York;Longmans, Gree...  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "### DEFACTORING INSPECTION\n",
    "metadata_file = pd.read_csv(classpath)\n",
    "#print(metadata_file.shape)\n",
    "print(\"The metadata files has {} rows and {} columns.\".format(\n",
    "    *metadata_file.shape))\n",
    "metadata_file.iloc[0:5,0:9] # display first 5 rows and 10 columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using *defactoring inspection* we can actually look at the metadata file and inspect the first five rows of metadata. By blending the original code with our inspection code and narrative the metadata becomes less of a conceptual abstraction and more of a tangible, material object that we can interrogate. Here we can see the file has 728 rows and 22 columns as well as the contents of the first five rows.\n",
    "\n",
    "> We worked with the HathiTrust Digital Library, which contains the aggregated collections of large public and university libraries: for 1820-1919, that gave us a collection of roughly 758,400 books in English, of which 53,200 include significant amounts of poetry\\...We gathered 360 reviewed and 360 random volumes, distributed in a similar way over the timeline (Underwood & Sellers 2016).\n",
    "\n",
    "This collection we have loaded into memory contains the metadata for the 720 poetry volumes out of the total of 53,200 identified poetry volumes in the HathiTrust Digital Library (Underwood 2014). What appears as a series of distance numbers in the MLQ article (Underwood & Sellers 2016, 325) are a tangible Python data structure capable of being directly inspected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "poemeta.csv\n",
      "DEFACTORING: Excluding volume with id loc.ark:/13960/t8sb4zz1q\n",
      "DEFACTORING: Excluding volume with id mdp.39015013402501\n",
      "DEFACTORING: Excluding volume with id mdp.39015011913525\n",
      "DEFACTORING: Excluding volume with id hardywessexpoems189.hardywessexpoems1898\n",
      "DEFACTORING: Excluding volume with id gerardmhopkins191.gerardmhopkins1918\n",
      "DEFACTORING: Excluding volume with id loc.ark:/13960/t3fx82c2q\n",
      "DEFACTORING: Excluding volume with id emilydickinso.emilydickinson\n",
      "DEFACTORING: Excluding volume with id ellisbell184.ellisbell1848\n",
      "We have 8 volumes in missing in metadata, and\n",
      "0 volumes missing in the directory.\n",
      "set()\n"
     ]
    }
   ],
   "source": [
    "### DEFACTORING FUNCTION \n",
    "### def get_metadata(classpath, volumeIDs, excludeif, excludeifnot, excludebelow, excludeabove):\n",
    "'''\n",
    "As the name would imply, this gets metadata matching a given set of volume\n",
    "IDs. It returns a dictionary containing only those volumes that were present\n",
    "both in metadata and in the data folder.\n",
    "\n",
    "It also accepts four dictionaries containing criteria that will exclude volumes\n",
    "from the modeling process.\n",
    "'''\n",
    "print(classpath)\n",
    "metadict = dict()\n",
    "\n",
    "with open(classpath, encoding = 'utf-8') as f:\n",
    "    reader = csv.DictReader(f)\n",
    "\n",
    "    anonctr = 0\n",
    "\n",
    "    for row in reader:\n",
    "        volid = dirty_pairtree(row['docid'])\n",
    "        theclass = row['recept'].strip()\n",
    "\n",
    "        # I've put 'remove' in the reception column for certain\n",
    "        # things that are anomalous.\n",
    "        if theclass == 'remove':\n",
    "            continue\n",
    "\n",
    "        bail = False\n",
    "        for key, value in excludeif.items():\n",
    "            if row[key] == value:\n",
    "                bail = True\n",
    "        for key, value in excludeifnot.items():\n",
    "            if row[key] != value:\n",
    "                bail = True\n",
    "        for key, value in excludebelow.items():\n",
    "            if forceint(row[key]) < value:\n",
    "                bail = True\n",
    "        for key, value in excludeabove.items():\n",
    "            if forceint(row[key]) > value:\n",
    "                bail = True\n",
    "\n",
    "        if bail:\n",
    "            print(\"DEFACTORING: Excluding volume with id \"+volid) ### DEFACTORING CODE\n",
    "            continue\n",
    "\n",
    "        birthdate = forceint(row['birth'])\n",
    "\n",
    "        pubdate = forceint(row['inferreddate'])\n",
    "\n",
    "        gender = row['gender'].rstrip()\n",
    "        nation = row['nationality'].rstrip()\n",
    "\n",
    "        #if pubdate >= 1880:\n",
    "            #continue\n",
    "\n",
    "        if nation == 'ca':\n",
    "            nation = 'us'\n",
    "        elif nation == 'ir':\n",
    "            nation = 'uk'\n",
    "        # I hope none of my Canadian or Irish friends notice this.\n",
    "\n",
    "        notes = row['notes'].lower()\n",
    "        author = row['author']\n",
    "        if len(author) < 1 or author == '<blank>':\n",
    "            author = \"anonymous\" + str(anonctr)\n",
    "            anonctr += 1\n",
    "\n",
    "        title = row['title']\n",
    "        canon = row['canon']\n",
    "\n",
    "        # I'm creating two distinct columns to indicate kinds of\n",
    "        # literary distinction. The reviewed column is based purely\n",
    "        # on the question of whether this work was in fact in our\n",
    "        # sample of contemporaneous reviews. The obscure column incorporates\n",
    "        # information from post-hoc biographies, which trumps\n",
    "        # the question of reviewing when they conflict.\n",
    "\n",
    "        if theclass == 'random':\n",
    "            obscure = 'obscure'\n",
    "            reviewed = 'not'\n",
    "        elif theclass == 'reviewed':\n",
    "            obscure = 'known'\n",
    "            reviewed = 'rev'\n",
    "        elif theclass == 'addcanon':\n",
    "            print(\"DEFACTORING: adding volume\") ### DEFACTORING CODE\n",
    "            obscure = 'known'\n",
    "            reviewed = 'addedbecausecanon'\n",
    "        else:\n",
    "            print(\"Missing class\" + theclass)\n",
    "\n",
    "        if notes == 'well-known':\n",
    "            obscure = 'known'\n",
    "        if notes == 'obscure':\n",
    "            obscure = 'obscure'\n",
    "\n",
    "        if canon == 'y':\n",
    "            if theclass == 'addcanon':\n",
    "                actually = 'Norton, added'\n",
    "            else:\n",
    "                actually = 'Norton, in-set'\n",
    "        elif reviewed == 'rev':\n",
    "            actually = 'reviewed'\n",
    "        else:\n",
    "            actually = 'random'\n",
    "\n",
    "        metadict[volid] = dict()\n",
    "        metadict[volid]['reviewed'] = reviewed\n",
    "        metadict[volid]['obscure'] = obscure\n",
    "        metadict[volid]['pubdate'] = pubdate\n",
    "        metadict[volid]['birthdate'] = birthdate\n",
    "        metadict[volid]['gender'] = gender\n",
    "        metadict[volid]['nation'] = nation\n",
    "        metadict[volid]['author'] = author\n",
    "        metadict[volid]['title'] = title\n",
    "        metadict[volid]['canonicity'] = actually\n",
    "        metadict[volid]['pubname'] = row['pubname']\n",
    "        metadict[volid]['firstpub'] = forceint(row['firstpub'])\n",
    "\n",
    "# These come in as dirty pairtree; we need to make them clean.\n",
    "\n",
    "cleanmetadict = dict()\n",
    "allidsinmeta = set([x for x in metadict.keys()])\n",
    "allidsindir = set([dirty_pairtree(x) for x in volumeIDs])\n",
    "missinginmeta = len(allidsindir - allidsinmeta)\n",
    "missingindir = len(allidsinmeta - allidsindir)\n",
    "print(\"We have \" \n",
    "      + str(missinginmeta) \n",
    "      + \" volumes in missing in metadata, and\")\n",
    "print(str(missingindir) + \" volumes missing in the directory.\")\n",
    "print(allidsinmeta - allidsindir)\n",
    "\n",
    "for anid in volumeIDs:\n",
    "    dirtyid = dirty_pairtree(anid)\n",
    "    if dirtyid in metadict:\n",
    "        cleanmetadict[anid] = metadict[dirtyid]\n",
    "\n",
    "# Now that we have a list of volumes with metadata, we can select the groups of IDs\n",
    "# that we actually intend to contrast. If we want to us more or less everything,\n",
    "# this may not be necessary. But in some cases we want to use randomly sampled subsets.\n",
    "\n",
    "# The default condition here is\n",
    "\n",
    "# category2sorton = 'reviewed'\n",
    "# positive_class = 'rev'\n",
    "# sizecap = 350\n",
    "# A sizecap less than one means, no sizecap.\n",
    "\n",
    "### DEFACTORING FUNCTION CALL\n",
    "### IDsToUse, classdictionary = metafilter.label_classes(metadict, category2sorton, positive_class, sizecap)\n",
    "\n",
    "### DEFACTORING NAMESPACE \n",
    "metadict = cleanmetadict  # put the data into the global namespace so execution can continue."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This above code cell is large due to a long `for` loop processing each row of the metadata file. At a high level, the code in this cell loads the metadata and determines which volumes to exclude in the analysis. It does this by loading the `poemeta.csv` file and excluding rows based upon the parameters specified in the `excludeif`, `excludeifnot`, `excludeabove`, and `excludebelow` variables. This process removed 8 volume designations from the (meta)data , which explains why there were 728 volumes in the metadata, but 720 volumes mentioned in the MLQ article. The resulting output immediately above is a mixture of the author's code and our own DEFACTORING inspection statements (marked with the comment \\#\\#\\# DEFACTORING). We have added a print statement so we can see the IDs of the volumes being excluded in the code.\n",
    "\n",
    "Beyond filtering out excluded (meta)data, this code also makes a series of normalizing decisions, that is, there are more conceptual distinctions being made (or unmade) in this code. First is the normalization of nationality, which is a clinical way of saying Underwood and Sellers lump Canada with the United States and Ireland with the UK. Nationality was briefly explored in the MLQ analysis, and it interesting to see the residue code here. However, the main thrust of the MLQ article focuses on word frequencies as features and leaves extant features such as gender and nationality as future directions to be explored more fully (Underwood & Sellers 2016).\n",
    "\n",
    "More importantly, this code cell splits the recept column of the metadata file into two columns, obscure and reviewed. These are important distinctions for the analysis because this is where the labels for training the computational model are specified. From the code and the comments, there are poems that were reviewed, and there are poems that were obscure. Lastly, the logic of the code indicates there are poems that are not in the set of reviewed poems but are nevertheless part of the canon. In the latter case the poems are set to be \"known\" `(obscure = \\'known\\')`. According to the author's comment this trumps the conflict when the author is known but not explicitly in the reviewed set.\n",
    "\n",
    "> To confirm that a given book had never been reviewed in any of \\[our selected\\] publications would have been a tedious task. It was more straightforward simply to select works at random from a very large collection while excluding authors already in our reviewed sample: in practice, this turned up mostly books that were rarely reviewed (Underwood & Sellers 2016).\n",
    "\n",
    "This is not the code that selects the sample of unreviewed works, that work had been done previously and was encoded in the metadata. Rather, this code is loading and preparing the data based upon that previous work while also \"removing authors already in our reviewed sample.\" Here we can see a connection between the prosaic description of methodology in the MLQ article with that same method made material in the code.\n",
    "\n",
    "What might be important to note here is how \"great debates\" in literary history about literary prestige, obscurity, and the canon are being ascribed in the code without much fanfare. There is a hard decision being made (in the blink of an eye) about the status of particular literary works. Most of these are, we suspect, fairly uncontroversial distinctions that accord with the broader community, but the code and computation enforce clear and unambiguous decisions for each and every volume. These hard decisions are pragmatically executed to get to more interesting analyses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "OCLC                                         NaN\n",
       "actualdate                                  1860\n",
       "author                     Drake, Joseph Rodman,\n",
       "birth                                       1795\n",
       "canon                                        NaN\n",
       "docid                             wu.89099921512\n",
       "enumcron                                     NaN\n",
       "firstpub                                    1860\n",
       "gender                                         m\n",
       "impaud                                       NaN\n",
       "imprint         New York;The Bradford Club;1860.\n",
       "inferreddate                                1860\n",
       "judge                                        NaN\n",
       "nationality                                   us\n",
       "notes                                        NaN\n",
       "othername                                    NaN\n",
       "pubname                                      NaN\n",
       "pubrev                                       NaN\n",
       "recept                                    random\n",
       "recordid                               8998599.0\n",
       "title                               The croakers\n",
       "yrrev                                        NaN\n",
       "Name: 534, dtype: object"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "### DEFACTORING INSPECTION\n",
    "# Examine the original metadata file \n",
    "\n",
    "defactoring_volume_id = 'wu.89099921512'\n",
    "\n",
    "the_croakers_metadata = (\n",
    "    metadata_file\n",
    "    .loc[metadata_file['docid'] == defactoring_volume_id]\n",
    "    .squeeze()\n",
    "    .sort_index()\n",
    ")\n",
    "the_croakers_metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'reviewed': 'not',\n",
       " 'obscure': 'obscure',\n",
       " 'pubdate': 1860,\n",
       " 'birthdate': 1795,\n",
       " 'gender': 'm',\n",
       " 'nation': 'us',\n",
       " 'author': 'Drake, Joseph Rodman,',\n",
       " 'title': 'The croakers',\n",
       " 'canonicity': 'random',\n",
       " 'pubname': '',\n",
       " 'firstpub': 1860}"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "### DEFACTORING INSPECTION\n",
    "# Examine the cleaned metadata\n",
    "cleanmetadict[defactoring_volume_id]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The inspection above shows the data expressed in the CSV file has been transformed into into a Python dictionary, `cleanmetadict`, with additional columns for expressing more granularity about the categorizations for each poetry volume. We also observe the raw metadata csv file has additional columns that are not reflected in the Python dictionary. What we see reflected in `cleanmetadict` is only the metadata necessary for the analysis with any dirty or unnecessary information removed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Sorting Training Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "### DEFACTORING FUNCTION\n",
    "### def label_classes(metadict, category2sorton, positive_class, sizecap):\n",
    "''' This takes as input the metadata dictionary generated\n",
    "by get_metadata. It subsets that dictionary into a\n",
    "positive class and a negative class. Instances that belong\n",
    "to neither class get ignored.\n",
    "'''\n",
    "\n",
    "all_instances = set([x for x in metadict.keys()])\n",
    "\n",
    "# The first stage is to find positive instances.\n",
    "\n",
    "all_positives = set()\n",
    "\n",
    "for key, value in metadict.items():\n",
    "    if value[category2sorton] == positive_class:\n",
    "        all_positives.add(key)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This block of code reads the metadata properties and puts a subset of all entries into a variable, `all_positives`, which will contain all of the volume ids for reviewed poems. If the `reviewed` column has a value of 'rev', then it is selected for inclusion. The name and value of the property are parameterized however, so technically it is more correct, but more opaque as well, to state: if a poem's metadata has the value 'rev' (specified by the `positive_class` variable) for the reviewed property (specified by the `category2sorton` variable) then it is labeled as a positive. Having thus collected all the reviewed poems into the set named `all_positives`, the next cell populates the variable `all_negatives` with all the instances not in the positive set by subtracting the set of positives from the set of all instances by applying a basic mathematical set operation (-)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_negatives = all_instances - all_positives\n",
    "iterator = list(all_negatives)\n",
    "for item in iterator:\n",
    "    if metadict[item]['reviewed'] == 'addedbecausecanon':\n",
    "        all_negatives.remove(item)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The negative labels are assigned to all instances that are not in the set of positive instances. This section includes additional code that filters out any item with `addedbecausecannon` set for the `reviewed` property, but this code should never execute because, as we have seen above, the volumes in the canon have already been removed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "360\n"
     ]
    }
   ],
   "source": [
    "if sizecap > 0 and len(all_positives) > sizecap:\n",
    "    positives = random.sample(all_positives, sizecap)\n",
    "else:\n",
    "    positives = list(all_positives)\n",
    "    print(len(all_positives))\n",
    "\n",
    "# If there's a sizecap we also want to ensure classes have\n",
    "# matching sizes and roughly equal distributions over time.\n",
    "\n",
    "numpositives = len(all_positives)\n",
    "\n",
    "if sizecap > 0 and len(all_negatives) > numpositives:\n",
    "    if not 'date' in category2sorton:\n",
    "        available_negatives = list(all_negatives)\n",
    "        negatives = list()\n",
    "\n",
    "        for anid in positives:\n",
    "            date = metadict[anid]['pubdate']\n",
    "\n",
    "            available_negatives = sort_by_proximity(available_negatives, \n",
    "                                                    metadict, date)\n",
    "            selected_id = available_negatives.pop(0)\n",
    "            negatives.append(selected_id)\n",
    "\n",
    "    else:\n",
    "        # if we're dividing classes by date, we obvs don't want to\n",
    "        # ensure equal distributions over time.\n",
    "\n",
    "        negatives = random.sample(all_negatives, sizecap)\n",
    "\n",
    "else:\n",
    "    negatives = list(all_negatives)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most of the code in the cell above does not actually execute because the number of entries in the `all_positives` and `all_negatives` lists are not greater than `sizecap`. The conditional statements on line 1 and line 12 will not be true so the the accompanying blocks of code never execute. If the `sizecap` variable was smaller, or the number of entries larger, this code would use random sampling to select smaller number of entries from the positives entries.\n",
    "\n",
    "Looking at the block of code for the negative entries is a bit more interesting. This block of code makes an unexecuted reference to a function `sort_by_proximity()` that samples from the negative elements with an equal distribution based upon some function of proximity. Because this code is not executing we are not going to spend more time and analytical attention to exactly how this function operates. Furthermore, we have not included the code for `sort_by_proximity()` in the notebook because it is not part of the execution path we are tracing. In the code's garden of forking paths, this is a path not taken.\n",
    "\n",
    "These issues of bespoke code that is not executed and functions that are not called point to properties of code that make it complex and therefore difficult to review or critique. Code has a textual and a processual dimension (Hiller 2015, Van Zundert 2016). The code-as-text is what we can see and read in the source code, the processual dimension of code is tied to its execution as a software program. Code critique moves in between these two modes of existence of code. We are here, as code critics, not simply looking at code-as-text. That is, we are in this case reviewing a *live* execution of the code. This is extremely significant and, we'd argue, distinguishes defactoring as more than analysing code-as-text; we are analyzing the code, the data, and their interaction in the computation. Leveraging the affordances of the Jupyter platform allows us the ability to interact with the methodology inscribed in the code. At each step of the incremental process we can ask questions by inspecting the state of variables (or even change them). This is more than simply treating the code as a text, the code is but one part of a complex assemblage we have been executing and inspecting. However, it is also not a complete inspection of all the ways in which the code can be possibly executed. As we *defactor* the Underwood and Sellers's code, we make choices about how much of the code to include for the argument we are trying to make (and for the sake of our time and attention). Thus we are dealing with a code-criticism conundrum: What is the required or adequate breadth and depth of the computational process and subsequently of the critique? The decision to include or not include `sort_by_proximity()` is a breadth issue. How broad should we be in including code that does not execute? Note that we are including code from a conditional block that does not execute, but are not going out the additional step to include non-executed functions defined elsewhere in the code. The decision to include or not include code from the standard library, code not written by the authors, is a depth issue. While there are many functions we are stepping over, like `len`, `list`, `append`, `pop`, `random.sample`, we argue there is no need to step into these functions because, following Hinsen's model, they are not part of the bespoke code of Pace of Change. Again, this raises the problematic issue of our decision to step over `sort_by_proximity()` even though it was written by the authors' for this particular project.\n",
    "\n",
    "Obviously the 'rules of the game' for defactoring are not yet clear established. Therefore we are more or less 'feeling' our way through an emerging methodology for code criticism. As we see vestiges of the authors' evolution in thinking in their code, this notebook is capturing the evolution of our thinking about defactoring as a practice.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "We have 360 positive, and\n",
      "360 negative instances.\n"
     ]
    }
   ],
   "source": [
    "# Now we have two lists of ids.\n",
    "\n",
    "IDsToUse = set()\n",
    "classdictionary = dict()\n",
    "\n",
    "print()\n",
    "print(\"We have \" + str(len(positives)) + \" positive, and\")\n",
    "print(str(len(negatives)) + \" negative instances.\")\n",
    "\n",
    "for anid in positives:\n",
    "    IDsToUse.add(anid)\n",
    "    classdictionary[anid] = 1\n",
    "\n",
    "for anid in negatives:\n",
    "    IDsToUse.add(anid)\n",
    "    classdictionary[anid] = 0\n",
    "\n",
    "for key, value in metadict.items():\n",
    "    if value['reviewed'] == 'addedbecausecanon':\n",
    "        print(\"DEFACTORING: Adding cannon supplement\") ### DEFACTORING CODE\n",
    "        IDsToUse.add(key)\n",
    "        classdictionary[key] = 0\n",
    "# We add the canon supplement, but don't train on it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this cell we are seeing yet another instance of metadata being shaped and transformed in preparation for analysis. The code first prints out the number of positive and negative instances by checking the length (using `len()`) of the volume ids stored in the `positives` and `negatives` variables. Two loops iterate over these lists and populate two more variables, `IDsToUse` and `classdictionary`. The first, `IDsToUse` contains a master list of all the volume identifiers to be used in the analysis. It is of the Python set datatype, meaning there will be no duplicate identifiers in the set list. The second, `classdictionary` is a dictionary that allows a simple lookup to see if a volume ID is in the positive or negative class–as indicated by a 0 or a 1. There is a final loop whose logic checks to see if any volumes have a specific metadata flag reviewed with a value of `addedbecausecanon`. We have added a defactoring statement to see if this logic is ever triggered. The output indicates the `if` statement's conditions were never satisfied, no volumes were added because of cannon.\n",
    "\n",
    "We have come to the end of the preparing metadata section. All of the code up to this point has focused on loading, normalizing, and transforming the metadata--namely the identifiers of the volumes to be analyzed. Based upon the values in the metadata fields and assumptions built into the logic of the code, the authors have assembled the list of volume ids and their associated binary labels: reviewed or random, 0 or 1. Because this is a *supervised* machine learning exercise, the authors need labeled data to train the model. All of the work in this section of the code was dedicated to assigning a class label (positive or negative) to the identifiers of the data files. The next section dives into the actual data itself."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transforming Words into Features\n",
    "\n",
    "Now that we know exactly which volumes of poetry the code will be analyzing, we can venture into the datafiles and begin the work of transforming the volume data files into a data structure suitable for analysis. The logistic regression requires the data to be in a specific shape, a matrix of binary features. This section does the work of *getting the data into shape*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "### DEFACTORING FUNCTION DEFINITIONS\n",
    "### We need to define the infer_date function\n",
    "\n",
    "def infer_date(metadictentry, datetype):\n",
    "    if datetype == 'pubdate':\n",
    "        return metadictentry[datetype]\n",
    "    elif datetype == 'firstpub':\n",
    "        firstpub = metadictentry['firstpub']\n",
    "        if firstpub > 1700 and firstpub < 1950:\n",
    "            return firstpub\n",
    "        else:\n",
    "            return metadictentry['pubdate']\n",
    "    else:\n",
    "        sys.exit(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This code cell defines a helper function, `infer_date()`, which is used in the code below to deal with differences in the `pubdate` and `firstpub` columns in the metadata. When `firstpub` falls between 1700 and 1950 the codes uses that as the date, otherwise it returns the value in `pubdate` (or it exists the script in the case of bad data)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# make a vocabulary list and a volsize dict\n",
    "wordcounts = Counter()\n",
    "\n",
    "volspresent = list()\n",
    "orderedIDs = list()\n",
    "\n",
    "positivecounts = dict()\n",
    "negativecounts = dict()\n",
    "\n",
    "for volid, volpath in zip(volumeIDs, volumepaths):\n",
    "    if volid not in IDsToUse:\n",
    "        continue\n",
    "    else:\n",
    "        volspresent.append((volid, volpath))\n",
    "        orderedIDs.append(volid)\n",
    "\n",
    "    date = infer_date(metadict[volid], datetype)\n",
    "    if date < pastthreshold or date > futurethreshold:\n",
    "        continue\n",
    "    else:\n",
    "        with open(volpath, encoding = 'utf-8') as f:\n",
    "            for line in f:\n",
    "                fields = line.strip().split('\\t')\n",
    "                if len(fields) > 2 or len(fields) < 2:\n",
    "                    # print(line)\n",
    "                    continue\n",
    "                word = fields[0]\n",
    "                if len(word) > 0 and word[0].isalpha():\n",
    "                    count = int(fields[1])\n",
    "                    wordcounts[word] += 1\n",
    "                    # for initial feature selection we use the number of\n",
    "                    # *documents* that contain a given word,\n",
    "                    # so it's just +=1.\n",
    "\n",
    "vocablist = [x[0] for x in wordcounts.most_common(numfeatures)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is an important section because it contains the code that opens the data files and selects the word-features used in the logistic regression. The main block of code in the cell above loops over each data file (representing a poem) and counts the number of instances of each word. Like the metadata file, we can use defactoring to insThis is an important section because it contains the code that opens the data files and selects the word-features used in the logistic regression. The main block of code in the cell above loops over each data file (representing a poem) and counts the number of instances of each word. This code performs the extremely important work of identifying the 3,200 word features used in the analysis. Underwood & Sellers only briefly mention this choice:\n",
    "\n",
    "The model we have trained uses thirty-two hundred variables---the frequencies of the thirty-two hundred words most common in the collection. A model this complex can encode a lot of information about a social boundary. But the complexity also creates room for a definition of literary prestige that can be about a lot of intersecting things at once; it is not required to map onto any single idea. Moreover, we may not be able to characterize the effect of any single word with great precision, since variables interact with each other in tricky ways. A list of the top ten words that individually have the largest effect on the model's predictions might not tell us very much (Underwood & Sellers 2016, 330).\n",
    "\n",
    "The decision to use 3,200 word features is never discussed in the MLQ article, but is rather situated in a broader discussion about the significance of individual words and how they can be marshaled as evidence of distinctions found via a statistical model. The 3,200 words selected here will eventually be transformed into a set of coefficients below with which the boundaries of statistical significance can be observed. The lofty philosophical discussion of the MLQ article is rooted in mundane materiality of a Pythonic function call.pect a datafile to get a sense of what these data look like before they are transformed by the code. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Word</th>\n",
       "      <th>Count</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>,</td>\n",
       "      <td>1453</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>the</td>\n",
       "      <td>1089</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>and</td>\n",
       "      <td>830</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>of</td>\n",
       "      <td>513</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>.</td>\n",
       "      <td>357</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>to</td>\n",
       "      <td>320</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>in</td>\n",
       "      <td>289</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>!</td>\n",
       "      <td>277</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>a</td>\n",
       "      <td>250</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>|'s|</td>\n",
       "      <td>229</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Word  Count\n",
       "0     ,   1453\n",
       "1   the   1089\n",
       "2   and    830\n",
       "3    of    513\n",
       "4     .    357\n",
       "5    to    320\n",
       "6    in    289\n",
       "7     !    277\n",
       "8     a    250\n",
       "9  |'s|    229"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "### DEFACTORING INSPECTION\n",
    "data_file = pd.read_csv(volumepaths[0], \n",
    "                        delimiter='\\t', \n",
    "                        names=[\"Word\", \"Count\"])\n",
    "data_file.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output above shows the first 10 lines of one of the poem data files. As we can plainly see, the volume has already been pre-processed into a list of words and their frequencies. This particular volume has 2,745 commas and 1445 instances of the word “the.” The authors’ code parses each of these files and assembles a vocabulary list of the 3200 most common words (as specified by the `numfeatures` variable)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Word      Count\n",
      "---------------\n",
      "the         720\n",
      "and         720\n",
      "of          720\n",
      "to          720\n",
      "in          720\n",
      "a           720\n",
      "that        720\n",
      "with        720\n",
      "all         720\n",
      "as          720\n",
      "on          720\n",
      "for         720\n",
      "is          720\n",
      "i           720\n",
      "but         720\n",
      "by          720\n",
      "at          720\n",
      "not         720\n",
      "his         719\n",
      "from        719\n"
     ]
    }
   ],
   "source": [
    "### DEFACTORING INSPECTION\n",
    "print(\"Word      Count\")\n",
    "print(\"---------------\")\n",
    "for word, count in wordcounts.most_common(n=20):\n",
    "    print(\"{:8}  {:5}\".format(word, count))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At first glance it might seem strange that the count is 720 for all of the top words in the corpus. However, when we dig deeper into the code we can see that the authors are not tabulating the total word frequencies across all volumes in the corpus, rather they are associating words and the number of volumes; a distinction not explicitly expressed in the MLQ article. The code loops over each file, opening it, and parses each line by splitting on the tab character (\"\\\\t\"). What is interesting is that Underwood and Sellers are only paying attention to the word and ignoring the frequency within the volume. They check to see if the word is longer than zero and use the `isalpha()` function to make sure the characters are alphabetic as opposed to punctuation. The comments in the code explain that the authors are just using the \"number of documents that contain a given word\".\n",
    "\n",
    "They are selecting their list of features (stored in the `vocablist` variable) by selecting words ranked by the number of documents in which they appear. The total number of documents we are working with is 720, so the table we generated above tells us that the top ten words appear in all of the documents. If we look at more than just the top ten, we can start to see the distribution of words in documents."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfcAAAFACAYAAAChujXqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAi80lEQVR4nO3df3ST9fn/8VfTULHUlrQphRYUC+0AraDCBASKkM0zddoPE9gYOsThBF0nbB4QN9TD1G5aimixjiGe4dHNndlujs0dY0fZ7NwKiCKKioLC4Ud/JLQU+oMm+f7B1xwrqaRN2yTvPh/neI55505y9SLlxXXfyX3H+Hw+nwAAgDEs4S4AAAB0L8IdAADDEO4AABiGcAcAwDCEOwAAhiHcAQAwDOEOAIBhrOEuoDsdPnz4rDW73a7a2towVGMOehg6ehg6ehga+he6SOxhenp6wHUmdwAADEO4AwBgGMIdAADDEO4AABiGcAcAwDCEOwAAhiHcAQAwDOEOAIBhCHcAAAzTK2eoW79+vXbu3KmkpCQVFhZKkjZv3qwdO3bIarUqLS1NS5Ys0YABAyRJpaWlKi8vl8Vi0W233aZx48b1RpkAABihVyb36dOna+XKle3WLrvsMhUWFurxxx/XkCFDVFpaKkk6dOiQKisrtWbNGt1///3auHGjvF5vb5QJAIARemVyHzNmjKqrq9utjR071v//2dnZevPNNyVJVVVVmjx5svr166dBgwZp8ODB2rdvn7Kzs3uj1HPyLLox4Hrshr/0ciUAAAQWEcfcy8vL/bveXS6XUlJS/PclJyfL5XKFqTIAAKJP2K8K9/LLLys2NlZTp06VJPl8vqAf63Q65XQ6JUkFBQWy2+1nbWO1WgOud9WxDta78zUiTXf3sC+ih6Gjh6Ghf6GLph6GNdy3bt2qHTt2aNWqVYqJiZEkpaSkqK6uzr+Ny+VScnJywMc7HA45HA7/7UCX4uutS/RF2mUAu1MkXuYw2tDD0NHD0NC/0EViDyPukq+7du3Sn//8Zy1fvlznnXeef338+PGqrKzU6dOnVV1drSNHjmjkyJHhKhMAgKjTK5P72rVr9d577+nEiRO68847NWfOHJWWlqqtrU2rV6+WJGVlZemOO+7QsGHDNGnSJC1btkwWi0W33367LJaI+GgAAABRoVfC/Z577jlrbcaMGR1uP2vWLM2aNasHKwIAwFyMxAAAGIZwBwDAMIQ7AACGIdwBADBM2E9iYwpOSwsAiBRM7gAAGIZwBwDAMIQ7AACGIdwBADAM4Q4AgGEIdwAADEO4AwBgGMIdAADDEO4AABiGcAcAwDCEOwAAhiHcAQAwDOEOAIBhCHcAAAxDuAMAYBjCHQAAwxDuAAAYhnAHAMAwhDsAAIYh3AEAMIw13AVEMs+iG8NdAgAAncbkDgCAYQh3AAAMQ7gDAGAYwh0AAMMQ7gAAGIZwBwDAMIQ7AACGIdwBADBMr5zEZv369dq5c6eSkpJUWFgoSWpsbFRRUZFqamqUmpqqpUuXKiEhQZJUWlqq8vJyWSwW3XbbbRo3blxvlNlrOjo5TuyGv/RyJQAAE/XK5D59+nStXLmy3VpZWZlycnK0bt065eTkqKysTJJ06NAhVVZWas2aNbr//vu1ceNGeb3e3igTAAAj9Eq4jxkzxj+Vf66qqkq5ubmSpNzcXFVVVfnXJ0+erH79+mnQoEEaPHiw9u3b1xtlAgBghLAdc6+vr5fNZpMk2Ww2NTQ0SJJcLpdSUlL82yUnJ8vlcoWlRgAAolHEXTjG5/MFva3T6ZTT6ZQkFRQUyG63n7WN1WoNuB6MY116VHuBXruj5+1qnT0tlB7iDHoYOnoYGvoXumjqYdjCPSkpSW63WzabTW63W4mJiZKklJQU1dXV+bdzuVxKTk4O+BwOh0MOh8N/u7a29qxt7HZ7wPXe0pnXDmedXyXcPTQBPQwdPQwN/QtdJPYwPT094HrYdsuPHz9eFRUVkqSKigpNmDDBv15ZWanTp0+rurpaR44c0ciRI8NVJgAAUadXJve1a9fqvffe04kTJ3TnnXdqzpw5ysvLU1FRkcrLy2W327Vs2TJJ0rBhwzRp0iQtW7ZMFotFt99+uywWvo4PAECweiXc77nnnoDrq1atCrg+a9YszZo1qwcrAgDAXIzEAAAYhnAHAMAwhDsAAIYh3AEAMAzhDgCAYQh3AAAMQ7gDAGAYwh0AAMMQ7gAAGIZwBwDAMBF3yVfTeBbdGO4SAAB9DJM7AACGIdwBADAM4Q4AgGEIdwAADEO4AwBgGMIdAADDEO4AABiGcAcAwDCEOwAAhiHcAQAwDOEOAIBhCHcAAAxDuAMAYBjCHQAAwxDuAAAYhnAHAMAw1nAXgK7xLLox4Hrshr/0ciUAgEjD5A4AgGGY3CNIoGmcSRwA0FlM7gAAGIZwBwDAMIQ7AACG4Zi7YSLluH2k1AEAfVHYw/2vf/2rysvLFRMTo2HDhmnJkiVqbW1VUVGRampqlJqaqqVLlyohISHcpQIAEBXCulve5XLp73//uwoKClRYWCiv16vKykqVlZUpJydH69atU05OjsrKysJZJgAAUSXsx9y9Xq9aW1vl8XjU2toqm82mqqoq5ebmSpJyc3NVVVUV5ioBAIgeYd0tn5ycrG9/+9tavHix4uLiNHbsWI0dO1b19fWy2WySJJvNpoaGhnCWCQBAVAlruDc2NqqqqkrFxcWKj4/XmjVrtG3btqAf73Q65XQ6JUkFBQWy2+1nbWO1WgOuB+NYlx7VvTqqvTO1dfXn/1xXehiovlDriGahvA9xBj0MDf0LXTT1MKzhvnv3bg0aNEiJiYmSpKuuukoffvihkpKS5Ha7ZbPZ5Ha7/fd/mcPhkMPh8N+ura09axu73R5wPVp0R+2hPkd39TCa/xxCFe3vw0hAD0ND/0IXiT1MT08PuB7WY+52u10fffSRWlpa5PP5tHv3bmVkZGj8+PGqqKiQJFVUVGjChAnhLBMAgKgS1sk9KytLEydO1PLlyxUbG6vhw4fL4XCoublZRUVFKi8vl91u17Jly8JZJgAAUSXs33OfM2eO5syZ026tX79+WrVqVZgqAgAguoX9q3AAAKB7Ee4AABiGcAcAwDBhP+aOrxboAiwAAHwVJncAAAxDuAMAYBjCHQAAwxDuAAAYhnAHAMAwhDsAAIYh3AEAMAzhDgCAYQh3AAAMQ7gDAGCYoMP9P//5T8D1N998s9uKAQAAoQs63EtKSgKuP/PMM91WDAAACN05Lxxz7NgxSZLX61V1dbV8Pl+7++Li4nquOgAA0GnnDPf8/Hz////4xz9ud9/AgQM1e/bs7q8KAAB02TnD/Q9/+IMk6YEHHtBDDz3U4wUBAIDQBH3MnWAHACA6nHNy/1x1dbVefPFFHThwQM3Nze3ue/rpp7u9MAAA0DVBh/sTTzyhtLQ03XrrrTrvvPN6siYAABCCoMP90KFDWr16tSwWznuDrvEsujHobWM3/KUHKwEAswWd1KNHj9aBAwd6sBQAANAdgp7cU1NT9fDDD+vrX/+6Bg4c2O6+uXPndnddAACgi4IO95aWFl155ZXyeDyqq6vryZoQZl/efX5MHe8m78yudgBA7wg63JcsWdKTdQAAgG4SdLh/fhraQNLS0rqlGAAAELqgw/2Lp6H9ss/PYgcAAMIv6HD/coAfP35cf/zjHzV69OhuLwqRh2PrABA9uvyl9YEDB2rBggV64YUXurMeAAAQoqAn90AOHz6slpaW7qoF8OtoTwEntwGAcws63FetWqWYmBj/7ZaWFh08eFA333xzjxQGAAC6JuhwnzFjRrvb/fv310UXXaQhQ4Z0e1HoXkzBANC3BB3u06dP75ECTp48qZKSEh08eFAxMTFavHix0tPTVVRUpJqaGqWmpmrp0qVKSEjokdcHAMA0QYd7W1ubXn75ZW3btk1ut1s2m03Tpk3TrFmzZLV2/dD9pk2bNG7cOP30pz9VW1ubWlpaVFpaqpycHOXl5amsrExlZWWaP39+l18DAIC+JOhPyz///PPavXu3Fi1apMcee0yLFi3Su+++q+eff77LL37q1Cm9//77/l3+VqtVAwYMUFVVlXJzcyVJubm5qqqq6vJrAADQ1wQ9cr/55pt67LHHdMEFF0iS0tPTdfHFF+vee+/VggULuvTi1dXVSkxM1Pr16/Xpp58qMzNTCxYsUH19vWw2myTJZrOpoaEh4OOdTqecTqckqaCgQHa7/axtrFZrwPVgdHxOPjN01JdI/rk782d57P8mB1xPK63srnKCFsr7EGfQw9DQv9BFUw+DDnefz9ftL+7xeLR//34tXLhQWVlZ2rRpk8rKyoJ+vMPhkMPh8N+ura09axu73R5wHYH7Fem6o+Zw/Ny8D0NHD0ND/0IXiT1MT08PuB70bvlJkybpV7/6lXbt2qVDhw5p165deuyxxzRx4sQuF5WSkqKUlBRlZWVJkiZOnKj9+/crKSlJbrdbkuR2u5WYmNjl1wAAoK8JenKfP3++/vSnP2njxo1yu91KTk7W1Vdfre985ztdfvGBAwcqJSVFhw8fVnp6unbv3q2hQ4dq6NChqqioUF5enioqKjRhwoQuvwbMwtf6AODczhnue/fu1fbt2zV//nzNnTtXc+fO9d/3/PPP65NPPlF2dnaXC1i4cKHWrVuntrY2DRo0SEuWLJHP51NRUZHKy8tlt9u1bNmyLj8/AAB9zTnDvbS0VNdee23A+y699FK9/PLLWrFiRZcLGD58uAoKCs5aX7VqVZefEwCAvuycx9wPHDigcePGBbwvJydH+/fv7+6aAABACM4Z7k1NTWprawt4n8fjUVNTU7cXBQAAuu6c4Z6RkaG333474H1vv/22MjIyur0oAADQdec85n799dfrN7/5jbxeryZMmCCLxSKv16uqqipt3LhRt956a2/UCQBA1Aj0zZ7e/FbPOcN9ypQpOn78uIqLi3X69GklJiaqoaFBcXFxmj17tqZMmdIbdQIAgCAF9T33G264QTNmzNCHH36oxsZGJSQkKDs7W/Hx8T1dHwAA6KSgT2ITHx/f4afmAQBA5Aj69LMAACA6EO4AABiGcAcAwDCEOwAAhiHcAQAwTNCflod5Orp8KgAgujG5AwBgGMIdAADDEO4AABiGY+7A/xfuCz0AQHdhcgcAwDBM7sBX6OgbBUz0ACIZkzsAAIZhcocRevs7+7090bMHAUBnMLkDAGAYwh0AAMMQ7gAAGIZj7uhz+uo59Ttz3J5j/EB0Y3IHAMAwhDsAAIZhtzwQJpzuFkBPYXIHAMAwhDsAAIYh3AEAMAzH3IE+rq9+NRAwGZM7AACGiYjJ3ev1asWKFUpOTtaKFSvU2NiooqIi1dTUKDU1VUuXLlVCQkK4ywR6HFM0gO4QEZP73/72N2VkZPhvl5WVKScnR+vWrVNOTo7KysrCVxwAAFEm7OFeV1ennTt3aubMmf61qqoq5ebmSpJyc3NVVVUVrvKAiOZZdONZ/wFA2MP9ueee0/z58xUTE+Nfq6+vl81mkyTZbDY1NDSEqzwAAKJOWI+579ixQ0lJScrMzNSePXs6/Xin0ymn0ylJKigokN1uP2sbq9UacD0Yx7r0KPRlHb3XAr0Pe+r91VEN3fF6Xf1d6g6h/C6D/nWHzvQw0O9bb/Y/rOH+wQcfaPv27XrrrbfU2tqqpqYmrVu3TklJSXK73bLZbHK73UpMTAz4eIfDIYfD4b9dW1t71jZ2uz3gOtATOnqv9eb7sCdfJ5y/S/wuh4b+hS7UHvZE/9PT0wOuhzXc582bp3nz5kmS9uzZo1deeUX5+fnavHmzKioqlJeXp4qKCk2YMCGcZQIAEFCkfs4l7MfcA8nLy9M777yj/Px8vfPOO8rLywt3SQAARI2I+J67JF1yySW65JJLJEkXXHCBVq1aFeaKAACIThE5uQMAgK4j3AEAMAzhDgCAYQh3AAAMEzEfqAPQPSL1qzkAeg+TOwAAhmFyB7pRoKk5dsNfwlAJgL6MyR0AAMMQ7gAAGIZwBwDAMBxzB3qYZ9GNXD4YQK9icgcAwDBM7gBCxrcEECk6Os9DX3s/MrkDAGAYJncAQKf01T01x/5vcrhLCBqTOwAAhiHcAQAwDLvlAfSI3v5gU2cvmNMXdiOboq8eBggFkzsAAIZhcgfQJ3VmGuzMXgEmSkQCJncAAAzD5A4AwBd09vMbkYjJHQAAwzC5A4CBvjx9fn7xIj4T0DcwuQMAYBgmdwDoBaZ/V5sLtkQWJncAAAzD5A4gaD31KeJImfoi/VPSfXX6D3XbvojJHQAAwzC5A+hV3TFxmT61McEiVEzuAAAYhnAHAMAw7JYHgD6OXfvmYXIHAMAwYZ3ca2trVVxcrOPHjysmJkYOh0PXXXedGhsbVVRUpJqaGqWmpmrp0qVKSEgIZ6kAAESNsIZ7bGysbrnlFmVmZqqpqUkrVqzQZZddpq1btyonJ0d5eXkqKytTWVmZ5s+fH85SAQCIGmENd5vNJpvNJkk6//zzlZGRIZfLpaqqKj344IOSpNzcXD344IOEOwB0A46v9w0Rc8y9urpa+/fv18iRI1VfX+8PfZvNpoaGhjBXBwBA9IiIT8s3NzersLBQCxYsUHx8fNCPczqdcjqdkqSCggLZ7faztrFarQHXg3Hs3JsA6EFfnDKj5fexo79vAtXfmW2jVWf+/u3JnzsS6uhqFnVF2MO9ra1NhYWFmjp1qq666ipJUlJSktxut2w2m9xutxITEwM+1uFwyOFw+G/X1taetY3dbg+4DgA9oTN/3/SFv5si5WeMhDp6oob09PSA62ENd5/Pp5KSEmVkZOiGG27wr48fP14VFRXKy8tTRUWFJkyYEMYqASB4nDo2Mpl+0Z0vC2u4f/DBB9q2bZsuvPBC3XvvvZKk733ve8rLy1NRUZHKy8tlt9u1bNmycJYJAEBUCWu4jxo1Si+99FLA+1atWtXL1QAA+hKT95xEzKflAQBA9yDcAQAwDOEOAIBhCHcAAAxDuAMAYBjCHQAAw4T9DHUAgL7F5K+gRQomdwAADEO4AwBgGMIdAADDcMwdANBjOL4eHkzuAAAYhnAHAMAwhDsAAIYh3AEAMAzhDgCAYQh3AAAMQ7gDAGAYwh0AAMMQ7gAAGIZwBwDAMIQ7AACGIdwBADAM4Q4AgGEIdwAADEO4AwBgGMIdAADDEO4AABiGcAcAwDCEOwAAhiHcAQAwDOEOAIBhCHcAAAxDuAMAYBhruAv4Krt27dKmTZvk9Xo1c+ZM5eXlhbskAAAiXsRO7l6vVxs3btTKlStVVFSkN954Q4cOHQp3WQAARLyIDfd9+/Zp8ODBSktLk9Vq1eTJk1VVVRXusgAAiHgRG+4ul0spKSn+2ykpKXK5XGGsCACA6BCxx9x9Pt9ZazExMe1uO51OOZ1OSVJBQYHS09MDPldH6+e0ZXvXHgcAQBhF7OSekpKiuro6/+26ujrZbLZ22zgcDhUUFKigoKDD51mxYkWP1dhX0MPQ0cPQ0cPQ0L/QRVMPIzbcR4wYoSNHjqi6ulptbW2qrKzU+PHjw10WAAARL2J3y8fGxmrhwoV6+OGH5fV6dc0112jYsGHhLgsAgIgXseEuSVdccYWuuOKKkJ7D4XB0UzV9Fz0MHT0MHT0MDf0LXTT1MMYX6JNrAAAgakXsMXcAANA1Eb1b/lxqa2tVXFys48ePKyYmRg6HQ9ddd50aGxtVVFSkmpoapaamaunSpUpISJAklZaWqry8XBaLRbfddpvGjRsX3h8izFpbW/XAAw+ora1NHo9HEydO1Jw5c+hhF3i9Xq1YsULJyclasWIFPeyku+66S/3795fFYlFsbKwKCgroYSedPHlSJSUlOnjwoGJiYrR48WKlp6fTwyAdPnxYRUVF/tvV1dWaM2eOcnNzo6+Hvijmcrl8H3/8sc/n8/lOnTrly8/P9x08eNC3efNmX2lpqc/n8/lKS0t9mzdv9vl8Pt/Bgwd9P/vZz3ytra2+Y8eO+e6++26fx+MJV/kRwev1+pqamnw+n893+vRp33333ef74IMP6GEXvPLKK761a9f6Hn30UZ/P56OHnbRkyRJffX19uzV62DlPPvmkz+l0+ny+M7/PjY2N9LCLPB6P74c//KGvuro6KnsY1bvlbTabMjMzJUnnn3++MjIy5HK5VFVVpdzcXElSbm6u/7S1VVVVmjx5svr166dBgwZp8ODB2rdvX9jqjwQxMTHq37+/JMnj8cjj8SgmJoYedlJdXZ127typmTNn+tfoYejoYfBOnTql999/XzNmzJAkWa1WDRgwgB520e7duzV48GClpqZGZQ+jerf8F1VXV2v//v0aOXKk6uvr/Se8sdlsamhokHTmlLZZWVn+xyQnJ3NKW53Znbx8+XIdPXpU1157rbKysuhhJz333HOaP3++mpqa/Gv0sPMefvhhSdI3vvENORwOetgJ1dXVSkxM1Pr16/Xpp58qMzNTCxYsoIdd9MYbb+jqq6+WFJ2/y0aEe3NzswoLC7VgwQLFx8d3uJ2PLwYEZLFY9Nhjj+nkyZN6/PHH9dlnn3W4LT08244dO5SUlKTMzEzt2bPnnNvTw8BWr16t5ORk1dfX65e//OVXnjaaHp7N4/Fo//79WrhwobKysrRp0yaVlZV1uD097FhbW5t27NihefPmfeV2kdzDqA/3trY2FRYWaurUqbrqqqskSUlJSXK73bLZbHK73UpMTJR09iltXS6XkpOTw1J3JBowYIDGjBmjXbt20cNO+OCDD7R9+3a99dZbam1tVVNTk9atW0cPO+nzHiQlJWnChAnat28fPeyElJQUpaSk+CfJiRMnqqysjB52wVtvvaWLL75YAwcOlBSdmRLVx9x9Pp9KSkqUkZGhG264wb8+fvx4VVRUSJIqKio0YcIE/3plZaVOnz6t6upqHTlyRCNHjgxL7ZGioaFBJ0+elHTmk/O7d+9WRkYGPeyEefPmqaSkRMXFxbrnnnt06aWXKj8/nx52QnNzs/+QRnNzs9555x1deOGF9LATBg4cqJSUFB0+fFjSmWPGQ4cOpYdd8MVd8lJ0ZkpUn8Rm7969WrVqlS688EL/FeO+973vKSsrS0VFRaqtrZXdbteyZcv8X1t4+eWX9c9//lMWi0ULFizQ5ZdfHs4fIew+/fRTFRcXy+v1yufzadKkSbr55pt14sQJetgFe/bs0SuvvKIVK1bQw044duyYHn/8cUlndi9PmTJFs2bNooeddODAAZWUlKitrU2DBg3SkiVL5PP56GEntLS0aPHixXrqqaf8h3mj8X0Y1eEOAADOFtW75QEAwNkIdwAADEO4AwBgGMIdAADDEO4AABiGcAcM8tJLL2ndunXd/rwPPvigXn/99W5/XgA9g3AHelBpaakeffTRdmv5+fkB1954443eLK3P2rNnj+68885wlwH0KMId6EGjR4/W3r175fV6JUnHjx+Xx+PRJ5980m7t6NGjGj16dKee2+PxdHu9AMwQ9eeWByLZyJEj5fF4dODAAWVmZuq9997TJZdcomPHjrVbS0tL819RasOGDdq7d68SEhJ00003yeFwSDqzy/3gwYPq16+fduzYoVtvvVU5OTkqLi7W/v37lZWV1e5iK62trSopKdGuXbvk9Xo1ZMgQLV++3H++7K/yzjvv6Nlnn5Xb7da0adPaXSDD6/WqtLRUr7/+ulpbWzVu3DgtXLjQfzavvXv36vnnn9ehQ4d0/vnna+7cuZo+fboefPBBTZ061X9Z3K1bt+r111/X6tWrJUlz5szR7bffri1btuj48eO67rrrNH36dD355JM6dOiQxo4dq/z8fFmtZ/7a2rFjh37/+9+rpqZGQ4cO1aJFi3TRRRdJku666y5de+212rZtm2pqajRu3Djddddd8nq9euSRR9TW1qZbbrlFkvTEE0/I5XLpt7/9rY4cOaK4uDhNmTJFP/jBD0L80wfCh3AHepDValVWVpbee+89ZWZm6v3339eoUaNks9narX0+tT/xxBMaNmyYnnnmGR0+fFirV69WWlqacnJyJEnbt2/X0qVLdffdd6utrU0PPfSQsrOz9fOf/1wfffSRCgoKNH78eElnzoF96tQpPf300+rXr58OHDiguLi4c9bc0NCgwsJCLV68WOPHj9err76q1157TdOmTZN0JpS3bt2qBx54QElJSXrqqae0ceNG/fjHP1Ztba0eeeQR3XHHHZo4caKampraXVjjXHbt2qWCggLV1dVp+fLl+vDDD5Wfn68LLrhA999/v/79739r+vTp+uSTT/T0009r+fLlGjFihLZt26Zf//rXWrt2rfr16ydJ+s9//qOVK1cqLi5Ov/jFL7R161Z985vf1MqVK/Xkk0+qpKTE/7qFhYW67rrrNG3aNDU3N3/llRGBaMBueaCHjR49Wu+//76kM1Pt6NGjz1obM2aMamtrtXfvXn3/+99XXFychg8frpkzZ2rbtm3+58rOztbXv/51WSwWNTQ06OOPP9bcuXPVr18/jRkzRldeeaV/29jYWDU2Nuro0aOyWCzKzMz8yksif+6tt97S0KFDNXHiRFmtVl1//fXtpv1///vfuuGGG5SWlqb+/ftr3rx5qqyslMfj0b/+9S/l5ORoypQpslqtuuCCCzR8+PCge3XTTTcpPj5ew4YN07Bhw3TZZZcpLS1N8fHxuvzyy3XgwAFJ0uuvvy6Hw6GsrCxZLBZNnz5dVqtVH330kf+5vvWtbyk5OVkJCQm68sor/Y8NxGq16ujRo2poaFD//v2VnZ0ddM1AJGJyB3rYmDFj9I9//EONjY1qaGjQkCFDlJSUpOLiYjU2Nuqzzz7TmDFj5Ha7lZCQoPPPP9//WLvdro8//th/OyUlxf//LpdLAwYMUP/+/f1rqampqq2tlSRNmzZNdXV1Wrt2rU6dOqWpU6fqu9/9rn+3dkfcbne714mJiWl32+12KzU1tV2NHo9H9fX1qqurU1paWhe6dMYX/xERFxd31u3jx49Lkmpra1VRUaFXX33Vf39bW5tcLleHz/XF+77szjvv1B/+8ActXbpUgwYN0s0339zuH0pAtCHcgR6WnZ2tU6dOyel06mtf+5okKT4+XjabTU6nU8nJyRo0aJAsFosaGxvV1NTkD/ja2toOrw9ts9l08uRJNTc3+wP+82CXzkyjs2fP1uzZs1VdXa1HH31U6enpmjFjxlfWO3DgwHa70n0+X7vbNptNNTU1/tu1tbWKjY1VUlKSUlJStG/fvoDPe95556mlpcV/+/Og7oqUlBTNmjVLs2bN6vRjP7+C5BcNGTJE99xzj7xer/73v/9pzZo12rhxY7t/OAHRhN3yQA+Li4vTiBEjtGXLFo0aNcq/PmrUKG3ZssV/vN1ut+trX/uaXnjhBbW2turTTz/VP//5T02dOjXg86ampmrEiBF66aWX1NbWpr1792rHjh3++99991199tln8nq9io+Pl9VqlcVy7l/5K664QgcPHtR///tfeTwe/f3vf28XxFdffbW2bNmi6upqNTc368UXX9SkSZMUGxurqVOnavfu3f7d9CdOnPDvDh8+fLj+97//qaWlRUePHlV5eXkXunnGzJkz9dprr+mjjz6Sz+dTc3Ozdu7c6b8m/FdJSkrSiRMndOrUKf/atm3b1NDQIIvF4j90EUyvgEjF5A70gjFjxujDDz88K9xfffXVdl+B+8lPfqINGzboRz/6kRISEjR79mxddtllHT5vfn6+iouLddtttyk7O1vTpk3TyZMnJZ2ZjDds2CCXy6X+/ftr0qRJHf5D4YsSExO1bNkybdq0SevXr9e0adP8exwk6ZprrpHb7dYDDzyg1tZWjR07VgsXLpR05h8o9913nzZv3qxnnnlG8fHxmjt3roYPH67rr79eH3/8sf9T7VOmTNHu3bs73UtJGjFihH70ox/p2Wef9X/CfdSoUUF9nTAjI0NXX3217r77bnm9Xq1Zs0a7du3S7373O7W0tCg1NVU/+clP/B8+vOWWW7Ry5cpOf1URCCeu5w4AgGHY7wQAgGEIdwAADEO4AwBgGMIdAADDEO4AABiGcAcAwDCEOwAAhiHcAQAwDOEOAIBh/h9keVMI+b8mnAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "### DEFACTORING INSPECTION\n",
    "plt.style.use('ggplot')\n",
    "ax = pd.Series([x[1] for x in wordcounts.most_common(n=3200)]).hist(\n",
    "    bins=72, \n",
    "    figsize=(8,5))\n",
    "ax.set_xlabel(\"Words  documents.\")\n",
    "ax.set_ylabel(\"Count\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot above shows a histogram of the top 3,200 words and how they are expressed across corpus. The spike on the right end of this chart shows there are nearly 60 words that appear in all 720 documents (as we can see in the text table above). As a whole, the higher bars on the left side of the chart indicate most of the words appear in a smaller number of documents. Here we use defactoring as a technique to investigate and even generate intermediate representations of the data, representations implicit in the data structures created by Underwood and Sellers, but not explicitly visualized in their MLQ narrative. For our purposes, this image is an interesting chapter in the story of the data precisely because it is in the middle of Underwood and Seller's analysis. These middle states are often glossed over in the hurried rush for analysis to generate a meaningful result. Defactoring is an effort to slow down, take a breather, and reflect upon the data-work that has happened up until this point in the code. The meandering step-by-step journey through the code sometimes reveals very interesting paths not taken, such as the commented out code block below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# vocablist = binormal_select(vocablist, positivecounts, negativecounts, totalposvols, totalnegvols, 3000)\n",
    "# Feature selection is deprecated. There are cool things\n",
    "# we could do with feature selection,\n",
    "# but they'd improve accuracy by 1% at the cost of complicating our explanatory task.\n",
    "# The tradeoff isn't worth it. Explanation is more important.\n",
    "# So we just take the most common words (by number of documents containing them)\n",
    "# in the whole corpus. Technically, I suppose, we could crossvalidate that as well,\n",
    "# but *eyeroll*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Underwood and Seller's code above does not actually perform any work as each line has been commented out, however we include it because it points towards an execution path not taken and an interesting rationale for why it was not followed. In the \"production\" code the heuristic for feature selection is to simply select the 3200 most common words by their appearance in the 720 documents. This is a simple and easy technique to implement and--more importantly--explain to a literary history and digital humanities audience. Selecting the top words is a well established practice in text analysis and it has a high degree of face validity. It is a good mechanism for removing features that have diminishing returns. However, the commented code above tells a different, and methodologically significant, story. The comment discusses an alternative technique for feature selection using binormal selection. Because this function is commented out and not used in the analysis, we have opted to not include it as part of the defactoring. Instead, we have decided to focus on the more interesting rationale about why binormal selection is not being used in the analysis as indicated in the comments:\n",
    "\n",
    "> There are cool things we could do with feature selection, but they'd improve accuracy by 1% at the cost of complicating our explanatory task. The tradeoff isn't worth it. Explanation is more important.\n",
    "\n",
    "This comment reveals much about the reasoning, the effort, and energy focused on the important, but in the humanities oft neglected, work of discussing methodology. As Underwood argued in *The literary uses of high-dimensional space* (Underwood 2015b), while there is enormous potential for the application of statistical methods in humanistic fields like literary history there is resistance to these methods because there is a resistance to methodology. Underwood has described the humanities disciplines relationship to methodology as an \"insistence on staging methodology as ethical struggle\" (Underwood 2013). In this commented code we can see the material manifestation of Underwood's methodological sentiment, in this case embodied by self-censorship in the decision to not use more statistically robust techniques for feature selection. We do not argue this choice compromises the analysis or final conclusions, rather we want to highlight the practical and material ways research methods are not a metaphysical abstraction, but rather have an tangible and observable reality. By focusing on a close reading of the code and execution environment, by *defactoring*, we illuminate methodology and its relation to the omnipresent *explanatory* task commensurate with the use of computational research methods in the humanities.\n",
    "\n",
    "In an algorithmic, data driven analysis, the selection of features is a *crucial* step because it affects the accuracy of the algorithm. In the digital humanities, feature selection is deeply embedded in the theory of the analysis and the context of the data. Claims made in and through this kind of analysis must attend to the representational configuration of the data. That is to say, we cannot take for granted how we have transformed data and what data are included or excluded from the analysis. Care, in the form of thorough documentation and thoughtful reflection, must be taken--especially at this unique moment in the development of digital humanities as we are still learning how algorithmic, data-driven techniques can be leveraged to better understand our objects of study."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Filtering Authors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "donttrainon = list()\n",
    "\n",
    "# Here we create a list of volumed IDs not to be used for training.\n",
    "# For instance, we have supplemented the dataset with volumes that\n",
    "# are in the Norton but that did not actually occur in random\n",
    "# sampling. We want to make predictions for these, but never use\n",
    "# them for training.\n",
    "\n",
    "for idx1, anid in enumerate(orderedIDs):\n",
    "    reviewedstatus = metadict[anid]['reviewed']\n",
    "    date = infer_date(metadict[anid], datetype)\n",
    "    if reviewedstatus == 'addedbecausecanon':\n",
    "        donttrainon.append(idx1)\n",
    "    elif date < pastthreshold or date > futurethreshold:\n",
    "        donttrainon.append(idx1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the comments describe, this block of code creates a list of volume IDs not to be used in the training. What that means *in code* is that any volume with the metadata label `addedbecauseofcanon` or with a date outside of the thresholds defined by `pastthreshold` and `futurethreshold` will be ignored. If we inspect the `donttrainon` variable we can see how many volumes satisfy these criteria."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The variable donttrainon contains 0 volume IDs\n"
     ]
    }
   ],
   "source": [
    "### DEFACTORING INSPECTION\n",
    "print(\"The variable donttrainon contains {} volume IDs\".format(\n",
    "    len(donttrainon)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It would appear there are no volumes to be filtered out by these criteria."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "authormatches = [list(donttrainon) for x in range(len(orderedIDs))]\n",
    "# For every index in authormatches, identify a set of indexes that have\n",
    "# the same author. Obvs, there will always be at least one.\n",
    "\n",
    "# Since we are going to use these indexes to exclude rows, we also add\n",
    "# all the ids in donttrainon to every volume\n",
    "\n",
    "for idx1, anid in enumerate(orderedIDs):\n",
    "    thisauthor = metadict[anid]['author']\n",
    "    for idx2, anotherid in enumerate(orderedIDs):\n",
    "        otherauthor = metadict[anotherid]['author']\n",
    "        if thisauthor == otherauthor and not idx2 in authormatches[idx1]:\n",
    "            authormatches[idx1].append(idx2)\n",
    "\n",
    "for alist in authormatches:\n",
    "    alist.sort(reverse = True)\n",
    "\n",
    "# I am reversing the order of indexes so that I can delete them from\n",
    "# back to front, without changing indexes yet to be deleted.\n",
    "# This will become important in the modelingprocess module."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this block of code Underwood and Sellers group the volumes by the same author. The list `authormatches` is a list of lists for each volume. Each sub-list contains the identifiers for all the volumes by the same author. Essentially this data structure represents the potential relations of each volume to other volumes, with that relation being \"other volumes by the same author.\" This raises the question, how many volumes share the same author in this corpus."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1    579\n",
       "2     82\n",
       "3     30\n",
       "5     15\n",
       "4      8\n",
       "6      6\n",
       "dtype: int64"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "### DEFACTORING INSPECTION\n",
    "# Tabular view of shared authorship\n",
    "pd.Series([len(x) for x in authormatches]).value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAFgCAYAAABHS1h8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAm6klEQVR4nO3de3CV9Z3H8c9JDiEQJLcToCFCJpA0REaRmxiuhjB2xEXXqghDncCWqiDQWl3RtWSpMkJragQTYhFDoR0X2y0WqN12I8q1bgVkuYNRruUSc7+YkJzL/sGY2TQk/M7lScLh/Zph5pznnPM8n3z/0M88V5vH4/EIAAAgSIR0dgAAAIBAotwAAICgQrkBAABBhXIDAACCCuUGAAAEFcoNAAAIKvbODtCeCxcudHaEoOJwOFRaWtrZMYIOcw08ZmoN5hp4zNQa7c01Pj7+ur9nzw0AAAgqlBsAABBUKDcAACCoUG4AAEBQodwAAICgQrkBAABBhXIDAACCCuUGAAAEFcoNAAAIKpQbAAAQVCg3AAAgqHTpZ0u55k5r8T50zeZOSgIAAG4U7LkBAABBhXIDAACCCuUGAAAEFcoNAAAIKpQbAAAQVCg3AAAgqFBuAABAUKHcAACAoEK5AQAAQaXD7lBcV1engoICnTt3TjabTU899ZRSUlI6avMAAOAm0WHlprCwUMOGDdOPf/xjOZ1OXblypaM2DQAAbiIdcljq66+/1rFjx5SRkSFJstvtioiI6IhNAwCAm0yH7LkpKSlR7969lZ+frzNnzigpKUlZWVkKDw9v8b2ioiIVFRVJkpYvX95qPQ6HoyPiBi273c4MLcBcA4+ZWoO5Bh4ztYa/c+2QcuNyuXTq1CnNmTNHycnJKiws1Pvvv6/HHnusxfcyMzOVmZnZ5npKS0utjhrUHA4HM7QAcw08ZmoN5hp4zNQa7c01Pj7+ur/vkMNSsbGxio2NVXJysiRpzJgxOnXqVEdsGgAA3GQ6pNxERUUpNjZWFy5ckCQdOnRICQkJHbFpAABwk+mwq6XmzJmjlStXyul0qk+fPpo3b15HbRoAANxEOqzcJCYmXvMkYQAAgEDiDsUAACCoUG4AAEBQodwAAICgQrkBAABBhXIDAACCCuUGAAAElQ67FNwXoWs2d3YEAABwg2HPDQAACCrXLTdut1vbtm1TU1NTR+QBAADwy3XLTUhIiNavX69u3bp1RB4AAAC/GB2WGjFihPbu3Wt1FgAAAL8ZnVDc1NSkX/ziF0pJSVFsbKxsNlvzZ08//bRl4QAAALxlVG5uvfVW3XrrrVZnAQAA8JtRuXnkkUeszgEAABAQxve5OXz4sHbs2KGKigpFR0drwoQJGjp0qJXZAAAAvGZ0QvGHH36o3NxcRUVFafTo0YqOjtYbb7yhoqIiq/MBAAB4xWjPzebNm/XSSy8pMTGxeVl6erpycnKUmZlpVTYAAACvGe25qampUUJCQotl8fHxqq2ttSQUAACAr4zKTWpqqtavX68rV65IkhoaGrRhwwalpKRYGg4AAMBbRoel5s6dq9zcXGVlZalXr16qra1VSkqKFi1aZHU+AAAArxiVm+joaC1dulRlZWXNV0vFxsZanQ0AAMBrXj0V3G6365ZbbpHT6dTly5d1+fJlq3IBAAD4xGjPzYEDB7R69WpVVla2+mzjxo2BzgQAAOAzo3Kzdu1affe739WkSZMUFhZmdSYAAACfGZWb2tpaTZkypcUDMwEAALoio3NuMjIy9NFHH1mdBQAAwG9t7rlZsmRJ854aj8ejDz74QH/4wx8UFRXV4ntLly61NCAAAIA32iw3GRkZ7b4HAADoitosN5MmTWp+/fnnnys5ObnVd4qLiy0JBQAA4Cujc25eeeWVay5ftmxZQMMAAAD4q92rpdxut6Sr59x88+8bly9fVmhoqLXpAAAAvNRuuZkxY0bz68cee6zFZyEhIfrnf/5na1IBAAD4qN1y8+abb8rj8ejf//3fW1wVZbPZ1Lt3b8tv6OeaO63VstA1my3dJgAAuLG1W27i4uIkSfn5+R0SBgAAwF9Gdyh+88032/zs6aefDlgYAAAAfxmVm759+7Z4X1lZqU8++UTjx4+3JBQAAICvjMrNI4880mpZRkaGfvvb3wY8EAAAgD+M7nNzLYmJiTp27FggswAAAPjNaM/N4cOHW7y/cuWKdu/erYSEBEtCAQAA+Mqo3KxevbrF+/DwcA0cOFCLFi2yJBQAAICvjMpNXl6e1TkAAAACwutzbjwej9xud/M/b7jdbv3rv/6rli9f7u1mAQAAjBjtuSkvL9fatWt17Ngx1dXVtfhs48aNxhv74IMP1L9/f9XX13uXEgAAwJDRnptf/vKXstvtWrJkicLDw7VixQqNHDlSc+fONd5QWVmZ9u/fr8mTJ/scFgAA4HqMys3Jkyf11FNPKTExUTabTYmJiXrqqae0detW4w2tW7dOs2bNks1m8zksAADA9RgdlgoJCVFoaKgkKSIiQtXV1erRo4fKy8uNNrJv3z5FRkYqKSlJR44cafN7RUVFKioqkqQ2z8txOBxG20Rrdrud+VmAuQYeM7UGcw08ZmoNf+dqVG4GDx6szz77TKNHj9Ydd9yh119/XWFhYRo0aJDRRk6cOKG9e/fqs88+U2Njo+rr67Vy5UotXLiwxfcyMzOVmZnZ7rpKS0uNtonWHA4H87MAcw08ZmoN5hp4zNQa7c01Pj7+ur83KjcLFiyQx+ORJGVlZWnLli2qr6/X1KlTjULOnDlTM2fOlCQdOXJEW7ZsaVVsAAAAAsGo3ERERDS/DgsL03e/+13LAgEAAPjDqNwE0m233abbbrutozcLAABuEj4/OBMAAKArotwAAICg4lW5cbvdqqiosCoLAACA34zOuamrq9Pbb7+tTz75RHa7XRs2bNDevXtVXFysxx57zOqMAAAAxoz23KxZs0Y9e/ZUfn6+7ParfSglJUV79uyxNBwAAIC3jPbcHDp0SG+99VZzsZGk3r17q6qqyrJgAAAAvjDac9OzZ0/V1NS0WFZaWqro6GhLQgEAAPjKaM/N5MmTlZOTo8cee0wej0cnT57Uu+++qylTplgaLnTNZkvXDwAAgo9RuXnggQfUrVs3rV27Vi6XS6tXr1ZmZqbuu+8+q/MBAAB4xajc2Gw2TZ061fhZUgAAAJ3F+PELJSUlOnv2rBoaGlosHzduXMBDAQAA+Mqo3GzatEm/+93vdOuttyosLKx5uc1mo9wAAIAuxajcbN26VStWrFBCQoLVeQAAAPxidCl4r169FBcXZ3UWAAAAvxntucnKytJbb72lqVOnKjIyssVnDofDkmAAAAC+MCo3TqdTBw8e1O7du1t9tnHjxoCHAgAA8JVRuXn77bc1Y8YMjR07tsUJxQAAAF2NUblxu9265557FBJidIoOAABApzFqK//0T/+k999/Xx6Px+o8AAAAfjHac/OnP/1JlZWV2rRpk3r16tXis9WrV1sSDAAAwBdG5WbBggVW5wAAAAgIo3KTlpZmdQ4AAICAMH621OnTp3Xs2DHV1NS0OPdm+vTplgQDAADwhVG5KSoq0q9+9SvdfvvtOnDggIYNG6aDBw9q5MiRVucDAADwitHVUn/4wx/04osv6rnnnlNYWJiee+45PfPMMwoNDbU6HwAAgFeMyk11dbWGDBki6eqTwN1ut+68807t27fP0nAAAADeMjosFRMTo5KSEvXp00ff+ta3tHfvXt1yyy2y241P2QEAAOgQRu3kgQce0N///nf16dNHDz/8sH7xi1/I6XRq9uzZVucDAADwilG5mTRpUvPrO++8U4WFhXI6nQoPD7cqFwAAgE+Mz7lpaGiQdPU5Uzt37tTf/vY3ud1uS8MBAAB4y6jcLF++XBcvXpQkvfvuu9qyZYu2bNmi9evXWxoOAADAW0bl5uLFi0pMTJQk7dy5Uy+++KKys7O1Z88eK7MBAAB4zeicm5CQEDmdTl28eFE9e/aUw+GQ2+1uPlQFAADQVRiVm2HDhun1119XTU2N0tPTJUnnz59XTEyMpeFcc6cZfzd0zWYLkwAAgBuFUbl58skntX37doWGhmrixImSpJqaGj3yyCOWhgMAAPCWUbnp1q2bMjMzWyy77bbbLAkEAADgD6MTigEAAG4UlBsAABBUKDcAACCoUG4AAEBQMTqhuKSkRO+++65Onz7d6t42q1evtiQYAACAL4zKzRtvvKG+ffvq8ccfV/fu3a3OBAAA4DOjcnP+/Hm9/PLLCgnx/SjW/PnzFR4erpCQEIWGhmr58uU+rwsAAKAtRuVmyJAhOn36tJKSkvzaWHZ2tnr37u3XOgAAANrTZrnZuHFj8+u4uDgtW7ZMo0ePVlRUVIvvTZ8+3bJwAAAA3mqz3JSVlbV4P2LECLlcrlbLvbFs2TJJ0pQpU1rd8ViSioqKVFRUJEleH7ZyOBw+57pZ2O125mQB5hp4zNQazDXwmKk1/J2rzePxeAKYp03l5eWKiYlRVVWVXnnlFc2ePVtpaWnt/ubc1JHG6+fBmdfncDhUWlra2TGCDnMNPGZqDeYaeMzUGu3NNT4+/rq/NzpDePbs2ddc/v3vf9/k55LU/ATxyMhIjRo1SsXFxca/BQAAMGVUblwuV6tlTqdTbrfbaCMNDQ2qr69vfn3w4EENGDDAi5gAAABm2r1aasmSJbLZbGpqalJ2dnaLz8rKypSSkmK0kaqqKr322muSrhalcePGadiwYb4lBgAAaEe75SYjI0OSVFxcrHvuuad5uc1mU2RkpIYOHWq0kb59++rnP/+5HzEBAADMtFtuJk2aJElKTk5W//79OyIPAACAX9osNzt27NCECRMkSSdOnNCJEyeu+b1v9u4AAAB0BW2Wm927dzeXm507d7a5AsoNAADoStosNy+88ELz6388mRgAAKCrMroU/IMPPtCZM2eszgIAAOA3owdnfvHFF9q6davq6+s1ZMgQpaWlKS0tTYmJiX49KRwAACDQjMrNggULJEklJSU6evSojh49qt/97neSpHXr1lkWDgAAwFtG5UaSLly4oKNHj+rIkSM6ceKEvvWtb1332VD+4nlRAADAW0blZu7cuQoPD9eYMWM0ceJE/eAHP1CPHj2szgYAAOA1o3IzYsQIHT9+XJ9++qnq6upUW1urIUOGKDY21up8AAAAXjEqN08++aQkqbKyUseOHdPRo0f19ttv65ZbbtGqVassDQgAAOAN43NuTp061XzOzbFjx9S9e3cNHjzYymwAAABeMyo3s2fPVs+ePTVkyBCNHDlSjz/+uPr162d1NgAAAK8ZlZsVK1aoT58+VmcBAADwm9Ed+Cg2AADgRsHthQEAQFCh3AAAgKBCuQEAAEHF6ITiw4cPq0+fPurTp48qKir0m9/8RiEhIZo5c6aioqIsjggAAGDOaM/N2rVrm5/+vX79erlcLtlsNr311luWhgMAAPCW0Z6b8vJyORwOuVwu/e///q/y8/Nlt9v1xBNPWJ0PAADAK0blpkePHqqsrNS5c+eUkJCg8PBwOZ1OOZ1Oq/MBAAB4xajcfOc739ELL7wgp9OprKwsSdLx48fVv39/K7MBAAB4zajcPPjggxo9erRCQkKaH7sQExPT/EBNAACArsL4UvBvrpTas2ePpKvlhjsXAwCArsZoz83Zs2e1YsUKdevWTWVlZUpPT9fRo0e1fft2/ehHP7I6IwAAgDGjPTdr1qzR9OnTlZubK7v9ah9KS0vT8ePHLQ0HAADgLaNyc/78eY0fP77FsvDwcDU2NloSCgAAwFdG5SYuLk5ffvlli2XFxcXNJxcDAAB0FUbn3EyfPl3Lly/XlClT5HQ6tWnTJv33f/83N/EDAABdjtGemxEjRuiFF15QdXW10tLS9NVXX+nZZ5/VHXfcYXU+AAAArxjtuZGkpKQkJSUlWZkFAADAb0blxuVyaffu3Tp16pQaGhpafGbloSnX3GmWrVuSQtdstnT9AACg4xmVm1WrVuns2bMaNmyYIiMjrc4EAADgM6Nyc+DAAa1evVo9evSwOg8AAIBfjE4oTkhIUG1trdVZAAAA/Ga052bBggUqKCjQHXfc0eqw1MSJEy0JBgAA4AujcvPxxx/r+PHjqqurU1hYWPNym81GuQEAAF2KUbn54IMPtGLFCiUkJFidBwAAwC9G59xERUXJ4XBYnQUAAMBvRntupk6dqpUrV+rBBx9sdc5N3759LQkGAADgC6Nys3btWknSvn37Wn22cePGwCYCAADwg1G58bfANDY2Kjs7W06nUy6XS2PGjNGjjz7q1zoBAACuxfjZUv7o1q2bsrOzFR4eLqfTqSVLlmjYsGFKSUnpiM0DAICbiFG5WbJkiWw22zU/W7p06XV/b7PZFB4eLunqc6pcLleb6wMAAPCHUbnJyMho8b6yslIfffSRxo8fb7wht9ut559/XpcuXdK9996r5OTkVt8pKipSUVGRJGn58uXG6/bVzXYFmN1uv+n+5o7AXAOPmVqDuQYeM7WGv3O1eTwejy8/vHTpkvLz8/XTn/7Uq9/V1dXptdde0+zZszVgwIB2v3tu6khfohm72Z4K7nA4VFpa2tkxgg5zDTxmag3mGnjM1BrtzTU+Pv66vze6z821xMTE6MyZM17/LiIiQmlpaTpw4ICvmwYAAGiT0WGpbdu2tXjf2Nio//mf/zE+Ibi6ulqhoaGKiIhQY2OjDh06pAceeMD7tAAAANdhVG527tzZ4n337t317W9/W1OnTjXaSEVFhfLy8uR2u+XxeHT33XdrxIgR3qcFAAC4DqNyk52d7ddGBg4cqJ/97Gd+rQMAAMBEm+Xm8uXLRivg8QsAAKArabPcLFy40GgFPH4BAAB0JW2WG0oLAAC4EXn1+IXS0lKVl5crJiaGmxYBAIAuyajcVFRUKDc3VydPntQtt9yimpoapaSkaNGiRYqJibE6IwAAgDGjm/itWbNGAwcOVGFhoX75y1+qsLBQiYmJWrNmjdX5AAAAvGJUbk6cOKHHH3+8+eGX4eHhmjVrlk6ePGlpOAAAAG8ZHZaKiIjQ+fPnlZiY2LzswoUL6tmzp1W5JN18z34CAAD+Myo306ZN08svv6yMjAzFxcXpq6++0scff6zp06dbnQ8AAMArRuUmMzNT/fr1065du3T27FlFR0dr0aJFGjp0qNX5AAAAvGJUbqqrqzV06FDKDAAA6PKMys28efN02223ady4cRo9erS6d+9udS4AAACfGF0tlZ+fr+HDh+svf/mL5s6dq9zcXO3du1cul8vqfAAAAF4x2nPTu3dv3Xvvvbr33ntVWlqqXbt26T/+4z+0evVqrV271uqMAAAAxoz23Px/lZWVqqysVE1NjSIiIqzIBAAA4DOjPTfnz5/Xrl27tHv3bjU2Nuruu+/Wc889p8GDB1udDwAAwCtG5eYnP/mJ7rrrLv3gBz/Q0KFDZbPZrM4FAADgE6Nys2bNGtntXj1AHAAAoFMYnXNDsQEAADcKr08oBgAA6MooNwAAIKi0WW7+7d/+rfn1b3/72w4JAwAA4K82y82FCxfU2NgoSdq6dWuHBQIAAPBHm2cKjxo1SosWLVKfPn3U2Nio7Ozsa35v6dKlloUDAADwVpvlZt68eTp+/LhKSkpUXFyse+65pyNzAQAA+KTda7xTU1OVmpoqp9OpSZMmdVAkAAAA3xndwCYjI0OHDx/Wjh07VFFRoejoaE2YMEFDhw61Oh8AAIBXjC4F//DDD5Wbm6uoqCiNHj1a0dHReuONN1RUVGR1PgAAAK8Y7bnZvHmzXnrpJSUmJjYvS09PV05OjjIzM63KBgAA4DWjPTc1NTVKSEhosSw+Pl61tbWWhAIAAPCVUblJTU3V+vXrdeXKFUlSQ0ODNmzYoJSUFEvDAQAAeMvosNTcuXOVm5urrKws9erVS7W1tUpJSdGiRYuszgcAAOAVo3ITHR2tpUuXqqysrPlqqdjYWKuzyTV3muXbuNGErtnc2REAAOjSjMrNN2JjYzuk1AAAAPiKp4IDAICgQrkBAABB5brlxu126/Dhw3I6nR2RBwAAwC/XLTchISH62c9+Jrvdq9NzAAAAOoXRYakhQ4bo5MmTVmcBAADwm9HumLi4OL366qsaOXKkYmNjZbPZmj+bPn26ZeEAAAC8ZVRuGhsbNWrUKElSeXm5pYEAAAD8YVRu5s2b59dGSktLlZeXp8rKStlsNmVmZuq+++7za50AAADXYnyW8Pnz5/XJJ5+oqqpK//Iv/6ILFy6oqalJAwcOvO5vQ0ND9b3vfU9JSUmqr6/X4sWLdfvtt7d6GCcAAIC/jE4o/utf/6rs7GyVl5drx44dkqT6+nqtX7/eaCPR0dFKSkqSJPXo0UP9+/fn8BYAALCE0Z6b9957Tz/5yU+UmJiov/71r5KkgQMH6vTp015vsKSkRKdOndLgwYO9/i0AAMD1GJWbqqqqVoefbDZbi6umTDQ0NCgnJ0dZWVnq2bNnq8+LiopUVFQkSVq+fLlX675ZOBwOn39rt9v9+j2ujbkGHjO1BnMNPGZqDX/nalRukpKStGPHDk2cOLF52e7du73a++J0OpWTk6Px48frrrvuuuZ3MjMzlZmZabzOm1FpaanPv3U4HH79HtfGXAOPmVqDuQYeM7VGe3ONj4+/7u+Nys3s2bP1yiuvaNu2bbpy5YqWLVumCxcu6KWXXjIK6fF4VFBQoP79++v+++83+g0AAIAvjMpN//79lZubq3379mnEiBGKjY3ViBEjFB4ebrSREydOaMeOHRowYICee+45SdKMGTM0fPhw35MDAABcg/Gl4N27d1dqaqrKy8sVExNjXGwkKTU1Ve+9955PAQEAALxhVG5KS0u1cuVKff7554qIiFBdXZ0GDx6shQsXKi4uzuqMAAAAxozuc5OXl6ekpCQVFhbq7bffVmFhoQYNGqS8vDyr8wEAAHjFqNx8+eWXmjVrVvOhqPDwcM2aNUtffvmlpeEAAAC8ZVRukpOTVVxc3GLZF198oZSUFEtCAQAA+KrNc242btzY/Lpv37569dVXNXz4cMXGxqqsrEyfffaZxo0b1yEhAQAATLVZbsrKylq8/+bGe9XV1erWrZtGjx6txsZGa9MBAAB4qc1yM2/evI7MAQAAEBDG97m5cuWKLl26pIaGhhbLv/3tbwc81DdC12y2bN0AACA4GZWb7du365133pHdbldYWFiLz1avXm1JMAAAAF8YlZtf//rX+vGPf6zbb7/d6jwAAAB+MboU3G63Ky0tzeosAAAAfjMqN9OnT9f69etVXV1tdR4AAAC/GB2Wio+P13vvvac///nPrT77//fDAQAA6GxG5WbVqlWaMGGC0tPTW51QDAAA0JUYlZva2lpNnz5dNpvN6jwAAAB+MTrnZtKkSdqxY4fVWQAAAPxmtOemuLhY//Vf/6Xf//73ioqKavHZ0qVLrcgFAADgE6NyM3nyZE2ePNnqLAAAAH4zKjeTJk2yOAYAAEBgGJWbbdu2tflZRkZGwMIAAAD4y6jc7Ny5s8X7yspKXbp0SampqZQbAADQpRiVm+zs7FbLtm3bpr///e8BDwQAAOAPo0vBr2XSpEntHq4CAADoDEZ7btxud4v3jY2N2rFjhyIiIiwJBQAA4CujcjNjxoxWy2JiYvTEE08EPBAAAIA/jMrNm2++2eJ99+7d1bt3b0sCAQAA+MOo3MTFxVmdAwAAICDaLTfXe7SCzWbTkiVLAhoIAADAH+2Wm/Hjx19zeXl5uf70pz/pypUrloQCAADwVbvl5h9v0FdTU6NNmzbpww8/VHp6uh5++GFLwwEAAHjL6Jybr7/+Wps3b9af//xnDR8+XCtWrFC/fv2szgYAAOC1dstNY2Oj/vjHP2rr1q1KS0vTT3/6U916660dlQ0AAMBr7Zab+fPny+12a9q0aRo0aJCqqqpUVVXV4jtDhw61NCAAAIA32i03YWFhkqS//OUv1/zcZrO1ugdOILnmTrNs3Tejy50dIEgx18BjptZgroHHTFsKXbO5syNIuk65ycvL66gcAAAAAeHzgzMBAAC6IsoNAAAIKpQbAAAQVCg3AAAgqFBuAABAUKHcAACAoEK5AQAAQcXo2VL+ys/P1/79+xUZGamcnJyO2CQAALhJdciem0mTJunFF1/siE0BAICbXIeUm7S0NPXq1asjNgUAAG5yHXJYylRRUZGKiookScuXL+/kNAAAwBsOhyMg67Hb7X6tq0uVm8zMTGVmZnZ2DAAA4IPS0tKArMfhcLS5rvj4+Ov+nqulAABAUKHcAACAoNIhh6Vyc3N19OhR1dTU6Mknn9Sjjz6qjIyMjtg0AAC4yXRIufnhD3/YEZsBAADgsBQAAAgulBsAABBUKDcAACCoUG4AAEBQodwAAICgQrkBAABBhXIDAACCSpd6ttQ/Cl2zubMjBJX2ntUB3zHXwGOm1mCugcdMuyb23AAAgKBCuQEAAEGFcgMAAIIK5QYAAAQVyg0AAAgqlBsAABBUKDcAACCoUG4AAEBQodwAAICgQrkBAABBhXIDAACCCuUGAAAEFZvH4/F0dggAAIBA6bJ7bhYvXtzZEYIOM7UGcw08ZmoN5hp4zNQa/s61y5YbAAAAX1BuAABAUOmy5SYzM7OzIwQdZmoN5hp4zNQazDXwmKk1/J0rJxQDAICg0mX33AAAAPiCcgMAAIKKvbMD/KMDBw6osLBQbrdbkydP1oMPPtjZkW4Y+fn52r9/vyIjI5WTkyNJqq2t1euvv66vvvpKcXFx+tGPfqRevXpJkjZt2qRt27YpJCREs2fP1rBhwzoxfddUWlqqvLw8VVZWymazKTMzU/fddx9z9UNjY6Oys7PldDrlcrk0ZswYPfroo8w0QNxutxYvXqyYmBgtXryYuQbA/PnzFR4erpCQEIWGhmr58uXM1U91dXUqKCjQuXPnZLPZ9NRTTyk+Pj5wM/V0IS6Xy/P00097Ll265GlqavI8++yznnPnznV2rBvGkSNHPF988YXnmWeeaV62YcMGz6ZNmzwej8ezadMmz4YNGzwej8dz7tw5z7PPPutpbGz0XL582fP00097XC5XZ8Tu0srLyz1ffPGFx+PxeL7++mvPwoULPefOnWOufnC73Z76+nqPx+PxNDU1eV544QXPiRMnmGmAbNmyxZObm+t59dVXPR4P/w0IhHnz5nmqqqpaLGOu/lm1apWnqKjI4/Fc/e9AbW1tQGfapQ5LFRcXq1+/furbt6/sdrvS09P16aefdnasG0ZaWlpzy/3Gp59+qokTJ0qSJk6c2DzPTz/9VOnp6erWrZv69Omjfv36qbi4uMMzd3XR0dFKSkqSJPXo0UP9+/dXeXk5c/WDzWZTeHi4JMnlcsnlcslmszHTACgrK9P+/fs1efLk5mXM1RrM1Xdff/21jh07poyMDEmS3W5XREREQGfapQ5LlZeXKzY2tvl9bGysPv/8805MdOOrqqpSdHS0pKv/o66urpZ0ddbJycnN34uJiVF5eXmnZLxRlJSU6NSpUxo8eDBz9ZPb7dbzzz+vS5cu6d5771VycjIzDYB169Zp1qxZqq+vb17GXANj2bJlkqQpU6YoMzOTufqhpKREvXv3Vn5+vs6cOaOkpCRlZWUFdKZdqtx4rnFVus1m64Qkwe9as0bbGhoalJOTo6ysLPXs2bPN7zFXMyEhIfr5z3+uuro6vfbaazp79myb32WmZvbt26fIyEglJSXpyJEj1/0+czX38ssvKyYmRlVVVXrllVcUHx/f5neZ6/W5XC6dOnVKc+bMUXJysgoLC/X++++3+X1fZtqlyk1sbKzKysqa35eVlTW3OPgmMjJSFRUVio6OVkVFhXr37i2p9azLy8sVExPTWTG7NKfTqZycHI0fP1533XWXJOYaKBEREUpLS9OBAweYqZ9OnDihvXv36rPPPlNjY6Pq6+u1cuVK5hoA38wlMjJSo0aNUnFxMXP1Q2xsrGJjY5v3xowZM0bvv/9+QGfapc65GTRokC5evKiSkhI5nU7t2bNHI0eO7OxYN7SRI0dq+/btkqTt27dr1KhRzcv37NmjpqYmlZSU6OLFixo8eHBnRu2SPB6PCgoK1L9/f91///3Ny5mr76qrq1VXVyfp6pVThw4dUv/+/Zmpn2bOnKmCggLl5eXphz/8oYYOHaqFCxcyVz81NDQ0H+ZraGjQwYMHNWDAAObqh6ioKMXGxurChQuSpEOHDikhISGgM+1ydyjev3+/fvWrX8ntduuee+7RQw891NmRbhi5ubk6evSoampqFBkZqUcffVSjRo3S66+/rtLSUjkcDj3zzDPNJx3//ve/10cffaSQkBBlZWXpzjvv7OS/oOs5fvy4lixZogEDBjQfIp0xY4aSk5OZq4/OnDmjvLw8ud1ueTwe3X333Xr44YdVU1PDTAPkyJEj2rJlixYvXsxc/XT58mW99tprkq4eThk3bpweeugh5uqn06dPq6CgQE6nU3369NG8efPk8XgCNtMuV24AAAD80aUOSwEAAPiLcgMAAIIK5QYAAAQVyg0AAAgqlBsAABBUKDcAACCoUG4ABJ358+fr4MGDnR0DQCeh3AAAgKDCTfwAWK60tFTr1q3TsWPH5PF4NHbsWM2ePVubNm3Shx9+qMbGRg0bNkxz5sxRz549deTIEa1atUoFBQXN65g/f76eeOIJ3X777Xrvvfd0/vx5hYWF6W9/+5scDofmz5+vQYMGadWqVdq1a5fsdrtCQkL08MMP64EHHujEvx5AR2PPDQBLud1urVixQg6HQ3l5eSooKNDYsWP18ccf6+OPP1Z2drbefPNNNTQ0aO3atcbr3bdvn9LT07Vu3TqNHDlS77zzjiRpwYIFcjgcev7557VhwwaKDXATotwAsFRxcbHKy8v1ve99T+Hh4QoLC1Nqaqp27dql+++/X3379lV4eLhmzpypPXv2yOVyGa03NTVVw4cPV0hIiCZMmKDTp09b+4cAuGFQbgBYqrS0VHFxcQoNDW2xvKKiQnFxcc3vHQ6HXC6XqqqqjNYbGRnZ/DosLExNTU3GxQhAcKPcALCUw+FQaWlpq+IRHR2tr776qvl9aWmpQkNDFRkZqe7du+vKlSvNn7ndblVXV3dYZgA3NsoNAEsNHjxY0dHR+s1vfqOGhgY1Njbq+PHjGjt2rP74xz+qpKREDQ0Nevfdd3X33XcrNDRU8fHxampq0v79++V0OvWf//mfampqMt5mVFSUSkpKLPyrAHRl9s4OACC4hYSE6Pnnn9c777yjefPmyWazaezYscrKylJFRYWys7PV2NioO+64Q3PmzJEk9ezZU9///vdVUFAgt9utadOmKTY21nibDz74oN555x39+te/1kMPPaRp06ZZ9ecB6IK4FBwAAAQVDksBAICgQrkBAABBhXIDAACCCuUGAAAEFcoNAAAIKpQbAAAQVCg3AAAgqFBuAABAUPk/+MRPWQP02ocAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "### DEFACTORING INSPECTION\n",
    "# Barchart of shared authorship\n",
    "ax = (pd.Series([len(x) for x in authormatches])\n",
    "      .value_counts()\n",
    "      .plot(kind=\"barh\",\n",
    "           figsize=(8,5)))\n",
    "ax.set_xlabel(\"count\")\n",
    "ax.set_ylabel(\"Number of volumes with same author\")\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This histogram tells us a majority of volumes are written by unique authors but that there are some authors who have written up to six volumes in the corpus. Note, we are generating this graph by counting the length of the list containing the volume IDs of other volumes by the same author. This means volumes written by the same author are counted twice. This is not an issue for the purposes of our inspection, just that the sum total number of volumes represented by this histogram is greater than 720."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Filtering Words"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "### DEFACTORING DEFINITIONS\n",
    "\n",
    "usedate = False\n",
    "# Leave this flag false unless you plan major\n",
    "# surgery to reactivate the currently-deprecated\n",
    "# option to use \"date\" as a predictive feature.\n",
    "\n",
    "def get_features(wordcounts, wordlist):\n",
    "    numwords = len(wordlist)\n",
    "    wordvec = np.zeros(numwords)\n",
    "    for idx, word in enumerate(wordlist):\n",
    "        if word in wordcounts:\n",
    "            wordvec[idx] = wordcounts[word]\n",
    "\n",
    "    return wordvec\n",
    "\n",
    "# In an earlier version of this script, we sometimes used\n",
    "# \"publication date\" as a feature, to see what would happen.\n",
    "# In the current version, we don't. Some of the functions\n",
    "# and features remain, but they are deprecated. E.g.:\n",
    "\n",
    "def get_features_with_date(wordcounts, wordlist, date, totalcount):\n",
    "    numwords = len(wordlist)\n",
    "    wordvec = np.zeros(numwords + 1)\n",
    "    for idx, word in enumerate(wordlist):\n",
    "        if word in wordcounts:\n",
    "            wordvec[idx] = wordcounts[word]\n",
    "\n",
    "    wordvec = wordvec / (totalcount + 0.0001)\n",
    "    wordvec[numwords] = date\n",
    "    return wordvec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This code cell defines two functions to be used below when opening and parsing the raw data files (in the `poems/` directory). The function `get_features()` simply takes the word counts from the parsed volume and filters out any words that are not part of `wordlist`, which contains the list of word features that had been selected for this analysis. We have also included a second function, `get_features_with_date()`, even though it is not executed. This residual code points to yet another path not taken, one that uses the volume's publication date as a feature. As Underwood and Seller's comment indicates, this was an experiment from an \"earlier version of this script...to see what would happen.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "volsizes = dict()\n",
    "voldata = list()\n",
    "classvector = list()\n",
    "\n",
    "for volid, volpath in volspresent:\n",
    "\n",
    "    with open(volpath, encoding = 'utf-8') as f:\n",
    "        voldict = dict()\n",
    "        totalcount = 0\n",
    "        for line in f:\n",
    "            fields = line.strip().split('\\t')\n",
    "            if len(fields) > 2 or len(fields) < 2:\n",
    "                continue\n",
    "\n",
    "            word = fields[0]\n",
    "            count = int(fields[1])\n",
    "            voldict[word] = count\n",
    "            totalcount += count\n",
    "\n",
    "    date = infer_date(metadict[volid], datetype)\n",
    "    date = date - 1700\n",
    "    if date < 0:\n",
    "        date = 0\n",
    "\n",
    "    if usedate:\n",
    "        features = get_features_with_date(voldict, \n",
    "                                          vocablist, \n",
    "                                          date, \n",
    "                                          totalcount)\n",
    "        voldata.append(features)\n",
    "    else:\n",
    "        features = get_features(voldict, vocablist)\n",
    "        voldata.append(features / (totalcount + 0.001))\n",
    "\n",
    "\n",
    "    volsizes[volid] = totalcount\n",
    "    classflag = classdictionary[volid]\n",
    "    classvector.append(classflag)\n",
    "    \n",
    "data = pd.DataFrame(voldata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is an important code block because we are now pulling the raw data files from the poems/ directory into memory, filtering out the unselected word features, and putting the data into a vectorized data structure. The code loops over the `volspresent` variable and parses each individual volume into the `voldict` dictionary. At this stage the code is reading in all the words of a volume including their frequencies, and it is tabulating the total number of words in that volume. Once all of the data for the volume has been read into memory, the code calls the `get_features` function that throws out the words not part of the selected word features stored in the `vocablist` variable. This is where the top 3200 words are foregrounded for each volume and the remaining, less commonly used words, are discarded.\n",
    "\n",
    "At this point, any prosaic resemblance left in the data is gone and now we are dealing entirely with textual data in a numeric form."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The vector representation of The croakers by Drake, Joseph Rodman,\n",
      "The vector has a length of 3200.\n",
      "The first 100 elements of the vector:\n",
      "[570. 254. 194. 167. 136. 123.  95.  86.  40.  56.  48.  61.  76.  89.\n",
      "  40.  30.  27.  50.  59.  63.  12.  20.  26.  32.  54.  34.  18.  13.\n",
      "  20.  44.  33.  26.  15.  41.  24.  16.  12.   9.  32.  37.  32.  17.\n",
      "  21.  39.  22.   3.   8.  26.  11.  18.  15.  26.  13.  21.  72.   7.\n",
      "  15.  21.   7.  14.  15.  47.  25.  36.  24.  11.   9.  43.   8.  13.\n",
      "   9.   9.  11.   4.  16.  13.  13.  13.   9.  10.   8.  20.   7.   9.\n",
      "   8.  21.  26.  19.  14.  25.  12.   6.   5.   6.  13.   6.   4.   9.\n",
      "   9.  17.]\n"
     ]
    }
   ],
   "source": [
    "### DEFACTORING INSPECTION\n",
    "print(\"The vector representation of {} by {}\".format(\n",
    "    metadict[defactoring_volume_id]['title'], \n",
    "    metadict[defactoring_volume_id]['author']))\n",
    "print(\"The vector has a length of {}.\".format(\n",
    "    len(features)))\n",
    "print(\"The first 100 elements of the vector:\")\n",
    "print(features[0:100])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The inspection above shows us the last volume processed by the loop, *The Croakers* by Joseph Rodman Drake. As we can see, the words for this volume of poetry are now represented as a list of numbers (representing word frequencies). However, this list of numbers still requires additional transformation in order to be consumable by logistic regression. The word frequencies need to be normalized so they are comparable across volumes. To do this Underwood and Sellers divide the frequency of each individual word (each number in the list above) by the total number of words in that volume (the `totalcount` variable. This makes volumes of different lengths comparable by turning absolute frequencies into relative frequencies. One thing we initially did not understand is why the value of 0.001 has been added to the `totalcount` variable. When we asked, it turned out this is a \"lazy\" way to prevent divide-by-zero errors.\n",
    "\n",
    "The end result of the code we have executed thus far in the notebook is a very neat and tidy table of numbers between zero and 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "      <th>5</th>\n",
       "      <th>6</th>\n",
       "      <th>7</th>\n",
       "      <th>8</th>\n",
       "      <th>9</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.058291</td>\n",
       "      <td>0.044428</td>\n",
       "      <td>0.027460</td>\n",
       "      <td>0.017129</td>\n",
       "      <td>0.015469</td>\n",
       "      <td>0.013382</td>\n",
       "      <td>0.011455</td>\n",
       "      <td>0.007440</td>\n",
       "      <td>0.006209</td>\n",
       "      <td>0.006049</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.048901</td>\n",
       "      <td>0.028083</td>\n",
       "      <td>0.022968</td>\n",
       "      <td>0.015976</td>\n",
       "      <td>0.011790</td>\n",
       "      <td>0.012650</td>\n",
       "      <td>0.008961</td>\n",
       "      <td>0.009504</td>\n",
       "      <td>0.004118</td>\n",
       "      <td>0.005137</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.046204</td>\n",
       "      <td>0.030856</td>\n",
       "      <td>0.019471</td>\n",
       "      <td>0.015614</td>\n",
       "      <td>0.014151</td>\n",
       "      <td>0.011757</td>\n",
       "      <td>0.007980</td>\n",
       "      <td>0.007820</td>\n",
       "      <td>0.003777</td>\n",
       "      <td>0.006969</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.057219</td>\n",
       "      <td>0.036929</td>\n",
       "      <td>0.023076</td>\n",
       "      <td>0.014638</td>\n",
       "      <td>0.013206</td>\n",
       "      <td>0.013127</td>\n",
       "      <td>0.010086</td>\n",
       "      <td>0.008340</td>\n",
       "      <td>0.006436</td>\n",
       "      <td>0.004592</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.052128</td>\n",
       "      <td>0.029098</td>\n",
       "      <td>0.016521</td>\n",
       "      <td>0.014094</td>\n",
       "      <td>0.013653</td>\n",
       "      <td>0.013404</td>\n",
       "      <td>0.007116</td>\n",
       "      <td>0.006592</td>\n",
       "      <td>0.003448</td>\n",
       "      <td>0.004496</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>715</th>\n",
       "      <td>0.048584</td>\n",
       "      <td>0.031039</td>\n",
       "      <td>0.023947</td>\n",
       "      <td>0.013898</td>\n",
       "      <td>0.009401</td>\n",
       "      <td>0.009725</td>\n",
       "      <td>0.010292</td>\n",
       "      <td>0.007861</td>\n",
       "      <td>0.004538</td>\n",
       "      <td>0.004295</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>716</th>\n",
       "      <td>0.041081</td>\n",
       "      <td>0.032973</td>\n",
       "      <td>0.017568</td>\n",
       "      <td>0.017297</td>\n",
       "      <td>0.012973</td>\n",
       "      <td>0.024865</td>\n",
       "      <td>0.009009</td>\n",
       "      <td>0.007477</td>\n",
       "      <td>0.004595</td>\n",
       "      <td>0.005405</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>717</th>\n",
       "      <td>0.036274</td>\n",
       "      <td>0.036804</td>\n",
       "      <td>0.017842</td>\n",
       "      <td>0.020672</td>\n",
       "      <td>0.010827</td>\n",
       "      <td>0.012596</td>\n",
       "      <td>0.009373</td>\n",
       "      <td>0.007172</td>\n",
       "      <td>0.006327</td>\n",
       "      <td>0.010493</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>718</th>\n",
       "      <td>0.058199</td>\n",
       "      <td>0.029858</td>\n",
       "      <td>0.029765</td>\n",
       "      <td>0.015161</td>\n",
       "      <td>0.013057</td>\n",
       "      <td>0.018038</td>\n",
       "      <td>0.008632</td>\n",
       "      <td>0.009715</td>\n",
       "      <td>0.003558</td>\n",
       "      <td>0.004115</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>719</th>\n",
       "      <td>0.059086</td>\n",
       "      <td>0.026329</td>\n",
       "      <td>0.020110</td>\n",
       "      <td>0.017311</td>\n",
       "      <td>0.014098</td>\n",
       "      <td>0.012750</td>\n",
       "      <td>0.009848</td>\n",
       "      <td>0.008915</td>\n",
       "      <td>0.004146</td>\n",
       "      <td>0.005805</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>720 rows × 10 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "            0         1         2         3         4         5         6  \\\n",
       "0    0.058291  0.044428  0.027460  0.017129  0.015469  0.013382  0.011455   \n",
       "1    0.048901  0.028083  0.022968  0.015976  0.011790  0.012650  0.008961   \n",
       "2    0.046204  0.030856  0.019471  0.015614  0.014151  0.011757  0.007980   \n",
       "3    0.057219  0.036929  0.023076  0.014638  0.013206  0.013127  0.010086   \n",
       "4    0.052128  0.029098  0.016521  0.014094  0.013653  0.013404  0.007116   \n",
       "..        ...       ...       ...       ...       ...       ...       ...   \n",
       "715  0.048584  0.031039  0.023947  0.013898  0.009401  0.009725  0.010292   \n",
       "716  0.041081  0.032973  0.017568  0.017297  0.012973  0.024865  0.009009   \n",
       "717  0.036274  0.036804  0.017842  0.020672  0.010827  0.012596  0.009373   \n",
       "718  0.058199  0.029858  0.029765  0.015161  0.013057  0.018038  0.008632   \n",
       "719  0.059086  0.026329  0.020110  0.017311  0.014098  0.012750  0.009848   \n",
       "\n",
       "            7         8         9  \n",
       "0    0.007440  0.006209  0.006049  \n",
       "1    0.009504  0.004118  0.005137  \n",
       "2    0.007820  0.003777  0.006969  \n",
       "3    0.008340  0.006436  0.004592  \n",
       "4    0.006592  0.003448  0.004496  \n",
       "..        ...       ...       ...  \n",
       "715  0.007861  0.004538  0.004295  \n",
       "716  0.007477  0.004595  0.005405  \n",
       "717  0.007172  0.006327  0.010493  \n",
       "718  0.009715  0.003558  0.004115  \n",
       "719  0.008915  0.004146  0.005805  \n",
       "\n",
       "[720 rows x 10 columns]"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "### DEFACTORING INSPECTION\n",
    "# Normalized perspective of the data\n",
    "data.iloc[:,0:10]# display first and last 5 rows and 10 columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last row in that table, 719, is the volume we have been tracking, *The Croakers* by Joseph Rodman Drake. It is just one of 720 relatively indistinguishable rows of numbers in this representation of 19th century poetry. This is a radical transformation of the original, prosaic representation literary historians are probably used to seeing (if the poem was of interest in the first place):\n",
    "\n",
    "![Screenshot of the Google Books site for The Croakers](notebook_resources/the-croakers.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "sextuplets = list()\n",
    "for i, volid in enumerate(orderedIDs):\n",
    "    listtoexclude = authormatches[i]\n",
    "    asixtuple = (data, \n",
    "                 classvector, \n",
    "                 listtoexclude, \n",
    "                 i, \n",
    "                 usedate, \n",
    "                 regularization)\n",
    "    sextuplets.append(asixtuple)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the last step before Underwood and Sellers's code moves away from the transformation of features and into the actual analysis of the data. This bit of code gathers all of the relevant data and metadata that has been cleaned and normalized in a structure suitable for performing the statistical analysis. The `sextuplets` variable is a list of 720 tuples containing six elements. Each item in the `sextuplets` list contains the necessary data structures to model each poem. The contents of each item in the list is as follows:\n",
    "\n",
    "- `data`: a normalized feature matrix. Word features are the columns and volumes are the rows with dimensions of 720 x 3200.\n",
    "- `classvector`: the classification or labels of volumes as either 'reviewed' (1) or 'random' (0).\n",
    "- `listtoexclude`: the list of poems to ignore because they are the same author.\n",
    "- ``i: the index of the volume\n",
    "- `usedate`: a flag indicating if date is a feature. It is false in this analysis.\n",
    "- `regularization`: a parameter for the logistic regression. This value was hardcoded at the beginning of the code in the Setting Parameters section.\n",
    "\n",
    "With all of the data assembled and in the right shape, a process we call *data fitness*, we can now venture into the algorithmic territory and perform the statistical analysis. As we can see, the *fitted* representation of features has traveled a great distance from the original poetry. One of the most important aspects of *distant reading* is the work of cleaning, preparing, and normalizing texts to be \"read\" by an algorithm. When considering distance, we should think not only of the perspective that we, the analyst, are reading from, but also the *distance traveled* in terms of successive transformations and representations of the data. If computational literary history is a triathalon, we have only completed the first endurance test."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Training Predictive Models\n",
    "\n",
    "We are now about to dive into the very heart of the analysis, training predictive models on each volume. In their MLQ article, Underwood and Sellers summarize this effort:\n",
    "\n",
    "> At the grittiest mathematical level, the predictive model we create is just an equation that translates word frequencies into a probability that a particular volume came from the reviewed sample. Although we may say that it models a boundary between the samples, probability is a continuum, and a probabilistic model allows us to treat social boundaries as fuzzy gradients. We \"train\" the model by showing a regression algorithm the volumes from all authors (except one) in both samples; the algorithm assigns each word a positive or negative weight in an effort to separate the samples. When we show the model a new volume, by the author it hasn't yet seen, it uses the weights assigned to different words to estimate the probability that this volume was reviewed. To avoid circularity, the model never makes predictions about one of an author's books by using information about others (Sculley and Pasanek 2008). So we actually train 636 slightly different models, each one excluding books by a different author. However, since each pair of models shares more than 99 percent of their evidence, we can describe them collectively as one model of the literary field (Underwood & Sellers 2016, 326).\n",
    "\n",
    "We can see a more *Pythonic* expression of this work below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "### DEFACTORING DEFINITION\n",
    "\n",
    "def sliceframe(dataframe, yvals, excludedrows, testrow):\n",
    "    numrows = len(dataframe)\n",
    "    newyvals = list(yvals)\n",
    "    for i in excludedrows:\n",
    "        del newyvals[i]\n",
    "        # NB: This only works if we assume that excluded rows\n",
    "        # has already been sorted in descending order !!!!!!!\n",
    "        # otherwise indexes will slide around as you delete\n",
    "\n",
    "    trainingset = dataframe.drop(dataframe.index[excludedrows])\n",
    "\n",
    "    newyvals = np.array(newyvals)\n",
    "    testset = dataframe.iloc[testrow]\n",
    "\n",
    "    return trainingset, newyvals, testset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function prepares the data for training one model by separating volumes used for training from a single volume held out to test the predictive power of the model. The function takes a dataframe containing the feature vectors, a list of the classifications for each volume, a list of volumes to exclude (because of shared authorship), and the the index of the specific volume to be held out. It returns the dataframe with the held out volume removed (`trainingset`), a list of the known classifications (`newyvals`) corresponding to the training set, and the held-out volume that will be classified once the model has been trained (`testset`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "### DEFACTORING DEFINITION\n",
    "\n",
    "def normalizearray(featurearray, usedate):\n",
    "    '''Normalizes an array by centering on means and\n",
    "    scaling by standard deviations. Also returns the\n",
    "    means and standard deviations for features.\n",
    "    '''\n",
    "\n",
    "    numinstances, numfeatures = featurearray.shape\n",
    "    means = list()\n",
    "    stdevs = list()\n",
    "    lastcolumn = numfeatures - 1\n",
    "    for featureidx in range(numfeatures):\n",
    "\n",
    "        thiscolumn = featurearray.iloc[ : , featureidx]\n",
    "        thismean = np.mean(thiscolumn)\n",
    "\n",
    "        thisstdev = np.std(thiscolumn)\n",
    "\n",
    "        if (not usedate) or featureidx != lastcolumn:\n",
    "            # If we're using date we don't normalize the last column.\n",
    "            means.append(thismean)\n",
    "            stdevs.append(thisstdev)\n",
    "            featurearray.iloc[ : , featureidx] = \\\n",
    "                (thiscolumn - thismean) / thisstdev\n",
    "        else:\n",
    "            print('FLAG')\n",
    "            means.append(thismean)\n",
    "            thisstdev = 0.1\n",
    "            stdevs.append(thisstdev)\n",
    "            featurearray.iloc[ : , featureidx] = \\\n",
    "                (thiscolumn - thismean) / thisstdev\n",
    "            # We set a small stdev for date.\n",
    "\n",
    "    return featurearray, means, stdevs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function standardizes the features by computing the z-score for the feature vectors. That is, it loops over each column of the data, subtracts the column mean from each value, and then divides that value by the standard deviation. This is an important step in the data preparation pipeline because it ensures all of the data values are on the same scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "### DEFACTORING DEFINITION\n",
    "\n",
    "def model_one_volume(data5tuple):\n",
    "    data, classvector, listtoexclude, i, usedate, regularization = \\\n",
    "        data5tuple\n",
    "    trainingset, yvals, testset = sliceframe(data, \n",
    "                                             classvector, \n",
    "                                             listtoexclude, \n",
    "                                             i)\n",
    "    newmodel = LogisticRegression(C = regularization)\n",
    "    trainingset, means, stdevs = normalizearray(trainingset, usedate)\n",
    "    newmodel.fit(trainingset, yvals)\n",
    "\n",
    "    testset = (testset - means) / stdevs\n",
    "    prediction = newmodel.predict_proba(testset.values.reshape(1, -1))[0][1]\n",
    "    if i % 50 == 0:\n",
    "        print(i)\n",
    "    # print(str(i) + \"  -  \" + str(len(listtoexclude)))\n",
    "    return prediction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In many respects, this is the most salient block of code in the entire document. The code above actually runs the logistic regression and does the computational work that generates a prediction about each individual volume. This function builds upon the two previous functions to assemble a normalized set of training data (`trainingset`) distinct from the single volume to be predicted (`testset`).\n",
    "\n",
    "There are three lines of code involved in the computational modeling of the data. First, Underwood and Sellers instantiate a model object with the regularization parameter (more on that below):\n",
    "```\n",
    "newmodel = LogisticRegression(C = regularization)\n",
    "```\n",
    "Then they fit the model using the normalized training data:\n",
    "```\n",
    "newmodel.fit(trainingset, yvals)\n",
    "```\n",
    "Once a model has been fit they can use that model to make predictions about unseen or held-out data. This is what they do with the `predict_proba()` function:\n",
    "```\n",
    "prediction = newmodel.predict_proba(testset.reshape(1, -1))[0][1]\n",
    "```\n",
    "Those three lines are all it takes to do the computational part of of the analysis, the rest of the code up until this point has all been data preparation, cleaning, normalization, and re-shaping. This ratio of standardized function calls to bespoke data management code indicates that claims of machines are eradicating scholars' jobs are greatly exaggerated (Basken 2017)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Regularization and Logistic Regression\n",
    "\n",
    "The three lines of code above hide a significant amount of intellectual and computational work. The call to the `newmodel.fit()` function is a crucial step in the analytical process. Underwood and Sellers are using an implementation of logistic regression from the 3rd party Python library scikit-learn.\n",
    "\n",
    "At a very high level, logistic regression is a machine learning algorithm for performing classification. Logistic regression works by estimating the parameters of a function, the *hypothesis representation*, that divides a multidimensional space into two parts (note, in this case we are talking about binomial or binary logistic regression, which classifies things into one of two bins). The hypothesis representation describes a line that winds its way through the space creating what is called the *decision boundary*. Every data point that lands on one side the boundary gets one label and every data point on the other side of the boundary gets the other label. Similar to linear regression, the goal is to find the best hypothesis representation, that is, the function that best draws a line dividing the space given the known data points. Once you have a good hypothesis representation, an appropriately *fit* model, you can begin to classify *new* data by dropping data points into the multidimensional space and seeing on which side of the decision boundary they land.\n",
    "\n",
    "The key to logistic regression is estimating the parameters of the hypothesis representation-the parameters to the function that draws a line through the multidimensional space. We can derive the parameters by using the *features* of existing data combined with their known labels; this is called *training data*. The modeling process is executed by the function call to `newmodel.fit(trainingset, yvals)`. In Underwood and Sellers’s c'de the model uses the training data---the matrix of word features in the data variable and known labels ('reviewed' or 'random') in the `classvector` variable---to \"learn\" the parameters through a process called *coordinate descent*. How so-called optimization functions  work is well beyond the scope of the discussion."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Overfitting\n",
    "\n",
    "One of the problems when fitting a logistic regression model is a tendency towards *overfitting*. Crudely this means the model, the function with the learned parameters, that you estimated have tailored themselves such that they are overly optimized to the particular training data you provided. As such, the model becomes less useful for prediction or classifying new data because they are outside the fitness of the model. An overfit model is like a snug pair of jeans, once you put on a few pounds (add new data) they don't fit. In Underwood and Sellers's case, they are fitting models on all volumes except one, which is held out. Then they test the predictive performance of the model by seeing if it correctly classifies the held-out volume. If they overfit the models, the model will to a terrible job guessing the status of the held out volumes.\n",
    "\n",
    "When Underwood and Sellers instantiated the model (`newmodel = LogisticRegression(C = regularization)`), they set a regularization parameter on the model. Regularization is a technique for logistic regression (and other machine learning algorithms) that smooths out the tendency toward overfitting with some more mathematical gymnastics. The diagram below shows how regularization can help with the fitness of the model:\n",
    "\n",
    "![Regularization](notebook_resources/regression_figures.png)\n",
    "*On the left side is a linear regression which doesn't quite fit the data. In the middle is an overfit logistic regression. On right side is a regularized logistic regression.*\n",
    "\n",
    "As the diagrams show, the regularized logistic expression (the right side) does have a bit of error, there are pink and blue dots on the wrong sides of the decision boundary, but as more data get added it will generally be more right than the overfitted model as represented by the middle diagram (the squiggly decision boundary)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Running the Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Beginning multiprocessing.\n",
      "0\n",
      "50\n",
      "100\n",
      "150\n",
      "200\n",
      "250\n",
      "300\n",
      "350\n",
      "450\n",
      "500\n",
      "400\n",
      "550\n",
      "600\n",
      "650\n",
      "700\n",
      "Multiprocessing concluded.\n"
     ]
    }
   ],
   "source": [
    "# Now do leave-one-out predictions.\n",
    "print('Beginning multiprocessing.')\n",
    "\n",
    "pool = Pool(processes = 4)\n",
    "res = pool.map_async(model_one_volume, sextuplets)\n",
    "\n",
    "# After all files are processed, write metadata, errorlog, and counts of phrases.\n",
    "res.wait()\n",
    "resultlist = res.get()\n",
    "\n",
    "assert len(resultlist) == len(orderedIDs)\n",
    "\n",
    "logisticpredictions = dict()\n",
    "for i, volid in enumerate(orderedIDs):\n",
    "    logisticpredictions[volid] = resultlist[i]\n",
    "\n",
    "pool.close()\n",
    "pool.join()\n",
    "\n",
    "print('Multiprocessing concluded.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This code automates the training of 720 volumes, holding out one volume, training the model on the remaining volumes, and then making a prediction for the held-out volume. As a very computation intensive process training or fitting a logistic regression model takes time, and training 720 different models obviously takes 720 times longer. Fortunately, this is a so called embarrassingly parallel computational task and so we can train the models using parallel processing instead of one after the other. Using Python's built in parallel processing modules, this code can speed up the process. Still, this block of code takes a bit of time to execute, around twenty minutes on a quad core MacBook Pro (late 2013 model). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 720 predictions.\n"
     ]
    }
   ],
   "source": [
    "### DEFACTORING INSPECTION\n",
    "print(\"There are {} predictions.\".format(len(logisticpredictions)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What emerges from the other side of this computationally intensive task are a series of predictions, 720 to be specific, one for each of the modeled volumes. These predictions, stored in the `logisticpredictions` variable, are the model's assertions of each volume's reviewed status. Additionally, because we already know the status of the modeled volumes we can compare the performance of the predictive model to the \"ground truth\" and see if the algorithm was able to detect patterns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "truepositives = 0\n",
    "truenegatives = 0\n",
    "falsepositives = 0\n",
    "falsenegatives = 0\n",
    "allvolumes = list()\n",
    "\n",
    "with open(outputpath, mode = 'w', encoding = 'utf-8') as f:\n",
    "    writer = csv.writer(f)\n",
    "    header = ['volid', \n",
    "              'reviewed', \n",
    "              'obscure', \n",
    "              'pubdate', \n",
    "              'birthdate', \n",
    "              'gender', \n",
    "              'nation', \n",
    "              'allwords', \n",
    "              'logistic', \n",
    "              'author', \n",
    "              'title', \n",
    "              'pubname', \n",
    "              'actually', \n",
    "              'realclass']\n",
    "    writer.writerow(header)\n",
    "    for volid in IDsToUse:\n",
    "        metadata = metadict[volid]\n",
    "        reviewed = metadata['reviewed']\n",
    "        obscure = metadata['obscure']\n",
    "        pubdate = infer_date(metadata, datetype)\n",
    "        birthdate = metadata['birthdate']\n",
    "        gender = metadata['gender']\n",
    "        nation = metadata['nation']\n",
    "        author = metadata['author']\n",
    "        title = metadata['title']\n",
    "        canonicity = metadata['canonicity']\n",
    "        pubname = metadata['pubname']\n",
    "        allwords = volsizes[volid]\n",
    "        logistic = logisticpredictions[volid]\n",
    "        realclass = classdictionary[volid]\n",
    "        outrow = [volid, \n",
    "                  reviewed, \n",
    "                  obscure, \n",
    "                  pubdate,\n",
    "                  birthdate, \n",
    "                  gender, \n",
    "                  nation, \n",
    "                  allwords, \n",
    "                  logistic, \n",
    "                  author, \n",
    "                  title, \n",
    "                  pubname, \n",
    "                  canonicity, \n",
    "                  realclass]\n",
    "        writer.writerow(outrow)\n",
    "        allvolumes.append(outrow)\n",
    "\n",
    "        if logistic > 0.5 and classdictionary[volid] > 0.5:\n",
    "            truepositives += 1\n",
    "        elif logistic <= 0.5 and classdictionary[volid] < 0.5:\n",
    "            truenegatives += 1\n",
    "        elif logistic <= 0.5 and classdictionary[volid] > 0.5:\n",
    "            falsenegatives += 1\n",
    "        elif logistic > 0.5 and classdictionary[volid] < 0.5:\n",
    "            falsepositives += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This code is a bit simpler than its predecessors. This block writes a CSV file to disk containing 720 rows of volume metadata, the predicted classification, and the actual classification.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Modeling Coefficients\n",
    "\n",
    "The code below generates a single logistic regression model, trained on all of the data with nothing held-out. The properties of this model, the coefficients of the hypothesis representation, are interrogated to better understand the influence of individual features, words, on reviewed or unreviewed volumes. Thus this individual model is not using computational modeling to predict a phenomena, it is using the computational model to explore and explain patterns and features of the phenomena."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "donttrainon.sort(reverse = True)\n",
    "trainingset, yvals, testset = sliceframe(data, \n",
    "                                         classvector, \n",
    "                                         donttrainon, \n",
    "                                         0)\n",
    "newmodel = LogisticRegression(C = regularization)\n",
    "trainingset, means, stdevs = normalizearray(trainingset, usedate)\n",
    "newmodel.fit(trainingset, yvals)\n",
    "\n",
    "coefficients = newmodel.coef_[0] * 100\n",
    "\n",
    "coefficientuples = list(zip(coefficients, \n",
    "                            (coefficients / np.array(stdevs)), \n",
    "                            vocablist + ['pub.date']))\n",
    "coefficientuples.sort()\n",
    "if verbose:\n",
    "    for coefficient, normalizedcoef, word in coefficientuples:\n",
    "        print(word + \" :  \" + str(coefficient))\n",
    "\n",
    "print()\n",
    "accuracy = (truepositives + truenegatives) / len(IDsToUse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This code functions much like the code in the `model_one_volume(`) function except it only trains a single model for the purposes of investigating the impact of particular words on the prediction. By inspecting the magnitude of the coefficients Underwood and Sellers can see how particular words influenced a positive or negative prediction. Looking at the code, specifically the call to `sliceframe()` reveals this model actually does have some hold-out data, the first volume at index zero. We suspect the cost of excluding a single volume is less than the effort of re-implementing the `sliceframe()` function. The code to instantiate and train the model is identical to the code above, but instead of predicting the status of the held-out data the code extracts the coefficients from the model and puts them in the `coefficients` variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 3200 coefficients in this model.\n"
     ]
    }
   ],
   "source": [
    "### DEFACTORING INSPECTION\n",
    "print(\"There are {} coefficients in this model.\".format(len(coefficientuples)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "First ten items in the list of coefficients.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[(-0.4518121129702824, -1531.7794821195973, \"i'll\"),\n",
       " (-0.42713135003767877, -7255.1107491467, 'mission'),\n",
       " (-0.41160936594270203, -862.5421235543723, 'tis'),\n",
       " (-0.37942470885804935, -1125.2718624911568, 'oft'),\n",
       " (-0.37837601538978427, -3337.887638344609, 'greet'),\n",
       " (-0.37397972301466276, -2558.840455180692, 'wondrous'),\n",
       " (-0.3713084261916791, -222.91243757658904, 'our'),\n",
       " (-0.36944484008913314, -4355.541309277771, 'anxious'),\n",
       " (-0.3592339045350108, -2088.196265707296, 'cheer'),\n",
       " (-0.35640117748763345, -1383.7344849648389, \"we'll\")]"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "### DEFACTORING INSPECTION\n",
    "print(\"First ten items in the list of coefficients.\")\n",
    "coefficientuples[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Last ten items in the list of coefficients.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[(0.3727659770624267, 4800.910150171263, 'brows'),\n",
       " (0.3812769536404176, 224.47954001488768, 'not'),\n",
       " (0.38249088903909245, 1149.0081244612627, 'half'),\n",
       " (0.38443007530124773, 398.82882492864155, 'what'),\n",
       " (0.3847771880614351, 6978.657839035639, 'whereon'),\n",
       " (0.39931483467532297, 3172.3564216702794, 'sake'),\n",
       " (0.4030191939031366, 2023.1474550567532, 'slow'),\n",
       " (0.4067213474693706, 3933.656549634385, 'hollow'),\n",
       " (0.4125900276423275, 1531.8208421175516, 'black'),\n",
       " (0.45244212603563333, 614.9111003777699, 'eyes')]"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "### DEFACTORING INSPECTION\n",
    "print(\"Last ten items in the list of coefficients.\")\n",
    "coefficientuples[-10:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The coefficients are a list of numbers, one per word feature, that determine the shape of the line through the multidimensional space. These results tell us the influence of particular words in the classification of volumes in the corpus. Looking at the normalized values helps us understand the degree to which particular words are more or less associated with reviewed or unreviewed poetry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "coefficientpath = outputpath.replace('.csv', '.coefs.csv')\n",
    "with open(coefficientpath, mode = 'w', encoding = 'utf-8') as f:\n",
    "    writer = csv.writer(f)\n",
    "    for triple in coefficientuples:\n",
    "        coef, normalizedcoef, word = triple\n",
    "        writer.writerow([word, coef, normalizedcoef])\n",
    "\n",
    "\n",
    "### DEFACTORING NAMESPACE\n",
    "rawaccuracy = accuracy "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This code generates the `mainmodelcoefficients.csv` output file, which contains the word, its coefficient, and its normalized coefficients."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting Results\n",
    "\n",
    "The final function of the analysis is to test the accuracy of the model(s). This code produces a plot giving a sense of how the model performed compared to the known classifications."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAABX0klEQVR4nO2dfXwU1b3/P5NNIBsIJCQlkA1pKbY+4EVpKQXqQykPt7a1RSjty1qfkIuiIk0CEiCQkBiIIqGKilqotP64Xh+I2uttbYQqtDwI1YsKXh+iKLBEKImBQDYh2Z3fH8ksZ2bPmTmzM7M72T3v1ysvyGQezjlz5nvO+Z7vgyTLsgyBQCAQJAUp8S6AQCAQCGKHEPoCgUCQRAihLxAIBEmEEPoCgUCQRAihLxAIBEmEEPoCgUCQRKTG8+HHjh1z9P65ubk4efKko89wG8lYZyA56y3qnDyQ9c7Pz7d0LzHTFwgEgiRCCH2BQCBIIoTQFwgEUeGtq8PgsWMxtKAAg8eOhbeuLt5F6lXEq/3iqtMXCAS9E29dHQbeey9SAgEAQKrfj4H33gsACEyfHs+i9Qri2X5cQn///v146qmnEAqFMGnSJEybNk3197a2Njz88MNoampCMBjEtddei4kTJzpRXoFA4AIya2rCAkshJRBAZk2NEPocxLP9DIV+KBTCxo0bUVZWhpycHCxevBhjxoxBQUFB+JxXX30VBQUFKC0txenTpzF//nxceeWVSE0VCwmBIBHxMCzvWMcFauLZfoY6/YaGBgwZMgR5eXlITU3FhAkTsG/fPtU5kiShvb0dsiyjvb0d/fv3R0qK2C4QCBKVIMNskHVcoCae7Wc4FW9ubkZOTk7495ycHHz88ceqc374wx/igQcewO23345AIICioiKq0N+6dSu2bt0KAKipqUFubq7V8uuSmprq+DPcRjLWGUjOese1ztXVkO+8E1JbW/iQnJEBVFc7WqaEec8m28/OehsKfVq4fUmSVL+/8847+OpXv4rly5fj+PHjqKqqwkUXXYSMjAzVeZMnT8bkyZPDvzvtZJGMjhzJWGcgOesd1zpPmQLv/fcjs6YGnmPHEMzPR2tpKQJTpgAOlilh3rPJ9rPTOctQ6Ofk5KCpqSn8e1NTE7Kzs1XnvP7665g2bRokScKQIUMwePBgHDt2DBdccIGlwgkEAvcSmD5dbNpaIF7tZ6h4HzFiBBobG3HixAl0dXVh165dGDNmjOqc3NxcvPfeewCAlpYWHDt2DIMHD3amxAKBQCCIGsOZvsfjwaxZs1BdXY1QKISJEydi2LBhqK+vBwBMnToVM2bMwGOPPYaSkhIAwA033IABAwY4W3KBQCAQmEaKZ45cEXDNfpKxzkBy1lvUOXkQAdcEAoFAEBVC6AsEAkESIYS+QCAQJBFC6AsEAkESIYS+QCAQxAg3hKMWEdEEAoEgBrglHLWY6QsEAkEM0AunHEuE0BcIBIIY4JZw1ELoCwQCQQxwSzhqIfQFAoEgBrSWliLk9aqOhbxetJaWxrQcYiNXIBAIYoCyWRsRTjnGkTbFTF8gEAiiQDG/fMLnU5lf6pllBqZPx4m9e9F49ChO7N0bl9DKYqYvEAgEJiHNLysBrOgxv0zbtw8Zzz8fd7NMPcRMXyAQxBU3OCyZhWV+2W/zZleYZeohhL5AIIgbyow51e+HJMvhmbHbBX+V3w8JgJI4Vvl/RTBIPT/WZpl6CKEvEAjihlsclsyyzOeDDEBJRqL8v8LjoZ5PM8uM1wpHCH2BQBATaELOLQ5LZmGZX5694QamWeaaNWvCx+K5whFCXyAQOA5LyIWysqjnx9phySyB6dNx6oEH0OXzYTmALp8Ppx54AKdXrQoflyUpfDwwfTpqa2vD18dzhSOsdwQCgeOwhFywb1+EvF7V3+LhsBQNgenTEZg+HbcDOEE5rkc8Vzhipi8QWCQa3SzPNb3RqoUFS5ilnDrFnBm7HZ73s2bNGvh8Pvh8PgAI/7+8f3/qPWOxwhFCXyCwQDS6WZ5reqtVCwu9uDMshyU3D3q876ekpAR+vx9+vx8Awv8vWrkybiEZhNAXCCwQjW6W5xrWOWuXLLGh1LGHtvEpSxI8fj9VoLOEasozz8Sy2Eys6uTJPYFYr3C4hP7+/fsxf/58zJs3Dy+99FLE3//0pz9h4cKFWLhwIUpKSvDLX/4SZ86csbusgihw82zJTUTbTtHoZnmuYZ1T1drKVS63oRJy6Bb4kixDAqizZJZQ9SxfHtuCM/D0zNx5jwNAcXGx6ndyhdNaWorMmpqIkA5OYCj0Q6EQNm7ciCVLlmDt2rXYuXMnjh49qjrnpz/9KVavXo3Vq1fj+uuvxyWXXIL+DJ2VIHYkmorAKaJV0QweOxaQZerf9XSzPCF23W69Eg2KkAv6fJA07aadJTMHzSNHnCwiPwx7fNpxpa88uHat4aqmEvRB0E4MhX5DQwOGDBmCvLw8pKamYsKECdi3bx/z/J07d+J73/uerYUUREdvdXyJNWbbSTVIUP4eSkvT1c3yhNglz6kAVN6fymZgVVVVRLl6w6qOZ6XDHPSGDXOiSOZheN4iGIwIvmY0oYj1d2postnc3IycnJzw7zk5Ofj444+p53Z0dGD//v247bbbqH/funUrtm7dCgCoqalBbm5uNGXmJjU11fFnuA2yznofV6K1i5V3bbad0lavhqT5SEkkSUJmZib65eYi5ZlnulUSR44Aw4YhWFmJ0Jw5CGVmQiKOhyor0e/669FPuQlxTvmRIyjvuTbtllvQ0dERrnPo6ae773/4MNCjMgG6Z4tZixYhMzMToeuvj6pdHGPYsO7yUo6H27u6GvKdd0Jqawv/Wc7IgFxdHdV7pr4HK+1SWEitg6KuUtreQ+krKYEAslavRr85cwB0h3So1NwDAJb7/VjaU1c7ZZmh0Jcpy1dJos1vgLfeegsXXnghU7UzefJkTJ48Ofz7yZMnecsZFbm5uY4/w22QdR6cn49Uio4xmJ+fcO1i5V2bbaehBioG6dw5YOlSnG1txcB77z3/0R8+jJS5c9Ha2tq9YTdlivpC7bOmTIk8B+e/m8fvvBMV9fXn76/5VqW2NmDpUpyk3COeeBcuVCUIB7pXOqcWLkRAaYMpU+C9//6I2PP9fvELnDx5Et66uvDfQllZgCwj5dQpaox6ZbbNfA821YFEamuD55Zb2Dc4ciT8Hpf5fFjR0/8knA/t0OXz4UTPOWT/zreo+jNU7+Tk5KCpqSn8e1NTE7Kzs6nn7ty5E1dccYWlAgnswy2ZetyO2Xbi0bd7jh1zZNlObgZWvfwyU+iEy+H3I2/kSORdeqlr1D68lit6ppykysTz5ZfwtLScV58UFzuuPtFuTNMgVXJayD4U6+/UUOiPGDECjY2NOHHiBLq6urBr1y6MGTMm4ry2tja8//771L8J4kM8zcJ6E2bbifaRagnm5zvidVlSUnJ+E5kDCYCnpQWeL7901Wa+Fdt8mhAnSensxIBly8K/O+X9qtSBFUqChSLQlbpm3XMPkJICWZKwHIDs8aBt5kzHvlND9Y7H48GsWbNQXV2NUCiEiRMnYtiwYaivrwcATJ06FQCwd+9eXHbZZUhPT3ekoILo4HEJF5hrJ1XaO79fpUsHztufw+OhbvgF8/NV6gkzafMemj0bD/zlL+HflZlkObo3fHkgZ7lkGdonTUL6tm3cZeKpg5l6kolJAHYCEj2zyHAdW1rC/w/qqO9sgaHuJgn3DkkCUlKQNW+eqt9IZ88CAFYAQDCIjOefR+d3vuPItyvJNKV9jDjmcJyJZNfpJxPxrHdYsNEGAKiX+CGvF20zZ6qyKynHeVZhg8eODQswUv8bfp4kAT3273rIAGRNzBtaWVll0gpo2vk857DqRtLl8+HE3r0Aut9zWkYGJJb1DFGXlnXrmO+Ft731UN57ld/fLaxtRlvvmOn0BQIn6S1mhnro2Z9L6F6uk2qj9G3botYxs1QScs9zzt54I4I9cV70b+SJKIN2oNArkxWvYtY9udUwBgIfAOR+/VRmtZIsQ5YkyIAtak6tbb0TOBV8TQh9QdyIt/OYlQHHTGx4hEIq3bUVHTOpkignjq8AIPWoBdonTdLdcwh5vVyCU69MVryKWcd5nNYAGA5qcloaKrq6Igc1WUawZ/ZsVW1itK9gB0456AmhL4gJNCEZT+cxKwMOMzb8wIHU8yOEFqdwo6F12lJQ1AspgQDSt21TbUwHs7IQzM5WrTa4VgM4v/+gfXdWvIpZx3mtWFhxfJRZfEttLap6fBm0sAYcMsEJD8x0iabuwiau1jsCgVVYQpIZv8SGZS1rFh+2mJg3L+oBhzVYQZIiBDLt422fNClSF99z3AgeU0HPsWMq65jjBw/i+IEDqtUGVXBq7hPyetE+aRL13dFWE3pexaxzmHXTWFEp7y0tPR2ZNTVomzlTdV7Lww9jQXGx4SyeNeCQCU54qPB46OkSTd3l/GAVzM5GMCsrJlZ2YiM3wXBjnVkbdLLHQ92QIzeweCHrzdpApG2gRpRJktCoiS2lZWhBQYTuXrm25eGHkVlTgzS/H50+H9VahWfDkocnLr4YladPRxxflpmJOz74wPB6rWUNzXons6aGWVbl71V+P5Yx6hqtlZL2Hjwbwj7G6kWxbNLbvPX5fGhWNn45yvqEz0fV5fNaUcnoVlPxtoedG7lC6CcYbqwzU0gi0oKEx6qCJkj6zZlz3hPZ5CCjLZPRx8gjtH0+XziGuha9QcNowFGV47XXkDJ3LlICgbAljxWrlAGLF6Pf5s3d+n6PB2dvuAH9nn7asKxagRkaOBCQJKS0tJgW9LR3qzfwkIMk2eY8QnzNmjXUGT7PIGFkRWWE2QFeCH1O3CgAncaNddYTkspHbcY2nDbrC61fjxM94Qb0Bhlji+rz9zRrrrj4+99X2dArFBcXo6SkJPy7XTP93NxcnH3yScOVBQ8DFi9Gvz/+UdU+Mro3RVM6OyPOD2Zl4fjBgwC6BWxQM3iT8A5ErHaVegY1LbIkYUFREVM1wxp0afgY6jLWOyHLyhL6MhC2yycnG4pZbbxm+kKnn+C4wSRST4fN8sxkwRNnnbkhygiHS/tg9fT7LN3z/A0bqFmSSIEP2Ot2r7RfMYc+W49+mzdHCFYJgEQR+ACwor1dlQbQ0yP8KijnWt4rYby3YH4+MzOVNnY9CyPvZtb+kjYxukxx0AqbiqakdG+kg55HYObVV3OV1S6E0E9g4m0SqZC+bRtVoKRv22b6Xjxx1llC9ewNN1CPm34WjMMIAGAOsoHp09E2c2a3/T7scbvXDiw8kBMCXhNOhRXt7Qh6vdwbmTyb88xzgkHTgyRPe5DfRznjHD2LKqUP3O73o+Xhh5mb6ymdnZAzMph5BHY1NIR/N2tFFA1C6Ccwbomnb2fsE54466yZ+OlVq6jHWeaLZu2ktUKENch66+qQ8fzzkILB7tlgj3299rwBixdjaGEhhvp8GFpYiAGLF5sqD6uMg8eOxVCfD1n33BOeEPCqvcJQnLv00LalKTPQnvdkFBupuCfQGu/Klvw+Kih/NxpYyGcNWLZMFQZai+fYMZW1Gu15gHkromgQOv0Eg6yzXRuGVrFLhw3w6fRp1+jtG5gNF8BCr56LfvnL8OyTpz1YOvazN92E06tWATDfv2n11EMG0HHFFejz1lu6evYK6FusKDrsUHZ2dwjklhZqWASe8BR6dTb7HvX2foz07Wbbssvng+eLL8K6fb1BlrYXIXT6Ai6sOAHZid06bNqsj5UQg0fFFU00UjMeuZ5jx1QzOJ6VD0vH3m/zZmaZjODxIlXUNLLHg7M33YTmZ581XB1VaK4PO4NBrcMOh0AGqGoOrVOZWXt11so2a9486qxfb2XB2h/R8/NgEe7rGhUaaedPouyVOKXqETP9BIPHXj0e4ZXtsNfWg/Wu7VxlKLDaVU5Ph+fLL6nPSiM2G3nKNNTno1usAGjsudZs/2bNbFll0IPWBlqrFJapJQueFahenY3qZzUgXDQrJcX09fSqVXj8oouYie1JU1Mx0xdEjZvi6Zu10rELJ2KpM61MZFm1oqlA9+w8jfigfT4fqvz+CGuPiJUPK1wvRxhfFkYrPNkgty+J1jNYmdGvQE+6wHvu4QqBbKZ8Vq/Xzvp5v49oZvcAIvZrilaujNj8Dnq9mHDBBeYqapGEFPqky7ZdZopuMH2MhngJWy1OtZ/Ru45GxWVUVtaAkXLqlMoqp9zjQetNN4VnbsoHvwLGUR/ljIzw/yuIZ5DHzWKU/MXskp8WXVSJAWR2c9iOWDM8yW1IU0lF8Ot9Hyr1YJTlUownyEEGOP/en9++PXwur6mpFRJOveOESsNNahIj3KjScqr9vHV1yCouVtmSy2lpaKmttXUJT9uIpAmAUL9+QChEfdagefNMOf+QqgrS+YdUgUTzrsnY/1T1kccDhEKmVHCsshpBc1IyUgMa1dmofiQ8qiyWKo5ZJ9A3aUm1HNBtmmnGzFaod3RwwkzRLaaPvRWn2m/A8uURzkNSZycGEI5aPEG8jCJ/ajciWQkzpLY2Zj2XM66xGmbYbqRg0LRPR3n//twRJ5XNSyVIWqPfH55hR+tXQr7HzJoatJaWomXdOsNZvyXfAYtoBX4sNQkJJ/Sd0OE6lWMzETFr1WKFFMqmKe04bQlvNvInCUvoVzAWzZ5jx7DMpC/A0lGjqIJ06ahRhuUD6O/BrKqCd2Cm6aq1jlrK/yUAoawsqiolmskB6z0CMIxGyjOAss6xUz0SayfKhBP6TsyQ3GL66HaYceYZiaPj2X5rlywx5fLPAyuDkqKmMGO2On/DBjSvW4fOnsGisyeA2PwNG1TnGQp34j0MWLbMdOIPnoFZq6umhSQgB8pKIn8tz7P0yqA3UCiDPW3Wz7uHwHpvct++htdqYc3iY61JSDihb6dNuJP3jDdOLCd5rVoAe9qPNZgA7BAICizTOQSD1FlcBcCtwlCVsaee0VhSKUILAHVmnPLMM9zCPSUQUCULJ8vNshcH+AdmMgaQXkgCQH9wNFsGnoHCihUb69q2X/6SGk9KTqGLVMU3QbFsGurzOb4SZpFwQt8JM0U3mT7agVPLST2rFjN6dVp5aeecrqpCKC0t4nzFQkP7cfFVwkNVfVRALRxlAKEe6xvaYFCOSKscRTAuKCoKC3EeBxyWRYdn+XJD4c6CnHkrZaYlUdEOzEbvq6SkJFzP8vR0atuwMDO5Csc4YqjUollFsupGUw+y4knJAwYgxBD84fM0AddivRJOOOsdEjdasjgNT52dcFiK5r48ljVG5/Baa+iFPibjpzPD+KLbXFIKBJAiywj1ON2kb9tGjauuuPLTEpMMmjfP0BmHB7MhpIPZ2ZDa2/XDAetY70RjCaVYV7HejTbsNI/1ztknn9R1kjLThwCE+w8tNAR5H9Lihmz7CpxfOcmShLM33ng+LwGMw3kHs7IgdXRwh58Q8fR1EEKfjlMxecwKBZ5BgncgMettykqywfIiJT/MCpwfJM59+9vo+49/MGfLWgEc8nrhCQQsCX1FMFb5/dRNZVK4k881GvSU8kKSqAI3msmCUtY0vx9dPSavqS0t4f+nnDplyjw0NzcXKV//Oj1JDugxc1jlprWTXt3Id5U3ciQ8PSsq8r0rKxUzeydkxjUeU1WrQj+V56T9+/fjqaeeQigUwqRJkzBt2rSIcw4ePIhNmzYhGAwiMzMTK1awbBwE8SaYn0//CEx0Jr0ZGW+4BWaOXOI4t75z2DDg8GHdMmuvCUyfrpoNKjM+mTLjgyRFRGRMCQTQd/fusFCnheclg5KtAICee5Cp/ZT/a2e9NMjZcyUiLYlCaWk4XdmtNde+h/nTp2M+8UzaECkB3YKZsIJR2siK7nk50O3f0NOuZLgK2rP0YD2vAsDtlMGHqXb88kvDWTizbgzPaKmjA1IoZHBXNcH8fFVfdBpDnX4oFMLGjRuxZMkSrF27Fjt37sRRzWzw7Nmz2LBhAxYtWoTa2tqYeJUJ6FRVVRmeY3VjWm9PwJQHMMtShjhutLkX1u0ePky1GqFdA6j15FpTRpq3LFNPTgTRqtB5dgXUewI8yVZo0HwTSJQW4HkPRnbsWgsSlu6ZdZxsV8ULGaAPjmasVVh9wuwGMQ/l/furksUo//9Bz4ARsbFPCPwK4j7hQHYG4TdcEU+/oaEBQ4YMQV5eHlJTUzFhwgTs27dPdc4//vEPfPe730Vubi4AYODAgc6UVmDIfffdZ3iO1Y1p20zMWIk7iON6AxRNqCjC2ujjIgUsyyGLjLZoNhuXEWY3rRUqWcJG+b2zk+s9FBcXqxO5MM5TzXRZ6jPGcVZUz4rIUwEAVZxqLlafMHu+bDDohbxeFK1cSR2gr05PNyynahXm8aCRTLbC+O5iEU8fsgG7d++W169fH/59+/bt8oYNG1TnPPXUU/Lvfvc7uby8XL733nvlN954g3qv1157TV60aJG8aNEiWZZluaOjw5Gfzk2b5FBhoRySJDlUWCh3btrk2LPc9gPA8WeEJEmWuz911U9Ikszdx+Oh3wdQvTfW+wwVFoavAXl9zznaa8rKyqKuT+emTXIoI0P994wMuWvOHDnUpw/1eqUu2mPlxPXLfvYzw2eQ/Ze8H1jP5HgPtGdR71VYGPV7Z53P+uHpu8FgUNUnlhOTaPJH+65p/SGUk2Oq/4HsW4MG0fse4/2EDOpGewat3h0dHUYi2xDDjdzdu3fjnXfewR133AEA2LFjBxoaGjBr1qzwORs3bsSnn36KZcuW4dy5cygrK0NpaanhhkNvib3jdtasWUOdIfDoiKOB3MQiIZNl88AKH6xg9N544tOQsDZNeTcoafsYADCwuFiVPFwpR9DnQ9fw4ei7ezfQkyVL7xlmQy6zLHDsiinDGyqZ9TxWP6Gh3eDWwtpDoiUo5/3eeQwayPtXgG7lRb6HctDVTMslCbczjCQemj2busF+7zXXhJ3xYhp7JycnB01NTeHfm5qakJ2dHXHOZZddhvT0dAwYMAAXX3wxPv/8c0sFi5ZkjJPDSg7thMAHYFvYX1aaQgWj98aK+VLev7+pchipkLRxXUg9eWZNjUrgK+VQBGWft94Kp0WkURXNpnUPNN04794M656k7llrT94+aZKpvaAV7e3MZwSzsxHMykI5utvLQ2xwaxOI6O0hWfneeZzBaCkVtZ7b5HuogHrvRvk/K0QHAFS/+y71mup33zWsQzQYCv0RI0agsbERJ06cQFdXF3bt2oUxY8aozhkzZgw++OADBINBdHR0oKGhQWWdEEtEnBznYW1q8jgFkfCEwtV7b6z45EUrV4bPWbNmDXUjjhQqrD0OAIZObHr9jSdTFTkr5PJIJRx/KohzZNBDNLPQyxrFSuBtNrNVFUPoA0Dg2muR0traPXPuCUEN0CcseoKdbH9S+PJ87zwGDTzJ2iuIa0OaCXH4dB156DqPXI/Hg1mzZqG6uhpFRUUYP348hg0bhvr6etTX1wMACgoKcPnll2PBggVYsmQJfvCDH6CwsNCRAhuR7HFyysrKHH+GXW2sTcRh9p6B6dPDm5FAjzpCkpB1zz3hTVDeVRDN2oUn/Z5eWxh9tFqBwyOEWDFfZK/XMMUfuTms9yw9IURrJ63FSdiiioUkdef/JRLD9/vjH5mn65WHbP8K4m88fZHHoIF5f0ay9tOVlZB7vMSVQUhOSYHU1sbcnCefUc44bicJ55yVjDp9Ersd0ph6bBfkLKDpW0m015t1hOJJv6eXzJulA2fpfYuLi7Hi5Mnz3pxEqj2jMrH2Mbi8UTV6crNOWGS7svTT5VA7gNHUXeUMvbdeeVpLSx393s3uGXjr6qh7PFoHPbNpG0U8fR0SLU5OPOEKWxvH+EY0fSuJVrdr1n+EJ/2ensqDNZv+zbp11NVH2fDhyHj+edUMWEm1Z1Qm1nGeKJRaW34rfhws/XSF4ZVsvbdeeZz+3lnZrlj3Z+3xkGj7ZaxlVsLN9EncEobBbJYcK9hZZ6di9NgFT+gFs6El1qxZg7Lhw5nxWGj313OhN4olQ86Sedrb7Ioo2pAbRuVmWYwtB5EyEZHWRSGvF1J7O71MHg8aGV7VRuVxCrOWcTx9EtBv/8OHPXjzzT74yU/a4fV230vE3uEkXkJfK+StBNUyi511dipGj13wmB2aHaB8Ph+CXi81XSI1iBlHsCwFmuAig6/xtveAxYt1VUAksRi4aW0GaAKR4bxFU9q+fd06feJcGcDZm25i1iNeQp+E5zvmTa/Y5fPhzWffwv/8jxevvJKO997rE3FOaelpzJt3BoBQ77iemHjVxQCmekOWXZEcnifRd/ukSdS/6bm707xzZcJKiLy/dO6c4WYvcD5mDqkqy9Kom3hUN966OkMVEEmsckHoed4qAl1RIZ1etQpnb7rpvDewx2Mo8GOZWcoKtPb+EN/ASizGaLwNCTIkyEjzH8UVV+Rh1aoBEQJ/0KAgbr75LG6//YwjZUxaoe90TkoeU0G3wxKqpO22XruRbZw3ciTyLr3U1vY2sv6RAKRv20a9lhyYte+KliCFFnJZQndeXBraNiJj5lT0/F3R/SrPXjpqlKGA5rEoIomFvpiV/1dBAtDvj39U50JYtQqNhw+j0e9H4+HDTIEPuMf3xmhPqKEhFSs/vwn/ln04LNwlyLgIH2EpVmI/RqvO/8pXgrj11jPYsuUkDh8+Br+/++e9945j5cpT6BM5+beFpFTvOGHho6f7q62t7ZXqHUCtSmDFmaeFs6W1MQnNgiSUlYVowu0CwBM+HzXMsFnvXFb0SZbFCes4SZfPx4z3LwGqchipMXgsilj92I69JVr5WFZKZsvHwm1qxo8+SsUrr6TjlVe8+PDDyCQ+WvLygvjJTwL4yU/aMWbMORjkWKEidPqcsASgU3pOMnZ4JyEIe6tO31tXh6z587lCxWo/Zh7dJk0frndPVhmVNjcKScCzKUfTT4e8Xsjp6apwwOE6cMRkV8rFI/SN4AltwGNayQs5ULAmSzSzVS0VOL9y4glvYTWOv1VkGfjgg1S88kq3zr2hwVi4Dx2qCPcAvvWtzqiEOwuh07eIEx5wpN4RUC/tnQo17bSKauCiRdyxwbXLbZ62TGlp0RUURkt4bZsrVPT8S4usaeSoVVxcTFeHXHttWHgr95fR7Vlq5GAGsMMPLzObYJsj1AVvtEoeyEGSpWZRma0ChonRyb5BvsOJDH29nfsSWvWqLAMHDqSipiYTV1wxGD5fPny+fBQU5GPy5MH47W8zIwR+QUEX7rjjDF555V84evQYAAl+/zH885/HUVFxGmPG2Cvw7UbM9AmszBzIe1aAPauxC9asK7R+PU5MmWLLM4wComkhl9s8M/1yRCYB0bunlscvuoiZ4DxkYNXCO+ulpWMkTRFZwdBUdQDQsm5dhNNOKC0Np2prdZ3PIlRfLS2G74RcPVgNxqfKMMZQUWnfkZk2I/sJec7yzEzc8cEHkffsWQ0oaSir/H4so6gXI8ooAwcOpOGHP9yAr31tIT77zDh/VGFhV1gtM2pUp2q8jXWQQ6He4SSWOv1Y6x1ZQlUuLETj7t22PENPiBlFjDTS6QPsKJGse0aUjxFlk8dzkke/zaoDK6KnXns1+v2mzA552o+kAucjQLKiVRoNdGQaRprHMOlVS8J6RzzRI/XeIa2s2hy52nfddt101NenY82aTBw8aKySAYCvfa1buF97bQAjR3aZjRuoalenfHKEesciZMwWxWSsbeZMSxYNsY75w1SfHDli30N0omkaLbe1ViPBrCwEs7MNhbzePRXCsV00g6w24qYnEMDaJUuo99Z+mDRVmdbjlxXRc8DixRjaE2uKZtbZccUVANTxfZQNUJZqjidYG/mMFUSdAfPWYtqENErwOuX+LK9aWZJQ5fdT6zB/wwY0r1uHzh6rqM6eVYMi8AHg+2lp1HbVI7OmBgi0ow7XAXgHEmR4Am0YNO9uFBTkY9asQQyB/yGA+wBcBkBCcXEJ/P5j2LnzBBYvbsWll5oX+Fp6g7l2Ugp9s7bOPMTKHlqBOZgMG2bbM87eeCNViJ298UYuM0BSyB0/eBDHDxxgC090B6YKeb2699SmNlQoBz2k7YozxrbONDvwrHvuUeXqraDcP+j1Yslll6mDhxHnKQK/+dlnDZ/HG72TJBwCWVMmgL1fwYKlrzd8tiyjEmwTXqUPAKAGhHtpzRpVpFSS7oGrAHPmbMP3v/8V+Hz56Nu3D9L8R+FBCDNQB2AUtXzjxnXg2WdPhnXu3T8XAVgGwN6Qxax2nTFjBvc9nN6fI0lo9c769esxd+7ciONOW+9o9Y5OeBHGQqcPmPP+5IGlw1UwUrPx7BWo1C+gm5Sy7lkB4zgxEhC2zsr6zW8gUdI+KiEFzJg48iRUIfkqAHrQAnMWQQBbPcnad5E9nnC9Wfp6Xr23t64OGaseQPqxK/ANz334OPh1EyX/G97AClyNHcwwDqT6xQlLOlY9Ab73EOuAawkt9FkvOBb692ijRppxNaed32/OnKhMNmPl5j5g8eKwC340mZ94YptUwDjiJuueRvsMIa8Xi7///fM6aaPNW0ofoDl6AeyMTaxyeAKBiH0MoFtQ325SsLEGGVaYiRWBAD1DFOPZyrfY1QW8+KIXtbWZOHzYeDNV4aqr2lFUdAY/+tEApPXtS90gVvZOAOuC2AiW7j6aAYZnEip0+haxqn/nWYqZ9SKMxtWcFSXRLFbd3M0sTdO3bQt/sOWMc/RUGyzTxwri/7T76rU9z3snk5SQOmlmYnSPh60yYVxDliNiT6Qn05Q2yQuNZVEkMGKpJ09XVVFVecsU88yec5X/K8/u7ASeeSYDY8d2m0ECMny+fHz1q/n4zW+ymQL/Bz9ox5/+9K+wvl35eeaZZowde667nYj6qeLPE8dZ5rlWzaeVvl5bW8vs6zRPfD1VT6yTqCTcTJ9nSWnFekeJoSKR8bLT0tCiMb0zu5qwS+UUjXOWlWfbFfWR97ks5yTtjI/HtFCBKwY8Q3VArlzCz0F3rJl+Tz9N7wPoTnoSTf+bMWMG9uzZE3H8qwA+I54djQouGuuirkAX+uIW5KMcx8A/2EyZ0o6iolZcdlmn8ckEWusdBb32s0ulwxNbn1wBiJl+jOBxwLESj4SMoaIgdXZiwHJ1BBLWjJR13Opor8xA0tLTTW8EeRgdk3WcxOyKxmhWbbT5nXLqlOH1rDZmPVtrZUK1VqHo7QGogocpm9Gy14t+Tz8NlocOK+sST//bsmULtayf9fxdL96QEXorx44O4MYbd2PUqDz4fPkYNO9ueAJt6ItzAJ7UEfh1AEbj/GZqt9XMpk3NpgU+WU4z1nfk7N5K7Ku1S5bAo0mI7gkEMI2QLdGYa8baCCThhD4vrA5upKpIobjiU48zZrMpX34ZEX2RZn6owKN6sByFUEdFYXgpMShVMI6T0Dq4LEncOV616etoZppLxo83/RGRVibU5+qoTJTgYZXoTmeYEgh0JxUPBiO8eMnkHyf27sWCoiJTqjmWFzKJFbVAezvw+9/3w8iRQ8LeqT5fPr7+9Xz87W8z0NRE7xM//nEAr712Ap2+AiLUmAQZMyBjPwB2qkqz0Kzv+j39NIb6fNRvlnyeFZPKFWfOUFVaO86do54/btw4rvvGOokK/05KL8RsvlitqoLMFGX2BbBmpGT0xbR9+3RjlvCO9kbZkQxhzGKZx8lT8vPDAmgFzgs31mCllCfaTWNtejxlia11SDoV5cb0vddcg9Abb0SoDnhnXSxHrhUAyokZqaJKqfX7cf+zz3KXj3zXVzPO4coPG5Dwxz9mYO3aTLS28s/9pk1rw/z5Z/DNb3apjnvr6pB5Sw3X6tAqtP6uqNGsfLNGkH2dhy1btnCfG5g+PWb5ARJ6pr9s2TJT5/OoKnjVNjyp9vpt3kwV+LyzXgU91RDXJitjRl/B4akSzdLUygY0OStaATA3NaN9xvwNG0zNurjCMiv/9viDDFi8mBmnyQjyXb9B+bu27VetWof16/vhG99Qz9wvuGAoKisHMgX+9Olt2L79BIqLS0CqZl56qR8mTsxTqUlYvhMkyzIzVb9bsUs3Wslov1m7wpxr+zq5yuxNodMTbiOXxOymJs/mKy3xMS2GCo8bvdkNRxZ65nYpZ88abzozzA55IkDGOgYJCblRxvLJsAse93ptWOYK0G3crwKwg3J8mSbeDA3Wuz6N/nhEugcrpAqcC/GFHwCAX/yiDffc04rhw41XdayNSSOfgmgSgbPIzc1Fyte/bjjjJs03eerA836dNgPVQ2zkOgSPKWdg+nScqq1VzwQpQbOMEnwA4DLd44E125Y6O7k2nfX01UbwbJzbCWvWpne+HfDqgsn3UAG6l/AbjOM83sP++ctRmboCqehUac4HohVL5WqmwL/++rPYteu4ygzS7z+GtWtbuAS+HqyZN2vFylpRs0JmaDHKmNZdKI+p1QTP+y0bPlzlPUx6QFsh1qsDLqG/f/9+zJ8/H/PmzcNLL70U8feDBw/i5ptvxsKFC7Fw4UK88MILdpczJrSWliKUpv5oQmlpEaoKXtWBcl7pNddQhfLZG26wZdeetREknT1LPV+76dw+aZJqw9HNy1bWIMNS5el9zHa4vmsHIcW642p0W5VQ4RjsW1okPPhgpkol4/Pl42v3zkF513IEGdtxd2A95vUbCVIlA0h45pn+eOEFtm0/Dywbd+ZkqcfkUPt9sAYJVsRULVwTqmCQatxw7zXXcD2DBmuw0qqueCC/p1jH6zEU+qFQCBs3bsSSJUuwdu1a7Ny5E0cpqoeLL74Yq1evxurVq/Hzn//ckcLGAkXYVWh+t8IDf/kLVSifXrXKtl17ZYDpbG83rScnHaYqQI8xUzZ8uOF9nMobwIsZIc5j8cSjC2YNQjsAtPz2t+FBXXEi0g72CzEIZagK501VhPvIkUOxdi1bmNxyy1l8jkKNpYyE9bgTD7X9nyOrL9b1Zvd17AhCqPR35irV46EK6Op3z8fdMavrJwcr0imMZ4WmJZ6B2QyFfkNDA4YMGYK8vDykpqZiwoQJ2LdvXyzKFnMya2rC6hBFFyt1dtqSi5O1OrDTq1Zrp8+76Wx2Y4wFj1Cxc8VADjIpzzxzPkqkLCPN78egefMA0D9mnk17q6orckZ6J76Ce/s/Ak+gDZl//AM8gTZIkLEaTagG28rsttvO4K23vohQy1RXn8KwQvqUxKnIrizM2s2Tg0QFrK0sWQMOy/KM7Otm36/WXJh23Axxy6EtG7B79255/fr14d+3b98ub9iwQXXOgQMH5FtvvVVesGCBXF1dLR8+fJh6r9dee01etGiRvGjRIlmWZbmjo8PRn2AwaOr8kCTJcrfFvIyef2VADkmS6WeXlZWRwRbDP2VlZVHXh3Vt56ZNcigjI1xeGZBDGRly15w5cigtTX08LU3u3LRJXe/CQtU5yk+5xTag/QBw5F2z6sB6HvmuZY568pT7yiuvlIHBMvAg7da6P7/5TZf8+ecm+/fTT1PfO/l+rfQ33h9W/9P2s4hrCgvlkCTJocJCuXPTJq42pn3TtHux+kOosDDq9xtNPSP7R6RMMFtvqxha7+zevRvvvPMO7rjjDgDAjh070NDQgFmzZoXPaWtrQ0pKCtLT0/H2229j06ZNePjhhw0HHLdZ77AyMfFYVeihtRiINriZWeuJrp5IkEbP4rE0sisDmFP5glmWVywLJLOhJ7Tl/uKLFDz2WH9s3NjfVDnvuqsVt99+Fjk5fGko9VBCEsQiUJ4edoUQ4ekbvN+0WQsh05nUGO1NWgFpLYJYwdh4krDE1HonJycHTU1N4d+bmpqQnZ2tOicjIwPp6ekAgG9961sIBoM4ffq0pYLFg6KVK6m780UrV9r2DMvesxT07PR51EfajTFtjlOrLuHR2Emb3mRl5BFg2Yd7/H6ueh47loKlSwdCCRim/Hz720MMBP59eO+9xgi1zJIlrbYIfAUr6kG7YrjbFTDMzj0hp7xcjdqb1NXz6u3JesdC128o9EeMGIHGxkacOHECXV1d2LVrF8aMGaM6p6WlBcqCoaGhAaFQCJlR7GjHG7KjAGwHKfJjyRs5EnmXXqr74ZAv1WysGh6BaUfWrnBn9vvR8vDDtn4svLrTcLv6fMi65x5TA2OwspKq2yUHbK0TkSTL4RAQn+Z9B7d9ay8GzbtbJdy/850h2LSpH/O5RUWt+PT+36HTV4CQlIJOXwGa1z2C4uIvMWhQ3FxgDLFz8sHT/3h01Xab+RoJaLuctvRgPYMM0eC0P4sWLuest99+G3/4wx8QCoUwceJETJ8+HfX19QCAqVOn4tVXX0V9fT08Hg/69OmDm266CRdeeKHhw92i3qEt2e47dIj6MoxUIbIkAbLMTNyhF389lJ0NyDJSTp2iLh1ZS1DepWysYubrYaYOWvTUBTyqjsFjx+KIPw3VWIqNmM1dZo9HRjC4DM3pG5DdflwVYbFt5sxwKI0KnM9T62TsFIVoIqoqMHMsezxAKGSqf/D0P7vUelbqrIeV8uk5bSkUFxejtraW6Rhm5OAokqhwwtNBzOr+eDIa6d1naGEhNdMSz/V6HdNIoDuRDD4a1qxZg7Lhw7mzSZHoeSpr3/WhQx48/HAmnnsug7tsffrIKC5uxS23nEVmpvqzIL1tVWGcOTJIOYUVAcgT4tpM/zDqf3YlD3ej0GfdR3tPnmewzhEeuTZiVt1iRk9JvQ+nwKddr6fzNLLTN1tPpygbPpyqVuAJ1BUaODDiWEODB/PmZaFv3z4qtcwVV+QxBX4/nMEDWIhW9IcMKaySOfuVYVh5fxZGTPpOhJrjKtBz+1ZwmAZahaVuqKqqMjyHBY/qz0z/oKlSWKoNNyYPj4WPSbz9WBQSRujTNqXIj4KF2dDALFtdnvsD5kMeaO2KoyXW2XlYmM0mpfA+LsYNZzdEeKhefXUe6urown3AgBCWLz+FhobuDdXmdY8g6M3AGWRiIR5Ef5xFyOtF+6RJhvrtNxQb9J7fyf9TE733V2/yWtETs4TkfffdZ3gOC65QBrDWP1h7OW7E7LfFep96gp3nGbEYGBJC6LM2pciPggUpxFcwjpOQHwstmJbe/bXX82CXo40dm73RQn4gTCESDCLk9eI9XIqZeE7jZypjJN7Hs530lHNZWSHU1nahoeG8lUzzukfQlFmIiqpsfPXq7pk7y6Ijfds241UQY0ZfAVAtvjquuEJ1nttmt9q2YIWMsLN/xM0ZyQH03qfb65kQQp85e+TArAu51sInmJWFYHY2t6mjkXkkbznMEuvsPCTkB6IIkf24DNehLizUUyDDE2jDKLyHFzCTep/BKf9CdXULPvlEbQZ58OAXuOuuEJTq6Vmm0NQQPKsgVl7W8sxMajYlWvpFM7BUIzNmzGAGnDMrYMhELmTICAU7+0dxcbHjgfncIFitenG7wmSzN6BV0Zhx6/7R5s3Uj/a+Q4eo569ZswaD5s1DmhILvaUFqV9+iQXFxdymjuGPTXMNLfm1XZussc7Oo/Duu2kA/hRWx6T5j0KCjNHYj5dwHfWaoQPPYH3a3ehAn/BcP+jNwAcPPYtbbmlDj0sIE7tSOJLHtaEDgPMmoSwVhhWTQJbw2LJlC/N50QrS2p4osU72D1aZ7BTUTgtMvYGY93o3kBDWOyyLGp548Hq77UZY3fF3wjvVKesGI95+Ow21tZl4/XUDiUxQiM+xHJW4CX9EGrozMSkWI4C57Fpkvc0mpbfL5JUV4lnPXI8Hn8+H5nXrIp49aN68qPuu9v529UM9j1TaOWbj2+v1b6e8vY2exbJM0vPIJYm1yWZCzPSZQZcEtrNvXx9cf/2giA3Va6/9io7A/xS//e2X+PxztVrmE9/3cBt+Hxb4gDrNo5GnKWvmZHb/gjdgmFFe5eUAOnsENNA94Whetw73P/ssAETt9XrvNddQ1VVXEb4wZjcAnXJM4vFI5VmJ8M7aY+FgZQay3LwrrljnpEgIoc9amrJy5PJ4yelBfuRmP2SnOiktyqYV9uzpg5//PCdCuE+blosdO+jCffjwLjzyyJc4fFgt3IERmDkzgNRUvk1dXosRlmAwu39BS7Sd8fzzXG340OzZ56N7Qp2jldxbKIe5tIgk1e++S1VX7fjww/DvZgVErAUNiZ3fQNnw4Qimpak30tPSuEKBm4Us37hx45h7K9prjOoa6wEqIdQ7LHhUHWaXhHY6Odm1HOUpE225LMvAzp19UFubiTff7Mv9vG98oxNFRa34yU/ajSwtI55N1tlqoC7yXtp3rahiJvn92MbwjlbgKQdL3UA6bZEsz8zE8gEDoq4f+TyzweTMOj9Z7Ye8Hqms9iNVI2bUHACQN3IkPC0tANTOccGsLBw/eDB8XzsGM1Y76an1eNQ7rgu41pvhsdPXosyYn/D5qDPmtUuWuMLJicRo41KWgdraA/jpT3NVs/aCgnz88pe5TIGfm3sMTz7ZjCNH1DP3N974F372Mz6BD9iXfANgz5y071pRxewADAOR8aw4zG4Srjhzhp0hikPAks8jbf4rYGyowFNWO2fbevb4ZlYQ0dj1V7a0UP0kKnsGAsCZDV5t+5Eo5d61a5fp+8ZipZXQQp/HTp/UhZLL8UrQl+OslG7ROLHY5YihPFsG8Gdcg7F4U5WJqaAgH8Bf8dZbfajXX3rpOfz+9804evQYACks3E+e9OHHP25HisVewhIw9x06ZNpixGy6RG05SJQBHozFrtbBilWfsHeu5lrWHkKlYUnVVLW2qiyHyNysRkLVLRYjJGSZeL4BPaFdAXq+4Yroi6dCLyez0QC1Z88e1e9kXeO5F5HQ6h2zy1Zymc+Ko8JazvPGrrcDWQZee60v1qzJxIEDdEFOZy+AFSgquggLFtCFhBVrJh7svCdLvTNjxoyID05BOd8owFvI64WH8Tdy2e7z+RD0eiNUa55AAM3r1lHVbp5AwJQVBwBD6x29a6ON9xINpHpixowZ2LJlS1TPY6kEI9Q7l14KT0/OZ/KbXZ6ejqr29oj7stRMPPCoZ8g687Yrz3lCvaODlRG0yu+nLhUn+f26MzteV36zhELAn/+cjsmTvxKhlrn11hymwB+PXfhzn5+i6eFHiM1UCX5/Afz+30UIfL3ZjJMzEKv3JGdOpHpHz5ZdgaYSA6Cy3tFeS/6fbCPFz6Mc58NxA2oDg3Kc9wEhr9fLt0ui+IYsKCoKq6tIQwW3hDwgBSpr4OWFR413urIScloagPNOc3JaGopWr47ZRjXZD/fs2UMtt9aWX8z0CezadAGcm+nTZl2sKJE8m3ahEPDKK+lYuzYTH32Uxl3esWM7UFzciiuuOIeMF/mjHOoRy5m+U7N+7XEay0EPo1HBOA7QZ8xkf6jqUQtq4d3Q097XqJ1Yhgq8G4sKdn5z2nKY2aTluQ+tzkah0e3qZzztFE3/jvVMP9XS1Q5QW1sb86QCCq2lpdTlOC2UglZtk3XPPdR7qlz5g8Cf/uTFmjWZOHSIv+nHj+9ASUkrxo8/xzxHKRNLEOjpTrUfJEtokOfH6x2ZZdy4cdTl9uMXXQRQ9mcqcH4FR1rH6M3AlD6zAt0Dhp4Khxcr+z1a9ZNROewOhUDrS0oYBt4ymYH2Pdb6fOF62bV3prdvwvp+yPq75ZtJOPUOEL0dPbkcXw565ixWByI37brgwdP4NUagASlyKKyWKSzMx913ZzMF/pVXdqCu7iS0KfZeeKFJV+Cz6k9GHNXrcCzVAKuuZq0heNRHvK7sPPclBTRNpwyoN0d5YLXfssxM03GfeIQQ+Txe/xG3YLdaxWjCYkQshK0V66VYh1x2hXrHytJPSyyThXR2Alu2eFFbmwm/n3/mPnFiO37zm1b87GeDbVefDH7tNaTMnRt1/c2qHszCWv5Gc0/yXZPZrHiSX6vUdD17OQBbvcPqiywb+nIAtzusutJT7xjZfTsFzzu1M4mKXSokK9jZp1kk3Eaund6BVpKFsGYN584Bmzdn4NvfzlNtqH7ta/koKclmCvypIz/DK6/8K2Lm/v/+XzPGjOk0XTcePMuXW/IjYM063ObuDkT3rsk6KJujpddcE7VJJMssc5mBiixWxDOKJasvOSWA4+llrEDW2S1JU7S4QujbiRXX/traR/CHP2TgssvUwn348Hzce28WvviC7o307/8ewKuvRgr3p+r7YPRotXB3XHgeOUI9zOtHwPpA7PqgyA+B5crO2xZkncoZx/WcaBS6Lr6Y6i/AAyv6Js3BjKZ2Y2G2n8RiUGap9fTUd05gNdqlk2Ug0X4b8Qy9QOIK9Q6J1eUoj0t9ezuweXM/1NZmoqWFf9z78Y8DmD+/FSNHdhmfzIET1jFDx4+HdPhwxHE7c7a6xX6ffNffB/BGz3FWXc0uvXn7omI9kub3o9NHD/lgRe1oRb1jJ06o/sx+73rqnXiptWJhGZdw6h0Sqy+qtbQ0bLcbQDrW4jfIwpdh71SfLx8jRuRj+fKBOgL/v3DTTasjZu5PPvmlbQLfKYKVlczQBmZmmnq4ZdlKzrK39xyzmvhDazfPgxLyAWCHfHBLjmKz8AYMi3aV4VQMfDdYyrjRGxpwodCPhvZ2oL6+L+bNy8LQojuQ0nkOEmRkIIBirMUpZFGvmz69Ddu3n9A4MB2D338VVq26wfFy8wpPM8I6dP31TFWFXc5jbviggMgsZkZhHHj0rXZZJmk/eCtqR7ODrJ2Dshv05Fri5UzIQq8vua2sAKd6Z//+/XjqqacQCoUwadIkTJs2jXpeQ0MDli5diqKiIi4zM7NhGAIBYPv2dLzySjpeecWLzk52qkGSm/AHlOZtwMC3n9U9z6llsRXMqgVYS36rES2dJpq2Z4VbIG3zY1EOnmudaP9YJ8yxS71jxcqGV73jBrTl6TXqnVAohI0bN2LJkiVYu3Ytdu7ciaOU7EOhUAibN2/G5ZdfbqlAABAISPif/0nH3LnZGDZsaFgtc8EF+bjttkF48cWMCIF/1TeP4JGsUjQiT5NSW8IfcAsuOrHT8LluUVuQ2KUWYEV21JtpOj0DiWazkoQVbmHLli2mVkesvLN21t/OHMV2507gxU4nJztXD26YPSu40cpNi6FxeUNDA4YMGYK8vDwAwIQJE7Bv3z4UFBSozvvLX/6C7373u/jkk0+iLkxDQyquvnqw7jlXXtmBn/wkgB/+sB25uSEAxsGzALZpHYlb1BYkVhONKFSCbn8eGjiQeY3T3tElJSWmPEe15WF5Qt57zTVY9cYb4f5AJjahrY727NkTDpqmsvcfPhz6blZq9ISi8lwlXMMyxoYvDXJTUtvXjepmJzx9IRYTJyvB5ewsA6s9eGbz8ZxgGs70m5ubkZOTE/49JycHzc3NEefs3bsXU6dOtVSYPn3Oa5q+//12PPhgC9577wvVZup//VcTfv3rtrDAB9jBsxSsbu7FE7Op/0wj8anIYgnvrIg1YyQzTVX0nGu0OrJjNWUkFJUN30oYx/gnIQVcLDaErcxKrZjwmnlGvIPLsQad2tra8yG7oY4IQLZNXCeYsgG7du2S169fH/59+/bt8saNG1XnrFmzRv7www9lWZblRx55RN69ezf1Xq+99pq8aNEiedGiRbIsy3JHRwf3T1lZGfNvIUmS5e6Iw6qfECCHCgvlzk2bTD0rXj+dmzbJocJCOSRJ4XJ3btokhzIy1PXKyGDWKRgMqtoM50OMh3/KyXtJUkQ7067Ra387fsj7AzBdHvIasj+AUVfWfVnn2/lDltXs+cy+bmNZzZaP9XPllVfa3nZk/6b1FSf7KE8boef7JPue3vcaTb2tYriR+9FHH+H555/H0qVLAQAvvvgiAOC6664Ln3PXXXeF/3/69Gn07dsXc+bMwdie0Y6FmY3caDbJ5MJCNO7ezf2MeKK3YQvwqwXWr1+PuXPnRhz3+dh5AFgbifHaGGM9V688ZL15oqVq76ucw3N+NJjdvGSdvywzE5WUQHFWyqpVVdj13p3oP9qN3Fja47Peybhx4wzDSC/LzMQdH3wQ9bNjupE7YsQINDY24sSJE+jq6sKuXbswZswY1TmPPvpo+GfcuHGYPXu2ocC3E9YmWfno0TErQzSQy2i9ZbsZtYBetjC7NhJ54N2MJY9b2QAjM2ctHTWKmhdh6ahRzOtZ59vl22B285J1ftHKlba/Rz3TQjdtQNKIpZqE9U60BgW0TF4rzpyJWTmNMBT6Ho8Hs2bNQnV1NYqKijB+/HgMGzYM9fX1qK+vd7RwvB2RtNcmbdOrXn7Z0fJZhZw12LVhy6K4uNh0akIr0R31dJ40eISinv6XTK4xf8MGNK9bh86eftPp6w6uNn/DBuq1xcXF1PNLf/Qj2xPjWIXV161u4tplUcOyhIpliAQ3YtsenA24LgwDi2iWim6z29VClo+lolqemUnNy0uqBZyKLuiEzboTrvxWn8c636xtPSs9oBazKgmzYRh40OsztbW1vUa9Ey9Y7/Ch2bPDlmM8kV95SegwDFaJ9TLV7PKfVb6lo0ZRl+1FK1dyzcDc4DWpFwjLzDux05zNijer2dUXb3pAnvcST3WLG/1V3AbrHc7fsMGUl3g86DUz/Wg2bJye6UcTREtvhlU2fHhUKQ/t9Eq0a9Vg18xbD56yWtnoMzvTj3WgM7tmvU6Vm2flw3o/rONumemz0OuTQPR7EHbO9HuN0I8GpzpzOCen/3zyDRJeSwqz5WN9COR9WNY70RAv9Y6dA7yVOvAM6k6EgtDiRqFvl9WM2ffmdqFPoheGgQeyjYV6h5OysjLb76kIglSGwAeMN19ZzhtGlJSUhK99ome5TwvkpBBPywuWioBHdeBU5EWz3HfoEHXT9L5Dh8LnTJgwwdIzeN6R0+qWaO7vlneUyDjVxgkh9FkfDmnGZxcs798K4v96O/UPzZ4dztgEAGl+PwbNm4eHZs82fDY54KwAPcMTWWerncaKsDE7C7Q6QJEDvBV9OHlObW1t2Fy28ejRsLks2a4sqxde/T7PO3J6b6a37v24nXHjxrnSFDYhhH4sZx2sWbwS18bIZrr63XepdrzV775r+GyWLb9TOCEM9Ew57bLTt2J+aGfSd+15vYV4ZeeaMGGCKwwS7GLPnj2m+yKrjUmTZKskhNCPJXqzeNryX4sVe3zWOcuJ/1dVVRl+mG4VQG754M3MQrXlpt1Hu7nHKzztcgwzixVfCh7cGKPfLbDaxk6tRa8V+vFaNmlzopKenIqqRm+2SA4aZF7X8v79DZ/Nk4R72bJlhh9UrPWxZpfzdsKjotLLo8vyurQiqHiEnkqV5xLHsFjSW81G9eSSW+rUa4V+vGYLpEdkuSSFvTfJMujBSqRNc8DSu1ahN0QQ5XVfV/7P+3HwDPA8/cGuqI3FxcWqa5V6KIOsmYkJqcpTVIdOplc0O4myS4Cx7tNbZ/16cikWEUh56LVCXw+nk0zQNvcAvmQOLDd63ueaccEnO43VlVEsVUK8H4fTKxbWR8cKT6Ett5XQEk6H5dASTWwgu55rBbeqKu3AqYEvIYQ++eGoTCodXhaTHU47y9P7aJRBY0FRUVglBPAJYtaAQ0Mbv9vKysguAWvFlNMptO9RgdU2PJY5vPVhPaO8f39qEDgeNWA0uNHKhIdoNt5jhVvUOVoSzjkrmlykYWcrhicsCyPHEjP3tcuRjMd5JZpnOe3dzAvL47GsrCxqpzSzdYuFMxPpGMaK4eKUo5Jb3jUNvRy5PLi5bnoI5ywdzC6LnVgZFBcXx3TFEU35eHB69hfNfZyybjCylLHaFmaX6vcdOgRPj8AHumf5nkBA1zIsWeB5F25focSThBP6pIVLBeM4Cblhppyvt2HG0+FKSkpMp7WL5VJQTwBp62HHZjnrA4ynV6f2PSoOcyxLmVgbDsTTrNGtagkzRGMimywkjNBXZmoevx9yT95XHocpcgWwgnGchPdjNLvisOtjturE4YQgdkrvatVjmHyPipNcRc/feS1l3CI47CyH2y1nzFpbCb8ANQkh9LXxcCRZDgt+IwsXpxKPO57QnIFe5iwr2DX745l18Q4STn+0rAGabAurg6RdsXdELBzrXt3JQkIIfa2KRgKQ0rM/rVjHsF781V6v6dR6CnofY2+yqedVWdlxTyB+nrfMGE2ZmdTjrAHazrK6IfZOb0b7DUZrIptMJITQJ2dkFVDnqDQSKrsaGkyl1iPR+xhVTlxwNpmCHZuMdgtiu2KO2DlLYwlYs3ln3aIj5gm5kQyY8eqO5QDq1veQECabLDNNCTA0zyJNuJwy54qlmVi05oROpFs0KhNPfgBerJqq2m22q4dd7U3W2c4+ZlesfCfQe89kG7ihDna+E2GyqYGlSmHBmqlZSQTem3F6oyuerva8s3IzTm9WcdvGorYtetP+AGs2HW+Bb4QIrWwRs6EN9GLB8MAT/dAJFYBbl4tGxDPmSLwGtFhjZcLSm4S8FrLsbngXvN99PNs8IdQ7JGaXzmaXYKSnZAW69xCM8uLatczTU5PYpZ7hyWvqRpzyRHYaK2oIu9Q7Pp8PxcXFjqr47MKqR24s0Sub2XIL9Y4OZmd2ZmcHsY5+yIOds1nejE+9EbPvWi+ZiF3wvCPW8yZPnhz1c7UzUkXgK20Ub5WTHm7ZSDeLW8rNNdPfv38/nnrqKYRCIUyaNAnTpk1T/X3fvn149tlnIUkSPB4PbrnlFlx00UWGD++NidGHFhRA6mkyJSYKAMiShMajR6nXaGdzZjYNnV652H19vHAiDo1RbKVYoVcOGmZn59r7u7kP9KaZvt4qztUz/VAohI0bN2LJkiVYu3Ytdu7ciaMa4fZv//ZvWL16NVavXo25c+fi8ccft1Qou3AiMXo00Q9ZyTF4YvKYncVHU2e3zEAE5nFiv8KuFVGy48ZVEsAh9BsaGjBkyBDk5eUhNTUVEyZMwL59+1TnpKenQ+rxgO3o6Aj/P944kRi9aOVKBL1eVY7boNeLopUrua7nifVj5SNatmyZ6TR7brMmiSduSdrNKsfYsWOZdunRlkkr5M2+91huSmrDjPRWi7t4bjobqnf27NmD/fv344477gAA7NixAx9//DFuu+021Xl79+7Ff/7nf+LUqVNYvHgxvvnNb0bca+vWrdi6dSsAoKamBufOnbOrHlRSU1PR1dVl+31TnnkGnuXLkXL4MEKFhQhWViJ0/fVc16alpzPVQ53t7QCAvn37oqOjI+Laqqoqw4Es9bnnIN1+O6S2tvAxOSMDwcce4yoj69lux4l3zWqLWLcRTzliXSayL8by2dpn9db+ahayf/fp08favYxOoI0JtJn82LFjMXbsWLz//vt49tlnqcJp8uTJqs0nJ2KBkzgVbxxTpgBTpqB4zRo0KrMizucMzs+nOpIF8/NVZaWVe+7cuYb1Gbp0qUrgA+j+felSnJwyxbB8xcXFjr8XJ3DqXbPuGes24ilHLMt03333qeI89e3bF0BsrH209eyN/dUsMdXp5+TkoKmpKfx7U1MTsrOzmedfcskl+OKLL3D69GlLBesNRNO5l44aRd0TuNrrtUd9cOQI9TBvmr1kVOmQ8OxvxHppznrelVdeaXitk/p2M5EurcB6J2IPKjoMhf6IESPQ2NiIEydOoKurC7t27cKYMWNU53zxxRfhFcGnn36Krq4uZDKCWCU78zdsoMb6eX77dnv06sOGUQ87Hd3TTVj5+Hn2N2IdHI71vKuvvtpQ+Nmpb9cLoufE8xRY70TsQUWHodD3eDyYNWsWqqurUVRUhPHjx2PYsGGor69HfX09gG69f0lJCRYuXIiNGzeiqKjINZu5bkRx+Qdgu8t/sLKy10T3NItWmLsxOYtd8NRh2bJljgg/vdAGtOfZGXZb4Dxczlnf+ta38NBDD2HdunWY3iOgpk6diqlTpwIApk2bhtraWqxevRrV1dVcNvqxxmyHikUHJD8Wu8wmQ9dfTw1J4WQsmVihFYROC3c3uPUD0fVdK30pmna1o+/yPFdrkuyWd9SbSMgwDMpMx4qbejydPqw827HNaxeg50S0fv16agIZsxuL8QxDwXLEA+jRYsl3bTZiqVXHITsjpJq5NpH7tx4iDIMOibC0F5yHdxMPsMdRKZ5hKMymAdReawbtd2J2dWCX/lw4BsaehBP6JGaTTLilA4ol63l4N/GccMSLN044XjkN63lm9wnEpqxzJIR6h7UsLisrw9y5cwH0LvWOFRJ5+aun3uFRdbCYMWMGdYY/bty4uKp6lDqw+qLZd80bx8mJvs9zT6HeYWOnesfQOas3UFJSQv1AkrWDJCraWaRdyVlIwW5W4DmVocmpexoNJPFErHBjQ0Krd0jMdqhE64CJoCPVCkI3qABisYfkFhWNWeK1T9DbcGVoZadwwnon2Wf6rDq7cWZnFZallhXMWu/Es12t1DnWOWTtaqdE/KbNqrWE9Y4gaXFils0j8N2y4W+FZJ1VCxJE6LM+Qm0Y1mSjtwont5cPcCaOvVl6U/9ONHWpVeL5bQr1ToJhVb0T62U/Db2y8lhqOQ1Zvt7qxNdbScRvWqh3BHHF7c5tLLvueNnp885gyRmc07M5q/fvDSstQfQknNCP5TLSyaxJduP25bXbVVE8ESb1IAdTKwMrTztZHbjdPvAnGrH+NhNOvUPi9FLQqeW1FVvxaOpsNvm600Sjiorlsj+a9+6ESoh1H6v3d7PaKBHVOzwI9Y5AhdWZWW91hXd7+fRWB3avZKyulJxeabllxSYQQt80Tn0cbldvxBK3q6K0IbFZ6AVQs2NgJcMMWx24nR74hcrIPQj1jgXiqd5xyorFDdY70RCvZT9vH3BCveOUI168Yu/wINQ7SaDe8dbVYfDYsRhaUIDBY8fCW1cX7yK5AqesWHqjwO8NkKsDp1cyVu8frzAMgtjg6pm+t64OA++9FymBQPhYyOvlzgTl9KzAqVmx2fsmu28CENt6u2Xjuze9azHTt4adM31XC/3BY8cildJRuny+cI5ZPZKlg8TLisVNuF294wS96V0LoW+NpFHveBiDAut4ssIzuxRLavOINrMPt2/OJxOuFvpBxojGOi5gk6zWE1YEN0+bCWHGh9grcg+uFvqtpaUIeb2qYyGvF62lpXEqkaC34fRgJ4SZoLfBJfT379+P+fPnY968eXjppZci/v73v/8dCxYswIIFC1BWVobPPvvMlsIFpk/HqQceQJfPB1mS0OXzcW/iCoT1RDSINhMkOoYbuaFQCPPnz0dZWRlycnKwePFizJ8/HwUFBeFzPvzwQ/h8PvTv3x//+7//i+effx4rV640fHhvt9N3I8mURIVEmyPXDusat7eZ6N/JQ0xz5DY0NGDIkCHIy8sDAEyYMAH79u1TCf0LL7ww/P9vfOMbaGpqslQogcAKbs8FKxDEE0Oh39zcjJycnPDvOTk5+Pjjj5nn/+1vf8Po0aOpf9u6dSu2bt0KAKipqUFubq7Z8poiNTXV8We4DVady8rKErot9N51tPV2e5uJ/p082FlvQ6FP0/5IkkQ998CBA3j99ddRWVlJ/fvkyZMxefLk8O9OL9OScSnIqvPcuXMTui1Y9S4uLo663m5vM9G/k4eY2unn5OSo1DVNTU3Izs6OOO/zzz/HE088gYULFyIzM9NSoQQCuxDWNfYhNrMTA0OhP2LECDQ2NuLEiRPo6urCrl27MGbMGNU5J0+exIMPPoi7777b8igkEAjcSbL6eiQahuodj8eDWbNmobq6GqFQCBMnTsSwYcNQX18PAJg6dSpeeOEFnDlzBhs2bAhfU1NT42zJBQKBQGAaV8fesUoy6v+Ssc5ActY7FnV2S3A5hWR8z0CMTTYFAkHyIsxfEw9Xh2EQCAQCgb0IoS8QCLgQweUSAyH0BQIBF8L8NTEQQl8gEAiSCCH0BQKBIIkQQl8gIBBep4JERwh9gYBAeJ0KEh0h9AUCgSCJEEJfkPSIbFmCZEJ45AqSHuF1KkgmxExfIBAIkggh9AUCAuF1Kkh0hNAXCAiE16kg0RFCXyAQCJKIuMbTFwgEAkFsSeiZfmlpabyLEHOSsc5ActZb1Dl5sLPeCS30BQKBQKBGCH2BQCBIIhJa6E+ePDneRYg5yVhnIDnrLeqcPNhZb7GRKxAIBElEQs/0BQKBQKBGCH2BQCBIInpVwLXHHnsMb7/9NgYOHBiOgPjZZ5/hd7/7Hc6dOwePx4PZs2fjggsuwLvvvovNmzejq6sLqampuPHGG3HppZcCAD799FM8+uijOHfuHEaPHo1bb70VkiTFs2q6mKm3wsmTJ1FUVISZM2fipz/9KYDeVW+zdf7888/x5JNPIhAIQJIkrFq1Cn369EnYOnd1deHxxx/HoUOHEAqFcNVVV+G6664D0LveM6Bf7/b2dnzlK1/BPffcg4yMDADAiy++iL/97W9ISUnBrbfeissvvxxA76q3mTrbLsvkXsTBgwflTz75RC4uLg4fq6qqkt9++21ZlmX5rbfeksvLy2VZluVPP/1UbmpqkmVZlj///HN5zpw54WtKS0vlDz/8UA6FQnJ1dXX4erdipt4Kq1evltesWSO//PLL4WO9qd5m6tzV1SWXlJTIhw4dkmVZlk+fPi0Hg0FZlhO3zn//+9/ltWvXyrIsy+3t7fKdd94pHz9+XJbl3lVnWabXu7S0VD548KAsy7K8bds2+ZlnnpFlWZaPHDkiL1iwQD537px8/Phx+e67706Yd82qs92yrFepdy655BL0799fdUySJAQCAQBAW1sbsrOzAQDDhw/HoEGDAADDhg1DZ2cnOjs78eWXXyIQCOCb3/wmJEnCVVddhX379sW2IiYxU28A2Lt3L/Ly8lBQUBA+1tvqbabO77zzDgoLC/G1r30NAJCZmYmUlJSErjMAtLe3IxgM4ty5c0hNTUVGRkavqzNAr/exY8dw8cUXAwBGjRqFN998EwCwb98+TJgwAWlpaRg8eDCGDBmChoaGXldvM3W2W5b1KqFP4+abb8bTTz+NuXPn4umnn8avfvWriHPefPNNDB8+HGlpaWhubkZOTk74bzk5OWhubo5lkW2BVe/29na8/PLLmDlzpur8RKg3q86NjY2QJAnV1dVYtGgRXn75ZQCJXedx48YhPT0dc+bMwZ133olrr70W/fv3T4g6A93C7Z///CcAYM+ePWhqagIQ+U4HDRqE5ubmhKg3q84kdsiyXi/06+vrcfPNN2P9+vW4+eab8fjjj6v+fuTIEWzevBn/8R//AQCQE8RClVXv5557Dj/+8Y+Rnp6uOj8R6s2qczAYxAcffIB58+ahsrISe/fuxXvvvZfQdW5oaEBKSgqeeOIJPPLII/jv//5vHD9+PCHqDABz587FX//6VyxatAiBQACpqd3bj6z6JUK9WXVWsEuW9Xqhv337dnz3u98FAIwfPx4NDQ3hvzU1NeHBBx/EXXfdhSFDhgDoHg3JEbSpqSm8dOpNsOrd0NCAzZs346677sKf//xnvPjii3j11VcTot6sOufk5OCSSy7BgAED0LdvX4wePRqHDh1K6Dr/4x//wOWXX47U1FQMHDgQF154IT755JOEqDPQncGsrKwM999/P773ve8hLy8PQOT329zcjEGDBiVEvVl1BuyVZb1e6A8aNAjvv/8+AODAgQPhBjl79ixqampw/fXX46KLLgqfn52dDa/Xi48++giyLGPHjh0YM2ZMXMpuBVa9Kysr8eijj+LRRx/Fj370I1x33XX44Q9/mBD1ZtX5sssuw+HDh9HR0YFgMIj/+7//Q0FBQULXOTc3FwcOHIAsy2hvb8fHH38Mn8+XEHUGgFOnTgEAQqEQ6urqMGXKFADAmDFjsGvXLnR2duLEiRNobGzEBRdckBD1ZtXZblnWqzxyf/vb3+L9999Ha2srBg4ciF/84hfIz8/HU089hVAohLS0NMyePRtf//rXsWXLFrz00kvhjwQAysrKMHDgQHzyySd47LHHcO7cOVx++eWYNWuWa027AHP1JnnuueeQnp4eNtnsTfU2W+cdO3bgpZdegiRJGD16NH79618DSNw6t7e347HHHsPRo0chyzImTpzYK98zQK93e3s7/vrXvwIAxo4di1/96lfhOtTV1eH1119HSkoKbrnlFowePRpA76q3mTrbLct6ldAXCAQCgTV6vXpHIBAIBPwIoS8QCARJhBD6AoFAkEQIoS8QCARJhBD6AoFAkEQIoS8QCARJhBD6AoFAkET8fwYhNMZ3JANhAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "If we divide the dataset with a horizontal line at 0.5, accuracy is:  0.7736111111111111\n",
      "Divided with a line fit to the data trend, it's  0.7902777777777777\n"
     ]
    }
   ],
   "source": [
    "### DEFACTORING FUNCTION PARAMETERS\n",
    "modeltype = 'linear'\n",
    "datelimits = []\n",
    "''' Takes a set of predictions produced by a model that knows nothing about date,\n",
    "and divides it along a line with a diachronic tilt. We need to do this in a way\n",
    "that doesn't violate crossvalidation. I.e., we shouldn't \"know\" anything\n",
    "that the model didn't know. We tried a couple of different ways to do this, but\n",
    "the simplest and actually most reliable is to divide the whole dataset along a\n",
    "linear central trend line for the data!\n",
    "'''\n",
    "\n",
    "\n",
    "\n",
    "listofrows = list()\n",
    "classvector = list()\n",
    "\n",
    "# DEPRECATED\n",
    "# if modeltype == 'logistic' and len(datelimits) == 2:\n",
    "#     # In this case we construct a subset of data to model on.\n",
    "#     tomodeldata = list()\n",
    "#     tomodelclasses = list()\n",
    "#     pastthreshold, futurethreshold = datelimits\n",
    "\n",
    "for volume in allvolumes:\n",
    "    date = volume[3]\n",
    "    logistic = volume[8]\n",
    "    realclass = volume[13]\n",
    "    listofrows.append([logistic, date])\n",
    "    classvector.append(realclass)\n",
    "\n",
    "    # DEPRECATED\n",
    "    # if modeltype == 'logistic' and len(datelimits) == 2:\n",
    "    #     if date >= pastthreshold and date <= futurethreshold:\n",
    "    #         tomodeldata.append([logistic, date])\n",
    "    #         tomodelclasses.append(realclass)\n",
    "\n",
    "y, x = [a for a in zip(*listofrows)]\n",
    "plt.axis([min(x) - 2, max(x) + 2, min(y) - 0.02, max(y) + 0.02])\n",
    "reviewedx = list()\n",
    "reviewedy = list()\n",
    "randomx = list()\n",
    "randomy = list()\n",
    "\n",
    "for idx, reviewcode in enumerate(classvector):\n",
    "    if reviewcode == 1:\n",
    "        reviewedx.append(x[idx])\n",
    "        reviewedy.append(y[idx])\n",
    "    else:\n",
    "        randomx.append(x[idx])\n",
    "        randomy.append(y[idx])\n",
    "\n",
    "plt.plot(reviewedx, reviewedy, 'ro')\n",
    "plt.plot(randomx, randomy, 'k+')\n",
    "\n",
    "if modeltype == 'logistic':\n",
    "    # all this is DEPRECATED\n",
    "    print(\"Hey, you're attempting to use the logistic-tilt option\")\n",
    "    print(\"that we deactivated. Go in and uncomment the code.\")\n",
    "\n",
    "    # if len(datelimits) == 2:\n",
    "    #     data = pd.DataFrame(tomodeldata)\n",
    "    #     responsevariable = tomodelclasses\n",
    "    # else:\n",
    "    #     data = pd.DataFrame(listofrows)\n",
    "    #     responsevariable = classvector\n",
    "\n",
    "    # newmodel = LogisticRegression(C = 100000)\n",
    "    # newmodel.fit(data, responsevariable)\n",
    "    # coefficients = newmodel.coef_[0]\n",
    "\n",
    "    # intercept = newmodel.intercept_[0] / (-coefficients[0])\n",
    "    # slope = coefficients[1] / (-coefficients[0])\n",
    "\n",
    "    # p = np.poly1d([slope, intercept])\n",
    "\n",
    "elif modeltype == 'linear':\n",
    "    # what we actually do\n",
    "\n",
    "    z = np.polyfit(x, y, 1)\n",
    "    p = np.poly1d(z)\n",
    "    slope = z[0]\n",
    "    intercept = z[1]\n",
    "\n",
    "plt.plot(x,p(x),\"b-\")\n",
    "plt.show(block = False)\n",
    "\n",
    "x = np.array(x, dtype='float64')\n",
    "y = np.array(y, dtype='float64')\n",
    "classvector = np.array(classvector)\n",
    "dividingline = intercept + (x * slope)\n",
    "predicted_as_reviewed = (y > dividingline)\n",
    "really_reviewed = (classvector == 1)\n",
    "\n",
    "accuracy = sum(predicted_as_reviewed == really_reviewed) / len(classvector)\n",
    "\n",
    "### DEFACTORING NAMESPACE\n",
    "tiltaccuracy = accuracy \n",
    "\n",
    "print('If we divide the dataset with a horizontal line at 0.5, accuracy is: ', \n",
    "      str(rawaccuracy))\n",
    "\n",
    "print(\"Divided with a line fit to the data trend, it's \", \n",
    "      str(tiltaccuracy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code above \"linear regression on the whole data set to identify a central trend\" (Underwood & Sellers 2016, 328). Then selects the predicted-as-reviewed volumes, which are the volumes falling above the dividing line. The accuracy of the dividing line is the ratio between the total number of volumes predicted-as-reviewed and the number of actually reviewed volumes. Note that this is the accuracy of the dividing line, not the accuracy of the model(s) to predict the reviewed status of the poems. A much more aesthetically pleasing figure was generated using the MLQ article, presumably using ggplot2 and R, which is a easier tool for generating figures for publication than matplotlib and Python. The published figure in the MLQ article below closely resembles the one generated in python above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Chart from the MLQ article](notebook_resources/mlq-chart.jpg)\n",
    "\n",
    "Underwood and Seller's discuss the tilt of this line at length in the MLQ article:\n",
    "\n",
    "> Technically, the upward drift is an error in the model. Volumes are not really more likely to be reviewed just because they were published later. But this is an error of an interesting kind, since the upward drift suggests that historical change across this century moved in the direction of the standards that govern reception for the century as a whole (Underwood & Sellers 2016, 328).\n",
    "\n",
    "The existence of this trend and that the pattern persists in temporal subsets is the significant discovery in the *Longue Duree of Literary Prestige:*\n",
    "\n",
    "> The words common in reviewed volumes (across the whole timeline) also tend to become more common in all volumes as we approach the end of the timeline, so the cloud of data points always tilts upward. We have compensated for this by allowing the dividing line between categories to tilt upward as well (Underwood & Sellers 2016, 335).\n",
    "\n",
    "What we \"discover\" in our defactoring is the path to their discovery. We follow their code-work and the traces of their methodology in the data and the code to actually bear witness to the material practices with which knowledge is produced.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}