CTC_Practice_Answer.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "DvgJT9G6dUrM"
   },
   "source": [
    "# CTC Forward, Backward and Gradient Practice. Compared with tf.ctc_loss.\n",
    "\n",
    "Credit 是此篇 [DingKe ipynb](https://github.com/DingKe/ml-tutorial/blob/master/ctc/CTC.ipynb) 的, 他完整呈現了 CTC loss 以及 gradient 的計算, 非常棒!\n",
    "\n",
    "此筆記加入自己的說明, 並且最後使用 tensorflow 來驗證.\n",
    "\n",
    "這篇另一個主要目的為改成可以練習的格式 (**#TODO tag**). 因為我相信最好的學習方式是自己造一次輪子. (這篇是參考答案)\n",
    "\n",
    "我們只專注在 CTC loss 的 forward, backwark and gradient. Decoding 部分請參考原作者的 [ipynb](https://github.com/DingKe/ml-tutorial/blob/master/ctc/CTC.ipynb). 最後使用 tf.nn.ctc_loss and tf.gradients 與我們的計算做對比"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "collapsed": true,
    "id": "J45YJwIgXnJ5"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import tensorflow as tf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "lHrFPI-FenD9"
   },
   "source": [
    "## 一些變數定義\n",
    "\n",
    "1. `Vocab = [0, 1, 2, 3, 4]`, 其中 **0 表示 'blank'**.\n",
    "\n",
    "2. `V = len(Vocab) = 5` 是字典大小, 譬如字典有 4 個 labels 加上一個 blank, 因此 V=5\n",
    "\n",
    "3. `l` 是正確答案 (已經安插 blanks 了), 例如 `l = [0, 3, 0, 3, 0, 4, 0]`\n",
    "\n",
    "4. `L = len(l) = 7`\n",
    "\n",
    "5. **後驗概率 `y`** (shape=`[V,T]`), 其中 `T` 表示 input sequence 長度, 所以 `y[k,t]` 表示時間點 `t`, label `k` 的後驗概率\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "collapsed": true,
    "id": "ZoDrr_F1d7mp"
   },
   "outputs": [],
   "source": [
    "Vocab = [0,1,2,3,4]\n",
    "l = [0, 3, 0, 3, 0, 4, 0]\n",
    "V, L = len(Vocab), len(l)\n",
    "T = 12\n",
    "logits = np.random.random([V,T])\n",
    "\n",
    "def softmax(logits):\n",
    "    max_value = np.max(logits, axis=0, keepdims=True)\n",
    "    exp = np.exp(logits - max_value)\n",
    "    exp_sum = np.sum(exp, axis=0, keepdims=True)\n",
    "    dist = exp / exp_sum\n",
    "    return dist\n",
    "  \n",
    "y = softmax(logits)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "YHAZYIwgpZc3"
   },
   "source": [
    "## Forward/Backward Dynamic Programming\n",
    "\n",
    "整個 CTC 關鍵就在計算 Dynamic Programming Table, 我們需要計算如下圖的 DP Table:\n",
    "<img src=\"https://i.imgur.com/lOPaABD.png\" height=\"70%\" width=\"70%\">\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "collapsed": true,
    "id": "DVIvRsnXiv6l"
   },
   "outputs": [],
   "source": [
    "def forward(y, label):\n",
    "  L = len(label)\n",
    "  V, T = y.shape\n",
    "  alpha = np.zeros([L,T])\n",
    "  # init first column\n",
    "  alpha[0,0] = y[label[0],0]  # TODO\n",
    "  alpha[1,0] = y[label[1],0]  # TODO\n",
    "  # run dp\n",
    "  for t in range(1,T):\n",
    "    for s in range(L):\n",
    "      k = label[s]\n",
    "      y_k_t = y[k,t]\n",
    "      alpha_tmp = alpha[s,t-1]\n",
    "      if s>0:\n",
    "        alpha_tmp += alpha[s-1,t-1]  # TODO\n",
    "      if s>1 and k!=0 and k!=label[s-2]:\n",
    "        alpha_tmp += alpha[s-2,t-1]  # TODO\n",
    "      alpha[s,t] = alpha_tmp*y_k_t\n",
    "  return alpha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "collapsed": true,
    "id": "y78iljgTAWlg"
   },
   "outputs": [],
   "source": [
    "def backward(y,label):\n",
    "  L = len(label)\n",
    "  V, T = y.shape\n",
    "  beta = np.zeros([L,T])\n",
    "  # init last column\n",
    "  beta[-1,-1] = y[label[-1],-1]  # TODO\n",
    "  beta[-2,-1] = y[label[-2],-1]  # TODO\n",
    "  # run dp\n",
    "  for t in range(T-2,-1,-1):\n",
    "    for s in range(L):\n",
    "      k = label[s]\n",
    "      y_k_t = y[k,t]\n",
    "      beta_tmp = beta[s,t+1]\n",
    "      if s<L-1:\n",
    "        beta_tmp += beta[s+1,t+1]  # TODO\n",
    "      if s<L-2 and k!=0 and k!=label[s+2]:\n",
    "        beta_tmp += beta[s+2,t+1]  # TODO\n",
    "      beta[s,t] = beta_tmp*y_k_t\n",
    "  return beta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 50
    },
    "colab_type": "code",
    "id": "kWg7qgg5EPbt",
    "outputId": "ea81766b-c6c3-4b43-c67e-f3a35273e20d"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "likelihood_by_forward = 1.2295022894015405e-05\n",
      "likelihood_by_backword = 1.2295022894015404e-05\n"
     ]
    }
   ],
   "source": [
    "# Forward and Backward likelihood should be very close\n",
    "\n",
    "alpha = forward(y,l)\n",
    "likelihood_by_forward =  alpha[-1,-1] + alpha[-2,-1]\n",
    "print('likelihood_by_forward = {}'.format(likelihood_by_forward))\n",
    "\n",
    "beta = backward(y,l)\n",
    "likelihood_by_backword =  beta[0,0] + beta[1,0]\n",
    "print('likelihood_by_backword = {}'.format(likelihood_by_backword))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "3nwPd5WCFf5A"
   },
   "source": [
    "## Gradient w.r.t. Posterior $y$\n",
    "\n",
    "[論文](https://www.cs.toronto.edu/~graves/icml_2006.pdf)的 (14) 如下:\n",
    "\n",
    "$$\n",
    "p(l|x)=\\sum_{s=1}^{|l'|}\\frac{\\alpha_t(s)\\beta_t(s)}{y_{l_s}^t}\n",
    "$$\n",
    "\n",
    "注意微分時 $\\alpha_t(s)$ 與 $\\beta_t(s)$ 是跟 $y_{l_s}^t$ 有關, 不能視為 constant, 利用它們的定義直接展開求微分即可. 譬如我們以 $t=2$ 當例子, 並且假設只有一條 alignment $\\pi=\\{\\pi_0,\\pi_1,...,\\pi_T\\}$, 其中 $y_{\\pi_2}^2=y_{l'_s}^2$, 就是先忽略 $\\sum_{\\pi\\in B^{-1}(l)}$, 不影響推導\n",
    "\n",
    "$$\n",
    "p(l|x)=\\sum_{s=1}^{|l'|}\\frac{\\alpha_2(s)\\beta_2(s)}{y_{l'_s}^2}\\\\\n",
    "=\\sum_{s=1}^{|l'|}\\frac{y_{\\pi_0}^0 y_{\\pi_1}^1 (y_{l'_s}^2)^2 y_{\\pi_3}^3 ... y_{\\pi_T}^T}{y_{l'_s}^2}\\\\\n",
    "=\\sum_{s=1}^{|l'|} {y_{\\pi_0}^0 y_{\\pi_1}^1 (y_{l'_s}^2) y_{\\pi_3}^3 ... y_{\\pi_T}^T}\\\\\n",
    "\\frac{\\partial p(l|x)}{\\partial y_{l'_s}^2}={y_{\\pi_0}^0 y_{\\pi_1}^1 y_{\\pi_3}^3 ... y_{\\pi_T}^T} \\mbox{,  只保留 index }s\\\\\n",
    "=\\frac{\\alpha_2(s)\\beta_2(s)}{(y_{l'_s}^2)^2}\n",
    "$$\n",
    "\n",
    "因此得到 Gradient:\n",
    "$$\n",
    "\\frac{\\partial p(l|x)}{\\partial y_k^t}=\\frac{1}{{y_k^t}^2}\\sum_{s\\in lab(l,k)}\\alpha_t(s)\\beta_t(s)\n",
    "$$\n",
    "\n",
    "解釋一下 $lab(l,k)$, 以我們上面的例子來說, `l = [0, 3, 0, 3, 0, 4, 0]`, `Vocabulary = [0,1,2,3,4]`, 而 $k$ 是指 Vocabulary 的哪一個值\n",
    "\n",
    "- $lab(l,0)=\\{0,2,4,6\\}$\n",
    "- $lab(l,1)=\\{\\}$\n",
    "- $lab(l,2)=\\{\\}$\n",
    "- $lab(l,3)=\\{1,3\\}$\n",
    "- $lab(l,4)=\\{5\\}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 67
    },
    "colab_type": "code",
    "id": "5syn0MgnblM7",
    "outputId": "53191cf3-378f-4a37-8176-b6ada28d2821"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0, 2, 4, 6]\n",
      "[]\n",
      "[1, 3]\n"
     ]
    }
   ],
   "source": [
    "def lab_l_k(l,k):\n",
    "  return [i for i,s in enumerate(l) if s==k]  # TODO\n",
    "print(lab_l_k(l,0))\n",
    "print(lab_l_k(l,1))\n",
    "print(lab_l_k(l,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "fRrwGv8AdnG_"
   },
   "source": [
    "因此我們的 gradient 計算公式如下\n",
    "$$\n",
    "\\frac{\\partial \\ln p(l|x)}{\\partial y_k^t}=\\frac{1}{p(l|x)}\\frac{\\partial p(l|x)}{\\partial y_k^t}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "collapsed": true,
    "id": "h8rxUftKFdZS"
   },
   "outputs": [],
   "source": [
    "def gradient(y, label):\n",
    "  L = len(label)\n",
    "  V, T = y.shape\n",
    "  alpha = forward(y,label)\n",
    "  beta = backward(y,label)\n",
    "  grad = np.zeros([V,T])\n",
    "  p = alpha[-1,-1] + alpha[-2,-1]\n",
    "  for t in range(T):\n",
    "    for k in range(V):\n",
    "      for s in lab_l_k(label,k):\n",
    "        grad[k,t] += alpha[s,t]*beta[s,t]  # TODO\n",
    "  grad /= (y*y*p)\n",
    "  return grad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 202
    },
    "colab_type": "code",
    "id": "DinVl4zLeGMh",
    "outputId": "cf7c2576-c113-44ba-e351-04e369ab2059"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[3.12717714 2.12614681 2.35360161 3.06430112 2.39102569 2.15188228\n",
      "  1.72376588 1.93891599 1.53795381 1.96067816 2.15547666 2.49423136]\n",
      " [0.         0.         0.         0.         0.         0.\n",
      "  0.         0.         0.         0.         0.         0.        ]\n",
      " [0.         0.         0.         0.         0.         0.\n",
      "  0.         0.         0.         0.         0.         0.        ]\n",
      " [2.73271271 2.55978634 2.57505254 3.09524208 2.46329795 2.30852701\n",
      "  2.43883506 1.52088175 1.59065849 0.57446128 0.21068645 0.        ]\n",
      " [0.         0.         0.         0.01646988 0.07262453 0.24797281\n",
      "  0.69042545 1.03487485 2.13589668 2.50074263 2.82067403 1.82337063]]\n",
      "(5, 12)\n"
     ]
    }
   ],
   "source": [
    "grad = gradient(y, l)\n",
    "print(grad)\n",
    "print(grad.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Aie-FDvceV4E"
   },
   "source": [
    "用數值方法確認 gradient 計算, $\\frac{df(x)}{dx}\\simeq \\frac{f(x+\\delta)-f(x-\\delta)}{2\\delta}$, 如果我們只用 $\\frac{df(x)}{dx}\\simeq \\frac{f(x+\\delta)-f(x)}{\\delta}$ 可以嗎? 是可以, 但是精確度較差, 原因可以由泰勒展開式知道:\n",
    "\n",
    "$$\n",
    "f(x+\\delta)=f(x)+f'(x)\\delta+\\frac{1}{2}f''(x)\\delta^2+o(|\\delta|^3)\\\\\n",
    "f(x-\\delta)=f(x)-f'(x)\\delta+\\frac{1}{2}f''(x)\\delta^2+o(|\\delta|^3)\\\\\n",
    "$$\n",
    "\n",
    "比較兩種方式的計算結果:\n",
    "\n",
    "$$\n",
    "\\frac{f(x+\\delta)-f(x-\\delta)}{2\\delta}=f'(x)+o(|\\delta|^2)\\\\\n",
    "\\frac{f(x+\\delta)-f(x)}{\\delta}=f'(x)+\\frac{1}{2}f''(x)\\delta+o(|\\delta|^2)=f'(x)+o(|\\delta|)\n",
    "$$\n",
    "\n",
    "可以發現 $\\frac{df(x)}{dx}\\simeq \\frac{f(x+\\delta)-f(x)}{\\delta}$ 估計的 error 高出一個 order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 437
    },
    "colab_type": "code",
    "id": "rgdDwpBeeVQE",
    "outputId": "c8b09942-93c3-4274-9e97-6a2cd258ad12"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1e-05\n",
      "1e-06\n",
      "[0, 0]：1.38e-06\n",
      "[0, 1]：1.33e-06\n",
      "[0, 3]：5.64e-06\n",
      "[0, 4]：4.98e-06\n",
      "[0, 5]：3.61e-06\n",
      "[0, 6]：1.91e-06\n",
      "[0, 7]：4.15e-06\n",
      "[0, 8]：1.82e-06\n",
      "[0, 9]：2.35e-06\n",
      "[0, 11]：4.31e-06\n",
      "[3, 1]：2.85e-06\n",
      "[3, 2]：1.26e-06\n",
      "[3, 3]：2.46e-06\n",
      "[3, 4]：2.72e-06\n",
      "[3, 6]：5.18e-06\n",
      "[3, 8]：1.99e-06\n",
      "[3, 9]：3.64e-06\n",
      "[3, 10]：1.65e-06\n",
      "[4, 3]：5.83e-06\n",
      "[4, 4]：1.82e-06\n",
      "[4, 5]：2.28e-06\n",
      "[4, 7]：4.00e-06\n",
      "[4, 8]：3.67e-06\n",
      "[4, 9]：3.65e-06\n",
      "[4, 10]：2.99e-06\n",
      "[4, 11]：2.51e-06\n"
     ]
    }
   ],
   "source": [
    "def check_grad(y, label, k=-1, t=-1, toleration=1e-3):\n",
    "  grad_1 = gradient(y, label)[k, t]\n",
    "\n",
    "  delta = 1e-10\n",
    "  original = y[k, t]\n",
    "\n",
    "  y[k, t] = original + delta\n",
    "  alpha = forward(y, label)\n",
    "  log_p1 = np.log(alpha[-1,-1] + alpha[-2,-1])\n",
    "\n",
    "  y[k, t] = original - delta\n",
    "  alpha = forward(y, label)\n",
    "  log_p2 = np.log(alpha[-1,-1] + alpha[-2,-1])\n",
    "\n",
    "  y[k, t] = original\n",
    "\n",
    "  grad_2 = (log_p1 - log_p2) / (2 * delta)  # TODO\n",
    "  if np.abs(grad_1 - grad_2) > toleration:\n",
    "    print('[%d, %d]：%.2e' % (k, t, np.abs(grad_1 - grad_2)))\n",
    "\n",
    "for toleration in [1e-5, 1e-6]:\n",
    "  print('%.e' % toleration)\n",
    "  for w in range(y.shape[0]):\n",
    "    for v in range(y.shape[1]):\n",
    "      check_grad(y, l, w, v, toleration)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "uNSuvN4ngghn"
   },
   "source": [
    "## Gradient w.r.t. Logits $u$\n",
    "\n",
    "<img src=\"https://i.imgur.com/1k0g7N2.png\" height=\"40%\" width=\"40%\">\n",
    "\n",
    "如上圖, 我們已經計算了 $\\frac{\\partial L}{\\partial y_i}$ (這邊假定固定的 $t$, 因此省略不寫), 我們希望計算 $\\frac{\\partial L}{\\partial u_k}$, 由 chain rule 知道如下關係\n",
    "\n",
    "$$\n",
    "\\frac{\\partial L}{\\partial u_k}=\\sum_i^V\\frac{\\partial L}{\\partial y_i}\\frac{\\partial y_i}{\\partial u_k}\n",
    "$$\n",
    "\n",
    "而 softmax 的 gradient 如下:\n",
    "\n",
    "$$\n",
    "\\frac{\\partial y_i}{\\partial u_k}=y_i(\\delta_{ik}-y_k)\n",
    "$$\n",
    "\n",
    "所以結合起來變成:\n",
    "\n",
    "$$\n",
    "\\frac{\\partial L}{\\partial u_k}=\\sum_i^V\\frac{\\partial L}{\\partial y_i}\\frac{\\partial y_i}{\\partial u_k}\\\\\n",
    "=\\sum_i^V \\left(\\frac{1}{p(l|x){y_i}^2}\\sum\\alpha\\beta\\right)y_i(\\delta_{ik}-y_k)\\\\\n",
    "=\\sum_i^V \\left(\\frac{\\sum\\alpha\\beta}{p(l|x)y_i}\\right)(\\delta_{ik}-y_k)\\\\\n",
    "=\\sum_i^V \\left(\\frac{\\sum\\alpha\\beta}{p(l|x)y_i}\\right)\\delta_{ik}-y_k\\sum_i^V \\left(\\frac{\\sum\\alpha\\beta}{p(l|x)y_i}\\right)\\\\\n",
    "=\\frac{\\sum\\alpha\\beta}{p(l|x)y_k}-y_k\n",
    "$$ \n",
    "\n",
    "最後一行的第一個 term 是因為只有 $i=k$時 $\\delta_{ik}=1$, 第二個 term 是因為 $\\sum_i^V\\frac{\\sum\\alpha\\beta}{y_i}=p(l|x)$, 請參考論文 $\\frac{\\alpha\\beta}{y_i}$的物理意義\n",
    "\n",
    "總結來說, 將原來 gradient 計算:\n",
    "$$\n",
    "\\frac{\\partial L}{\\partial y_i}=\\frac{\\sum\\alpha\\beta}{p(l|x){y_i}^2}\n",
    "$$\n",
    "改成如下計算:\n",
    "$$\n",
    "\\frac{\\partial L}{\\partial u_k}=\\frac{\\sum\\alpha\\beta}{p(l|x)y_k}-y_k\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "fZNt3xUJgXXo",
    "outputId": "86b62f58-2bdf-4ba8-ec1a-f9c7ad3e077f"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.2903241643264e-15\n"
     ]
    }
   ],
   "source": [
    "def gradient_logits_naive(y, label):\n",
    "  '''\n",
    "  gradient by back propagation\n",
    "  '''\n",
    "  y_grad = gradient(y, label)\n",
    "\n",
    "  sum_y_grad = np.sum(y_grad * y, axis=0, keepdims=True)\n",
    "  u_grad = y * (y_grad - sum_y_grad) \n",
    "\n",
    "  return u_grad\n",
    "\n",
    "def gradient_logits(y, label):\n",
    "  L = len(label)\n",
    "  V, T = y.shape\n",
    "  alpha = forward(y,label)\n",
    "  beta = backward(y,label)\n",
    "  grad = np.zeros([V,T])\n",
    "  p = alpha[-1,-1] + alpha[-2,-1]\n",
    "  for t in range(T):\n",
    "    for k in range(V):\n",
    "      for s in lab_l_k(label,k):\n",
    "        grad[k,t] += alpha[s,t]*beta[s,t]  # TODO\n",
    "  grad /= (y*p)\n",
    "  grad -= y  # TODO\n",
    "  return grad\n",
    "\n",
    "grad_l = gradient_logits_naive(logits, l)\n",
    "grad_2 = gradient_logits(logits, l)\n",
    "\n",
    "print(np.sum(np.abs(grad_l - grad_2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "k7mqs3_3iiZn"
   },
   "source": [
    "## Tensorflow Gradient w.r.t. Logits $u$\n",
    "\n",
    "使用 tensorflow [`tf.nn.ctc_loss`](https://www.tensorflow.org/api_docs/python/tf/nn/ctc_loss) 來驗證我們的計算\n",
    "\n",
    "tf.nn.ctc_loss(\n",
    "    **labels,\n",
    "    inputs,\n",
    "    sequence_length**,\n",
    "    preprocess_collapse_repeated=False,\n",
    "    ctc_merge_repeated=True,\n",
    "    ignore_longer_outputs_than_inputs=False,\n",
    "    time_major=True\n",
    ")\n",
    "\n",
    "很多參數, 我們就看主要的前三個\n",
    "- `labels`: An int32 SparseTensor. labels.indices[i, :] == [b, t] means labels.values[i] stores the id for (batch b, time t). labels.values[i] must take on values in [0, **num_labels**). See core/ops/ctc_ops.cc for more details.\n",
    "- `inputs`: 3-D float Tensor. If time_major == False, this will be a Tensor shaped: [batch_size, max_time, **num_classes**]. If time_major == True (default), this will be a Tensor shaped: [max_time, batch_size, **num_classes**]. The logits.\n",
    "- `sequence_length`: 1-D int32 vector, size [batch_size]. The sequence lengths.\n",
    "\n",
    "\n",
    "其中有四點要注意:\n",
    "1. `labels` 用的是 num_labels, 而 `inputs` 用的是 num_classes, 它們的關係為: num_classes = num_labels + 1. 所以 num_labels 表示的是 true labels, 而 num_classes 是多了一個 'blank' 的數量\n",
    "2. 'blank' 的 index 不像我們上面一樣用 0 表示, 會是 num_classes-1.\n",
    "3. tf 的 ctc_loss 用的 input 直接就是 logits $u$, 不像我們上面算 gradient 時給的是後驗概率 $y$\n",
    "4. tf 是 loss, 不是 likelihood, 所以算的 gradient 跟我們上面的正負號相反\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "zmjQfYfF4KN3"
   },
   "source": [
    "#### 1. 建立 inputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "collapsed": true,
    "id": "KSFarG8a4J6Z"
   },
   "outputs": [],
   "source": [
    "# logits.shape = y.shape = [V,T]\n",
    "perm_idx = list(range(1,V))+[0]\n",
    "inputs = logits[perm_idx,:].T.reshape([T,1,V]) # [max_time, batch_size, num_classes]\n",
    "inputs_tensor = tf.constant(inputs, dtype=tf.float32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "bwEF7QTY4fTo"
   },
   "source": [
    "#### 2. 建立 labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "collapsed": true,
    "id": "-ku-Nf3S4lf-"
   },
   "outputs": [],
   "source": [
    "# 原來 l = [0, 3, 0, 3, 0, 4, 0], 在這邊要變成 [4, 2, 4, 2, 4, 3, 4], true labels = {0,1,2,3}, blank={4}\n",
    "# blank 不用自己填, 因此給 tensorflow 只需要 [2,2,3] 這樣的 array\n",
    "labels_tensor = tf.SparseTensor(indices=[[0, 0], [0, 1], [0, 2]], values=[2, 2, 3], dense_shape=[1, 3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "i_iccAlu5dt2"
   },
   "source": [
    "### 3. 建立 ctc_loss, 和計算 gradient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 67
    },
    "colab_type": "code",
    "id": "HKoPZasV5gmt",
    "outputId": "bcb78732-3ffd-4e2d-fd9d-33951fc271d9"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "inputs_tensor=Tensor(\"Const:0\", shape=(12, 1, 5), dtype=float32)\n",
      "ctc_loss=Tensor(\"CTCLoss:0\", shape=(1,), dtype=float32)\n",
      "(1, 12, 1, 5)\n"
     ]
    }
   ],
   "source": [
    "ctc_loss = tf.nn.ctc_loss(labels_tensor, inputs_tensor, [12])\n",
    "print('inputs_tensor={}\\nctc_loss={}'.format(inputs_tensor,ctc_loss))\n",
    "tf_grad = tf.gradients(ys=ctc_loss,xs=inputs_tensor)\n",
    "with tf.Session() as sess:\n",
    "  tf_grad_results = sess.run(tf_grad)\n",
    "print(np.array(tf_grad_results).shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Hcw1Vd095tmx"
   },
   "source": [
    "### 4. 比較 tf 算的 gradient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "_OHiPJLP5z4_",
    "outputId": "2edc72c5-f780-482b-e0c2-753e52fcbcd7"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.170273107219122e-07\n"
     ]
    }
   ],
   "source": [
    "# tf 是 loss, 我們之前算的是 likelihood, 因此差一個正負號\n",
    "tf_grad_results = -np.array(tf_grad_results).reshape([T,V]).T\n",
    "\n",
    "# tf 的 'blank' 是最大的 index, 將他挪回 index 0, true index 全部加 1\n",
    "perm_idx = [V-1]+list(range(0,V-1))\n",
    "tf_grad_results = tf_grad_results[perm_idx,:]\n",
    "\n",
    "# print(tf_grad_results)\n",
    "# print(np.array(tf_grad_results).shape)\n",
    "\n",
    "# 我們之前計算的 gradient\n",
    "grad = gradient_logits(y, l)\n",
    "# print(grad)\n",
    "# print(grad.shape)\n",
    "\n",
    "# compare\n",
    "diff_with_tf = np.mean(np.abs(tf_grad_results - grad))\n",
    "print(diff_with_tf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "guASq8jE08Xd"
   },
   "source": [
    "## Log-scale Forward/Backward/Gradient\n",
    "\n",
    " 採用 Log scale 的前向後向以及 Gradient 計算會穩定很多"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "collapsed": true,
    "id": "djTvHNhTrcMV"
   },
   "outputs": [],
   "source": [
    "ninf = -np.float('inf')\n",
    "\n",
    "def _logsumexp(a, b):\n",
    "    '''\n",
    "    np.log(np.exp(a) + np.exp(b))\n",
    "\n",
    "    '''\n",
    "    \n",
    "    if a < b:\n",
    "        a, b = b, a\n",
    "        \n",
    "    if b == ninf:\n",
    "        return a\n",
    "    else:\n",
    "        return a + np.log(1 + np.exp(b - a)) \n",
    "    \n",
    "def logsumexp(*args):\n",
    "    '''\n",
    "    from scipy.special import logsumexp\n",
    "    logsumexp(args)\n",
    "    '''\n",
    "    res = args[0]\n",
    "    for e in args[1:]:\n",
    "        res = _logsumexp(res, e)\n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "hT5PT9DG1HgB",
    "outputId": "1a19ed77-8299-4a9c-d1a6-4e22a04f8356"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.346686827687929e-17\n"
     ]
    }
   ],
   "source": [
    "def forward_log(log_y, label):\n",
    "  L = len(label)\n",
    "  V, T = log_y.shape\n",
    "  alpha = np.ones([L,T])*ninf\n",
    "  # init first column\n",
    "  alpha[0,0] = log_y[label[0],0]\n",
    "  alpha[1,0] = log_y[label[1],0]\n",
    "  # run dp\n",
    "  for t in range(1,T):\n",
    "    for s in range(L):\n",
    "      ls = label[s]\n",
    "      log_y_ls_t = log_y[ls,t]\n",
    "      alpha_tmp = alpha[s,t-1]\n",
    "      if s>0:\n",
    "        alpha_tmp = logsumexp(alpha_tmp,alpha[s-1,t-1])  # TODO\n",
    "      if s>1 and ls!=0 and ls!=label[s-2]:\n",
    "        alpha_tmp = logsumexp(alpha_tmp,alpha[s-2,t-1])  # TODO\n",
    "      alpha[s,t] = alpha_tmp + log_y_ls_t\n",
    "  return alpha\n",
    "\n",
    "log_alpha = forward_log(np.log(y), l)\n",
    "alpha = forward(y, l)\n",
    "print(np.sum(np.abs(np.exp(log_alpha) - alpha)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "Ivd945tX2GWE",
    "outputId": "b885b8fc-67be-4bce-b94f-e0bec24c22ec"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.2692450278827943e-16\n"
     ]
    }
   ],
   "source": [
    "def backward_log(log_y, label):\n",
    "  L = len(label)\n",
    "  V, T = log_y.shape\n",
    "  beta = np.ones([L,T])*ninf\n",
    "  # init last column\n",
    "  beta[-1,-1] = log_y[label[-1],-1]\n",
    "  beta[-2,-1] = log_y[label[-2],-1]\n",
    "  # run dp\n",
    "  for t in range(T-2,-1,-1):\n",
    "    for s in range(L):\n",
    "      ls = label[s]\n",
    "      log_y_ls_t = log_y[ls,t]\n",
    "      beta_tmp = beta[s,t+1]\n",
    "      if s<L-1:\n",
    "        beta_tmp = logsumexp(beta_tmp,beta[s+1,t+1])\n",
    "      if s<L-2 and ls!=0 and ls!=label[s+2]:\n",
    "        beta_tmp = logsumexp(beta_tmp,beta[s+2,t+1])\n",
    "      beta[s,t] = beta_tmp + log_y_ls_t\n",
    "  return beta\n",
    "\n",
    "log_beta = backward_log(np.log(y), l)\n",
    "beta = backward(y, l)\n",
    "print(np.sum(np.abs(np.exp(log_beta) - beta)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "Ou_iq-Vr2oAM",
    "outputId": "e1f319ad-d8e1-4d06-a83b-866ab3f5e35a"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.5844270340181765e-13\n"
     ]
    }
   ],
   "source": [
    "def gradient_log(log_y, label):\n",
    "  L = len(label)\n",
    "  V, T = log_y.shape\n",
    "  log_alpha = forward_log(log_y,label)\n",
    "  log_beta = backward_log(log_y,label)\n",
    "  log_grad = np.ones([V,T])*ninf\n",
    "  log_p = logsumexp(log_alpha[-1,-1], log_alpha[-2,-1])\n",
    "  for t in range(T):\n",
    "    for k in range(V):\n",
    "      for s in lab_l_k(label,k):\n",
    "        log_grad[k,t] = logsumexp(log_grad[k,t],log_alpha[s,t]+log_beta[s,t])  # TODO\n",
    "  log_grad -= (2*log_y+log_p)\n",
    "  return log_grad\n",
    "\n",
    "log_grad = gradient_log(np.log(y), l)\n",
    "grad = gradient(y, l)\n",
    "print(np.sum(np.abs(np.exp(log_grad) - grad)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "6oOggns39hxO"
   },
   "source": [
    "Scaling 方法跳過"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "thdZuTmKu9jZ"
   },
   "source": [
    "## Tensorflow ctc_decoder\n",
    "\n",
    "使用上沒啥好說的, 有一點要注意的是如同上面說的 tf 會將最大的 label 視為 'blank', 所以下面 decode 會看到 label 0 但是沒有 label 4, 另外, 預設 merge_repeated=True 會去掉 blank 和 重複"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "pXDDZaYMypU_"
   },
   "source": [
    "tf.nn.ctc_greedy_decoder(\n",
    "    inputs,\n",
    "    sequence_length,\n",
    "    merge_repeated=True\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 101
    },
    "colab_type": "code",
    "id": "6mnBUFH5vQXk",
    "outputId": "26135a0d-48a6-49d1-86dd-67d02e6b72af"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[SparseTensorValue(indices=array([[0, 0],\n",
      "       [0, 1],\n",
      "       [0, 2],\n",
      "       [0, 3],\n",
      "       [0, 4],\n",
      "       [0, 5],\n",
      "       [0, 6],\n",
      "       [0, 7]], dtype=int64), values=array([0, 2, 1, 2, 3, 0, 2, 0], dtype=int64), dense_shape=array([1, 8], dtype=int64))]\n",
      "[[-10.676455]]\n"
     ]
    }
   ],
   "source": [
    "decoded, neg_sum_logits = tf.nn.ctc_greedy_decoder(inputs_tensor, [12])\n",
    "with tf.Session() as sess:\n",
    "  decoded_out, neg_sum_logits_out = sess.run([decoded, neg_sum_logits])\n",
    "print(decoded_out)\n",
    "print(neg_sum_logits_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "oWB5DXp_yRm4"
   },
   "source": [
    "tf.nn.ctc_beam_search_decoder(\n",
    "    inputs,\n",
    "    sequence_length,\n",
    "    beam_width=100,\n",
    "    top_paths=1,\n",
    "    merge_repeated=True\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 101
    },
    "colab_type": "code",
    "id": "8yWBVuD3yWZc",
    "outputId": "0ebfc825-aa9d-4371-80ab-2507b7f31ef9"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[SparseTensorValue(indices=array([[0, 0],\n",
      "       [0, 1],\n",
      "       [0, 2],\n",
      "       [0, 3],\n",
      "       [0, 4],\n",
      "       [0, 5],\n",
      "       [0, 6],\n",
      "       [0, 7]], dtype=int64), values=array([0, 2, 1, 2, 3, 0, 2, 0], dtype=int64), dense_shape=array([1, 8], dtype=int64))]\n",
      "[[-10.676455]]\n"
     ]
    }
   ],
   "source": [
    "decoded_beam, neg_sum_logits_beam = tf.nn.ctc_greedy_decoder(inputs_tensor, [12])\n",
    "with tf.Session() as sess:\n",
    "  decoded_beam_out, neg_sum_logits_beam_out = sess.run([decoded_beam, neg_sum_logits_beam])\n",
    "print(decoded_beam_out)\n",
    "print(neg_sum_logits_beam_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "RB0ryUsVzirA"
   },
   "source": [
    "更多 decoding 的算法細節, 請參考 [DingKe](https://github.com/DingKe/ml-tutorial/blob/master/ctc/CTC.ipynb) 的文章"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "collapsed": true,
    "id": "7aqnGyac4cxq"
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "CTC_Practice_Answer.ipynb",
   "provenance": [],
   "toc_visible": true,
   "version": "0.3.2"
  },
  "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.5.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}