02_tutorial_basic_kernels.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pystencils.session import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 02: Basic Kernel generation with *pystencils*\n",
    "\n",
    "Now that you have an [overview of pystencils](01_tutorial_getting_started.ipynb), \n",
    "this tutorial shows in more detail how to formulate, optimize and run stencil kernels.\n",
    "\n",
    "## 1) Kernel Definition\n",
    "\n",
    "### a) Defining kernels with assignment lists and the  `kernel` decorator \n",
    "\n",
    "*pystencils* gets a symbolic formulation of the kernel. This can be either an `Assignment` or a sequence of `Assignment`s that follow a set of restrictions. \n",
    "\n",
    "Lets first create a kernel that consists of multiple assignments:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "src_arr = np.zeros([20, 30])\n",
    "dst_arr = np.zeros_like(src_arr)\n",
    "\n",
    "dst, src = ps.fields(dst=dst_arr, src=src_arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left [ grad_{x} \\leftarrow \\frac{{{src}_{E}}}{2} - \\frac{{{src}_{W}}}{2}, \\quad grad_{y} \\leftarrow \\frac{{{src}_{N}}}{2} - \\frac{{{src}_{S}}}{2}, \\quad {{dst}_{C}} \\leftarrow grad_{x} + grad_{y}\\right ]$$"
      ],
      "text/plain": [
       "⎡         src_E   src_W            src_N   src_S                         ⎤\n",
       "⎢gradₓ := ───── - ─────, grad_y := ───── - ─────, dst_C := gradₓ + grad_y⎥\n",
       "⎣           2       2                2       2                           ⎦"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grad_x, grad_y = sp.symbols(\"grad_x, grad_y\")\n",
    "\n",
    "symbolic_description = [\n",
    "    ps.Assignment(grad_x, (src[1, 0] - src[-1, 0]) / 2),\n",
    "    ps.Assignment(grad_y, (src[0, 1] - src[0, -1]) / 2),\n",
    "    ps.Assignment(dst[0, 0], grad_x + grad_y),\n",
    "]\n",
    "kernel = ps.create_kernel(symbolic_description)\n",
    "symbolic_description"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We created subexpressions, using standard sympy symbols on the left hand side, to split the kernel into multiple assignments. Defining a kernel using a list of `Assignment`s is quite tedious and hard to read. \n",
    "To simplify the formulation of a kernel, *pystencils* offers the `kernel` decorator, that transforms a normal Python function with `@=` assignments into an assignment list that can be passed to `create_kernel`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left [ grad_{x} \\leftarrow \\frac{{{src}_{E}}}{2} - \\frac{{{src}_{W}}}{2}, \\quad grad_{y} \\leftarrow \\frac{{{src}_{N}}}{2} - \\frac{{{src}_{S}}}{2}, \\quad {{dst}_{C}} \\leftarrow grad_{x} + grad_{y}\\right ]$$"
      ],
      "text/plain": [
       "⎡         src_E   src_W            src_N   src_S                         ⎤\n",
       "⎢gradₓ := ───── - ─────, grad_y := ───── - ─────, dst_C := gradₓ + grad_y⎥\n",
       "⎣           2       2                2       2                           ⎦"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@ps.kernel\n",
    "def symbolic_description_using_function():\n",
    "    grad_x @= (src[1, 0] - src[-1, 0]) / 2\n",
    "    grad_y @= (src[0, 1] - src[0, -1]) / 2\n",
    "    dst[0, 0] @= grad_x + grad_y\n",
    "symbolic_description_using_function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The decorated function can contain any Python code, only the `@=` operator, and the ternary inline `if-else` operator have different meaning. \n",
    "\n",
    "### b) Ternary 'if' with `Piecewise`\n",
    "\n",
    "The ternary operator maps to `sympy.Piecewise` functions, that can be used to introduce branching into the kernel. Piecewise defined functions must give a value for every input, i.e. there must be a 'otherwise' clause in the end that is indicated by the condition `True`. Piecewise objects are standard sympy terms that can be integrated into bigger expressions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$${{src}_{E}} + \\begin{cases} 1.0 & \\text{for}\\: {{src}_{N}} > 0 \\\\0.0 & \\text{otherwise} \\end{cases}$$"
      ],
      "text/plain": [
       "        ⎛⎧1.0  for src_N > 0⎞\n",
       "src_E + ⎜⎨                  ⎟\n",
       "        ⎝⎩0.0    otherwise  ⎠"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.Piecewise((1.0, src[0,1] > 0), (0.0, True)) + src[1, 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Piecewise objects are created by the `kernel` decorator for ternary if-else statements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left [ grad_{x} \\leftarrow \\begin{cases} \\frac{{{src}_{E}}}{2} - \\frac{{{src}_{W}}}{2} & \\text{for}\\: {{src}_{W}} > 0 \\\\0.0 & \\text{otherwise} \\end{cases}\\right ]$$"
      ],
      "text/plain": [
       "⎡         ⎧src_E   src_W               ⎤\n",
       "⎢         ⎪───── - ─────  for src_W > 0⎥\n",
       "⎢gradₓ := ⎨  2       2                 ⎥\n",
       "⎢         ⎪                            ⎥\n",
       "⎣         ⎩     0.0         otherwise  ⎦"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@ps.kernel\n",
    "def kernel_with_piecewise():\n",
    "    grad_x @= (src[1, 0] - src[-1, 0]) / 2 if src[-1, 0] > 0 else 0.0\n",
    "kernel_with_piecewise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### c) Assignment level optimizations using `AssignmentCollection`\n",
    "\n",
    "When the kernels get larger and more complex, it is helpful to organize the list of assignment into a more structured way. The `AssignmentCollection` offers optimizating transformation on a list of assignments. It holds two assignment lists, one for subexpressions and one for the main assignments. Main assignments are typically those that write to an array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>Subexpressions:</div><table style=\"border:none; width: 100%; \"><tr style=\"border:none\"> <td style=\"border:none\">$$a \\leftarrow {{src}_{N}} + {{src}_{W}}$$</td>  </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$b \\leftarrow 2 {{src}_{E}} + {{src}_{S}}$$</td>  </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$c \\leftarrow - {{src}_{C}} + 2 {{src}_{E}} + {{src}_{N}} + {{src}_{S}} + {{src}_{W}}$$</td>  </tr> </table><div>Main Assignments:</div><table style=\"border:none; width: 100%; \"><tr style=\"border:none\"> <td style=\"border:none\">$${{dst}_{C}} \\leftarrow a + b + c$$</td>  </tr> </table>"
      ],
      "text/plain": [
       "Equation Collection for dst_C"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@ps.kernel\n",
    "def somewhat_longer_dummy_kernel(s):\n",
    "    s.a @= src[0, 1] + src[-1, 0]\n",
    "    s.b @= 2 * src[1, 0] + src[0, -1]\n",
    "    s.c @= src[0, 1] + 2 * src[1, 0] + src[-1, 0] + src[0, -1] - src[0,0]\n",
    "    dst[0, 0] @= s.a + s.b + s.c\n",
    "    \n",
    "ac = ps.AssignmentCollection(main_assignments=somewhat_longer_dummy_kernel[-1:], \n",
    "                             subexpressions=somewhat_longer_dummy_kernel[:-1])\n",
    "ac"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'adds': 8, 'muls': 2, 'divs': 0}"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ac.operation_count"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `pystencils.simp` submodule offers several functions to optimize a collection of assignments.\n",
    "It also offers functionality to group optimization into strategies and evaluate them. \n",
    "In this example we reduce the number of operations by reusing existing subexpressions to get rid of two unnecessary floating point additions. For more information about assignment collections and simplifications see the [demo notebook](demo_assignment_collection.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>Subexpressions:</div><table style=\"border:none; width: 100%; \"><tr style=\"border:none\"> <td style=\"border:none\">$$a \\leftarrow {{src}_{N}} + {{src}_{W}}$$</td>  </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$b \\leftarrow 2 {{src}_{E}} + {{src}_{S}}$$</td>  </tr> <tr style=\"border:none\"> <td style=\"border:none\">$$c \\leftarrow - {{src}_{C}} + a + b$$</td>  </tr> </table><div>Main Assignments:</div><table style=\"border:none; width: 100%; \"><tr style=\"border:none\"> <td style=\"border:none\">$${{dst}_{C}} \\leftarrow a + b + c$$</td>  </tr> </table>"
      ],
      "text/plain": [
       "Equation Collection for dst_C"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "opt_ac = ps.simp.subexpression_substitution_in_existing_subexpressions(ac)\n",
    "opt_ac"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'adds': 6, 'muls': 1, 'divs': 0}"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "opt_ac.operation_count"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### d) Ghost layers and iteration region\n",
    "\n",
    "When creating a kernel with neighbor accesses, *pystencils* automatically restricts the iteration region, such that all accesses are safe. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>.highlight .hll { background-color: #ffffcc }\n",
       ".highlight  { background: #f8f8f8; }\n",
       ".highlight .c { color: #408080; font-style: italic } /* Comment */\n",
       ".highlight .err { border: 1px solid #FF0000 } /* Error */\n",
       ".highlight .k { color: #008000; font-weight: bold } /* Keyword */\n",
       ".highlight .o { color: #666666 } /* Operator */\n",
       ".highlight .ch { color: #408080; font-style: italic } /* Comment.Hashbang */\n",
       ".highlight .cm { color: #408080; font-style: italic } /* Comment.Multiline */\n",
       ".highlight .cp { color: #BC7A00 } /* Comment.Preproc */\n",
       ".highlight .cpf { color: #408080; font-style: italic } /* Comment.PreprocFile */\n",
       ".highlight .c1 { color: #408080; font-style: italic } /* Comment.Single */\n",
       ".highlight .cs { color: #408080; font-style: italic } /* Comment.Special */\n",
       ".highlight .gd { color: #A00000 } /* Generic.Deleted */\n",
       ".highlight .ge { font-style: italic } /* Generic.Emph */\n",
       ".highlight .gr { color: #FF0000 } /* Generic.Error */\n",
       ".highlight .gh { color: #000080; font-weight: bold } /* Generic.Heading */\n",
       ".highlight .gi { color: #00A000 } /* Generic.Inserted */\n",
       ".highlight .go { color: #888888 } /* Generic.Output */\n",
       ".highlight .gp { color: #000080; font-weight: bold } /* Generic.Prompt */\n",
       ".highlight .gs { font-weight: bold } /* Generic.Strong */\n",
       ".highlight .gu { color: #800080; font-weight: bold } /* Generic.Subheading */\n",
       ".highlight .gt { color: #0044DD } /* Generic.Traceback */\n",
       ".highlight .kc { color: #008000; font-weight: bold } /* Keyword.Constant */\n",
       ".highlight .kd { color: #008000; font-weight: bold } /* Keyword.Declaration */\n",
       ".highlight .kn { color: #008000; font-weight: bold } /* Keyword.Namespace */\n",
       ".highlight .kp { color: #008000 } /* Keyword.Pseudo */\n",
       ".highlight .kr { color: #008000; font-weight: bold } /* Keyword.Reserved */\n",
       ".highlight .kt { color: #B00040 } /* Keyword.Type */\n",
       ".highlight .m { color: #666666 } /* Literal.Number */\n",
       ".highlight .s { color: #BA2121 } /* Literal.String */\n",
       ".highlight .na { color: #7D9029 } /* Name.Attribute */\n",
       ".highlight .nb { color: #008000 } /* Name.Builtin */\n",
       ".highlight .nc { color: #0000FF; font-weight: bold } /* Name.Class */\n",
       ".highlight .no { color: #880000 } /* Name.Constant */\n",
       ".highlight .nd { color: #AA22FF } /* Name.Decorator */\n",
       ".highlight .ni { color: #999999; font-weight: bold } /* Name.Entity */\n",
       ".highlight .ne { color: #D2413A; font-weight: bold } /* Name.Exception */\n",
       ".highlight .nf { color: #0000FF } /* Name.Function */\n",
       ".highlight .nl { color: #A0A000 } /* Name.Label */\n",
       ".highlight .nn { color: #0000FF; font-weight: bold } /* Name.Namespace */\n",
       ".highlight .nt { color: #008000; font-weight: bold } /* Name.Tag */\n",
       ".highlight .nv { color: #19177C } /* Name.Variable */\n",
       ".highlight .ow { color: #AA22FF; font-weight: bold } /* Operator.Word */\n",
       ".highlight .w { color: #bbbbbb } /* Text.Whitespace */\n",
       ".highlight .mb { color: #666666 } /* Literal.Number.Bin */\n",
       ".highlight .mf { color: #666666 } /* Literal.Number.Float */\n",
       ".highlight .mh { color: #666666 } /* Literal.Number.Hex */\n",
       ".highlight .mi { color: #666666 } /* Literal.Number.Integer */\n",
       ".highlight .mo { color: #666666 } /* Literal.Number.Oct */\n",
       ".highlight .sa { color: #BA2121 } /* Literal.String.Affix */\n",
       ".highlight .sb { color: #BA2121 } /* Literal.String.Backtick */\n",
       ".highlight .sc { color: #BA2121 } /* Literal.String.Char */\n",
       ".highlight .dl { color: #BA2121 } /* Literal.String.Delimiter */\n",
       ".highlight .sd { color: #BA2121; font-style: italic } /* Literal.String.Doc */\n",
       ".highlight .s2 { color: #BA2121 } /* Literal.String.Double */\n",
       ".highlight .se { color: #BB6622; font-weight: bold } /* Literal.String.Escape */\n",
       ".highlight .sh { color: #BA2121 } /* Literal.String.Heredoc */\n",
       ".highlight .si { color: #BB6688; font-weight: bold } /* Literal.String.Interpol */\n",
       ".highlight .sx { color: #008000 } /* Literal.String.Other */\n",
       ".highlight .sr { color: #BB6688 } /* Literal.String.Regex */\n",
       ".highlight .s1 { color: #BA2121 } /* Literal.String.Single */\n",
       ".highlight .ss { color: #19177C } /* Literal.String.Symbol */\n",
       ".highlight .bp { color: #008000 } /* Name.Builtin.Pseudo */\n",
       ".highlight .fm { color: #0000FF } /* Name.Function.Magic */\n",
       ".highlight .vc { color: #19177C } /* Name.Variable.Class */\n",
       ".highlight .vg { color: #19177C } /* Name.Variable.Global */\n",
       ".highlight .vi { color: #19177C } /* Name.Variable.Instance */\n",
       ".highlight .vm { color: #19177C } /* Name.Variable.Magic */\n",
       ".highlight .il { color: #666666 } /* Literal.Number.Integer.Long */</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div class=\"highlight\"><pre><span></span><span class=\"n\">FUNC_PREFIX</span> <span class=\"kt\">void</span> <span class=\"nf\">kernel</span><span class=\"p\">(</span><span class=\"kt\">double</span> <span class=\"o\">*</span> <span class=\"n\">RESTRICT</span> <span class=\"n\">fd_dst</span><span class=\"p\">,</span> <span class=\"kt\">double</span> <span class=\"o\">*</span> <span class=\"n\">RESTRICT</span> <span class=\"k\">const</span> <span class=\"n\">fd_src</span><span class=\"p\">)</span>\n",
       "<span class=\"p\">{</span>\n",
       "   <span class=\"k\">for</span> <span class=\"p\">(</span><span class=\"kt\">int</span> <span class=\"n\">ctr_0</span> <span class=\"o\">=</span> <span class=\"mi\">2</span><span class=\"p\">;</span> <span class=\"n\">ctr_0</span> <span class=\"o\">&lt;</span> <span class=\"mi\">18</span><span class=\"p\">;</span> <span class=\"n\">ctr_0</span> <span class=\"o\">+=</span> <span class=\"mi\">1</span><span class=\"p\">)</span>\n",
       "   <span class=\"p\">{</span>\n",
       "      <span class=\"kt\">double</span> <span class=\"o\">*</span> <span class=\"n\">RESTRICT</span> <span class=\"n\">fd_dst_C</span> <span class=\"o\">=</span> <span class=\"mi\">30</span><span class=\"o\">*</span><span class=\"n\">ctr_0</span> <span class=\"o\">+</span> <span class=\"n\">fd_dst</span><span class=\"p\">;</span>\n",
       "      <span class=\"kt\">double</span> <span class=\"o\">*</span> <span class=\"n\">RESTRICT</span> <span class=\"k\">const</span> <span class=\"n\">fd_src_2E</span> <span class=\"o\">=</span> <span class=\"mi\">30</span><span class=\"o\">*</span><span class=\"n\">ctr_0</span> <span class=\"o\">+</span> <span class=\"n\">fd_src</span> <span class=\"o\">+</span> <span class=\"mi\">60</span><span class=\"p\">;</span>\n",
       "      <span class=\"kt\">double</span> <span class=\"o\">*</span> <span class=\"n\">RESTRICT</span> <span class=\"k\">const</span> <span class=\"n\">fd_src_W</span> <span class=\"o\">=</span> <span class=\"mi\">30</span><span class=\"o\">*</span><span class=\"n\">ctr_0</span> <span class=\"o\">+</span> <span class=\"n\">fd_src</span> <span class=\"o\">-</span> <span class=\"mi\">30</span><span class=\"p\">;</span>\n",
       "      <span class=\"k\">for</span> <span class=\"p\">(</span><span class=\"kt\">int</span> <span class=\"n\">ctr_1</span> <span class=\"o\">=</span> <span class=\"mi\">2</span><span class=\"p\">;</span> <span class=\"n\">ctr_1</span> <span class=\"o\">&lt;</span> <span class=\"mi\">28</span><span class=\"p\">;</span> <span class=\"n\">ctr_1</span> <span class=\"o\">+=</span> <span class=\"mi\">1</span><span class=\"p\">)</span>\n",
       "      <span class=\"p\">{</span>\n",
       "         <span class=\"n\">fd_dst_C</span><span class=\"p\">[</span><span class=\"n\">ctr_1</span><span class=\"p\">]</span> <span class=\"o\">=</span> <span class=\"n\">fd_src_2E</span><span class=\"p\">[</span><span class=\"n\">ctr_1</span><span class=\"p\">]</span> <span class=\"o\">+</span> <span class=\"n\">fd_src_W</span><span class=\"p\">[</span><span class=\"n\">ctr_1</span><span class=\"p\">];</span>\n",
       "      <span class=\"p\">}</span>\n",
       "   <span class=\"p\">}</span>\n",
       "<span class=\"p\">}</span>\n",
       "</pre></div>\n"
      ],
      "text/plain": [
       "FUNC_PREFIX void kernel(double * RESTRICT fd_dst, double * RESTRICT const fd_src)\n",
       "{\n",
       "   for (int ctr_0 = 2; ctr_0 < 18; ctr_0 += 1)\n",
       "   {\n",
       "      double * RESTRICT fd_dst_C = 30*ctr_0 + fd_dst;\n",
       "      double * RESTRICT const fd_src_2E = 30*ctr_0 + fd_src + 60;\n",
       "      double * RESTRICT const fd_src_W = 30*ctr_0 + fd_src - 30;\n",
       "      for (int ctr_1 = 2; ctr_1 < 28; ctr_1 += 1)\n",
       "      {\n",
       "         fd_dst_C[ctr_1] = fd_src_2E[ctr_1] + fd_src_W[ctr_1];\n",
       "      }\n",
       "   }\n",
       "}"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kernel = ps.create_kernel(ps.Assignment(dst[0,0], src[2, 0] + src[-1, 0]))\n",
    "ps.show_code(kernel)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When no additional ghost layer information is given, *pystencils* looks at all neighboring field accesses and introduces the required number of ghost layers **for all directions**. In the example above the largest neighbor accesses was ``src[2, 0]``, so theoretically we would need 2 ghost layers only the the end of the x coordinate. \n",
    "By default *pystencils* introduces 2 ghost layers at all borders of the domain. The next cell shows how to change this behavior. Be careful with manual ghost layer specification, wrong values may lead to SEGFAULTs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>.highlight .hll { background-color: #ffffcc }\n",
       ".highlight  { background: #f8f8f8; }\n",
       ".highlight .c { color: #408080; font-style: italic } /* Comment */\n",
       ".highlight .err { border: 1px solid #FF0000 } /* Error */\n",
       ".highlight .k { color: #008000; font-weight: bold } /* Keyword */\n",
       ".highlight .o { color: #666666 } /* Operator */\n",
       ".highlight .ch { color: #408080; font-style: italic } /* Comment.Hashbang */\n",
       ".highlight .cm { color: #408080; font-style: italic } /* Comment.Multiline */\n",
       ".highlight .cp { color: #BC7A00 } /* Comment.Preproc */\n",
       ".highlight .cpf { color: #408080; font-style: italic } /* Comment.PreprocFile */\n",
       ".highlight .c1 { color: #408080; font-style: italic } /* Comment.Single */\n",
       ".highlight .cs { color: #408080; font-style: italic } /* Comment.Special */\n",
       ".highlight .gd { color: #A00000 } /* Generic.Deleted */\n",
       ".highlight .ge { font-style: italic } /* Generic.Emph */\n",
       ".highlight .gr { color: #FF0000 } /* Generic.Error */\n",
       ".highlight .gh { color: #000080; font-weight: bold } /* Generic.Heading */\n",
       ".highlight .gi { color: #00A000 } /* Generic.Inserted */\n",
       ".highlight .go { color: #888888 } /* Generic.Output */\n",
       ".highlight .gp { color: #000080; font-weight: bold } /* Generic.Prompt */\n",
       ".highlight .gs { font-weight: bold } /* Generic.Strong */\n",
       ".highlight .gu { color: #800080; font-weight: bold } /* Generic.Subheading */\n",
       ".highlight .gt { color: #0044DD } /* Generic.Traceback */\n",
       ".highlight .kc { color: #008000; font-weight: bold } /* Keyword.Constant */\n",
       ".highlight .kd { color: #008000; font-weight: bold } /* Keyword.Declaration */\n",
       ".highlight .kn { color: #008000; font-weight: bold } /* Keyword.Namespace */\n",
       ".highlight .kp { color: #008000 } /* Keyword.Pseudo */\n",
       ".highlight .kr { color: #008000; font-weight: bold } /* Keyword.Reserved */\n",
       ".highlight .kt { color: #B00040 } /* Keyword.Type */\n",
       ".highlight .m { color: #666666 } /* Literal.Number */\n",
       ".highlight .s { color: #BA2121 } /* Literal.String */\n",
       ".highlight .na { color: #7D9029 } /* Name.Attribute */\n",
       ".highlight .nb { color: #008000 } /* Name.Builtin */\n",
       ".highlight .nc { color: #0000FF; font-weight: bold } /* Name.Class */\n",
       ".highlight .no { color: #880000 } /* Name.Constant */\n",
       ".highlight .nd { color: #AA22FF } /* Name.Decorator */\n",
       ".highlight .ni { color: #999999; font-weight: bold } /* Name.Entity */\n",
       ".highlight .ne { color: #D2413A; font-weight: bold } /* Name.Exception */\n",
       ".highlight .nf { color: #0000FF } /* Name.Function */\n",
       ".highlight .nl { color: #A0A000 } /* Name.Label */\n",
       ".highlight .nn { color: #0000FF; font-weight: bold } /* Name.Namespace */\n",
       ".highlight .nt { color: #008000; font-weight: bold } /* Name.Tag */\n",
       ".highlight .nv { color: #19177C } /* Name.Variable */\n",
       ".highlight .ow { color: #AA22FF; font-weight: bold } /* Operator.Word */\n",
       ".highlight .w { color: #bbbbbb } /* Text.Whitespace */\n",
       ".highlight .mb { color: #666666 } /* Literal.Number.Bin */\n",
       ".highlight .mf { color: #666666 } /* Literal.Number.Float */\n",
       ".highlight .mh { color: #666666 } /* Literal.Number.Hex */\n",
       ".highlight .mi { color: #666666 } /* Literal.Number.Integer */\n",
       ".highlight .mo { color: #666666 } /* Literal.Number.Oct */\n",
       ".highlight .sa { color: #BA2121 } /* Literal.String.Affix */\n",
       ".highlight .sb { color: #BA2121 } /* Literal.String.Backtick */\n",
       ".highlight .sc { color: #BA2121 } /* Literal.String.Char */\n",
       ".highlight .dl { color: #BA2121 } /* Literal.String.Delimiter */\n",
       ".highlight .sd { color: #BA2121; font-style: italic } /* Literal.String.Doc */\n",
       ".highlight .s2 { color: #BA2121 } /* Literal.String.Double */\n",
       ".highlight .se { color: #BB6622; font-weight: bold } /* Literal.String.Escape */\n",
       ".highlight .sh { color: #BA2121 } /* Literal.String.Heredoc */\n",
       ".highlight .si { color: #BB6688; font-weight: bold } /* Literal.String.Interpol */\n",
       ".highlight .sx { color: #008000 } /* Literal.String.Other */\n",
       ".highlight .sr { color: #BB6688 } /* Literal.String.Regex */\n",
       ".highlight .s1 { color: #BA2121 } /* Literal.String.Single */\n",
       ".highlight .ss { color: #19177C } /* Literal.String.Symbol */\n",
       ".highlight .bp { color: #008000 } /* Name.Builtin.Pseudo */\n",
       ".highlight .fm { color: #0000FF } /* Name.Function.Magic */\n",
       ".highlight .vc { color: #19177C } /* Name.Variable.Class */\n",
       ".highlight .vg { color: #19177C } /* Name.Variable.Global */\n",
       ".highlight .vi { color: #19177C } /* Name.Variable.Instance */\n",
       ".highlight .vm { color: #19177C } /* Name.Variable.Magic */\n",
       ".highlight .il { color: #666666 } /* Literal.Number.Integer.Long */</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div class=\"highlight\"><pre><span></span><span class=\"n\">FUNC_PREFIX</span> <span class=\"kt\">void</span> <span class=\"nf\">kernel</span><span class=\"p\">(</span><span class=\"kt\">double</span> <span class=\"o\">*</span> <span class=\"n\">RESTRICT</span> <span class=\"n\">fd_dst</span><span class=\"p\">,</span> <span class=\"kt\">double</span> <span class=\"o\">*</span> <span class=\"n\">RESTRICT</span> <span class=\"k\">const</span> <span class=\"n\">fd_src</span><span class=\"p\">)</span>\n",
       "<span class=\"p\">{</span>\n",
       "   <span class=\"k\">for</span> <span class=\"p\">(</span><span class=\"kt\">int</span> <span class=\"n\">ctr_0</span> <span class=\"o\">=</span> <span class=\"mi\">0</span><span class=\"p\">;</span> <span class=\"n\">ctr_0</span> <span class=\"o\">&lt;</span> <span class=\"mi\">18</span><span class=\"p\">;</span> <span class=\"n\">ctr_0</span> <span class=\"o\">+=</span> <span class=\"mi\">1</span><span class=\"p\">)</span>\n",
       "   <span class=\"p\">{</span>\n",
       "      <span class=\"kt\">double</span> <span class=\"o\">*</span> <span class=\"n\">RESTRICT</span> <span class=\"n\">fd_dst_C</span> <span class=\"o\">=</span> <span class=\"mi\">30</span><span class=\"o\">*</span><span class=\"n\">ctr_0</span> <span class=\"o\">+</span> <span class=\"n\">fd_dst</span><span class=\"p\">;</span>\n",
       "      <span class=\"kt\">double</span> <span class=\"o\">*</span> <span class=\"n\">RESTRICT</span> <span class=\"k\">const</span> <span class=\"n\">fd_src_2E</span> <span class=\"o\">=</span> <span class=\"mi\">30</span><span class=\"o\">*</span><span class=\"n\">ctr_0</span> <span class=\"o\">+</span> <span class=\"n\">fd_src</span> <span class=\"o\">+</span> <span class=\"mi\">60</span><span class=\"p\">;</span>\n",
       "      <span class=\"kt\">double</span> <span class=\"o\">*</span> <span class=\"n\">RESTRICT</span> <span class=\"k\">const</span> <span class=\"n\">fd_src_W</span> <span class=\"o\">=</span> <span class=\"mi\">30</span><span class=\"o\">*</span><span class=\"n\">ctr_0</span> <span class=\"o\">+</span> <span class=\"n\">fd_src</span> <span class=\"o\">-</span> <span class=\"mi\">30</span><span class=\"p\">;</span>\n",
       "      <span class=\"k\">for</span> <span class=\"p\">(</span><span class=\"kt\">int</span> <span class=\"n\">ctr_1</span> <span class=\"o\">=</span> <span class=\"mi\">1</span><span class=\"p\">;</span> <span class=\"n\">ctr_1</span> <span class=\"o\">&lt;</span> <span class=\"mi\">30</span><span class=\"p\">;</span> <span class=\"n\">ctr_1</span> <span class=\"o\">+=</span> <span class=\"mi\">1</span><span class=\"p\">)</span>\n",
       "      <span class=\"p\">{</span>\n",
       "         <span class=\"n\">fd_dst_C</span><span class=\"p\">[</span><span class=\"n\">ctr_1</span><span class=\"p\">]</span> <span class=\"o\">=</span> <span class=\"n\">fd_src_2E</span><span class=\"p\">[</span><span class=\"n\">ctr_1</span><span class=\"p\">]</span> <span class=\"o\">+</span> <span class=\"n\">fd_src_W</span><span class=\"p\">[</span><span class=\"n\">ctr_1</span><span class=\"p\">];</span>\n",
       "      <span class=\"p\">}</span>\n",
       "   <span class=\"p\">}</span>\n",
       "<span class=\"p\">}</span>\n",
       "</pre></div>\n"
      ],
      "text/plain": [
       "FUNC_PREFIX void kernel(double * RESTRICT fd_dst, double * RESTRICT const fd_src)\n",
       "{\n",
       "   for (int ctr_0 = 0; ctr_0 < 18; ctr_0 += 1)\n",
       "   {\n",
       "      double * RESTRICT fd_dst_C = 30*ctr_0 + fd_dst;\n",
       "      double * RESTRICT const fd_src_2E = 30*ctr_0 + fd_src + 60;\n",
       "      double * RESTRICT const fd_src_W = 30*ctr_0 + fd_src - 30;\n",
       "      for (int ctr_1 = 1; ctr_1 < 30; ctr_1 += 1)\n",
       "      {\n",
       "         fd_dst_C[ctr_1] = fd_src_2E[ctr_1] + fd_src_W[ctr_1];\n",
       "      }\n",
       "   }\n",
       "}"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gl_spec = [(0, 2),   # 0 ghost layers at the left, 2 at the right border\n",
    "           (1, 0)]   # 1 ghost layer at the lower y, one at the upper y coordinate\n",
    "kernel = ps.create_kernel(ps.Assignment(dst[0,0], src[2, 0] + src[-1, 0]), ghost_layers=gl_spec)\n",
    "ps.show_code(kernel)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2 ) Restrictions\n",
    "\n",
    "\n",
    "### a) Independence Restriction\n",
    "\n",
    "*pystencils* only works for kernels where each array element can be updated independently from all other elements. This restriction ensures that the kernels can be easily parallelized and also be run on the GPU. Trying to define kernels where the results depends on the iteration order, leads to a ValueError."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Field dst is written at two different locations\n"
     ]
    }
   ],
   "source": [
    "invalid_description = [\n",
    "    ps.Assignment(dst[1, 0], src[1, 0] + src[-1, 0]),\n",
    "    ps.Assignment(dst[0, 0], src[1, 0] - src[-1, 0]),\n",
    "]\n",
    "try:\n",
    "    invalid_kernel = ps.create_kernel(invalid_description)\n",
    "    assert False, \"Should never be executed\"\n",
    "except ValueError as e:\n",
    "    print(e)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The independence restriction makes sure that the kernel can be safely parallelized by checking the following conditions: If a field is modified inside the kernel, it may only be modified at a single spatial position. In that case the field may also only be read at this position. Fields that are not modified may be read at multiple neighboring positions.\n",
    "\n",
    "Specifically, this rule allows for in-place updates that don't access neighbors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "valid_kernel = ps.create_kernel(ps.Assignment(src[0,0], 2*src[0,0] + 42))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If a field stores multiple values per cell, as in the next example, this restriction only applies for accesses with the same index."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "v = ps.fields(\"v(2): double[2D]\")\n",
    "valid_kernel = ps.create_kernel([ps.Assignment(v[0,0](1), 2*v[0,0](1) + 42),\n",
    "                                 ps.Assignment(v[0,1](0), 2*v[1,0](0) + 42)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### b) Static Single Assignment Form\n",
    "\n",
    "All assignments that don't write to a field must be in SSA form\n",
    "1. Each sympy symbol may only occur once as a left-hand-side (fields can be written multiple times)\n",
    "2. A symbol has to be defined before it is used. If it is never defined it is introduced as function parameter\n",
    "\n",
    "The next cell demonstrates the first SSA restriction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Assignments not in SSA form, multiple assignments to a\n"
     ]
    }
   ],
   "source": [
    "@ps.kernel\n",
    "def not_allowed():\n",
    "    a, b = sp.symbols(\"a b\")\n",
    "    a @= src[0, 0]\n",
    "    b @= a + 3\n",
    "    a @= src[-1, 0]\n",
    "    dst[0, 0] @= a + b\n",
    "try:\n",
    "    ps.create_kernel(not_allowed)\n",
    "    assert False\n",
    "except ValueError as e:\n",
    "    print(e)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, for right hand sides that are Field.Accesses this is allowed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "KernelFunction kernel([<double * RESTRICT fd_dst>, <double * RESTRICT const fd_src>])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@ps.kernel\n",
    "def allowed():\n",
    "    dst[0, 0] @= src[0, 1] + src[1, 0]\n",
    "    dst[0, 0] @= 2 * dst[0, 0]\n",
    "ps.create_kernel(allowed)"
   ]
  }
 ],
 "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.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}