01_tutorial_getting_started.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pystencils.session import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 01: Getting Started\n",
    "\n",
    "\n",
    "## Overview\n",
    "\n",
    "\n",
    "*pystencils* is a package that can speed up transformations on *numpy* arrays. All computations are carried out fully parallel on CPUs (single node with OpenMP, multiple nodes with MPI) or on GPUs.\n",
    "\n",
    "It is suited for applications that run the same operation on *numpy* arrays multiple times like for example:\n",
    "- numerical simulations (finite difference, lattice Boltzmann)\n",
    "- image processing\n",
    "\n",
    "As the name suggests, *pystencils* was developed for **stencil codes**, i.e. operations that update an element of a numpy array using only a local neighborhood. \n",
    "\n",
    "<img width=\"25%\" src='../../img/pystencils_stencil_four_points_with_arrows.svg'>\n",
    "\n",
    "It generates C code, compiles it behind the scenes, and lets you call the compiled C function as if it was a native Python function. \n",
    "\n",
    "But lets not dive too deep into the concepts of *pystencils* here, they are covered in detail in the following tutorials. Lets instead look at a simple example, that computes the average neighbor values of a *numpy* array.\n",
    "\n",
    "We create two rather large arrays for in- and output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_arr = np.random.rand(1024, 1024)\n",
    "output_arr = np.zeros_like(input_arr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first implement a version of this algorithm using pure numpy and benchmark it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def numpy_kernel():\n",
    "    output_arr[1:-1, 1:-1] = input_arr[2:, 1:-1] + input_arr[:-2, 1:-1] + \\\n",
    "                             input_arr[1:-1, 2:] + input_arr[1:-1, :-2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7.13 ms ± 254 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit \n",
    "numpy_kernel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets see how to run the same algorithm with *pystencils*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm8AAAAZCAYAAABn0sl9AAAABHNCSVQICAgIfAhkiAAADHBJREFUeJztnXmQHUUdxz+BNQkSBAUJoiSLQEANFAQEQRI2oJyGq4QgCiwUoiIGVApRRBeUIEQlBlSUEiKHBQbkUBDQWGswSAKRI3IYLokkEgjIFSFAXP/49tSbndfzXs/xjt39fapevX3d0z395vubfr/59bFgGIZhGIZhDFrGAH3Ab1rdkDakE12b2SXVNw14CHjN1XtKSfUahmEYhjGAWSvj8RPc+99ynu8ryBE5Mmf5VvIwcF+TznUE8CPgdWAmcBZwF+U7iIYRZzSwBphV45ijkA32Acc3o1FG06inv2nfnphuA5NC/W1HxpNFztuijOXKKt9Krge+DmwOPNngc30i9r48lt7Z4PMaQ5uD0APd9Sn5mwEXAq8Co5rVKKNp1NLftG9fTLeBSaH+NmvkbUf3njfytqNryKM5y7eS6AIf0oRzberel9c8yjDK5RDgeWCeJ28YcJnLv7iEc3Wjp8muEuoyyiFNf9O+vTHdBiaF+luf89YBnAw8gOZbPQWc5iqbACwDVsSOn4jmwD2OhvmeBRYC02PHnIdE3wZ5kGuohAI/U/crtgd3A0+T33nbBbgWeAZ4A/gX8DMqjhpAD7omk93nvtirh0rE75hEXnesjm7gOuAJpN/LwHz817mTyjDsOOAapN//0A1aLz/icGSAL7lzLkZRyhGJ841y331+In0dZDt9KEwc50SXfpyn/T5C7BF0jfuAGcDOwI3ACy7tQ7HjutD3fhpYDfwbuA04OLA9A4X1gT2B36L7M8k0l38ssKqJ7UrDdC6XWvq3m/YQpv9Q0N50G3y6QYB2yWHT4cDvgI+j+V0XARuiOVdbAZu4k0V8AzgHWIou2ko0jrsTsK/LBw2T/hI5HXcCf4jV8ed637KNuAE5E6Pp78DW41jgEmRgNyHHbSs0hj0F+Ai6hr3u+G5gLLruEb3ABsixvt+1JSI+F++naKHDPGTMGwL7A1cAWwNnetq3BbAAWAJchZyplwPzpyNHbSXwKxRZ3c+l74Ns6U137KvoRt0FWA94xaV/lIqjt5dra8Se7n2up91JQu0RKkP449GNcgtypscAj7i8meh6rwRuRtdzLLA3sBv9NRjoHIDuf99ipA8A30PzMOdR0aRVmM7lk6Z/u2kP4foPBe1Nt8GlG+TU7hLk3Z6JIm0Rk+gfAQJd9LeAO1wjkmyU+HyCK39CSEPalOiJwPcdOvEvJhiHok2PAe9N5O2JvO7kmHevqyv0HHG28KQNR87Pm4k2RPX1UR2xCMnf1eUtRY59RAdy8vvo/0MKcLZLPyCWdi6ypT8hxzZiLXRzPu45d5Ks9niVa8cryHlOMt3lXwusm8gbhTqKwcQc5FyPTKR3APcA/0BOO1QixEUmP3eTbwjGdG4MPv3bTXvIpv9Q0N50689A1g1yarezOyAtEvaQy5/iPu/hPv8isLEXu+N38uS9E0WyfI5HI7kWrYANZW3kTPzek9eJ37G6gGpnJc71yKjXi6X1kt95S+NQV/ZoT33PUD3EGZIfOfs+Z3YcckyfSKRHdvPDWNpCFNn7ossb59InuM8/934jf72h9vgw6cOxO7i2342/o0nSKvvNQi1bH4k6yjmevLPRtdg1ltZD634ITOd85NG/3bSHbPoPBu1Nt/4MZt0gp3ZXuAP2SMm/w+VHkZuNgBdd2k3AVHSh0liIIlA+J2AGGlZNciKa5/U6GnqdWKN+H/XKb4fG0tfPUOdlaPgzWaYTv2N1l0s/H4mQfM13+TvGyvSS33kbA/wYhZf/S/+5cX1oiDNZ320pddXLX+Tyt0zJf8rlbxBLG+7adb/7vD5yXqejcHEf8AWXd6r7fERK/XGy2OO66AZZgX/eZ3Qv7BVwXki331AmoUjlcnfeT+aoo4itH4h/C5+dkTbnJ9J7yPZD8E+q7bDWa3aNutpN56J9VBmE2E9W/dtRewjXv1n3eBH9TbdqmqHboOtvV6Co0rCU/GUoAhNnPPIeV7nK30JRqQmJ4zrQRHbfPmlvR6Imv/xUNMz3WfSjPguFGcekfYGc5RehiE8oaRe+E78RP0qY8ced5l7yOW/vRxqucXXMAr6DxJ/tyvZ46ktzPOrlP+byk6HriMhxHZtIvx0tetgYLZeO34zL0FMLyJai40IItcfdXP6lKfU8i268kNXYafabhf2A71KJjmbtTIrauu+BpAOF7h+i+oGrh2w/BKdQ/dByAxVbTubVm3DcLjoX7aPKItR+QvVvZ+0hTP9m3ONF9TfdWqPboOpvR7rMtC1Aogt6S0r+cOBjwK/dcSsTDdiOdEEOQ2IkncYFaFguzqNoflQIoeW/TfXqx1qkhTw78TtW97j0d2Q4Ry/5nLeLXH63J+9TpDtvafXVy48ib2khaV/kDeB0KhG1WejpJRr7vxItjx6BbojFKXXXop49nuTSP+cpG90L9waeK81+85KnMyli62lTATYg/Kl7Zsb2QjnbDrRa56J9VCOoZT+h+g8E7aG2/s24x8vU33QTze6bB2x/G602XeNeaRGOaNVjmnP3BvBH97oD2B1NUFzq8rd3776LvjsVJyBiOBpG/H7i2NuRI1mPLOUXAGegyYGvBdT9OnAr8t5Hus+1uMu1ZSJaFVOEaEnx2in50fDldZ68PTxpRbkXPTV1Ub2oYEvgfSi0/GIiL1o5uhca259P5TrOBT6Nhk7XJWyVaZJ69hg96d3jKRvd6KHRPp/9NpOitj4JrUhOLppZTfr8lAlo7slf0NPiX/M0vARaqXPRPqoVhOo/ELSH2vo3+h5vpv6mm2iHvrkt+9tokuCURPrXqHh80R5nO+CPtmwJ/AdFXOJhzS+78sd4ytxI9bDcpu74SYn0b7nG1yNL+SgqmGVC45FUX6tO/FGqbZCxLqEyET/OcKrDur34DW4UGkZMW1QSLQpJargPClGXHXmLIrJPAu+Opa9NJcR+hqfcWshOnnXHxFekRv8/d4V7PzDl3HGy2uN96GZJm/C62J37cE/eOPo7zz77LULWJ8Gitn4heigYneGcPWQbgvHRTfan+HbSuWgf1Shq2U8Z+vfQGu0hm/6NvsfL1t90E83umwdsfxvf5+1c9EWvA65G89u6gG3R9g2bUYm8TUOO2ELgQfQjvDmVH9rjkJMREf07rHPQuPcqV24OtaNXSQdmmCetFiHlI494HcK5GTlkh9J/3zsfj6DrcSn6zrciR+5tyFGZCDyHnLx6vIq8+YloOfUSZAg3oU2Vf4L2lJuDdFyGrve+KEQ9NfD7hXInmlx5GvB3NFdtFYpKjkdPCjM85SIH9CD3OR5dW4qieFug7xayD2AWexwBfBBdrzdS6jsdXdOrXb0PorD29ug+eE/s2DT77UFh81pMprK3X1Hy2vrBSMcs+xa2inbUuWgfFdFD4+1lqOjfLO2hPP1rYbpVKEu3orS8v407b5ej1R8no/lRL6Af38+jEN/zyEMGebQdaHXEYegiLUebtJ5H9b+/mgd8KfYagVYXzkHjv8lVJ9Gk+00S6RsT9qWzlH+Xe38uoN6Il9CeZFOQp+/bITnOlWh15VdRB7w3cnKWI4fnmgznPgptP7Iv0mkY2mH6AfeajCZj7o80uh85mS9SvvMGiszei+YqHI2c0seBbwI/IP0mnIuct5epDpHPRc7bInSt65HFHse7Ntb6/7o3oweX09Gw7t7IphajjRPj+OwXNP/w6jrtXlonP4Qitv5hNLR9QQntaAbtpHPRPipJM+xlqOjfjHu8bP1rYbpVKKNvLsJQ6m/rciqK2iRZQPX+XkvItmAhpPzxKEKVlWjT4cn1DjQGNWn2m5e8E2jz2Hq02eXmGc83FPHpXLSPagS17Mf0z0ejfqPimG7lE9I3W39bgG2RJ5vcGX0qitocj5bjzkTDhmMD6w0tfznhG37GGe3afWGOssbgIc1+szAKhf23Rzf3qe7vLNvi5LH1h/Fv32NU49O5aB9VFqH2Y/rno1G/UaZbY0nTzfrbErkTDbslORFtFrgahVKTkwS70cXvTKm3Xvnof3T6/gWHYYSSZr+hdOFfEj47dkw3Zuutxqdz0T6qDLqobz+mfzEa8RvVhenWaHy6dWH9bWnsg8KOaVtgpHEWmrDYUe/AFE5CS3wNowh57TcLZuutJ4/ORXUrC9O/GPYbNTAx3ZrANLIPN9xNsTlnJwBbFyhvGBF57DcLZuvtQVadi+pWFqZ/cew3amBiuhmGYRiGYRiGYRiGYRiGYRiGYRiGYRiGYRiGYRhGC/k/Z0xkHfv+t+EAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$${{dst}_{(0,0)}} \\leftarrow \\frac{{{src}_{(-1,0)}}}{4} + \\frac{{{src}_{(0,-1)}}}{4} + \\frac{{{src}_{(0,1)}}}{4} + \\frac{{{src}_{(1,0)}}}{4}$$"
      ],
      "text/plain": [
       "         src_W   src_S   src_N   src_E\n",
       "dst_C := ───── + ───── + ───── + ─────\n",
       "           4       4       4       4  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "src, dst = ps.fields(src=input_arr, dst=output_arr)\n",
    "\n",
    "symbolic_description = ps.Assignment(dst[0,0], \n",
    "                                     (src[1, 0] + src[-1, 0] + src[0, 1] + src[0, -1]) / 4)\n",
    "symbolic_description"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANoAAADTCAYAAADnEg0TAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAClpJREFUeJzt3W9sVWcdwPHv0z+X0tLC4I7hUGZchluiiTFxMdmcm7qRS7Y5/JNg3AvFxAjTSDRzJMZk0akkRoPGzMRIfLGgkTmXbHEdL/ybEaPGF1OjYwwWUNDJZUBLKbDePr5o6coKtL1rf4eefj9v6D33nPVHeb65596de5tyzkiaXS1FDyDNB4YmBWgregDNjNRbvw24DVgKpIvsloGjwK9zrfq7qNkEyedoc1/qrX8R+PQ0D3s416rfnY15NJGhzXGpt94N/AFon+ahZ4B351r11MxPpdfyOdrcdwPTjwxgAXD9DM+iizC0uW9BQcdqGgytjPqPtXDfratYt+o6Xni2UvQ4MrRy6uga5ms/O8SNd/QXPYpGGFoZtVdg6VWNosfQqwxNCmBoUgBDkwJ4CVZZbVm3kgPPdXB4f4U19x7nzg19RY80nxlaWW19/FDRI+hVnjpKAQxNCuCpYxmtvXL1Re976sjzgZNolKGV0bmYdu3oZvuDy9m5d1/BE817njqW1XADdj/ZzbIVQ0WPIkMrr107erjprn7Sxd5srUiGVkaNIXjmiW5uX+9FxZcJQyujpx/p4ea7+2lpLXoSjTK0Mjq4p8KvH+3h/ruu478HF7Bt8/KiR5rvfNWxjDZurTM40MnL/72ahz6R+Mw3PIUsmI9oZXXy+BXk3MKXfzzytQplaGXUGGrl7OnOsdtnBrtoDPlvXSB/+HPfxM8LPHli8YRtA30Tt8HwLMyjCzC0ue/YebdyhoG+K8g5jduWOHniCiZ+hufLsz+ewNDK4B/Af8ZunT7VyXBj4uv6w402zgwuHLfl37lW3Tv74wkMbc7LtWoG7gNG3n/W3n6Wjs6TdHSeHNvp3O229ldGt/wL+Gz0rPOZHwleIqm3vhpYxrlfcvHFtbshNfj2L28Z3WUYOOojWTxDK7GUUgaGcs7NfGS4ZpCnjlIAQ5MCGJoUwNCkAIYmBTA0KYChSQEMTQpgaFIAQ5MCGJoUwNCkAIYmBTA0KYChSQEMTQpgaFIAQ5MCGJoUwNCkAIYmBTA0KYChSQEMTQpgaFIAQ5MCGJoUwNCkAIYmBTA0KYChSQEMTQpgaFIAQ5MCGJoUwNCkAIYmBTA0KYChSQEMTQpgaFIAQ5MCGJoUwNCkAIYmBTA0KYChSQEMTQpgaFIAQ5MCGJoUwNCkAIYmBTA0KYChSQEMTQrQFvnNUm+9A3gL0H6J3V4BXsy16mDMVNLkUm+9E3gzl167Z4H9uVY9M+H4nPMsjfaab9Rb/zzwSWDhFHYfBH4CfCvXqjEDllBKKQNDOedLLQ5dQuqtJ2ALsB7omMIhp4DtuVb9/viNIaeOqbdeAzYxtcgY3e9TwD2zNpQ0NR8FPsHUIgPoBD6XeusfGL8x6jna7U0ed8eMTiFNX7Nr97zjokK7ssnjqjM6hTR9za7B846LCi1N2NJ/rIX7bl3FulXX8cKzlYsc56uiKtrENdjE2i1uIXd0DfO1nx3ixjv6C5tBakYTa7e40NorsPSqRmHfX2pWE2vXUzMpgKFJAQxNChB6CdYEW9at5MBzHRzeX2HNvce5c0PfVA5LKS3OOZ+Y7fE0f0x7TU1z7RYb2tbHD01115RSG1ADvgC8N6V0T875iVmbTfNGSunDwM6U0m+B7wBP55wv/WLHNNYuzIFTx5TSm1JKDwEvATuAWxm5nmxxkXOpVBYDp4H3AT8FXkopfTWltHKmvkGxj2iXMjjQyZ6/vAfYO7plwbh7E7AhpbQqfrC5J6X05aJnuMytGfd19+ifXwLu56+7T7L6HYN0dL2ud5MUF9raK1df9L6njjzPwIklDPQtAjITryzpYOSR7dbZGq9E2oCHih5iDhh+ze0FwDADfVVO9vWfF9pka/cCigvt3EC7dnSz/cHl7Ny777z7l73hMDfceIiRU8Z7GPlBdI7eewrYlHN+JG7guce3yUxNSmkD8D2ga3TTANAKPMYN73oTi5etOO+AydbuBRT7HG24Abuf7GbZiqEJ96UES5cfyzl/DFgBPAC8yMgPYcGE/aXXp8LI2toH3A9clXO+lyXV46SJl+pecu1eQLGh7drRw0139V/wLzJOzvlEzvn7wLXA+4HtwDMBE2p++B3wI+A24Lqc8w9yzpf+X01TXLvnFBdaYwieeaKb29dP+cLMPOKPOeeNOecXZ3M8zR85530550055z/nqXzkQBNrt7jQnn6kh5vv7qeltbARpKY0sXaLC+3gngq/ebSHBz74Rl462M62zcsLm0WajibWbnGvOm7cWh/7etMt17B52/8Km0WajibW7uVxZcjDvz9Q9AhSU6a4di+P0KSSMzQpQFRoZ4OPk2bKjKzdqND+3uRxf5vRKaTpa3btnndcVGg7gaPTPOY4Ix8LLhVpBzClNySPcwR4bPyGyM/eXwl8BHgbk/+Si38Aj+Va9WDEbGXlRcUzI/XWr2Fk7V7P5L/k4u/Az3Otevi8/0ZUaIpnaJcPX3WUAhiaFMDQpACGJgUwNCmAoUkBDE0KYGhSAEOTAhiaFMDQpACGJgUwNCmAoUkBDE0KYGhSAEOTAhiaFMDQpACGJgUwNCmAoUkBDE0KYGhSAEOTAhiaFMDQpACGJgUwNCmAoUkBDE0KYGhSAEOTAhiaFMDQpACGJgUwNCmAoUkBDE0KYGhSAEOTAhiaFMDQpACGJgUwNCmAoUkBDE0KYGhSAEOTAhiaFMDQpACGJgUwNCmAoUkBDE0KYGhSAEOTAqScc9Ez6HVKvfUE3Am8n5yXMtDXBcCfdn0IgBvX/AKArp4BUjoK/Ar4Za5V/ccP0lb0AJoRXwE+DsDgQBfHj6wkpczqd47ce6K+jpwTra3/ZuGiU8Ba4O3ANwuad97x1HGOS731JcD6sQ2VjtOklMk5jW3LOZFSptJxZtyhH0+99e64Sec3Q5v73gq0jt1qa2tQ6Tg1Ya/KwgFa2xrjtrQDq2d9OgGGVgaVCVsWLTlGSsNjt1PKdC85NqVjNSsMrYxeOXOar29o4XPvg0P7ILU0WLBwsOix5jNDK6OFi4bZ8sOjvPO2kUezRYuPkdLkx2nWGFoZtVfg6muPj93u6ukrcBrhy/vl1dbWIKUh2ipnX/MiiApgaGXW2naKniteLnoMeeoohfARray2rFvJgec6OLy/wpp7j3PnBp+nFcjQymrr44eKHkGv8tRRCmBoUgBPHcto7ZUXv4bxqSPPB06iUYZWRudi2rWjm+0PLmfn3n0FTzTveepYVsMN2P1kN8tWDBU9igytvHbt6OGmu/q9xvHyYGhl1BiCZ57o5vb1/UWPohGGVkZPP9LDzXf309I6+b4KYWhldHBPhd882sMDH3wjLx1sZ9vm5UWPNN/5qmMZbdxaH/t60y3XsHnb/wqcRviIVn4P//5A0SPI0KQQhjb3vZ4PQR2efBfNBEOb++qT7zIrx2oaDG2Oy7Xqc8DBJg7dn2tVL80KYmjlsBHYM439/wlsmqVZdAH+kosSSb31q4GlwMWuuxoGjuVa9XDcVAJDk0J46igFMDQpwP8BiJCQfoddbocAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 216x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(3,3))\n",
    "ps.visualize_stencil_expression(symbolic_description.rhs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we first have created a symbolic notation of the stencil itself. This representation is built on top of *sympy* and is explained in detail in the next section. \n",
    "This description is then compiled and loaded as a Python function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "kernel = ps.create_kernel(symbolic_description).compile()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This whole process might seem overly complicated. We have already spent more lines of code as for the *numpy* implementation and don't have anything running yet! However, this multi stage process of formulating the algorithm symbolically, and just in the end actually running it, is what makes *pystencils* faster and more flexible than other approaches.\n",
    "\n",
    "Now finally lets benchmark the *pystencils* kernel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pystencils_kernel():\n",
    "    kernel(src=input_arr, dst=output_arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.36 ms ± 51.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "pystencils_kernel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [