{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# theano的scan函数\n",
    "`scan`函数是theano中一个十分重要的概念，利用它我们可以处理时序数据，从而完成许多复杂的计算过程。其原理有点类似一个高度封装过的for循环，每个时刻都调用相同的回调函数处理该时刻的数据，最后再将处理的结果按照时间顺序堆叠、汇总。  \n",
    "`scan`相对于for循环的好处：\n",
    "* 迭代次数可以作为符号图的一部分\n",
    "* 最小化GPU数据传输\n",
    "* 可以按时序计算梯度\n",
    "* 比python在for循环中调用编译好的函数速度更快\n",
    "* 可以通过检测实际需要使用的内存量来减少总的内存使用\n",
    "\n",
    "它的缺点是过于复杂，难于调试，因此也成为笔者学习theano遇到的第一道坎。  \n",
    "\n",
    "scan函数的定义：\n",
    "```py\n",
    "theano.scan(fn, sequences=None, outputs_info=None, non_sequences=None, n_steps=None, truncate_gradient=-1, go_backwards=False, mode=None, name=None, profile=False, allow_gc=None, strict=False)\n",
    "```\n",
    "* fn：供scan在每个时刻调用的函数句柄，fn的输入参数顺序应该与传给scan函数的参数顺序一致  \n",
    "fn读取参数的顺序为：sequence->outputs_info->non_sequences\n",
    "* sequences：存放作为输入的时序数据，类型可以是列表或字典  \n",
    "举个例子，假设sequences=[a, b]，那么执行scan时fn会依次读取该列表中每个变量第t时刻的数据a[t]，b[t]。需要注意的是传入时要通过dimshuffle把时间维放在axis0\n",
    "* outputs_info：scan函数输出的初始值，同时通过outputs_info我们可以实现递归计算，其中taps是scan实现递归运算的关键。默认情况下，如果不指定taps参数，则outputs_info输入参数的前一时刻的值将会参与当前时刻的计算，如果outputs_info为None，表示不使用任何tap，此时不具有递归结构。如果使用多个时刻的taps，需要有额外的维度，比如taps=[-5,-2,-1]，而输出的shape为(6,3)，那么初始状态至少需要(5,)+(6,3)=(5,6,3)才能足够容纳scan函数的输出。如果类型为dict，需要指定initial（初值）和taps（参与计算的时刻），比如taps=-2表示计算过程需要用到t-2时刻的数据\n",
    "* non_sequences：作为输入的非时序数据\n",
    "* n_steps：循环执行fn的次数，可以是符号变量\n",
    "* truncate_gradient：用于限定BPTT中梯度回传的时刻数，如果为-1则表示不适用梯度截断\n",
    "* go_backwards：默认为False，如果为True，则从序列的末尾后往前计算各时刻输出，该参数一般用在双向的递归结构，比如BiGRU和BiLSTM的实现中\n",
    "* 返回值：scan函数的返回值是fn在每个时刻上处理的结果按照时间维度的堆叠\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们以累加器为例，来研究scan函数到底是怎么工作的：\n",
    "$$ sum(n) = \\sum_{i=0}^n i$$\n",
    "这个关系可以表示为如下的递归关系：\n",
    "$$sum(n)=sum(n-1)+n$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  0.   1.   3.   6.  10.  15.]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "DEBUG: nvcc STDOUT mod.cu\r\n",
      "   ���ڴ����� C:/Users/hschen/AppData/Local/Theano/compiledir_Windows-10-10.0.14393-Intel64_Family_6_Model_60_Stepping_3_GenuineIntel-2.7.12-64/tmpiviokf/265abc51f7c376c224983485238ff1a5.lib �Ͷ��� C:/Users/hschen/AppData/Local/Theano/compiledir_Windows-10-10.0.14393-Intel64_Family_6_Model_60_Stepping_3_GenuineIntel-2.7.12-64/tmpiviokf/265abc51f7c376c224983485238ff1a5.exp\r\n",
      "\n",
      "Using gpu device 0: GeForce GTX 960M (CNMeM is disabled, cuDNN 5103)\n",
      "C:\\Anaconda2\\lib\\site-packages\\theano\\sandbox\\cuda\\__init__.py:600: UserWarning: Your cuDNN version is more recent than the one Theano officially supports. If you see any problems, try updating Theano or downgrading cuDNN to version 5.\n",
      "  warnings.warn(warn)\n"
     ]
    }
   ],
   "source": [
    "import theano\n",
    "import theano.tensor as T\n",
    "import numpy as np\n",
    "\n",
    "n = T.iscalar()\n",
    "acc_out, updates = theano.scan(lambda i, acc_sum: acc_sum + i, sequences=T.arange(n+1), \n",
    "                               outputs_info=T.constant(np.float64(0)))\n",
    "accumulate_sum = theano.function([n], acc_out)\n",
    "print accumulate_sum(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "为了获得累加结果，我以0作为outputs_info的initial_state的初始值。接着传给sequences一个以0起始、增量为1的整数序列，T.arange是theano版本的np.arange。随后我定义了一个匿名函数，这个函数的第一个参数`i`是序列的第i个元素；第二个参数`acc_sum`是上一个时刻的输出值，即$sum(i-1)$；而这个函数的返回值是$sum(i)=sum(i-1)+i$。scan函数按照时间顺序计算每个时刻的输出，并将结果按照时间顺序堆叠成一个np.ndarray数组：$[sum(0), sum(1),...,sum(n)]$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面结合一些scan函数的examples进行讲解，以加深理解"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 平方和累加\n",
    "给定一个正整数$n$，我们要通过scan函数求解如下的式子：\n",
    "$$\\sum_{i=1}^n i^2$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "55\n"
     ]
    }
   ],
   "source": [
    "n = T.iscalar()\n",
    "acc_out, updates = theano.scan(lambda i, acc_sum: acc_sum + i**2, sequences=T.arange(n+1), \n",
    "                               outputs_info=T.constant(np.int64(0)))\n",
    "acc_out = acc_out[-1]\n",
    "accumulate_square_sum = theano.function([n], acc_out)\n",
    "print accumulate_square_sum(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 斐波那契数列\n",
    "斐波那契数列的递推式为\n",
    "$$ x(n)=x(n-1)+x(n-2),(n\\geq 2)$$\n",
    "其中$x(0)=0, x(1)=1$  \n",
    "其scan版本的实现如下："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1  2  3  5  8 13 21 34 55 89]\n"
     ]
    }
   ],
   "source": [
    "n = T.iscalar()\n",
    "x0 = T.ivector()\n",
    "fib_out, updates = theano.scan(lambda xtm1, xtm2: xtm1 + xtm2, \n",
    " outputs_info=[dict(initial=x0, taps=[-1,-2])], n_steps=n)\n",
    "fib_out = fib_out\n",
    "fib = theano.function([x0, n], fib_out)\n",
    "print fib(np.int32([0, 1]), 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 计算圆周率\n",
    "圆周率可以通过下面的积分式计算\n",
    "$$\\pi=2\\int_{-1}^1 \\sqrt{1-x^2}dx$$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.14156113248\n"
     ]
    }
   ],
   "source": [
    "inp=T.dvector()\n",
    "dx = T.dscalar()\n",
    "pi, updates = theano.scan(lambda xt, pi_sum: pi_sum+2.*T.sqrt(1-xt**2)*dx, sequences=[inp], \n",
    "                          outputs_info=T.constant(np.float64(0.)))\n",
    "pi = pi[-1]\n",
    "cal_pi = theano.function([inp, dx], pi)\n",
    "n_interval = 100000\n",
    "print cal_pi(np.linspace(-1, 1, n_interval)[1:-1], 2. / n_interval)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 计算泰勒展开\n",
    "我们用scan实现$e^x$的泰勒展开式：\n",
    "$$ e^x=1+x+\\frac{1}{2!}x^2+\\frac{1}{3!}x^3+\\cdots = \\sum_{n=0}^\\infty \\frac{1}{n!}x^n$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "calc_exp(1)=2.718282\n",
      "calc_exp(0)=1.000000\n",
      "calc_exp(0)=0.367879\n",
      "whether calc_exp(1) equals np.exp(1):True\n"
     ]
    }
   ],
   "source": [
    "n = T.iscalar()\n",
    "x = T.dscalar()\n",
    "#factorial = T.cumprod(T.arange(1, n + 1)))\n",
    "factorial = T.gamma(n+1)\n",
    "\n",
    "def fn(n, power, exp_sum, x):\n",
    "    power = power*x\n",
    "    return  power, exp_sum + 1./T.gamma(n+1)*power\n",
    "result, updates = theano.scan(fn, sequences=T.arange(1, n),\n",
    "                         outputs_info=[T.constant(np.float64(1.)), T.constant(np.float64(1.))],\n",
    "                         non_sequences=x)\n",
    "exp_ = result[1][-1]\n",
    "calc_exp = theano.function([n, x], exp_)\n",
    "print \"calc_exp(1)=%f\"%calc_exp(15, 1)\n",
    "print \"calc_exp(0)=%f\"%calc_exp(15, 0)\n",
    "print \"calc_exp(0)=%f\"%calc_exp(15, -1)\n",
    "print \"whether calc_exp(1) equals np.exp(1):%s\"%np.allclose(calc_exp(15, 1), np.exp(1))"
   ]
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}