02_decoding_symbols/DecodingTutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<title>Decoding</title>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy.stats import multivariate_normal\n",
    "from tkinter import *\n",
    "import numpy as np\n",
    "import numpy.matlib\n",
    "import math\n",
    "%run performance_utils.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<script>\n",
       "code_show=true; \n",
       "function code_toggle() {\n",
       " if (code_show){\n",
       " $('div.input').hide();\n",
       " } else {\n",
       " $('div.input').show();\n",
       " }\n",
       " code_show = !code_show\n",
       "} \n",
       "$( document ).ready(code_toggle);\n",
       "</script>\n",
       "The raw code for this IPython notebook is by default hidden for easier reading.\n",
       "To toggle on/off the raw code, click <a href=\"javascript:code_toggle()\">here</a>."
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import HTML\n",
    "HTML('''<script>\n",
    "code_show=true; \n",
    "function code_toggle() {\n",
    " if (code_show){\n",
    " $('div.input').hide();\n",
    " } else {\n",
    " $('div.input').show();\n",
    " }\n",
    " code_show = !code_show\n",
    "} \n",
    "$( document ).ready(code_toggle);\n",
    "</script>\n",
    "The raw code for this IPython notebook is by default hidden for easier reading.\n",
    "To toggle on/off the raw code, click <a href=\"javascript:code_toggle()\">here</a>.''')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Probabilistic Decoding for Intelligent Text Entry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3><em>CHI 2019 Course on Computational Interaction</em></h3>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"http://pokristensson.com\">Per Ola Kristensson</a>, Department of Engineering, University of Cambridge, United Kingdom.\n",
    "\n",
    "Copyright (c) 2017-2019 Per Ola Kristensson. All rights reserved."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An intelligent text entry method is a text entry method that _infers_ users' intended text. We have known how to design such methods for quite some time. One of the more interesting examples is a 12th century shorthand system called _Nova ars notaria_ (\"the new note art\") by the monk of John of Tilbury. What is interesting with this system is that the _design principles_ are known:\n",
    "\n",
    "* Simplify letters to line marks.\n",
    "* Compress common word stems into sequences of simple line marks and dots.\n",
    "* Identify common word stems by frequency analysis\n",
    "\n",
    "In other words, two fundamental principles underpin the design of any efficient text entry method:\n",
    "\n",
    "* Minimise users' effort in articulating their input.\n",
    "* Exploit the redundancies in natural languages using statistical language modelling.\n",
    "\n",
    "However, _how_ can we both minimise users' articulation effort and leverage the redundancies in natural languages? The solution is to perform _statistical decoding_. A statistical decoder is a generative probabilistic model that is capable of searching a vast hypothesis space in order to identify users' intended text given noisy observations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Statistical Decoding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The problem we are adressing is a user attempting to communciate information over some form of channel. In Human-Computer Interaction (HCI) we typically model a user transmitting _information_ to a computer system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The term _information_ in HCI usually refers to characters (for example, typing on a keyboard), words (using speech recognition), or commands (for instance by using keyboard shortcuts or touchscreen gestures)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More formally, we intend to transmit a message $y$ via some form of signal $x$. In a perfect world this would be trivially achieved via a lookup-table. Unfortunately we live in an imperfect world and as a consequence our signal will be perturbed by noise in our neuromuscular system, device sensor imprecision, cognitive errors by the user, etc. Due to this inherent uncertainty it makes sense to model the problem _probabilistically_."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then wish to compute the probability of the message $y$ _given_ the signal $x$. This can be written mathematically as $P(y | x)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>\n",
    "The probability $P(y|x)$ is known as a _conditional probability_ and it is (by either the definition of conditional probability or as an axiom of probability):\n",
    "</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation}\n",
    "P(y|x) = \\frac{P(x \\cap y)}{P(x)}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(We are going to assume $P(x)\\neq0$ and $P(y)\\neq 0$).\n",
    "\n",
    "The above equation states that the conditional probability of $y$ given $x$ is identical to the ratio of the _joint probability_ of $x$ and $y$ and the probability of $x$. The joint probability of $x$ and $y$ can also be written as $P(x,y)$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can rewrite the conditional probability $P(y|x)$ as follows:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "P(y|x) &= \\frac{P(x \\cap y)}{P(x)}\\\\\n",
    "P(x|y) &= \\frac{P(y \\cap x)}{P(y)} = \\frac{P(x \\cap y)}{P(y)}\\\\\n",
    "\\Rightarrow P(x \\cap y) &= P(x|y)P(y) = P(y|x)P(x)\\\\\n",
    "\\Rightarrow P(y|x) &= \\frac{P(x|y)P(y)}{P(x)}\\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This last expression is known as _Bayes' rule_ (or theorem). Usually we have many possible messages that we wish to decode and $P(y|x)$ will then become the _posterior_ probability distribution, assigning a probability to every possible message. Our objective is to compute this posterior probability distribution and select the most probable message."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, since we are usually only interested in the most probable message $\\hat{y}$ we can write:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation}\n",
    "\\hat{y}=\\underset{y}{\\arg\\max}\\left[P(y|x)\\right]\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have already seen that the conditional probability $P(y|x)$ can be written using Bayes' rule:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation}\n",
    "\\hat{y}=\\underset{y}{\\arg\\max}\\left[\\frac{P(x|y)P(y)}{P(x)}\\right]\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, as we are only interested in the message that maximises the conditional probability of the message given the signal, the denominator $P(x)$ will be invariant in the search and can therefore be dropped:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation}\n",
    "\\hat{y}=\\underset{y}{\\arg\\max}\\left[P(x|y)P(y)\\right]\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$P(x|y)$ is the likelihood of the signal $x$ given a particular hypothesis for what the message $y$ could be. $P(y)$ is the _prior_ probability of the message, that is, without taking any signal into account. For instance, if a system can only recognise two messages, $x_1$ and $x_2$, and both are equally likely in the absence of any additional information, then the prior probability of either $x_1$ or $x_2$ is $0.5$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Identifying the highest probable message $y$ is a _search problem_. We search by consulting a model of the likelihood of a signal $x$ given a message $y$ under consideration and by consulting a model of the prior probability of a message $y$ without any consideration to any signal. This search will generate _hypotheses_ and these hypotheses will have probabilities assigned to them. Usually, the hypothesis with the highest probability assigned to it is our preferred hypothesis:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation}\n",
    "\\hat{\\text{hypothesis}}=\\underset{\\text{hypotheses}}{\\arg\\max}\\left(\\text{likelihood model}\\cdot\\text{prior model}\\right)\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Sequence Decoding</h2>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us assume we have a simple touchscreen keyboard which can generate $k$ distinct letters $L=\\left\\{l_1, l_2, \\ldots, l_k\\right\\}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We view our signal as a discrete <em>observation sequence</em>, which, for ease of notation, we will denote $O$. When a user is typing on the touchscreen keyboard we are provided with a touch point coordinate $(x,y)$, where $x$ and $y$ are the horizontal and vertical components of the touch point respectively. A set of $n$ touch points then form $n$ observations:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation}\n",
    "O = \\left\\{o^{(1)}, o^{(2)}, \\ldots, o^{(n)}\\right\\}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since subscripts are often used to denote specific members of a set, we use a superscript within parantheses to denote the index of an individual observation in an sequence (indexing starts at 1). This means that for instance the sequence $l_1, l_2$ means two letters selected from the set of letters that can be generated by the keyboard (which, since we have defined the keyboard above to generate distinct letters, means $l_1 \\neq l_2$). In contrast, $l_1^{(1)},l_1^{(2)}$ means we have generated the _same_ letter $l_1$ twice. Similarly, $l_1^{(1)},l_2^{(2)}$ means we have generated the letter sequence $l_1, l_2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are now given a set of observations $O = \\left\\{o^{(1)}, o^{(2)}, \\ldots, o^{(n)}\\right\\}$ and we wish to find a _hypothesis_ for $O$, which is a corresponding letter sequence ${L} = \\left\\{l^{(1)}, l^{(2)}, \\ldots, l^{(m)}\\right\\}$. Note that the cardinality of $O$ and $L$ are not necessarily identical, that is, it is not neccesarily true that $n = m$ ($|O|=|L|$). There are obviously many possible letter sequences that can match the observation sequence. However, we are only interested in the most _probable_ letter sequence $\\hat{L}$, and this can be computed using the expression we derived earlier:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation}\n",
    "\\hat{L}=\\underset{L}{\\arg\\max}\\left[P(O|L)P(L)\\right]\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$P(L)$ is the prior. What would be a suitable model for the prior probability of a letter sequence? This particular problem is known as _statistical language modelling_ and it is a research area that has received considerable attention. A very simple model is the unigram probability model, which in this case simply assigns a probability to each individual letter $l_i$. The probabilities assigned over the letters $l_1, l_2, \\ldots, l_k$ (remember that our keyboard could generate $k$ distinct letters) then must sum to one and this is our prior probability distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$P(O|L)$ is the likelihood of the observation sequence given a particular letter sequence hypothesis. For a touchscreen keyboard this means that we need a model that assigns a probability to an individual $(x,y)$ touch point given a particular letter key hypothesis, that is $P(o_i|l_i)$. A simple way to achieve this is to centre a two-dimensional Gaussian distribution at each key centre and assume each 2D Gaussian has equal variance and that there is no correlation between the horizontal and vertical axes (that is, the covariance matrix is diagonal). \n",
    "\n",
    "<img src=\"imgs/key_with_2d_gaussian.svg\">\n",
    "\n",
    "To compute our likelihood distribution we need to compute the probability of observing a specific touch point for all possible $k$ letters that can be generated by the keyboard. For example, if our keyboard has two keys, say $A$ and $B$ generating the letters $a$ and $b$ respectively, then our likelihood distribution would need to be computed by computing the probability of the touch point for key $A$ and key $B$. If the touch point is centred at key $A$ the probability will be maximised for $A$ and smaller for $B$ (we assume keys do not overlap spatially).\n",
    "\n",
    "<img src=\"imgs/ab_keys_with_2d_gaussian.svg\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>\n",
    "The search for the most probable letter sequence given an observation sequence can be modelled as a (discrete-time) <em>Markov chain</em>. A Markov chain assumes a future state is conditional on the previously known state and independent of any other previous or future states.\n",
    "</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Statistical Decoding using Token Passing</h2>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>\n",
    "We will now explain how to statistically decode the observation sequence $O$ into the most probable letter sequence $L$. There are several algorithms for learning a hidden Markov model and performing inference. For our particular decoding problem a flexible and efficient method is known as <em>token passing</em>. A particular advantage is that it can efficiently search more complex models and it can be easily parallelised.\n",
    "</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>Token</h3>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>\n",
    "Fundamental to token passing is the notion of a <em>token</em>. A token is a data structure that at a minimum contains the following information: 1) the hypothesis generated so far; and 2) the accumulated probability of the hypothesis. For a touchscreen keyboard the hypothesis generated so far is a letter sequence (a string). For bookkeeping reasons we will also let each token store which observation index it belongs to.\n",
    "</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Token:\n",
    "    def __init__(self, hypo, acc_prob, obs_index):\n",
    "        self.hypo = hypo\n",
    "        self.acc_prob = acc_prob\n",
    "        self.obs_index = obs_index"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>Observation</h3>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>\n",
    "An observation is simply an $(x,y)$ touch point.\n",
    "</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Observation:\n",
    "    def __init__(self, x, y):\n",
    "        self.x = x\n",
    "        self.y = y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>Observation Sequence</h3>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>\n",
    "The observation sequence $O = \\left\\{o^{(1)}, o^{(2)}, \\ldots, o^{(n)}\\right\\}$ is simply a series of observations. They could be kept in a list but it is often convenient to represent them as a separate class when we later improve the model's performance.\n",
    "</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class ObservationSequence:\n",
    "    def __init__(self, seq):\n",
    "        self.seq = seq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>Key</h3>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>\n",
    "A key describes the geometry of the keyboard. In this case our likelihood model is very simple. All we need is information on where to centre the 2D Gaussian, which means we need the centre coordinate of the key, and the letter generated by the key. Below we have generated a simple keyboard layout with four keys: A, B, C and D.\n",
    "</p>\n",
    "\n",
    "<img src=\"imgs/key_layout.svg\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Key:\n",
    "    def __init__(self, centre_x, centre_y, letter):\n",
    "        self.centre_x = centre_x\n",
    "        self.centre_y = centre_y\n",
    "        self.letter = letter\n",
    "                \n",
    "a_key = Key(0.0,0.0,'a')\n",
    "b_key = Key(0.0,1.0,'b')\n",
    "c_key = Key(1.0,0.0,'c')\n",
    "d_key = Key(1.0,1.0,'d')\n",
    "keys = [a_key, b_key, c_key, d_key]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prior"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The prior is in this instance the simple unigram probability of the letter. It can be implemented using a simple hash table. To make the demonstration simple we will assume a unigram letter model and all of the letters will be equally probable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "lm = {'a':0.25,'b':0.25,'c':0.25,'d':0.25}\n",
    "def get_prior(context, lm):\n",
    "    c = context[-1]\n",
    "    return lm[c]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Likelihood"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The likelihood is a 2D Gaussian probability density function for the key corresponding to the letter under consideration, evaluated at the coordinate specified by the observation. To make the demonstration simple we will assume a symmetrical 2D Gaussian with unit variance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_likelihood(x, y, key):\n",
    "    key_x = key.centre_x\n",
    "    key_y = key.centre_y\n",
    "    p = multivariate_normal.pdf([x,y], mean=[key_x,key_y], cov=[[1,0],[0,1]]);\n",
    "    return p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Substitution-Only Decoder"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first decoder we will look at is a _substitution-only_ decoder. It takes as input:\n",
    "\n",
    "* A set of $k$ symbols that can be recognised.\n",
    "* A set of keys (a centre coordinate and a letter, the letter must correspond to a symbol in the symbol set).\n",
    "* A language model (in this case a simple mapping of a single letter to a probability).\n",
    "* An observation sequence."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our decoder will receive a series of touch points, that is $O = \\left\\{o^{(1)}, o^{(2)}, \\ldots, o^{(n)}\\right\\}$.\n",
    "\n",
    "To be able to decode such a model we start with a single initial token. The process now works as follows. We take any token at observation index $i \\in [0,n)$ and we propagate $k$ tokens to observation index $i+1$. We do this until we are at the last observation and can no longer propagate any tokens. These final tokens represent all hypotheses for the observation sequence. The token with the most probable hypothesis is the most likely letter sequence corresponding to the observation sequence."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model of the substitution-only decoder is shown below. An initial start state leads to a state each for every observation. States are linked by edges. Each edge links two states and generates a symbol. A series of edges generate a _path_. As is evident in the figure, several paths are possible from the start state to the end state. The subsitution-only decoder below will search _all_ such possible paths and substitute an observation for a symbol every time a particular edge is traversed when transitioning from one state to another.\n",
    "\n",
    "<img src=\"imgs/simple_chain.svg\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code below provides several possible observation sequences. The default observation sequence is empty and generates the expected result that the decoder predicts an empty hypothesis (no letter generated) given an empty observation sequence. By commenting and uncommenting the observation sequences it is possible to explore a few observation sequences:\n",
    "\n",
    "* No observations. This generates a single-state search space.\n",
    "* A noisy single touch in the vicinity of the A key. This generates a two-state search space.\n",
    "* Two noisy touches, both in the vicinity of the A key. This generates a three-state search space.\n",
    "* Two noisy touches near the B key followed by the A key. This generates a three-state search space.\n",
    "* Three noisy touches near the B key followed by the A key, followed by the D key. This generates a four-state search space (illustrated in the figure above)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "<strong>REFLECT:</strong> \n",
    "<ul>\n",
    "    <li>\n",
    "        How does the number of tokens relate to the number of observations?\n",
    "    </li>\n",
    "    <li>\n",
    "        Why is the letter sequence probability the same for 'b' and 'c' in the sequence with a single observation?\n",
    "    </li>\n",
    "</ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "letter sequence\n",
      "nth best, hypothesis, probability\n",
      "1, bad, 0.0000612\n",
      "2, bcd, 0.0000431\n",
      "3, dad, 0.0000410\n",
      "4, bac, 0.0000390\n",
      "5, bab, 0.0000353\n",
      "6, bbd, 0.0000336\n",
      "7, aad, 0.0000336\n",
      "8, dcd, 0.0000289\n",
      "9, bcc, 0.0000275\n",
      "10, dac, 0.0000262\n",
      "11, bcb, 0.0000249\n",
      "12, bdd, 0.0000237\n",
      "13, acd, 0.0000237\n",
      "14, dab, 0.0000237\n",
      "15, baa, 0.0000225\n",
      "...\n",
      "50, cdd, 0.0000087\n",
      "51, caa, 0.0000083\n",
      "52, dba, 0.0000083\n",
      "53, adc, 0.0000083\n",
      "54, cbc, 0.0000079\n",
      "55, adb, 0.0000075\n",
      "56, cbb, 0.0000071\n",
      "57, aba, 0.0000068\n",
      "58, cca, 0.0000058\n",
      "59, dda, 0.0000058\n",
      "60, cdc, 0.0000056\n",
      "61, cdb, 0.0000050\n",
      "62, ada, 0.0000048\n",
      "63, cba, 0.0000045\n",
      "64, cda, 0.0000032\n"
     ]
    }
   ],
   "source": [
    "symbols = ['a','b','c','d']\n",
    "#obs_seq = ObservationSequence([]) #No observations\n",
    "#obs_seq = ObservationSequence([Observation(0.1,0.1)]) #a + noise\n",
    "#obs_seq = ObservationSequence([Observation(0.1,0.1),Observation(0.15,-0.1)]) #aa + noise\n",
    "#obs_seq = ObservationSequence([Observation(0.1,1.1),Observation(0.15,-0.1)]) #ba + noise\n",
    "obs_seq = ObservationSequence([Observation(0.1,1.1),Observation(0.15,-0.1),Observation(1.05,0.95)]) #bad + noise\n",
    "completed_tokens = []\n",
    "\n",
    "def get_key_for_symbol(symbol):\n",
    "    for key in keys:\n",
    "        if key.letter == symbol:\n",
    "            return key\n",
    "    return None\n",
    "\n",
    "# Convenience function for nicely printing token passing results\n",
    "def print_sorted_tokens(sorted_tokens):\n",
    "    results = sorted_tokens\n",
    "    results = list(zip(range(1,len(results) + 1), results))\n",
    "    shortened_id = 15\n",
    "    shortened = False\n",
    "    if (len(results) > 32):\n",
    "        results = results[:shortened_id] + results[-shortened_id:]\n",
    "        shortened = True\n",
    "    print(\"letter sequence\")\n",
    "    print(\"nth best, hypothesis, probability\")\n",
    "    for i, result in results:\n",
    "        if (result.acc_prob < 1e-7):\n",
    "            print(\"%d, %s, %.3g\" % (i, result.hypo, result.acc_prob))\n",
    "        else:\n",
    "            print(\"%d, %s, %.7f\" % (i, result.hypo, result.acc_prob))\n",
    "        if (i == shortened_id and shortened):\n",
    "            print(\"...\")\n",
    "\n",
    "def propagate(token, obs_seq, completed_tokens):\n",
    "    # Which observation is the token in?\n",
    "    ix = token.obs_index\n",
    "    # Are there more observations?\n",
    "    next_ix = ix + 1\n",
    "    if next_ix == len(obs_seq.seq):\n",
    "        # No more observations, add token to the list of completed tokens\n",
    "        completed_tokens.append(token)\n",
    "    else:\n",
    "        # There is another observation, propagage a token per symbol to the next observation\n",
    "        next_obs = obs_seq.seq[next_ix]\n",
    "        for symbol in symbols:\n",
    "            # Find the key that generates the symbol\n",
    "            key = get_key_for_symbol(symbol)\n",
    "            # Calculate the prior probability of our symbol given previous symbols we have generated\n",
    "            prior = get_prior(token.hypo+symbol, lm)\n",
    "            # Calculate the likelihood of the observation (touch coordinate) given our current symbol hypothesis\n",
    "            likelihood = get_likelihood(next_obs.x, next_obs.y, key)\n",
    "            # Calculate an accumulated probability for our posterior probability\n",
    "            acc_prob = token.acc_prob * prior * likelihood\n",
    "            # Create a new token for the next observation index that contains the accumulated hypothesis and probability\n",
    "            new_token = Token(token.hypo+symbol, acc_prob, next_ix)\n",
    "            # Propagate our token\n",
    "            propagate(new_token, obs_seq, completed_tokens)\n",
    "\n",
    "# To start the process we need to create an initial token which we can then propagate to the first observation index\n",
    "seed_token = Token(\"\", 1.0, -1)\n",
    "# Propagating the first token starts the search process\n",
    "propagate(seed_token, obs_seq, completed_tokens)\n",
    "# Sort the list of completed tokens according to the probabilities of their hypotheses\n",
    "results = sorted(completed_tokens, key=lambda token: token.acc_prob, reverse=True)\n",
    "\n",
    "# Print out sequence hypotheses and their probabilities\n",
    "print_sorted_tokens(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>\n",
    "<div class=\"alert alert-block alert-success\"> \n",
    "<h3>Exercise 1</h3> \n",
    "<br>\n",
    "Write a test to independently compute the probabilities of the length 1 and 2 observation sequences examined above by completing the followings stubs.  Ensure you understand what the recursive token propagation method is achieving.\n",
    "</div>\n",
    "</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\"> \n",
    "<strong>Task 1:</strong> Complete the following stub code to ensure the test passes for the single observation case\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def TestSingleObservation(observation_sequence, sequence_probabilities, tolerance):\n",
    "    # This flag will stay True if your test passes\n",
    "    test_passes = True\n",
    "\n",
    "    # Pull out the single observation\n",
    "    obs1 = observation_sequence.seq[0]\n",
    "    \n",
    "    token_count = 0\n",
    "    \n",
    "    # Loop through possible symbols (letter keys)\n",
    "    for symbol in symbols:\n",
    "        # Find the key that generates the symbol\n",
    "        key = get_key_for_symbol(symbol)\n",
    "\n",
    "        # TODO compute the prior, likelihood and accumulated probability after the single observation\n",
    "        prior = get_prior(key.letter, lm)\n",
    "        likelihood = get_likelihood(obs1.x, obs1.y, key)\n",
    "        acc_prob = prior * likelihood\n",
    "\n",
    "        # Compare accumulated probability to that reported from our recursive token passing method\n",
    "        if not (abs(obs_seq1_prob[token_count] - acc_prob) < tol):\n",
    "            test_passes = False\n",
    "        token_count += 1\n",
    "        \n",
    "    return test_passes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the test to ensure it passes (no code changes required)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PASSED: well done!\n"
     ]
    }
   ],
   "source": [
    "# Do not change the tolerance\n",
    "tol = 0.000001\n",
    "\n",
    "# A sequence with a single Observation\n",
    "obs_seq1 = ObservationSequence([Observation(0.1,0.1)]) #a + noise\n",
    "# The probabilities associated with the four possible token hypotheses, i.e. 'a','b','c','d'\n",
    "obs_seq1_prob = [0.039393,0.026406,0.026406,0.017700]\n",
    "\n",
    "if (TestSingleObservation(obs_seq1,obs_seq1_prob,tol)):\n",
    "    print(\"PASSED: well done!\")\n",
    "else:\n",
    "    print(\"FAILED: try again!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\"> \n",
    "<strong>Task 2:</strong> Complete the following stub code to ensure the test passes for two observations in the sequence\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def TestTwoObservations(observation_sequence, sequence_probabilities, tolerance):\n",
    "    # This flag will stay True if your test passes\n",
    "    test_passes = True\n",
    "\n",
    "    # Pull out the two observations in the sequence\n",
    "    obs1 = observation_sequence.seq[0]\n",
    "    obs2 = observation_sequence.seq[1]\n",
    "\n",
    "    token_count = 0\n",
    "    \n",
    "    # Loop through possible symbols (letter keys) for the first observation\n",
    "    for symbol_obs1 in symbols:\n",
    "        # Find the key that generates the symbol\n",
    "        key1 = get_key_for_symbol(symbol_obs1)\n",
    "\n",
    "        # TODO compute the prior, likelihood and accumulated probability after the first observation\n",
    "        obs1_prior = get_prior(key1.letter, lm)\n",
    "        obs1_likelihood = get_likelihood(obs1.x, obs1.y, key1)\n",
    "        acc_prob_obs1 = obs1_prior * obs1_likelihood\n",
    "\n",
    "        # Loop through possible symbols (letter keys) for the second observation\n",
    "        for symbol_obs2 in symbols:\n",
    "            # Find the key that generates the symbol\n",
    "            key2 = get_key_for_symbol(symbol_obs2)\n",
    "            \n",
    "            # TODO compute the prior, likelihood and accumulated probability after the second observation\n",
    "            obs2_prior = get_prior(key2.letter, lm)\n",
    "            obs2_likelihood = get_likelihood(obs2.x, obs2.y, key2)\n",
    "            acc_prob = acc_prob_obs1 * obs2_prior * obs2_likelihood\n",
    "\n",
    "            # Compare accumulated probability to that reported from our recursive token passing method\n",
    "            if not (abs(obs_seq2_prob[token_count] - acc_prob) < tol):           \n",
    "                test_passes = False\n",
    "            token_count += 1\n",
    "    return test_passes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the test for the case of two observations in a sequence (no code changes required)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PASSED: well done!\n"
     ]
    }
   ],
   "source": [
    "# Do not change the tolerance\n",
    "tol = 0.000001\n",
    "\n",
    "# A sequence with a two Observations\n",
    "obs_seq2 = ObservationSequence([Observation(0.1,0.1),Observation(0.15,-0.1)]) #aa + noise\n",
    "# The probabilities associated with the 16 possible token hypotheses, i.e. 'aa','ab','ac','ad', 'ba', 'bb', ... ,'dd'\n",
    "obs_seq2_prob = [0.001542,0.000846,0.001087,0.000596,\\\n",
    "                 0.001034,0.000567,0.000728,0.000400,\\\n",
    "                 0.001034,0.000567,0.000728,0.000400,\\\n",
    "                 0.000693,0.000380,0.000488,0.000268,]\n",
    "\n",
    "if (TestTwoObservations(obs_seq2,obs_seq2_prob,tol)):\n",
    "    print(\"PASSED: well done!\")\n",
    "else:\n",
    "    print(\"FAILED: try again!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Decoding with Substitutions and Deletions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our first decoder has some limitations. In particular, it cannot identify a hypothesis which contains more symbols than the number of observations. However, we can easily modify the decoder to handle both substitutions and deletions of observations. We do this by modifying our <em>propagate</em> function so that instead of generating $k$ transitions ($k$ tokens), each transition modelling a substituion of an observation with a letter, we generate $k+1$ transitions ($k+1$ tokens). The additional transition (token) models an $\\epsilon$-transition (see figure below).\n",
    "\n",
    "<img src=\"imgs/simple_epsilon_chain.svg\">\n",
    "\n",
    "Incorporating the $\\epsilon$-transition, we can see below that the decoder is able to ignore noisy observations. For instance, for the observation sequence modelling letter keys: B, A, (noise), D, the decoder is able to delete the noisy observation. If we had a more sophisticated language model, the decoder would now be able to correct misspellings and accidental touches on letter keys."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "<strong>REFLECT:</strong> \n",
    "<ul>\n",
    "    <li>\n",
    "        Compare the deletion penalty value used here to the probabilities determined for the single observation sequence examined above.\n",
    "    </li>\n",
    "    <li>\n",
    "        Roughly, at what distance does the spurious third observation in the four observation sequence below start to be interpreted as a legitimate touch for the best hypothesis?\n",
    "    </li>\n",
    "</ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "letter sequence\n",
      "nth best, hypothesis, probability\n",
      "1, bad, 0.0000012\n",
      "2, bcd, 0.0000009\n",
      "3, dad, 0.0000008\n",
      "4, bac, 0.0000008\n",
      "5, bab, 0.0000007\n",
      "6, bbd, 0.0000007\n",
      "7, aad, 0.0000007\n",
      "8, bd, 0.0000006\n",
      "9, ad, 0.0000006\n",
      "10, ba, 0.0000006\n",
      "11, dcd, 0.0000006\n",
      "12, bcc, 0.0000006\n",
      "13, dac, 0.0000005\n",
      "14, bcb, 0.0000005\n",
      "15, bdd, 0.0000005\n",
      "...\n",
      "611, bba, 2.95e-12\n",
      "612, cbb, 2.93e-12\n",
      "613, db, 2.85e-12\n",
      "614, ba, 2.75e-12\n",
      "615, ccba, 2.75e-12\n",
      "616, ddba, 2.75e-12\n",
      "617, cb, 2.73e-12\n",
      "618, cdbc, 2.61e-12\n",
      "619, cdbb, 2.36e-12\n",
      "620, adba, 2.25e-12\n",
      "621, cbba, 2.14e-12\n",
      "622, dba, 2.08e-12\n",
      "623, cdb, 2.07e-12\n",
      "624, cba, 1.99e-12\n",
      "625, cdba, 1.51e-12\n"
     ]
    }
   ],
   "source": [
    "symbols = ['a','b','c','d','ɛ'] #Epsilon represents a delete observation event (i.e. do not generate a hypothesis)\n",
    "#obs_seq = ObservationSequence([]) #No observations\n",
    "#obs_seq = ObservationSequence([Observation(0.1,1.1),Observation(0.15,-0.1),Observation(1.05,0.95)]) #bad + noise\n",
    "obs_seq = ObservationSequence([Observation(0.1,1.1),Observation(0.15,-0.1),Observation(3.0,-2.7),Observation(1.05,0.95)]) #b - a - (Noisy input far away from any key) - d + noise\n",
    "completed_tokens = []\n",
    "deletion_penalty = 0.02\n",
    "\n",
    "def get_key_for_symbol(symbol):\n",
    "    for key in keys:\n",
    "        if key.letter == symbol:\n",
    "            return key\n",
    "    return None\n",
    "\n",
    "def propagate(token, obs_seq, completed_tokens):\n",
    "    # Which observation is the token in?\n",
    "    ix = token.obs_index\n",
    "    # Are there more observations?\n",
    "    next_ix = ix + 1\n",
    "    if next_ix == len(obs_seq.seq):\n",
    "        # No more observations, add token to the list of completed tokens\n",
    "        completed_tokens.append(token)\n",
    "    else:\n",
    "        # There is another observation, propagage a token per symbol to the next observation\n",
    "        next_obs = obs_seq.seq[next_ix]\n",
    "        for symbol in symbols:\n",
    "            if symbol == 'ɛ':\n",
    "                # Delete observation\n",
    "                acc_prob = token.acc_prob * deletion_penalty\n",
    "                # Create a new token for the next observation index but don't accumulate a hypothesis\n",
    "                new_token = Token(token.hypo, acc_prob, next_ix)\n",
    "                propagate(new_token, obs_seq, completed_tokens)\n",
    "            else:\n",
    "                # Substitute observation for symbol\n",
    "                key = get_key_for_symbol(symbol)\n",
    "                prior = get_prior(token.hypo+symbol, lm)\n",
    "                likelihood = get_likelihood(next_obs.x, next_obs.y, key)\n",
    "                acc_prob = token.acc_prob * prior * likelihood\n",
    "                new_token = Token(token.hypo+symbol, acc_prob, next_ix)\n",
    "                propagate(new_token, obs_seq, completed_tokens)\n",
    "\n",
    "seed_token = Token(\"\", 1.0, -1)\n",
    "propagate(seed_token, obs_seq, completed_tokens)\n",
    "results = sorted(completed_tokens, key=lambda token: token.acc_prob, reverse=True)\n",
    "\n",
    "# Print out sequence hypotheses and their probabilities\n",
    "print_sorted_tokens(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Beam Pruning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our prior decoders will generate exponential growth in the number of tokens as a function of the number of observations. At observation index 1 we have propagated $k$ tokens, at observation index 2 we have propaged $k^2$ tokens, and so on.\n",
    "\n",
    "To keep the search complexity tractable we observe that many tokens that we are propagated are highly unlikely to contain the most probable hypotheses. These tokens can be filtered out using <em>beam pruning</em>. For every observation index we store the highest accumulated probability we have generated so far for that particular observation. We then impose a condition before we propagate a token: only propagate the token if the token's accumulated probability is within a certain distance to the highest accumulated probability for the next observation. This distance is known as the <em>beam width</em> and it is a parameter that can be used to set the operating point (or trade-off) between accuracy and time when searching for a hypothesis. The narrower the beam the more greedy the search becomes.\n",
    "\n",
    "We can see that the decoder generates less hypotheses as some hypotheses are eliminated in the search."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "<strong>REFLECT:</strong> \n",
    "<ul>\n",
    "    <li>\n",
    "        Examine how increasing the beam width to around 0.7 affects the resulting hypotheses.\n",
    "    </li>\n",
    "    <li>\n",
    "        What is a potential downside of introducing beam pruning?\n",
    "    </li>\n",
    "</ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "letter sequence\n",
      "nth best, hypothesis, probability\n",
      "1, bad, 0.0000012\n",
      "2, bcd, 0.0000009\n",
      "3, dad, 0.0000008\n",
      "4, bac, 0.0000008\n",
      "5, bab, 0.0000007\n",
      "6, bbd, 0.0000007\n",
      "7, aad, 0.0000007\n",
      "8, bd, 0.0000006\n",
      "9, ad, 0.0000006\n",
      "10, ba, 0.0000006\n",
      "11, dcd, 0.0000006\n",
      "12, bcc, 0.0000006\n",
      "13, dac, 0.0000005\n",
      "14, bcb, 0.0000005\n",
      "15, bdd, 0.0000005\n",
      "...\n",
      "69, aca, 0.0000002\n",
      "70, a, 0.0000002\n",
      "71, adc, 0.0000002\n",
      "72, adb, 0.0000001\n",
      "73, aba, 0.0000001\n",
      "74, aacd, 4.72e-09\n",
      "75, aacc, 3.01e-09\n",
      "76, aacb, 2.73e-09\n",
      "77, aac, 2.38e-09\n",
      "78, aaca, 1.74e-09\n",
      "79, aaad, 3.88e-10\n",
      "80, aaac, 2.47e-10\n",
      "81, aaab, 2.24e-10\n",
      "82, aaa, 1.95e-10\n",
      "83, aaaa, 1.43e-10\n"
     ]
    }
   ],
   "source": [
    "class ObservationSequence:\n",
    "    def __init__(self, seq):\n",
    "        self.seq = seq\n",
    "        self.beam = [0.0]*len(seq) \n",
    "\n",
    "symbols = ['a','b','c','d','ɛ'] #Epsilon represents a delete observation event (i.e. do not generate a hypothesis)\n",
    "#obs_seq = ObservationSequence([]) #No observations\n",
    "#obs_seq = ObservationSequence([Observation(0.1,1.1),Observation(0.15,-0.1),Observation(1.05,0.95)]) #bad + noise\n",
    "obs_seq = ObservationSequence([Observation(0.1,1.1),Observation(0.15,-0.1),Observation(3.0,-2.7),Observation(1.05,0.95)]) #b - a - (Noisy input far away from any key) - d + noise\n",
    "completed_tokens = []\n",
    "deletion_penalty = 0.02\n",
    "beam_width = 0.2\n",
    "\n",
    "def get_key_for_symbol(symbol):\n",
    "    for key in keys:\n",
    "        if key.letter == symbol:\n",
    "            return key\n",
    "    return None\n",
    "\n",
    "def beam_prune(p, obs_seq, obs_seq_index):\n",
    "    bp = obs_seq.beam[obs_seq_index]\n",
    "    if p > bp:\n",
    "        obs_seq.beam[obs_seq_index] = p\n",
    "        return False\n",
    "    else:\n",
    "        if p < beam_width * bp:\n",
    "            return True\n",
    "        else:\n",
    "            return False\n",
    "\n",
    "def propagate(token, obs_seq, completed_tokens):\n",
    "    ix = token.obs_index\n",
    "    next_ix = ix + 1\n",
    "    if next_ix == len(obs_seq.seq):\n",
    "        completed_tokens.append(token)\n",
    "    else:\n",
    "        next_obs = obs_seq.seq[next_ix]\n",
    "        for symbol in symbols:\n",
    "            if symbol == 'ɛ':\n",
    "                acc_prob = token.acc_prob * deletion_penalty\n",
    "                new_token = Token(token.hypo, acc_prob, next_ix)\n",
    "                if not beam_prune(acc_prob, obs_seq, next_ix):\n",
    "                    propagate(new_token, obs_seq, completed_tokens)\n",
    "            else:\n",
    "                key = get_key_for_symbol(symbol)\n",
    "                prior = get_prior(token.hypo+symbol, lm)\n",
    "                likelihood = get_likelihood(next_obs.x, next_obs.y, key)\n",
    "                acc_prob = token.acc_prob * prior * likelihood\n",
    "                new_token = Token(token.hypo+symbol, acc_prob, next_ix)\n",
    "                if not beam_prune(acc_prob, obs_seq, next_ix):\n",
    "                    propagate(new_token, obs_seq, completed_tokens)\n",
    "\n",
    "seed_token = Token(\"\", 1.0, -1)\n",
    "propagate(seed_token, obs_seq, completed_tokens)\n",
    "results = sorted(completed_tokens, key=lambda token: token.acc_prob, reverse=True)\n",
    "\n",
    "# Print out sequence hypotheses and their probabilities\n",
    "print_sorted_tokens(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Decoding with Substitutions, Deletions and Insertions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our latest decoder can substitute observations into either letters or empty strings ($\\epsilon$-transitions). However, it cannot model the fact that a user may have accidentally failed to touch a letter key. To correct this behaviour we need to model <em>insertion</em> transitions. This can be done by adding <em>self-loops</em> to each state, see the figure below. Each self-loop can generate any symbol."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To model such insertions we modify our <em>propagate</em> function again. Instead of simply propagating $k+1$ tokens to observation index $i+1$, we also propagate $k$ tokens to observation index $i$, the observation index of the token under consideration."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our search space is no longer increasing exponentially: it is infinite. The reason is that the decoding process can stay at observation index $i$ indefinately and keep generating new tokens. Therefore insertions need to carry an <em>insertion penalty</em>. If the beam width and the insertion penalty are set correctly, the search becomes tractable again.\n",
    "\n",
    "<img src=\"imgs/simple_self-loop_chain.svg\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "<strong>REFLECT:</strong> \n",
    "<ul>\n",
    "    <li>\n",
    "        Examine how adjusting the insertion penalty affects the length of the resulting hypotheses (advise not going above ~3).\n",
    "    </li>\n",
    "    <li>\n",
    "        How does the concept of an insertion penalty loosely relate to the way we have modelled the likelihood of an observation?\n",
    "    </li>\n",
    "</ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "letter sequence\n",
      "nth best, hypothesis, probability\n",
      "1, bad, 0.0000012\n",
      "2, bcd, 0.0000009\n",
      "3, dad, 0.0000008\n",
      "4, bac, 0.0000008\n",
      "5, bab, 0.0000007\n",
      "6, bbd, 0.0000007\n",
      "7, aad, 0.0000007\n",
      "8, bd, 0.0000006\n",
      "9, ad, 0.0000006\n",
      "10, ba, 0.0000006\n",
      "11, dcd, 0.0000006\n",
      "12, bcc, 0.0000006\n",
      "13, dac, 0.0000005\n",
      "14, bcb, 0.0000005\n",
      "15, bdd, 0.0000005\n",
      "...\n",
      "213, aaabd, 1.16e-10\n",
      "214, aaacd, 1.16e-10\n",
      "215, aaadd, 1.16e-10\n",
      "216, aaaca, 7.42e-11\n",
      "217, aaacb, 7.42e-11\n",
      "218, aaacc, 7.42e-11\n",
      "219, aaacd, 7.42e-11\n",
      "220, aaaba, 6.71e-11\n",
      "221, aaabb, 6.71e-11\n",
      "222, aaabc, 6.71e-11\n",
      "223, aaabd, 6.71e-11\n",
      "224, aaaaa, 4.28e-11\n",
      "225, aaaab, 4.28e-11\n",
      "226, aaaac, 4.28e-11\n",
      "227, aaaad, 4.28e-11\n"
     ]
    }
   ],
   "source": [
    "class ObservationSequence:\n",
    "    def __init__(self, seq):\n",
    "        self.seq = seq\n",
    "        self.beam = [0.0]*len(seq) \n",
    "\n",
    "symbols = ['a','b','c','d','ɛ'] #Epsilon represents a delete observation event (i.e. do not generate a hypothesis)\n",
    "#obs_seq = ObservationSequence([]) #No observations\n",
    "#obs_seq = ObservationSequence([Observation(0.1,1.1),Observation(0.15,-0.1),Observation(1.05,0.95)]) #bad + noise\n",
    "obs_seq = ObservationSequence([Observation(0.1,1.1),Observation(0.15,-0.1),Observation(3.0,-2.7),Observation(1.05,0.95)]) #b - a - (Noisy input far away from any key) - d + noise\n",
    "completed_tokens = []\n",
    "deletion_penalty = 0.02\n",
    "insertion_penalty = 1.2\n",
    "beam_width = 0.2\n",
    "\n",
    "def get_key_for_symbol(symbol):\n",
    "    for key in keys:\n",
    "        if key.letter == symbol:\n",
    "            return key\n",
    "    return None\n",
    "\n",
    "def beam_prune(p, obs_seq, obs_seq_index):\n",
    "    bp = obs_seq.beam[obs_seq_index]\n",
    "    if p > bp:\n",
    "        obs_seq.beam[obs_seq_index] = p\n",
    "        return False\n",
    "    else:\n",
    "        if p < beam_width * bp:\n",
    "            return True\n",
    "        else:\n",
    "            return False\n",
    "\n",
    "def propagate(token, obs_seq, completed_tokens):\n",
    "    # Which observation is the token in?\n",
    "    ix = token.obs_index\n",
    "    # Are there more observations?\n",
    "    next_ix = ix + 1\n",
    "    if next_ix == len(obs_seq.seq):\n",
    "        # No more observations, add token to the list of completed tokens\n",
    "        completed_tokens.append(token)\n",
    "    else:\n",
    "        # There is another observation, propagage tokens for all symbols to the next observation index\n",
    "        next_obs = obs_seq.seq[next_ix]\n",
    "        for symbol in symbols:\n",
    "            if symbol == 'ɛ':\n",
    "                acc_prob = token.acc_prob * deletion_penalty\n",
    "                new_token = Token(token.hypo, acc_prob, next_ix)\n",
    "                if not beam_prune(acc_prob, obs_seq, next_ix):\n",
    "                    propagate(new_token, obs_seq, completed_tokens)\n",
    "            else:\n",
    "                key = get_key_for_symbol(symbol)\n",
    "                prior = get_prior(token.hypo+symbol, lm)\n",
    "                likelihood = get_likelihood(next_obs.x, next_obs.y, key)\n",
    "                acc_prob = token.acc_prob * prior * likelihood\n",
    "                new_token = Token(token.hypo+symbol, acc_prob, next_ix)\n",
    "                if not beam_prune(acc_prob, obs_seq, next_ix):\n",
    "                    propagate(new_token, obs_seq, completed_tokens)\n",
    "    # Propagate tokens for all symbols (except epsilon) within the same observation index\n",
    "    if ix >= 0:\n",
    "        for symbol in symbols:\n",
    "            if not symbol == 'ɛ':\n",
    "                key = get_key_for_symbol(symbol)\n",
    "                prior = get_prior(token.hypo+symbol, lm)\n",
    "                obs = obs_seq.seq[ix]\n",
    "                acc_prob = token.acc_prob * prior * insertion_penalty\n",
    "                # Create a new token but don't increment the observation index\n",
    "                new_token = Token(token.hypo+symbol, acc_prob, ix)\n",
    "                if not beam_prune(acc_prob, obs_seq, ix):\n",
    "                    propagate(new_token, obs_seq, completed_tokens)\n",
    "\n",
    "seed_token = Token(\"\", 1.0, -1)\n",
    "propagate(seed_token, obs_seq, completed_tokens)\n",
    "results = sorted(completed_tokens, key=lambda token: token.acc_prob, reverse=True)\n",
    "\n",
    "# Print out sequence hypotheses and their probabilities\n",
    "print_sorted_tokens(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see from the output that our new decoder is capable of generating longer hypotheses than there are observations. However, the most likely hypothesis remains the same as before. The reason is our parameter settings and our language model. Since our language model simply distributes the probability mass evenly over our four possible letters and does not model probabilities over letter sequences, there is no structure that allows the decoder to sophisticatedly decide to insert a letter. For insertions to be used effectively, a high-quality language model is required."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Log probabilities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The previous code works but will not scale very well due to numerical underflow. To avoid this, statistical decoding is usually working with probabilities on a logarithmic scale (_log-probs_).\n",
    "\n",
    "To transform a probability $P$ to the logarithmic scale we compute the natural logarithm, $\\ln(P)$. To transform a log-prob back into a regular probability we compute $\\exp(\\ln(P))$. Adding log-probs has the same effect as multiplying regular probabilties.\n",
    "\n",
    "An example: \n",
    "\n",
    "Let $P_1 = P_2 = 0.5$. \n",
    "\n",
    "Then $\\ln(P_1) = \\ln(P_2) \\approx -0.693$. \n",
    "\n",
    "$\\ln(P_1) + \\ln(P_2) \\approx -1.386$ and $\\exp\\left[\\ln(P_1) + \\ln(P_2)\\right] = 0.25$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Multiplied Probabilities - Summed Log Probabilities\n",
      "0.25 - 0.25\n"
     ]
    }
   ],
   "source": [
    "P1 = 0.5\n",
    "P2 = 0.5\n",
    "\n",
    "joint_P1_P2 = P1 * P2\n",
    "\n",
    "# Take log\n",
    "log_P1 = math.log(P1)\n",
    "log_P2 = math.log(P2)\n",
    "\n",
    "# Sum log probs\n",
    "log_joint_P1_P2 = log_P1 + log_P2\n",
    "\n",
    "# Take exp(log(P))\n",
    "e_log_joint_P1_P2 = math.exp(log_joint_P1_P2)\n",
    "\n",
    "print('Multiplied Probabilities - Summed Log Probabilities')\n",
    "print('%.3g - %.3g' % (joint_P1_P2,e_log_joint_P1_P2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>\n",
    "<div class=\"alert alert-block alert-success\"> \n",
    "<h3>Exercise 2: Use Log Probabilities</h3> \n",
    "<br>\n",
    "Extend the decoder code below to use log probabilities throughout. Check that the results returned are the same as those above.\n",
    "</div>\n",
    "</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "letter sequence\n",
      "nth best, hypothesis, probability\n",
      "1, bad, 0.0000012\n",
      "2, bcd, 0.0000009\n",
      "3, dad, 0.0000008\n",
      "4, bac, 0.0000008\n",
      "5, bab, 0.0000007\n",
      "6, bbd, 0.0000007\n",
      "7, aad, 0.0000007\n",
      "8, bd, 0.0000006\n",
      "9, ad, 0.0000006\n",
      "10, ba, 0.0000006\n",
      "11, dcd, 0.0000006\n",
      "12, bcc, 0.0000006\n",
      "13, dac, 0.0000005\n",
      "14, bcb, 0.0000005\n",
      "15, bdd, 0.0000005\n",
      "...\n",
      "213, aaabd, 1.16e-10\n",
      "214, aaacd, 1.16e-10\n",
      "215, aaadd, 1.16e-10\n",
      "216, aaaca, 7.42e-11\n",
      "217, aaacb, 7.42e-11\n",
      "218, aaacc, 7.42e-11\n",
      "219, aaacd, 7.42e-11\n",
      "220, aaaba, 6.71e-11\n",
      "221, aaabb, 6.71e-11\n",
      "222, aaabc, 6.71e-11\n",
      "223, aaabd, 6.71e-11\n",
      "224, aaaaa, 4.28e-11\n",
      "225, aaaab, 4.28e-11\n",
      "226, aaaac, 4.28e-11\n",
      "227, aaaad, 4.28e-11\n"
     ]
    }
   ],
   "source": [
    "class ObservationSequence:\n",
    "    def __init__(self, seq):\n",
    "        self.seq = seq\n",
    "        self.beam = [0.0]*len(seq) \n",
    "\n",
    "symbols = ['a','b','c','d','ɛ'] #Epsilon represents a delete observation event (i.e. do not generate a hypothesis)\n",
    "#obs_seq = ObservationSequence([]) #No observations\n",
    "#obs_seq = ObservationSequence([Observation(0.1,1.1),Observation(0.15,-0.1),Observation(1.05,0.95)]) #bad + noise\n",
    "obs_seq = ObservationSequence([Observation(0.1,1.1),Observation(0.15,-0.1),Observation(3.0,-2.7),Observation(1.05,0.95)]) #b - a - (Noisy input far away from any key) - d + noise\n",
    "completed_tokens = []\n",
    "deletion_penalty = 0.02\n",
    "insertion_penalty = 1.2\n",
    "beam_width = 0.2\n",
    "\n",
    "def get_key_for_symbol(symbol):\n",
    "    for key in keys:\n",
    "        if key.letter == symbol:\n",
    "            return key\n",
    "    return None\n",
    "\n",
    "def beam_prune(p, obs_seq, obs_seq_index):\n",
    "    bp = obs_seq.beam[obs_seq_index]\n",
    "    if p > bp:\n",
    "        obs_seq.beam[obs_seq_index] = p\n",
    "        return False\n",
    "    else:\n",
    "        if p < beam_width * bp:\n",
    "            return True\n",
    "        else:\n",
    "            return False\n",
    "\n",
    "def propagate(token, obs_seq, completed_tokens):\n",
    "    # Which observation is the token in?\n",
    "    ix = token.obs_index\n",
    "    # Are there more observations?\n",
    "    next_ix = ix + 1\n",
    "    if next_ix == len(obs_seq.seq):\n",
    "        # No more observations, add token to the list of completed tokens\n",
    "        completed_tokens.append(token)\n",
    "    else:\n",
    "        # There is another observation, propagage tokens for all symbols to the next observation index\n",
    "        next_obs = obs_seq.seq[next_ix]\n",
    "        for symbol in symbols:\n",
    "            if symbol == 'ɛ':\n",
    "                acc_prob = token.acc_prob * deletion_penalty\n",
    "                new_token = Token(token.hypo, acc_prob, next_ix)\n",
    "                if not beam_prune(acc_prob, obs_seq, next_ix):\n",
    "                    propagate(new_token, obs_seq, completed_tokens)\n",
    "            else:\n",
    "                key = get_key_for_symbol(symbol)\n",
    "                prior = get_prior(token.hypo+symbol, lm)\n",
    "                likelihood = get_likelihood(next_obs.x, next_obs.y, key)\n",
    "                acc_prob = token.acc_prob * prior * likelihood\n",
    "                new_token = Token(token.hypo+symbol, acc_prob, next_ix)\n",
    "                if not beam_prune(acc_prob, obs_seq, next_ix):\n",
    "                    propagate(new_token, obs_seq, completed_tokens)\n",
    "    # Propagate tokens for all symbols (except epsilon) within the same observation index\n",
    "    if ix >= 0:\n",
    "        for symbol in symbols:\n",
    "            if not symbol == 'ɛ':\n",
    "                key = get_key_for_symbol(symbol)\n",
    "                prior = get_prior(token.hypo+symbol, lm)\n",
    "                obs = obs_seq.seq[ix]\n",
    "                acc_prob = token.acc_prob * prior * insertion_penalty\n",
    "                # Create a new token but don't increment the observation index\n",
    "                new_token = Token(token.hypo+symbol, acc_prob, ix)\n",
    "                if not beam_prune(acc_prob, obs_seq, ix):\n",
    "                    propagate(new_token, obs_seq, completed_tokens)\n",
    "\n",
    "seed_token = Token(\"\", 1.0, -1)\n",
    "propagate(seed_token, obs_seq, completed_tokens)\n",
    "results = sorted(completed_tokens, key=lambda token: token.acc_prob, reverse=True)\n",
    "\n",
    "# Print out sequence hypotheses and their probabilities\n",
    "print_sorted_tokens(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Statistical Language Modelling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A _statistical language model_ is a model that assigns probabilities to sequences of letters or words."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Information, Entropy and Perplexity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To understand how to model languages we need to understand what a _message_ is and how we can quantify the information in a message. A message $m$ in some message space $M$ is sent from a _sender_ to a _receiver_. Importantly, for a message $m$ to carry information $I(m)$ from the sender to the recevier the information must not already be known to the receiver. If the information is already known $I(m)=0$. $I(m)$ measures the amount of new information, or _surprisal_, there is in receiving a message. The more surprising the message is to the receiver, the more information content in the message.\n",
    "\n",
    "$I(m)$ is called _self-information_ and is defined as:\n",
    "\n",
    "\\begin{equation}\n",
    "I(m) = \\log_2\\left(\\frac{1}{P(m)}\\right)=-\\log_2\\left(P\\left(m\\right)\\right),\n",
    "\\end{equation}\n",
    "\n",
    "where $P(m)$ is the probability of the message in the message space $M$.\n",
    "\n",
    "The above equation is better understood by considering a message as communicating an event (such as the outcome of flipping a fair coin: heads or tails). Let information about one event $m_1$ be $I(m_1)$ and another event $m_2$ be $I(m_2)$.\n",
    "\n",
    "Then, _assuming these events are independent_, the information about _both_ these events would be $I_{both} = I(m_1\\cap m_2)= I(m_1) + I(m_2)$.\n",
    "\n",
    "If we now think about the _probabilities_ of these events we have $P(m_1)$ and $P(m_2)$ and $P_{both} = P(m_1\\cap m_2)= P(m_1) \\cdot P(m_2)$.\n",
    "\n",
    "In other words, self-information of independent events _adds_ while the probabilities of the events _multiply_. This is exactly how log-probs operate, which explains the logarithm in the above equation. Since logarithms scale with a constant factor we can set this constant factor to be negative in order to ensure information is a positive quantity. The base of the logarithm determines the unit. By convention the base is usually 2 ($\\log_2)$ and thus the unit is in _bits_.\n",
    "\n",
    "<h4>Example 1:</h4>\n",
    "\n",
    "Consider a Bernoulli trial involving flipping a fair coin. The probability of either outcome is $0.5$. To communicate the specific outcome of an individual event to the sender is the equivalent of sending a message $m$ with $I(m) = \\log_2\\left(\\frac{1}{0.5}\\right) =  \\log_2(2) = 1$ bit.\n",
    "\n",
    "<h4>Example 2:</h4>\n",
    "\n",
    "Consider we wish to communicate the outcome of having flipped a fair coin four times as either $m_a$ (four heads in a row) or $\\bar{m_a}$ (_not_ observing four heads in a row). The probability $P(m_a)$ of observing four heads in a row is $\\frac{1}{16}$ and $I(m_a) = -\\log_2\\left(P\\left(m_a\\right)\\right) = -\\log_2\\left(\\frac{1}{16}\\right) = 4$ bits. Communicating the other possible outcome, $\\bar{m_a}$, we first find the probability $P(\\bar{m_a}) = 1 - P(m_a) = 1 - \\frac{1}{16} = \\frac{15}{16}$. Then $I(\\bar{m_a}) = -\\log_2\\left(\\frac{15}{16}\\right) \\approx 0.093$ bits.\n",
    "\n",
    "Entropy $H$ is a measure of the _uncertainty_ in the message space $M$.\n",
    "\n",
    "\\begin{equation}\n",
    "H(M) = \\sum_{m \\in M} P(m)I(m)=-\\sum_{m \\in M}P(m)\\log_2\\left(P(m)\\right)\n",
    "\\end{equation}\n",
    "\n",
    "Entropy is the _average_ self-information of all messages in the message space. Entropy measures disorder. An entropy of zero means there is no disorder and everything is completely predictable.\n",
    "\n",
    "Entropy tells us the average bits we need to communicate a message, assuming we have found an optimal coding scheme. The difference between the average bits actually used for communication and the optimal bits truly necessary to encode all messages is a measure of _redundancy_. Natural langauges have high redundancy.\n",
    "\n",
    "A measure related to entropy is _perplexity_, which measures how well a probability model is at prediction. Perplexity is defined as $2^H$, where $H$ is the entropy of the model.\n",
    "\n",
    "<h4>Example from Jurafsky and Martin 2000; _Speech and Language Processing_:</h4>\n",
    "\n",
    "We want to bet on eight horses via a middleman. How do we communicate the horse we want to bet on?. One way is to encode this message using binary representation for each horse, horse 1 becomes 001, horse 2 becomes 010, etc. If we spend the whole day betting and each horse is coded with 3 bits then on average we are sending 3 bits per race.\n",
    "\n",
    "Can we do better? Assume each horse is not equally likely but instead distributed like this: horse 1: $\\frac{1}{2}$; horse 2: $\\frac{1}{4}$; horse 3: $\\frac{1}{8}$; horse 4: $\\frac{1}{16}$; and horses 5--8: $\\frac{1}{64}$. Then the entropy of the random variable $X$ that ranges over horses gives us a lower bound on the number of bits: $H(X) = -\\sum^8_{i=1}P(i)\\log_2\\left(P(i)\\right)=2$ bits. Such a code can be built by using shorter encoding for more probable horses and longer encodings for less probable horses (for example, the most likely horse can be encoded as 0 while the least likely horse can be coded as 111111. If the horses have equal probability then if we repeat the above calculation we will find that the entropy would be 3.\n",
    "\n",
    "Perplexity is the weighted average number of choices a random variable has to make. Choosing between 8 equally likely horses ($H = 3$ bits) has a perplexity of $2^3 = 8$. Choosing between the biased horses above ($H = 2$ bits) has a perplexity of $2^2 = 4$.\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Modelling language"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since natural languages are highly redundant it follows that the vast majority of letter and word combinations are highly improbable. We can capture probable sequences of letters and words by counting their occurance in some large representative dataset, a _corpus_. A _language model_ assigns probabilities to sequences of letters and words.\n",
    "\n",
    "We will demonstrate these ideas by using an example from: Shannon, C.E. 1948. A mathematical theory of communication. _Bell System Technical Journal_ 27: 379–423, 623–656.\n",
    "\n",
    "Consider a discrete source with five letters \\{A,B,C,D,E\\}. Chosen with probability of 0.2 (uniformaly sampled) and each choice independent of previous choices we might expect a typical outcome such as: B D C B C E C C C A D C B D D A A E C E E A A B B D A E E C A C E E B A E E C B C E A D\n",
    "\n",
    "We can build a generative model that generates increasingly English-like language output, as follows:\n",
    "\n",
    "<b>Zero-order approximation (uniform random sampling from {A-Z}:</b>\n",
    "\n",
    "<em>XFOML RXKHRJFFJUJ ZLPWCFWKCYJ FFJEYVKCQSGHYD QPAAMKBZAACIBZLHJQD</em>\n",
    "\n",
    "<b>First-order approximation (sampling according to the probabilities of English letters occurring in common texts):</b>\n",
    "\n",
    "<em>OCRO HLI RGWR NMIELWIS EU LL NBNESEBYA TH EEI ALHENHTTPA OOBTTVA NAH BRL</em>\n",
    "\n",
    "<b>Second-order approximation (sampling so that the every letter is dependent on its previous letter according to letter probabilities in English):</b>\n",
    "\n",
    "<em>ON IE ANTSOUTINYS ARE T INCTORE ST BE S DEAMY ACHIN D ILONASIVE TUCOOWE AT TEASONARE FUSO TIZIN ANDY TOBE SEACE CTISBE</em>\n",
    "\n",
    "<b>Third-order approximation (sampling similar to second-order approximation, except each letter is now dependent on two previous letters):</b>\n",
    "\n",
    "<em>IN NO IST LAT WHEY CRATICT FROURE BIRS GROCID PONDENOME OF DEMONSTURES OF THE REPTAGIN IS REGOACTIONA OF CRE</em>\n",
    "\n",
    "<h4>Zipf's law</h4>\n",
    "\n",
    "Another striking example of language redudandancy is Zipf's law, which estimates the probability $P_r$ of occurrence of a word in a corpus to be $P_r \\propto \\frac{1}{r^\\alpha}$, where $r$ is the statistical rank of the word in decreasing order and $\\alpha$ is close to unity. As an example of this relationship, the 100 most frequently used words in the British National Corpus comprise close to 46% of the entire corpus.\n",
    "\n",
    "<h4><em>n</em>-gram models</h4>\n",
    "\n",
    "A _unigram model_ assumes any word will follow any other word at an equal probability. The simplest unigram model possible believes every word is equally probable. So if a language consists of $1000$ words then the probability of occurrence of a word is exactly $\\frac{1}{1000}$.\n",
    "\n",
    "The word \"the\" occurs 69,971 times in the Brown corpus while the word \"rabbit\" only occurs 11 times. However, what is the most reasonable choice of word to complete the following sentence fragment: \"Just then, the white...\". Clearly a unigram model makes too simplistic assumptions about natural languages.\n",
    "\n",
    "We can handle the case with the white rabbit by probabilistically modelling every single possible word combination that can occur in the world. However, it is clearly intractable to model such a large space, which is countably infinite. Therefore, language models make a _Markov assumption_: the probability of a word only depends on a finite set of previous words.\n",
    "\n",
    "A bigram model approximates the probability of a word given previous words by the conditional probability of a single preceding word. So instead of computing this probability: $P(\\text{rabbit}|\\text{Just the other day I saw a})$, a bigram model computes the following approximate probability: $P(\\text{rabbit}|\\text{a})$.\n",
    "\n",
    "A trigram model is a minor extension to a bigram model: it has two words of prior context. In general, an $n$-gram model has $n-1$ words of context an is said to be of order $n-1$.\n",
    "\n",
    "The basic way of training a language model is by counting instances of $n$-grams in corpora (frequency estimation). These counts have to be normalised in order to create a probabilistic model. Fundamentally, one takes the count of say a particular bigram and divide it by the sum of all the bigrams that share the first word (which is identical to the unigram count of the first word). This is a form of _parameter estimation_. In general, we estimate an $n$-gram probability by dividing the observed frequency of a particular sequence by the observed frequency of a prefix. Using such relative frequency to estimate probabilities is an example of maximum likelihood estimation (the likelihood of the training set given the model is maximised). Typically language models use log-probs.\n",
    "\n",
    "The typical way to evaluate a language model is to train it on a training set and test it on a separate test set. The lower the perplexity the better the language model is in predicting text on the test set."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Smoothing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our previous language modelling cannot account for $n$-grams we have not previously observed in our training data. Fundamentally, the problem is that our language model is overestimating the probabilities of observed sequences and underestimting unobserved sequences by assigning them zero probability mass. The solution is to distribute some of the probability mass in the language model to unobserved sequences.\n",
    "\n",
    "The simplest method is to add one to every count of every $n$-gram in the model (this is known as the Laplace method):\n",
    "\n",
    "\\begin{equation}\n",
    "c^*=c+1,\n",
    "\\end{equation}\n",
    "\n",
    "where $c$ is the original count and $c^*$ is the modified count. It does solve the problem of zero counts but in practice it results in poor language modelling performance.\n",
    "\n",
    "A slightly more sophisticated version is known as Good-Turing discounting and it is defined as:\n",
    "\n",
    "\\begin{equation}\n",
    "c^*=(c+1)\\frac{N_{c+1}}{N_c},\n",
    "\\end{equation}\n",
    "\n",
    "where $c$ is the original count and $c^*$ is the modified count as before and $N_c$ and $N_{c+1}$ are the number of $n$-grams that occur exactly $c$ and $c+1$ times. An obvious problem is when $N_c = 0$. One solution is to plot $\\log(c)$ against $\\log(Z_c)$, where $Z_c$ is:\n",
    "\n",
    "\\begin{equation}\n",
    "Z_c = \\frac{2N_c}{c''-c'},\n",
    "\\end{equation}\n",
    "\n",
    "where $c''$ and $c'$ are the nearest higher and lower sample counts such that $N_{c''}$ and $N_{c'}$ are positive. A linear fit between $\\log(c)$ against $\\log(Z_c)$ can then be obtained by ordinary least squares (OLS) regression and used to infer $N_c$ that are unobservable.\n",
    "\n",
    "See <a href=\"https://www.grsampson.net/AGtf.html\">Good–Turing Frequency Estimation Without Tears</a> for a complete explanation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>\n",
    "<div class=\"alert alert-block alert-success\"> \n",
    "<h3>Exercise 3: Improve the Statistical Decoder</h3> \n",
    "<br>\n",
    "The code for the statistical decoder you were provided with has a few limitations. In this exercise you will work in your group and introduce a better language model.  Change the uniform unigram model to a character-level language model that learns probabilities for sequences of characters based on a training text and uses smoothing to estimate the probabilities for unseen character sequences.\n",
    "</div>\n",
    "</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\"> \n",
    "    <strong>Task 1:</strong> Construct a character-level bigram language model based on training text. It is useful to represent the start of a word as a psuedo character.  \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  1.62400000e+03   1.08300000e+03   1.20400000e+03   1.09300000e+03]\n",
      " [  3.04000000e+02   1.36800000e+03   5.95000000e+02   1.20100000e+03]\n",
      " [  1.52700000e+03   3.18000000e+02   1.73000000e+02   1.00000000e+00]\n",
      " [  7.55000000e+02   1.00000000e+00   1.47000000e+02   2.95000000e+02]\n",
      " [  1.09300000e+03   1.00000000e+00   2.84000000e+02   1.41000000e+02]]\n",
      "[[  3.24540368e-01   2.16426859e-01   2.40607514e-01   2.18425260e-01]\n",
      " [  8.76585928e-02   3.94463668e-01   1.71568627e-01   3.46309112e-01]\n",
      " [  7.56315007e-01   1.57503715e-01   8.56859832e-02   4.95294700e-04]\n",
      " [  6.30217028e-01   8.34724541e-04   1.22704508e-01   2.46243740e-01]\n",
      " [  7.19552337e-01   6.58327847e-04   1.86965109e-01   9.28242265e-02]]\n"
     ]
    }
   ],
   "source": [
    "# Load the training text\n",
    "f_training_text = open('training_text.txt', 'r')\n",
    "training_text = f_training_text .read()\n",
    "f_training_text.close()\n",
    "\n",
    "# Initalize a structure to count character-to-character transitions, include a start-of-word one-way transition denoted '_'\n",
    "transition_counts = np.matlib.zeros((5,4))\n",
    "from_ids = {'_':0,'a':1,'b':2,'c':3,'d':4}\n",
    "to_ids = {'a':0,'b':1,'c':2,'d':3}\n",
    "\n",
    "# Split the training text into words\n",
    "training_words = training_text.split(' ')\n",
    "\n",
    "# Count character transitions\n",
    "for word in training_words:\n",
    "    prev_character = '_'\n",
    "    for character in word:\n",
    "        # TODO: Update the transition_counts matrix, make use of from_ids and to_ids\n",
    "        from_id = from_ids[prev_character] # REMOVE\n",
    "        to_id = to_ids[character] # REMOVE        \n",
    "        transition_counts[from_id,to_id] += 1 # REMOVE\n",
    "        \n",
    "        # Update previous character\n",
    "        prev_character = character\n",
    "\n",
    "# TODO: Apply smoothing (if necessary) and convert the transition counts into a useful bigram\n",
    "\n",
    "smoothed_transition_counts = transition_counts + np.matlib.ones((5,4)) # REMOVE\n",
    "print(smoothed_transition_counts) # REMOVE\n",
    "transition_probs = smoothed_transition_counts/smoothed_transition_counts.sum(axis=1) # REMOVE\n",
    "print(transition_probs) # REMOVE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\"> \n",
    "    <strong>Task 2:</strong> Modify the get_prior method to use use the computed bigram.  \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# REMOVE\n",
    "lm = {'_a':3.24540368e-01,'_b':2.16426859e-01,'_c':2.40607514e-01,'_d':2.18425260e-01,\\\n",
    "'aa':8.76585928e-02,'ab':3.94463668e-01,'ac':1.71568627e-01,'ad':3.46309112e-01,\\\n",
    "'ba':7.56315007e-01,'bb':1.57503715e-01,'bc':8.56859832e-02,'bd':4.95294700e-04,\\\n",
    "'ca':6.30217028e-01,'cb':8.34724541e-04,'cc':1.22704508e-01,'cd':2.46243740e-01,\\\n",
    "'da':7.19552337e-01,'db':6.58327847e-04,'dc':1.86965109e-01,'dd':9.28242265e-02}\n",
    "\n",
    "# Populate the bigram with computed probabilities\n",
    "#lm = {'_a':0.25,'_b':0.25,'_c':0.25,'_d':0.25,\\\n",
    "#      'aa':0.25,'ab':0.25,'ac':0.25,'ad':0.25,\\\n",
    "#      'ba':0.25,'bb':0.25,'bc':0.25,'bd':0.25,\\\n",
    "#      'ca':0.25,'cb':0.25,'cc':0.25,'cd':0.25,\\\n",
    "#      'da':0.25,'db':0.25,'dc':0.25,'dd':0.25}\n",
    "\n",
    "# Modify get_prior to return prior for specified context\n",
    "#def get_prior(context, lm):\n",
    "#    c = context[-1]\n",
    "#    return lm[c]\n",
    "\n",
    "# REMOVE\n",
    "def get_prior(context, lm):\n",
    "    if (len(context) == 1):\n",
    "        context = '_' + context\n",
    "    c = context[-2:]\n",
    "    return lm[c]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\"> \n",
    "    <strong>Task 3:</strong> Examine the effect of the new language model on decoder performance.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "letter sequence\n",
      "nth best, hypothesis, probability\n",
      "1, bad, 0.0002221\n",
      "2, bada, 0.0001917\n",
      "3, bab, 0.0001459\n",
      "4, dad, 0.0001429\n",
      "5, baba, 0.0001324\n",
      "6, dada, 0.0001234\n",
      "7, babad, 0.0000954\n",
      "8, babad, 0.0000954\n",
      "9, dab, 0.0000939\n",
      "10, badab, 0.0000908\n",
      "11, abad, 0.0000865\n",
      "12, daba, 0.0000852\n",
      "13, baba, 0.0000844\n",
      "14, badaba, 0.0000824\n",
      "15, babada, 0.0000824\n",
      "...\n",
      "217, aababa, 0.0000054\n",
      "218, aababa, 0.0000054\n",
      "219, aabad, 0.0000052\n",
      "220, aaca, 0.0000051\n",
      "221, aadab, 0.0000050\n",
      "222, aadc, 0.0000048\n",
      "223, aabada, 0.0000045\n",
      "224, aadaba, 0.0000045\n",
      "225, aaa, 0.0000020\n",
      "226, aaab, 0.0000009\n",
      "227, aaaba, 0.0000008\n",
      "228, aaad, 0.0000008\n",
      "229, aaada, 0.0000007\n",
      "230, aaac, 0.0000004\n",
      "231, aaabab, 0.0000004\n"
     ]
    }
   ],
   "source": [
    "#obs_seq = ObservationSequence([]) #No observations\n",
    "#obs_seq = ObservationSequence([Observation(0.1,1.1)])\n",
    "obs_seq = ObservationSequence([Observation(0.1,1.1),Observation(0.15,-0.1),Observation(1.05,0.95)]) #bad + noise\n",
    "#obs_seq = ObservationSequence([Observation(0.1,1.1),Observation(0.15,-0.1),Observation(3.0,-2.7),Observation(1.05,0.95)]) #b - a - (Noisy input far away from any key) - d + noise\n",
    "completed_tokens = []\n",
    "deletion_penalty = 0.02\n",
    "insertion_penalty = 1.2\n",
    "beam_width = 0.2\n",
    "\n",
    "def propagate(token, obs_seq, completed_tokens):\n",
    "    # Which observation is the token in?\n",
    "    ix = token.obs_index\n",
    "    # Are there more observations?\n",
    "    next_ix = ix + 1\n",
    "    if next_ix == len(obs_seq.seq):\n",
    "        # No more observations, add token to the list of completed tokens\n",
    "        completed_tokens.append(token)\n",
    "    else:\n",
    "        # There is another observation, propagate tokens for all symbols to the next observation index\n",
    "        next_obs = obs_seq.seq[next_ix]\n",
    "        for symbol in symbols:\n",
    "            if symbol == 'ɛ':\n",
    "                acc_prob = token.acc_prob * deletion_penalty\n",
    "                new_token = Token(token.hypo, acc_prob, next_ix)\n",
    "                if not beam_prune(acc_prob, obs_seq, next_ix):\n",
    "                    propagate(new_token, obs_seq, completed_tokens)\n",
    "            else:\n",
    "                key = get_key_for_symbol(symbol)\n",
    "                prior = get_prior(token.hypo+symbol, lm)\n",
    "                likelihood = get_likelihood(next_obs.x, next_obs.y, key)\n",
    "                acc_prob = token.acc_prob * prior * likelihood\n",
    "                new_token = Token(token.hypo+symbol, acc_prob, next_ix)\n",
    "                if not beam_prune(acc_prob, obs_seq, next_ix):\n",
    "                    propagate(new_token, obs_seq, completed_tokens)\n",
    "    # Propagate tokens for all symbols (except epsilon) within the same observation index\n",
    "    if ix >= 0:\n",
    "        for symbol in symbols:\n",
    "            if not symbol == 'ɛ':\n",
    "                key = get_key_for_symbol(symbol)\n",
    "                prior = get_prior(token.hypo+symbol, lm)\n",
    "                obs = obs_seq.seq[ix]\n",
    "                acc_prob = token.acc_prob * prior * insertion_penalty\n",
    "                # Create a new token but don't increment the observation index\n",
    "                new_token = Token(token.hypo+symbol, acc_prob, ix)\n",
    "                if not beam_prune(acc_prob, obs_seq, ix):\n",
    "                    propagate(new_token, obs_seq, completed_tokens)\n",
    "\n",
    "seed_token = Token(\"\", 1.0, -1)\n",
    "propagate(seed_token, obs_seq, completed_tokens)\n",
    "results = sorted(completed_tokens, key=lambda token: token.acc_prob, reverse=True)\n",
    "\n",
    "# Print out sequence hypotheses and their probabilities\n",
    "print_sorted_tokens(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extensions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Several extensions are possible to gradually augment the above models. Some of these are explained below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adaptive language modelling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We previously developed an $n$-gram model to use as a model for our prior beliefs in the hypotheses we generate during the search. Our $n$-gram model was trained on general text. However, not all users are interested in typing the same text. This observation leads to the idea of _adaptation_. An adaptive language model learns sequences of letters (or words) that the user prefers to type. How to best adapt to users' typing is a challenging and ongoing area of research.\n",
    "\n",
    "One of the simplest ways is to train a _background_ model and a _foreground_ model. The background model is trained on a large corpus and the foreground model is trained online on text actually typed by the user (by updating the counts in the language model online). Since background models tend to be large it would take a very long time for the users' updated counts to be reflected in the prior. A faster way is to use a mixture model where the log-prob estimates for the background and foreground models are linearly combined with a mixture weight $\\alpha$:\n",
    "\n",
    "\\begin{equation}\n",
    "\\ln\\left(P_{mixture}\\right)=\\alpha\\ln\\left(P_{background}\\right)+(1-\\alpha)\\ln\\left(P_{foreground}\\right)\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Confusion networks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Our current decoder only searches for the highest probable hypothesis. However, it is possible to generate a hypothesis space. One useful representation of such a hypotheiss space is a _confusion network_. A confusion network is a set of hypotheses and corresponding posterior probabilities for each observation index. The posterior probabilities are normalized so that they sum to unity. Each observation index has a corresponding _confusion cluster_. A series of confusion clusters form the confusion network.\n",
    "\n",
    "It is relatively easy to change the decoder code to generate confusion networks. These confusion networks can then be outputted from the decoder and used to for example offer the user a selection of next-best choices if the most probable hypothesis was the uninteded one."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parallelisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Token passing lends itself to parallelisation as we propagate new tokens and mostly read from an immutable observation sequence. This means tokens can be propagated in multiple cores with minimal synchronisation. The beam width is however updated during a search and needs to either be an atomic read/write operation, have a mutual exclusion lock, or be separate for each core."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\"> \n",
    "<h1>Further Exploration: Grid Search</h1>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we have designed a probabilistic model for our problem and implemented a token-passing decoder to infer hypotheses from observations, we need to set the parameters. There are many ways to learn parameters. The simplest method is _grid search_, which fixes all parameters except one and then searches for the optimal parameter value under the fixed constraints. The process can be iterated.\n",
    "\n",
    "A more advanced method is to use the expectation-maximisation (EM) algorithm to learn the parameter set. This is indeed the foundation of the Baum-Welch algorithm for hidden Markov models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Metric"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To perform grid search, we need a suitable metric to optimise.  A useful metric for comparing the desired letter sequence to the letter sequence actually reproduced is the character error rate (CER).  Here we define character error rate as the number of changes required to transform the received letter sequence (response) into the desired letter sequence (stimulus) divided by the length of the desired letter sequence.\n",
    "\n",
    "The Levenshtein distance is the number of single-character edits (can be any of: insertion, deletion or substitution) required to transform string A into string B.  The Levenshtein distance provides the numerator when computing character error rate. An example is included below.\n",
    "\n",
    "When performing grid search we thus want to minimise the average character error rate across the test set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "stimulus: coffee, has 6 characters\n",
      "2 single-character edits are required to transform covfefe into coffee\n",
      "CER is 0.33\n"
     ]
    }
   ],
   "source": [
    "def character_error_rate(stimulus,response):\n",
    "    cer = 0.0;\n",
    "    ld = float(levenshtein(stimulus,response))\n",
    "    cer = ld / len(stimulus)\n",
    "\n",
    "    return cer\n",
    "\n",
    "# Credit to: https://en.wikibooks.org/wiki/Algorithm_Implementation/Strings/Levenshtein_distance#Python\n",
    "def levenshtein(s1, s2):\n",
    "    if len(s1) < len(s2):\n",
    "        return levenshtein(s2, s1)\n",
    "\n",
    "    # len(s1) >= len(s2)\n",
    "    if len(s2) == 0:\n",
    "        return len(s1)\n",
    "\n",
    "    previous_row = range(len(s2) + 1)\n",
    "    for i, c1 in enumerate(s1):\n",
    "        current_row = [i + 1]\n",
    "        for j, c2 in enumerate(s2):\n",
    "            insertions = previous_row[j + 1] + 1 # j+1 instead of j since previous_row and current_row are one character longer\n",
    "            deletions = current_row[j] + 1       # than s2\n",
    "            substitutions = previous_row[j] + (c1 != c2)\n",
    "            current_row.append(min(insertions, deletions, substitutions))\n",
    "        previous_row = current_row\n",
    "    \n",
    "    return previous_row[-1]\n",
    "\n",
    "stimulus = 'coffee'\n",
    "response = 'covfefe'\n",
    "\n",
    "print('stimulus: %s, has %d characters' % (stimulus,len(stimulus)))\n",
    "\n",
    "ld = levenshtein(stimulus,response)\n",
    "\n",
    "print('%d single-character edits are required to transform %s into %s' % (ld,response,stimulus))\n",
    "\n",
    "cer = character_error_rate(stimulus,response)\n",
    "\n",
    "print('CER is %.2f' % (cer))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\"> \n",
    "<h3>Exercise 4: Learn Parameters from User Data</h3> \n",
    "    \n",
    "Collect sample data and implement grid search to train your statistical decoder.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\"> \n",
    "    <strong>Task 1:</strong> Use the ABCD toy keyboard below and generate sample data. The toy keyboard will present a set of 10 stimulus words randomly choosen from a larger stimulus set.  Clicking next on the 10th stimulus word will close the keyboard and write a complete stimulus-response log to file.  You can repeat this task as many times as you want but the log file will overwrite. \n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "logging touches to: stimulus_response.log\n"
     ]
    }
   ],
   "source": [
    "from toy_keyboard import ToyKeyboard\n",
    "\n",
    "ToyKeyboard()\n",
    "%gui tk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\"> \n",
    "    <strong>Task 2:</strong> Implement grid search to find parameter values that minimise character error rates across the sample data you collected. You can explore the effect of modifying any of the parameters you think might be relevant.  You may want to manually inspect reasonable bounds for parameter values before starting the grid search. Code to parse and extract the stimulus-response log is provided below. You should be able to reduce CER to 0. \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sequence for: abad\n",
      "0.150, -0.300\n",
      "0.020, 0.820\n",
      "0.180, 0.030\n",
      "0.800, 0.890\n",
      "S-R: abad-abad\n",
      "CER: 0.00\n",
      "\n",
      "Sequence for: dac\n",
      "0.790, 1.240\n",
      "0.240, -0.210\n",
      "0.980, 0.150\n",
      "S-R: dac-dac\n",
      "CER: 0.00\n",
      "\n",
      "Sequence for: cdc\n",
      "0.900, 0.060\n",
      "0.180, 0.900\n",
      "0.920, 1.140\n",
      "S-R: cdc-cbd\n",
      "CER: 0.67\n",
      "\n",
      "Sequence for: bad\n",
      "0.180, 0.980\n",
      "0.250, 0.000\n",
      "0.850, 1.120\n",
      "S-R: bad-bad\n",
      "CER: 0.00\n",
      "\n",
      "Sequence for: baba\n",
      "0.310, 0.950\n",
      "0.190, -0.110\n",
      "-0.060, 1.030\n",
      "0.140, -0.120\n",
      "S-R: baba-baba\n",
      "CER: 0.00\n",
      "\n",
      "Sequence for: bac\n",
      "0.230, 1.010\n",
      "0.300, 0.150\n",
      "0.960, 0.010\n",
      "S-R: bac-bac\n",
      "CER: 0.00\n",
      "\n",
      "Sequence for: ca\n",
      "0.930, -0.110\n",
      "0.240, -0.150\n",
      "S-R: ca-ca\n",
      "CER: 0.00\n",
      "\n",
      "Sequence for: cc\n",
      "1.010, 0.070\n",
      "1.010, 0.070\n",
      "S-R: cc-cc\n",
      "CER: 0.00\n",
      "\n",
      "Sequence for: ac\n",
      "0.120, 0.050\n",
      "0.840, 0.110\n",
      "S-R: ac-ac\n",
      "CER: 0.00\n",
      "\n",
      "Sequence for: bab\n",
      "0.100, 1.180\n",
      "0.130, 0.140\n",
      "0.050, 0.990\n",
      "S-R: bab-bab\n",
      "CER: 0.00\n",
      "\n",
      "Mean CER across training set: 0.07\n",
      "0.0666666666667\n",
      "0.1\n",
      "0.005\n",
      "1.0\n",
      "0.05\n"
     ]
    }
   ],
   "source": [
    "from parse_sr_log import ParseStimulusResponseLog\n",
    "\n",
    "# Extract a matched array of stimuli and observation sequences from the stimulus-response.log\n",
    "stimulus_set, observation_sequence_set = ParseStimulusResponseLog('stimulus_response.log')\n",
    "\n",
    "# Default parameter values\n",
    "deletion_penalty = 0.01 #0.02\n",
    "insertion_penalty = 1.2\n",
    "beam_width = 0.2\n",
    "sigma = 0.25\n",
    "\n",
    "def get_likelihood(x, y, key):\n",
    "    key_x = key.centre_x\n",
    "    key_y = key.centre_y\n",
    "    p = multivariate_normal.pdf([x,y], mean=[key_x,key_y], cov=[[sigma,0],[0,sigma]]);\n",
    "    return p\n",
    "\n",
    "# Keep track of the character error rates\n",
    "cer_set = []\n",
    "\n",
    "# Print out stimuli and observation sequences\n",
    "stimulus_count = 0\n",
    "for observation_sequence in observation_sequence_set:\n",
    "    stimulus = stimulus_set[stimulus_count]\n",
    "    \n",
    "    # Print out observation points for stimulus\n",
    "    print(\"Sequence for: \" + stimulus)\n",
    "    for observation in observation_sequence.seq:\n",
    "        print('%.3f, %.3f' % (observation.x, observation.y))\n",
    "    stimulus_count += 1\n",
    "    \n",
    "    completed_tokens = []\n",
    "    \n",
    "    seed_token = Token(\"\", 1.0, -1)\n",
    "    propagate(seed_token, observation_sequence, completed_tokens)\n",
    "    results = sorted(completed_tokens, key=lambda token: token.acc_prob, reverse=True)\n",
    "\n",
    "    # Get best hypothesis from the decoder result list\n",
    "    best_response = results[0].hypo\n",
    "    \n",
    "    # Compute character error rate\n",
    "    cer = character_error_rate(stimulus,best_response)\n",
    "    cer_set.append(cer)\n",
    "    \n",
    "    print('S-R: %s-%s' % (stimulus,best_response))\n",
    "    print('CER: %.2f\\n' % (cer))\n",
    "    \n",
    "print('Mean CER across training set: %.2f' % (np.mean(cer_set)))\n",
    "\n",
    "# TODO: implement grid search to minimize CER\n",
    "\n",
    "min_mean_cer = 1.0\n",
    "best_dp = -1\n",
    "\n",
    "sigma_values = np.arange(0.1,1.0,0.1)\n",
    "for sigma in sigma_values:\n",
    "    def get_likelihood(x, y, key):\n",
    "        key_x = key.centre_x\n",
    "        key_y = key.centre_y\n",
    "        p = multivariate_normal.pdf([x,y], mean=[key_x,key_y], cov=[[sigma,0],[0,sigma]]);\n",
    "        return p\n",
    "    bw_values = np.arange(0.05,0.25,0.05)\n",
    "    for beam_width in bw_values:\n",
    "        ip_values = np.arange(1.0,1.5,0.25)\n",
    "        for insertion_penalty in ip_values:\n",
    "            dp_values = np.arange(0.005,0.02,0.005)\n",
    "            for deletion_penalty in dp_values:\n",
    "                # Keep track of the character error rates\n",
    "                cer_set = []\n",
    "\n",
    "                # Print out stimuli and observation sequences\n",
    "                stimulus_count = 0\n",
    "                for observation_sequence in observation_sequence_set:\n",
    "                    stimulus = stimulus_set[stimulus_count]\n",
    "                    stimulus_count += 1\n",
    "\n",
    "                    completed_tokens = []\n",
    "\n",
    "                    seed_token = Token(\"\", 1.0, -1)\n",
    "                    propagate(seed_token, observation_sequence, completed_tokens)\n",
    "                    results = sorted(completed_tokens, key=lambda token: token.acc_prob, reverse=True)\n",
    "\n",
    "                    # Get best hypothesis from the decoder result list                \n",
    "                    best_response = \"\"\n",
    "                    if (len(results) > 0):\n",
    "                        best_response = results[0].hypo\n",
    "\n",
    "                    # Compute character error rate\n",
    "                    cer = character_error_rate(stimulus,best_response)\n",
    "                    cer_set.append(cer)\n",
    "\n",
    "                # Compute mean CER across test set\n",
    "                mean_cer = np.mean(cer_set)\n",
    "                if (mean_cer < min_mean_cer):\n",
    "                    best_bw = beam_width\n",
    "                    best_dp = deletion_penalty\n",
    "                    best_ip = insertion_penalty\n",
    "                    best_sigma = sigma\n",
    "                    min_mean_cer = mean_cer\n",
    "\n",
    "print(min_mean_cer)\n",
    "print(best_sigma)\n",
    "print(best_dp)\n",
    "print(best_ip)\n",
    "print(best_bw)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\"> \n",
    "    <strong>Task 3:</strong> Use a noisy toy keyboard and generate sample data. Perform grid search again. Examine what noise level you can accomodate by adjusting sigma.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "logging touches to: stimulus_noisy_response.log\n"
     ]
    }
   ],
   "source": [
    "from noisy_toy_keyboard import NoisyToyKeyboard\n",
    "\n",
    "NoisyToyKeyboard(sigma=30)\n",
    "%gui tk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sequence for: add\n",
      "0.279, 0.153\n",
      "0.651, 1.316\n"
     ]
    },
    {
     "ename": "NameError",
     "evalue": "name 'Token' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-1-9ea1f7473177>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     33\u001b[0m     \u001b[0mcompleted_tokens\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     34\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 35\u001b[0;31m     \u001b[0mseed_token\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mToken\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1.0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     36\u001b[0m     \u001b[0mpropagate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mseed_token\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mobservation_sequence\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcompleted_tokens\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     37\u001b[0m     \u001b[0mresults\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msorted\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcompleted_tokens\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkey\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mlambda\u001b[0m \u001b[0mtoken\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mtoken\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0macc_prob\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mreverse\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'Token' is not defined"
     ]
    }
   ],
   "source": [
    "from parse_sr_log import ParseStimulusResponseLog\n",
    "\n",
    "# Extract a matched array of stimuli and observation sequences from the stimulus-response.log\n",
    "#stimulus_set, observation_sequence_set = ParseStimulusResponseLog('stimulus_noisy_response.log')\n",
    "stimulus_set, observation_sequence_set = ParseStimulusResponseLog('test_data.log')\n",
    "\n",
    "# Default parameter values\n",
    "deletion_penalty = 0.005 #0.02\n",
    "insertion_penalty = 1.0 #1.2\n",
    "beam_width = 0.05 #0.2\n",
    "sigma = 0.3 #1.0\n",
    "\n",
    "def get_likelihood(x, y, key):\n",
    "    key_x = key.centre_x\n",
    "    key_y = key.centre_y\n",
    "    p = multivariate_normal.pdf([x,y], mean=[key_x,key_y], cov=[[sigma,0],[0,sigma]]);\n",
    "    return p\n",
    "\n",
    "# Keep track of the character error rates\n",
    "cer_set = []\n",
    "\n",
    "# Print out stimuli and observation sequences\n",
    "stimulus_count = 0\n",
    "for observation_sequence in observation_sequence_set:\n",
    "    stimulus = stimulus_set[stimulus_count]\n",
    "    \n",
    "    # Print out observation points for stimulus\n",
    "    print(\"Sequence for: \" + stimulus)\n",
    "    for observation in observation_sequence.seq:\n",
    "        print('%.3f, %.3f' % (observation.x, observation.y))\n",
    "    stimulus_count += 1\n",
    "    \n",
    "    completed_tokens = []\n",
    "    \n",
    "    seed_token = Token(\"\", 1.0, -1)\n",
    "    propagate(seed_token, observation_sequence, completed_tokens)\n",
    "    results = sorted(completed_tokens, key=lambda token: token.acc_prob, reverse=True)\n",
    "\n",
    "    # Get best hypothesis from the decoder result list\n",
    "    best_response = results[0].hypo\n",
    "    \n",
    "    # Compute character error rate\n",
    "    cer = character_error_rate(stimulus,best_response)\n",
    "    cer_set.append(cer)\n",
    "    \n",
    "    print('S-R: %s-%s' % (stimulus,best_response))\n",
    "    print('CER: %.2f\\n' % (cer))\n",
    "    \n",
    "print('Mean CER across training set: %.2f' % (np.mean(cer_set)))\n",
    "\n",
    "# TODO: implement grid search to minimize CER\n",
    "\n",
    "min_mean_cer = 1.0\n",
    "best_dp = -1\n",
    "\n",
    "sigma_values = np.arange(0.1,1.0,0.1)\n",
    "for sigma in sigma_values:\n",
    "    def get_likelihood(x, y, key):\n",
    "        key_x = key.centre_x\n",
    "        key_y = key.centre_y\n",
    "        p = multivariate_normal.pdf([x,y], mean=[key_x,key_y], cov=[[sigma,0],[0,sigma]]);\n",
    "        return p\n",
    "    bw_values = np.arange(0.05,0.25,0.05)\n",
    "    for beam_width in bw_values:\n",
    "        ip_values = np.arange(1.0,1.5,0.25)\n",
    "        for insertion_penalty in ip_values:\n",
    "            dp_values = np.arange(0.005,0.02,0.005)\n",
    "            for deletion_penalty in dp_values:\n",
    "                # Keep track of the character error rates\n",
    "                cer_set = []\n",
    "\n",
    "                # Print out stimuli and observation sequences\n",
    "                stimulus_count = 0\n",
    "                for observation_sequence in observation_sequence_set:\n",
    "                    stimulus = stimulus_set[stimulus_count]\n",
    "                    stimulus_count += 1\n",
    "\n",
    "                    completed_tokens = []\n",
    "\n",
    "                    seed_token = Token(\"\", 1.0, -1)\n",
    "                    propagate(seed_token, observation_sequence, completed_tokens)\n",
    "                    results = sorted(completed_tokens, key=lambda token: token.acc_prob, reverse=True)\n",
    "\n",
    "                    # Get best hypothesis from the decoder result list                \n",
    "                    best_response = \"\"\n",
    "                    if (len(results) > 0):\n",
    "                        best_response = results[0].hypo\n",
    "\n",
    "                    # Compute character error rate\n",
    "                    cer = character_error_rate(stimulus,best_response)\n",
    "                    cer_set.append(cer)\n",
    "\n",
    "                # Compute mean CER across test set\n",
    "                mean_cer = np.mean(cer_set)\n",
    "                if (mean_cer < min_mean_cer):\n",
    "                    best_bw = beam_width\n",
    "                    best_dp = deletion_penalty\n",
    "                    best_ip = insertion_penalty\n",
    "                    best_sigma = sigma\n",
    "                    min_mean_cer = mean_cer\n",
    "\n",
    "print(min_mean_cer)\n",
    "print(best_sigma)\n",
    "print(best_dp)\n",
    "print(best_ip)\n",
    "print(best_bw)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\"> \n",
    "    <strong>Task 4:</strong> Evaluate your refined decoder on the test data generated with the same noisy toy keyboard (sigma=30). You will be judged based on your mean CER across the set of observation sequences. The file (test_data.log) will be committed to the repository at the end of the allocated time.  Ensure all your parameter values are configured as desired before executing. Do not modify the following code.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from parse_sr_log import ParseStimulusResponseLog\n",
    "\n",
    "# Extract a matched array of stimuli and observation sequences from the stimulus-response.log\n",
    "stimulus_set, observation_sequence_set = ParseStimulusResponseLog('test_data.log')\n",
    "\n",
    "# Default parameter values\n",
    "deletion_penalty = 0.02\n",
    "insertion_penalty = 1.2\n",
    "beam_width = 0.2\n",
    "sigma = 1.5\n",
    "\n",
    "def get_likelihood(x, y, key):\n",
    "    key_x = key.centre_x\n",
    "    key_y = key.centre_y\n",
    "    p = multivariate_normal.pdf([x,y], mean=[key_x,key_y], cov=[[sigma,0],[0,sigma]]);\n",
    "    return p\n",
    "\n",
    "# Keep track of the character error rates\n",
    "cer_set = []\n",
    "\n",
    "# Print out stimuli and observation sequences\n",
    "stimulus_count = 0\n",
    "for observation_sequence in observation_sequence_set:\n",
    "    stimulus = stimulus_set[stimulus_count]\n",
    "    stimulus_count += 1\n",
    "    \n",
    "    completed_tokens = []\n",
    "    \n",
    "    seed_token = Token(\"\", 1.0, -1)\n",
    "    propagate(seed_token, observation_sequence, completed_tokens)\n",
    "    results = sorted(completed_tokens, key=lambda token: token.acc_prob, reverse=True)\n",
    "\n",
    "    # Get best hypothesis from the decoder result list\n",
    "    best_response = results[0].hypo\n",
    "    \n",
    "    # Compute character error rate\n",
    "    cer = character_error_rate(stimulus,best_response)\n",
    "    cer_set.append(cer)\n",
    "    \n",
    "    print('S-R: %s-%s' % (stimulus,best_response))\n",
    "    print('CER: %.2f\\n' % (cer))\n",
    "    \n",
    "print('Mean CER across training set: %.4f' % (np.mean(cer_set)))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}