BA10J - Forward-Backward - Soft Decoding problem.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## HMM beyond Viterbi learning:\n",
    "\n",
    "https://www.youtube.com/watch?v=yUZ8CBdeJRs&list=PLQ-85lQlPqFPnk31Uut2ajVkBvlFmMtdx&index=9\n",
    "\n",
    "Also, because Viterbi learning is dependent on the initial guess for Parameters, it may become stuck in a local optimum. Like other heuristics, it is often run many times, retaining the best choice of Parameters.\n",
    "\n",
    "**Exercise Break**: Apply Viterbi learning to learn parameters for an HMM modeling CG-islands as well as for the profile HMM for the gp120 HIV alignment.\n",
    "\n",
    "---\n",
    "\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Recall the outcome likelihood problem: \n",
    "\n",
    "Given arbitrary HMM and emission string $x$, we calculate the sum of the probabilities of the hidden paths entering a node. We progressively fill in the viterbi matrix from left to right (hence *forward algorithm*), using a similar reccurence to the Viterbi algorithm, except in the Forward algo we take the sum of all the incoming edges instead of the max.\n",
    "\n"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "def outcome_likelihood(x, T, E, states): #outcome-likelihood\n",
    "    \"\"\"returns likelihood of emission string, given HMM\"\"\"\n",
    "\n",
    "    S = len(states); n = len(x)\n",
    "    # viterbi graph |states|*|time of emission| matrix\n",
    "    viterbi = np.zeros(shape = (S, n)) \n",
    "\n",
    "    # init first column of viterbi with Pr_emission & 1/States\n",
    "    for state in range(S):\n",
    "        viterbi[state][0] = 1/S * E[state][emission[0]]\n",
    "\n",
    "    # Fill viterbi: sum incoming edges probabilities for each node viterbi[state, emission]\n",
    "    # probability for an edge is the product of p(state -> emission), p(transition between states), and p(hidden path to parent node leading to transition)\n",
    "    for i in range(1,n):\n",
    "        for state in range(S):\n",
    "            em = E[state][emission[i]]\n",
    "            for prev in range(S):\n",
    "                trans = T[prev][state]\n",
    "                viterbi[state][i] += trans * em * viterbi[prev][i-1] \n",
    "    return sum(viterbi[s][n-1] for s in range(S))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "---\n",
    "\n",
    "<br>\n",
    "\n",
    "### Soft Decoding Problem:\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "The Soft Decoding Problem\n",
    "In a previous chapter, we introduced a “soft” clustering algorithm, based on the more general expectation maximization algorithm, that relaxed the Lloyd algorithm’s rigid assignment of points to clusters. Analogously, by generating a single optimal hidden path, the Viterbi algorithm provides a rigid “yes” or “no” answer to the question of whether an HMM was in state k at time i. But how certain are we that this was the case?\n",
    "\n",
    "Returning to the crooked casino analogy once more, say that the i-th coin flip is heads. If this flip occurs in the middle of ten consecutive heads, then you should be relatively confident that the biased coin was used. But what if, of the ten flips surrounding the i-th flip, six are heads and four are tails? In this case, you should be less certain that the biased coin was used.\n",
    "\n",
    "In the case of an **arbitrary HMM, we would like to compute the conditional probability Pr(πi = k|x) that the HMM was in state k at time i given that it emitted string x**.\n",
    "\n",
    "\n",
    "\n",
    "Find the **probability** that an HMM was in a **particular state at a particular moment** given its emitted string.\n",
    "\n",
    "**Input**:\n",
    "   -  A string x = x1 ... xn emitted by an HMM.  \n",
    "   \n",
    "**Output**:\n",
    "   - The **conditional probability** Pr(πi = k|x) that the **HMM was in state k at step i** given that it emitted x.\n",
    "\n",
    "<br>\n",
    "\n",
    "### *another learning problem, this time with soft decisions.*\n",
    "\n",
    "*in the first Viterbi learning algorithm, we assigned a state for each position in the hidden path; instead we should assign a weight to each state for each position in the path. This allows for more probabilistic decisions instead of a hard assignment scenario.*\n",
    "\n",
    "*The idea is analogous to Lloyd algorithm vs soft k-means clustering from gene expression chapter, except asking for a specific outcome's probability*\n",
    "\n",
    "---\n",
    "\n",
    "<br>\n",
    "\n",
    "The **unconditional probability that a hidden path will pass through state $k$ at time $i$ and emit $x$** can be written as the **sum**:\n",
    "\n",
    "<br>\n",
    "\n",
    "\\begin{aligned} \\mathrm{Pr}(\\pi_i=k, x) = \\displaystyle\\sum_{\\text{all paths }\\pi\\text{ with }\\pi_i =k}\\mathrm{Pr}(x, \\pi) \\end{aligned}\n",
    "\n",
    "\n",
    "<br> \n",
    "\n",
    "##### The **conditional probability $Pr(π_i = k|x)$** is equal to the **proportion of paths that pass through** state k at time i and emit x with respect to all paths emitting x:\n",
    "\n",
    "\\begin{aligned} \\mathrm{Pr}(\\pi_i =k|x) & = \\dfrac{\\mathrm{Pr}(\\pi_i=k, x)}{\\mathrm{Pr}(x)}\\\\ & = \\dfrac{\\sum_{\\text{all paths }\\pi\\text{ with }\\pi_i =k}\\mathrm{Pr}(x, \\pi)}{\\sum_{\\text{all paths } \\pi} \\mathrm{Pr}(x, \\pi)}\\,. \\end{aligned} \n",
    "\n",
    "<br>\n",
    "\n",
    "**STOP and Think**: If the Viterbi algorithm for the crooked casino emits a path $π = π_1, π_2, ...π_n$ with $π_i = B$, is the dealer more likely to have used a biased coin at step $i$? \n",
    "\n",
    "*not necessarily, Viterbi algorithm finds the optimal path through the entire graph; it could be that another state was more likely at this particular time, but the paths passing through that state were not the optimal path*\n",
    "\n",
    "Is it possible that $π_i = B$ but that $Pr(π_i = B|x)$ is smaller than $Pr(π_i = F|x)$? (*yes, perhaps less likely*)\n",
    "\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "#### The forward-backward algorithm\n",
    "\n",
    "\n",
    "\n",
    "We note that $Pr(π_i = k, x)$ is equal to the sum of product weights Pr(π, x) of all paths π through the Viterbi graph for x that pass through the node $(k, i)$. As shown in the figure below, we can break each such path into a blue subpath from source to (k, i), which we denote $π$, and a subpath from $(k, i)$ to $sink$, which we denote $π$. Writing $Weight(π_{blue})$ and $Weight(π_{red})$ as the respective product weights of these subpaths yields the recurrence:\n",
    "\n",
    "\n",
    "\\begin{aligned} \\mathrm{Pr}(\\pi_i = k, x) & = \\sum_{\\text{all paths }\\pi \\text{ with }\\pi_i =k} \\mathrm{Pr}(x, \\pi)\\\\ & = \\sum_{\\color{blue}{\\text{all paths }\\pi_\\text{blue}}} \\sum_{\\color{red}{\\text{all paths }\\pi_\\text{red}}} \n",
    "\\textit{Weight}(\\pi_\\text{blue}) \\cdot \\textit{Weight}(\\pi_\\text{red})\\\\ & = \\sum_{\\color{blue}{\\text{all paths }\\pi_\\text{blue}}} \\textit{Weight}(\\pi_\\text{blue}) \\cdot \\sum_{\\color{red}{\\text{all paths }\\pi_\\text{red}}} \\textit{Weight}(\\pi_\\text{red}) \\end{aligned}\n",
    "\n",
    "<br>\n",
    "\n",
    "#### The Forward algorithm solves the Outcome Likelihood Problem & is similar to the recurrence in the Viterbi algorithm\n",
    "We have already computed the **sum of product weights of all blue subpaths**; it is just $forward(k)$, i , which we encountered when solving the **Outcome Likelihood Problem**. Now we would like to compute the sum of product weights of all red subpaths, which we denote as backwardk, i , so that the preceding equation becomes\n",
    "\n",
    "\n",
    "\\begin{aligned} \\mathrm{Pr}(\\pi_i=k, x) = {\\color{blue}{\\textit{forward}_{k,i}}} \\cdot {\\color{red}{\\textit{backward}_{k,i}}} \\end{aligned}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Backward algorithm\n",
    "The name of $backward(k, i)$ derives from the fact that to compute this value, we can **simply reverse the directions of all edges in the Viterbi graph** (see figure below) and apply the same dynamic programming algorithm used to compute forwardk, i . Since the reversed edge connecting (l, i + 1) to (k, i) has weight $Weight_i(k, l) = transition(k, l) \\cdot emission_l(x_{i+1})$, we have that\n",
    "\n",
    "<br>\n",
    "\n",
    "\\begin{align} {\\color{red}{\\textit{backward}_{k,i}}} = \\displaystyle\\sum_{\\text{all states }l} {\\color{red}{\\textit{backward}_{l,i+1}}} \\cdot \\textit{Weight}_i(k,l) \\end{align}\n",
    "\n",
    "<br>\n",
    "STOP and Think: How should this recurrence be initialized?\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### *Initializing the Backward algorithm* (I think...)\n",
    "\n",
    "When we entered the first column of |states| for the first position in $emissions$ in the Forward algorithm, we simply gave the transition probabilities as being equal for each state: $T(source \\rightarrow \\pi_0) = \\dfrac{1}{|states|}$\n",
    "\n",
    "<br> \n",
    "\n",
    "For the Backward algorithm need to consider the edges from nodes in the last column of the Viterbi graph $(state, n)$ that always lead to the sink.\n",
    "\n",
    "So instead of Transition probability = 1/|states| it should be just 1 since there is only one edge from each final state to the sink, so the transition prob is 1.\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def forward(X, T, E, S): \n",
    "    \"\"\"returns likelihood of emission string, given HMM\"\"\"\n",
    "    # F: |states|*|time of emission| matrix\n",
    "    # first column of F with P(state,emission) * (1/States; naive probability initials)\n",
    "    F = np.zeros(shape = (S, len(X))) \n",
    "    for state in range(S):\n",
    "        F[state][0] = E[state][X[0]] / S\n",
    "    # Fill forward matrix: p(node) =  p(emission) * sum(p(transition) * p(previous node))\n",
    "    for i in range(1,len(X)):\n",
    "        for state in range(S):\n",
    "            F[state][i] = sum(T[k][state]*F[k][i-1] for k in range(S))\n",
    "            F[state][i] *= E[state][X[i]]\n",
    "    return F\n",
    "\n",
    "# backwards outcome likelihood\n",
    "def backward(X, T, E, S): \n",
    "    B = np.zeros(shape = (S, len(X)))\n",
    "    # last col: p(trans(sink,last state)) is p=1 instead of 1/S as in forward algo??\n",
    "    for state in range(S):\n",
    "        B[state][len(X)-1] =  1\n",
    "    for i in range(len(X)-2, -1, -1):\n",
    "        for state in range(S):\n",
    "            B[state][i] = sum(T[state][k] * B[k][i+1] * E[k][X[i+1]] for k in range(S))\n",
    "    return B\n",
    "\n",
    "def forward_backward(X,T,E,S):\n",
    "    \"\"\" return responsibility matrix (smoothing)\"\"\"\n",
    "    F = forward(X, T, E, S)\n",
    "    B = backward(X, T, E, S)\n",
    "    F_sink = sum(F[state][len(X)-1] for state in range(S))\n",
    "    return np.multiply(F, B)/F_sink\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/jasonmoggridge/Dropbox/Rosalind/Coursera_textbook_track/Course6/data\n",
      "/Users/jasonmoggridge/Dropbox/Rosalind/Coursera_textbook_track/Course6/results\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>A</th>\n",
       "      <th>B</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.5438</td>\n",
       "      <td>0.4562</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.6492</td>\n",
       "      <td>0.3508</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.9647</td>\n",
       "      <td>0.0353</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.9936</td>\n",
       "      <td>0.0064</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.9957</td>\n",
       "      <td>0.0043</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.9891</td>\n",
       "      <td>0.0109</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.9154</td>\n",
       "      <td>0.0846</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0.9640</td>\n",
       "      <td>0.0360</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>0.8737</td>\n",
       "      <td>0.1263</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0.8167</td>\n",
       "      <td>0.1833</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        A       B\n",
       "0  0.5438  0.4562\n",
       "1  0.6492  0.3508\n",
       "2  0.9647  0.0353\n",
       "3  0.9936  0.0064\n",
       "4  0.9957  0.0043\n",
       "5  0.9891  0.0109\n",
       "6  0.9154  0.0846\n",
       "7  0.9640  0.0360\n",
       "8  0.8737  0.1263\n",
       "9  0.8167  0.1833"
      ]
     },
     "execution_count": 92,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "%cd /Users/jasonmoggridge/Dropbox/Rosalind/Coursera_textbook_track/Course6/data\n",
    "with open(\"ba10j.sample.txt\",'r') as file:\n",
    "    inputs = file.readlines()\n",
    "    alpha = inputs[2].strip().split()\n",
    "    X = list(alpha.index(x) for x in inputs[0].strip())\n",
    "    states = inputs[4].strip().split()\n",
    "    S = len(states)\n",
    "    # HMM Transition, emission matrices,\n",
    "    T = np.array([line.split()[1:] for line in inputs[7: 7+S]], float)\n",
    "    E = np.array([line.split()[1:] for line in inputs[7+S+2:]], float)\n",
    "    del(inputs)\n",
    "\n",
    "##\n",
    "R = np.transpose(forward_backward(X,T,E,S))\n",
    "R = np.around(R, 4)\n",
    "R = pd.DataFrame(R, columns=states)\n",
    "%cd /Users/jasonmoggridge/Dropbox/Rosalind/Coursera_textbook_track/Course6/results\n",
    "\n",
    "R.to_csv(\"ba10j.sample.out.txt\", index=False)\n",
    "R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/jasonmoggridge/Dropbox/Rosalind/Coursera_textbook_track/Course6/data\n",
      "/Users/jasonmoggridge/Dropbox/Rosalind/Coursera_textbook_track/Course6/results\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>A</th>\n",
       "      <th>B</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.6283</td>\n",
       "      <td>0.3717</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.8035</td>\n",
       "      <td>0.1965</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.8047</td>\n",
       "      <td>0.1953</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.6160</td>\n",
       "      <td>0.3840</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.6458</td>\n",
       "      <td>0.3542</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.6136</td>\n",
       "      <td>0.3864</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.8228</td>\n",
       "      <td>0.1772</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0.5910</td>\n",
       "      <td>0.4090</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>0.8075</td>\n",
       "      <td>0.1925</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0.7976</td>\n",
       "      <td>0.2024</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        A       B\n",
       "0  0.6283  0.3717\n",
       "1  0.8035  0.1965\n",
       "2  0.8047  0.1953\n",
       "3  0.6160  0.3840\n",
       "4  0.6458  0.3542\n",
       "5  0.6136  0.3864\n",
       "6  0.8228  0.1772\n",
       "7  0.5910  0.4090\n",
       "8  0.8075  0.1925\n",
       "9  0.7976  0.2024"
      ]
     },
     "execution_count": 95,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "%cd /Users/jasonmoggridge/Dropbox/Rosalind/Coursera_textbook_track/Course6/data\n",
    "with open(\"dataset_26261_5.txt\",'r') as file:\n",
    "    inputs = file.readlines()\n",
    "    alpha = inputs[2].strip().split()\n",
    "    X = list(alpha.index(x) for x in inputs[0].strip())\n",
    "    states = inputs[4].strip().split()\n",
    "    S = len(states)\n",
    "    # HMM Transition, emission matrices,\n",
    "    T = np.array([line.split()[1:] for line in inputs[7: 7+S]], float)\n",
    "    E = np.array([line.split()[1:] for line in inputs[7+S+2:]], float)\n",
    "    del(inputs)\n",
    "\n",
    "##\n",
    "\n",
    "\n",
    "R = np.transpose(forward_backward(X,T,E,S))\n",
    "R = np.around(R, 4)\n",
    "\n",
    "%cd /Users/jasonmoggridge/Dropbox/Rosalind/Coursera_textbook_track/Course6/results\n",
    "R = pd.DataFrame(R, columns=states)\n",
    "R.to_csv(\"ba10j.problem.out.txt\", index=False, sep='\\t')\n",
    "R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/jasonmoggridge/Dropbox/Rosalind/Coursera_textbook_track/Course6/data\n",
      "/Users/jasonmoggridge/Dropbox/Rosalind/Coursera_textbook_track/Course6/results\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>A</th>\n",
       "      <th>B</th>\n",
       "      <th>C</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.4851</td>\n",
       "      <td>0.1737</td>\n",
       "      <td>0.3412</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.5028</td>\n",
       "      <td>0.1895</td>\n",
       "      <td>0.3077</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.3179</td>\n",
       "      <td>0.4505</td>\n",
       "      <td>0.2316</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.3489</td>\n",
       "      <td>0.4598</td>\n",
       "      <td>0.1913</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.3605</td>\n",
       "      <td>0.3580</td>\n",
       "      <td>0.2815</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.4952</td>\n",
       "      <td>0.1946</td>\n",
       "      <td>0.3102</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.3631</td>\n",
       "      <td>0.4414</td>\n",
       "      <td>0.1955</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0.5272</td>\n",
       "      <td>0.1792</td>\n",
       "      <td>0.2936</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>0.3471</td>\n",
       "      <td>0.3502</td>\n",
       "      <td>0.3027</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0.5379</td>\n",
       "      <td>0.1915</td>\n",
       "      <td>0.2706</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        A       B       C\n",
       "0  0.4851  0.1737  0.3412\n",
       "1  0.5028  0.1895  0.3077\n",
       "2  0.3179  0.4505  0.2316\n",
       "3  0.3489  0.4598  0.1913\n",
       "4  0.3605  0.3580  0.2815\n",
       "5  0.4952  0.1946  0.3102\n",
       "6  0.3631  0.4414  0.1955\n",
       "7  0.5272  0.1792  0.2936\n",
       "8  0.3471  0.3502  0.3027\n",
       "9  0.5379  0.1915  0.2706"
      ]
     },
     "execution_count": 96,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# rosalind data rosalind_ba10j\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "%cd /Users/jasonmoggridge/Dropbox/Rosalind/Coursera_textbook_track/Course6/data\n",
    "with open(\"rosalind_ba10j.txt\",'r') as file:\n",
    "    inputs = file.readlines()\n",
    "    alpha = inputs[2].strip().split()\n",
    "    X = list(alpha.index(x) for x in inputs[0].strip())\n",
    "    states = inputs[4].strip().split()\n",
    "    S = len(states)\n",
    "    # HMM Transition, emission matrices,\n",
    "    T = np.array([line.split()[1:] for line in inputs[7: 7+S]], float)\n",
    "    E = np.array([line.split()[1:] for line in inputs[7+S+2:]], float)\n",
    "    del(inputs)\n",
    "\n",
    "##\n",
    "\n",
    "\n",
    "R = np.transpose(forward_backward(X,T,E,S))\n",
    "R = np.around(R, 4)\n",
    "\n",
    "%cd /Users/jasonmoggridge/Dropbox/Rosalind/Coursera_textbook_track/Course6/results\n",
    "R = pd.DataFrame(R, columns=states)\n",
    "R.to_csv(\"ba10j.rosalind.out.txt\", index=False, sep='\\t')\n",
    "R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Finished Forward-Backward algorithm for soft-decoding problem\n",
    "# matrix answers what was likelihood of state,i given Xi\n",
    "# now hidden path is probabilistic in nature instead of hard assignment\n",
    "\n",
    "def forward_backward(X,T,E,S):\n",
    "    \n",
    "    \n",
    "    def forward(X, T, E, S): \n",
    "        \"\"\"returns likelihood of emission string, given HMM\"\"\"\n",
    "        # F: |states|*|time of emission| matrix\n",
    "        # first column of F with P(state,emission) * (1/States; naive probability initials)\n",
    "        F = np.zeros(shape = (S, len(X))) \n",
    "        for state in range(S):\n",
    "            F[state][0] = E[state][X[0]] / S\n",
    "        # Fill forward matrix: p(node) =  p(emission) * sum(p(transition) * p(previous node))\n",
    "        for i in range(1,len(X)):\n",
    "            for state in range(S):\n",
    "                F[state][i] = sum(T[k][state]*F[k][i-1] for k in range(S))\n",
    "                F[state][i] *= E[state][X[i]]\n",
    "        return F\n",
    "\n",
    "    # backwards outcome likelihood\n",
    "    def backward(X, T, E, S): \n",
    "        B = np.zeros(shape = (S, len(X)))\n",
    "        # last col: p(trans(sink,last state)) is p=1 instead of 1/S as in forward algo??\n",
    "        for state in range(S):\n",
    "            B[state][len(X)-1] =  1\n",
    "        for i in range(len(X)-2, -1, -1):\n",
    "            for state in range(S):\n",
    "                B[state][i] = sum(T[state][k] * B[k][i+1] * E[k][X[i+1]] for k in range(S))\n",
    "        return B\n",
    "    \n",
    "    \n",
    "    \"\"\" return responsibility matrix (smoothing)\"\"\"\n",
    "    F = forward(X, T, E, S)\n",
    "    B = backward(X, T, E, S)\n",
    "    F_sink = sum(F[state][len(X)-1] for state in range(S))\n",
    "    return np.multiply(F, B)/F_sink\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/jasonmoggridge/Dropbox/Rosalind/Coursera_textbook_track/Course6/data\n",
      "/Users/jasonmoggridge/Dropbox/Rosalind/Coursera_textbook_track/Course6/results\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>A</th>\n",
       "      <th>B</th>\n",
       "      <th>C</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.4851</td>\n",
       "      <td>0.1737</td>\n",
       "      <td>0.3412</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.5028</td>\n",
       "      <td>0.1895</td>\n",
       "      <td>0.3077</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.3179</td>\n",
       "      <td>0.4505</td>\n",
       "      <td>0.2316</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.3489</td>\n",
       "      <td>0.4598</td>\n",
       "      <td>0.1913</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.3605</td>\n",
       "      <td>0.3580</td>\n",
       "      <td>0.2815</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.4952</td>\n",
       "      <td>0.1946</td>\n",
       "      <td>0.3102</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.3631</td>\n",
       "      <td>0.4414</td>\n",
       "      <td>0.1955</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0.5272</td>\n",
       "      <td>0.1792</td>\n",
       "      <td>0.2936</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>0.3471</td>\n",
       "      <td>0.3502</td>\n",
       "      <td>0.3027</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0.5379</td>\n",
       "      <td>0.1915</td>\n",
       "      <td>0.2706</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        A       B       C\n",
       "0  0.4851  0.1737  0.3412\n",
       "1  0.5028  0.1895  0.3077\n",
       "2  0.3179  0.4505  0.2316\n",
       "3  0.3489  0.4598  0.1913\n",
       "4  0.3605  0.3580  0.2815\n",
       "5  0.4952  0.1946  0.3102\n",
       "6  0.3631  0.4414  0.1955\n",
       "7  0.5272  0.1792  0.2936\n",
       "8  0.3471  0.3502  0.3027\n",
       "9  0.5379  0.1915  0.2706"
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# rosalind data rosalind_ba10j\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "%cd /Users/jasonmoggridge/Dropbox/Rosalind/Coursera_textbook_track/Course6/data\n",
    "with open(\"rosalind_ba10j.txt\",'r') as file:\n",
    "    inputs = file.readlines()\n",
    "    alpha = inputs[2].strip().split()\n",
    "    X = list(alpha.index(x) for x in inputs[0].strip())\n",
    "    states = inputs[4].strip().split()\n",
    "    S = len(states)\n",
    "    # HMM Transition, emission matrices,\n",
    "    T = np.array([line.split()[1:] for line in inputs[7: 7+S]], float)\n",
    "    E = np.array([line.split()[1:] for line in inputs[7+S+2:]], float)\n",
    "    del(inputs)\n",
    "\n",
    "##\n",
    "\n",
    "\n",
    "R = np.transpose(forward_backward(X,T,E,S))\n",
    "R = np.around(R, 4)\n",
    "\n",
    "%cd /Users/jasonmoggridge/Dropbox/Rosalind/Coursera_textbook_track/Course6/results\n",
    "R = pd.DataFrame(R, columns=states)\n",
    "R.to_csv(\"ba10j.rosalind.out.txt\", index=False, sep='\\t')\n",
    "R"
   ]
  }
 ],
 "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}