gistfile1.txt

{
 "metadata": {
  "name": "Theano Scan"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "raw",
     "metadata": {},
     "source": "This notebook will introduce you to the use of scan in theano, at first it can seem a bit obscure, but you will quickly get the idea.\n\nScan is here to help us add loops into theano natively. You can of course just make a loop where at every iteration you call a theano function, but the former can be faster. \n\nYou can also use scan to iterate over a tensor, the elements of a vector, rows of a column and so on.\n\nRather than explain exactly how it works, I think it is less confusing to just dive into some simple examples."
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "import theano as th\nfrom theano import tensor as T\nimport numpy as np",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "raw",
     "metadata": {},
     "source": "Let's try counting to ten."
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "count=0\nfor i in range(1,11):\n    count+=1\nprint count",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "10\n"
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "raw",
     "metadata": {},
     "source": "Now we try it theano style."
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "i = T.iscalar('i')\nresults, updates = th.scan(fn=lambda previous_count: previous_count+1,\n                            outputs_info=0,\n                            n_steps=i)",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "raw",
     "metadata": {},
     "source": "Let's try and understand some of this.\n\nFirstly the integer variable i is there to tell scan how many times to iterate.\n\nNext we pass scan a function which it applies on each iteration, the first argument to the function is always the output from the previous call of the function.\n\nWe also have to tell scan what the outputs of the function should look like and initialize the count to zero.\n\nWhat we create are a variable results that will hold all of the function outputs, and updates, which are used as updates in a theano function."
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "f = th.function(inputs=[i], outputs=results, updates=updates)\nprint f(10)",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "[ 1  2  3  4  5  6  7  8  9 10]\n"
      },
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": "/usr/local/lib/python2.7/dist-packages/theano/scan_module/scan_perform_ext.py:85: RuntimeWarning: numpy.ndarray size changed, may indicate binary incompatibility\n  from scan_perform.scan_perform import *\n"
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "raw",
     "metadata": {},
     "source": "Notice that it returns the result at each iteration. If you only care about the end result you could create a different function."
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "f = th.function(inputs=[i], outputs=results[-1], updates=updates)\nprint f(10)",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "10\n"
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "raw",
     "metadata": {},
     "source": "This is worth doing since Theano is smart and will optimize by discarding the intermediate values once they are no longer needed.\n\nNext we will try a slightly trickier form of iteration. You have seen that scan stores the results of all the iterations, and so the function to be applied at each stage can accept as arugments any of the previous results.\n\nThis allows us to compute any recurrence relation we like, let's try everybody's favourite: the Fibonacci sequence."
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "i = T.iscalar('i') #Number of iterations.\nx0 = T.ivector('x0') #Initializes the recurrence, since we need the previous \n                    #two terms in the Fibonacci sequence, we need an inital vector\n                    #with two terms.\nresults, updates = th.scan(fn=lambda f_m_1,f_m_2: f_m_1+f_m_2,\n                            outputs_info=[{'initial':x0, 'taps':[-2,-1]}],\n                            n_steps=i)",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 64
    },
    {
     "cell_type": "raw",
     "metadata": {},
     "source": "In order to set this up, we need to provide scan with a bit more information. \n\nIn the outputs_info argument we provide a list of dictionaries. W\n\nWe need to say what the inital variable is, in this case x0, and what arguments will be passed to the function at each step. The key 'taps' is a list of indices to the results to be used as arguments, -1 and -2 mean the last and penultimate entries.\nBy default, taps is set to [-1]."
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "f=th.function(inputs=[i,x0], outputs=results, updates=updates)\nprint f(50, np.asarray([0,1], dtype=np.int32))\n",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "[          1           2           3           5           8          13\n          21          34          55          89         144         233\n         377         610         987        1597        2584        4181\n        6765       10946       17711       28657       46368       75025\n      121393      196418      317811      514229      832040     1346269\n     2178309     3524578     5702887     9227465    14930352    24157817\n    39088169    63245986   102334155   165580141   267914296   433494437\n   701408733  1134903170  1836311903 -1323752223   512559680  -811192543\n  -298632863 -1109825406]\n"
      }
     ],
     "prompt_number": 68
    },
    {
     "cell_type": "raw",
     "metadata": {},
     "source": "I deliberately set the number of iterations high, you can see that when we get too large for an int32 to handle, we get overflow. So next let's try and some flow control to our loops, so we can end early if a condition is met, like say overflow."
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "def fib(f_m_1, f_m_2):\n    ret= f_m_1+ f_m_2\n    return ret, th.scan_module.until(ret <0)\ni = T.iscalar('i') #Number of iterations.\nx0 = T.ivector('x0') #Initializes the recurrence, since we need the previous \n                    #two terms in the Fibonacci sequence, we need an inital vector\n                    #with two terms.\nresults, updates = th.scan(fn=fib,\n                            outputs_info=[{'initial':x0, 'taps':[-2,-1]}],\n                           \n                            n_steps=i)",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 85
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "f=th.function(inputs=[i,x0], outputs=results, updates=updates)\nprint f(50, np.asarray([0,1], dtype=np.int32))",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "[          1           2           3           5           8          13\n          21          34          55          89         144         233\n         377         610         987        1597        2584        4181\n        6765       10946       17711       28657       46368       75025\n      121393      196418      317811      514229      832040     1346269\n     2178309     3524578     5702887     9227465    14930352    24157817\n    39088169    63245986   102334155   165580141   267914296   433494437\n   701408733  1134903170  1836311903 -1323752223]\n"
      }
     ],
     "prompt_number": 86
    },
    {
     "cell_type": "raw",
     "metadata": {},
     "source": "This conditional ending is a bit obscure, but you see that the key is to return two values from the function, the second being a th.scan_module.until(termination_condition)."
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Now let's do something a little less trival and implement the Newton-Raphson method for finding the zeros of a function with continuous second derivative. You can read about the method here: http://en.wikipedia.org/wiki/Newton%27s_method."
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "#First we create a symbolic variable.\nx=T.dscalar('x')\nf = 6*x**3 - 2*x**2 +9*x+ 1 + T.cos(x)\nf_prime = T.grad(f,x)\nf_prime_2 = T.grad(f_prime, x)\n\n#Then the compiled theano functions for plotting.\nF = th.function(inputs=[x], outputs=f)\nF_prime = th.function(inputs=[x], outputs=f_prime)\nF_prime_2 = th.function(inputs=[x], outputs=f_prime_2)\n\n#Now let's make a plot.\nxs = np.linspace(-1,1,1000)\ny1 = [F(z) for z in xs]\ny2 = [F_prime(z) for z in xs]\ny3 = [F_prime_2(z) for z in xs]\n\nimport matplotlib.pyplot as plt\n%matplotlib inline\nplt.plot(xs, y1, label='f')\nplt.plot(xs, y2, label='f\\'')\nplt.plot(xs, y3, label='f\\'\\'')\nplt.plot(xs, [0 for z in xs], label='zero')\nplt.legend()",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 128,
       "text": "<matplotlib.legend.Legend at 0x7fef820b21d0>"
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcVXX+x/HXZXE3EU1AEFEWCRfADTNNXHBLUVscpcxG\nnV+To1aWW9a4pIJZ41hmzmSlNS3qVOKUGlhi5RLuqaigoiICLoALqGzn98c3LppgeO+Fu32ej8d5\n3IV77/l4ub7v4Xu+i07TNA0hhBA2zcHcBQghhKh6EvZCCGEHJOyFEMIOSNgLIYQdkLAXQgg7IGEv\nhBB2wCRhX1xcTGhoKIMHDwYgOzubiIgIAgIC6Nu3L7m5uabYjRBCCAOZJOyXLFlCUFAQOp0OgJiY\nGCIiIkhOTqZ3797ExMSYYjdCCCEMZHTYnz17lg0bNjBu3DhKx2etX7+e0aNHAzB69GjWrVtn7G6E\nEEIYweiwf/HFF1m0aBEODmUvlZWVhZubGwBubm5kZWUZuxshhBBGMCrsv/nmG5o0aUJoaCgVzbqg\n0+n0zTtCCCHMw8mYJ2/fvp3169ezYcMGbty4wZUrVxg1ahRubm5kZmbi7u5ORkYGTZo0ueO5fn5+\nnDhxwpjdCyGE3fH19eX48eP3/kTNRBISErRBgwZpmqZpU6ZM0WJiYjRN07To6Ght2rRpdzzehLsW\nmqbNmjXL3CXYFHk/TUveT9MxNDtN2s++tLlm+vTpxMfHExAQwA8//MD06dNNuRshhBD3yKhmnFv1\n6NGDHj16AODq6srmzZtN9dJCCCGMJCNobUR4eLi5S7Ap8n6alryf5qf7rQ2o+nes01XYg0cIIUT5\nDM1OkzXjCCGEsVxdXcnJyTF3GRahYcOGZGdnm+z15MheCGExJBfKVPReGPoeSZu9EELYAQl7IYSw\nAxL2QghhByTshRCiEo4dO0ZISAj33XcfS5cuNXc590x64wghRCW88cYb9O7dm/3795u7FIPIkb0Q\nQlTC6dOnCQoKMncZBpOul0IIi2GpudCrVy9+/PFHnJ2dcXZ2Zu/evfj5+VXpPk3d9VLCXghhMSw5\nF3r27MmoUaMYM2ZMtezP1GEvbfZCCKthqnWQDP0+sdQvosqQsBdCWA1zZ601r7onJ2iFEMIOmDXs\nky4kmXP3QghxT6y5GceosL9x4wZhYWGEhIQQFBTEjBkzAMjOziYiIoKAgAD69u1Lbm5uuc/vsbIH\nk7+bzOUbl40pQwghqoVZm3GKi2HlSoOfbnRvnPz8fOrUqUNRURHdunXjzTffZP369TRu3JipU6ey\ncOFCcnJyiImJuX3HOh1Z17KY+f1Mvkn5hvm95vNMyDM46KRlSQh7Zcm9caqb/r3QNNi4EaZNgwYN\n0G3bZt6ul/n5+fTo0YOVK1fy2GOPsXXrVtzc3MjMzCQ8PJyjR4+W/w8Bdp/bzcSNEykuKeadAe8Q\n5hVmipKEEFZGwr6MTqdD27ULpk6FjAyIiYHISHQODuaZ4rikpISQkBDc3Nzo2bMnrVu3JisrCzc3\nNwDc3NzIysq662t0bNqRbWO2MaHzBIatHsafY/9M5rVMY0sTQgjrNmQIjBgBBw+q60Y0Ixnd9dLB\nwYH9+/dz+fJl+vXrx5YtW277uU6nq7Cda/bs2frr4eHhPB3+NEMDhzLvx3m0WdaGGd1mMDFsIjUc\naxhbphBCWJ/kZBJ27SJh3jyjX8qkI2hff/11ateuzYoVK0hISMDd3Z2MjAx69ux512ac8hy7eIzJ\ncZM5dvEYC/ss5NEHHrXqPq5CiD8mzThlLGqlqosXL+p72ly/fp34+HhCQ0OJjIxk1apVAKxatYqh\nQ4fe82u3atyKb6O+Zfmg5cz9cS4Pr3yYxPREY8oVQgi7ZdSR/cGDBxk9ejQlJSWUlJQwatQopkyZ\nQnZ2NsOHD+fMmTP4+PiwZs0aXFxcbt/xPXw7FZcU8/GBj3lty2s83PxhontH09yluaFlCyEslBzZ\nl7HridDyCvJ4c/ubvJ34Nn9p/xdmdJtBg1oNqqhCIUR1k7AvY1HNONWtbo26zAqfxcHnDnI+7zyt\nlrbi3cR3KSwuNHdpQghh0awq7Es1rd+UD4d8yHdPfUfssVgeePcBPj/4OSVaiblLE0LYqFuXJXzn\nnXfMXc49s8qwLxXsHkzcqDj+Pfjf/POXf9L+X+3ZmLJR/gwUQphc6bKEV65cYe/evfpOKNbCqsO+\nVK8Wvdg5diezeszipbiX6LGyB9vObDN3WUIIG2LtyxLaRNiDOmkx7IFhHHzuIGNCxxD1VRSRn0dy\nMOuguUsTQli5Xr16kZCQwIQJE6hfvz6FhYVWN+7Hqnrj3IubRTdZvns50T9HE+Ebwawes/Bzrdo1\nI4UQxrHk3jiyLKGFqulUk+e7PM+Y0DEs3rmYBz94kEEBg3i1+6v4uvqauzwhhAF0c0xzNK3NMuwL\nxVK/iCrDZsO+VP2a9fl7j78zKWwSS3YuIWxFGJGtIpnZfaaEvhBWxtCQNhVra7q5lc202f8Rl1ou\nzAqfxfFJx/Fu4E3YijDGxo7lZM5Jc5cmhBBVzm7CvpRLLRdmh88mZWIKXvd50fn9zoxbP47UnFRz\nlyaEsHDW3Ixjd2FfqmHthszpOYfkick0rd+Uju93ZPS60Ry5cMTcpQkhLJQ1N+PYbG+ce5VzPYdl\nu5bxduLbPNTsIWZ0m0Enz07mLksIu2JpuWBOdj0RWnXIL8xnxd4VvLn9TVo1bsWMbjPo6dPTqr/R\nhbAWlpoL5iBhX00Kigv47OBnxPwcQ4NaDZjRbQaRrSJlQXQhqpCl50J1krCvZsUlxaw7uo7on6O5\nXnSdlx98mai2UdR0qmnu0oSwOdaSC9Xh9+9FUUkRaw6v4cl2T0rYVyVN0/g+9Xve2vEW+zP3M77j\neP7a8a/cX/d+c5cmhM2wtlyoSqXvxZWbV/hg7wcs+WUJzV2a8+Off6z++ezT0tLo2bMnrVu3pk2b\nNrz99tsAZGdnExERQUBAAH379tUvXWjNdDodfVr2YeOTG9k8ajNnLp8hYGkAz/7vWY5ePPrHLyCE\nEPfo5biXabGkBYnnEln7xFq2PrPV4Ncy6sg+MzOTzMxMQkJCuHbtGh06dGDdunV89NFHNG7cmKlT\np7Jw4UJycnKIiYm5fcc28A1+Pu88y3Yt473d79GpaScmPzhZTuYKYQRbyAVT0el0vPTdS0wKm4R3\nA+/b7jd7M87QoUOZMGECEyZMYOvWrbi5uZGZmUl4eDhHj95+9GtLv9Trhdf59OCn/GPHP6jhWIO/\ndfobUW2jqFujrrlLE8Kq2FIuGMtiT9CeOnWKHj16cOjQIby9vcnJyQFUW7erq6v+trEFW7ISrYT4\nE/G8u+tdtqdt5+ngpxnfabzMtilEJdliLhjKIme9vHbtGo899hhLliyhfv36dxRWUbPG7Nmz9dfD\nw8MJDw83RTlm46BzoJ9fP/r59eNU7imW715O1w+60t6jPRM6T2CA3wAcHRzNXaYQwgDHjh3jT3/6\nEydPnmT+/PlMnDixWvabkJBAQkKC0a9j9JF9YWEhgwYNYsCAAbzwwgsABAYGkpCQgLu7OxkZGfTs\n2dOmm3Hu5kbRDdYcXsPSxKVcyL/Acx2fY0zoGBrXaWzu0oSwOJacC2PHjsXFxYW33nqLP//5z4SH\nhzN69GhWrlzJ1q1b+eijj0y6P1Mf2RvVG0fTNMaOHUtQUJA+6AEiIyP16zOuWrWKoUOHGrMbq1bL\nqRZPBz9N4l8SWfP4GpIuJOH3th8j/juCzSc3yyLpQliJ3y9LaG0dMYwK+23btvGf//yHLVu2EBoa\nSmhoKJs2bWL69OnEx8cTEBDADz/8wPTp001Vr1Xr5NmJlUNXkvp8Kt29u/Ny3Mv4vu3LvB/nkX4l\n3dzlCSEqUN6yhKXu1lRtSWRQlRlpmsaejD2s2LuCNYfX8JD3Q4wLHcdA/4E4Ozqbuzwhqp0l54K1\nL0soE72YkU6no2PTjiwftJy0F9N47IHHWLR9Ed7/9GZK3BRZLF2I39PpTLMZyFK/iCpDwt5C1K1R\nl2dCnuHnMT/zw9M/4OTgxMDPBhKyPIS3tr9FxtUMc5cohPlpmmk2A1lDc01FJOwt0AP3P0B0n2hO\nv3Caxf0Wc+jCIYKWBdH/P/359NdPySvIM3eJQggrI2FvwRx0DvRs0ZOPhnxE+uR0RgeP5tODn+L5\nD09GrxvNxpSNFBQXmLtMIeyGNOOIKlfHuQ4j245kw5MbODrhKO3d2/P6j6/T9K2mjFs/jrgTcRQW\nF/7xCwkhDGbNzTjSG8fKnbl8hv8m/Zc1h9dwIucEwwKHMbz1cMJ9wnFyMMkAaSGqjeRCGYudG+ee\ndyy/VJM7lXuKtYfXsiZpDadzTzMscBhDA4fSq0UvWWxFWLb9+2H6dHTffSe58BsJe1EpJ3NO8mXS\nl8Qei+XQ+UNE+EYwpNUQHvF/hIa1G5q7PCGUU6fgtdcgPh5eew3dhAmSC7+RsBf37Hzeeb5J/obY\nY7FsSd1CJ89ODGk1hCGthtDcpbm5yxP26NIlmD8fVq2CCRPg5Zehfn3JhVvc+l6UlMCuXfDll7Bo\nkYS9qIT8wnziT8Sz7tg6vkn+Bo96HvT3609/v/481Owhae4RVSs/H5YsgbfeguHD4e9/B3d3/Y8l\nF8rodDri4jRiYyE2FurVg8ceg/nzJezFPSouKWbXuV1sOr6JTcc3ceTiEXo076EP/5YNW5q7RGEr\niopg5UqYPRu6doV58yAg4I6HSS6U0el0dOmiMWQIREZC6Rxs0owjjHYp/xLxJ+P14d+gVgP6+faj\nT8s+PNz8YVxquZi7RGFtNA3Wr4cZM6BJE3jjDejcucKHSy6UkTZ7US1KtBJ+zfqVTcc38UPqD+w4\nu4NWjVrRq0Uvevr0pJt3N+rXrP/HLyTs17ZtMG0aXL4MCxfCgAF/OC+N5EIZCXthFjeLbpKYnsiW\nU1vYcmoLu8/tpk2TNvTy6UUPnx6EeYbRoFYDc5cpLMGRI+pIfu9eeP11eOopcKzcCm2SC2Uk7IVF\nuF54nR1nd7AldQs/nfmJPRl7aNmwJV29uvKQ90N0bdaVFi4trHrEobhH586pNvmvv1ZH9BMmQK1a\nlXpqiVbC0YtHad2kteTCbyTshUUqLC5kf+Z+tqVtY3vadralbaNEK6Frs6509epKJ89OhLqHStOP\nLbp8WbXFL18O48bB9OnQ8O5jOS7fuMzuc7vZcXYH29O2s+PsDlxru3Ly+ZOSC7+xuLAfM2YM3377\nLU2aNOHgQTX/enZ2Nn/60584ffo0Pj4+rFmzBheX20/uSdjbNk3TOHP5DNvStrEjbQe7M3bza9av\neDfwpmPTjnTw6EDHph0JcQ+hXo165i5XGOLmTXjvPYiOhkcegTlzoFmzOx52+cZl9mbsZU/GHnaf\n282ejD1kXM0g2D2Yrl5d6dqsKw82exD3eu42nQtFRUU4OVV+ChOLC/uffvqJevXq8fTTT+vDfurU\nqTRu3JipU6eycOFCcnJyiImJMUnBwnoVlRSRdCFJ/Yc/t4fdGbs5dP4QPi4+hLiH0Ob+NrR1a0vb\nJm3xbuAtTUCWqqQEPv8cXn0V2rRRYd+mjf4L/tD5Qxw8f5ADWQfYfW43GVczCHEPoYNHBzo07UAH\njw4ENg7E0eHOdnxLzYXVq1czbtw4/e2CggK6du3Kd999xyuvvMLatWu5efMmw4YNY/HixdSqVYuE\nhASeeuopJk2axOLFi+nbty/vv/8+U6dOZe3atQAMHz6chQsXUqNGjTv2aXFhD3Dq1CkGDx6sD/vA\nwEC2bt2Km5sbmZmZhIeHc/ToUZMULGxLYXEhhy8c5kDmAQ6eP6i2rIPkFebRpkkb2jZpS5smbWjT\npA0BjQLwqOchXwLmomkQF4c2bRqFzo4cenkUP/s46MP98PnD3FfzPv3vq51bu7sGe3msIReuXr1K\nWFgYL774IklJSaSmprJy5UqcnJyIioqiTZs2LFiwgISEBCIiInj55ZeZO3cuxcXFLFiwgM2bN7N+\n/XoAhgwZQu/evZk7d+4d+7GKsG/YsCE5OTmA+nPe1dVVf9vYgoV9uJR/SR8iB7MOcvjCYVKyU8gv\nzMff1Z+ARgF3bDIOwLRyrueQfCmZlOwUru7YSrd3/8d953OZ2VvHxna1af1bqJd+Kbdu0hrX2q5G\n7fOPckGXkGDU65fSwsMNel5JSQmRkZE0b96cd999l3r16vHrr7/SsqUagLhjxw6efPJJTp48SUJC\nAv369ePq1av6I3c/Pz+WLl1K//79AYiLi+PZZ58lNTX1jn2ZOuyrfA7cu628Pnv2bP318PBwwg38\nBQjb06hOI3r49KCHT4/b7s+9kUvKpRSSLyWTfCmZb1O+ZfHOxSRfSsZR54iPiw/NXZrTvEFzdb1B\nc5q7qOuNajeSvwp+o2kaF/Mvcvryac5cPsPpXHV55oq6fir3FAXFBYSXeDN941XaHM0m6a+Pc33s\nON52f8DoUDe4bjNnxMyZM8nLy+Ptt9/m/Pnz5Ofn06FDB/3PNU2jpKREf/v++++/rYnm3LlzNG9e\nNh+Vt7c3586du+s+ExISSDDBl1yVhH1p8427uzsZGRk0adKk3MfdGvZCVIZLLRc6eXaik2en2+7X\nNI3s69mcyj3F6cunOZ17mtOXT/PTmZ/UfbmnyS/Mx72eOx71PXCv54573bLrHvU8cKvnhmttVxrW\naohLLZdKNz1YCk3TuHLzCjk3cjifd57Ma5lkXcsi81qmup6Xpb9Mv5JOLadaNHdpjncDb5o3UJdd\nm3XFu4E3PoV1afKPf6H77DN44QXY9CJd6tn3ifQvvviC1atXs2vXLhwdHWncuDG1a9cmKSkJDw+P\ncp/z+4OLpk2bcurUKR544AEAzpw5Q9OmTe+6398fCM+ZM8eg+qsk7CMjI1m1ahXTpk1j1apVDB06\ntCp2I4SeTqejUZ1GNKrTiA5NO5T7mPzCfLKuZZFxLYPMa5lkXFWXu9J3kXEtg6y8LLKvZ5NzPYcr\nN69Qr0Y9XGu7qi+A2g1pWKsh9WvUp26NutR1rksd5zr663VrqNu1nWrj5OBU4ebo4EhxSTHFWnGF\nlzeLbpJfmE9+YT7Xi67rr5du1wqukXsjl5wbOeRczyHnRg65N3K5fOMytZ1r41LLBbe6brjXc9df\ntmrcih4+PfT3Na3ftPxusHl58I9/qMnKnnxSDZCq4GDNnuzbt4+JEyeyefNmGjVqBICDgwN/+ctf\neOGFF1i6dCn3338/6enpHD58mL59+5b7OiNHjmTevHl06qQOVubOncuoUaOq5d9gdNiPHDmSrVu3\ncvHiRZo1a8bcuXOZPn06w4cP54MPPtB3vRTC3Oo416FFwxa0aNjiDx9bXFLM5ZuXybmeo74AbqjL\nawXXyCvII68wj7yCPLKvZ6vrhXkqnAuvU6wVU1RSVO5WXFKMg84BRwdHHHWO5V7WdKxJHec6t221\nnWpTx7kODWs1pG6Nuvq/Pkq/hFxqueBSywVnR2fD3pzCQlixQo147dEDfvkFfH0Ney0btH79enJz\nc+nWrZv+vocffpivvvqKuXPn0qVLFy5evIinpyfjx4/Xh/3vj+xfffVVrly5Qrt27QDVG+fVV1+t\nln+DDKoSwp6VlMDataobZcuWqhtl+/ZmK0dyoYzVnaAVQliozZvVaFdQo1979zZvPaJKSdgLYW/2\n7FEhf+qUWi3q8cfBwaFaSyguhhMn4OBB+PXXsktRdaQZRwh7cfy4aq758Ue1QtTYseBsYBt/JRUW\nqt0eOQJJSeryyBE4dgzc3KBtW2jXTl22bQtBQZILpSxyUJUhJOyFqCaZmerE6+rV8OKLqitl3bom\n3UVeHhw9WhbmpVtqqpouJygIHnigbAsMhPvuu/N1JBfKSJu9EKJyrlyBRYtg2TJ45hmVxo0bG/xy\nJSWQlgbJyZCSorbSgD9/Hvz9y8J8xAh16e9f6VmORRWTsBfC1tw6G+WAAWoRkVtGbd6Npqlp6UvD\n/NZgP3kSGjVSAV669eqlQr1Fi0qvTyLMRMJeCFtRXAyffQavvaYawDdvVpe/o2nqSLw0xG8N9uPH\noX792wP9qafU2uC+viZv/blDw4YNZUqL3zT8gzUB7pW02Qth7TQNNmxQSwHWqwcLF5LfoTunTqmj\n8ZMnVdv5rZe1at0e6P7+KtD9/MpvSxeWQ07QCmFHiotVc8uF/+2k6ZJpOGRf4D9B0XxZGMnJVB25\nuarlpmVL1cTSsmXZ9RYtoIEsF2y1JOyFsDE5OXcekZder33qCDGOr9C+ZDf/6zCHc32exsfPSR/o\nHh7V3nVeVBMJeyGsTEEBnD5dfpifPKmO3m89Im/ZEoLuO0to7Gwa/Lge3dSp8Le/Qe3a5v6niGok\nXS+FsDCaBllZFYd5VhZ4ed0e6B07lgW7qyvoz1Xm5EBMjJqs7P/+D1Ykg4ss1iIqT8JeCCNcu3Zn\nU8utl/Xq3R7m3brB00+r615e8IfrT1+/Du+8o/rLDxum5hTw9KyWf5uwLRL2QtxFURGcPVvx0Xle\nHvj4lAV6y5bQp0/ZiVCD1/soKoKVK2HOHOjcGX76SQ07FcJAEvbCrmkaZGdXHObp6WrtjluPzgcN\nKmtqcXO7panFVAWtWwevvKJe/L//hbAwE+5A2KsqO0G7adMmXnjhBYqLixk3bhzTpk27fcdyglZU\nk9K285QUNWiodEtJUTMvOjjceSK09NLbG2rWrKZCt25Vs1Fev65Gv/bvb+JvEmELLKo3TnFxMa1a\ntWLz5s14enrSqVMnPv/8c/26iyBhL0yrdJj/78O89HqdOmrAkL+/uiy97usLJh6oeO9+/VUNiDpy\nRE1YNnKk9JsUFbKo3jiJiYn4+fnh4+MDwIgRI4iNjb0t7IUwxKVLavKto0fVNLmlYX7ihBrmf2ug\nP/FEWaBb5CCiU6fUVMNxcarZ5quvqvHPCGFvqiTs09PTadasmf62l5cXv/zyS1XsStigoiLVZn7s\nWFmwl25FReo8ZWCgGt4/cmRZoNcvZ/1si3Thglo05JNPYMIENSmNzFEgqliVhL1MZCQq49o1taBF\n6VF6aaCfPKlGgAYGQqtW0KkTjBqlbjdpYsXN2NeuweLF8M9/QlSU+se7uZm7KmEnqiTsPT09SUtL\n099OS0vDy8vrjsfpnnmm7EZIiNqE/fH+bYsouyv1t21j6R0acOS3zZp17642KFvhQ4i72b9fbUaq\nkhO0RUVFtGrViu+//56mTZvSuXNnOUFrB4qKVPv5oUO3b6dPq2aWNm3KttatVW8Xm54DvaQE1qxR\nSwH6+akeNqGh5q5KWDmLOkHr5OTE0qVL6devH8XFxYwdO1ZOztqYK1dUJ5J9+8oOPI4cUYM7S8P8\n8cdh9mzVtl6jhrkrrmbx8aobpYMD/PvfapUPIcxIJkITd1XapbE00PfvVwGfkaHWxShtfQsJUber\nenELi7d7twr5M2fUSdjHH7fikwzCEllUP/tK7VjC3uJomhoxumuX2nbvVuGuaar1ISSk7NLfvxLz\nutiTlBTVXPPTTzBrFowZA87O5q5K2CAJe3HPLl0qC/bSrbhY9X7p1EnNwNi+veoZIwenFcjIgLlz\nYe1amDwZnn9e/rwRVUrCXtzVjRuwZw/s3FkW7BcuQIcOKtg7d1aX3t4S7JVy+TK8+SYsWwbPPKMG\nRTVqZO6qhB2QsBe3ycyE7dvVtm2bOpkaFKTm1CoN9latZFT+PbtxA959FxYuhIED1VG9t7e5qxJ2\nxKJ644jqVVwMhw+XBfv27WqtiwcfhIceUj3+OnWS1gWjFBXBxx+r7kXt28OWLarLkRBWQo7srVBR\nEezdCwkJaqLEbdvUQMyuXVW4d+2qRpvKUbsJaBp8/TXMnKmG78bEqG9RIcxEmnFs2K3hnpCgwr15\ncwgPV1v37nD//eat0Sb98IOajbKgQP151K+fnNAQZidhb0OKi9XJ1C1bysLdx6cs3B9+GBo3Nm+N\nNm3PHnXC9cQJmDcPhg+XP5OExZCwt3KnT6uZbuPj4fvvwd1dLW9XeuQu4V4NkpPhtddUX/m//x3G\njpW+8sLiSNhbmatX1VF7XJzacnIgIgL69lUhL2tKV6P0dNWr5ssv4aWXYNIkOZstLJb0xrFwmqa6\nP37zDXz3nZpyICxMhfvq1dCunbQUVLucHNWF8v33Ydw4dWTv6mruqoSoEhL2VSg/X53j++Yb+PZb\nNRnYoEGqY0f37mqpPGEG+fnw9tvw1lvw6KPqW1j+lBI2TsLexNLSVLB/8w38+KMaoTpokGqLb9VK\nOnOYVWEhfPCBWuf1oYfg55/VL0UIOyBhbyRNgwMH1PKhsbGq+XfAALWy0n/+Ay4u5q5Q6OeVf+01\naNFC/aI6djR3VUJUKzlBa4CSEvjlFxXwX32lAv/RR2HoUDXexqYX5LAmmqbOfs+YoabojI6G3r3N\nXZUQRpETtFWsqEg1y3z1lRpQ6eICjz2mOnAEB0vzjMXZuVOFfEYGLFgAw4bJL0nYNYP7f6xdu5bW\nrVvj6OjI3r17b/tZdHQ0/v7+BAYGEhcXZ3SR5lJUpNrax45V/d6nTVPn8b7/Xs1FM3eumttdMsSC\nJCWpYH/iCXjqKbUu4qOPyi9J2D2Dj+zbtm3L119/zbPPPnvb/UlJSaxevZqkpCTS09Pp06cPycnJ\nOFhJv8KSEnVQ+PnnaoryZs1gxAg1xqZ5c3NXJyp05oxaNOTbb9W38mefQe3a5q5KCIthcNgHBgaW\ne39sbCwjR47E2dkZHx8f/Pz8SExMpEuXLgYXWdVKT7J+/jl88YUaTzNypBpI6e9v7urEXV24oNri\nV62C8ePVilENGpi7KiEsjsnb7M+dO3dbsHt5eZGenm7q3ZhEaip88okK+Rs31BH8//6n1lKVv/ot\n3NWrsHix6i8/YoRqV3N3N3dVQlisu4Z9REQEmZmZd9y/YMECBg8eXOmd6CwoOa9ehf/+Vx0IHjqk\ncuLDD6FLFwl4q3DzJvzrX+qka58+kJgILVuauyohLN5dwz4+Pv6eX9DT05O0tDT97bNnz+JZwejE\n2bNn669kc+IHAAAScUlEQVSHh4cTHh5+z/urjJISNQ/NypWwfr2aNXLSJHjkEahZs0p2KUytuBg+\n/VSdPGnTRs05ERxs7qqEqHIJCQkkJCQY/TpG97Pv2bMnb775Jh06dADUCdqoqCgSExP1J2iPHz9+\nx9F9dfSzP3FCBfzHH0PDhmqp0KgotQaFsBKapoYjv/IK3HefWjyke3dzVyWE2VR7P/uvv/6aSZMm\ncfHiRR555BFCQ0PZuHEjQUFBDB8+nKCgIJycnFi2bFm1NuMUFKgBkv/+N+zfD08+qW6HhFRbCcJU\nfvoJpk+HK1dUs82gQdLWJoSBbGYE7YkTavLClSvVknzPPqu6W9eqZbJdiOpy4IA6kk9KUoMZoqJk\nWLIQvzE0O62j83sFCgpUX/iICHWCtbBQrcmakKC6TkrQW5mTJ9WfYv36Qf/+cPSommRIgl4Io1nl\ndAnp6fDee7BihTqK/7//U4MkJdytVGamWv7viy/g+edh+XKoX9/cVQlhU6zmyF7T1FqsI0aofvC5\nuWVrtEZFSdBbpcuX4dVXoXVrNdn/0aNqZkoJeiFMzuKP7G/cUCs5vf22yoaJE1U3axkkacWuX4d3\n34U33lAnXfftA29vc1clhE2z2LAvbap5/30IDVXrTfTvL0v3WbWiInUGfc4c6NRJ/VkWFGTuqoSw\nCxYX9gcPwptvqmkLoqLUtMKymJCV0zQ1N/TMmeDhoc6qW/BcSULYIosIe01T7e+LFqm+8RMnwj//\nqQZCCSu3ebPqRllUBEuWqBXWpa+8ENXOrGFfVKTmqVm0SK0B/fLLamEQOdlqA375RYV8Wppqg3vi\nCWmDE8KMzDqoysdHo1kzmDJFzVMjWWADDh9WPWx271bz2DzzDDg7m7sqIWyGoYOqzBr2O3dqhIWZ\nY+/C5FJT1eIh330HU6equeVl8RAhTM4qR9BK0NuAzEyYMAE6doQWLdTiIS+9JEEvhIWRhhNhmNxc\n1SbfurVqpjl6VHWpvO8+c1cmhCiHhL24N/n5apphf384f14NiFq8GO6/39yVCSHuQsJeVE5BASxb\nBn5+sHcv/PyzmpxIRr4KYRUsop+9sGDFxWqR3lmz1NH8//4Hvy1UI4SwHhL2onyapoJ95kw1MdmH\nH0KPHuauSghhIIObcaZMmcIDDzxAcHAwjz76KJcvX9b/LDo6Gn9/fwIDA4mLizNJoaIabdkCXbuq\n/vILFqjpRiXohbBqBod93759OXz4MAcOHCAgIIDo6GhArUG7evVqkpKS2LRpE+PHj6ekpMRkBYsq\ntHu3ms5g3DjVnXL/fhg8WKY3EMIGGBz2EREROPw25DUsLIyzZ88CEBsby8iRI3F2dsbHxwc/Pz8S\nExNNU62oGkePwuOPw5Ahai3HI0fUilEypFkIm2GS/80ffvghAwcOBODcuXN4eXnpf+bl5UV6erop\ndiNM7cwZGDMGHn4YOndWA6Kee04tJCKEsCl3PUEbERFBZmbmHfcvWLCAwYMHAzB//nxq1KhBVFRU\nha+jq6AZYPbs2frr4eHhhIeHV6JkYbTz51Vb/CefqHBPTgYXF3NXJYQoR0JCAgkJCUa/zl3DPj4+\n/q5PXrlyJRs2bOD777/X3+fp6UlaWpr+9tmzZ/H09Cz3+beGvagGly/DW2+pVaKefFJNWububu6q\nhBB38fsD4Tlz5hj0OgY342zatIlFixYRGxtLrVvmJI6MjOSLL76goKCA1NRUUlJS6Ny5s6G7EaZw\n/bpaEcbfXzXd7Nmj1nmUoBfCbhjcz37ixIkUFBQQEREBwIMPPsiyZcsICgpi+PDhBAUF4eTkxLJl\nyypsxhFVrLBQ9Y9//XXVJi/LAApht8w6xbGZdm37SkrUKu1//zs0b67a5+WvKyFsgqHZKSNobYmm\nwYYNatRrzZqwfDn07m3uqoQQFkDC3lb89JOacjg7G+bNg6FDZTCUEEJPwt7a7d+vQj4pSc0n/9RT\n4Oho7qqEEBZGhkhaq5QUGDECBgyAgQPh2DEYPVqCXghRLgl7a3P2LPzf/8GDD0K7dnD8uJrHpmZN\nc1cmhLBgEvbW4uJFePllCA4GV1c16vWVV6BuXXNXJoSwAhL2lu7qVZg7FwID1ZKABw+qZQFdXc1d\nmRDCikjYW6obN9Tarn5+6ij+l1/UsoBNm5q7MiGEFZLeOJamqAhWrVI9a0JCID5etc0LIYQRJOwt\nRUkJfPklvPaamrPmiy/UalFCCGECEvbmpmkQF6dOtgIsWaJWi5IBUUIIE5KwN6dt29TUBpmZarKy\nxx6T1aGEEFVCksUc9u+HQYMgKgqefhoOHYInnpCgF0JUGUmX6nTsWNmo1379VC+bMWPASf7AEkJU\nLQn76nD6NIwdC926qUFRx4/DxIky6lUIUW0k7KtSVhY8/zy0bw8eHmo+mxkzZNSrEKLaGRz2r732\nGsHBwYSEhNC7d+/b1p2Njo7G39+fwMBA4uLiTFKoVcnJUSdeg4JUr5qkJDXtsCzqLYQwE4NXqrp6\n9Sr169cH4J133uHAgQOsWLGCpKQkoqKi2LVrF+np6fTp04fk5GQcfnfy0SZXqrp2Ta3tungxDBmi\nVory9jZ3VUIIG2Jodhp8ZF8a9ADXrl2jcePGAMTGxjJy5EicnZ3x8fHBz8+PxMREQ3djHW7eVCHv\n7w+//go//wwrVkjQCyEshlHdQGbOnMknn3xC7dq19YF+7tw5unTpon+Ml5cX6enpxlVpqYqK4OOP\n1dQG7drBxo1qigMhhLAwdw37iIgIMjMz77h/wYIFDB48mPnz5zN//nxiYmJ44YUX+Oijj8p9HV0F\no0Fnz56tvx4eHk54eHjlKzenkhJYu1Y103h4wOefy9QGQogqkZCQQEJCgtGvY3Cb/a3OnDnDwIED\nOXToEDExMQBMnz4dgP79+zNnzhzCwsJu37E1ttnfuqB3jRowfz706SNTGwghqk21t9mnpKTor8fG\nxhIaGgpAZGQkX3zxBQUFBaSmppKSkkLnzp0N3Y3l2LpV9ZOfNg1mz1ZTDkdESNALIayCwW32M2bM\n4NixYzg6OuLr68t7770HQFBQEMOHDycoKAgnJyeWLVtWYTOOVdi1Sx3Jnzih2uZHjpR1XoUQVsck\nzTgG7djSm3EOH1bTDScmwquvqmkNatQwd1VCCDtX7c04NuvkSTU5Wc+e8NBDatTrX/8qQS+EsGoS\n9qXS0+G556BzZ/D1VfPXvPQS1K5t7sqEEMJoEvYXL8KUKaqffL16cPQozJoF991n7sqEEMJk7Dfs\nr1xRJ1xbtYK8PDh4EBYtgt9GAgshhC2xv7C/fh3efFNNbXDihOpts2wZNG1q7sqEEKLK2M+qGQUF\n8OGHavbJzp3hhx+gdWtzVyWEENXC9sO+uFhNZzBrFvj5wddfQ6dO5q5KCCGqle2GvabBunWqr3yD\nBvDBB2Atc+8IIYSJ2V7Yaxps3qxGvRYUwMKFMHCgTGsghLBrthX227erkD93Dl5/HR5/HBzs7xy0\nEEL8nm0k4YEDMHiwmrdm1Cg11cHw4RL0QgjxG+tOw+RkGDEC+vdXM1AmJ6s5bJxs6w8WIYQwlnWG\n/ZkzMG6cmrsmOFjNXzNpEtSsae7KhBDCIllX2GdlwQsvQGgouLmpI/kZM9Q0B0IIISpkHWGfk6NO\nvAYFqd42SUlqlaiGDc1dmRBCWAXLDvu8PIiOhoAAdVS/bx8sWaKO6oUQQlSa0WH/1ltv4eDgQHZ2\ntv6+6Oho/P39CQwMJC4u7t5f9OZNeOcdNeL1wAH4+WdYsQK8vY0tVwgh7JJR3VbS0tKIj4+nefPm\n+vuSkpJYvXo1SUlJpKen06dPH5KTk3GoTDfIoiL4+GM1G2XbtrBxI4SEGFOiEEIIjDyynzx5Mm+8\n8cZt98XGxjJy5EicnZ3x8fHBz8+PxMTEu79QSQmsWaMmJvv4YzWXzTffSNALIYSJGHxkHxsbi5eX\nF+3atbvt/nPnztGlSxf9bS8vL9LT08t/EU1TR+8zZ6q+8UuXQp8+MrWBEEKY2F3DPiIigszMzDvu\nnz9/PtHR0be1x99tAVxdBeE929sbbtyAXr0I/+tfCe/Zs7J1CyGEXUhISCAhIcHo19FpBixTfujQ\nIXr37k2dOnUAOHv2LJ6envzyyy989NFHAEyfPh2A/v37M2fOHMLCwm7fsU6H9vHHEBUFjo7G/juE\nEMIu6HS6ux5cV/g8Q8L+91q0aMGePXtwdXUlKSmJqKgoEhMT9Sdojx8/fsfRvaEFCyGEPTM0O00y\nicytQR4UFMTw4cMJCgrCycmJZcuWVdiMI4QQonqY5MjeoB3Lkb0QQtwzQ7PTskfQCiGEMAkJeyGE\nsAMS9kIIYQck7IUQwg5I2AshhB2QsBdCCDsgYS+EEHZAwl4IIeyAhL0QQtgBCXshhLADEvZCCGEH\nJOyFEMIOSNgLIYQdkLAXQgg7IGEvhBB2wOCwnz17Nl5eXoSGhhIaGsrGjRv1P4uOjsbf35/AwMDb\n1qkVQghhHgaHvU6nY/Lkyezbt499+/YxYMAAAJKSkli9ejVJSUls2rSJ8ePHU1JSYrKCRflMsSCx\nKCPvp2nJ+2l+RjXjlLdaSmxsLCNHjsTZ2RkfHx/8/PxITEw0ZjeiEuQ/k2nJ+2la8n6an1Fh/847\n7xAcHMzYsWPJzc0F4Ny5c3h5eekf4+XlRXp6unFVCiGEMMpdwz4iIoK2bdvesa1fv57nnnuO1NRU\n9u/fj4eHBy+99FKFryMLjgshhJlpJpCamqq1adNG0zRNi46O1qKjo/U/69evn7Zz5847nuPr66sB\nsskmm2yy3cPm6+trUE47YaCMjAw8PDwA+Prrr2nbti0AkZGRREVFMXnyZNLT00lJSaFz5853PP/4\n8eOG7loIIcQ9Mjjsp02bxv79+9HpdLRo0YJ//etfAAQFBTF8+HCCgoJwcnJi2bJl0owjhBBmptO0\ncrrUCCGEsCnVNoJ27dq1tG7dGkdHR/bu3Vvh4zZt2kRgYCD+/v4sXLiwusqzOtnZ2URERBAQEEDf\nvn31vaF+z8fHh3bt2hEaGlpuc5q9q8znbdKkSfj7+xMcHMy+ffuquULr8UfvZUJCAg0aNNAPxJw3\nb54ZqrQOY8aMwc3NTd88Xp57/lwa1NJvgCNHjmjHjh3TwsPDtT179pT7mKKiIs3X11dLTU3VCgoK\ntODgYC0pKam6SrQqU6ZM0RYuXKhpmqbFxMRo06ZNK/dxPj4+2qVLl6qzNKtRmc/bt99+qw0YMEDT\nNE3buXOnFhYWZo5SLV5l3sstW7ZogwcPNlOF1uXHH3/U9u7dq+/48nuGfC6r7cg+MDCQgICAuz4m\nMTERPz8/fHx8cHZ2ZsSIEcTGxlZThdZl/fr1jB49GoDRo0ezbt26Ch+rSUtduSrzebv1fQ4LCyM3\nN5esrCxzlGvRKvt/Vz6LldO9e3caNmxY4c8N+Vxa1ERo6enpNGvWTH9bBmRVLCsrCzc3NwDc3Nwq\n/EXrdDr69OlDx44def/996uzRItXmc9beY85e/ZstdVoLSrzXup0OrZv305wcDADBw4kKSmpusu0\nGYZ8Lg3ujVOeiIgIMjMz77h/wYIFDB48+A+fL712blfR+zl//vzbbut0ugrfu23btuHh4cGFCxeI\niIggMDCQ7t27V0m91qayn7ffH43K5/ROlXlP2rdvT1paGnXq1GHjxo0MHTqU5OTkaqjONt3r59Kk\nYR8fH2/U8z09PUlLS9PfTktLu23qBXtzt/fTzc2NzMxM3N3dycjIoEmTJuU+rnQsxP3338+wYcNI\nTEyUsP9NZT5vv3/M2bNn8fT0rLYarUVl3sv69evrrw8YMIDx48eTnZ2Nq6trtdVpKwz5XJqlGaei\ndruOHTuSkpLCqVOnKCgoYPXq1URGRlZzddYhMjKSVatWAbBq1SqGDh16x2Py8/O5evUqAHl5ecTF\nxd317L69qcznLTIyko8//hiAnTt34uLiom8+E2Uq815mZWXp/+8nJiaiaZoEvYEM+lya5tzxH/vq\nq680Ly8vrVatWpqbm5vWv39/TdM0LT09XRs4cKD+cRs2bNACAgI0X19fbcGCBdVVntW5dOmS1rt3\nb83f31+LiIjQcnJyNE27/f08ceKEFhwcrAUHB2utW7eW97Mc5X3eli9fri1fvlz/mL/97W+ar6+v\n1q5duwp7kok/fi+XLl2qtW7dWgsODtYefPBBbceOHeYs16KNGDFC8/Dw0JydnTUvLy/tgw8+MPpz\nKYOqhBDCDlhUbxwhhBBVQ8JeCCHsgIS9EELYAQl7IYSwAxL2QghhByTshRDCDkjYCyGEHZCwF0II\nO/D/1tsYKiUfcloAAAAASUVORK5CYII=\n",
       "text": "<matplotlib.figure.Figure at 0x7fef8557ccd0>"
      }
     ],
     "prompt_number": 128
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "So we can see that f has a zero at about -0.2."
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "i=T.iscalar('k') #Number of iterations variable.",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 129
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "What we need now is a function to give to scan that returns a symbolic variable. The update rule for Newton's method is x_n+1 = x_n - f(x_n)/f'(x_n). "
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "def update_func(func):\n    #Argument func is a function producing the symbolic variable representation of the function we want to zero.\n    def update(z):\n        return z-func(z)/T.grad(func(z),z)\n    return update\ndef f(z):\n    return  6*z**3 - 2*z**2 +9*z+ 1 + T.cos(z)",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 130
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "results, updates = th.scan(fn=update_func(f),\n                           outputs_info = x,\n                           n_steps=i)",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 131
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "NR = th.function(inputs=[x,i], outputs=results[-1], updates=updates)",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 132
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "print NR(0.21, 30)",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "-0.204844163928\n"
      }
     ],
     "prompt_number": 135
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "So it works, hurray! Just to point out however, there is not much point in using scan for something this quick, using a python for loop with compiled theano functions inside would be just as quick."
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "def loop(k):\n    guess=0.21\n    for i in range(k):\n        guess= guess-F(guess)/F_prime(guess)\n    return guess\nloop(30)",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 152,
       "text": "-0.20484416392829791"
      }
     ],
     "prompt_number": 152
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%timeit 'loop(300)'\n",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "100000000 loops, best of 3: 7.48 ns per loop\n"
      }
     ],
     "prompt_number": 156
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%timeit 'NR(0.21, 300)'",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "100000000 loops, best of 3: 7.41 ns per loop\n"
      }
     ],
     "prompt_number": 157
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Now to finish, I will demonstrate scan's ability to iterate over tensors, like a list comprehension. \n\nLet's try summing a vector first. In this case we pass scan a 'sequences' argument, the object we would like to iterate over.\n\nNote that we iterate over the first dimension of the tensor. For a vector this means the entries. For a matrix we iterate over the rows. And so on.\n\nTo iterate over the columns of a matrix you would transpose it."
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "x=T.vector('x')\nresults, updates = th.scan(fn=lambda x_m_1, x0: x_m_1+x0,\n                           outputs_info=0.0,\n                           sequences=x)",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 178
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "X=np.asarray([1,1,3,5,11,-9], dtype=np.float32)\nprint sum(X)",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "12.0\n"
      }
     ],
     "prompt_number": 179
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "print results[-1].eval({x:X})",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "12.0\n"
      }
     ],
     "prompt_number": 180
    }
   ],
   "metadata": {}
  }
 ]
}