hmm/hmm_guide.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Implementation of HMM algorithms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Choosing the appropriate data structures\n",
    "You can choose any data structure to represent your HMM as long as your implementation works. The important thing to bear in mind is that your model needs to contain: \n",
    "* N states  \n",
    "  Represented by characters, strings or numbers. Can include the first and the last state explicitly but can also work without them. \n",
    "* A set of transition probabilities {a<sub>kj</sub>}  \n",
    "  a<sub>kj</sub> = P(π(i) = j | π(i-1) = k)  \n",
    "  The probability of being in the state j at step i given that at step i-1 we were in the state k.\n",
    "* A set of emission probabilities {e<sub>k</sub>(c)}  \n",
    "  e<sub>k</sub>(c) = P(s<sub>i</sub> = c | π(i) = k)  \n",
    "  The probability of observing the symbol c at the i-th position of the sequence given that at the i-th step we are in the state k.\n",
    "  \n",
    "One possible representation is described below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### States\n",
    "States are represented as a list of characters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "states = ['b','y','n','e']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transitions\n",
    "Transitions are represented as a dictionary where keys are tuples and values are transition probabilities. The first element of the tuple is the state from which we transition and the second element is the state we transition to."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "transitions = {('b','y') : 0.2,\n",
    "               ('b','n') : 0.8,\n",
    "               ('y','y') : 0.7,\n",
    "               ('y','n') : 0.2,\n",
    "               ('y','e') : 0.1,\n",
    "               ('n','n') : 0.8,\n",
    "               ('n','y') : 0.1,\n",
    "               ('n','e') : 0.1\n",
    "    }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Emissions\n",
    "Emissions are represented as a dictionary where keys are states and values are dictionaries. The internal dictionaries contain emitted symbols as keys and emission probabilities as values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "emissions = {'y' : {'A':0.1, 'C':0.4, 'G':0.4, 'T':0.1},\n",
    "             'n' : {'A':0.25, 'C':0.25, 'G':0.25, 'T':0.25}\n",
    "    }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sequence\n",
    "The sequence is simply a string."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "sequence = 'ATGCG'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initializing matrices\n",
    "We can write a simple utility function to initialize matrices that will be useful afterwards. It takes the desired number of rows, the desired number of columns and the value to fill the matrix with. If the third parameter is not provided, the matrix is filled with zeros by default."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def initialize_matrix(dim1,dim2,value=0):\n",
    "    F = []\n",
    "    for i in range(0,dim1):\n",
    "        F.append([])\n",
    "        for j in range(0,dim2):\n",
    "            F[i].append(value)\n",
    "    return F"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing matrices\n",
    "For visualization purposes only we implement some printing functions. They are not important for the sake of the solution and the details of the implementation can be ignored."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_matrix(matrix,axis1,axis2):\n",
    "    w = '{:<10}'\n",
    "    print w.format('') + w.format('0') + ''.join([w.format(char) for char in axis2]) + w.format('0')\n",
    "    for i, row in enumerate(matrix):\n",
    "        print w.format(axis1[i]) + ''.join(['{:<10.2e}'.format(item) for item in row])\n",
    "        \n",
    "def print_matrix_p(matrix,axis1,axis2):\n",
    "    w = '{:<10}'\n",
    "    print w.format('') + w.format('0') + ''.join([w.format(char) for char in axis2]) + w.format('0')\n",
    "    for i, row in enumerate(matrix):\n",
    "        print w.format(axis1[i]) + ''.join(['{:<10s}'.format(item) for item in row])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Forward algorithm\n",
    "We start by initializing matrix F that will contain as many rows as there are states and as many columns as there are symbols in the sequence plus two additional columns (the first and the last). The probability of starting from the begin state is 1 so we set the first element of the matrix to 1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          0         A         T         G         C         G         0         \n",
      "b         1.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n",
      "y         0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n",
      "n         0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n",
      "e         0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n"
     ]
    }
   ],
   "source": [
    "F = initialize_matrix(len(states),len(sequence)+2)\n",
    "F[0][0] = 1\n",
    "\n",
    "print_matrix(F,states,sequence)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we calculate the values for the first symbol. For each state, it is the probability of transitioning from the begin state to the current state times the probability of emitting the first symbol of the sequence in the current state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          0         A         T         G         C         G         0         \n",
      "b         1.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n",
      "y         0.00e+00  2.00e-02  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n",
      "n         0.00e+00  2.00e-01  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n",
      "e         0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n"
     ]
    }
   ],
   "source": [
    "for i in range(1,len(states)-1):\n",
    "        F[i][1] = transitions[(states[0],states[i])]*emissions[states[i]][sequence[0]]\n",
    "        \n",
    "print_matrix(F,states,sequence)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For all of the other symbols, from the second to the last, we calculate the values as the sum of probabilities. For each state, it is the sum of probabilities of transitioning from any previous state to the current state times the probability of emitting the corresponding symbol of the sequence in the current state times the probability of the previous state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          0         A         T         G         C         G         0         \n",
      "b         1.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n",
      "y         0.00e+00  2.00e-02  3.40e-03  2.59e-03  1.06e-03  3.69e-04  0.00e+00  \n",
      "n         0.00e+00  2.00e-01  4.10e-02  8.37e-03  1.80e-03  4.14e-04  0.00e+00  \n",
      "e         0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n"
     ]
    }
   ],
   "source": [
    "for j in range(2,len(sequence)+1): # loops on the symbols\n",
    "        for i in range(1,len(states)-1): # loops on the states\n",
    "            p_sum = 0\n",
    "            for k in range(1,len(states)-1): # loops on all of the possible previous states\n",
    "                p_sum += F[k][j-1]*transitions[(states[k],states[i])]*emissions[states[i]][sequence[j-1]]\n",
    "            F[i][j] = p_sum\n",
    "            \n",
    "print_matrix(F,states,sequence)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we need to calculate the final value. It is the sum of probabilities of transitioning from any previous state to the end state times the probability of the previous state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          0         A         T         G         C         G         0         \n",
      "b         1.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n",
      "y         0.00e+00  2.00e-02  3.40e-03  2.59e-03  1.06e-03  3.69e-04  0.00e+00  \n",
      "n         0.00e+00  2.00e-01  4.10e-02  8.37e-03  1.80e-03  4.14e-04  0.00e+00  \n",
      "e         0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  7.83e-05  \n"
     ]
    }
   ],
   "source": [
    "p_sum = 0\n",
    "for k in range(1,len(states)-1):\n",
    "    p_sum += F[k][len(sequence)]*transitions[(states[k],states[-1])]\n",
    "F[-1][-1] = p_sum\n",
    "\n",
    "print_matrix(F,states,sequence)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now put everything into a function that takes states, transitions, emissions and a sequence and returns a matrix F."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def forward(states,transitions,emissions,sequence):\n",
    "    F = initialize_matrix(len(states),len(sequence)+2)\n",
    "    F[0][0] = 1\n",
    "    for i in range(1,len(states)-1):\n",
    "        F[i][1] = transitions[(states[0],states[i])]*emissions[states[i]][sequence[0]]\n",
    "    for j in range(2,len(sequence)+1):\n",
    "        for i in range(1,len(states)-1):\n",
    "            p_sum = 0\n",
    "            for k in range(1,len(states)-1):\n",
    "                p_sum += F[k][j-1]*transitions[(states[k],states[i])]*emissions[states[i]][sequence[j-1]]\n",
    "            F[i][j] = p_sum\n",
    "    p_sum = 0\n",
    "    for k in range(1,len(states)-1):\n",
    "        p_sum += F[k][len(sequence)]*transitions[(states[k],states[-1])]\n",
    "    F[-1][-1] = p_sum\n",
    "    return F"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Backward algorithm\n",
    "We start by initializing matrix F that will contain as many rows as there are states and as many columns as there are symbols in the sequence plus two additional columns (the first and the last). The probability of starting from the end state is 1 so we set the last element of the matrix to 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          0         A         T         G         C         G         0         \n",
      "b         0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n",
      "y         0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n",
      "n         0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n",
      "e         0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  1.00e+00  \n"
     ]
    }
   ],
   "source": [
    "F = initialize_matrix(len(states),len(sequence)+2)\n",
    "F[-1][-1] = 1\n",
    "\n",
    "print_matrix(F,states,sequence)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we calculate the values for the last symbol. For each state, it is the probability of transitioning from the current state to the end state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          0         A         T         G         C         G         0         \n",
      "b         0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n",
      "y         0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  1.00e-01  0.00e+00  \n",
      "n         0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  1.00e-01  0.00e+00  \n",
      "e         0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  1.00e+00  \n"
     ]
    }
   ],
   "source": [
    "for i in range(1,len(states)-1):\n",
    "    F[i][-2] = transitions[(states[i],states[-1])]\n",
    "    \n",
    "print_matrix(F,states,sequence)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For all of the other symbols, from the second to last to the first, we calculate the values as the sum of probabilities. For each state, it is the sum of probabilities of transitioning from the current state to any successive state times the probability of emitting the corresponding symbol of the sequence in the successive state times the probability of the successive state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          0         A         T         G         C         G         0         \n",
      "b         0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n",
      "y         0.00e+00  3.08e-04  3.23e-03  1.04e-02  3.30e-02  1.00e-01  0.00e+00  \n",
      "n         0.00e+00  3.61e-04  1.64e-03  6.12e-03  2.40e-02  1.00e-01  0.00e+00  \n",
      "e         0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  1.00e+00  \n"
     ]
    }
   ],
   "source": [
    "for j in range(len(sequence)-1,0,-1): # loops on the symbols\n",
    "    for i in range(1,len(states)-1): # loops on the states\n",
    "        p_sum = 0\n",
    "        for k in range(1,len(states)-1): # loops on all of the possible successive states\n",
    "            p_sum += F[k][j+1]*transitions[(states[i],states[k])]*emissions[states[k]][sequence[j]]\n",
    "        F[i][j] = p_sum\n",
    "\n",
    "print_matrix(F,states,sequence)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we need to calculate the final value. It is the sum of probabilities of transitioning from the begin state to any successive state times the probability of emitting the first symbol of the sequence in the successive state times the probability of the successive state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          0         A         T         G         C         G         0         \n",
      "b         7.83e-05  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n",
      "y         0.00e+00  3.08e-04  3.23e-03  1.04e-02  3.30e-02  1.00e-01  0.00e+00  \n",
      "n         0.00e+00  3.61e-04  1.64e-03  6.12e-03  2.40e-02  1.00e-01  0.00e+00  \n",
      "e         0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  1.00e+00  \n"
     ]
    }
   ],
   "source": [
    "p_sum = 0\n",
    "for k in range(1,len(states)-1):\n",
    "    p_sum += F[k][1]*transitions[(states[0],states[k])]*emissions[states[k]][sequence[0]]\n",
    "F[0][0] = p_sum\n",
    "\n",
    "print_matrix(F,states,sequence)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now put everything into a function that takes states, transitions, emissions and a sequence and returns a matrix F."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def backward(states,transitions,emissions,sequence):\n",
    "    F = initialize_matrix(len(states),len(sequence)+2)\n",
    "    F[-1][-1] = 1\n",
    "    for i in range(1,len(states)-1):\n",
    "        F[i][-2] = transitions[(states[i],states[-1])]\n",
    "    for j in range(len(sequence)-1,0,-1): \n",
    "        for i in range(1,len(states)-1):\n",
    "            p_sum = 0\n",
    "            for k in range(1,len(states)-1):\n",
    "                p_sum += F[k][j+1]*transitions[(states[i],states[k])]*emissions[states[k]][sequence[j]]\n",
    "            F[i][j] = p_sum\n",
    "    p_sum = 0\n",
    "    for k in range(1,len(states)-1):\n",
    "        p_sum += F[k][1]*transitions[(states[0],states[k])]*emissions[states[k]][sequence[0]]\n",
    "    F[0][0] = p_sum\n",
    "    return F"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Viterbi algorithm\n",
    "We start by initializing matrices F and FP that will contain as many rows as there are states and as many columns as there are symbols in the sequence plus two additional columns (the first and the last). Matrix FP is filled with 'b' symbols that indicate the begin state. The probability of starting from the begin state is 1 so we set the first element of the matrix F to 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          0         A         T         G         C         G         0         \n",
      "b         1.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n",
      "y         0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n",
      "n         0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n",
      "e         0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n",
      "\n",
      "          0         A         T         G         C         G         0         \n",
      "b         b         b         b         b         b         b         b         \n",
      "y         b         b         b         b         b         b         b         \n",
      "n         b         b         b         b         b         b         b         \n",
      "e         b         b         b         b         b         b         b         \n"
     ]
    }
   ],
   "source": [
    "F = initialize_matrix(len(states),len(sequence)+2)\n",
    "FP = initialize_matrix(len(states),len(sequence)+2,states[0])\n",
    "F[0][0] = 1\n",
    "\n",
    "print_matrix(F,states,sequence)\n",
    "print\n",
    "print_matrix_p(FP,states,sequence)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we calculate the values for the first symbol. For each state, it is the probability of transitioning from the begin state to the current state times the probability of emitting the first symbol of the sequence in the current state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          0         A         T         G         C         G         0         \n",
      "b         1.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n",
      "y         0.00e+00  2.00e-02  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n",
      "n         0.00e+00  2.00e-01  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n",
      "e         0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n"
     ]
    }
   ],
   "source": [
    "for i in range(1,len(states)-1):\n",
    "    F[i][1] = transitions[(states[0],states[i])]*emissions[states[i]][sequence[0]]\n",
    "    \n",
    "print_matrix(F,states,sequence)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to choose the maximum probability at each of the following steps we will need to find the maximum in a list of values. It will also be useful to know the index of the state that has this maximum probability so that we can set the corresponding pointer. Therefore, we implement a simple function that, given a list of values, returns the maximum value and the index of the maximum value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_max_val_ind(values):\n",
    "    max_val = values[0]\n",
    "    max_ind = 0\n",
    "    for ind, val in enumerate(values):\n",
    "        if val>max_val:\n",
    "            max_val = val\n",
    "            max_ind = ind\n",
    "    return max_val, max_ind"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For all of the other symbols, from the second to the last, we calculate the values as the maximum of probabilities. For each state, it is the maximum of probabilities of transitioning from any previous state to the current state times the probability of emitting the corresponding symbol of the sequence in the current state times the probability of the previous state. We also set the corresponding pointers accordingly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          0         A         T         G         C         G         0         \n",
      "b         1.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n",
      "y         0.00e+00  2.00e-02  2.00e-03  1.60e-03  4.48e-04  1.25e-04  0.00e+00  \n",
      "n         0.00e+00  2.00e-01  4.00e-02  8.00e-03  1.60e-03  3.20e-04  0.00e+00  \n",
      "e         0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n",
      "\n",
      "          0         A         T         G         C         G         0         \n",
      "b         b         b         b         b         b         b         b         \n",
      "y         b         b         n         n         y         y         b         \n",
      "n         b         b         n         n         n         n         b         \n",
      "e         b         b         b         b         b         b         b         \n"
     ]
    }
   ],
   "source": [
    "for j in range(2,len(sequence)+1): # loops on the symbols\n",
    "        for i in range(1,len(states)-1): # loops on the states\n",
    "            values = []\n",
    "            for k in range(1,len(states)-1): # loops on all of the possible previous states\n",
    "                values.append(F[k][j-1]*transitions[(states[k],states[i])]*emissions[states[i]][sequence[j-1]]) # appends the value to a list\n",
    "            max_val, max_ind = get_max_val_ind(values) # finds the maximum and the index of the maximum in the list\n",
    "            F[i][j] = max_val # sets the probability to the maximum probability\n",
    "            FP[i][j] = states[max_ind+1] # sets the corresponding pointer to the appropriate previous state\n",
    "\n",
    "print_matrix(F,states,sequence)\n",
    "print\n",
    "print_matrix_p(FP,states,sequence)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we need to calculate the final value. It is the maximum of probabilities of transitioning from any previous state to the end state times the probability of the previous state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          0         A         T         G         C         G         0         \n",
      "b         1.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  \n",
      "y         0.00e+00  2.00e-02  2.00e-03  1.60e-03  4.48e-04  1.25e-04  0.00e+00  \n",
      "n         0.00e+00  2.00e-01  4.00e-02  8.00e-03  1.60e-03  3.20e-04  0.00e+00  \n",
      "e         0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  3.20e-05  \n",
      "\n",
      "          0         A         T         G         C         G         0         \n",
      "b         b         b         b         b         b         b         b         \n",
      "y         b         b         n         n         y         y         b         \n",
      "n         b         b         n         n         n         n         b         \n",
      "e         b         b         b         b         b         b         n         \n"
     ]
    }
   ],
   "source": [
    "values = []\n",
    "for k in range(1,len(states)-1):\n",
    "    values.append(F[k][len(sequence)]*transitions[(states[k],states[-1])])\n",
    "max_val, max_ind = get_max_val_ind(values)\n",
    "F[-1][-1] = max_val\n",
    "FP[-1][-1] = states[max_ind+1]\n",
    "\n",
    "print_matrix(F,states,sequence)\n",
    "print\n",
    "print_matrix_p(FP,states,sequence)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now put everything into a function that takes states, transitions, emissions and a sequence and returns a matrix F and a matrix FP of pointers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def viterbi(states,transitions,emissions,sequence):\n",
    "    F = initialize_matrix(len(states),len(sequence)+2)\n",
    "    FP = initialize_matrix(len(states),len(sequence)+2,states[0])\n",
    "    F[0][0] = 1\n",
    "    for i in range(1,len(states)-1):\n",
    "        F[i][1] = transitions[(states[0],states[i])]*emissions[states[i]][sequence[0]]\n",
    "    for j in range(2,len(sequence)+1):\n",
    "        for i in range(1,len(states)-1):\n",
    "            values = []\n",
    "            for k in range(1,len(states)-1):\n",
    "                values.append(F[k][j-1]*transitions[(states[k],states[i])]*emissions[states[i]][sequence[j-1]])\n",
    "            max_val, max_ind = get_max_val_ind(values)\n",
    "            F[i][j] = max_val\n",
    "            FP[i][j] = states[max_ind+1]\n",
    "    values = []\n",
    "    for k in range(1,len(states)-1):\n",
    "        values.append(F[k][len(sequence)]*transitions[(states[k],states[-1])])\n",
    "    max_val, max_ind = get_max_val_ind(values)\n",
    "    F[-1][-1] = max_val\n",
    "    FP[-1][-1] = states[max_ind+1]\n",
    "    return F, FP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to traceback the path we need to start from the end state and move through the FP matrix according to the pointers. Therefore, we implement a function that takes states and the FP matrix, starts from the last element of the matrix and moves through the matrix appending, at each step, the current state to the path variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def traceback(states,FP):\n",
    "    path = ['e'] # the last element of the path is the end state\n",
    "    current = FP[-1][-1] # the current state is the one written in the last cell of the matrix\n",
    "    for i in range(len(FP[0])-2,0,-1): # loops on the symbols\n",
    "        path = [current] + path # appends the current state to the path\n",
    "        current = FP[states.index(current)][i] # finds the index of the current state in the list of states and moves to the corresponing row of FP \n",
    "    path = ['b'] + path # the first element of the path is the begin state\n",
    "    return ' '.join(path) # transforms the list into a string where elements are separated by spaces"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now visualize the most likely path we obtained."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "- A T G C G -\n",
      "b n n n n n e\n"
     ]
    }
   ],
   "source": [
    "path = traceback(states,FP)\n",
    "print '- '+' '.join(sequence)+' -'\n",
    "print path"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}