quora-question-similarity/quora-question-similarity.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Doc2Vec for Kaggle Quora competition\n",
    "\n",
    "Here we're tackling Kaggle's Quora question pairs competition mostly using the Paragraph Vector idea by [Le and Mikolov](https://arxiv.org/pdf/1405.4053v2.pdf). \n",
    "\n",
    "As training data we get ~400K question pairs labeled 1 if the questions in the pair are duplicates/have the same intent, and labeled 0 if they are considered to be different.\n",
    "\n",
    "We then have to classify a test set of ~2.3M question pairs. The classification should be as correct as possible since the score is log loss, which punishes how wrong guesses are.\n",
    "\n",
    "TL;DR: we get __best results with tuned XGBoost__ (surprise!) against a very basic, untuned deep net. Also, Doc2Vec vectorization sets a performance ceiling.\n",
    "\n",
    "![Ask More Questions poster](../notebooks-img/ask_more_questions_part.jpeg?raw=True \"Ask More Questions poster\")\n",
    "\n",
    "Unlike many approaches currently discussed in the competition [Kernels](https://www.kaggle.com/c/quora-question-pairs/kernels), which take into account the [distribution of labels](https://www.kaggle.com/anokas/data-analysis-xgboost-starter-0-35460-lb) [in the test data](https://www.kaggle.com/davidthaler/how-many-1-s-are-in-the-public-lb), for the code below I don't take into account any information other than the training data.\n",
    "\n",
    "See [the Kaggle Quora competition webiste](https://www.kaggle.com/c/quora-question-pairs) for full details and for the data!\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start by importing libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "import sys, os\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from numpy.linalg import norm\n",
    "from six.moves import cPickle as pickle\n",
    "from collections import Counter\n",
    "from gensim.models.doc2vec import Doc2Vec\n",
    "from gensim import matutils\n",
    "import multiprocessing\n",
    "from scipy.stats import expon, bernoulli, geom\n",
    "from itertools import izip\n",
    "from scipy.spatial.distance import cdist\n",
    "import time\n",
    "from xgboost import XGBClassifier\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import log_loss\n",
    "from skopt import gbrt_minimize\n",
    "from skopt.space import Categorical\n",
    "from sklearn.cross_validation import cross_val_score\n",
    "import tensorflow as tf\n",
    "\n",
    "# note that I'm using a special gensim that allows use of Doc2Vec with pretrained word embeddings\n",
    "# you can get it at https://github.com/jhlau/gensim\n",
    "# most of the code is the same if using regular gensim, except for small syntactical fixes\n",
    "# also, if using regular gensim, don't use the pretrained_emb parameter in Doc2Vec. the difference isn't huge.\n",
    "import gensim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data management\n",
    "\n",
    "\n",
    "\n",
    "We then define the function to read in the csv file and process the question pairs. We also have another function that takes a dataframe with question pairs with tags and reduces them to a list of unique questions for further processing.\n",
    "\n",
    "Note the pickling of variables along the way. I favor this approach in case I need to step away and resume processing later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def processRows(tokens_only=False, train_size=None):\n",
    "\n",
    "    if tokens_only: # test data\n",
    "        pickle_filename = \"pickles/test_processed_df.pickle\"\n",
    "    else: # train data\n",
    "        # can extend this if-statement in case we want to \n",
    "        # work with fractions of the data\n",
    "        if train_size is None: # full data, ~404K\n",
    "            pickle_filename = \"pickles/train_processed_df.pickle\"\n",
    "        elif train_size == \"1000\": # 0.25%, ~1000 rows\n",
    "            pickle_filename = \"pickles/train_processed_df_1000.pickle\"\n",
    "    \n",
    "    # load data from pickle if we've previously processed it\n",
    "    if os.path.exists(pickle_filename):\n",
    "        \n",
    "        print('%s already present - skipping pickling.' % pickle_filename)\n",
    "        with open(pickle_filename, 'rb') as f:\n",
    "            data_df = pickle.load(f)\n",
    "\n",
    "    # otherwise, process and pickle\n",
    "    else:\n",
    "\n",
    "        if tokens_only:\n",
    "            data_df = pd.read_csv(\"test.csv\")\n",
    "        else:\n",
    "            data_df = pd.read_csv(\"train.csv\")\n",
    "\n",
    "        print(\"Processing training rows\" if not tokens_only else \"Processing test rows\")\n",
    "        for index, row in data_df.iterrows():\n",
    "            \n",
    "            # lowercases, tokenizes\n",
    "            q1 = gensim.utils.simple_preprocess(str(row.question1))\n",
    "            q2 = gensim.utils.simple_preprocess(str(row.question2))\n",
    "\n",
    "            if tokens_only:\n",
    "                data_df.set_value(index, \"question1\", q1)\n",
    "                data_df.set_value(index, \"question2\", q2)\n",
    "            else:\n",
    "                # for training data, add tags\n",
    "                # using doc2vec.TaggedDocument for Doc2Vec later on\n",
    "                data_df.set_value(index, \"question1\", gensim.models.doc2vec.TaggedDocument(q1, [row.qid1]))\n",
    "                data_df.set_value(index, \"question2\", gensim.models.doc2vec.TaggedDocument(q2, [row.qid2]))\n",
    "\n",
    "\n",
    "        try:\n",
    "            with open(pickle_filename, 'wb') as f:\n",
    "                pickle.dump(data_df, f, pickle.HIGHEST_PROTOCOL)\n",
    "        except Exception as e:\n",
    "            print('Unable to save data to', pickle_filename, ':', e)\n",
    "\n",
    "    return data_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def getUniqueTrainQuestions(train_processed_df=None):\n",
    "    # since there's ~808K questions but only ~537K are unique, \n",
    "    # consider only unique ones to make training faster\n",
    "\n",
    "    print(\"getting unique questions for data\")\n",
    "\n",
    "    if train_processed_df is None:\n",
    "        train_processed_df = processRows()\n",
    "\n",
    "    qid_set = set()\n",
    "    q_list = []\n",
    "\n",
    "    for q1, q2 in izip(train_processed_df.question1, train_processed_df.question2):\n",
    "        if q1.tags[0] not in qid_set:\n",
    "            q_list.append(q1)\n",
    "            qid_set.add(q1.tags[0])\n",
    "        if q2.tags[0] not in qid_set:\n",
    "            q_list.append(q2)\n",
    "            qid_set.add(q2.tags[0])\n",
    "\n",
    "    return q_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Doc2Vec (+ random search for near-optimal parameters)\n",
    "\n",
    "With very lightly processed text, we can start vectorizing sentences. \n",
    "\n",
    "There's many ways of doing this. Using TF-IDF, for instance, is an obvious choice. There's a few other newer vectorizing methods that I've been meaning to try, like Doc2Vec and its Paragraph Vector idea by [Le and Mikolov](https://arxiv.org/pdf/1405.4053v2.pdf). One of its main insights is that \"The paragraph token can be thought of as another word. It acts as a memory that remembers what is missing from the current context – or the topic of the paragraph.\"\n",
    "\n",
    "Doc2Vec is a Paragraph Vector implementation in gensim. It has many parameters that I used random search to comb through, according to this paper from [Bergstra and Bengio](http://www.jmlr.org/papers/v13/bergstra12a.html) on how, compared to grid search, \"random search over the same domain is able to find models that are as good or better within a small fraction of the computation time.\" (Admittedly, after having spent some time doing grid search.)\n",
    "\n",
    "Here's a rough, brief and very illuminating image (from the paper itself) explaining why. \n",
    "\n",
    "![Random Search vs. Grid Search](../notebooks-img/bergstra_bengio_grid_vs_random.png?raw=True \"Random Search vs. Grid Search\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below is a somewhat inelegant function to perform random search for the best parameters for Doc2Vec, as well as a sanity check to assess the quality of the models, taken from [gensim's own code](https://github.com/RaRe-Technologies/gensim/blob/develop/docs/notebooks/doc2vec-lee.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def sanity(data, model, iterations):\n",
    "\n",
    "    print(\"generating sanity check on\", len(data), \"docs\")\n",
    "    ranks = []\n",
    "    second_ranks = []\n",
    "\n",
    "    # check the same 1000 random documents to have comparable model assessments\n",
    "    np.random.seed(0)\n",
    "    doc_id_sample = np.random.choice(xrange(len(data)), 1000)\n",
    "\n",
    "    # iterate over documents\n",
    "    for i, doc_id in enumerate(doc_id_sample):\n",
    "        \n",
    "        sys.stdout.write(\"\\rProcessing doc #{:d} (doc_id = {:d}) of 1000 sampled docs\".format(i, doc_id))\n",
    "        sys.stdout.flush()\n",
    "        \n",
    "        # infer vector from document's words\n",
    "        inferred_vector = model.infer_vector(data[doc_id].words, steps=iterations)\n",
    "        \n",
    "        # take most similar documents to inferred vector\n",
    "        sims = model.docvecs.most_similar([inferred_vector], topn=len(model.docvecs))\n",
    "        \n",
    "        # index of current doc on sorted list of documents most similar to inferred vector\n",
    "        # should be 0 or close to 0\n",
    "        rank = [docid for docid, sim in sims].index(doc_id)\n",
    "        \n",
    "        # collect index\n",
    "        ranks.append(rank)\n",
    "\n",
    "    print()\n",
    "    print(Counter(ranks).most_common(5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that in ``randomSearch()`` we do 100 iterations, which may be too many. According to [Alice Zheng from Dato](https://stats.stackexchange.com/a/209409) (then Turi, which was acquired by Apple), \"if the close-to-optimal region of hyperparameters occupies at least 5% of the grid surface, then random search with 60 trials will find that region with high probability.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def randomSearch(data, train_size, n_iter=100):\n",
    "\n",
    "    n_jobs = multiprocessing.cpu_count()\n",
    "\n",
    "    # only train in unique questions\n",
    "    print(\"pre-unique questions:\", data.shape[0]*2)\n",
    "    data = getUniqueTrainQuestions(data)\n",
    "    print(\"post-unique questions:\", len(data))\n",
    "    \n",
    "    # sampling at different orders of magnitude. manually tuned.\n",
    "    dimensions = np.ceil(expon.rvs(scale=1.5, size=n_iter)).astype(int)*100\n",
    "    min_count = np.ceil(expon.rvs(scale=1.5, size=n_iter)).astype(int)\n",
    "    iterations = np.ceil(expon.rvs(scale=23, size=n_iter)).astype(int)*10\n",
    "    window_size = geom.rvs(0.085, loc=2, size=100)\n",
    "    sampling_threshold = expon.rvs(scale=2e-4, loc=1e-8, size=n_iter)\n",
    "    dm_concat = bernoulli.rvs(p=0.5, size=n_iter)\n",
    "    dm = bernoulli.rvs(p=0.5, size=n_iter)\n",
    "    dbow_words = bernoulli.rvs(p=0.5, size=n_iter)\n",
    "    negative = np.ceil(expon.rvs(scale=5, loc=-1, size=n_iter)).astype(int)\n",
    "\n",
    "    # visually ensure we have enough range\n",
    "    # I was personally looking for sampling_threshold the \n",
    "    # min in the e-07 range, and max in the e-03 range.\n",
    "    print(\"min(sampling_threshold)\", min(sampling_threshold))\n",
    "    print(\"max(sampling_threshold)\", max(sampling_threshold))\n",
    "\n",
    "    # create search space list\n",
    "    random_search_list = [(dimensions[i], min_count[i], iterations[i], window_size[i], \n",
    "                           sampling_threshold[i], dm_concat[i], dm[i], dbow_words[i], \n",
    "                           negative[i]) for i in xrange(n_iter)]\n",
    "\n",
    "    # iterate over search space list\n",
    "    for dimensions, min_count, iterations, window_size, sampling_threshold, \\\n",
    "        dm_concat, dm, dbow_words, negative in random_search_list:\n",
    "\n",
    "        print(\"generating model\")\n",
    "        print(\"train_size\", train_size, end=\" | \")\n",
    "        print(\"dimensions\", dimensions, end=\" | \")\n",
    "        print(\"min_count\", min_count, end=\" | \")\n",
    "        print(\"iterations\", iterations, end=\" | \")\n",
    "        print(\"window_size\", window_size, end=\" | \")\n",
    "        print(\"sampling_threshold\", sampling_threshold, end=\" | \")\n",
    "        print(\"dm_concat\", dm_concat, end=\" | \")\n",
    "        print(\"dm\", dm, end=\" | \")\n",
    "        print(\"dbow_words\", dbow_words, end=\" | \")\n",
    "        print(\"negative\", negative, end=\"\\n\")\n",
    "\n",
    "\n",
    "        model = Doc2Vec(size=dimensions, min_count=min_count, iter=iterations, \n",
    "                        window=window_size, workers=n_jobs, sample=sampling_threshold, \n",
    "                        dm_concat=dm_concat, dm=dm, dbow_words=dbow_words, negative=negative)\n",
    "\n",
    "        \n",
    "        model.build_vocab(data)\n",
    "        print(\"training started on model\")\n",
    "        # for regular gensim:\n",
    "        #    model.train(data, total_examples=model.corpus_count, epochs=model.iter)\n",
    "        model.train(data)\n",
    "        \n",
    "        sanity(data, model, iterations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a demonstration, let's run the code so far on only 1K documents (as opposed to ~400K)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pickles/train_processed_df_1000.pickle already present - skipping pickling.\n",
      "pre-unique questions: 2020\n",
      "getting unique questions for data\n",
      "post-unique questions: 2014\n",
      "min(sampling_threshold) 1.18748078328e-05\n",
      "max(sampling_threshold) 0.000437384587252\n",
      "generating model\n",
      "train_size 1000 | dimensions 200 | min_count 4 | iterations 520 | window_size 7 | sampling_threshold 5.3497446733e-05 | dm_concat 1 | dm 1 | dbow_words 1 | negative 0\n",
      "training started on model\n",
      "generating sanity check on 2014 docs\n",
      "Processing doc #999 (doc_id = 707) of 1000 sampled docss\n",
      "[(1, 104), (0, 40), (2, 28), (3, 17), (8, 12)]\n",
      "generating model\n",
      "train_size 1000 | dimensions 100 | min_count 4 | iterations 180 | window_size 8 | sampling_threshold 0.000215620894966 | dm_concat 1 | dm 1 | dbow_words 1 | negative 0\n",
      "training started on model\n",
      "generating sanity check on 2014 docs\n",
      "Processing doc #999 (doc_id = 707) of 1000 sampled docss\n",
      "[(1, 39), (2, 20), (0, 19), (4, 10), (5, 9)]\n",
      "generating model\n",
      "train_size 1000 | dimensions 200 | min_count 3 | iterations 110 | window_size 5 | sampling_threshold 2.06616947804e-05 | dm_concat 0 | dm 0 | dbow_words 0 | negative 2\n",
      "training started on model\n",
      "generating sanity check on 2014 docs\n",
      "Processing doc #999 (doc_id = 707) of 1000 sampled docss\n",
      "[(3, 16), (0, 14), (4, 10), (1, 9), (2, 9)]\n",
      "generating model\n",
      "train_size 1000 | dimensions 300 | min_count 1 | iterations 20 | window_size 18 | sampling_threshold 4.63524612343e-05 | dm_concat 1 | dm 1 | dbow_words 0 | negative 6\n",
      "training started on model\n",
      "generating sanity check on 2014 docs\n",
      "Processing doc #999 (doc_id = 707) of 1000 sampled docss\n",
      "[(6, 4), (43, 4), (162, 4), (367, 4), (463, 4)]\n",
      "generating model\n",
      "train_size 1000 | dimensions 200 | min_count 1 | iterations 140 | window_size 20 | sampling_threshold 0.000100941050888 | dm_concat 1 | dm 0 | dbow_words 1 | negative 0\n",
      "training started on model\n",
      "generating sanity check on 2014 docs\n",
      "Processing doc #999 (doc_id = 707) of 1000 sampled docss\n",
      "[(1, 287), (2, 30), (3, 19), (0, 12), (4, 11)]\n",
      "generating model\n",
      "train_size 1000 | dimensions 200 | min_count 1 | iterations 160 | window_size 36 | sampling_threshold 0.000107799872776 | dm_concat 0 | dm 1 | dbow_words 0 | negative 0\n",
      "training started on model\n",
      "generating sanity check on 2014 docs\n",
      "Processing doc #999 (doc_id = 707) of 1000 sampled docss\n",
      "[(1, 19), (0, 6), (3, 6), (2, 4), (21, 4)]\n",
      "generating model\n",
      "train_size 1000 | dimensions 200 | min_count 1 | iterations 50 | window_size 35 | sampling_threshold 0.000437384587252 | dm_concat 1 | dm 0 | dbow_words 1 | negative 2\n",
      "training started on model\n",
      "generating sanity check on 2014 docs\n",
      "Processing doc #999 (doc_id = 707) of 1000 sampled docss\n",
      "[(1, 301), (2, 32), (3, 19), (0, 10), (4, 9)]\n",
      "generating model\n",
      "train_size 1000 | dimensions 200 | min_count 3 | iterations 430 | window_size 8 | sampling_threshold 1.18748078328e-05 | dm_concat 1 | dm 1 | dbow_words 1 | negative 4\n",
      "training started on model\n",
      "generating sanity check on 2014 docs\n",
      "Processing doc #999 (doc_id = 707) of 1000 sampled docss\n",
      "[(1, 102), (0, 35), (2, 31), (3, 26), (4, 16)]\n",
      "generating model\n",
      "train_size 1000 | dimensions 100 | min_count 3 | iterations 210 | window_size 23 | sampling_threshold 0.000255110214606 | dm_concat 0 | dm 0 | dbow_words 1 | negative 2\n",
      "training started on model\n",
      "generating sanity check on 2014 docs\n",
      "Processing doc #999 (doc_id = 707) of 1000 sampled docss\n",
      "[(1, 228), (2, 36), (3, 31), (0, 25), (4, 10)]\n",
      "generating model\n",
      "train_size 1000 | dimensions 300 | min_count 1 | iterations 100 | window_size 5 | sampling_threshold 2.5518819046e-05 | dm_concat 1 | dm 1 | dbow_words 1 | negative 3\n",
      "training started on model\n",
      "generating sanity check on 2014 docs\n",
      "Processing doc #999 (doc_id = 707) of 1000 sampled docss\n",
      "[(1, 49), (2, 15), (0, 13), (3, 10), (4, 9)]\n"
     ]
    }
   ],
   "source": [
    "train_size = \"1000\"\n",
    "data = processRows(train_size = train_size)\n",
    "randomSearch(data, train_size, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When I ran this on the training set, I got the following best parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dimensions = 300 \n",
    "min_count = 1 \n",
    "iterations = 300 \n",
    "window_size = 17 # average of top performers\n",
    "sampling_threshold = 0.000232331506591 # median of top performers\n",
    "dm_concat = 1 \n",
    "dm = 0 \n",
    "dbow_words = 1\n",
    "negative = 4 # median of top performers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Doc2Vec and word embeddings\n",
    "\n",
    "Then I wanted to see the role of pretrained word embeddings on training. I used the following function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def explorePretrained(data):\n",
    "\n",
    "    print(\"pre-unique questions:\", data.shape[0]*2)\n",
    "    data = getUniqueTrainQuestions(data)\n",
    "    print(\"post-unique questions:\", len(data))\n",
    "\n",
    "    n_jobs = multiprocessing.cpu_count()\n",
    "\n",
    "    # fetched mostly through https://github.com/3Top/word2vec-api\n",
    "    for pretrained_emb in [ \"word_emb/toy_pretrained_word_embeddings.txt\",\n",
    "                            \"word_emb/glove.840B.300d.modified_header.txt\", \n",
    "                            \"word_emb/glove.twitter.27B/glove.twitter.27B.100d.modified_header.txt\", \n",
    "                            \"word_emb/glove.twitter.27B/glove.twitter.27B.200d.modified_header.txt\", \n",
    "                            \"word_emb/glove.twitter.27B/glove.twitter.27B.25d.modified_header.txt\", \n",
    "                            \"word_emb/glove.twitter.27B/glove.twitter.27B.50d.modified_header.txt\", \n",
    "                            \"word_emb/glove.wikiplusgigaword5.6B/glove.6B.100d.modified_header.txt\", \n",
    "                            \"word_emb/glove.wikiplusgigaword5.6B/glove.6B.200d.modified_header.txt\", \n",
    "                            \"word_emb/glove.wikiplusgigaword5.6B/glove.6B.300d.modified_header.txt\", \n",
    "                            \"word_emb/glove.wikiplusgigaword5.6B/glove.6B.50d.modified_header.txt\",\n",
    "                            \"word_emb/freebase-vectors-skipgram1000-en-without-en.txt\",\n",
    "                            \"word_emb/GoogleNews-vectors-negative300.txt\"]:\n",
    "\n",
    "        try:\n",
    "\n",
    "            print(\"pretrained_emb\", pretrained_emb)\n",
    "\n",
    "            # note that the values of the parameters are assigned above e.g. dimensions = 300\n",
    "            model = Doc2Vec(size=dimensions, min_count=min_count, iter=iterations, \n",
    "                    window=window_size, workers=n_jobs, sample=sampling_threshold, \n",
    "                    dm_concat=dm_concat, dm=dm, dbow_words=dbow_words, negative=negative,\n",
    "                    pretrained_emb=pretrained_emb)\n",
    "\n",
    "            sanity(data, model, iterations)\n",
    "\n",
    "        except Exception as e:\n",
    "            print(\"failed with\", pretrained_emb, \"with exception\", e)\n",
    "            \n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "explorePretrained(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "My results on the whole data indicated that ``glove.twitter.27B/glove.twitter.27B.50d.modified_header.txt`` yielded the best results.\n",
    "\n",
    "![Word Embedding experiments results](../notebooks-img/word_embed_results.png?raw=True \"Word Embedding experiments results\")\n",
    "\n",
    "Note that the results aren't spectacular, however. The ``sanity()`` function picks 1000 random documents and tries to make sure that each document's already-known vector is close to the vector inferred by the model purely from the document's words. \n",
    "\n",
    "Here, as with the ``sanity()`` assessment for the random search, we see about 20-30% of the 1000 docs taking up one of the top 5 most similar spots to their own inferred vector. Ideally, we'd like close to 100% of the documents taking the single top similar spot. This indicated that perhaps Doc2Vec vectorization, as powerful as it is, has its limits as currently implemented."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training Doc2Vec\n",
    "\n",
    "Now that we have the right hyper-parameters to train our model, we can use ``trainModel()`` below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def trainModel():\n",
    "\n",
    "    data = processRows() # get all training data\n",
    "    data = getUniqueTrainQuestions(data) # unique training data\n",
    "\n",
    "    n_jobs = multiprocessing.cpu_count()\n",
    "\n",
    "    dimensions = 300 \n",
    "    min_count = 1 \n",
    "    iterations = 300 \n",
    "    window_size = 17 \n",
    "    sampling_threshold = 0.000232331506591\n",
    "    dm_concat = 1 \n",
    "    dm = 0 \n",
    "    dbow_words = 1\n",
    "    negative = 4\n",
    "    pretrained_emb = \"word_emb/glove.twitter.27B/glove.twitter.27B.50d.modified_header.txt\"\n",
    "\n",
    "    try:\n",
    "\n",
    "        model = Doc2Vec(size=dimensions, min_count=min_count, iter=iterations, \n",
    "                window=window_size, workers=n_jobs, sample=sampling_threshold, \n",
    "                dm_concat=dm_concat, dm=dm, dbow_words=dbow_words, negative=negative,\n",
    "                pretrained_emb=pretrained_emb)\n",
    "\n",
    "\n",
    "        start_time = time.time()    \n",
    "        model.build_vocab(data)\n",
    "        print(time.time()-start_time, \"secs to build vocab\")\n",
    "\n",
    "        start_time = time.time()    \n",
    "        print(\"training started on model\")\n",
    "        model.train(data)\n",
    "        print(time.time()-start_time, \"secs to train\")\n",
    "\n",
    "        pickle_filename = \"pickles/best_model.pickle\"\n",
    "\n",
    "        try:\n",
    "            with open(pickle_filename, 'wb') as f:\n",
    "                pickle.dump(model, f, pickle.HIGHEST_PROTOCOL)\n",
    "        except Exception as e:\n",
    "            print('Unable to save data to', pickle_filename, ':', e)\n",
    "            \n",
    "        return model\n",
    "\n",
    "    except Exception as e:\n",
    "        print(\"failed with exception\", e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pickles/train_processed_df.pickle already present - skipping pickling.\n",
      "getting unique questions for data\n",
      "23.2628879547 secs to build vocab\n",
      "training started on model\n",
      "28809.9361 secs to train"
     ]
    }
   ],
   "source": [
    "model = trainModel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## (Efficiently) vectorizing millions of sentences\n",
    "\n",
    "Now that we've trained our Doc2Vec model to (likely) close to its best possible performance in vectorizing documents, we want to vectorize all test questions.\n",
    "\n",
    "This process can be _extremely slow_. The process starts really fast but then seems to be hindered by memory. In order to avoid that, I split the data in chunks of 10K question pairs at a time and processed them separately. This brought down an expected processing time of several days to about ~7h on my laptop!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def partitionData(data):\n",
    "\n",
    "    partitions = np.split(data, xrange(10000,data.shape[0],10000), axis=0)\n",
    "\n",
    "    for i, p in enumerate(partitions):\n",
    "\n",
    "        pickle_filename = \"pickles/test_processed/test_processed_df_\"+str(i+1)+\"-\"+str(len(partitions))+\".pickle\"\n",
    "\n",
    "        try:\n",
    "            with open(pickle_filename, 'wb') as f:\n",
    "                pickle.dump(p, f, pickle.HIGHEST_PROTOCOL)\n",
    "        except Exception as e:\n",
    "            print('Unable to save data to', pickle_filename, ':', e)\n",
    "\n",
    "    return len(partitions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "num_partitions = partitionData(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I processed these text vectors in parallel manually. Each line defining a for-loop in ``inferTestVectorsDivideAndConquer()`` below used its own processor.\n",
    "\n",
    "Making this code nicer as a proper python-multiprocessing function is a project for a future afternoon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def inferTestVectorsDivideAndConquer(num_partitions, model):\n",
    "\n",
    "    dimensions = 300 \n",
    "\n",
    "    for number in [1, 8, 15, 22, 29, 36, 43, 50, 57, 64, 71, 78, 85, 92, 99, 106, 113, 120, 127, 134, 141, 148, 155, 162, 169, 176, 183, 190, 197, 204, 211, 218, 225, 232, 239, 246]:\n",
    "    # for number in [2, 9, 16, 23, 30, 37, 44, 51, 58, 65, 72, 79, 86, 93, 100, 107, 114, 121, 128, 135, 142, 149, 156, 163, 170, 177, 184, 191, 198, 205, 212, 219, 226, 233]:\n",
    "    # for number in [3, 10, 17, 24, 31, 38, 45, 52, 59, 66, 73, 80, 87, 94, 101, 108, 115, 122, 129, 136, 143, 150, 157, 164, 171, 178, 185, 192, 199, 206, 213, 220, 227, 234]:\n",
    "    # for number in [4, 11, 18, 25, 32, 39, 46, 53, 60, 67, 74, 81, 88, 95, 102, 109, 116, 123, 130, 137, 144, 151, 158, 165, 172, 179, 186, 193, 200, 207, 214, 221, 228, 235]:\n",
    "    # for number in [5, 12, 19, 26, 33, 40, 47, 54, 61, 68, 75, 82, 89, 96, 103, 110, 117, 124, 131, 138, 145, 152, 159, 166, 173, 180, 187, 194, 201, 208, 215, 222, 229]:\n",
    "    # for number in [6, 13, 20, 27, 34, 41, 48, 55, 62, 69, 76, 83, 90, 97, 104, 111, 118, 125, 132, 139, 146, 153, 160, 167, 174, 181, 188, 195, 202, 209, 216, 223, 230]:\n",
    "    # for number in [7, 14, 21, 28, 35, 42, 49, 56, 63, 70, 77, 84, 91, 98, 105, 112, 119, 126, 133, 140, 147, 154, 161, 168, 175, 182, 189, 196, 203, 210, 217, 224, 231]:\n",
    "\n",
    "\n",
    "        pickle_filename = \"pickles/test_processed/test_processed_df_\" + str(number) +\"-\" + num_partitions + \".pickle\"\n",
    "        output_pickle = \"pickles/inferred/inferred_\" + str(number) +\"-\" + num_partitions + \".pickle\"\n",
    "\n",
    "        with open(pickle_filename, 'rb') as f:\n",
    "            data_df = pickle.load(f)\n",
    "\n",
    "        print(\"inferring test vectors for \" + output_pickle)\n",
    "\n",
    "        start_time = time.time()\n",
    "\n",
    "        q1_vec = pd.DataFrame(columns=range(0,dimensions))\n",
    "        q2_vec = pd.DataFrame(columns=range(0,dimensions))\n",
    "\n",
    "        for index, row in data_df.iterrows():\n",
    "\n",
    "            sys.stdout.write(\"\\rProcessing doc #{:d} of {:d} docs\".format(index, data_df.shape[0]))\n",
    "            sys.stdout.flush()\n",
    "                \n",
    "            q1_vec.loc[index] = model.infer_vector(row.question1, steps=300)\n",
    "            q2_vec.loc[index] = model.infer_vector(row.question2, steps=300)\n",
    "\n",
    "        print()\n",
    "        print(\"file, q1_vec.shape, q2_vec.shape\")\n",
    "        print(output_pickle, q1_vec.shape, q2_vec.shape)\n",
    "\n",
    "        print(time.time()-start_time, \"sec to infer test vectors for\", output_pickle)\n",
    "\n",
    "        q_vec = {\"q1_vec\": q1_vec, \"q2_vec\": q2_vec}\n",
    "\n",
    "        start_time = time.time()\n",
    "        try:\n",
    "            with open(output_pickle, 'wb') as f:\n",
    "                pickle.dump(q_vec, f, pickle.HIGHEST_PROTOCOL)\n",
    "        except Exception as e:\n",
    "            print('Unable to save data to', output_pickle, ':', e)\n",
    "\n",
    "        print(time.time()-start_time, \"sec to write inferred test vectors for\", output_pickle)\n",
    "\n",
    "    return q1_vec, q2_vec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "inferTestVectorsDivideAndConquer(num_partitions, model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given that the test data is now partitioned, we have to put it back together:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def uniteTestData():\n",
    "    \n",
    "    output_pickle = \"pickles/inferred_test_vectors.pickle\"\n",
    "    \n",
    "    if os.path.exists(output_pickle):\n",
    "        \n",
    "        print('%s already present - skipping pickling.' % pickle_filename)\n",
    "        with open(output_pickle, 'rb') as f:\n",
    "            data_df = pickle.load(f)\n",
    "            \n",
    "            q1_vec = data[\"q1_vec\"]\n",
    "            q2_vec = data[\"q2_vec\"]\n",
    "    \n",
    "    else:\n",
    "        \n",
    "        start_time = time.time()\n",
    "    \n",
    "        for number in xrange(1,236):\n",
    "\n",
    "            print(number, end=\"/235...\")\n",
    "\n",
    "            input_pickle = \"pickles/inferred/inferred_\" + str(number) +\"-235.pickle\"\n",
    "\n",
    "            with open(input_pickle, 'rb') as f:\n",
    "                data_df = pickle.load(f)\n",
    "\n",
    "                if number == 1:\n",
    "                    q1_vec = data_df[\"q1_vec\"].copy()\n",
    "                    q2_vec = data_df[\"q2_vec\"].copy()\n",
    "                else:\n",
    "                    q1_vec = pd.concat([q1_vec, data_df[\"q1_vec\"]])\n",
    "                    q2_vec = pd.concat([q2_vec, data_df[\"q2_vec\"]])\n",
    "\n",
    "        \n",
    "        print(time.time()-start_time, \"sec to gather inferred test vectors into\", output_pickle)\n",
    "                    \n",
    "        print()\n",
    "\n",
    "        q_vec = {\"q1_vec\": q1_vec, \"q2_vec\": q2_vec}\n",
    "\n",
    "        start_time = time.time()\n",
    "\n",
    "        try:\n",
    "            with open(output_pickle, 'wb') as f:\n",
    "                pickle.dump(q_vec, f, pickle.HIGHEST_PROTOCOL)\n",
    "        except Exception as e:\n",
    "            print('Unable to save data to', output_pickle, ':', e)\n",
    "\n",
    "        print(time.time()-start_time, \"sec to write inferred test vectors for\", output_pickle)\n",
    "                \n",
    "    return q1_vec, q2_vec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1/235...2/235...3/235...4/235...5/235...6/235...7/235...8/235...9/235...10/235...11/235...12/235...13/235...14/235...15/235...16/235...17/235...18/235...19/235...20/235...21/235...22/235...23/235...24/235...25/235...26/235...27/235...28/235...29/235...30/235...31/235...32/235...33/235...34/235...35/235...36/235...37/235...38/235...39/235...40/235...41/235...42/235...43/235...44/235...45/235...46/235...47/235...48/235...49/235...50/235...51/235...52/235...53/235...54/235...55/235...56/235...57/235...58/235...59/235...60/235...61/235...62/235...63/235...64/235...65/235...66/235...67/235...68/235...69/235...70/235...71/235...72/235...73/235...74/235...75/235...76/235...77/235...78/235...79/235...80/235...81/235...82/235...83/235...84/235...85/235...86/235...87/235...88/235...89/235...90/235...91/235...92/235...93/235...94/235...95/235...96/235...97/235...98/235...99/235...100/235...101/235...102/235...103/235...104/235...105/235...106/235...107/235...108/235...109/235...110/235...111/235...112/235...113/235...114/235...115/235...116/235...117/235...118/235...119/235...120/235...121/235...122/235...123/235...124/235...125/235...126/235...127/235...128/235...129/235...130/235...131/235...132/235...133/235...134/235...135/235...136/235...137/235...138/235...139/235...140/235...141/235...142/235...143/235...144/235...145/235...146/235...147/235...148/235...149/235...150/235...151/235...152/235...153/235...154/235...155/235...156/235...157/235...158/235...159/235...160/235...161/235...162/235...163/235...164/235...165/235...166/235...167/235...168/235...169/235...170/235...171/235...172/235...173/235...174/235...175/235...176/235...177/235...178/235...179/235...180/235...181/235...182/235...183/235...184/235...185/235...186/235...187/235...188/235...189/235...190/235...191/235...192/235...193/235...194/235...195/235...196/235...197/235...198/235...199/235...200/235...201/235...202/235...203/235...204/235...205/235...206/235...207/235...208/235...209/235...210/235...211/235...212/235...213/235...214/235...215/235...216/235...217/235...218/235...219/235...220/235...221/235...222/235...223/235...224/235...225/235...226/235...227/235...228/235...229/235...230/235...231/235...232/235...233/235...234/235...235/235...5867.43962097 sec to gather inferred test vectors into pickles/inferred/inferred.pickle\n",
      "\n",
      "Unable to save data to pickles/inferred/inferred.pickle : error return without exception set\n",
      "69.5035700798 sec to write inferred test vectors for pickles/inferred/inferred.pickle\n"
     ]
    }
   ],
   "source": [
    "q1_vec, q2_vec = uniteTestData()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expressed [here](http://stackoverflow.com/a/38246020/583834), ``cPickle`` has its limits when it comes to large files. We use instead the pandas builtin HDF5 format output function ``to_hdf()``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1658.35788798 sec to write inferred test vectors for pickles/inferred_test_vectors.h5\n"
     ]
    }
   ],
   "source": [
    "output_pickle = \"pickles/inferred_test_vectors.h5\"\n",
    "\n",
    "start_time = time.time()\n",
    "\n",
    "q_vec = pd.concat([q1_vec, q2_vec], ignore_index=True, axis=1)\n",
    "\n",
    "try:\n",
    "    q_vec.to_hdf(output_pickle, 'q_vec', mode='w')\n",
    "    del q1_vec\n",
    "    del q2_vec\n",
    "except Exception as e:\n",
    "    print('Unable to save data to', output_pickle, ':', e)\n",
    "\n",
    "print(time.time()-start_time, \"sec to write inferred test vectors for\", output_pickle)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A quick note that this will be a bit heavy on your memory..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2345796, 600)\n",
      "<class 'pandas.core.frame.DataFrame'>\n",
      "Int64Index: 2345796 entries, 0 to 2345795\n",
      "Columns: 600 entries, 0 to 599\n",
      "dtypes: float64(600)\n",
      "memory usage: 10.5 GB\n",
      "None\n"
     ]
    }
   ],
   "source": [
    "print(q_vec.shape)\n",
    "print(q_vec.info())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## (Unsupervised) cosine similarity predictions\n",
    "\n",
    "From here we can submit a prediction simply by obtaining a measure of similarity between the first and the second question for each pair. We calculate it [the same way gensim's ``doc2vec`` does](https://github.com/RaRe-Technologies/gensim/blob/develop/gensim/models/doc2vec.py), through a cosine similarity. \n",
    "\n",
    "We implement cosine similarity code because Doc2Vec doesn't yet support calculating similarity with provided vectors. Also, it's surprising that row-wise cosine distance for matrices is not implemented in any python library (that I'm aware of). [``einsum``](https://docs.scipy.org/doc/numpy-1.12.0/reference/generated/numpy.einsum.html) is used, thanks to [this post](http://stackoverflow.com/a/15622926/583834). [This other post](http://stackoverflow.com/a/33641428/583834) is a good primer on ``einsum``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "245.740324974 sec to generate 2345796 predictions\n"
     ]
    }
   ],
   "source": [
    "q1_vec = q_vec.loc[:,:299]\n",
    "q2_vec = q_vec.loc[:,300:]\n",
    "\n",
    "start_time = time.time()\n",
    "\n",
    "# cosine similarity\n",
    "cosine_predictions = np.einsum('ij,ij->i', q1_vec, q2_vec) / (norm(q1_vec, axis=1) * norm(q2_vec, axis=1))\n",
    "\n",
    "# replace all values less than 0 for 0\n",
    "cosine_predictions[cosine_predictions < 0] = 0\n",
    "\n",
    "# make dataframe for output\n",
    "cosine_predictions = pd.concat([pd.Series(np.arange(cosine_predictions.shape[0])), pd.Series(cosine_predictions)], axis=1)\n",
    "\n",
    "# rename dataframe's columns\n",
    "cosine_predictions.columns = [\"test_id\",\"is_duplicate\"]\n",
    "\n",
    "# write csv submission\n",
    "cosine_predictions.to_csv(\"similarity_predictions.csv\", index=False)\n",
    "\n",
    "print(time.time()-start_time, \"sec to generate\", cosine_predictions.shape[0], \"predictions\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results are not optimal, since a good score would be closer to 0.4 and below:\n",
    "\n",
    "![Paragraph Vector unsupervised similarity score: 0.70604](../notebooks-img/doc2vec_unsupervised_score_wide.png?raw=True \"Paragraph Vector unsupervised similarity score: 0.70604\")\n",
    "\n",
    "This score would likely be improved by running Doc2Vec over both training and test data, since we only have ~400K question pairs in training and ~2.34M pairs in testing.\n",
    "\n",
    "## Supervised training methods\n",
    "\n",
    "## XGBoost (untuned)\n",
    "\n",
    "We can also improve the score by using XGBoost on the training data and predicting the test data from it. In order to switch to a supervised training solution, we want to treat each question pair as a single vector. For this, we use the 600-dimension concatenation of the 300-dimension vectors of both questions in each pair. Some other operations such as addition or subtraction on the two 300-dimenstional vectors _could_ work, but usually vectors are concatenated to avoid data loss.\n",
    "\n",
    "Here's a basic XGBoost application. First, we load the data (because I restarted the server since I ran the code above...):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pickles/train_processed_df.pickle already present - skipping pickling.\n",
      "(404290, 600)\n",
      "(404290,)\n"
     ]
    }
   ],
   "source": [
    "pickle_filename = \"pickles/best_model.pickle\"\n",
    "\n",
    "# training data and model\n",
    "train_data = processRows()\n",
    "labels = train_data.is_duplicate\n",
    "\n",
    "with open(pickle_filename, 'rb') as f:\n",
    "    model = pickle.load(f)\n",
    "\n",
    "# vectorized questions\n",
    "q1_vec = model.docvecs[train_data.qid1]\n",
    "q2_vec = model.docvecs[train_data.qid2]\n",
    "q_vec = pd.concat([pd.DataFrame(q1_vec), pd.DataFrame(q2_vec)], ignore_index=True, axis=1)\n",
    "\n",
    "print(q_vec.shape)\n",
    "print(labels.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we perform XGBoost training and extract predictions: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "874.805408001 sec to train XGBoost model\n",
      "Log Loss: 10.880096\n"
     ]
    }
   ],
   "source": [
    "test_size = 0.25\n",
    "X_train, X_test, y_train, y_test = train_test_split(q_vec, labels, test_size=test_size)\n",
    "\n",
    "start_time = time.time()\n",
    "\n",
    "xgb_model = XGBClassifier(nthread=-1)\n",
    "xgb_model.fit(X_train, y_train)\n",
    "\n",
    "print(time.time()-start_time, \"sec to train XGBoost model\")\n",
    "\n",
    "# pickle model\n",
    "output_pickle = \"pickles/xgboost_model.pickle\"\n",
    "try:\n",
    "    with open(output_pickle, 'wb') as f:\n",
    "        pickle.dump(xgb_model, f, pickle.HIGHEST_PROTOCOL)\n",
    "except Exception as e:\n",
    "    print('Unable to save data to', output_pickle, ':', e)\n",
    "\n",
    "# make and report predictions scores for test data\n",
    "predictions = xgb_model.predict(X_test)\n",
    "logloss_score = log_loss(y_test, predictions)\n",
    "print(\"Log Loss: {:f}\".format(logloss_score))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This result is _extremely_ bad, but ameliorated by the fact that ``model.predict`` outputs label values, while ``model.predict_proba`` outputs probabilities. Since the ``log_loss`` metric penalizes wrong results, getting at least close to the right result is important, and probabilities help us with that. We use ``predict_proba()`` this time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Log Loss: 0.571065\n"
     ]
    }
   ],
   "source": [
    "predictions = xgb_model.predict_proba(X_test)\n",
    "logloss_score = log_loss(y_test, predictions)\n",
    "print(\"Log Loss: {:f}\".format(logloss_score))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This result is much more reasonable. Now we can generate predictions for the test values and submit them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2345796, 600)\n",
      "529.756315231 sec to generate 2345796 predictions\n"
     ]
    }
   ],
   "source": [
    "qtest_vec = pd.read_hdf(\"pickles/inferred_test_vectors.h5\")\n",
    "\n",
    "print(qtest_vec.shape)\n",
    "\n",
    "start_time = time.time()\n",
    "predictions = xgb_model.predict_proba(qtest_vec)\n",
    "print(time.time()-start_time, \"sec to generate\", predictions.shape[0], \"predictions\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to use only the second column of the prediction result because, [as the source reveals](http://xgboost/python-package/xgboost/sklearn.py), a horizontal stack of ``classzero_probs`` and ``classone_probs`` is returned."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# prepare for submission\n",
    "predictions = pd.concat([pd.Series(np.arange(predictions.shape[0])), pd.Series(predictions[:,1])], axis=1)\n",
    "predictions.columns = [\"test_id\",\"is_duplicate\"]\n",
    "predictions.to_csv(\"xgboost_predictions.csv\", index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result of the submission is a bit worse than that of the similarity predictions:\n",
    "\n",
    "![Vanilla (untuned) XGBoost score: 0.82654](../notebooks-img/vanilla_xgboost_score_wide.png?raw=True \"Vanilla (untuned) XGBoost score: 0.82654\")\n",
    "\n",
    "## Cross-validated XGBoost\n",
    "\n",
    "We can improve this by tuning the XGBoost parameters with cross validation. The code below is taken from Bernstein and and Potts' excellent [Optimizing the hyperparameter of which hyperparameter optimizer to use](http://roamanalytics.com/2016/09/15/optimizing-the-hyperparameter-of-which-hyperparameter-optimizer-to-use/) [\\[github\\]](https://github.com/roaminsight/roamresearch/tree/master/BlogPosts/Hyperparameter_tuning_comparison).\n",
    "\n",
    "We first define the cross validation functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cross_validated_scorer(X_train, y_train, model_class, params, loss, kfolds=5):\n",
    "    \"\"\"\n",
    "    The scoring function used through this module, by all search\n",
    "    functions.\n",
    "    \"\"\"\n",
    "    mod = model_class(**params)\n",
    "    cv_score = -1 * cross_val_score(\n",
    "        mod,\n",
    "        X_train,\n",
    "        y=y_train,\n",
    "        scoring=loss,\n",
    "        cv=kfolds,\n",
    "        n_jobs=multiprocessing.cpu_count()).mean()\n",
    "    return cv_score\n",
    "\n",
    "def skopt_search(X_train, y_train, model_class, param_grid, loss, skopt_method, n_calls=100):\n",
    "    \"\"\"\n",
    "    General method for applying `skopt_method` to the data.\n",
    "    \"\"\"\n",
    "    param_keys, param_vecs = zip(*param_grid.items())\n",
    "    param_keys = list(param_keys)\n",
    "    param_vecs = list(param_vecs)\n",
    "\n",
    "    def skopt_scorer(param_vec):\n",
    "        params = dict(zip(param_keys, param_vec))\n",
    "        err = cross_validated_scorer(\n",
    "            X_train, y_train, model_class, params, loss)\n",
    "        return err\n",
    "    outcome = skopt_method(skopt_scorer, list(param_vecs), n_calls=n_calls, verbose=True)\n",
    "    results = []\n",
    "    for err, param_vec in zip(outcome.func_vals, outcome.x_iters):\n",
    "        params = dict(zip(param_keys, param_vec))\n",
    "        results.append({'loss': err, 'params': params})\n",
    "    return results\n",
    "\n",
    "def skopt_gbrt_search(\n",
    "        X_train, y_train, model_class, param_grid, loss, n_calls=100):\n",
    "    \"\"\"\n",
    "    `skopt` according to the gradient-boosting-tree search method.\n",
    "    \"\"\"\n",
    "    return skopt_search(\n",
    "        X_train, y_train, model_class, param_grid, loss, gbrt_minimize, n_calls=n_calls)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that \"``scikit-optimize`` asks you to specify just the upper and lower bounds of the space to be searched.\"\n",
    "\n",
    "We define these bounds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "skopt_grid = {\n",
    "    'max_depth': (4, 12),\n",
    "    'learning_rate': (0.01, 0.5),\n",
    "    'n_estimators': (20, 200),\n",
    "    'objective' : Categorical(('binary:logistic',)),\n",
    "    'gamma': (0, 0.5),\n",
    "    'min_child_weight': (1, 5),\n",
    "    'subsample': (0.1, 1),\n",
    "    'colsample_bytree': (0.1, 1)}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also note that, at least for my AWS machine, I had an issue with the package multiprocessing v0.70a1 and pandas v0.20.1 similar to the one [here](https://github.com/scikit-learn/scikit-learn/issues/7981). I reverted to pandas v0.19.2, which I knew was working in another machine, and everything ran well. For the output, I had to tweak [``gbrt.py``](https://github.com/scikit-optimize/scikit-optimize/blob/master/skopt/optimizer/gbrt.py) due to a [small bug](https://github.com/scikit-optimize/scikit-optimize/issues/326) in passing the ``verbose`` parameter. It should be resolved by scikit-optimize v0.4.\n",
    "\n",
    "We now run the cross validation (and wait...)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration No: 1 started. Evaluating function at random point.\n",
      "Iteration No: 1 ended. Evaluation done at random point.\n",
      "Time taken: 4465.8883\n",
      "Function value obtained: 0.5287\n",
      "Current minimum: 0.5287\n",
      "Iteration No: 2 started. Evaluating function at random point.\n",
      "Iteration No: 2 ended. Evaluation done at random point.\n",
      "Time taken: 544.0995\n",
      "Function value obtained: 0.5615\n",
      "Current minimum: 0.5287\n",
      "Iteration No: 3 started. Evaluating function at random point.\n",
      "Iteration No: 3 ended. Evaluation done at random point.\n",
      "Time taken: 1703.3487\n",
      "Function value obtained: 0.5177\n",
      "Current minimum: 0.5177\n",
      "Iteration No: 4 started. Evaluating function at random point.\n",
      "Iteration No: 4 ended. Evaluation done at random point.\n",
      "Time taken: 1754.8015\n",
      "Function value obtained: 0.5727\n",
      "Current minimum: 0.5177\n",
      "Iteration No: 5 started. Evaluating function at random point.\n",
      "Iteration No: 5 ended. Evaluation done at random point.\n",
      "Time taken: 1233.1887\n",
      "Function value obtained: 0.5124\n",
      "Current minimum: 0.5124\n",
      "Iteration No: 6 started. Evaluating function at random point.\n",
      "Iteration No: 6 ended. Evaluation done at random point.\n",
      "Time taken: 3502.9155\n",
      "Function value obtained: 0.5019\n",
      "Current minimum: 0.5019\n",
      "Iteration No: 7 started. Evaluating function at random point.\n",
      "Iteration No: 7 ended. Evaluation done at random point.\n",
      "Time taken: 208.8112\n",
      "Function value obtained: 0.5315\n",
      "Current minimum: 0.5019\n",
      "Iteration No: 8 started. Evaluating function at random point.\n",
      "Iteration No: 8 ended. Evaluation done at random point.\n",
      "Time taken: 514.5910\n",
      "Function value obtained: 0.4954\n",
      "Current minimum: 0.4954\n",
      "Iteration No: 9 started. Evaluating function at random point.\n",
      "Iteration No: 9 ended. Evaluation done at random point.\n",
      "Time taken: 516.5554\n",
      "Function value obtained: 0.5347\n",
      "Current minimum: 0.4954\n",
      "Iteration No: 10 started. Evaluating function at random point.\n",
      "Iteration No: 10 ended. Evaluation done at random point.\n",
      "Time taken: 252.5022\n",
      "Function value obtained: 0.5520\n",
      "Current minimum: 0.4954\n",
      "Iteration No: 11 started. Searching for the next optimal point.\n",
      "Iteration No: 11 ended. Search finished for the next optimal point.\n",
      "Time taken: 747.3653\n",
      "Function value obtained: 0.5880\n",
      "Current minimum: 0.4954\n",
      "Iteration No: 12 started. Searching for the next optimal point.\n",
      "Iteration No: 12 ended. Search finished for the next optimal point.\n",
      "Time taken: 447.4338\n",
      "Function value obtained: 0.5019\n",
      "Current minimum: 0.4954\n",
      "Iteration No: 13 started. Searching for the next optimal point.\n",
      "Iteration No: 13 ended. Search finished for the next optimal point.\n",
      "Time taken: 1051.1582\n",
      "Function value obtained: 0.4938\n",
      "Current minimum: 0.4938\n",
      "Iteration No: 14 started. Searching for the next optimal point.\n",
      "Iteration No: 14 ended. Search finished for the next optimal point.\n",
      "Time taken: 1178.8049\n",
      "Function value obtained: 0.4925\n",
      "Current minimum: 0.4925\n",
      "Iteration No: 15 started. Searching for the next optimal point.\n",
      "Iteration No: 15 ended. Search finished for the next optimal point.\n",
      "Time taken: 412.5154\n",
      "Function value obtained: 0.4981\n",
      "Current minimum: 0.4925\n",
      "Iteration No: 16 started. Searching for the next optimal point.\n",
      "Iteration No: 16 ended. Search finished for the next optimal point.\n",
      "Time taken: 1067.9992\n",
      "Function value obtained: 0.4678\n",
      "Current minimum: 0.4678\n",
      "Iteration No: 17 started. Searching for the next optimal point.\n",
      "Iteration No: 17 ended. Search finished for the next optimal point.\n",
      "Time taken: 1076.8978\n",
      "Function value obtained: 0.5003\n",
      "Current minimum: 0.4678\n",
      "Iteration No: 18 started. Searching for the next optimal point.\n",
      "Iteration No: 18 ended. Search finished for the next optimal point.\n",
      "Time taken: 2339.6743\n",
      "Function value obtained: 0.4500\n",
      "Current minimum: 0.4500\n",
      "Iteration No: 19 started. Searching for the next optimal point.\n",
      "Iteration No: 19 ended. Search finished for the next optimal point.\n",
      "Time taken: 855.9505\n",
      "Function value obtained: 0.4688\n",
      "Current minimum: 0.4500\n",
      "Iteration No: 20 started. Searching for the next optimal point.\n",
      "Iteration No: 20 ended. Search finished for the next optimal point.\n",
      "Time taken: 2343.4761\n",
      "Function value obtained: 0.4690\n",
      "Current minimum: 0.4500\n",
      "Iteration No: 21 started. Searching for the next optimal point.\n",
      "Iteration No: 21 ended. Search finished for the next optimal point.\n",
      "Time taken: 2010.8827\n",
      "Function value obtained: 0.7040\n",
      "Current minimum: 0.4500\n",
      "Iteration No: 22 started. Searching for the next optimal point.\n",
      "Iteration No: 22 ended. Search finished for the next optimal point.\n",
      "Time taken: 1047.5082\n",
      "Function value obtained: 0.7326\n",
      "Current minimum: 0.4500\n",
      "Iteration No: 23 started. Searching for the next optimal point.\n",
      "Iteration No: 23 ended. Search finished for the next optimal point.\n",
      "Time taken: 3835.9812\n",
      "Function value obtained: 0.4876\n",
      "Current minimum: 0.4500\n",
      "Iteration No: 24 started. Searching for the next optimal point.\n",
      "Iteration No: 24 ended. Search finished for the next optimal point.\n",
      "Time taken: 1067.2297\n",
      "Function value obtained: 0.5205\n",
      "Current minimum: 0.4500\n",
      "Iteration No: 25 started. Searching for the next optimal point.\n",
      "Iteration No: 25 ended. Search finished for the next optimal point.\n",
      "Time taken: 2496.8922\n",
      "Function value obtained: 0.4419\n",
      "Current minimum: 0.4419\n",
      "Iteration No: 26 started. Searching for the next optimal point.\n",
      "Iteration No: 26 ended. Search finished for the next optimal point.\n",
      "Time taken: 611.4180\n",
      "Function value obtained: 1.0712\n",
      "Current minimum: 0.4419\n",
      "Iteration No: 27 started. Searching for the next optimal point.\n",
      "Iteration No: 27 ended. Search finished for the next optimal point.\n",
      "Time taken: 1571.0556\n",
      "Function value obtained: 2.1142\n",
      "Current minimum: 0.4419\n",
      "Iteration No: 28 started. Searching for the next optimal point.\n",
      "Iteration No: 28 ended. Search finished for the next optimal point.\n",
      "Time taken: 655.3671\n",
      "Function value obtained: 0.6329\n",
      "Current minimum: 0.4419\n",
      "Iteration No: 29 started. Searching for the next optimal point.\n",
      "Iteration No: 29 ended. Search finished for the next optimal point.\n",
      "Time taken: 4351.4282\n",
      "Function value obtained: 0.5632\n",
      "Current minimum: 0.4419\n",
      "Iteration No: 30 started. Searching for the next optimal point.\n",
      "Iteration No: 30 ended. Search finished for the next optimal point.\n",
      "Time taken: 185.3983\n",
      "Function value obtained: 0.5635\n",
      "Current minimum: 0.4419\n",
      "Iteration No: 31 started. Searching for the next optimal point.\n",
      "Iteration No: 31 ended. Search finished for the next optimal point.\n",
      "Time taken: 479.1749\n",
      "Function value obtained: 0.7071\n",
      "Current minimum: 0.4419\n",
      "Iteration No: 32 started. Searching for the next optimal point.\n",
      "Iteration No: 32 ended. Search finished for the next optimal point.\n",
      "Time taken: 3992.0674\n",
      "Function value obtained: 0.5514\n",
      "Current minimum: 0.4419\n",
      "Iteration No: 33 started. Searching for the next optimal point.\n",
      "Iteration No: 33 ended. Search finished for the next optimal point.\n",
      "Time taken: 2283.0408\n",
      "Function value obtained: 0.5263\n",
      "Current minimum: 0.4419\n",
      "Iteration No: 34 started. Searching for the next optimal point.\n",
      "Iteration No: 34 ended. Search finished for the next optimal point.\n",
      "Time taken: 276.8985\n",
      "Function value obtained: 0.6364\n",
      "Current minimum: 0.4419\n",
      "Iteration No: 35 started. Searching for the next optimal point.\n",
      "Iteration No: 35 ended. Search finished for the next optimal point.\n",
      "Time taken: 970.2984\n",
      "Function value obtained: 0.6137\n",
      "Current minimum: 0.4419\n",
      "Iteration No: 36 started. Searching for the next optimal point.\n",
      "Iteration No: 36 ended. Search finished for the next optimal point.\n",
      "Time taken: 2100.6165\n",
      "Function value obtained: 0.5099\n",
      "Current minimum: 0.4419\n",
      "Iteration No: 37 started. Searching for the next optimal point.\n",
      "Iteration No: 37 ended. Search finished for the next optimal point.\n",
      "Time taken: 4367.9989\n",
      "Function value obtained: 0.5174\n",
      "Current minimum: 0.4419\n",
      "Iteration No: 38 started. Searching for the next optimal point.\n",
      "Iteration No: 38 ended. Search finished for the next optimal point.\n",
      "Time taken: 1208.8153\n",
      "Function value obtained: 1.2260\n",
      "Current minimum: 0.4419\n",
      "Iteration No: 39 started. Searching for the next optimal point.\n",
      "Iteration No: 39 ended. Search finished for the next optimal point.\n",
      "Time taken: 691.7738\n",
      "Function value obtained: 0.5213\n",
      "Current minimum: 0.4419\n",
      "Iteration No: 40 started. Searching for the next optimal point.\n",
      "Iteration No: 40 ended. Search finished for the next optimal point.\n",
      "Time taken: 2973.6736\n",
      "Function value obtained: 0.4886\n",
      "Current minimum: 0.4419\n",
      "Iteration No: 41 started. Searching for the next optimal point.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration No: 41 ended. Search finished for the next optimal point.\n",
      "Time taken: 1120.7123\n",
      "Function value obtained: 0.4993\n",
      "Current minimum: 0.4419\n",
      "Iteration No: 42 started. Searching for the next optimal point.\n",
      "Iteration No: 42 ended. Search finished for the next optimal point.\n",
      "Time taken: 353.5213\n",
      "Function value obtained: 0.5569\n",
      "Current minimum: 0.4419\n",
      "Iteration No: 43 started. Searching for the next optimal point.\n",
      "Iteration No: 43 ended. Search finished for the next optimal point.\n",
      "Time taken: 2070.9181\n",
      "Function value obtained: 2.4857\n",
      "Current minimum: 0.4419\n",
      "Iteration No: 44 started. Searching for the next optimal point.\n",
      "Iteration No: 44 ended. Search finished for the next optimal point.\n",
      "Time taken: 1378.0682\n",
      "Function value obtained: 0.5624\n",
      "Current minimum: 0.4419\n",
      "Iteration No: 45 started. Searching for the next optimal point.\n",
      "Iteration No: 45 ended. Search finished for the next optimal point.\n",
      "Time taken: 811.8483\n",
      "Function value obtained: 1.1005\n",
      "Current minimum: 0.4419\n",
      "Iteration No: 46 started. Searching for the next optimal point.\n",
      "Iteration No: 46 ended. Search finished for the next optimal point.\n",
      "Time taken: 393.4375\n",
      "Function value obtained: 0.5848\n",
      "Current minimum: 0.4419\n",
      "Iteration No: 47 started. Searching for the next optimal point.\n",
      "Iteration No: 47 ended. Search finished for the next optimal point.\n",
      "Time taken: 772.5899\n",
      "Function value obtained: 0.9117\n",
      "Current minimum: 0.4419\n",
      "Iteration No: 48 started. Searching for the next optimal point.\n",
      "Iteration No: 48 ended. Search finished for the next optimal point.\n",
      "Time taken: 1118.7057\n",
      "Function value obtained: 2.4966\n",
      "Current minimum: 0.4419\n",
      "Iteration No: 49 started. Searching for the next optimal point.\n",
      "Iteration No: 49 ended. Search finished for the next optimal point.\n",
      "Time taken: 379.4622\n",
      "Function value obtained: 0.5513\n",
      "Current minimum: 0.4419\n",
      "Iteration No: 50 started. Searching for the next optimal point.\n",
      "Iteration No: 50 ended. Search finished for the next optimal point.\n",
      "Time taken: 2131.0316\n",
      "Function value obtained: 4.6486\n",
      "Current minimum: 0.4419\n",
      "73925.8363371 sec to optimize hyperparameters for XGBoost with cross-validation\n"
     ]
    }
   ],
   "source": [
    "LOG_LOSS = 'neg_log_loss'\n",
    "\n",
    "start_time = time.time()\n",
    "\n",
    "n_calls = 50\n",
    "# temporary redirection in case the kernel tends to interrupt and mess up the output\n",
    "# stdout = sys.stdout\n",
    "# with open('output.txt', 'w') as sys.stdout:\n",
    "#     res = skopt_gbrt_search(q_vec, labels, XGBClassifier, skopt_grid, LOG_LOSS, n_calls=n_calls)\n",
    "# sys.stdout = stdout\n",
    "\n",
    "res = skopt_gbrt_search(q_vec, labels, XGBClassifier, skopt_grid, LOG_LOSS, n_calls=n_calls)\n",
    "\n",
    "\n",
    "print(time.time()-start_time, \"sec to optimize hyperparameters for XGBoost with cross-validation\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "50 calls took 73,925.84 seconds $\\approx$ 20.5 hours to run.\n",
    "\n",
    "We _could_ make more calls and potentially get a better set of parameters. On AWS, each extra call increases the running time on by ~25 minutes. I decided against it for now. We can then take a look at the 5 best parameter sets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0.44190591057460116, {'n_estimators': 177, 'subsample': 0.9547615647987954, 'colsample_bytree': 0.14751315421666406, 'gamma': 0, 'objective': 'binary:logistic', 'learning_rate': 0.20374312390199137, 'max_depth': 12, 'min_child_weight': 1})\n",
      "(0.450003499736888, {'n_estimators': 173, 'subsample': 0.7724272099539228, 'colsample_bytree': 0.14301412377293404, 'gamma': 0, 'objective': 'binary:logistic', 'learning_rate': 0.20035580997038052, 'max_depth': 12, 'min_child_weight': 1})\n",
      "(0.4678065646592991, {'n_estimators': 172, 'subsample': 0.8083830575540177, 'colsample_bytree': 0.1466420600454964, 'gamma': 0, 'objective': 'binary:logistic', 'learning_rate': 0.3236444877922251, 'max_depth': 8, 'min_child_weight': 1})\n",
      "(0.46877660968947515, {'n_estimators': 175, 'subsample': 0.9154769880774636, 'colsample_bytree': 0.14608616503688332, 'gamma': 0, 'objective': 'binary:logistic', 'learning_rate': 0.3480749063995038, 'max_depth': 7, 'min_child_weight': 1})\n",
      "(0.46903095714082, {'n_estimators': 191, 'subsample': 0.9213832212741234, 'colsample_bytree': 0.12177194107091244, 'gamma': 0, 'objective': 'binary:logistic', 'learning_rate': 0.2972309168523623, 'max_depth': 12, 'min_child_weight': 1})\n"
     ]
    }
   ],
   "source": [
    "res_tuples = [(x['loss'], x['params']) for x in res]\n",
    "for r in sorted(res_tuples)[:5]:\n",
    "    print(r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now we can train a model based on the best parameters over the whole training data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n_estimators 177\n",
      "subsample 0.954761564799\n",
      "colsample_bytree 0.147513154217\n",
      "gamma 0\n",
      "objective binary:logistic\n",
      "learning_rate 0.203743123902\n",
      "max_depth 12\n",
      "min_child_weight 1\n",
      "238.549133062 sec to train XGBoost model\n"
     ]
    }
   ],
   "source": [
    "best_params = sorted(res_tuples)[0][1]\n",
    "for param in best_params.keys():\n",
    "    print(param, best_params[param])\n",
    "\n",
    "start_time = time.time()\n",
    "\n",
    "xgb_model = XGBClassifier(nthread=multiprocessing.cpu_count(), \n",
    "                          colsample_bytree = best_params['colsample_bytree'],\n",
    "                          learning_rate = best_params['learning_rate'],\n",
    "                          min_child_weight = best_params['min_child_weight'],\n",
    "                          n_estimators = best_params['n_estimators'],\n",
    "                          subsample = best_params['subsample'],\n",
    "                          objective = best_params['objective'],\n",
    "                          max_depth = best_params['max_depth'],\n",
    "                          gamma = best_params['gamma'])\n",
    "xgb_model.fit(q_vec, labels)\n",
    "\n",
    "print(time.time()-start_time, \"sec to train XGBoost model\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we generate and export predictions over the training data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "864.285174847 sec to generate 2345796 predictions\n"
     ]
    }
   ],
   "source": [
    "# make and report predictions scores for test data\n",
    "start_time = time.time()\n",
    "predictions = xgb_model.predict_proba(qtest_vec)\n",
    "print(time.time()-start_time, \"sec to generate\", predictions.shape[0], \"predictions\")\n",
    "\n",
    "predictions = pd.concat([pd.Series(np.arange(predictions.shape[0])), pd.Series(predictions[:,1])], axis=1)\n",
    "predictions.columns = [\"test_id\",\"is_duplicate\"]\n",
    "predictions.to_csv(\"xgboost_cv_predictions.csv\", index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results are the following:\n",
    "\n",
    "![Cross-validation-tuned XGBoost score: 0.53307](../notebooks-img/cv_xgboost_score_wide.png?raw=True \"Cross-validation-tuned XGBoost score: 0.53307\")\n",
    "\n",
    "This is a significant improvement over any of the previous classifiers! However, as mentioned before, it seems that even with optimized parameters and a powerful classifier such as XGBoost, Doc2Vec may not be producing the best vectorized sentences given our potentially-too-small training data.\n",
    "\n",
    "## Deep neural net (untuned)\n",
    "\n",
    "We can also try a deep neural net. Once again, this is adapted from a somewhat user-friendly ``DeepClassifier`` class using tensorflow by Roam's [Dingwall, Potts and Senaratna](https://roamanalytics.com/2016/09/13/prescription-based-prediction/#Deep-classifiers) [\\[github\\]](https://github.com/roaminsight/roamresearch/tree/master/BlogPosts/Prescription_based_prediction). You can see more details about the parameters and return values of the functions in ``DeepClassifier`` in [github](https://github.com/roaminsight/roamresearch/tree/master/BlogPosts/Prescription_based_prediction).\n",
    "\n",
    "For more nuts and bolts of implementing a net like this in tensorflow, check out [this repo](https://github.com/arturomp/udacity-deep-learning/blob/master/3_regularization.ipynb) or [this other one](https://github.com/arturomp/kaggle/blob/master/digit-recognition/cnn-digit-recognition.ipynb) where I use deep nets for optical character recognition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class DeepClassifier:\n",
    "    \"\"\"Defines a feed-forward neural network with two hidden layers.\n",
    "    Roughly,\n",
    "    h1 = f(xW1 + b1)\n",
    "    h2 = g(h1W2 + b2)\n",
    "    y = softmax(h2W3 + b3)\n",
    "    where drop-out is applied to h1 and h2.    \n",
    "    \"\"\"\n",
    "    def __init__(self,\n",
    "            hidden_dim1=200,\n",
    "            hidden_dim2=100,\n",
    "            activation1=tf.nn.relu,\n",
    "            activation2=tf.nn.relu,\n",
    "            keep_prob1=0.7,\n",
    "            keep_prob2=0.7,\n",
    "            eta=0.01,\n",
    "            max_iter=100,\n",
    "            tol=1e-05,\n",
    "            verbose=True):\n",
    "        \n",
    "        self.hidden_dim1 = hidden_dim1\n",
    "        self.hidden_dim2 = hidden_dim2\n",
    "        self.activation1 = activation1\n",
    "        self.activation2 = activation2\n",
    "        self.keep_prob1 = keep_prob1\n",
    "        self.keep_prob2 = keep_prob2\n",
    "        self.eta = eta\n",
    "        self.max_iter = max_iter\n",
    "        self.tol = tol\n",
    "        self.verbose = verbose\n",
    "        self.params = ('hidden_dim1', 'hidden_dim2',\n",
    "                       'activation1', 'activation2',\n",
    "                       'keep_prob1', 'keep_prob2',                       \n",
    "                       'eta', 'max_iter', 'tol',\n",
    "                       'verbose')\n",
    "        \n",
    "    def fit(self, X, y):\n",
    "        \"\"\"Specifies the model graph and performs training.\n",
    "        \"\"\"\n",
    "        # Set-up the dataset:\n",
    "        self.input_dim = X.shape[1]        \n",
    "        self.classes_ = sorted(set(y))\n",
    "        self.output_dim = len(self.classes_)\n",
    "        y_ = self.onehot_encode(y)        \n",
    "        # Begin the tf session:\n",
    "        tf.reset_default_graph()\n",
    "        self.sess = tf.InteractiveSession()        \n",
    "        # Model:\n",
    "        self._build_graph()        \n",
    "        # Optimization:\n",
    "        cost = tf.reduce_mean(\n",
    "            tf.nn.softmax_cross_entropy_with_logits(self.model, self.outputs))\n",
    "\n",
    "        # from previous deep nets\n",
    "        global_step = tf.Variable(0)  # count the number of steps taken.\n",
    "        decay_steps = 100\n",
    "        decay_rate = 0.9\n",
    "        learning_rate = tf.train.exponential_decay(self.eta, global_step, decay_steps, decay_rate)\n",
    "        \n",
    "        optimizer = tf.train.AdagradOptimizer(learning_rate).minimize(cost)\n",
    "        # Initialization:\n",
    "        init = tf.initialize_all_variables()\n",
    "        self.sess.run(init)                        \n",
    "        for i in range(1, self.max_iter+1):\n",
    "            # Training step:\n",
    "            _, loss = self.sess.run([optimizer, cost],\n",
    "                feed_dict={self.inputs: X,\n",
    "                           self.outputs: y_,\n",
    "                           self.keepprob_holder1: self.keep_prob1,\n",
    "                           self.keepprob_holder2: self.keep_prob2})\n",
    "            # Progress report:\n",
    "            self._progressbar(loss, i)\n",
    "            if loss <= self.tol:\n",
    "                sys.stderr.write('Stopping criteria reached.')\n",
    "        if self.verbose:\n",
    "            sys.stderr.write(\"\\n\")\n",
    "\n",
    "    def _build_graph(self):\n",
    "        \"\"\"Builds the core computation graph.\"\"\"\n",
    "        # Inputs and outputs:\n",
    "        self.inputs = tf.placeholder(tf.float32, [None, self.input_dim])\n",
    "        self.outputs = tf.placeholder(tf.float32, [None, self.output_dim])\n",
    "        # Layer 1:        \n",
    "        W1 = self._weight_init(self.input_dim, self.hidden_dim1, name='W1')\n",
    "        b1 = self._bias_init(self.hidden_dim1, name='b1')\n",
    "        hidden1 = self.activation1(tf.matmul(self.inputs, W1) + b1)\n",
    "        # Dropout 1:\n",
    "        self.keepprob_holder1 = self._dropout_init('keep_prob1')\n",
    "        dropout_layer1 = tf.nn.dropout(hidden1, self.keepprob_holder1)\n",
    "        # Layer 2:\n",
    "        W2 = self._weight_init(self.hidden_dim1, self.hidden_dim2, name='W2')\n",
    "        b2 = self._bias_init(self.hidden_dim2, name='b2')\n",
    "        hidden2 = self.activation2(tf.matmul(dropout_layer1, W2) + b2)\n",
    "        # Dropout 2:\n",
    "        self.keepprob_holder2 = self._dropout_init('keep_prob2')\n",
    "        dropout_layer2 = tf.nn.dropout(hidden2, self.keepprob_holder2)\n",
    "        # Output layer:\n",
    "        W3 = self._weight_init(self.hidden_dim2, self.output_dim, name='W3')\n",
    "        b3 = self._bias_init(self.output_dim, name='b3')\n",
    "        # No softmax here; that's handled by the cost function.\n",
    "        self.model = tf.matmul(dropout_layer2, W3) + b3\n",
    "\n",
    "    def predict(self, X, prob=False):\n",
    "        \"\"\"Predict method that mimics `sklearn` by accepting\n",
    "        an np.array and returning a vector of predictions.\n",
    "        \"\"\"\n",
    "        predictions = self.sess.run(self.model,\n",
    "            feed_dict={self.inputs: X,\n",
    "                       self.keepprob_holder1: 1.0,\n",
    "                       self.keepprob_holder2: 1.0})\n",
    "        # changed to output probabilities\n",
    "        if prob:\n",
    "            return tf.nn.softmax(predictions)\n",
    "        else:\n",
    "            return self._predictionvecs2class(predictions)                               \n",
    "    \n",
    "    def _weight_init(self, m, n, name):\n",
    "        \"\"\"Weight initialization according to the heuristic\n",
    "        that the values should be uniformly distributed around\n",
    "        \"\"\"\n",
    "        x = np.sqrt(6.0/(m+n))\n",
    "        with tf.name_scope(name) as scope: \n",
    "            return tf.Variable(\n",
    "                tf.random_uniform(\n",
    "                    [m, n], minval=-x, maxval=x), name=name)\n",
    "\n",
    "    def _bias_init(self, dim, name, constant=0.0):\n",
    "        \"\"\"Bias initialization, by default as all 0s.\n",
    "        \"\"\"\n",
    "        with tf.name_scope(name) as scope:            \n",
    "            return tf.Variable(\n",
    "                tf.constant(constant, shape=[dim]), name=name)\n",
    "\n",
    "    def _dropout_init(self, name):\n",
    "        \"\"\"Initialize a placeholder for a dropout value.\"\"\"        \n",
    "        with tf.name_scope(name) as scope:\n",
    "            return tf.placeholder(tf.float32, name=name)\n",
    "                \n",
    "    def onehot_encode(self, y, on_value=1.0):\n",
    "        \"\"\"Turns the list of class labels `y` into a matrix of\n",
    "        one-hot encoded vectors. This could be replaced by\n",
    "        `tf.one_hot`, but this native version does the job.\n",
    "        \"\"\"        \n",
    "        classmap = dict(zip(self.classes_, range(self.output_dim)))        \n",
    "        y_ = np.zeros((len(y), self.output_dim))\n",
    "        for i, cls in enumerate(y):\n",
    "            y_[i][classmap[cls]] = on_value            \n",
    "        return y_\n",
    "\n",
    "    def _progressbar(self, loss, index):\n",
    "        \"\"\"Overwriting progress bar for feedback on training process.\n",
    "        Prints to standard error.\n",
    "        \"\"\"        \n",
    "        if self.verbose:        \n",
    "            sys.stderr.write('\\r')\n",
    "            sys.stderr.write(\"Iteration {}: loss is {}\".format(index, loss))\n",
    "            sys.stderr.flush()\n",
    "\n",
    "    def get_params(self, deep=True):\n",
    "        \"\"\"Gets the hyperparameters for the model, as given by the\n",
    "        `self.params` attribute. This is called `get_params` for\n",
    "        compatibility with sklearn. `deep=True` is ignored, but is\n",
    "        needed for sklearn.\n",
    "        \"\"\"\n",
    "        return {p: getattr(self, p) for p in self.params}\n",
    "\n",
    "    def set_params(self, **params):\n",
    "        \"\"\"Use the params dict to set attribute values. This\n",
    "        is needed for sklearn `GridSearchCV` compatibility.\n",
    "        \"\"\"        \n",
    "        for key, val in six.iteritems(params):\n",
    "            setattr(self, key, val)\n",
    "        return self\n",
    "\n",
    "    def _predictionvecs2class(self, predictions):\n",
    "        \"\"\"\n",
    "        Convert the matrix of prediction probabilities into classes.\n",
    "        In cases of ties, a random choices is made to avoid spurious\n",
    "        patterns resulting from guessing classes that are earlier\n",
    "        in ordering in case of ties.              \n",
    "        \"\"\"\n",
    "        maxprobs = predictions.max(axis=1)\n",
    "        cats = []\n",
    "        for row, maxprob in zip(predictions, maxprobs):\n",
    "            i = np.random.choice([i for i, val in enumerate(row)\n",
    "                                  if val==maxprob])\n",
    "            cats.append(self.classes_[i])        \n",
    "        return cats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After defining the class, we train a model based on it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Iteration 3000: loss is 0.563651680946\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "38095.938349 sec to train deep net\n",
      "422.033621073 sec to generate 2345796 predictions\n"
     ]
    }
   ],
   "source": [
    "start_time = time.time()\n",
    "\n",
    "mod = DeepClassifier(\n",
    "           hidden_dim1=200,\n",
    "           hidden_dim2=100,\n",
    "           keep_prob1=0.7,\n",
    "           keep_prob2=0.7,\n",
    "           activation1=tf.nn.relu,\n",
    "           activation2=tf.nn.relu,\n",
    "           verbose=True,\n",
    "           max_iter=3000,\n",
    "           eta=0.01)\n",
    "\n",
    "mod.fit(q_vec, labels)\n",
    "\n",
    "print(time.time()-start_time, \"sec to train deep net\")\n",
    "\n",
    "start_time = time.time()\n",
    "# Test predictions and scoring:\n",
    "predictions = tf.nn.softmax(mod.predict(qtest_vec, prob=True)).eval()[:,1:]\n",
    "print(time.time()-start_time, \"sec to generate\", predictions.shape[0], \"predictions\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Training this deep net took 38095.94 seconds $\\approx$ 11.5 hours to run.\n",
    "\n",
    "Ideally, we'd like to cross-validate different parameters for our particular architecture and, even further, try out a couple of different architectures with their own cross-validations. However, that would be a bit too resource-intensive for a cursory view at a deep net's performance, as is evident from running only one of them, and as Roam's Dingwall, Potts and Senaratna themselves acknowledge.\n",
    "\n",
    "Now we export the predictions to csv and submit them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "deep_predictions = pd.concat([pd.Series(np.arange(predictions.shape[0])), pd.Series(predictions.ravel())], axis=1)\n",
    "deep_predictions.columns = [\"test_id\",\"is_duplicate\"]\n",
    "deep_predictions.to_csv(\"deep_predictions.csv\", index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is better than untuned XGBoost, but not better than just plain cosine similarity, and considerably worse than XGBoost tuned with scikit-optimize.\n",
    "\n",
    "![Deep (untuned) neural net score: 0.53307](../notebooks-img/deep_score_wide.png?raw=True \"Deep (untuned) neural net score: 0.53307\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Opportunities and conclusions\n",
    "\n",
    "Taken together, these results point toward an opportunity to better vectorize our sentences. Further experiments can easily incorporate TFIDF-vectorization. Another possibility is to use the test data, or even non-uniqued training data, to train the Doc2Vec model. At the same time, it may be the case that Doc2Vec is not the best way to vectorize documents. \n",
    "\n",
    "Once we have a vector representaion for sentences, we can achieve good results with cosine similarity, tuned XGBoost or optimized deep nets classifiers. Coupled with a better way of representing sentences, these classifiers should be able to improve the current best score on Kaggle data of 0.53307."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}