{
 "metadata": {
  "language": "Julia",
  "name": "",
  "signature": "sha256:9ffa086197a683b322822329bf4c22052357265c43076a8ea7a2c46e4f811b81"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Metaprogramming in Julia"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Metaprogramming** is writing **code that writes code**.\n",
      "\n",
      "Inspired by several other languages, notably Scheme, Julia provides *built-in* facilities for metaprogramming:\n",
      "\n",
      "* Julia provides access to its parser and [abstract syntax tree](http://en.wikipedia.org/wiki/Abstract_syntax_tree): You can get a *symbolic representation* of any Julia code.\n",
      "\n",
      "* You can manipulate these symbolic representations to transform and generate Julia code at runtime, and *evaluate* it to run the resulting code.\n",
      "\n",
      "* Julia provides **symbolic macros**: these are essentially functions evaluated at *parse time* which take the *syntax tree* of the code, perform arbitrary transformations, and insert new code to be later compiled.\n",
      "\n",
      "Julia macros, inspired by Scheme's [hygienic macros](http://en.wikipedia.org/wiki/Hygienic_macro), effectively allow you to both **extend the syntax of Julia** with arbitrary parse-time code generation."
     ]
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Symbolic expressions in Julia"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ex = x - 2y"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "LoadError",
       "evalue": "x not defined\nwhile loading In[1], in expression starting on line 1",
       "output_type": "pyerr",
       "traceback": [
        "x not defined\nwhile loading In[1], in expression starting on line 1",
        ""
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Using `:(.....)` or `quote .... end` produces a **symbolic expression** of type `Expr`, which contains the **parsed syntax tree** of a Julia expression."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ex = :(x - 2y)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 2,
       "text": [
        ":(x - 2y)"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The `dump` function uses **introspection** to print the contents of any data structure:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dump(ex)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Expr "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "  head: Symbol call\n",
        "  args: Array(Any,(3,))\n",
        "    1: Symbol -\n",
        "    2: Symbol x\n",
        "    3: Expr \n",
        "      head: Symbol call\n",
        "      args: Array(Any,(3,))\n",
        "        1: Symbol *\n",
        "        2: Int64 2\n",
        "        3: Symbol y\n",
        "      typ: Any\n",
        "  typ: Any\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Macros in Julia"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Essentially **functions evaluated at parse-time**, which take a **symbolic expression** as input and produce **another expression** as output, which is **inserted into the code** before compilation:"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "parse \u2192 expressions \u2192 macro \u2192 new expr. \u2192 compile "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "A simple macro example: **reverse the order of function arguments**"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "macro reverse(ex)\n",
      "    if isa(ex, Expr) && ex.head == :call\n",
      "        return Expr(:call, ex.args[1], reverse(ex.args[2:end])...)\n",
      "    else\n",
      "        return ex\n",
      "    end\n",
      "end"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# equivalent to 3 - 1\n",
      "@reverse 1 - 4"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "3"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "A useful macro: Polynomial evaluation by Horner"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The following macro evaluates the polynomial\n",
      "\n",
      "* $p(x) = c_0 + c_1 x + \\cdots + c_n x^n$\n",
      "\n",
      "by **Horner's rule**\n",
      "\n",
      "* $c_0 + x \\cdot (c_1 + x \\cdot (c_2 + x \\cdot (c_3 + \\cdots)))$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "macro horner(x, c...)\n",
      "    ex = esc(c[end])\n",
      "    for i = length(c)-1:-1:1\n",
      "        ex = :($(esc(c[i])) + t * $ex)\n",
      "    end\n",
      "    return Expr(:block, :(t = $(esc(x))), ex)\n",
      "end"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Special-function evaluation"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Fast inline polynomial evaluation is very useful for **special functions**.  For example, evaluating the inverse $\\mathrm{erf}^{-1}(x)$ of the **error function**:\n",
      "\n",
      "* $\\mathrm{erf}(x) = \\frac{2}{\\sqrt{pi}} \\int_0^x e^{-t^2} dt$\n",
      "\n",
      "via **rational approximants** (ratios of polynomials) [Blair (1976)]:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "function my_erfinv(x::Float32) # specialized for single-precision args\n",
      "    a = abs(x)\n",
      "    if a >= 1.0f0\n",
      "        if x == 1.0f0\n",
      "            return inf(Float32)\n",
      "        elseif x == -1.0f0\n",
      "            return -inf(Float32)\n",
      "        end\n",
      "        throw(DomainError())\n",
      "    elseif a <= 0.75f0 # Table 10 in Blair et al.                               \n",
      "        t = x*x - 0.5625f0\n",
      "        return x * @horner(t, -0.13095_99674_22f2,\n",
      "                               0.26785_22576_0f2,\n",
      "                              -0.92890_57365f1) /\n",
      "                   @horner(t, -0.12074_94262_97f2,\n",
      "                               0.30960_61452_9f2,\n",
      "                              -0.17149_97799_1f2,\n",
      "                               0.1f1)\n",
      "    elseif a <= 0.9375f0 # Table 29 in Blair et al.                             \n",
      "        t = x*x - 0.87890625f0\n",
      "        return x * @horner(t, -0.12402_56522_1f0,\n",
      "                               0.10688_05957_4f1,\n",
      "                              -0.19594_55607_8f1,\n",
      "                               0.42305_81357f0) /\n",
      "                   @horner(t, -0.88276_97997f-1,\n",
      "                               0.89007_43359f0,\n",
      "                              -0.21757_03119_6f1,\n",
      "                               0.1f1)\n",
      "    else # Table 50 in Blair et al.                                             \n",
      "        t = 1.0f0 / sqrt(-log(1.0f0 - a))\n",
      "        return @horner(t, 0.15504_70003_116f0,\n",
      "                          0.13827_19649_631f1,\n",
      "                          0.69096_93488_87f0,\n",
      "                         -0.11280_81391_617f1,\n",
      "                          0.68054_42468_25f0,\n",
      "                         -0.16444_15679_1f0) /\n",
      "              (copysign(t, x) *\n",
      "               @horner(t, 0.15502_48498_22f0,\n",
      "                          0.13852_28141_995f1,\n",
      "                          0.1f1))\n",
      "    end\n",
      "end\n",
      "@vectorize_1arg Real my_erfinv"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "my_erfinv (generic function with 4 methods)"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This is precisely how `erfinv` is implemented in Julia (in single and double precision), and is ***3\u20134\u00d7 faster*** than the compiled (Fortran?) code in Matlab, and ***2\u20133\u00d7 faster*** than the compiled (Fortran Cephes) code used in SciPy.\n",
      "\n",
      "The difference (at least in Cephes) seems to be mainly that **they have explicit arrays of polynomial coefficients** and call a **subroutine** for Horner's rule, versus inlining it via a macro."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "using PyCall\n",
      "@pyimport scipy.special as s\n",
      "x = rand(10^7);\n",
      "@time erfinv(x);\n",
      "@time s.erfinv(x);\n",
      "norm(erfinv(x) - s.erfinv(x)) / norm(erfinv(x))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "elapsed time: 0."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "2604303 seconds (80000160 bytes allocated)\n",
        "elapsed time: "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.464521205 seconds (80003272 bytes allocated, 4.88% gc time)\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 9,
       "text": [
        "1.2616095260735443e-15"
       ]
      }
     ],
     "prompt_number": 9
    }
   ],
   "metadata": {}
  }
 ]
}