{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "def traceLcs(D, x, y):\n",
    "    ''' Backtrace for LCS; returns LCS as string. '''\n",
    "    i, j = len(x), len(y) # start in lower-right\n",
    "    st = []\n",
    "    while i > 0 and j > 0:\n",
    "        # get the three contributions\n",
    "        diag, vert, horz = 0, 0, 0\n",
    "        if i > 0 and j > 0:\n",
    "            delt = -1 if x[i-1] == y[j-1] else 1\n",
    "            diag = D[i-1, j-1] + delt\n",
    "        if i > 0: vert = D[i-1, j]\n",
    "        if j > 0: horz = D[i, j-1]\n",
    "        if diag <= vert and diag <= horz:\n",
    "            # diagonal is best, thus, this char is part of LCS\n",
    "            st.append(x[i-1])\n",
    "            i -= 1; j -= 1 # move up and left\n",
    "        elif vert <= horz: i -= 1 # vertical is best; move up\n",
    "        else: j -= 1 # horizontal is best; move left\n",
    "    # reverse it, then return string-ized LCS\n",
    "    return (''.join(st))[::-1]\n",
    "\n",
    "import numpy\n",
    "def lcsDp(x, y):\n",
    "    ''' Return LCS of x and y.  Uses bottom-up dynamic programming\n",
    "        approach. '''\n",
    "    D = numpy.zeros((len(x)+1, len(y)+1), dtype=int)\n",
    "    for i in range(1, len(x)+1):\n",
    "        for j in range(1, len(y)+1):\n",
    "            delt = -1 if x[i-1] == y[j-1] else 1\n",
    "            D[i, j] = min(D[i-1, j-1] + delt, D[i-1, j], D[i, j-1])\n",
    "    return traceLcs(D, x, y), D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'TCTA'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lcs, D = lcsDp('ATCTGAT', 'TGCATA') # example from Jones & Pevzner 6.5\n",
    "lcs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  0,  0,  0,  0,  0,  0],\n",
       "       [ 0,  0,  0,  0, -1, -1, -1],\n",
       "       [ 0, -1, -1, -1, -1, -2, -2],\n",
       "       [ 0, -1, -1, -2, -2, -2, -2],\n",
       "       [ 0, -1, -1, -2, -2, -3, -3],\n",
       "       [ 0, -1, -2, -2, -2, -3, -3],\n",
       "       [ 0, -1, -2, -2, -3, -3, -4],\n",
       "       [ 0, -1, -2, -2, -3, -4, -4]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D"
   ]
  }
 ],
 "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.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}