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": [
    "This benchmark shows that *pystencils* is a lot faster than pure *numpy*, especially for large arrays. \n",
    "If you are interested in performance details and comparison to other packages like Cython, have a look at [this page](demo_benchmark.ipynb).\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Short *sympy* introduction\n",
    "\n",
    "In this tutorial we continue with a short *sympy* introduction, since the symbolic kernel definition is built on top of this package. If you already know *sympy* you can skip this section. \n",
    "You can also read the full [sympy documentation here](http://docs.sympy.org/latest/index.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "sp.init_printing()  # enable nice LaTeX output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*sympy* is a package for symbolic calculation. So first we need some symbols:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "sympy.core.symbol.Symbol"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = sp.Symbol(\"x\")\n",
    "y = sp.Symbol(\"y\")\n",
    "type(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The usual mathematical operations are defined for symbols:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKoAAAAZBAMAAABTBqhqAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEHarIkSJZt3NVLsyme8Q6PJIAAAACXBIWXMAAA7EAAAOxAGVKw4bAAACe0lEQVRIDa2Vv2sTYRjHv5deLtdcEo9WBLdrBLfSUgqut4ggqMVBHIRonRyKxUUdlCAuXeSwVDx1iOIgTi26ieg/II2LUhTM6mQq6iBCfH88z937Nj0Ummd4n1/f93N3712eANqWXocU7dnlKC8eS/aM0wAD5ce1HyOiGih/wf05KqqJqm6PiAoYKD9S1GP/z67EBVpCye6SkrjrBcod5bvhuw7mdxQ51SiZVbqqNh4r98/l2+AXcHB3GaFk8wzOSXdDLpk9ySIOqpGO7s3FgEcJbB2hhLD24PmG1N+WS2a2WpaZOisTpydXYZaOUaLeGAz6wrltseRmqVXZomKNpJZOoyYOHT0cUbskTtdpzp9e1DmpnekIJ0NVYuqRK59Fzq+FqTnKOV+55SdEDZaBA7gWrVhUr9rDjK4wdQOnFoDrtI2vnqO80NuuxtT2O8AFXA6f6pzUx+sJXukKU4GgDTzTRT5XA+VgfJ2a4njFDYR4qPJLaXonTVdFHDYifBfeS9P7X9N0VrXrYm5syijXWah9kZLJRVKhCNLzPaAV1v6oAn8DYsB5Ym4oqqGDgdp8q3eIVZ4AmJBTt1DqiYYwOoFGgrqYG2d1Mbs6MlQ1nEEppLZ8W49LfUzpnN/tCwT6sZlaFzeQDL0tA9XqvEE2Urwuyr+Dvivg0pj6CS2q0L26XdyMgQ9KlekM1GRz4iNtAcptOHNT01dtNSYXvwiGNP4GTlx8L7JHqpZRLRS1lCv6xW6RiKkqddtU5WeidMjZ02W/7o8lDv9ROLGxJZsupDNadhiYu7jlLwcRx6YvmISmRMe7Tu1yU/7oh+3lcKmgkn0OBX2jXPgPA/wF1cGLzQch/MoAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$x^{2} \\left(x + y + 5\\right) + x^{2}$$"
      ],
      "text/plain": [
       " 2                2\n",
       "x ⋅(x + y + 5) + x "
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr = x**2 * ( y + x + 5) + x**2\n",
    "expr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can do all sorts of operations on these expressions: expand them, factor them, substitute variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAIYAAAAYBAMAAADNO5iRAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEHarIkSJZt3NVLsyme8Q6PJIAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAB6ElEQVQ4EY2UMUvDQBTH/2lrEtOkBnVxqUHBTXRycsgigoM4OQktFQQHoR9B/ABSdNCIQ3F0UhylmG+gICKiYhzdrFgEEfTs3eVeyglmuHvv/373z7u7EIA/tWtfRJlpo6WVCaMIr+6dkIIMzTDfkLF+JoRXtV41kB267xqZSBlC24e9VOiQBZowQ5RjDQE4ba1MREXkdohMQjsgiTYkhHYvwIZ2HdYmYlkgROFTinS2rmiWxiNLnrwvRZTjwgcnnCAlWbCMFZoeiaSF4okIFWGHuWcuZjzcg+NTjUefunBCuLXVUONR+v7OfDaij2KSGnNicHxuIkg1yD6UbEwGWPS7hPAYOKzdAISoWlu2PCAGCg9DyaaTYIq/RHiM7qHUpIRvtp2QI7+j8DCVvMDu4JwD0qMN75kQBvrlEQNmFO2/RNE0QGS/FOCNeaxH0XYU7bJoYBNWhxBMCfhL+CjPg8gV3/3iRdEH27vFPgWy8DLWeij5Drkk45FPWB9ASjj+FHL81Luc6IPKZyiy3f0+og+XnUdCiErzAvOc6I7Cg8r3qNQ5ITwwg3KTEENjg7eCoB5UHqo9hVkP83EWlODVdEzPNFVYcCcS2QetaWIj7BXzDUP+zIZ7a//N7Xox+C/7F9c39vBXies/0sx5RKtmUccAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$x^{3} + x^{2} y + 6 x^{2}$$"
      ],
      "text/plain": [
       " 3    2        2\n",
       "x  + x ⋅y + 6⋅x "
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr.expand()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAHsAAAAZBAMAAADj1UwdAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEHarIkSJZt3NVLsyme8Q6PJIAAAACXBIWXMAAA7EAAAOxAGVKw4bAAACE0lEQVQ4EaVUPYgTQRT+Jj+bzW4uDndgKXsR7MIdYi3bHIrFESyshMTUBlMI2sgFsbHbRjDanGB11Yml6F117QWukKDotVrlRMETIb75ebsZXQUvDzLz3ve+b+btm8wAxnqvpfVOMnlxMTmJzmr8uPZ1Hnmr9G0OORAczSX3Iy2/9B+L3J7h9rRf2p6B/uV63ZsoximjMtJuNUPSVJ4jHmEMkR3WNVxXtHsON4ickAJGwhFeAV1O155svVD+Awb0zOQMZKQdK+wMJ+rT6YT80oCGzJj8J/JGQ9UWsHh27VxkCQXqgGhcuGqrsnLRjLAuNYUXPN5rktQbQXQq9/3EysM+cBp3o4cmtmQvOMSKg4jjPj4DlQSe9I6C2Mr9TeAGbsnnDvnKQqIapYzrmUp8kAgGEKhumxSNdapI4qmOveHw8afhcJWQeoQvhGUIvgMbLdQmhJ6KNF0NSg5NVTPvhbas/dRAitDXbfSNfH/XpGhUxYOpmXyMwiElyLh1F/Xu5QECuYKCNEmo1j0rTLBsYia/RLjqIvv07bp17c0dpPeETqL8I5yUaBVlLH+H9m9IsSOo8wsdLDUW39ocQNWI88vNO0adype6H2MD8YLoHewCVfo59rc/7diyUrmOLztaCtwrI2KdLyaCnyKLWFl6ZWyM0Ag4NLPfDyMXMdHMhbXp3Oei3Hifp87bKz2EXIEL0mP1C2cKaKjr4RNsAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$x^{2} \\left(x + y + 6\\right)$$"
      ],
      "text/plain": [
       " 2            \n",
       "x ⋅(x + y + 6)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr.factor()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANoAAAAZBAMAAAClTy3yAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEHarIkSJZt3NVLsyme8Q6PJIAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADG0lEQVRIDbVWPWgUQRh9e9ns/eQuWRIJWKibCBELMYSA7TaijXhYiIiQmPiDRVAkoBbqISLY6KEYXBU8RVBsvARb8RAEG5NTJIJYmEoUxIuoRVDiNzPf7M79JTY3xc733vfme7Mzw+wCqk08czlqWRdZOH5bvmU2qrBhkfDTP1vsZlgksvavVruZFqnFFrsBhkXCk247/stTqyYbqeN+I5Y4thDZCSmxi7Jb5RGq2vwa5TX3dQHDNaSGykKgeFlySV92qzxClVV7kL8v/wbWNh7OFiK5FwdEd1Y8wnY/jKqDSDXGiZSngqkhH3AYQLMqxxYE0rceTwvuokrws5lbpFrPSl13UGDrUw0robYg0Lm8XKHOztEjak3c7FwoSWZVWOWGG5zXrITKonvj9gGP07Eyzax/eA8vkHZbd+QDrLGBLF483ExKQ+XQANF03W0nPxLSx0GzkYU1Gr+QyMshQMcJoBenvcsKs5v9FbP+Kx8zVhEJyhiqOI/UdaexOwuc4XLMGhaO6yymfE4nCsBBHHcfKMxuyTyelC4Bcy8XEKOMoUrllFK70VSIeaRI/caGhYVkkZO0fTQxF7clPhYEV4LgOsUbynSMloCu0aszPmFDla4QdoLg5pcgGKQQyNB9OyeCiK2y6PJEUjZRB/ihAMDvNuLRxtA12lXsvfuXcoZKuhHF70YfEoeE0i1iaSRV4DZX0pFcI6T/aMxu4bu9c3Ge5iNWklXtOaVlt848MnTf7lOkngPZlzTjbkXMZSD2/16sgj6Fo31D4Rswst9DGykMVc0pydBU8nWnJBVZjBSeQ1+yoAPdvtRRsamkaOxmf4bjz5YwlVpAh087Eqkyo0IXrqRdxjkfmJdcyBoWPf3d77k4QAtjDfVtOcVqdkPP20Ow58ez7Uc3HaaUoUqWlFSfyV3jb4i4o0i9klUWnJJddCcJqN1MhYgj1U5OaTcJ7Vwjljmzi+5bwa4xU0YcqcaYtXwjHd7KVawh0KHYltVbqKr74sixTb449XXD72R9ymBCVWhrJIGnVWglEB7PlUThIZ5spGr6pwD8A4ssuzZjg3jWAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$x^{2} \\left(x + \\cos{\\left (x \\right )} + 5\\right) + x^{2}$$"
      ],
      "text/plain": [
       " 2                     2\n",
       "x ⋅(x + cos(x) + 5) + x "
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr.subs(y, sp.cos(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also built equations and solve them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANAAAAAZBAMAAACybb07AAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEHarIkSJZt3NVLsyme8Q6PJIAAAACXBIWXMAAA7EAAAOxAGVKw4bAAACxklEQVRIDa1VTWsTURQ9k04m6eTDoRHB3TSCu9JSCt3JbEQQ1OJCXAjRutFFsbhRF0pQF7qRoCiOHxDFhbhq0Z2I/gFpRFCKgnHpylRUUIT4Pu6dvEkzSaB5i7n3nnveOe9NXt4Aeiy98ijbcugr5QRjtS07aIH+Utkg/3NERv2lsgv2r1EZDZByN0ZkBPSXyvrKaN/wdpkggUtSCd0lhdsrCe0u+Kb3to65LpBLLcVVV8w0FDAedOEJ5ff2b2Bn7yZJmU33a1QdwTGZX4wAmTyOVbJwfQ3dmQ0Ah4ouHklponqW3kcHIH/v2arErhrtfkYzkmc1iR1bEEtRT4VUZFRst1sCsqvi0RkxAQXzjpQRbhM1xiOpjorIlNHErr27fYJTDbHO8tzhRV2TgDXl46CnIDaaP/tZ1Pyrs5EpRYo6SCPreOZytkZwbhnYgfP+dV2TgOM2Ma0RNlrFoQXgAk3jBZlSOB3K8VBRpJHjORtuQDOydeAEznhPdE0C+ws1vNQIGwG5KvBUg3wYYlLU0kHtCOMrEVgUy/RwX9ViRTfC8JbIvaKPHyI6YXj3WxjOqHZB3I1rMuvwLFNKkfihfiNs87mGNIISVRC/+4qX/6c5tCNx1TviblRGosE8U0rz6amN1t5EoHx1YNGOwDpSTc0ho2INBfHaj2qwY2RI4dQ1Oa4oijRyvWmkPJohD8OjVAuTuuaVPkdOvy/+wxYAeYC6DkNMihR1kEaV+mtEN6nTQPpvrmULPznY6BMqhNCO7AYuBcAHxYp4MSlq6SCNSuWJj6QCpKuwZienzsUFUFr8EmiIT92Bk+8E8CDOi0lRS4XCnj/zZi3ypCtonXhspEq7SijvnMphQvxS3a6njNUs/vxagaESXarEM1oD05wpxOzscs7n3IwJnwmTkpj3/PCly/Ji2zxebIaGR6IDOHhK4qe859T/bS2hHg5dmHoAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$x^{2} \\left(x + y + 5\\right) + x^{2} = 1$$"
      ],
      "text/plain": [
       " 2                2    \n",
       "x ⋅(x + y + 5) + x  = 1"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eq = sp.Eq(expr, 1)\n",
    "eq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAH8AAAAyBAMAAACHVRmSAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMA74lUMhDN3XarIkRmu5l0i/HRAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABq0lEQVRIDWOQ//+JgVzA9P+/AIOwiyu5+hlYXZwFGETI1g7SyEKmAcyiUGvJNCDy5VfKDGDgGDVgNAyAaWjA0wHfxG9zIGmZcF5g0lCBJnusFEEDWLMZnmHVCRUkaADvA4a7ELXMDdgMImjA/gUwbWQacBGmn4EoAyI6Vzc1wLWAGN+c+jZABIgxgFWBp4I9AaIeQrJ+O8AgTLwBTAFMX5kXoBjwP4DhYQBYBJ8LWNPLgaDMgZWB7QJYcRiIX54CZH9nYDgP9ANTWlq6WFraBJDsfwTALJX5G0BKkEAZ0IADYD4+FyA0+Dsg2GDWRIgLgGxiDGAO6GfggPgYZo4/MAwgbHQD9O6BVKImpP0G9xnWwLRCaE4FVuyxwLSAMwHDgMiOiDcHIBrhpN5LBwgbzQXsC1h+ASVQXQDXhY2BbsAGrt+kGcC6AM1YZlDtRIIL0LQzMLA3UGiAHshIClzA84BsAyB5dhPDFjJdAMmzLMnGV8k0AJJnGf///0CmAfA8S3YYMCDlWfJiASnPkmMASp4lxwCUPEuOASh5lhwDQIEPB8PDAIqb+5R2OACieHsa0A6UPAAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\left [ - x - 6 + \\frac{1}{x^{2}}\\right ]$$"
      ],
      "text/plain": [
       "⎡         1 ⎤\n",
       "⎢-x - 6 + ──⎥\n",
       "⎢          2⎥\n",
       "⎣         x ⎦"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.solve(sp.Eq(expr, 1), y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A *sympy* expression is represented by an abstract syntax tree (AST), which can be inspected and modified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKoAAAAZBAMAAABTBqhqAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEHarIkSJZt3NVLsyme8Q6PJIAAAACXBIWXMAAA7EAAAOxAGVKw4bAAACe0lEQVRIDa2Vv2sTYRjHv5deLtdcEo9WBLdrBLfSUgqut4ggqMVBHIRonRyKxUUdlCAuXeSwVDx1iOIgTi26ieg/II2LUhTM6mQq6iBCfH88z937Nj0Ummd4n1/f93N3712eANqWXocU7dnlKC8eS/aM0wAD5ce1HyOiGih/wf05KqqJqm6PiAoYKD9S1GP/z67EBVpCye6SkrjrBcod5bvhuw7mdxQ51SiZVbqqNh4r98/l2+AXcHB3GaFk8wzOSXdDLpk9ySIOqpGO7s3FgEcJbB2hhLD24PmG1N+WS2a2WpaZOisTpydXYZaOUaLeGAz6wrltseRmqVXZomKNpJZOoyYOHT0cUbskTtdpzp9e1DmpnekIJ0NVYuqRK59Fzq+FqTnKOV+55SdEDZaBA7gWrVhUr9rDjK4wdQOnFoDrtI2vnqO80NuuxtT2O8AFXA6f6pzUx+sJXukKU4GgDTzTRT5XA+VgfJ2a4njFDYR4qPJLaXonTVdFHDYifBfeS9P7X9N0VrXrYm5syijXWah9kZLJRVKhCNLzPaAV1v6oAn8DYsB5Ym4oqqGDgdp8q3eIVZ4AmJBTt1DqiYYwOoFGgrqYG2d1Mbs6MlQ1nEEppLZ8W49LfUzpnN/tCwT6sZlaFzeQDL0tA9XqvEE2Urwuyr+Dvivg0pj6CS2q0L26XdyMgQ9KlekM1GRz4iNtAcptOHNT01dtNSYXvwiGNP4GTlx8L7JHqpZRLRS1lCv6xW6RiKkqddtU5WeidMjZ02W/7o8lDv9ROLGxJZsupDNadhiYu7jlLwcRx6YvmISmRMe7Tu1yU/7oh+3lcKmgkn0OBX2jXPgPA/wF1cGLzQch/MoAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$x^{2} \\left(x + y + 5\\right) + x^{2}$$"
      ],
      "text/plain": [
       " 2                2\n",
       "x ⋅(x + y + 5) + x "
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.38.0 (20140413.2041)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"422pt\" height=\"260pt\"\n",
       " viewBox=\"0.00 0.00 422.00 260.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 256)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-256 418,-256 418,4 -4,4\"/>\n",
       "<!-- Add(Pow(Symbol(x), Integer(2)), Mul(Pow(Symbol(x), Integer(2)), Add(Integer(5), Symbol(x), Symbol(y))))_() -->\n",
       "<g id=\"node1\" class=\"node\"><title>Add(Pow(Symbol(x), Integer(2)), Mul(Pow(Symbol(x), Integer(2)), Add(Integer(5), Symbol(x), Symbol(y))))_()</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"135\" cy=\"-234\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"135\" y=\"-230.3\" font-family=\"Times,serif\" font-size=\"14.00\">Add</text>\n",
       "</g>\n",
       "<!-- Pow(Symbol(x), Integer(2))_(0,) -->\n",
       "<g id=\"node2\" class=\"node\"><title>Pow(Symbol(x), Integer(2))_(0,)</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-162\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"99\" y=\"-158.3\" font-family=\"Times,serif\" font-size=\"14.00\">Pow</text>\n",
       "</g>\n",
       "<!-- Add(Pow(Symbol(x), Integer(2)), Mul(Pow(Symbol(x), Integer(2)), Add(Integer(5), Symbol(x), Symbol(y))))_()&#45;&gt;Pow(Symbol(x), Integer(2))_(0,) -->\n",
       "<g id=\"edge1\" class=\"edge\"><title>Add(Pow(Symbol(x), Integer(2)), Mul(Pow(Symbol(x), Integer(2)), Add(Integer(5), Symbol(x), Symbol(y))))_()&#45;&gt;Pow(Symbol(x), Integer(2))_(0,)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M126.65,-216.765C122.288,-208.283 116.853,-197.714 111.959,-188.197\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"114.99,-186.439 107.304,-179.147 108.765,-189.641 114.99,-186.439\"/>\n",
       "</g>\n",
       "<!-- Mul(Pow(Symbol(x), Integer(2)), Add(Integer(5), Symbol(x), Symbol(y)))_(1,) -->\n",
       "<g id=\"node5\" class=\"node\"><title>Mul(Pow(Symbol(x), Integer(2)), Add(Integer(5), Symbol(x), Symbol(y)))_(1,)</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"171\" cy=\"-162\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"171\" y=\"-158.3\" font-family=\"Times,serif\" font-size=\"14.00\">Mul</text>\n",
       "</g>\n",
       "<!-- Add(Pow(Symbol(x), Integer(2)), Mul(Pow(Symbol(x), Integer(2)), Add(Integer(5), Symbol(x), Symbol(y))))_()&#45;&gt;Mul(Pow(Symbol(x), Integer(2)), Add(Integer(5), Symbol(x), Symbol(y)))_(1,) -->\n",
       "<g id=\"edge2\" class=\"edge\"><title>Add(Pow(Symbol(x), Integer(2)), Mul(Pow(Symbol(x), Integer(2)), Add(Integer(5), Symbol(x), Symbol(y))))_()&#45;&gt;Mul(Pow(Symbol(x), Integer(2)), Add(Integer(5), Symbol(x), Symbol(y)))_(1,)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M143.35,-216.765C147.712,-208.283 153.147,-197.714 158.041,-188.197\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"161.235,-189.641 162.696,-179.147 155.01,-186.439 161.235,-189.641\"/>\n",
       "</g>\n",
       "<!-- Symbol(x)_(0, 0) -->\n",
       "<g id=\"node3\" class=\"node\"><title>Symbol(x)_(0, 0)</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">x</text>\n",
       "</g>\n",
       "<!-- Pow(Symbol(x), Integer(2))_(0,)&#45;&gt;Symbol(x)_(0, 0) -->\n",
       "<g id=\"edge3\" class=\"edge\"><title>Pow(Symbol(x), Integer(2))_(0,)&#45;&gt;Symbol(x)_(0, 0)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M84.4297,-146.834C74.2501,-136.938 60.4761,-123.546 48.9694,-112.359\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"51.4055,-109.846 41.7957,-105.385 46.5259,-114.865 51.4055,-109.846\"/>\n",
       "</g>\n",
       "<!-- Integer(2)_(0, 1) -->\n",
       "<g id=\"node4\" class=\"node\"><title>Integer(2)_(0, 1)</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"99\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">2</text>\n",
       "</g>\n",
       "<!-- Pow(Symbol(x), Integer(2))_(0,)&#45;&gt;Integer(2)_(0, 1) -->\n",
       "<g id=\"edge4\" class=\"edge\"><title>Pow(Symbol(x), Integer(2))_(0,)&#45;&gt;Integer(2)_(0, 1)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M99,-143.697C99,-135.983 99,-126.712 99,-118.112\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"102.5,-118.104 99,-108.104 95.5001,-118.104 102.5,-118.104\"/>\n",
       "</g>\n",
       "<!-- Pow(Symbol(x), Integer(2))_(1, 0) -->\n",
       "<g id=\"node6\" class=\"node\"><title>Pow(Symbol(x), Integer(2))_(1, 0)</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"171\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"171\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">Pow</text>\n",
       "</g>\n",
       "<!-- Mul(Pow(Symbol(x), Integer(2)), Add(Integer(5), Symbol(x), Symbol(y)))_(1,)&#45;&gt;Pow(Symbol(x), Integer(2))_(1, 0) -->\n",
       "<g id=\"edge5\" class=\"edge\"><title>Mul(Pow(Symbol(x), Integer(2)), Add(Integer(5), Symbol(x), Symbol(y)))_(1,)&#45;&gt;Pow(Symbol(x), Integer(2))_(1, 0)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M171,-143.697C171,-135.983 171,-126.712 171,-118.112\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"174.5,-118.104 171,-108.104 167.5,-118.104 174.5,-118.104\"/>\n",
       "</g>\n",
       "<!-- Add(Integer(5), Symbol(x), Symbol(y))_(1, 1) -->\n",
       "<g id=\"node9\" class=\"node\"><title>Add(Integer(5), Symbol(x), Symbol(y))_(1, 1)</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"279\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"279\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">Add</text>\n",
       "</g>\n",
       "<!-- Mul(Pow(Symbol(x), Integer(2)), Add(Integer(5), Symbol(x), Symbol(y)))_(1,)&#45;&gt;Add(Integer(5), Symbol(x), Symbol(y))_(1, 1) -->\n",
       "<g id=\"edge6\" class=\"edge\"><title>Mul(Pow(Symbol(x), Integer(2)), Add(Integer(5), Symbol(x), Symbol(y)))_(1,)&#45;&gt;Add(Integer(5), Symbol(x), Symbol(y))_(1, 1)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M189.812,-148.807C207.002,-137.665 232.618,-121.062 251.993,-108.504\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"253.916,-111.429 260.403,-103.053 250.108,-105.555 253.916,-111.429\"/>\n",
       "</g>\n",
       "<!-- Symbol(x)_(1, 0, 0) -->\n",
       "<g id=\"node7\" class=\"node\"><title>Symbol(x)_(1, 0, 0)</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"99\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x</text>\n",
       "</g>\n",
       "<!-- Pow(Symbol(x), Integer(2))_(1, 0)&#45;&gt;Symbol(x)_(1, 0, 0) -->\n",
       "<g id=\"edge7\" class=\"edge\"><title>Pow(Symbol(x), Integer(2))_(1, 0)&#45;&gt;Symbol(x)_(1, 0, 0)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M156.43,-74.8345C146.25,-64.9376 132.476,-51.5462 120.969,-40.3591\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"123.405,-37.8461 113.796,-33.3847 118.526,-42.865 123.405,-37.8461\"/>\n",
       "</g>\n",
       "<!-- Integer(2)_(1, 0, 1) -->\n",
       "<g id=\"node8\" class=\"node\"><title>Integer(2)_(1, 0, 1)</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"171\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"171\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">2</text>\n",
       "</g>\n",
       "<!-- Pow(Symbol(x), Integer(2))_(1, 0)&#45;&gt;Integer(2)_(1, 0, 1) -->\n",
       "<g id=\"edge8\" class=\"edge\"><title>Pow(Symbol(x), Integer(2))_(1, 0)&#45;&gt;Integer(2)_(1, 0, 1)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M171,-71.6966C171,-63.9827 171,-54.7125 171,-46.1124\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"174.5,-46.1043 171,-36.1043 167.5,-46.1044 174.5,-46.1043\"/>\n",
       "</g>\n",
       "<!-- Integer(5)_(1, 1, 0) -->\n",
       "<g id=\"node10\" class=\"node\"><title>Integer(5)_(1, 1, 0)</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"243\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"243\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">5</text>\n",
       "</g>\n",
       "<!-- Add(Integer(5), Symbol(x), Symbol(y))_(1, 1)&#45;&gt;Integer(5)_(1, 1, 0) -->\n",
       "<g id=\"edge9\" class=\"edge\"><title>Add(Integer(5), Symbol(x), Symbol(y))_(1, 1)&#45;&gt;Integer(5)_(1, 1, 0)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M270.65,-72.7646C266.288,-64.2831 260.853,-53.7144 255.959,-44.1974\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"258.99,-42.4395 251.304,-35.1473 252.765,-45.6409 258.99,-42.4395\"/>\n",
       "</g>\n",
       "<!-- Symbol(x)_(1, 1, 1) -->\n",
       "<g id=\"node11\" class=\"node\"><title>Symbol(x)_(1, 1, 1)</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"315\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"315\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x</text>\n",
       "</g>\n",
       "<!-- Add(Integer(5), Symbol(x), Symbol(y))_(1, 1)&#45;&gt;Symbol(x)_(1, 1, 1) -->\n",
       "<g id=\"edge10\" class=\"edge\"><title>Add(Integer(5), Symbol(x), Symbol(y))_(1, 1)&#45;&gt;Symbol(x)_(1, 1, 1)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M287.35,-72.7646C291.712,-64.2831 297.147,-53.7144 302.041,-44.1974\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"305.235,-45.6409 306.696,-35.1473 299.01,-42.4395 305.235,-45.6409\"/>\n",
       "</g>\n",
       "<!-- Symbol(y)_(1, 1, 2) -->\n",
       "<g id=\"node12\" class=\"node\"><title>Symbol(y)_(1, 1, 2)</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"387\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"387\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">y</text>\n",
       "</g>\n",
       "<!-- Add(Integer(5), Symbol(x), Symbol(y))_(1, 1)&#45;&gt;Symbol(y)_(1, 1, 2) -->\n",
       "<g id=\"edge11\" class=\"edge\"><title>Add(Integer(5), Symbol(x), Symbol(y))_(1, 1)&#45;&gt;Symbol(y)_(1, 1, 2)</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M297.812,-76.8069C315.002,-65.6653 340.618,-49.0622 359.993,-36.5043\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"361.916,-39.4294 368.403,-31.0533 358.108,-33.5553 361.916,-39.4294\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.files.Source at 0x7f72368cf1d0>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ps.to_dot(expr, graph_style={'size': \"9.5,12.5\"} )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Programatically the children node type is acessible as ``expr.func`` and its children as ``expr.args``.\n",
    "With these members a tree can be traversed and modified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "sympy.core.add.Add"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr.func"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAL4AAAAcBAMAAAAtjhhLAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAMkS7zRCZdiKJ71Rmq90icBAQAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADFklEQVRIDbVVz2sTQRR++bHJbpKuRcRr4kk8aApS8SKuGhARYQ5C9VCMtEUEwb20RUG6KHhQoT1o60XMf9CgiEJB9yAKRaWg0IIeqqgtHqQWi0oP8c2vnX12tVXqg8x73zfffp3OzrwFiMLp6YtqWozuYpT4J3QTXiU/5wRtXvLMmtiNSjUP5SDxATdIfU+cWI30uSDVoWTvYLiiSprcir1ImTWi7VyXNXv7yJTUobhE8RqRE6LwlBHPmZJWrvhP4Rpl/4RKAZ/9AGCbpTl1ziXFqCBtvZNJkhh3kG1uwgwnBgBipk9iGlKW6gJmA5FWHcZbPwCOctkwQKGp9emRdLuuad4GWzhB//59qkGkmUPTAa6bb6kbQplp3aaJu1GtOZFThyc6eTFPWO1mSM1Mccpq4JDrgElMMsZbLV3SXGi1FpCxZwmt3QypGeEP+3Ei34BeTNbgzNYxo4tVXRcunvMVzteJUrlZAz5cYUKi/d88PYmYn4n0EuzEdBxe+O+F5JfBGinNuZ4iMyFRKjen2IAhqdD+neKiPkYy9QWOYLoBt9kDqaGjw5ylYqA4t0mUyu1qzhNrRJH2B8jMAjxEBu/9R0wMDuCYEBZkOyK6UDHK67Xap1qNr40VfPiG2TAIctivypjtZbF+EArEK2OD3n08yRU+zb1E6NVWGe5CnGnzwMF+JfwXpb9WSF18LLdHiO8P31EV2r+XnxIRiil4kMOmsAM57LsvMd3LL0C/1NCxyIYgzxTH329Mqf33QGZKKhSTw3vlAfD3i+fnPI7LmQUbH8awvoqkh2pzr2lq2Eliyuht9kA1lHrlb9dhMgA4jWSpAVUG1nT/wDMpgc96tQJ3D3adCdUMpGeJUq+/e+xsICWauXyrD4l9+MP7K9+aFOAYncaIMcXv+kOvkmh/AcVtz/qQGzEGWDkEUUD72zE52eZZ+uOmGMmL/naCQapBTLoIoiATUCyQG2b8BFr25+c485rMhgRRkPh9SQ/yZrMydnPqLf7Eh2DlfBLzt9/HEt/7Yj3Jal24S8Llzrp4JZnwY/qf4yeXiLcck/e0SwAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\left ( x^{2}, \\quad x^{2} \\left(x + y + 5\\right)\\right )$$"
      ],
      "text/plain": [
       "⎛ 2   2            ⎞\n",
       "⎝x , x ⋅(x + y + 5)⎠"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr.args"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using *pystencils* \n",
    "\n",
    "\n",
    "### Fields\n",
    "\n",
    "*pystencils* is a module to generate code for stencil operations. \n",
    "One has to specify an update rule for each element of an array, with optional dependencies to neighbors.\n",
    "This is done use pure *sympy* with one addition: **Fields**.\n",
    "\n",
    "Fields represent a multidimensional array, where some dimensions are considered *spatial*, and some as *index* dimensions. Spatial coordinates are given relative (i.e. one can specify \"the current cell\" and \"the left neighbor\") whereas index dimensions are used to index multiple values per cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_field = ps.Field.create_generic('f', spatial_dimensions=2, index_shape=(3,))\n",
    "# or equivalently:\n",
    "my_field = ps.fields(\"f(3) : double[2D]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Neighbors are labeled according to points on a compass where the first coordinate is west/east, second coordinate north/south and third coordinate top/bottom. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAACkAAAAdBAMAAADFpVh+AAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEJl2IquJVETdZu8yu83OyatpAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABHUlEQVQoFaWQPUvEQBCGnxjMFyEGQbirTOEPEIurFzsL8UARxCZcZ3eNhVpcGktBEBtBPPAnWNgGW1GvsfcneF1KJ7vJyRI7B3beeZ9ZhmHAxHujlrwcWbY1Z21h6f/o6paiO+FwMu5Sp+rRpf63bHNhbSTGu4bth30bR8e3mU20W9n8A7JXClW/nViXE8kemKsNBvhFjU/kHWCuFuX+kI2afkKcNxuHBXPWavoEy2VD10ueSRQ4FYTS1Jf4UNzhTiGYw86Cptwgo/FmcLWgSgYGM04T+X/ZUplbEX9x3xfy1lK9QzDksd6u7sjVdlM3k6FRLlYiTPXVXqeMRilLhYa4Y6NGEmNA9pDwlZZznSWZItLWybRIan5p66b8AFnOOMJB1gZeAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$${{f}_{(1,0)}^{1}}$$"
      ],
      "text/plain": [
       "f_E__1"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "field_access = my_field[1, 0](1)\n",
    "field_access"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result of indexing a field is an instance of ``Field.Access``. This class is a subclass of a *sympy* Symbol and thus can be used whereever normal symbols can be used. It is just like a normal symbol with some additional information attached to it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isinstance(field_access, sp.Symbol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Building our first stencil kernel\n",
    "\n",
    "Lets start by building a simple filter kernel. We create a field representing an image, then define a edge detection filter on the third pixel component which is blue for an RGB image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "img_field = ps.Field.create_generic('img', \n",
    "                                    spatial_dimensions=2, # 2D image \n",
    "                                    index_dimensions=1)   # multiple values per pixel: e.g. RGB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAysAAAAsBAMAAACXlDpKAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAZnbNRO8QMqsimd27VInIquLFAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAH70lEQVR4Ae1aXYhbRRQ+ySY3f3uTaBGhD21stZRql+g+WB+UgBQLSrv2QRGqXhSh2ocsVF9Ebbo+aFEw+qa1bGylqKtlrVURoUYRqaUPUSxUbem1vugqui2ldYV2PTN3Zu7c3Lk/ycZNupuB3XvO/HxzzndmMndOAtAvPcZA9v39PWZR3xxk4Hd4vs9D7zHwGeyo9p5Vi96iQzBeWvQk9CIBp41etGrR23SQMhDJL3oieoIArWKZka3R50pL6//vNgN3WgYspY/Yfd02pz+/xUC0Qp4xMzaKj8fIv37pBQZuJkZsP3XSwMeLRPYqnbtyTh01vCYJW3/lQwTQOU4oenh2Fv8nZvxY6diVM1tNFv0mCtG2ACAC6ExWBA0DtijqbKFjV85M1T/+9pSe0gKACKAzcUE4P14VokLo2JUzU0qdV+C3UrUAIILovEXwMSEktdCxK6duLwX1TMG1CwDCn85rBAebhaQWrCunuq2l2kyhpe6qzgsAwp9OkQ1LXZT9T9ZljcjsytlcHai/7Oox5aoJqLjyIVqmM24yTiLTMjnJvKwR2bpyNtcG6zubu0RqzTVB+pUP0TKd4m01PenLDrty+vYJ17gR7grX0btXaIjvPTHmE0JhRBCdUb5JMnXFaLuKXTntinalxO2nvml3LBsXHsIzLPMKoXA3iE6dXyJzNcVou8q6ctp621J8dpavhHYxwkN4hmVeIRR+BtGZuMQG5QrS6Ef2Quy311eMxc58jbU/rTgptSlFHBAp6u5cp7amABt3vba6sPGmEoC+d21dOZxU/h8QLCxzsaIFCDLNEjhuuFwMQ6fTfe0sAynLh7y5HrKDX0H5ejiO72D11L+umZoqzOMwWImdh0TR2ZCNjYB5rKCdK+kjAH/CeM3ZLmkcojnh1QaEsIJxOk8QWfTxCGwquRJeYejk7jPbzzFmxqWwpBr74NlMEcZNWA6AFx+ypaoShbKYQIUMyOXh28Q7k3ITwJ50MdH4CFKXITINeAnMNbCzswtqMoQr4dU6hG0FC8s8QIyiGzgNnIXBanPCK4hO2X1u+2VG0g4kDB79nJRbNTgHRrlEAoJ3hi00i5ll/ZwbAmOrlwA02Ex+plEDnYaFweAII1fQtPMkIPjqHq/T33KYHGk7FYaGHBAs4dUShNqKX4eHPx4eXudpBVuXljVzgYieMeEHMk28gPm+NLCEl3DBg07LffXEwD+fduQt++h/ktw4DTABSHfiLESR7F2smYeTqlR5F8XYDNwAYLKwsK7kgZHVpyE9CfEC4BQ3AiSrrPnpl4gQMfUKSBDuhFcghLcVbLcorQjpSAiIZXAC0sSrTdXoCGTBlfBS0mm572U7D4v8IQYDFQP2ARxGurXIDAzWIGHitLRYG0JStqKMvB/BjzJ3WNZpgFs7Y6LJxvIqHAV4jg0F+JlImRJcABmCbC1HCYSgs8oQwkTOqRJC9CKzEaVdiEOQM+j37xN4wMIeRMOVJBc1ndR9j4k1/iFWLkhA5YYJDwF++qSLe3Cx44mQHOXNLm8GcKGgOYdxnbjCol3KAiYS8ANxX6qB+w8PqRc4kBWW8ig8CDIENCW8giHorDKEMJGFRQ0hehGDiNImhHYR4g24DUGOkPVnonAQ/6SiptMOi3viFH8Ty5kSzqZUHtPvkWkYME1tBl6pQka0uryJ5gGixdgWeNIdltj0q+RHaPjicEcWytUk3pLQf1aoXXgkfQgyBDTlzIIhKKcyhDCRhUUNIXoRc4jSJgSeJ5k8XIcgK2F9PVFy5w/VdNphcU8s7i3xCuLykrwKYpOA3/Trb43ChqE/AHbzJuqAQ8GzAeC9M1uHMDwOV0mvnXm4G2A1Xn6uhdjYX5OQag6LgWtLhnDlzAIh6KwyhLCCf4gpIUQvYidR2oTA76xwg59AkPSBVWO/4HMpQZSKmk47LO6J8a5hlcGihNMk4olwj6hyeROrizZ3WOwmIg3WQOrNd8tmuRICslUKCDqrhGtbwcOitCKcI8EQuFvwfcYmiP9kxTmppDE67bC4bY/w8zU7Ig2UxQQ5Y2ADwMljWAzms6TIL5q6FwpBHCVnTKwCOgHCxBi1C88Wxy3UP1ulhKAmKa14nEzbVDhESEeCIfBsyZXAet2lc/kmvASdlvsetotoNL//CGf+NvA8h2eE7l5k2MxKauKfApddz+wFPBMhURcNNCwZ8iYWq4hK32yVGoK6JkGEtKIDjlgQ5E2MrFtefBNeNp3SbqnwocBsz/AakYURPZhw4uoPUNomqnFDpEQOmOyOiCnafAV955oG3jwnRSdM7LxpRGt6ZY4Q0J4VHXDEglhDDpZlwi1/waaTuu9heznPUXAl+5SMwRpJOBP3SgqkS0wL9/iEd3vjizF4Kg9TU8YcIegKa9mKDjjCILIH9gO8zd0K9UQ6Lfc9bD9d5TD2t/q8Rn5GcaHbRbdFfN2XlWB5raOLBds1iA44wiCKDreCFE6nl/t4l2el7CCe14pnXUgo/Cgr38lKsDwod9GrVOsaRAccsSD0vOxWsFynXTzd5x9H5EboC+ZgTrZBq/mOczU6HIjQ5u5BdMARC+IJl5/+FRadXu5jak2U+4WkElhg3U1Rw13nW+NYobTnAoBY5euyu1Gm0+1+XFouS9yD+zXdYWClNO1AQVL6YhcZSD0gTa59Kil9sYsMbGvIk++uylpf7hoDXzpmTo061L7SJQa0kpj4P197oU2OwIdhAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\left(- {{img}_{(-1,-1)}^{2}} w_{1} - {{img}_{(-1,0)}^{2}} w_{2} - {{img}_{(-1,1)}^{2}} w_{1} + {{img}_{(1,-1)}^{2}} w_{1} + {{img}_{(1,0)}^{2}} w_{2} - {{img}_{(1,1)}^{2}} w_{1}\\right)^{2}$$"
      ],
      "text/plain": [
       "                                                                              \n",
       "(-img_SW__2⋅w₁ - img_W__2⋅w₂ - img_NW__2⋅w₁ + img_SE__2⋅w₁ + img_E__2⋅w₂ - img\n",
       "\n",
       "          2\n",
       "_NE__2⋅w₁) "
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w1, w2 = sp.symbols(\"w_1 w_2\")\n",
    "color = 2\n",
    "sobel_x = (-w2 * img_field[-1,0](color) - w1 * img_field[-1,-1](color) - w1 * img_field[-1, +1](color) \\\n",
    "           +w2 * img_field[+1,0](color) + w1 * img_field[+1,-1](color) - w1 * img_field[+1, +1](color))**2\n",
    "sobel_x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have mixed some standard *sympy* symbols into this expression to possibly give the different directions different weights. The complete expression is still a valid *sympy* expression, so all features of *sympy* work on it. Lets for example now fix one weight by substituting it with a constant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzUAAAAsBAMAAACu84sRAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAZnbNRO8QMqsimd27VInIquLFAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAJcklEQVR4Ae1abYhcVxl+Zmf2zs73WAmCP5oxaonRxsH9YaW0zB+xoCTbFCxCbQc/ILbgDFT/FHUn8YdGC46Cgqllx0b8io1rUz+KVAeRUGN+jGIhakquKWJcxaZx2WaFur7vOffcOffecz/m7mAWMufHPe97zvs89/2Ye865dxeYtZ2YgeoPvr0T3Zr5RBn4Oz47y8MOzcAzWO7vUNdueLeexkrrhk/Cjk3AxfaOde2Gd+wpkYFM/YZPxM5JwGHHlepACG/cOZ7NPMk1ZQ5eL7rc+2cZ2UEZ2CuLYue6JHyUL7O2UzKQFTvMwy9caJNHn0/u1dpZBqRs037VnQZfSo6UMEPigkyFDTZ7YGuLrvlNA8Q8VO0vOKuheT56dNqvutPgS8mREmbIj4FJe1ayPQPEPFTqT1DIAMW0X3WnwZeSIyUskBLAwLTcdu1W+q4YJ5RahfU4m/D5ab/qToMvJUdKmCE3BqZsw7U76UoJhIpYDBMYGk2m/ao7Db6UHClhhrQEmDJLrtUhV0oglBoJjEJN5Ktu6PTEE9PgS8mREmYIMcA0XpsKr2j21qmbu6xmb8k/CSwMtSkprgVGQgeCaOdVNxThn4jzJgnfNDgMiUhyaxnOg3/5mxDe0X5kBDzqjxEGpvuVUeaKkqjPjqxnWJ3f2rKpNnWW9ZYZ6Fq0HETLV91olD4b500SvmlwGBKR5NYylDvwUJ+lB7b+S9cjclC7Gpjcv9wUVzXDXcCLrGYvv0UbHYsH8J6xMqmUk6+6yWEx3kTyVRryPtvhwO/DfI28tQdEm0e1ySN3HhUl8kySYmJ6uzIqDZVE/XOA2JoWtDFNzN/xwhlNnVB0XnWTo2K8ieRTtdkOR3htIm/tCbDcgzw/9TzDrmJiuqiqWBu4dsA14OCI9JDa0FJ3RbOeUJSvuhOAYryJ5FO12Q5HeG0ib+2JsNZD7lUe6XmGXcXEtMwl4FZriE5crP9QbVhfeM2+Lj70BHKXv77naO7Sr2noT3suCJukl+Rossw0K8EPrtvyxqnNtjhUbZKEEhZEZ4DcVc7ZkZu+A2tfAweOfWVv48DbWkDliVuHPBNonboz5AqkF2i7WrZJKLfpIGe/G9Xyr9B5E54HqsMCFW6ClhxtP49yL7cO/8e6bXnj1GZbHKo2SUIJC2LFRp4PAThDf2Su5pZgn2tYV1uVJeCfWBkYE+o+LiuqSGTG6Re1Ifm3hdEJfLrUBLHvFtsQP5oFIxkNdn0TRnTfZyRgbFmr4zeBj3Xx3hj55D1UbWIjiuBwahMXSp7uGBoEpV/Whn7wq48Xm/nRaX4E6HRM21BtRFDNAWaiVrNFR7Xg+Q//nNs73RWAhn74NVxFu9Piw8GjwH3ym6hC4WEHLjpS/iAE63bmOd0FLAO6qjD5Jktzl2yGWTjE/1MycD7WaRSUV7HCkq3RGx8f+yTdqC4u3n7P4mKPyBNwmEPBi4uLP15cvC00FAnbvx+VVkQQ7poGFDfbtYZlrXNVFoaYH4p/paEgJBMlhZmouUeA5brQ5eUaPSBcKzphv9Ti88VF4CQoefmXMbdK+1DfMX7kixqKlaKaUeMG9DFnLv9dogJuxnkBy23izYAd/FgX542Pz+OG89zw6SYyomOICMU5QxtDkbCMXenhexRMWBA1OqfxgkMf8Kvr/DuvXEFxFfMNUOLfShPKAZEUZtJqo69pOAfsbtPsPfRLbWd7bZwAfkH3tTKbKA+Az9CcbH9WAvekWD19gOQgOm8rk4qozdOotRlGrj5HiwJtjr6PdTHe+Pk8bqjaJOAID8WpjTkUASu1sIHDFFdYEOUhMhzWfBPFDdxmodhEycbBfnt3H2fpLzS2SB9ZcFKYCejYoiOh4Qjc0ZsaVQMYAKfRGdm4H/QMFpuPU7V5cfwcT4rmDwjvUjNOH0QvdJWJqI31CuZHDCv3qP68PpUaykD2Md74+GSMyg1VmwQc4aE4tTGHImCdLj6IbD88CHr3pGpQ7Si8pvVqFfM2aKM4URjRkkQPFAchHeCkMBPgPi41m1WnlUbWveg0cRyZD+BgoU6nNdq1srZtbeJLhCO40wIBvUHNOH0QXXItRG3oT0GlOgg218zdh0/Q5JprIIUYb3x8MkblhqpNAo7wUJzamEMRMNoof4S5ekQQz9Lf/K1/ozDArn7uypf5XwLpaHVnFZ3+wiYVjEId14aZnCMAC/M9vjrNurS3i1IPhVNfaNNLDnKrmKMF85td3LX/H3QaiajNecXh9EH0Y66FrM0GPygM+/6lw/vJqcDHuhhvfHwyRuWGqk0CDpkaxzlWFIc6Q5tDEbDlNp5CpRcRxMf30GLzM+D4Zfrn8yN1vBfYS2+Lr0Pu6L8onRyEdICTIphwsEWD3MpN2cdez9KGN3StJJ+jsvI+d8ooEHpsoZ4b2hHHg5joY12QT8ao+CoNoxueQYcjPBTnufGAwLsE30XA6Lk5pKdloiDKA41J7DcywWLL55tWl/ga1/K87yBHz9A5avRRTXimK56TqI9Oou8CLjC6LdwA7Te1lnYWn+BjnZFP+qTcsPo+H/yqyxERysf8INq7RSIoFAmj/WYd8pVA2E4QRJf3HbhMIimS6a/qtv7DkRr39i+1yz1ybOiOBn5sfJOwJtGfcqfFcwM+pwnf5PgEH+vMfCJbUW6492dhzJEyFAEr8Tkt1yM+2ZIHUd3gL8vgpEgHOCm5Hqm8AspmvaykqP78TU/StEVwp9FnnG+0NYXeVsKbRD/kGlSWUDiDfbywR8Fce79g5uNPS8n5xhwpQxGwuQFtERnb72ACvXJk34jMOCnSAUqKZNIqwtVL2n6iDI//8ig+WXc0VvAtNRXal9rOVOHktUb+blRP0e4YD0vMl9ANDx/5lDIUB7a21kax5eGcSHEd4KRIpgod35y2SwkJ+ls9Nlx3tzVdKUyY89hXpFk8LIwO0+BTHB7X4n3ywrKhLsZPGJiK4/t3PH5Fs5X16Upf0yp1TQkRh/r4H4WSBKajPPJQ11LySY6UoTiw3+l+TCp7HBBM8+NUzjWT03lSmdFxD+pKiOyJQTqQBBbCBkyDT3KkDEXCrEGohwkmdAckE31cdNu9rhQvyF+nwe4Ww5h/yPPrlJNJYH4aV58GX0oOHTbXdj1KIQSZ7tZYXqvJM/G6Z8CzjmUb192fmQPjDHy1P5Zh/VRTZuJ1zoD1rMeBx/RKeWZmyv89A2V5SFL3LXSVNOuvewY+Mvbgf+qq3TZRfHcjAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\left(- 0.5 {{img}_{(-1,-1)}^{2}} - {{img}_{(-1,0)}^{2}} w_{2} - 0.5 {{img}_{(-1,1)}^{2}} + 0.5 {{img}_{(1,-1)}^{2}} + {{img}_{(1,0)}^{2}} w_{2} - 0.5 {{img}_{(1,1)}^{2}}\\right)^{2}$$"
      ],
      "text/plain": [
       "                                                                              \n",
       "(-0.5⋅img_SW__2 - img_W__2⋅w₂ - 0.5⋅img_NW__2 + 0.5⋅img_SE__2 + img_E__2⋅w₂ - \n",
       "\n",
       "              2\n",
       "0.5⋅img_NE__2) "
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sobel_x = sobel_x.subs(w1, 0.5)\n",
    "sobel_x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets built an executable kernel out of it, which writes the result to a second field. Assignments are created using *pystencils* `Assignment` class, that gets the left- and right hand side of the assignment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABA4AAAAeCAYAAACsR1zIAAAABHNCSVQICAgIfAhkiAAAFNJJREFUeJztnXm0HFWdxz8QDDAEAWVYBZ4BgkRAFllEhIcYgfF4cJxB1AHmMQfjgGwH0ZFFCfEIhx0cYAgCRkZk8ABhGWUCAhk2WWSRYREUfBFk3weTEJK8+eN773S96lvVVdXVXdXdv885ffq9e+tW3arfr36361e/+7tgGIZhGIZhGIZhGIZhGAPIisBQ1Z0YII4F7gfeBl4BbgC2qLRHRrcw2Rt1xvTTKIrpTn9h8jS6Sd30bTKwXItt6tZnw+gKywGzgDWr7sgAMRc4CBmYLYE5wIvAB6rslNEVTPZGnTH9NIpiutNfmDyNblI3fdsAOL7FNqX2eUNgDLimSOM+Zwhdm9kl7e8I4HFgodvvUSXtd1A4Bti/6k50kKORXnyl6o6kMAlYCny+6o70IXWXv8l+sDH9NIpQd70B0508mDyNbmL6lo3DgANybD+uz8vnPNi27vvBnO08XqhfLdi+Sp4AHu7Ssb4MnAssAs4BTgLuoXznRL8yBPwd8NMOHuNDwKXA88C7wCiS1Ro59zOKZBr6vJjS7uPu+4Gcx+smqyIb83rVHelD6i5/k321mH1Kx/QzjOlNa0x3smPyHDyqtCGmb9m4AL2cfn/G7cf1udU8hzgzge8CfwPcmLMt6EHuH4CPAE8WaF8lJ6N5H5OBPwbqh1z5T4CRNo/lr9P66ObrxDH6mcuQ4Ti3Q/vfGLgbWAu4DvgdsAOwO9LrTwKvZdzXKLA6Mqxx3gHOSGi3IfBX7nhjGY/Vba4EpiBjvrTivvQbdZe/yb46zD61xvSzGdObbJjuZMfkOVhUbUNM37IzAzkDvpdh27b6/AskjLXzNnQ8Afwv+SMd6sD26NyPTqgforxogFsJK32Zx+hX1kfTOzqZ22AuksPhsfKzXPmFOfY16j79xunAC8AmVXfE6Dom+2ox+5ROP+vnCJLxcIG2pjet6WfdiTNCcV3qFQZJnt3AbEg6ddK3DyMHTKtIkEx9XgE4EngEPYDNB76NohNeAJ6Lbf8plPPgaRRa/zJwH3pD7zmV5JCTXpqH/ixwR0LdEOkP9TsCV6EQm8VuX7OA9SLbzCD5OqXVjUT2MQJcDTyD5Pc2cBfh6xzt8xTkVXoZWIYGi1b1ni8BtwNvuWP+D4rOWDF2vEnu3O+Kla+MdGeM5nk3h7ryfwr0P8S3kO52ismoP3+k2QG2KroR/wKsknF/o+Q3jru7PpyWUH4GmlZ0LQotegslN1nHbTcV+BmS5VvAfyJPbZyVgO+gXBuLkC04DpiAHIBJ1/lM4CV3HKN8QvIvW/ZQTP69KvvPoOs3M1a+Iw07OxSruxzZws063bkc1NU+mX52hxGKPezVVW+i5TaudZcR2nMcmDxb0y/jDlRvQ0zf8vMsMD2lPlOfJwI3oYv8EPI0XIoeBn/kyq+PbH+cK5sPXIScBZcAv2V8HoQvoYfPMfTQOCPy2aDFidWJf0VhGqGIiyGSHQcHAUvQTXMFUuw5bl/P01DOYXRNRmk4C/xnGIXsjKFcC9G6rSPHWojC9GcDpyC5POfafT+hz3cAbwD3Amcjr+C2GepBMh9DS3b8G9KZR13ZPOB9sWPeCbyHDInHG8/Q9bvKlW9ENu5AUzk6xcGoP7MS6r3HdY+M+xtFDrn90f10JDJ0E1LafJNwAhhffj2wADmQTkeOnDHglyi5yTtum9PRkiteVlFWAX7t6n6DnH+zkaH8d1ceus4/pH7Gsd8Iyb9M2UMx+fey7LdF53R2rPwaGrZp20j5esgJel1Xepeduton08/uMEKxh7266k203Ma17jJCe44Dk2dr+mXcgeptiOlbfq5HDpIQmfvsnQPfZXz+g10Z/+Yb9PC8BD2oTQzsKx4qPt21T/Nu1B3vuQqdwxDhB98p6Eb/Awqjj/Jp5DyYEyuf5/aV9RhRNg6UTQRuQQ/s0T74/Y0xPkIka/0nXN2faHjtQFErN7i642JtZrryz0XKTkG6dCvygHmWB15F0SxZmIhu4E6uQHE66v83E+rPc/WHZNzfKOEokmeA3RLaXO62iXucffmLwFaR8lWR42cJMsSfiNRNRPJbhjypnp8StgV7RPp4ZOz4F6AIl08jffCfSQnnYRQjJP8yZQ/55d/rsv8wOqeLI2Wb0rDP8R88P3Blu3Srgxmpq30y/ewOIxR72Kur3kTLbVzrLiO05zgwebamX8YdqN6GmL7l5yR0DeL5DjP3eQd00v+dcIDHXb1fQmI39/8lGTt4odv+4602rDET0INsKDHkEOGH+rNpflCOMgcpbvQN/DyKOw6S+KJre2Bgfy/SPK0gS713NIUcKVOQ8XsmVu715qxI2X0omuEbrm6KK/fe2IuCZ9TMVLf9fhm3L8JF7hgHJ9R7w35sxv2diG7OtVFSly3QvbIMeUc/FmiTlCvkCXfsvQJtHnJ1oSkrN7u6D7r/d3L/X5vQ5/mu/lOx8pCRjzobjXIIyb8s2UMx+fe67FdH/b0yUnYR8CYNB+m+rnwlFGF1Tzc7mJG62ifTz+4wQrGHvbrqjS+3ca37jNCe48Dk2Zp+GXegehti+paff0Z9WS9WnrnPPswi5MkBRRaM0XhjvSZS7jEU7rAf6UkW7kNv3kMPoGugkIjQ2/JOchXJyQ6T+DFaYmS1WPkQ4Yf6e1z5aYyfXuA/d7n67SJt5lHccbAhcD7KZrqAZsFHb1q/v7kJ+2pV/4CrT0qa4W+k1SNlE12/fuv+Xw05Tk4GNnfbe4/kMe7/LyfsP87ebvu9U7YZJfmmCH3iSzq2Mo5+6sZ3MvY5iTPcfuLRKKsgh0w814YvjztqPK+jbLbxqSMATyHvosfbgu0T9vUgMt6rJtQbDUZpT9/ihORfpuxhMOW/HLqGv3T/r42il05Fy0uNAV93dV9z//99CccdpVz9qKN9Mv3sDKPk053ZKfuqo95Ey21c6yyjlKdLYPLMSlXjDvTX2GP6VoyvonPaKU+jFSJ/fxZdwNsTtp2MHu7/7P5/FYXLnIiWZ/w8EtDNwPGMz3GwArAlilp4N7Dv44D/ojkk/VCU7G5d4DEUgp6UnDBEq/Yz0UP6JSgZRhbmIE/s51AijVZ4b9W3WmxXRtjKZOSgWQOd503ovJYiJ8A/EnbcpK2rnFbvnScvJNS/gBwZqyEnE8h5dCfKa7AW8qpOQFMpnkA5H/ZA+RJ8+M+tLfrn8WuSxn9oRvFJPLPyfOx/rydxx1G8D1n1KYkLUcjXrrHyjyGP6oMJ5TcH9jWEdOIaNF0lyiTksIsmrPS24DcJfVsX+D3y7hrptKtvcULyL1P2MJjyH0P3rB/0j0Q/6s6lcZ5rROqeRtfUcyyK6toMjXH3uLJHWxy3bP2oo30y/ewM5zDeKQ/Kd7QPmmc7Gqt7OGVfddSbaLmNa52lTF0Ck2dW2h13wMYeMH0rin9Wen/qVjG842Al9BDnwzbi7IxCGeIh+o+iMJqJSIjT3f/bo8gE7ySY6o4RFyooBOVrNKZAePZDN8+h6EHzEHf8qWjuSSuytH8EZQDdH72lz8JNKInG35LNcRC9mdIeaMvgaOSoOIhmj/BXkOMgREjmWer9ua1DOA/BurHtPLcC01AY0s5IT/wNehuKGFgRhfw8hrKcZsHn2gg5pzxZE7Mk8aT7npJQv6n7fqrN4/hzjmeg9YlyHkooDxm17VLqtmG8wfW24GHCcv8Ikve8QN0aKNJlZ7LnpehFrkJrFZ/VakPa17c4IfmXJXsoLv9ekX2a7N5AP+AmoRC+n9H4sbQU+AD68fBR4DD0dsEzjOYF3o9++M0EfoXGm9dT+lO2ftTRPtVBP6E3dDSPbQmtbz6CHvZmE74GSdRRb6LlNq7lpypdApNnnE6NO2BjD5i+hchy/3vHUShPYSJ+LshS91krYbuT3HfowR/0FvlXaPWEO9HDa3TlAZ/1Py5U0Jv7Za5dlKORwfoRehN9BHqDnTWxRtb216NwjawsQtERe9OcsCmEn48Un/dShKXuOymrqJ8ycHWgLmkKSjt4eQ4n9OVDyDHzZqzuFve9B3Ie3EVDgW9BhvIQZBhuITsL3HcnQ4Vuc9+fJbzkzCfRyhbtzkPziVriIVbeEMbvRV/+QGBf26XUbRPb3xL3SZp29O2E40Ny5FBWdkVJNZ+n3JC8svswEziBZM96JwnJvyzZQ3H5tyt7KEf+hyKbswidc9zupsnO/4Cbjt6+nRGpexvZpaPQW4cfx9ru6coeRdmaDwD+GtmDblJH+1QH/YTq7VOdbUsd9SZaXqdxrQ7jVJZ+1G2cipb3mzyrGnfAxh7oLX2D+owlPtLgL3kO7oX7HgqxWJ/mN///gsLKoXGBtyGcj2ATlMDiT2gJQI8P1w+9cd+Fxlx5z0Qk1Jti296EvDatyNP+XhQhsXKG/XrmoIfaaRm2PQ9d37MJe+Imkt2p8Aa6TknrWo+67+FY+Z4kzztqh0vd9wnIUHkmIOO3POHkmQ8gZ8I+yIsadQ74v30uhqzTFEDTZyBn2E1OnkZ6NISSOUY5CenFZTTfiBsjr2R0LtVH0aAQZyOkN9A8r2xbNDg9HihfTDg8zRvIkFGLG9wlyCO8Ec0e6W+gaJbQvnzk0MUUZxWU++KwNvbRLln6EI1U6jYh+Zcleygm/zJkD+3L30eZnYzGqLtQlFnUXqbJ7g30w+AoNOf0sUjdW2ht7b3Q250FTa3Hsyqyf2lvfDpBHe1T1foJ9bBPdbYtddQbX163ca0O41SWftRtnPLl/SbPOo070FtjT8h+QH4b0kv6BvUZS/yz0qsp26RyIHooXYyEexqaL78QOQLG0EUDebeWIc/RJWg5vf9ASr2A5gvrl3N8Di3ZMYNGptDraF73cj23fXwe3PdohMOkkaf9Vm7bPIkZV0Ph8FHv3xDJSWP2R9f1PeQlOhOtkXkt8iL+Lrb9PJKnB/waXfvLUX6JE2gsMbKV69ciV38aMkTLkHzGGJ8VM63PWepBSVzGUP6L890x/fqoSUt1gs7dJ13ZMVb3B1e+hHzecp9M5sBWG7bJxuh8x9B5nIIcHGNIvz4YaDPq6ociZTOQrG5Eg8KpKLxoodv2F4y/fisiPbovtm9fnjT36hWa5yh6HkH3bDSKZX93/EXo3jwVDYYvo+idZTR7XvdFA1V8WZeiVPkmJ0sfTqR57nWnCcm/bNlDfvmXLXsoJv97UYRZlN+j+zNKkux+TsMmxSO0/DS+hSRH5kW50rVJW7e+U9TJPtVBP6F+9qlTtmXE7Xu4QNs66U20vM7jWh3GKUjuRxW6NGjyrNO4A7019ozSbD8gnw3pZX2DaseSI1z7+NSxXByJwj/eQ8K/Gnle5jPeI/EFlGHySRRFsBgJ4GIa81jiHOa2X+Q6+gNXPheYFdvWP/jH38SfSPNDdog87Td1226RYb9RbkTXxCvWEOkP2Vu6uvno4f515AWbhcL1o8wj2XGwCXI+vIYUdQwZeM/O6EZ9AyXpuBPJa5jOOA5Aqx7c6Y63CHlMjyd9KsfhNBLDxA3cLFd3b4vjhngJ+H6BdnnZADmOXkD6Px95nUNeUggbyN2AK5BOvonuu1dQIpcDaTY227l9XJhQHr+PQM6+McLTV1Z0xwyFjh2OPMiLkePwfLev1wh7aM8lnICmKHX4QZbWh73QtckTqdQuIfl3QvaQT/5lyx7yy38icjTuGys/n+YlhpNk5+1O/AcINGxylqVhT0d2IWm1mW5QF/tUB/2E+tmnTtmWEYo7DqA+ehMtr/O4VodxCpL7UYUuDZI86zTuQO+NPaOEHQd5bEgv6xtUO5acR3OURk9wOfK4RclzM4bI035HdOHXJh/TXbvdc7YzOsvPkTfSKJ8DkM4fE6gLRQ61Qx1+kKX1oUikUq+TJP+yZQ/55d/pKLOsnImcl1M7sG8jnV6yT2Zb6kM7elOHcQqS+zGIutRNedZl3AEbe6qi3XGnyrHkbrIvDPD/xBNYVMFDNCv6YjQPPp5DYBo60Vbkab8lSjDxUpbORrgOvfH/Ys52Rme5kkZiEyM/EwiHxH0GLZH5LM2eXVB0SXxpnxk0QvCSPsMl9DlEp4+90H13M+KgGxSRf0j2UI3845FaywXKOiW7H6Jwxt3pUS9+D9Av9imNfrUtVVKm3pTJDGycKkLd5FnluAM29nSauulbVtJ0bnnkWIi/uG/JCq036Thz0VyRNRk/HeIsNB3iPjRH4+vIuxcSTois7XdFGS/z8hLVzCEy0rkBJWbcgXDol5HOVLSsz1wUkvU+tCrKLuj+3ActRxrnVZrnd52HcmukkWVp1SJ0+tg+7O6VNvZRR4rIPyR76K78X0WrzqwTK1+LZqdwJ2R3Afrh9gU0Dc334x3C94tRjH6xT2n0q22pkjL1pkxsnCpGXeRZ9bgDNvZ0g7roW17SdG4asi23d6875XI34eyQh6I5MO+iCIJ4ONAI4fkxWduvjHI07JS7x0adOQAZUyM/mwHXAH9GntIFKGfFqaRP5zmGcObaotQhBDStDweja9RvFJF/2bKHYvK/l+a5oE/RnKSqE7JLels4o+TjDDr9Yp8G0bZUSaf0pg7jFCT3o191qU7yrHLcARt7ukEnx52qxpIrqGbFldLYE91oed/gn4SEVzRy4jCal2w0ep/lUVKZpIRORvlsiTzva7axj0nIi7s1jTljW5O8/GgnyNqHywgvNTqIlCF7aF/++6FpagcDmwPnoLcAG8W2M9kNHnWwT2Zbeo+Q3tRhnMraD9Ol8XRCnjbuGEkkjTtVjyXrAbdR7kpYlXAEzTdaK+6nveSE05Enyeg/NgfOrroTA0ZS5FBWhgl7zmdHthkhPcqoXbL0wSKVmmlX9lCO/C3KzEiiavuUpb3pZ/2I680w1Y9TWfphuhSmE/K0ccdIIjTuDJOucyOk61ur9pCucz8heQVEwxhopmGJErtJ0cihPLQbZVQGFqnUTDdkDxZlZhSnF+yT6Wf9KKI3Nk7VlyrkabIYXKqyH0k6tytaptEwjAQsgWV3KRI5lId2o4zKwCKVwnRa9mBRZkZ71N0+mX7Wk7x6Y+NUvem2PE0Wg00V9iNJ5+yZyDAMwzAMwzAMwzAMwzAMwzAMwzAMwzAMwzAMwzAMwzAMwzAMwzAMwzAMwzAMwzCMnuD/AAhkq0BVIUB0AAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$${{dst}_{(0,0)}} \\leftarrow \\left(- 0.5 {{img}_{(-1,-1)}^{2}} - {{img}_{(-1,0)}^{2}} w_{2} - 0.5 {{img}_{(-1,1)}^{2}} + 0.5 {{img}_{(1,-1)}^{2}} + {{img}_{(1,0)}^{2}} w_{2} - 0.5 {{img}_{(1,1)}^{2}}\\right)^{2}$$"
      ],
      "text/plain": [
       "                                                                              \n",
       "dst_C := (-0.5⋅img_SW__2 - img_W__2⋅w₂ - 0.5⋅img_NW__2 + 0.5⋅img_SE__2 + img_E\n",
       "\n",
       "                       2\n",
       "__2⋅w₂ - 0.5⋅img_NE__2) "
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dst_field = ps.fields('dst: [2D]' )\n",
    "update_rule = ps.Assignment(dst_field[0,0], sobel_x)\n",
    "update_rule"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we can see *pystencils* in action which creates a kernel for us."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pystencils import create_kernel\n",
    "ast = create_kernel(update_rule, cpu_openmp=False)\n",
    "compiled_kernel = ast.compile()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This compiled kernel is now just an ordinary Python function. \n",
    "Now lets grab an image to apply this filter to:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6IAAAEfCAYAAABbM3sFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XuQXOd53/nfMxcMMAAG1wEwuAMkeLUjUoJlyVI5tGnZkrwVOlkrKznJcr3ccLdKduxNaiPZla1VqnZr5Yod7aUSVTGRYu6WI0uWraXWKztSKDEXK6IIirQoihdAIAgMrgNwMAAGwFzf/QPN9zzd6LfnnOnuMz09308VCk/3vKfP22dO98yZ/p3nWAhBAAAAAACUpWepJwAAAAAAWFk4EAUAAAAAlIoDUQAAAABAqTgQBQAAAACUigNRAAAAAECpOBAFAAAAAJSKA1EAAAAAQKmaOhA1sw+a2WtmdszMPtmqSQEAAAAAupeFEBa3oFmvpNclfUDSqKTnJH0shPDD1k0PAAAAANBt+ppY9t2SjoUQjkuSmf2hpEckJQ9Et27dGvbv39/EKgEAAAAAner555+/GEIYXmhcMweiuySdcrdHJf1kowX279+vI0eONLFKAAAAAECnMrM384xr5hxRq3PfbTlfM3vczI6Y2ZGxsbEmVgcAAAAA6AbNHIiOStrjbu+WdKZ2UAjhiRDC4RDC4eHhBT+hBQAAAAB0uWYORJ+TdMjMDpjZKkkflfTV1kwLAAAAANCtFn2OaAhh1sx+TdK/kdQr6fMhhJdbNjMAAAAAQFdqplmRQghfk/S1Fs0FAAAAALACNBPNBQAAAACgMA5EAQAAAACl4kAUAAAAAFAqDkQBAAAAAKXiQBQAAAAAUCoORAEAAAAApeJAFAAAAABQKg5EAQAAAACl4kAUAAAAAFAqDkQBAAAAAKXiQBQAAAAAUCoORAEAAAAApeJAFAAAAABQKg5EAQAAAACl4kAUAAAAAFAqDkQBAAAAAKXiQBQAAAAAUCoORAEAAAAApepb6gkAAIClMzs7G+u+Pn4tAACUg09EAQAAAACl4kAUAAAAAFAqMjgAACxjPlo7Ojoa67m5uVivXbs21tu3b69a/rXXXov1/fff344pAgBwmwU/ETWzz5vZBTP7gbtvs5l9w8yOVv7f1N5pAgAAAAC6RZ5o7u9L+mDNfZ+U9HQI4ZCkpyu3AQAAAABY0ILR3BDCvzez/TV3PyLpoUr9pKRnJH2ihfPCAq5dux7rk6ezKNZbN+ZjfeHKTFa/NVG1/GW3/OSNqVhPz2QRr7n5EOu+3uxvFqtXZbvN+jUDsd44lEW/RrZsiPXmgSweNrxpY6zvOLBXAIDmjI+Px3r//v2xnp/Pfh7cuHGj7v0SnXIBAEtjsc2KtocQzkpS5f9trZsSAAAAAKCbtb1rrpk9bmZHzOzI2NhYu1cHAAAAAOhwi83jnDezkRDCWTMbkXQhNTCE8ISkJyTp8OHDITUOmevXswjVX7x0LNbfPpbFry6ELPp69vJkrGdms00cgq+ro1ghDLgbq7JyPrW8qzVfd7zCVTf+St1l164+H+udg1mnxp+5u7rf1Yffc0+s169bJwBAfT526zvlTk5O1hsOAEBHWOwnol+V9GilflTSU62ZDgAAAACg2+W5fMsXJP0nSXeb2aiZPSbp05I+YGZHJX2gchsAAAAAgAXl6Zr7scSXHm7xXFa0F149Eevffy47l/ZklnDV3ExvrGenL7qlXWw2Fa1VdSraf02ug6Ifl4zm+jhuavx8/XVPTGXrmpiwWL92cbpqfn/28qVY/8bPH4r1g/ceFAAgY5a9l6bet73U/QAAlKntzYoAAAAAAPA4EAUAAAAAlIqrWC+hV4+fifVnns26G45fchcen5uN9dRM1g1xzsdpXcyqvzeLaLlSqo1iJWK7fliovuEXrjsmGQUOqcd3z+fGtarpHbuZTf4fPXU01p9ykd+fuP8OAQAyRHMBAMsFn4gCAAAAAErFgSgAAAAAoFREc0vmI1F/9P2sM+z4W+Oxnp/PIqszs1n9wL5NsX7XgS2xHruSRXn/7UtZ3HfyZtaJtrfH53Rr41v+/nnV+0IygpuK7xYdX9vV10WPxy9l2+n3vpF1Dv5nI9k22LJ5owAAAFaq69evx3pmZibW/f39sZ6dzU756u3NfqcaGBiIte/E7cc04h+3r4/DC+TDJ6IAAAAAgFJxIAoAAAAAKBWfnZfs3LkLsT72VhZjmJ+t3x33r967Pda/9ciPxbqvp/7fEA4f3Brrf/zl52M9OzdfNa4qqFswUpuK1xYfr/pjbrud1afOX471//Mfvh/rxx75aQFAu8270waOHs06em/evDnWw8PDpc7Jo2susHK98cYbsfbvSdu2bYv12NhYrH0E14/x73M+slsb0z116lSs9+zZs9hpYwXjE1EAAAAAQKk4EAUAAAAAlIpobsnOjV2M9eXJqVj7oJRvcPufv3tvrFNxXO+n7sqiFfftyjrJPvejC1XjVvdl8Yp0pNYt0ExH3KLx3UbLzGbb7NvHs5jur7oYSU+O7QQAtXwc7cyZrAP59HTWgXxoaCjWd999d6yfe+65WC9lNBfAynX//fcvOMa/zx04cCDWr7/+eqzvuuuuWJ8+fTrWu3btqnqswcHBWN+8eTPWq1evzjljrHT8xg4AAAAAKBUHogAAAACAUhHNLdnUjIvjVsVPs6hEcEHd3iZipr2uG9rtyddicdx2d9CtnWByeRcpuTKbXaD52rXJWA8NrRcAFHXlypVY7969O9ZTU9n7to+f+fvzXvS93eiaC2Ax+vrqHxL497kLF6pP8/Kddq9fv96eiaGr8YkoAAAAAKBUHIgCAAAAAErFgSgAAAAAoFScI7qUqs6RzMrZ2ew8yP/ve6diffcvLtyW+wenxmP9w9GsXtVrVePyXWplvuD4Zs4jrX4eedY3Pzsb6zl37igALEbqXNDJycl6wzvmXEt/OQbOEQWwGHNzc7H27ykbNmyIdX9/f9Uy58+fj/X69fTnQHF8IgoAAAAAKBUHogAAAACAUhHNLZtLkKbipwN92d8HvvbCyVhfuzkd6/ceylpmnx3PYmNPHTkR68mpmVj3VSdzq+NbyhGp9cumxvt8sXLEcZVatmaZRuMAAACQS+p3O3/5lomJibr1wYMHq5a5fPlyrDdv3tyqKWIFWfATUTPbY2bfMrNXzOxlM/uNyv2bzewbZna08v+m9k8XAAAAALDc5Ynmzkr6ByGEeyW9R9LHzew+SZ+U9HQI4ZCkpyu3AQAAAABoaMFobgjhrKSzlfqqmb0iaZekRyQ9VBn2pKRnJH2iLbPsUnliqj0uUvtvXzod66//ZdZN1z9Ov1vAJXxvi2Lki9cm5poaXzSCm+qM22h+OeYKAM1art1nl+u8AbTPm2++GevBwcG6Y3zXcF/PuisUvPDCC1XL3HnnnbG+ePFirEdGRhY/WawohZoVmdl+SQ9KelbS9spB6tsHq9sSyzxuZkfM7MjY2FhzswUAAAAALHu5D0TNbJ2kP5b0myGEK3mXCyE8EUI4HEI4PDw8vJg5AgAAAAC6SK6uuWbWr1sHoX8QQviTyt3nzWwkhHDWzEYkXWjXJLuKuaxtKkKViE35brohmKurWvHWfczbormp9eVZvuD41LKN57fwXKueNwC00HKNuC7XeQNon3379i045p577ql7/44dO3KtY/369YXmBEj5uuaapM9JeiWE8E/dl74q6dFK/aikp1o/PQAAAABAt8nziej7JP0dSS+Z2YuV+35b0qclfcnMHpN0UtJH2jNFAAAAAEA3ydM19z9KssSXH27tdLpfb29vrFe7zmX9vXkirj6C6x60aDxWkuZT3Wez+sb1G7GenZlJrLoN8d2cy/vnAAAAAGD5KNQ1FwAAAACAZnEgCgAAAAAoVa6uuWidn3jwr8T694ZPx7rXUunnVskbY83m8U++/O1YvzC2KtbzM9mFjot2xM0V5W24jBvjnlOYp4PuSnT16tVY+wtw9/Vlb21r1qyJ9cDAQDkTWyJ+G3T7c22Hnp76f5vN02XW73NLia65y9f169fr1v77tXr16ljTpbS41DaWqrezf/8cGhpq/8SAFYpPRAEAAAAApeJAFAAAAABQqs7IEnWgS2+9FetTZ87H+ubNm9mgRcSbrCeLvq7q71/c5Fou1C2np2fdGPc3i2Y64ubphnvb8v4L9TsH+1jcpUtjsT516lisey3H7u5W1lMTtRvZsSvWm7duX/ixkMvExETd+saNrGvzjOvafOedd8baRykHXRfqlAsXLlTdHh8fj7WP8O7du3fBxyrbpUuXYv2We3/y+765iP/+/ftLmVcR/nV98eLFWPtIsee/v9PT07FezHPz32u/b/k5rVqVnYJQG9urx49Zu3ZtrE+fzk67aBSD9c/P73+bNm1acN2dwm/X8+ezn5Xr1q2L9bVr12I9706j8M/fxx9rY87btm1rzWSX0BtvvBFrvy9v3Lixbp16P/Od9/32rn1v89tzx44dse62yP7Jkydj7d//xsay3wN27cp+dufZxlL1fup/7zt16lTd8X4deX4WAbiFT0QBAAAAAKXiQBQAAAAAUKoVH8098sPjsf7ai1mc6tRs1o3u4uUsxjWbO1o6X/9rfpn5PLHWHI9TcHzj9WX13Ez2d4owd9ONqf+weeK4ueK7DVaSetzBwSzWdunMkVg/8I5nsvGTWRxqUXqzSM/pN7Ko0+qhX4j1li07hPpOnDgRax972rlzZ6w3b95cd9m5ublY+/huUT7WJlXHqfzXfIxu69atsW53l8rJycmq2z5ul4rX+S6ar7zyShtn1zz/+vXbtfZ5v83HXZ999tlYLyaa6+Ohu3fvrrtuH8Hz+1yKH+OfT+1+lpKK8zYTzc3TNbfRMp6PxI+OjsbaR6n998JvV8/HIfOo7V7sX4/+sToxwuy3k//++tev32/8Pnf58uVC6/L72YEDB5Jfu3LlSqz9+6c/zaET+desj/L797zh4eFY++3q941mtnEtHxP329jH+v3Puk48RQLoJHwiCgAAAAAoFQeiAAAAAIBSrZhoro8S/YuvvxTrZ0azSNL0TBazmJvKulIWjZzetoxv7zpfPzqbGu/ju8oxJlR1lU3M4bZlqtrSJsao7viqx01tGyXqRtHhxLjqdSfmNJ91zquK407OqCkhi0nu2prtH9NTWaRu9OTfjPXuvZ0de2oH30HTd3mVqiO4vmtkKirl41c+AuU7EvouiV5VxDwRg5Oqo1+ej3v5eKLvUrl9e+s7J9d2at2zZ0+sfTTN8881tT06kY/Epr4P/rnNzs7WHZOXj4On1rGYWGsrlm2lPPOo3ZY+cu5jhb7zrd/f/evXv+Z9t1v/+vWKxp+l6tejX3e7X495HT+eneLjO277bZM6pcBvp6Idbf32S72GpOr3T/8+fPXq1Vj7aGnZnXX999THsA8ePBjrLVu2xNrvQ/45+M7Tnt/Gqa7NUr6fG77jtue3/8jISKxfein7ffPHf/zH6y4LrGR8IgoAAAAAKBUHogAAAACAUnV1NNfHj/7Pb2Zxj784mUUzZm5k8bBk5DQRP01FTm+77eO4qZhqKuKaGj+fXnfd8Q2ir0rFYAvOtXgEt8H2axArLk2ojfK529eyfWtVfxb7Hl77lVhPT/+9bMyq7rqIuHfu3LlY+663tTFRH01LRW37+/tjfebMmVj7/cFHq3yk0K/PR1x9zKy2K6+Pdfkolr/fL++jhz6ulYohNqs2SozF89HDV199NdZ+f/KdNv33OhXr9dE+f6F7f3+jmK5/HRTtLNuM2i7FGzZsiLWfr48h+uih39+Hhobqjnn99ddj7bef71TtX79+2dr3jtTr0W9nH4Mto5uuf3779u2L9djYWN3x/nvtI88+iurn7U8n8qcH+P3Sf698V16p+nvsfxfy78N+W6Zi0u2K6fqfG/5930eH/ekIqe3nn48/HcQ/Z7+f+feB2teB3wZ++/uord9P/Xi/z771Vnbqzl133RVrH+H2sWNgJeMTUQAAAABAqTgQBQAAAACUqqujuX/0F1l05j8ez+IUsz6OmydymoqfNuj6moysJrrPJteXYx55lm04v6LLJ8er7v2F19Vgfe3h11W/6/At8/XHTWXjBoayDrpnRp+O9c6DH252kh3FR498N8NUZ1epOlrlo3qnT2fbzHeJveOOO5qeZ61GFzP380tFtHwk7OzZs7GuvZh8q+TpgLqUHVqbUfZz8/E6X3s+Vunjp6mOpL5Lp5/r/fffv+h5NivPdm0UffXWrl0bax9j9DFYP8Z3rk11sfUdUn2c2Ud2/ePXSr0eL1zIupq3I5rr36ek6u64/v3Q8zHQVIfbVu0rPuoqVUd4/Xub3/5+W/q5+lMCWhXNPXnyZNVt//32kVq/bfy8/b7sX6c+Ft0uftv6n11+G6eizX7f2L17d6yvXLkSax9vB1YaPhEFAAAAAJSKA1EAAAAAQKm6Lpo7eia7wPWfHc3iJdVxXKeZCGmDCFQqstqq9RWN496efF3882vV/BpFh9seN6zqiJuK49Z0yqyab2L5a1k8Z13fs7Gemno41mVfLLwdfEfC1AW+a/moo4+13nvvva2b2AJqO5P6SJiPWfoIr48x+gjjtm3bYn3ixIlY79+/vxVTxRLw8bo870F+jH9NdLra55bq3uvjlD6K2sx7mI+A+vi97x7r4/5SuqOufz36KLCP6frXaVH+fWB4eLjqa/40BD+nVMTfP+92dEiu7Zrr47WpuGvqFATfxfbNN9+MdTMx2Nroqo96p2Kt/meG395ld5z129bHdH0s3XfvrY2+v81/T3xkl2guVjI+EQUAAAAAlGrBA1EzW21m3zWzvzSzl83sH1fuP2Bmz5rZUTP7opmtWuixAAAAAADI84nolKSfDSG8Q9IDkj5oZu+R9DuSPhNCOCRpXNJj7ZsmAAAAAKBbLHiOaLh1MsnbJxj0V/4FST8r6Vcq9z8p6VOSPtv6KRbz9ReOx/rKzew4O3VeaOqSI02dm9lg+VyXZmniXND0c6s+37Ed53a26jzSvNKXfKi65crUeaF5zxHNcV7pbHb/0FB2ntLYpaOxHt75Y/WmvWzluVxE7ddS59CUzZ+LNjo6GuvBwcFY+/N6PH9OVWpMs/ylAvJc4qQbLuWS5/4yLKdL5+S+LFaFP09Qqr4Mhd+XDx061KopLshfAqn2MlCp8xo9f3/e89YX4s/j86/FWn47+/MG/WVomjlXdTFS51f681P9JUQ8f06pv9TR9evXY+3fI/No3Edj4deXv+zMUvLni7722mux9ucQp34e+OeT2o+BlSbXOaJm1mtmL0q6IOkbkn4k6XII4e3uDKOSdiWWfdzMjpjZEX/tJwAAAADAypTrQDSEMBdCeEDSbknvllSvzWXdP2mFEJ4IIRwOIRyu7ToHAAAAAFh5Cl2+JYRw2cyekfQeSRvNrK/yqehuSWfaML+884r1C+dmYj03vUSR04bLz9cb0uSlWQrOdRHLF4/vFlvXYvT3Z7uvj5bNzGRt4atjtz4Kk4jWVsVva7dZjmiuHzOX7Ys3xl/M7iea2zGRRm/dunWx9vtTKubn79+9e3esffKj2T++Lad4aDPyPLdG0ch2WK7bPs+8fWxTqr4Mit+Xl4qPtErVsVYfG/X8c/XLF73kiL+kk3/9pmKsUnV0uNO2pVQ9j2PHjsV68+bNsc7zPuejyj5KnUft67fo66tTTufw/M8Jvw/keT7L6XJPQDvl6Zo7bGYbK/UaST8n6RVJ35L0y5Vhj0p6ql2TBAAAAAB0jzyfiI5IetLMenXrwPVLIYQ/NbMfSvpDM/ufJb0g6XNtnCcAAAAAoEvk6Zr7fUkP1rn/uG6dL7rkLk9kkZmz17MPeefnsq6Wqe647egY23j5qkGF1tGqOO3tXyqvI26zUTaz7Pu7ZbW5+11sZ+5cVs+6OG6errnJLrt5x7n6RrbuTUMnhc7mu0meO5ftQz097j0lEQ/1MSsfX+O8eHSi2vfhVNx1qaxZs6bqdtFYdjNdrH2sd+vWrbmW8VFnH83tRP5779/bUnynXP88/fckz+N0o2Z+t1mp2wyoxSsBAAAAAFAqDkQBAAAAAKUq1DW3U427COT41azD20BfdpydL3KqBcfnjZw2tXwzHXHzRIIlBRWcn3KsLzGmlZ0le1YNxPqh+7ILS/t17N13OVvgchbPLtwptzaam4rzVj0/d/9sVq/fkHVivHolm9/6oSwOulwtpmtup/Odb/fs2RPrPBHGVsYci3aWXE7buNO70nb6/FIWM+9O60i6YcOGqts+Btru18Hg4GCsfcfYRo/p5+c7qXaivXv3xnp6ejrWqec3N5f9frV+/fpYnzmTXShhMd2Bl+vry/Pf66L7aNldwIFOxSeiAAAAAIBScSAKAAAAAChVZ2dIcrrorsNcHYkoGDlNjU+MScZVa9adjKwqEd/INX7xc5Uk+RhJnmUS27Vq/PzCz20xfKfcg1uzbor/xQcOx3r0+HdivWfrqJtqE51yb4vm5uiUm4r8TmVRz7NXTsWaaG5n8hd6zxO58p0lt23bVnf8YuKP3RBfSykaZSvbct32y3XeXm9vb9VtHw/N85yKvtZ8p9t169bFemJiItfy/f39sfadZTvRwEB2esvrr78e65GRkVj7SLLnvw9TU1N1x+TVDfup123PBygLn4gCAAAAAErFgSgAAAAAoFRdEc29cHky1snOsnkip0U7wKZiwDnnkas7btF4cY5uvbcvU3/ehbv3JqO8xfk47rZtw7H+H3/5gVjPz2cxoS1rv5ktfMNFhkI2pqlobaNxubrpzsRqTc+40Nl859s1a9Y0GHmLj6xt3JjFrf3j+I6TADrHxYsXYz00NJRrGR8fvnLlSqx37drVuom1mY/H18ahiywLAIvFJ6IAAAAAgFJxIAoAAAAAKFVXRHPHr9e/KHOeSG16vOren2fZRsvkWj413nWllesK2NO3yg2vHw2dr72AecHYrY+c5lm2KB/FlaR37tsU6//hb7wr1of27Yz19Fufi/Vgb3ZxbU35OG6e2GyObrq3LZ+j624iytunMS13qY6nebvmdnr3QN91s+i8Z2ayGPbNmzdj3Ww0t9s6MRbtMlnG8yy6X3fKtm923p3yPFKKPqeisdHJyez0Hh/NbbRdfHfc73//+7FeTtFcf9pB0S7WvlP4YnTi678Zy+n9AugkfCIKAAAAACgVB6IAAAAAgFJ1RTT35vRsdqOJSGwzcdrbIhdNRF/93wd6+rJv0baNg7EeXuO63U1mHf8WEylMdxGuW1bNNU/QxD/kQH/2fEaGN8f6oQcOVS3z3gfvjvX161lHwumLT8R6VfjLbIGpLA7ZTGw2OaZ2XDLy2yDaW3Fj8q2696NzrFrl4+7F4lQ+4tbsRd8BtJ/vep2Xf1/o61uev0r19Cz+s4jFbDMAqMUnogAAAACAUnEgCgAAAAAo1fLMk9SYc91kfcw0TyS2mW66yRhwznX4GGdPf9aB7/4da2P9S+/YGuv3u7jq4GDW7a4bXJm4WHV77NSfxXp447PZF2bOZfWsjwY1EZvNE9ltdnl399xMc90GO81iuuYuJ810bpydnV140CLmkbq/G7ZxnvvLsFy7YC7XeefVjufUTJdsqTtiqkX3G7/Nylhfp+u2LsBAWfhEFAAAAABQKg5EAQAAAACl6opobp5YTXMdcd3KQr4Lh6cjv9nyq9ZkF7j/6I9lHXH/7t94SPWcfOM7sZ669kqsZ25kXWVDrj62bVK16ux5+gTP/FzW3fbQvuzvIENDl6sfa3AiqyemEytpUew2OWYRXXNT667eidRNujGa66N2ReNUvmtuM10p8657OW3XlDzPzW/XMizXbb9c591Iu5/TwMBAoXVJ0vXr2SkWBw4cWPS6l9LMTPbzuOg27u/vb2rd3bifvq3bng/QTrl/SzKzXjN7wcz+tHL7gJk9a2ZHzeyLZrZqoccAAAAAAKDIn+t/Q9Ir7vbvSPpMCOGQpHFJj7VyYgAAAACA7pQrmmtmuyX9oqT/RdLft1tZ2J+V9CuVIU9K+pSkz7Zhjgta0++Op5vqiLv4+O7tHVbrj+tzEaD/7qe2xPojP/eTsT598vlYDw9+M9Z7d7j46ly5MbWk4LoFFo2u3nCR28u1XQcLRmebis3mifXmXX7hOG5f32p1k7xxo+UUS5qamor1+vXrG4y8nT9VYPXq5r7XebbZctquXqc/t06fX8pynXcjeWLZzTynoaGhwo/jO2Jv3Lgx1j7W39vbu+g5leHGjRuxLrrfDA4ONhi5sG7bT7vt+QBlyfuJ6P8m6R8q+y17i6TLIYS334lHJe1q8dwAAAAAAF1owQNRM/vPJF0IITzv764ztO6feszscTM7YmZHxsbGFjlNAAAAAEC3yBPNfZ+kv2ZmH5a0WtKQbn1CutHM+iqfiu6WdKbewiGEJyQ9IUmHDx9uSy5hw2DWvS0Zx1WizjM+T3e3muNwP66nL5vfh+7aEOuqOO6J/xDrXXuezh5ofDKrryxlR9w88dMWxWabXb5wbDbPvJtcd8j+5rN6cKu6STd2zfWR2qJdc32nXB/5W4xu7izp45ad+NyW67ZfrvNupN3PafPmzbH23XAXE9MdHx+P9datnf1e79+rUp3CPX+/P31hMbptP+225wOUZcFPREMIvxVC2B1C2C/po5K+GUL4W5K+JemXK8MelfRU22YJAAAAAOgazVzk7hO61bjomG6dM/q51kwJAAAAANDNcnXNfVsI4RlJz1Tq45Le3fopAQAAAAC6WaED0U61bePaWFfl9HOcC6oc547Wnv9Zd0xN9t/fHt64LtZ/9xffFeuzo8djvWvPv8sWvnit7vrKkbp8Seo8yDznbybGp8bkfdy2X7KlwZwKn3vq9pX+nVruFnMOzHI6P8af8+WlnoO/TMPVq1djvX///qbmUfS8o+W0jTv9nKrluu07fbsuRtHnlOdyL96Au6zaK69kl0vftm1brG/evJlc3p9fOTExEetOPEfUb5t9+/bF2p/bmuLf54aHhxe9Xqk799O3ddvzAdqpmWguAAAAAAAVii0ZAAAbzklEQVSFcSAKAAAAAChVV0Rzh1e74+lU3CNxfzIekmvZ9KU+/CVb3r0z28ybNmSXczh97gvZAteWMI5bNKbaVEQ1x5jlPKfUY/Wsc/Xyj+Z2u8nJ7LJJ69atazDyFh/tO3v2bKybjeYWNT09Xer6kCF21x38JZf85U3ySsX6O8XFixdjnSfC7OO4MzMzsR4ZGWntxACsSHwiCgAAAAAoFQeiAAAAAIBSdUU0d9Ngdjy9ad3qWF+/MRXrZEfcFsV3a1NZZhbr9x7aHutrV6/Eetfet7IFrrQ7zlMzwU7rSlu7AdsRr21m3ouZU092/5XrWYxp+x171E3ydg7t9Oii73a7d+/eWF++fHnBZX2Eb8uWLYueg3/fkKq7caa2n+/muWPHjkWveyl1SpdJ/30s2tVzzZo17ZtYAd3cjVRq/3PavHlzrH3UvdHj+9fg4OBgrP3r10dcl9I1dxqQn2vq+fkxly5dinWz0dxu20+77fkAZeETUQAAAABAqTgQBQAAAACUqiuiuZs3bYz1jsEsGnl8Kou5hblUh9VmOuiq/nhVH+GPbMoiWxfOj8Z63Y4raq9UJLb2dgd0pa2dXzsiv8kxeR5TOefuxqzO9r/LEw/EeqhDIlpLqVNiap6P4NZGZOvxMc4bN27Eupl47Pr166tu+7iwf4/x8/NdOn3HzwsXLsR627Zti55TK/ltlmcbl81HLH0n5DyaiWSjc2zcmP0+8eqrr8a6tnv21FR26k/qNXj8+PFYHzp0qKXzLCL13uHft7zUe9u+ffvaMDsAKxmfiAIAAAAASsWBKAAAAACgVF0RzfVRkwdHsjjVG1fdcfbs9VgW74hbML4rqcelznp7s808fWNGbVVGV9q2x2YXs+52z7vB46aeh63Kyg0/rW6Vt2uu7yA5OTnZ1jnlNTY2FuvVq7OO2xMTEwsu67tJnj9/PtbNxNdq46A+Guijnz5C6vnt2ildXD0ff56Zyd4L83SZLKPjpN9+q1Zlr9888/Mx6rK7F+f6GZUYv9yU+ZyGh4dj7eOqUr7XYF9f9rPfR1zLfm0ePXo01j7+n3of9u9t/jXbynl3237abc8HKAufiAIAAAAASsWBKAAAAACgVF0RzfV+5h1ZLO5Pf/SjWF8v3BG3aHw33WE1+LhnyxpFNhNRbbR8J3TTbTAu15xaFdldxJx6svrG/P2x3rXvbnWrvPFJH03zXVxfe+21WG/atKnumFbxnWQl6fr1LLI/P599T1PPw3f79V0z77333lZNscqGDRti3d/fH+tUt0vfvdPP7/XXX4/1nj17Yt2uiOC1a9di/eabb8Z6+/btdceklB1l89vGRy9T8/BjfJR3KXVjRHCpnpOPwx87dqzqaz6+6iOu/jXo9wkf9/fvI+3ab0ZHsw79mzdvjvX4+Hjd8al53H13e352ddt+2m3PBygLn4gCAAAAAErFgSgAAAAAoFRdF82999AdsX7/nVkE5et/6brD+Qie6scpqiMUeeK71fNoSwSjcPQ1bwfYotHZZrrS5njM3Otu1byb3WauXpe9pC5O/VKs9/TwNx/PxzJ9PNTHY1966aVY+9iYr/3F4/1r7sqVK3XXW3tRet/J19e+E7d/XB+V9fFYH5ttJd+B18eK/fqmp6dj7efto7lr166NtY/s+efmO3z68f5+qfp5p6LAPmLtu3S+9dZbdefaKfz+5Dsh54lh+m125MiRWN93331V6/CPheXjzjvvrLrto7r+fcF3mfWvCf9aO3fuXN117Ny5M9a1r7t6/D5aG3X376v+def591IfF/avWQBoJ347BgAAAACUigNRAAAAAECpui6a6/3t9+yO9Utns+6Yp89kcZaqaK7qxy3DfP0x6WWbjeYmuvq2Kn5au0zLIr9FY7OJWG/ex21Zl97FzMnpz/6eM2M/Getdd/zV+uO7TN6uud7q1avr3u87S/q4m7+YvI90+liqX7ePnPm4r4/NNeLX5+fhu6Tu3r1b7eaf68aNG2Ptn5OP6aYuUO/v9xG8gYGBuuN95NSvq3ZOPj7ot5nvzOnH+Fiqjy36WLS3mH2rVUZGRmLto41+3/L7n9/GvuNzbWTy6NGjsfavg+Hh4Vj7LqdFdWP3zqLPqXafbYc77shOA/IxXR9r9V25/f7uXxP+Nej3sxMnTtR9nL1798bax29rTw9IdaX2EXLPz7uM+Hi37adFn89yem5AO+U6EDWzE5KuSpqTNBtCOGxmmyV9UdJ+SSck/c0QQv2+4AAAAAAAVBSJ5v5MCOGBEMLhyu1PSno6hHBI0tOV2wAAAAAANNRMNPcRSQ9V6iclPSPpE03Op6X27toR63/4gSwW87/+eRbbOX8x+xB3btp1uEtFcKuioakI7SJiF3m6z+aKqOZYtuHyBTvitiM2m3tcni69LdqWtQazeON0yDqbTqz69VgP0ym3io90pmJj/gLyMzMzsfaR2FR32zxqX5c+puoja37dPua3a9euQutrJR/h27p1a6x9hM93ffVRwNT285G/vPw29NFSvy39PHxkNRUR9rHeVKzSjy+D77DsI7h+W/rviR/jOwv7WqqOWPs4pI+lNxPNzROrrB1z9erVRa+vHWr3Ab+f+di4559TGfuKf+85dOhQrF999dVY+/0jz2vC19u3b4+1fz5+//PdwWufs98eqfj+nj176o4pQ9H9tNP2Uan6vdDHpFO//6W6bwMrWd7flIOkr5vZ82b2eOW+7SGEs5JU+X9bcmkAAAAAACryfiL6vhDCGTPbJukbZvbqgktUVA5cH5eqT7IHAAAAAKxMuQ5EQwhnKv9fMLOvSHq3pPNmNhJCOGtmI5IuJJZ9QtITknT48OElaxP24N1ZfPJ312Uxn//7my/F+nvns9jEhfEsNpHsdOZinPM1UYxZy2Iu83OJDn5NRVybiKjmXXePqwd8BNJ9kF7VaHi+iTG18/O33bqDn0dI3O/HFxzTKCTQl8XlTo2/M9aDO/+bWA9vHdFKkHpNNIqk++iSj7X5C8VfuJC9jfj4ru9U6iO0qZheSm38zEdwfUddH7XLc2H5svnn4efqY7o+Euqjnn6bpWLOvu6piZj7dZ88eTLWvruwjwn676/3ne98J9Y+8uzjwn5/qp1HmVJR2ePHj8fa7687dmSnhfh9TKqOnfrt7Pfxovz2fv3112Od6kZc+zpYysh5PbXz89vZx1H99vPboOwYt3fPPffE2sdJ33zzzVj714p/rv61mXpf9XFfv2xtzN5HP/3+e/DgwRzPovVqvydF99NO20el6vj+j370o1j7mL7nt0Ez8Xugmyz4k93M1prZ+rdrST8v6QeSvirp0cqwRyU91a5JAgAAAAC6R54/9W+X9JXKXx77JP3rEMKfm9lzkr5kZo9JOinpI+2bJgAAAACgW1iZF9U9fPhwOHLkSGnrK+r8+YuxPjl6OtY3XTfdvMzFyN71V+6P9fFjP4j1/fu+lC0w6dbRjk65ebvSrs6iTqPjD8d6fOo+dRofcQvzMw1GxlGFHr+3rzoetn3H7lhv2767dviKcubMmVj7OFmjLqy+Y6Af56NseYyNjcX64sXsNVvbnbTeeoeHh6u+5rv0djMfD/Xx59Q2811KfRRakjZt2tSSOX3729+O9e7d2evJ7xv+e+djwO9///tbMod28a+JS5cuVX3Nx/Z83NhH9XyMHd3NR439vuK74/rf03yE278WV8p7GYDlwcyed5f8TOL6EgAAAACAUnEgCgAAAAAoFQeiAAAAAIBSdd51CZbQ9u1b69atNF91audcok6cv9nMJVtqzxFNnnuanSPat2oo1juH97vhicvRlMGtu7fPnyvD+TErhT/Ps/acT9TXiZdBqL2syUL8eaudzl/yZ/v27Us4E3S6oaGhujUArAR8IgoAAAAAKBUHogAAAACAUhHNXUqpSG3hS7MUHHPb11x9PRu3o/ePs/tn/t8c60itL0+MOG902N0/NBLLa1d/J9br1m8QyuMvLdDoclBlXioKnc9fhmJuLjs1IbU/LadoLgAAWBifiAIAAAAASsWBKAAAAACgVERzl5TvlFs07tpMN90G4/zjzs5m9cxkYt055pRr3jXzyxM37smiemEDsc+lQjQXefl9YP/+/bEeHx9fcDwdRQEA6C58IgoAAAAAKBUHogAAAACAUhHNXUpV8dU5/4X6Y5rqlFsbfW1Hl95ULPime5jZ+mNk1fOz/sS6XT3P31GA5eTSpUuxnvXR/wQfzc0zHgAALB/8Jg8AAAAAKBUHogAAAACAUhHNLZ2Pr/o4bqrjbBtis4tavuCc5l0cd/VdWb3hF+ove/lr1fO7/nJW2yqhc9E1F3lNTEzEure3N9apfWPVquy1T9dcAAC6C5+IAgAAAABKxYEoAAAAAKBURHOXUiqC21Q0dxFdc5PdbvOs290/P5XVAwezeu/vZnXvBtW14QPVt9/4eFbfPJrVxHSBZaunJ/vb540bN+qOGRgYiPXY2Fisd+3a1b6JAQCA0vGJKAAAAACgVByIAgAAAABKRTS3ZGaW3aiKwc7Vv79VsdnQKJqbI86bJ/Lro7lDD2d1Ko7r9W6svu276954Navd5lN/9ncU6+FvKkuFrrlo5Lvf/W6st27dGuu5ubl6w9Xf3193PAAA6C65fns3s41m9mUze9XMXjGz95rZZjP7hpkdrfy/qd2TBQAAAAAsf3k/RvrfJf15COEeSe+Q9IqkT0p6OoRwSNLTldsAAAAAADS0YDTXzIYk/bSk/0qSQgjTkqbN7BFJD1WGPSnpGUmfaMcku8mWrduyGz2uA2xwHSTbEZutjeY2FflNrK/qcWbUlFA/tqfeLJt75eq6WK/bsa7eaJSAaC4k6cKFC7G+efNmrH28dnJysu6ya9eurTvmvvvua+UUAQBAB8nziehBSWOS/pWZvWBm/9LM1kraHkI4K0mV/7c1ehAAAAAAAKR8B6J9kt4p6bMhhAclTapADNfMHjezI2Z2xF8TDgAAAACwMuXpmjsqaTSE8Gzl9pd160D0vJmNhBDOmtmIpAv1Fg4hPCHpCUk6fPjwis/mbd+xJ9ZTYwdjPdDzfDZoLhW7bVFsttG4ZiK/viPw5T/L6o2/mNX9iQ/OZ85X355wy5vbTQd7Y3nl2k/FeoiuuUBb1P4B8erVq7H2MdqhoaFYX7t2Ldap7rg+jrtqVXaawp49e+oNBwAAXWbB395DCOcknTKzuyt3PSzph5K+KunRyn2PSnqqLTMEAAAAAHSVvNcR/XVJf2BmqyQdl/SrunUQ+yUze0zSSUkfac8UAQAAAADdJNeBaAjhRUmH63zp4dZOp/v1uAjpNXsk1gNDr2eDLr3llmhVp9yaVHTRyG+uWLDbnabeyOoTv57VG35edU18vfr2Tbf86jWxnJ49EOvNe/96/cdCqeiau7ycO3cu1r7T7cTERKx9p1sfoZWk+fnsNe8jtX55b2BgoO54H9nduXNnrHuI2QMAsCLwEx8AAAAAUCoORAEAAAAApcp7jijaYMu2Q7E+f+pXY719y+ezQeMXs3o21R23YDddqUFH3aKR3wTL4niaOpHV5/95YoGaXXH9YCyn5/fHeqzvt2O9a92QsPQGB7Pv1c2bN2Od6paKpXXmzJlYDw8Pxzr1fRwfH08+Vm9v1sXaL+/vv3z5cqzXrMli9gcOZDF7AACw8vCJKAAAAACgVByIAgAAAABKRTS3Q2zf895Ynx7NOlYOD3wx1qvWPJ8tMO3isTcKRmul5jrlFmX99esB93eQgSzWJ0lnrmYNmdfu+q9jvWvDlsXPAy0zNJTFok+dOhVrH8n0MUxJWr16dayvXLnSxtmhkfXr18faf798V2Pf6baW73x7/fr1WM/MzMTax3QffPDBxU8WAAB0LT4RBQAAAACUigNRAAAAAECpOBAFAAAAAJSKc0Q70K7dh9ytfxSrM28+F+u5ia/Hes/2Y9nwqbNuWXde51zNOZ7J80pzXJrF67WsHsjxd43V2SUbzo5n547ND/xC1bBd991TbB4o1bp162J977331h1z/Pjxqtv+vFB/DiHKtXHjxliPjo7Guqcne/3680j9pVyk6vOD9+3b144pAgCAFYBPRAEAAAAApeJAFAAAAABQKqK5y8jOfT/hbmX15bcuxvrqtVdjfWPi+7Gev/py1WPds/dqdmP6clbPT7lRPqabXeZBfVksc2Iyu0THG6PbY71pJIvdzvbflz2HHVnkdmSk+vIe6C4HDx5c6imgjuHh4bo1AABAmfhEFAAAAABQKg5EAQAAAAClIprbBTZu3urq97uvvP/2wfWE2ayen6s/xtzfLCzrlLvBsl3ogXxrAwAAALDC8YkoAAAAAKBUHIgCAAAAAEpFNBeSi9eql10CAAAAQHvxiSgAAAAAoFQciAIAAAAASsWBKAAAAACgVAseiJrZ3Wb2ovt3xcx+08w2m9k3zOxo5f9NZUwYAAAAALC8LXggGkJ4LYTwQAjhAUnvknRd0lckfVLS0yGEQ5KertwGAAAAAKChoi1SH5b0oxDCm2b2iKSHKvc/KekZSZ9o3dQgSSGElowpm5m1ZAwAAACA7lP0HNGPSvpCpd4eQjgrSZX/t9VbwMweN7MjZnZkbGxs8TMFAAAAAHSF3AeiZrZK0l+T9EdFVhBCeCKEcDiEcHh4eLjo/AAAAAAAXaZINPdDkr4XQjhfuX3ezEZCCGfNbETShdZPr/OkYrC19/vb8/Pzde8vOqaZupHUuKLx2lbVktTT07PguDxjUnMFAAAAsHSKRHM/piyWK0lflfRopX5U0lOtmhQAAAAAoHvlOhA1s0FJH5D0J+7uT0v6gJkdrXzt062fHgAAAACg2+SK5oYQrkvaUnPfJd3qorus5YnKFq0bfW1ubq4l62gm1lvv9kKKxmtTsVl/f6pu9LXe3t5Cj1W0zhPxBQAAANCcol1zAQAAAABoCgeiAAAAAIBSFemau6wVjb6mIrT+/lSdd1zRumjENxXfrVW0a27R2G0qTuvv93Wjr7WqLhr3XehrAAAAAPLjt2kAAAAAQKk4EAUAAAAAlKqro7ntiN36enZ2tu79tV9Ljcszpui6F9PVN080t2j3WR937evrq3t/o2iuXya1fNEx/jn7Mf75p+5vhJguAAAAUAy/QQMAAAAASsWBKAAAAACgVF0dzW2HVFwzb4wz77h641P1Ysbneawy5pFH0fkVXRYAAABAufhEFAAAAABQKg5EAQAAAACl6upobp5upr4zbJ7aP6bvVlvb9dV3tU11cS3apTdP599UN9zarrn+a/755Rnjt0Gezrr++afur91+ebrrFu3Mm2fdqeew0NcAAAAA5Mdv0wAAAACAUnEgCgAAAAAoFQeiAAAAAIBSdfU5ol7q/L7UpUVS513mqaV853DmGVP0/M/FXL4lj2bOpc1zHmmj8zHznMOZZ0yeOvV8AAAAALQOn4gCAAAAAErFgSgAAAAAoFQrJpqbkida6jWKu3rNxGiLRm3zxnFTUuPyRFOLRnbzbu9mIr95HgcAAADA0uETUQAAAABAqTgQBQAAAACUasVHc4vKG/VMRUXzyBOpLdr1tuj4RopGXItGfAEAAAB0t1xHS2b235vZy2b2AzP7gpmtNrMDZvasmR01sy+a2ap2TxYAAAAAsPwteCBqZrsk/T1Jh0MIPyapV9JHJf2OpM+EEA5JGpf0WDsnCgAAAADoDnmjuX2S1pjZjKRBSWcl/aykX6l8/UlJn5L02VZPcCUiygoAAACgmy34iWgI4bSk35V0UrcOQCckPS/pcghhtjJsVNKuesub2eNmdsTMjoyNjbVm1gAAAACAZStPNHeTpEckHZC0U9JaSR+qM7RuN5wQwhMhhMMhhMPDw8PNzBUAAAAA0AXyNCv6OUlvhBDGQggzkv5E0k9J2mhmb0d7d0s606Y5AgAAAAC6SJ4D0ZOS3mNmg3brxMSHJf1Q0rck/XJlzKOSnmrPFAEAAAAA3STPOaLPSvqypO9JeqmyzBOSPiHp75vZMUlbJH2ujfMEAAAAAHQJC6HuqZ3tWZnZmKRJSRdLWylWqq1iP0P7sZ+hDOxnKAP7GcrAfrYy7AshLNgcqNQDUUkysyMhhMOlrhQrDvsZysB+hjKwn6EM7GcoA/sZvDzniAIAAAAA0DIciAIAAAAASrUUB6JPLME6sfKwn6EM7GcoA/sZysB+hjKwnyEq/RxRAAAAAMDKRjQXAAAAAFCqUg9EzeyDZvaamR0zs0+WuW50NzM7YWYvmdmLZnakct9mM/uGmR2t/L9pqeeJ5cXMPm9mF8zsB+6+uvuV3fJ/VN7fvm9m71y6mWM5SexnnzKz05X3tBfN7MPua79V2c9eM7NfWJpZYzkxsz1m9i0ze8XMXjaz36jcz/sZWqbBfsb7Geoq7UDUzHol/TNJH5J0n6SPmdl9Za0fK8LPhBAecG3BPynp6RDCIUlPV24DRfy+pA/W3Jfarz4k6VDl3+OSPlvSHLH8/b5u388k6TOV97QHQghfk6TKz82PSrq/ssw/r/x8BRqZlfQPQgj3SnqPpI9X9iXez9BKqf1M4v0MdZT5iei7JR0LIRwPIUxL+kNJj5S4fqw8j0h6slI/KemXlnAuWIZCCP9e0ls1d6f2q0ck/V/hlu9I2mhmI+XMFMtZYj9LeUTSH4YQpkIIb0g6pls/X4GkEMLZEML3KvVVSa9I2iXez9BCDfazFN7PVrgyD0R3STrlbo+q8c4JFBEkfd3Mnjezxyv3bQ8hnJVuvTlK2rZks0M3Se1XvMeh1X6tEov8vDu1gP0MTTGz/ZIelPSseD9Dm9TsZxLvZ6ijzANRq3MfLXvRKu8LIbxTt+JEHzezn17qCWHF4T0OrfRZSXdIekDSWUm/V7mf/QyLZmbrJP2xpN8MIVxpNLTOfexnyKXOfsb7Geoq80B0VNIed3u3pDMlrh9dLIRwpvL/BUlf0a1ox/m3o0SV/y8s3QzRRVL7Fe9xaJkQwvkQwlwIYV7Sv1AWV2M/w6KYWb9uHRz8QQjhTyp3836Glqq3n/F+hpQyD0Sfk3TIzA6Y2SrdOjn5qyWuH13KzNaa2fq3a0k/L+kHurV/PVoZ9qikp5Zmhugyqf3qq5L+y0q3yfdImng78gYUVXM+3l/Xrfc06dZ+9lEzGzCzA7rVTOa7Zc8Py4uZmaTPSXolhPBP3Zd4P0PLpPYz3s+Q0lfWikIIs2b2a5L+jaReSZ8PIbxc1vrR1bZL+sqt9z/1SfrXIYQ/N7PnJH3JzB6TdFLSR5ZwjliGzOwLkh6StNXMRiX9T5I+rfr71dckfVi3mi1cl/SrpU8Yy1JiP3vIzB7QrZjaCUn/rSSFEF42sy9J+qFudaj8eAhhbinmjWXlfZL+jqSXzOzFyn2/Ld7P0Fqp/exjvJ+hHguBKDYAAAAAoDxlRnMBAAAAAOBAFAAAAABQLg5EAQAAAACl4kAUAAAAAFAqDkQBAAAAAKXiQBQAAAAAUCoORAEAAAAApeJAFAAAAABQqv8frAYk2paRUysAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import requests\n",
    "import imageio\n",
    "from io import BytesIO\n",
    "\n",
    "response = requests.get(\"https://www.python.org/static/img/python-logo.png\")\n",
    "img = imageio.imread(BytesIO(response.content)).astype(np.double)\n",
    "img /= img.max()\n",
    "plt.imshow(img);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6IAAAEfCAYAAABbM3sFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XmwXNV99vtnGSSE5lmABJIAWUgIJIzADAZkYwhymFI2t7Apvwom5o0LEju4cu33rXJMKkmVU5W8vr6phATH3JCKHdvEZrCNzSCQGQxYs4VmMQjNE5qFBHLW/UOHpafbvc/Zfbp7n+4+30+Vi99pdfdevXvt3b3dz/7tEGMUAAAAAABF+UBPDwAAAAAA0LtwIAoAAAAAKBQHogAAAACAQnEgCgAAAAAoFAeiAAAAAIBCcSAKAAAAACgUB6IAAAAAgELVdCAaQrguhLA6hLAuhPDVeg0KAAAAANC+Qoyxew8M4QRJayRdI2mjpPmSPh1jXFG/4QEAAAAA2s2JNTz2YknrYoyvS1II4fuSbpKUeSAaQujeUS8AAAAAoBXsjDGO6upOtURzx0raYH9v7LgNAAAAANA7rc9zp1p+EQ0VbvudXzxDCHdKurOG5QAAAAAA2kgtB6IbJZ1uf4+TtLn8TjHG+yXdL5VGc2fNmlXDogEAAAAAzWDevHlVP6aWaO58SZNCCBNDCH0l3SrpsRqeDwAAAADQC3T7F9EY49EQwt2SnpB0gqQHYozL6zYyAAAAAEBbqiWaqxjj45Ier9NYAAAAAAC9QC3RXAAAAAAAqsaBKAAAAACgUByIAgAAAAAKxYEoAAAAAKBQHIgCAAAAAArFgSgAAAAAoFAciAIAAAAACsWBKAAAAACgUByIAgAAAAAKxYEoAAAAAKBQHIgCAAAAAArFgSgAAAAAoFAciAIAAAAACsWBKAAAAACgUByIAgAAAAAKdWJPDwAAABRn3759JX+ffvrpFe+3d+/eIoYDAOil+EUUAAAAAFAoDkQBAAAAAIUimgsAQItZsmRJqq+55ppU33777aleuHBhqteuXZvquXPnljzX9ddfX3EZzz33XM3jBAAgC7+IAgAAAAAKxYEoAAAAAKBQRHNbyPLly1N97rnnpvq2225L9dtvv53qdevWpfrIkSMlz/WBDxz//yBCCBWXF2Ps8j5Zz7lt27ZUDxkyJNX33HNPqn/yk5+ket68eV0+PwDgmAkTJqT66NGjqf72t7+d6h07dqS6T58+qT506FDJc23fvj3Vvt8HAKCRuvxFNITwQAhhewjhVbtteAjhqRDC2o7/DmvsMAEAAAAA7SJPNPffJF1XdttXJc2NMU6SNLfjbwAAAAAAuhTyxHBCCBMk/TTGOK3j79WSZsUYt4QQTpU0L8Y4OcfzpIXNmjWrm0NuT8uWLUv1tddem+qPfOQjqfao7ZYtW1LtMas8Edq898v7XF3xyK5HyN59991UHzx4sOQxHuFdunRpqp999tm6jAkAWtlLL72U6v/+7/9O9XvvvZfqgQMHpnr8+PGpXrFiRclzXXHFFRWX4ftuAAA6U3aa3cIY48yuHtPdT5kxMcYtktTx39HdfB4AAAAAQC/T8GZFIYQ7Jd3Z6OUAAAAAAFpDdw9Et4UQTrVo7vasO8YY75d0v1QazYX01ltvpdqjykOHDk31k08+mWq/OLl3pfUo1gknnJDqfv36pbo8YlVLZ8Ssx3qU17v0ehz3pJNOSvW4ceNSPX369JLn+ta3vpXq2bNnp9q7BXsXYQDoTU499dRUT5w4MdUe0/3tb3+bau+au2rVqpLnIoILAOgJ3f30eUzSnI56jqRH6zMcAAAAAEC7y3P5lv+U9JKkySGEjSGEOyR9Q9I1IYS1kq7p+BsAAAAAgC51Gc2NMX4645+urvNYegWP41588cWpXrx4carXr1+fao/dDht2/HKtY8eOTbVHcJ1HaMvjtB7fqlZWNDdP3NejYvv370/1T3/605L79e/fP9UrV65M9Re+8IVUe+R39erVXS4bANqF7299v+rdyJ2ftgEAQDPgxBAAAAAAQKE4EAUAAAAAFKrhl2+BtGTJklT/8R//caq98613mfWuucOHD0+1d6U9cOBAqvfs2ZPqE088sWKdV554bS338du9g+6YMWNK7nfw4MFUe0z3n//5n1P9l3/5l6kePfr4pWyff/75LscHAO0iz6kWtXRKB9C8/AoCM2bMSPXMmTNT7dF8/+61b9++VO/duzfV/v3qBz/4QaofeeSRkmWvW7cu1bfcckvF8fmpZ0A5fhEFAAAAABSKA1EAAAAAQKGI5hagb9++qV64cGGqPU51+PDhVI8cObLi83iEwu/v3WMPHTqUao+0ZnXWlRofx83inR7LDRgwINUe0/W48VNPPZVq7ygMAPXkp1fMnj071ZMmTUr1FVdckWo/PeC5555r8OjyqaVTOoDm5d9/Bg8enGq/MsGyZctS7d+9Tj/99FT798rzzjsv1R7Tveiii0qWPWLEiFT7VSF27dqV/wWgV+MXUQAAAABAoTgQBQAAAAAUimhuAbZv357qVatWpdrjER6d/cAHjv//A35xco+levzCo7ne9cwf26dPn5Ix5bm4eU9FdqXSdePR5t27d6f60UcfTfXnP//5VN98882pLu/wBgDOY7fXXXddqj/72c+mev78+alevXp1qh944IFUv/fee40aIgBk8u9FDz/8cMX7eIT27bffTvULL7zQ5fPfd999qfZTECRpwYIFqfbvlf699Mwzz+xyGei9+EUUAAAAAFAoDkQBAAAAAIUimluAEEKqvRusRxe8o6FHcE8++eSKz+lRXo+EeQz2yJEjVY81K0Zbbby21oun++N93XgnYI+BeHzXuwsDQGfGjx+fat9nfvvb3071s88+m+oDBw6k2rtVfve73031hAkTUu2nWhSh1n0vgNbi+6Gs74zexfaTn/xkqjdv3pxq/67lpyz46Qj+XUuSLrvssor38zF1doUEgF9EAQAAAACF4kAUAAAAAFAoorkF8+hDVoTKo7xDhw6teH/viOv3P+mkk1LtnXI9BtzZsmtRSwfdzh7r/+axZe9G7BeQ9/t4hBkAynkn85UrV6baT3kYMmRIqv1C736R+OnTp6c6z36+Vnk+PwC0v9NOO63L+2zcuDHVW7duTfX69etTPXHixFQPGjQo1R7H3bZtW8nzjho1KtUjR45MNV3EkRe/iAIAAAAACsWBKAAAAACgUByIAgAAAAAKxTmiPSjrHB+/hICfC+r5e7+EgJ+PNHz48FQfPnw41eWXNCk/Z7Qa9bqUS97zQp1fssXPhfLzGfIuAwBOOeWUVJ955pmp9ksOZF1eyzXLJQrY5wEo59+XvPbvVP490fuTvPPOO6nesWNHyfP6OaMXXXRRqv18+6x9JiDxiygAAAAAoGAciAIAAAAACsXv5QXziFfW7R7H3b9/f6q9hbZHUf2SLTt37ky1t8/OGxvrToy2Hvev9XmJowHoDt93+D6zlS4/4K/Bay7lAvRee/bsSbVHcP2SK2vWrEm1f/ccM2ZMxdojt5J05ZVXpnrw4MGp9svFeLQXKNflL6IhhNNDCM+GEFaGEJaHEL7YcfvwEMJTIYS1Hf8d1vjhAgAAAABaXZ5o7lFJX44xTpF0iaS7QghTJX1V0twY4yRJczv+BgAAAACgU11Gc2OMWyRt6aj3hxBWShor6SZJszru9qCkeZK+0pBRtimP437gA8f/PwGPR3i0wm/ftWtXqr0DrkcgvPtuEdHVPMsgWgugmWTFWom4Amhl3vnWr6jwwgsvpNq/M/opX6+99lqqt27dmuq9e/eWLOOxxx6ruOwzzjgj1f79FihX1TmiIYQJki6Q9IqkMR0HqYoxbgkhjM54zJ2S7qxtmAAAAACAdpH7QDSEMFDSjyR9Kca4L+//QxxjvF/S/R3Pwc9eAAAAANDL5ToQDSH00bGD0O/GGH/ccfO2EMKpHb+Gnippe6MG2U6yIgpZXW29c+Pu3bsr3qc7EbJqI7JZz1tLHDevrLic8/ER8wVQT3n2Qc2oVccNoHZ+atfBgwcr3se/Y+7YsSPVfsUG/37lp4JJpd14BwwYkGqPAns3XaBcnq65QdJ3JK2MMf4f+6fHJM3pqOdIerT+wwMAAAAAtJs8v4heLumzkpaFEJZ03Pa/JX1D0g9DCHdIekvSLY0ZIgAAAACgneTpmvuCpKy859X1HU57yrp48EknnZTqK664ItUe083T0dG77zqPUHh8QsqO2o4dOzbVO3fuTLV3UPPuvY2K7BYR+UVrWrFiRaqnTp2aar+wts/LG264IdUPP/xwyXO9+OKLjRhioV599dVU/+mf/mnF+zzzzDNFDaelZe3PsvaxeU6LaNR+iv1fe1i+fHmqp02blmr/TuCn9Fx//fWp/tGPfpTql156qVFDbHl51rFUug3feOONqf7xj3+c6lb6zPCuuRdeeGEPjgTIRk9lAAAAAEChOBAFAAAAABSqquuItrus+MZHP/rRVHuHMY9x5L1g77hx41LtXcz69OlT8XnzdLfNio35c5Y/T9Zj/HVUexH3ojvo+vgWL16c6gsuuCDVs2fPTnXWe+c66w73R3/0R6l+8MEHU/2LX/wi19hx3MqVK1P9sY99LNXTp09P9erVq1Pt0fIhQ4aket++fanOikzNnTs31R4rL//7tttuS/WmTZtSvWTJEvUUX0++H7r88stTPWbMmFS//vrrqS6/8HhP8fX3+c9/PtWnnnpqqv29O/HE4x9L3vXxoosuSvWWLVtS7acNlNuwYUOqfZ/uHR1PPvnkVB8+fDjVPrf69++faj91wrtB+ueHXxi+PBLn+5i333471Zdeemmq/fWtXbtWzWbVqlWpPv/881N96623pnrp0qWp9u0p63PT7/MHf/AHJcvbuHFjqn1f3+zWrVuXao+BTp48OdW+HezatSvVzz//fKp9PT399NOp9v2Xz2OpdC7/+Z//ear91IZWjez7Np+1L/TvV/69a/v24xd48HVc7oUXXkj1gQMHUu3bv++r7rrrrlT7OiYyDXSOX0QBAAAAAIXiQBQAAAAAUKhQZOe9EEJa2KxZswpbbrm33nor1VOmTEn1Oeeck2qPCXnEw2NFeaOrfr9FixalevPmzan2uEcj5B2rR3Y92utRNr89K+Lr6tlB1yNKHl+bMGFCqj0Kd8YZZ6TaL9CcN0rtPOrkz3XPPfek+rnnnkv1U089VfUy2onH0iTp6quPN9k+66yzUr1gwYJU+3rN4nFVvwB31lz0WO8pp5xS8m++/XsU0ztGe+TKl1evyNWyZctS/Yd/+Icl/zZp0qRUP/7446k+dOhQql9++eVUjx49uuIyfDsomkdzv/71r6fao+8/+9nPUu37F1/fZ599dqp9n+BxxnIeffV16ftb72Tu+wUfn3efdB6FztqnnHvuuSV/+/2yXp/zueg8nnj66aenujx+/j7v0l4+d6+66qpU+/vlUdvLLrss1fPnz0/1m2++mepa97Hv8yikVLoO7r777lR7zLInI5D+Xvz+7/9+qj2i7evVTxPx/Zm/zjz7M4/1+vcXqXQe7N69O9WXXHJJqv0zzb8X9aSsOLPvq32d+fbv69V5VDbPOpayPze8U7vvh7Zt25Zq3+b9Plu3bs1cHtAO5s2b538ujDHO7Oox/CIKAAAAACgUB6IAAAAAgEK1dTTXL/R+8803p9q7J3qcxzs3eoSivHtqLartRFs0nw8e5fJomneK8/tXO5fy3j9PNNdjNK6ekWfvSjh+/PhUe+TPY7o+t5588sm6jaPZeFfPj3/846kujxp6nO+VV15JtcfwPE7m3XG9HjhwYMU6q6PykSNHKi5LKo1TeSTU32t/3ttvv73i83o8tloeM/vc5z5X8m+//OUvU+3RPo+U9e3bN9XnnXdeqvv165fqRkf/O+NRHY+seRT6nXfeSbVvy9751+P3WXOgnEfyPCbpz+vL8+fyOef7PJfVZddP5yiP//k+zMfh8VpftnfmddVGc32eeBxeKo14eudq70L6xhtvpNo/Wz3+7J+VPv+8dr4uPG4+YMCAkvt512LfNn1/49tjEZ2GPRbs251va75f8Lni9xk7dmyq/XVnzTmPn/r8Ke+S7RFU/yz3Zfjc/7M/+7OKY21UZ13fn82ZMyfVfgrHf/zHf6Ta99UeB/d549vdsGHDUu3btX+X8X2+VLqteu3rw2PO/j3R17F//nziE59INTFdtDuiuQAAAACApseBKAAAAACgUG0XzfW4xxe+8IVUe/To/vvvT7XHLDza41EOjxVlra/O1mOezrJZql1eLcsqf7zXeSK4tXbHzXO/rGiux86yYpW1yorqeOTP6wcffDDVfnHsuXPn1m1MPSUrlubr26PuUmmkzKNLfrFxj6xlReInTpyYat+ufQ54l9gNGzakeunSpSXP5du8j2/9+vUVl+3xy3vvvTfVPh86u0h6JR6/9c6QUumc8/2Qx848Munrpsh9e2c8qvORj3wk1f7afLv2uJy/D95p2edGZ6/TY8/+eI+7egzWo9vejTNrLvr77u/PiBEjUl0e//O//b2bOfN4gsm3L59zzqO5Hu/MGqvHMH255f+WtV/1bdO3Wd/n+fvoEUiPDnuU0l+n8/UqlUYx/f329+trX/taxWXUq5tu+Vh9v+dRZf9O4evDo+X+nnoE1/dHPi+9e6zvO727/8aNG0vG5zFVj0/78vx937NnT6q/9KUvpdpj87XEdH1eSdKMGTNS7dFr/3zMiriPGjUq1T5Hfe57p1s/lca/H/h2I5V2ffbPDa99bvq4fR17TNcj0ln7v/J1A7QqorkAAAAAgKbHgSgAAAAAoFBtF831C4RnxW49tuPRPo+sZEVUs3QnEluvKGstnWtrHUcjXkNnsqK5eTqe5lled7oaZ0V1Zs+enWqPAC1cuLDqZTSDzZs3p/rDH/5wqhcvXpxqf0/KL2jv68DjXh6nuuuuu1Lt0UPvslv+vJX49uixvvL9js+VNWvWVHwuv9C7dyT1+KVHasu7V3blueeeS/Vll11W8m++brKirB7P8/hps/CoziWXXJLqrAip74drfW0el3MeM/XPBo/w+fr299T3Nd6d1Tt/+v7lV7/6Vcmy/b3zKGBW5+/hw4dXvL3aaG7WnCkfR9ay/XQV3w68I713AvWOxV5/8pOfTLXPh4ceeihzPL7O/TPe970e1/T3sdrt0WWdgiCVrkP/ruGfIb7O/Lk+9KEPpdpPCfDnzNrPZe3brrrqqpL7+b4tq4uwzyHn0VLvdJ3VkTmLfx5cfPHFJf/m+wXviOuv+8wzz0y1v+7zzz8/1T4H/DMqz+dEOV+Gj2ny5Mmp9nXj98nq9u1div19+PKXv1xxuX6qBtBqiOYCAAAAAJoeB6IAAAAAgEJVvmJyi1m1alWqvRuidzH06Ix38PNIhMe18kRFmyWOW6/l1uMx1Ty26K6e1cZx847Pozre5e/RRx9N9Re/+MVUe5TtkUceybWMZuDdCT1i5PEzj8eNHDmy5PH+mOuvvz7Vvv4feOCBVHus0mNtHj30x/r75bd7984nn3yyZEw33XRTqr2rpY/V9xcehfPOo/5as2LLeZTHKvPELJHNT8PwuKdvp85jdN4NNovHKv1zyOdueSzS56ZHHT3W78vOiuZWyz+vyruF+rbt25qvpzvuuCPVHld/+umnU+3dYL2jrX8uexdbj3F5/NHvX/68Hon3rsrLly9Pte+TLrroolTn2R7XrVuXao8R+7YvZce+fdn+XLfddluqfR3Pnz+/4u1++oLz6Pro0aNTXb5vu/HGG1Ptpw54h1/f5/nr8w66fvuFF16Y6jydXv2z0TtYS6Xz2qPX/pp82f6Z4dvNokWLUu1RWd/efZvzuvw7XNa6XbZsWaq9u7W/Bn/v/DPHx+Hb9csvv1xxWUBvw7cZAAAAAEChOBAFAAAAABSqLaK5HoXxiFxW7LbajriuO3HcPFopjtvoDrqdPaYR8dq8Y8pzP48YeQzHO3N6fKiV+IW5Pfbu8VGP1HlsUSqNU/lzeezJ417+XH4h8Gp5PK58TA8//HCqp02blmrvOOnv3WmnnZZqjyF7pNMju+hZ48aN6/I+Pj/8ffdunB6F9n2NR1e99v3A9OnTqxhxfrWc2lAe+d6wYUOqPQLp8dp//Md/TLXvwzzO6NvBOeeck+olS5ak2iOt/hp8HZfH+rN413uPQ/q+d+PGjbmeq9KyX3zxxcz7+XcKf4x3vr300ksrPpfv/zy+6rH+PLzrra97Sfr+97+f6ltuuSXVZ599dqq9U7jHdP07lUfLvRuxb1seL3YeVy3fFn078n2p832vrzN/rHer9ehwrfyzyD83vIO7r0tf/x6r9nnpkXZfr/4++KkgQG/AL6IAAAAAgEJ1eSAaQugXQvh1CGFpCGF5COEvO26fGEJ4JYSwNoTwgxBC366eCwAAAACAPL+IHpH0sRjjdEkzJF0XQrhE0t9K+maMcZKk3ZLu6OQ5AAAAAACQlOMc0XjsJI73TxTp0/G/KOljkj7TcfuDku6VdF/9h/i7/HwTSZo6dWqq/dyGBQsWpNrbbFd7bk2t54U2+pzKWi+D0uhl1PNc1WrPBc26vEctY+jsfr68HTt2pPqZZ56p+Nhhw4ZVNaae5Oe9+eUenG8rfu6YVNrG3s+h8fNszjjjjJrHWc6fv/xSGn4Okp8HNHHixFT7ZQf8fC4/Z8mf99prr031lVdemernnnuu6rFXq5b5XrSs7SZLM742H5Of49iMl9fxc+n8XDWp9DJGW7duTbVfmsX5+XB+bnWWGTNmVLzdzyH0+pRTTim5n39++7l7fq66X+rD91U+V7LORfTzhP1cSd+flZ/77edL+v7Mz4X08wmz9o2+vGqdddZZqV69enXJv/k6+8lPfpLq8847L9VZ58L7e+Hr0l9P1rp0nX1m+Dbi59v6usn6zPD31M8Trid/Xp9z/j76udb+We7ni/tr8+0m73nQQLvL9WkZQjghhLBE0nZJT0l6TdKeGOP7W+FGSWMzHntnCGFBCGFBpX8HAAAAAPQuuQ5EY4y/jTHOkDRO0sWSplS6W8Zj748xzowxzqz07wAAAACA3qWqy7fEGPeEEOZJukTS0BDCiR2/io6TtLkB46toyJAhJX979GHfvn2p9thJtRoVx+2p+3fn8Y2+lEve+3v8zd9rj8Ls2bMn1T4HPAqTJ/7X2TiqfU0eS/LIuL+GVuKv37ePvOvFH+OXxnCNjjROmjSp5G+Pl3kMbPHixan2y8j4ezdixIhUe+zR3/f169dXNb7ydVn0fqGnFHEKQrXaeV2WX8bIo+W+n/T4oMdA62X06NGp9kjwBRdcUHK/PJe/8ku2ZMVGs9aNX2LDt1/f3j3qWu7QoUOp9n2Yx279M8ovA1UvkydPLvnb31Nftr/3ftrCq6++mmrf5/nldvz0gs985jOpvvnmm1P9yCOPpDrrM6P876z3xe/j0dzOTgFphKyY7pQpx3+T8cuTebTZ+dzK+n7aqvsdoLvydM0dFUIY2lGfLOnjklZKelbSpzruNkfSo40aJAAAAACgfeT5v5JOlfRgCOEEHTtw/WGM8achhBWSvh9C+GtJiyV9p4HjBAAAAAC0iTxdc38j6YIKt7+uY+eLFs4jPFJ2VMc773ncIU8MpFrd6bDaiNhdd6KyPRX/y/s8HhXz+Mvu3btT7VEiv09WN756RgHzdOz12NM999yTao/nPP7447mW1wy6MweyOoz2JO98u3fv3lR7tM/j3T7n/L177bXXUr1s2bKKy/JumvXUSlGuZo8at9K6dHnG7dufx9Cl7C7YjYjjZvF4q8cfpdJ9ep7O6dW+j97h3GOsnXVn9e8Lu3btSvUVV1yRat/PlXe1bTQ/DcH3bR579lMkfKz+2jye7fFY3//568/SndMOmrErtc9FX6/lcfdKapmjQLtqjm+DAAAAAIBegwNRAAAAAEChGt9urAG8w52UHdvxTpbNEsct8v71jH40ustuOX8f/T06cOBAqmfOPH5FIO9g6jHJTZs21W1M1b6+rLnoF2FHz/KukZs3H2/8PXbs8csiZ0WufF5698+NGzem+u677061R3O9s2StfJ55TfQLWcrnhnel9X8rMgLpn+Plp9/4Npin+7nL2ievWLEi1b/3e7+X6qNHj6Z60aJFqS5fF/63x4r98c1yCoJfaeChhx5KtXca9g6//nr88+rgwYOpfuKJJ1Lt69Wfpx35PPXTwvx9z4N9NXBMc+wlAQAAAAC9BgeiAAAAAIBCtWQ01+MkUmmsqNFxh2aP49ZTkReZL49b+d8eB/Ko0/r161O9atWqVGd19qu2Y10932vvtvrwww9XfGzWBa6bXRHR8CJ4zOqDH/xgqj2KlXVR9f79+6fao2z+vvvFzPPKs25btRNjVmSyWV5Dq87r7ow7K97dU8rH4Pv9es13X8bAgQNT7TFg35bLT+Pxjux+asgzzzxT8T5+GkDR/DuS1/6Z4/snf63eUdnXR3f2Z65Vt68s1b6eVn2dQL3xiygAAAAAoFAciAIAAAAACtWS0dwTTywdtsdI8kS8qu2UW3SMM4+srrLdiYE0Q4y4/D3x+NC5556b6l//+tep9ouQO48S5YlxFdE11/lrraVrc9GyXn/eWGWeddOTcSWP0Q0YMCDVHkHbuXNnxcd6ZG3Pnj2pXrBgQap9v1V+esH72j2u1YwxtVrndU+pZUzlj22G7bGzfXUt+8ms1+afH/Pnz0+1dz/1bufl9u/fn2qPDnvtUf6eNG3atFS/8cYbqfZO3lmnwPh74fvFkSNHdrncWudPs39mNON+AWg1/CIKAAAAACgUB6IAAAAAgEK1ZDS3nMd2si4g3Yxx3GpjnB4Z8hjhSSedlGqPpeZV7WutVzTXX4N3KZWkX/3qV6n2C497d8MDBw6k2tefRyBriQLWM7Lr8zKrky961uHDh1O9ZcuWVGddqNzfO59zHucbPnx4l48FUKro7cO32bzfFfzzy08f8X29x3ebhe/PvEOwdwfP+h7l3y+yuryzbwNQDX4RBQAAAAAUigNRAAAAAECh2iKa6xrdhbSecVyPv3hQR1uLAAAgAElEQVRc5tChQ6n2yOnNN9+c6pUrV1asu9PBtNr71Ct649Fa75IrSWeeeWaXj/HYY1YH3XpGcOv1vFmxp1aKNPlYPaKW9do6e3yz8H2Hb4O+beZ9fe/LirRnRfE7G1O7yXNKRdFqndc9xcfdnfXajNujy7MdNPo1lHfQ9Si/d5xtVXnmSp4O/Vk664ScZ9mtNEfb4TMeKEpzf7oCAAAAANoOB6IAAAAAgEK1XTS3Fo3qsOr/5vEej5lOnDgx1fv27Uv12rVrU/3II4+ketSoUanu379/qssjrt0Zb3fvn+fC8P76vfZuuJK0devWVHvMxdeNdwuu11hrfUy18UuiOu0hz9yvVZ74Wt5ofitq59dWT80eI65V1uurZVvLE232qLYkDRgwINUe5ffHl3+utbp67s/aQbWnTvj6Y38GHNPen1gAAAAAgKbDgSgAAAAAoFAtGc2tZ4yhlufK+1iPYHgcd9KkSan+zW9+k+rOIqvvO/HE42/dDTfckGsc9ZInfponwvPCCy+ketOmTSX/5uvJL6Ld6DhurV2Ru/NcrSgrklRLV8VmkhX/a3QX286eP2vdZtXNqNp5U/Trabd53Q7j7ky95kp3tmv/DF62bFmq/+RP/iTVp512Wqr91JqelOfzO8+8qfZ0hO7s21pVO+yrgaLk/kU0hHBCCGFxCOGnHX9PDCG8EkJYG0L4QQihb1fPAQAAAABANdHcL0paaX//raRvxhgnSdot6Y56DgwAAAAA0J5yRXNDCOMk/b6kv5F0TziWvfiYpM903OVBSfdKuq8BY6y7WmKVeW/3eMq7776b6gkTJqR64cKFqT7llFNSPW3atFRv2bIl1Rs2bEj16tWrU71mzZqKY2pGeTroSlK/fv1SXW2cpREdgbuz7N7eYbDRMdZ68wj4rl27Uu3dMocOHdrl8+SJrNUay2q1dVuNrNfWLFG2Zln3WZ1e88SLW1me15qHrz/f9t3Ro0czl+WfUX4KjW/nvh9pFv5Z6132s7oF++379+9PtX9ncdu2batpfO0yT9/Xbq8HqLe8v4j+P5L+b0nvb1EjJO2JMb6/l94oaWydxwYAAAAAaENdHoiGEK6XtD3GuNBvrnDXiv93dQjhzhDCghDCgm6OEQAAAADQRvJEcy+XdGMI4ROS+kkarGO/kA4NIZzY8avoOEmbKz04xni/pPslKYRQl2xVeeTRow+1XDC41jhu1hi9A+zrr7+e6ptvvjnVixcvTvWiRYtSvXfv3lR7fMijQM0SWat2HeeN3BbZEbc76zLP/PDao1F+e/kF05tNs8yzRvH3xeP0We+Lx9f8/h5Z8wi9P8+UKVNqGyzqpl7zOiuK3+7bTTPK89k/evToVC9dujTVAwYMSPXZZ5+d6oMHD5Y83rvm7ty5M9U///nPU+3fTUaOHJlr7I2wZMmSVN96662pfu+991L9y1/+MtVDhgxJte/nPKqc1dHfo7nM/Wz1PFUDaGVd/iIaY/xfMcZxMcYJkm6V9EyM8TZJz0r6VMfd5kh6tGGjBAAAAAC0jWq65pb7io41LlqnY+eMfqc+QwIAAAAAtLNcXXPfF2OcJ2leR/26pIvrPyQAAAAAQDur6kC0Wfj5WFLpuRp+PoOf/5ClUecHeVv0iRMnpnrTpk2pfuKJJ1I9bNiwVO/bty/Vfl5otZcx6Uy9LjFR9LmZjbg0Sz3PC81z/6xL05Sfg4RiDRo0qGKd59xdPw8861IQhw4dqmF02fJcMohzkBqD9dqc8rwvvl37uY/Ov1uUP6ef/9m/f/9UHz58OPc4i+Lnvfr5rL4O/LU6/741fPjwVD/55JOp9nUxadKk2gbbS7DvAI6pJZoLAAAAAEDVOBAFAAAAABSqJaO5HoOTpNNPPz3VHk3bsWNHxcc34lIu5TzOsmrVqlR75M/ju5s3H7/6zcknn1z18iopInJay/PUM47bLJHdrMdnxSeb/ZItrh2jRGvWrEn1Rz7ykVR7rH/58uWp9ui/v4++Lc+cOTPVfkkYj8RlKV/HWfPGb/dYoMe7PS6XdamFojXjHPL16uvMa3/fnd/+zjvvVLyPR/GRTxGfXdOmTUu1X77l4ouPt77wz/Hy2Lsvz6Ovfom2z372s6mePn16qh955JGqxlor/y7kl2nxOe7fo9zu3btTfdJJJ6Xa577fB9macf8H9DR+EQUAAAAAFIoDUQAAAABAoVoymrt3796Sv6dOnZpqj7V6FM4jKFnqGZvwiNyECRNSPWXKlFRv2LAh1R4RLCKO2w7dcXsqspv3MX77yJEjU/3SSy+l2qO5Ht1qd83S3XXUqFGp3rhxY8VxZHWT9Pdu+/btqT7rrLNS7d04uxPNzepqmTW+devWpfr6669P9ZgxY1I9f/78LsdRT0OHDk31ypUrU+375HPOOafQMTmPYXuM1mOZWV1VPaqY9f6OGzeu1iGiwU477bRUL1y4MNVDhgxJtc9jqbQ7rneT7dOnT6p37dqVao//FmHJkiWpvvrqq1Ptc9w/izxq6/udPXv2pPruu+9OtW+//jwAUA1+EQUAAAAAFIoDUQAAAABAoVoymusxM0l64YUXUu3dIf1+HpX1CEq9OuiWRw09tuJd9Pbt25dqv8C9d9esVtb4Oovf1tIVuKdisz297FqW0SxR1HrJ6gicdR+pdI7v37+/4mMa3d117dq1JX97p0jfZt94442K93HeKfKCCy5ItcfUvFulx/+cx4O9w7ZU2hHXI7/eodXjoWeeeWaqt23bVnEcRfN93hlnnJHqok+XyOKfDaeeemqqPZr71ltvVXysz2l/r/K8tiK06v6lXKNfh5/u4xFr3zZPOeWUksd4DN63QT/l5qc//WmqvRvvrFmzUv3yyy93c9Sd8/H6/tZfq+9v/XuRR4o9Nv/ss8+m2k9HqDV+3qrzNM/nIIDO8YsoAAAAAKBQHIgCAAAAAArVktFcj3pJpR0NvQOid9D1+FUR0c2sLotZFz3PI08EN8/9u3O/Vu2aW+T9y3k875JLLql4u8/LdtBZHNw7UL722mup9m32L/7iL1K9ZcuWVC9evLjbY/II3ezZs0v+7ciRI6letmxZqkeMGFHxubybpMfXPNLp277H17KiuVnPL5XGPT1G5/Omb9++qfbOvx4R9NMA/vqv/zrVvo5rjQh6R9yPf/zjqf6Xf/mXVP/sZz9L9ebNm1Odtb/M2rfVk2+PBw4cSLVHc7Pu7/FH/7wpfx97ShHrr1nU8lrPPvvsVPu2v3Xr1lSXd8/2Lro+b/yUIO+gu3r16lTfcMMNqfbPhlq2QR+rVBrHffvttyvePnny5FT7ftjntZ8u4Pt0f221RnPbbZ5mvZ52e51APfCLKAAAAACgUByIAgAAAAAK1ZLR3BkzZpT87R0NPQrnkTyP22R1NKxnHDQrolhtNKNecdxmjOl25/6tFMcdPXp0ql955ZVUe2zvvPPOq2l5zcbXX/l25lEuj5l6tPSxxx5Ltcc1PRbnMbrx48enesCAAan2qKx3wP33f//3kjF5ZNXfl7Fjx6batzWPk3qkzsfn+53yfVUlU6ZMSXV5NM8vRO9dLb2zp4/bY4G+Pvw1PPTQQ6k+fPhwqn2d+Tr2SLVU2kXYY37eOfjNN99M9d/8zd9UHKu/vz7uzuLdjTBo0KBUP//886n2qK13LPY4t4/VY87+2v7+7/++ZHke0Zw3b143R53dsbPWz4BW0oiupb5t+Tbn75tUut1mdU/2uKzv5/7u7/4u1T63vvKVr6T6N7/5Tap9Xvr+5aqrrkr1Rz/60ZLxeQzeT2caPHiwKvEovz+vn1Lgpy/kieN29v608zytdtsEejN+EQUAAAAAFIoDUQAAAABAoVoymlvOO0h69zqP444cOTLVHqHy6E2WPBGSvLGTapdRy/2bMY7bnWhKM3TQzftcHkfNmltZ3ThbSXeigB799Eiob49ZF1v3qKjf36OyWZFOj9R5Z0ipNELp43M7duxI9fnnn5/qJUuWpNq7YdfSTbL8/t/73vdSPXPmzFRPnz491YsWLar4XN4p0yOCPj6Pnw4fPjzVK1asyBzjSy+9lGpfzx7b8/2w38fj6r4832485lgEnzf+Pk6YMCHV55xzTqr9ffd1OWzYsFSfdNJJqf6Hf/iHkuV5dNPjmldccUWqvftxeSS0knaL/9X6eVqLSZMmpdq3G98PSKXv98UXX5zqJ554ItXefdafy7c7f95vfvObqfY5NHHixFSfccYZFR/rkXupdJvy2verflrTJz/5SVXiEeFauuP25HvaKO223QE9gV9EAQAAAACF4kAUAAAAAFCotojmeve6UaNGpfr1119Ptcf5PBrlsd48Md288nR+zNMRt+jusz0Vg21UZLfWZeR5Lo89ZUWXPIqap5NqOyjv7Lxt27ZUe6dI77bqMUTv4jh37txUe5w2a1vxbdkjdB6Vk0ojch7L3LVrV6ovuuiiVC9fvjzVHiH1i9t7tK9a3rFTKu3G++KLL6baO196h0vvzuxzzuN83inYu4w/9dRTqe6s03fWRe193+tdcP298FMk/LEeyS464uafHz4/XnvttVR7F2XvdP3qq6+m2j9jfL3665RK54ovb82aNanOE9/398XrrK7w3rG4fB37NuXzutou7/Xi814qHXtW7a/bX4/Pv2o/4/396devX8m/+WkB//RP/5TqqVOnptr3PW7Pnj2p9v2ix+l9nvj6WLt2bap9n+XrQird1nzd+H719ttvT7XvY9evX59qn7/+2vLImqPlY3JZ87TZ5qhUOp/yvJ5GzFGg1eU6EA0hvClpv6TfSjoaY5wZQhgu6QeSJkh6U9L/FWPc3ZhhAgAAAADaRTX/18tHY4wzYozvd8z4qqS5McZJkuZ2/A0AAAAAQKdqiebeJGlWR/2gpHmSvpJ156J4tM3jLB6j8Quye2TNu9pV29m01os1NyKOW0TX3CLisUVGcDt7bFaHPO9o6JEmj+R4zKodZEWuOot0egzeI14vv/xyqr2j4yc+8YlUexdRj3F6vCmLx6HKu+Z6LO7IkSOp/qu/+qtUz5s3r+JzdRb57a7y2PbKlStT7V00f/WrX6X6hRdeSPWFF16Yao+cZnVePfPMM1PtkbDO4m7eEdfjir7+rrnmmlT7vtfjxePHj6/4PB4pLlpWd2Gfo97p1j9L/L3zbqYet5RKX9/SpUsr3m/IkCGpHjRoUMWx+vr27cC3If8M9DhjeczU56/HRvNsX42wb9++kr99vL5u/PV55Nw7Hu/efTyo5fugar355pslf/v76PuF+fPnp/ryyy9P9fXXX5/qxYsXp9pPA/Bty8fqpwT46QQevy1/r3x8N954Y6o3bdqU6n/9139Ntc8Pn+O+H6lW1hyVqp+nzTZHpex94YABA1KdNUd9/+fbvsetgd4g7y+iUdKTIYSFIYQ7O24bE2PcIkkd/x2d+WgAAAAAADrk/UX08hjj5hDCaElPhRBW5V1Ax4HrnV3eEQAAAADQK4RqY4whhHslHZD0eUmzYoxbQginSpoXY5zcxWPTwmbNmlX1YKvlF2v2CIXHI/x2j/x4V88sHscp73TmEU1fnsc3ymOCXaln99ms+/lr8mipR3W8znNB5+7Eemt5XpfVyS6P8uf3GOMbb7yRao8GZUU3zz///G6Po1l4vPO0005Ltb/+rM61krRixYpUz5w5M9X+eI+p+bbiHbA9Tu8dd/PE3j0yKknPPPNMqhcuXJhqj/96tMrjsbVc3L07PJ7n8S2fc/6+eCS22nXmt1933XUl//bkk0+m2iOGvp78PfWIocckfazeBdgf691nJ08+/vFSRGddn+9+qobvF/01eHTT54bHpaXSee1R4M997nOp9g6mq1ZV/v99/TQUX2f+ueIRQeedU6XSOKTXtcQyq+Xzu7wDrP/tkUbn88xj1T4vPcpaT76teSTW3xef4/79wLfNiy++uMtl+XP6tvn000+X3G/RokWp9oisj9X30R75Pffcc7scRx5Zc1Qq3aY8put8njbDHJVKT5fI6uDskX3Xk3MUKIKfyiRpofUVytRlNDeEMCCEMOj9WtK1kl6V9JikOR13myPp0SrHCwAAAADohfJEc8dIerjjV4UTJX0vxviLEMJ8ST8MIdwh6S1JtzRumAAAAACAdtHlgWiM8XVJ0yvcvkvS1Y0YVL145MVjVh5l80iE3+7RyzzxtfIuk/v370+1x048muaxqawulbXEUhvVobYRF5AuX27WMvw98viQx2Wq7YyY9/X4HPJokEdy2i2O67Li0v76O7swt0ezfLvzSKO/dx4x9LiSv1+vvvpq7vFL0pIlS0r+9giud9D1sfv2myey3yhZcTm/wL1H4XyO+rpctmxZqvNs7x6ZlEpPNfDa920eP/X1l1V7TM/nRp6IfqN4FNh5B2fvdtm/f/9U+zbh67ucr7+vfe1rqc4T1fPHepdZ34b8eVz5PtIj51lx3kbz97e8c6hHS/2z2fn+xV+DR2IbxSPavk/y7wEeMx81alSqfX/mUXeX9V3D7+/zQSrt3OwxUP+MmjRpUsXl1UvWHJWqn6fNMEel0u8g/vqafY4Czaqa64gCAAAAAFAzDkQBAAAAAIXiQBQAAAAAUKi81xFteVnn+/j5T1nnV3WHn/Pgrfz9nDRXr3M+a718i8s6F6L8XJSe4pfb8POznN9e7bmt5ff3S/0MGjQo1dOn/84p1G0v67I4eS+X4+dgfvCDH+zy/n7pl02bNuVaRh5+rtaVV15Zt+ctUtZ5Xr5vyzpnqTv8nMWpU6dW9dh169al2s8LzTpvvxllXbbHz1PPu759H+PrYPjw4anOOkfUeyC0g2nTpvX0EOoiz+vw/dmaNWvqstzyefLhD3+4Ls9bi3abo1L9Lm0D4Bh+EQUAAAAAFIoDUQAAAABAoZo7A1WAImIW8+bNS7VfsqDRl2ap9XIHRV8uoVoe4fvCF76Q6tGjR6f68ccfL3RMaIxqI6BongiZXzLnU5/6VKr9UgbPPfdcqvPE6Tu7dFZPmTJlSsUaKMf+DACO4RdRAAAAAEChOBAFAAAAABSq10dzi+aRsmqjr0V0zfWYW9++ffMNrMJz1hqdy/Nc7733Xqq9g26zR4p7i7wddNHe+vTpk+p33nkn1b/97W9T7V2Uffv1TtUeuXfbt2+vyzgBAECx+EUUAAAAAFAoDkQBAAAAAIUimtuEGtERt9bIbrX3P+GEE1Lt3W0PHTqU6s6imx7Vy4rt5RkHgJ61e/fuVD/zzDOp9q65Y8aMSbVHdj1+P3jw4IrPTzQXAIDWxC+iAAAAAIBCcSAKAAAAACgU0dyC1RIt7ck4bp77eBx34MCBqfZonsdxPXLrF7GXsiO8/hgAzW/AgAGpPnLkSKp9n5IV3/fIrsd6/XkmTZpUv8ECAIDC8K0eAAAAAFAoDkQBAAAAAIUimtuD6hWJbUSstzvP27dv31R750uP3fnF7Z1fuL78ud59991UZ0VzPRbsfBwAirF169ZUn3XWWan2mP3OnTsrPtajub4f8G3fo7kAAKA18YsoAAAAAKBQHIgCAAAAAApFNLdgHmXtTly2msfWs2tunsd6PHb8+PGp9jhe1kXply5dWvL3/v37u1yeCyF0eR8A9bVmzZpUX3/99am+9NJLU/3888+n+p133kn18OHDU+1x3FNOOSXVvk/xx44bN66WYQMAgCaQ6xfREMLQEMJ/hRBWhRBWhhAuDSEMDyE8FUJY2/HfYY0eLAAAAACg9eWN5n5L0i9ijOdImi5ppaSvSpobY5wkaW7H3wAAAAAAdKrLaG4IYbCkKyX9oSTFGN+V9G4I4SZJszru9qCkeZK+0ohBthOPkHoXyKzurj3ZEbfa+3sE16O5AwcOTPX27dtTffDgwVR7Z11J6t+/f8Vle7fMYcOO/wi/aNGiiuOYMWNGl+NGPlmxcq99fqN1rVixItVXXXVVyb9dfvnlqV6/fn2q33zzzVS//PLLqfau10ePHk11v379Uu37hXPOOafLMXmsFwAAtKY8v4ieKWmHpP8vhLA4hPCvIYQBksbEGLdIUsd/RzdwnAAAAACANpHnQPRESR+SdF+M8QJJB1VFDDeEcGcIYUEIYUE3xwgAAAAAaCN5uuZulLQxxvhKx9//pWMHottCCKfGGLeEEE6VtL3Sg2OM90u6X5JCCL2ynemIESNS/fbbb6faI40DBgxItcfXsjQiWlvrY957771Ue2Tv6quvTvXo0cd/OPeo3eTJk0ueyy92v2XLlopjKo/zvs+juShWvbpCo7527NiRao+0e/dZv33ixImp9m1Rkn7+85+netWqVan2CK6fauBdcH3f9tZbb6X6zjvvTLVH9l988cWKYwUAAK2vy19EY4xbJW0IIbx/pHC1pBWSHpM0p+O2OZIebcgIAQAAAABtJe91RP9E0ndDCH0lvS7pdh07iP1hCOEOSW9JuqUxQwQAAAAAtJNcB6IxxiWSZlb4p6sr3IYyWR1xPbK7a9euVJ988sldPjZLERHcLB7N9TjtggXHTw++6KKLUu0dbb3rrVTagdMvZD927NhUe4fWw4cPV3xe1I+vb6+z5mhWdBrFO+GEE1I9ZMiQVHsX6nXr1qXaY/N79+4teS5/zKhRo1I9ePDgVHsEd9u2bam+4YYbKt5+3333pdrnTd++fVNNp1wAANpL3uuIAgAAAABQFxyIAgAAAAAKRXauAB4V9cibd9AdOnRoqvfs2ZPqgQMHptq7UtYap63l8VmP9W61HjX22O1rr71W8bH+miXp0KFDqT7jjDNSvWHDhorL846faIz9+/en2t+T008/PdW7d+9O9YEDB0oe7/PGo5seG/XIL+rHu+b6tuXvqb8ngwYNSvWECRNKnsvfR4/peoR+9uzZqfaY/tNPP53qTZs2pfoDHzj+/4n26dMn1eeee26llwMAANoAv4gCAAAAAArFgSgAAAAAoFBEcws2derUVK9cuTLVHtN1HoPzaG5WZ12Pq9ZTtVFeH4dfoN7jmv6cHvWUSl9rVhzX48znn39+VeND9caPH59qn68evfTbp02bVvJ4f7+82+q+fftSXW2XaOTj2+DkyZNT7RFcj1V7lNdjs5I0Z86cVL/88supXrp0aap/8YtfpHrnzp2p9o64o0ePTvWkSZNyvAoAANBO+EUUAAAAAFAoDkQBAAAAAIXiQBQAAAAAUKhQ62VAqlpYCGlhs2bNKmy5rWDZsmWp9nO1/P3xesSIEan2S6X45V7ee++9zOX5uZbVzoGs+/ulN/w+fvs555yT6nXr1qXaLwMhlV7Sw88b9PMMp0+fXs2wUaPVq1en2i/74fPM56XPRUnatm1bqv1cQa+nTJlSn8GixPbt21Pt51973b9//1R3tk/w7dEvBePnB/ft2zfVfh4q54ICANCe5s2b538ujDHO7Oox/CIKAAAAACgUB6IAAAAAgEJx+ZYmcd5551W8ffny5an2CG7W5V788iid8ehrLTFd5xHcLP56yi8L4QYPHpxqLs3SHPyyH1nWrl2b6q1bt5b8m1+mxWOgXqMxfN/h+4jDhw+n2i8J5bcfPXq05Ll83+FR7EsuuaQ+gwUAAL0Cv4gCAAAAAArFgSgAAAAAoFBEc5vcueeeW/H2JUuWpHrPnj2p7tevX6q9u2V55Darq22eeG3W/b32+J7XBw8eTLV305wxY0aXy0Xzoytqc6IbMQAAaDb8IgoAAAAAKBQHogAAAACAQhHNbVG1Rlk9zusXpffao7Ye5c3qeDp06NCaxgQAAACgd+AXUQAAAABAoTgQBQAAAAAUimhuE/KOuB7Bvfnmm3tiOHX1xS9+MdXf+ta3Uv3II4/0xHAAAAAA9IAufxENIUwOISyx/+0LIXwphDA8hPBUCGFtx3+HFTFgAAAAAEBr6/JANMa4OsY4I8Y4Q9KFkg5JeljSVyXNjTFOkjS3428AAAAAADoVvBtql3cO4VpJX48xXh5CWC1pVoxxSwjhVEnzYoyTu3h8WtisWbO6OeTitXNUtt0RBQYAAAAaa968ef7nwhjjzK4eU+05ordK+s+OekyMcYskdRyMjq70gBDCnZLurHI5AAAAAIA2lbtrbgihr6QbJT1UzQJijPfHGGfmOSoGAAAAALS/3NHcEMJNku6KMV7b8XfLRnObPWrrsVEf6549e3piOHU1dOjQVDfjus+K8krEeQEAAIBKuhPNreY6op/W8ViuJD0maU5HPUfSo1U8FwAAAACgl8p1IBpC6C/pGkk/tpu/IemaEMLajn/7Rv2HBwAAAABoN1V1za15YQVHc+sVwW3nqGy7q2cUmA68AAAAwO9qdDQXAAAAAICacSAKAAAAAChUtdcRbXp54rhEbXsPf089MlAWH0iyorzlPKbriOkCAAAAXeMXUQAAAABAoTgQBQAAAAAUqu2iuc5jt0RwkUdWlFcqnUMAAAAAuo9fRAEAAAAAheJAFAAAAABQqLaO5gKNcu+991a8vbNOuwAAAACO4RdRAAAAAEChOBAFAAAAABSqraO53gF16NChqZ41a1aq6aYL5/Oks/mQNbcAAAAAdI1fRAEAAAAAheJAFAAAAABQKA5EAQAAAACFartzRP3yGVnnf/rtLutcvzzPieaUdc5nnvM6y++T9Xgu2QIAAABUh19EAQAAAACF4kAUAAAAAFCotovmuqzIZFa8NiuuWW2Utzt6Y/y3s0ulNOKSKHmes7P54Jf9AQAAANB9/CIKAAAAACgUB6IAAAAAgEK1dTQ3S7VdTquN8nZHEfHfZtao10mnWwAAAKD55PpFNITwZyGE5SGEV0MI/xlC6BdCmBhCeCWEsDaE8IMQQt9GDxYAAAAA0Pq6PBANIYyV9KeSZsYYp0k6QdKtkv5W0jdjjJMk7ZZ0RyMHCgAAAABoD3mjuSdKOjmE8J6k/pK2SPqYpM90/PuDku6VdF+9B9gMiohxFhH/ba5zesoAAAVASURBVDadvU6iswAAAED76vIX0RjjJkl/J+ktHTsA3StpoaQ9McajHXfbKGlspceHEO4MISwIISyoz5ABAAAAAK0sTzR3mKSbJE2UdJqkAZJmV7hrrPT4GOP9McaZMcaZtQwUAAAAANAe8kRzPy7pjRjjDkkKIfxY0mWShoYQTuz4VXScpM3VLHjevHlVDrV38vhqOyt/ncwPAAAAoH3l6Zr7lqRLQgj9QwhB0tWSVkh6VtKnOu4zR9KjjRkiAAAAAKCd5DlH9BVJ/yVpkaRlHY+5X9JXJN0TQlgnaYSk7zRwnAAAAACANhFirHhqZ2MWFsIOSQcl7SxsoeitRop5hsZjnqEIzDMUgXmGIjDPeofxMcZRXd2p0ANRSQohLKBxERqNeYYiMM9QBOYZisA8QxGYZ3B5zhEFAAAAAKBuOBAFAAAAABSqJw5E7++BZaL3YZ6hCMwzFIF5hiIwz1AE5hmSws8RBQAAAAD0bkRzAQAAAACFKvRANIRwXQhhdQhhXQjhq0UuG+0thPBmCGFZCGFJCGFBx23DQwhPhRDWdvx3WE+PE60lhPBACGF7COFVu63ivArH/L8d+7ffhBA+1HMjRyvJmGf3hhA2dezTloQQPmH/9r865tnqEMLv9cyo0UpCCKeHEJ4NIawMISwPIXyx43b2Z6ibTuYZ+zNUVNiBaAjhBEn/KGm2pKmSPh1CmFrU8tErfDTGOMPagn9V0twY4yRJczv+Bqrxb5KuK7sta17NljSp4393SrqvoDGi9f2bfneeSdI3O/ZpM2KMj0tSx+fmrZLO7XjMP3V8vgKdOSrpyzHGKZIukXRXx1xif4Z6yppnEvszVFDkL6IXS1oXY3w9xviupO9LuqnA5aP3uUnSgx31g5Ju7sGxoAXFGJ+T9HbZzVnz6iZJ/x6PeVnS0BDCqcWMFK0sY55luUnS92OMR2KMb0hap2Ofr0CmGOOWGOOijnq/pJWSxor9Geqok3mWhf1ZL1fkgehYSRvs743qfHIC1YiSngwhLAwh3Nlx25gY4xbp2M5R0ugeGx3aSda8Yh+Heru7Ixb5gJ1awDxDTUIIEyRdIOkVsT9Dg5TNM4n9GSoo8kA0VLiNlr2ol8tjjB/SsTjRXSGEK3t6QOh12Mehnu6TdJakGZK2SPr7jtuZZ+i2EMJAST+S9KUY477O7lrhNuYZcqkwz9ifoaIiD0Q3Sjrd/h4naXOBy0cbizFu7vjvdkkP61i0Y9v7UaKO/27vuRGijWTNK/ZxqJsY47YY429jjP8t6ds6HldjnqFbQgh9dOzg4Lsxxh933Mz+DHVVaZ6xP0OWIg9E50uaFEKYGELoq2MnJz9W4PLRpkIIA0IIg96vJV0r6VUdm19zOu42R9KjPTNCtJmsefWYpP/R0W3yEkl734+8AdUqOx/vD3RsnyYdm2e3hhBOCiFM1LFmMr8uenxoLSGEIOk7klbGGP+P/RP7M9RN1jxjf4YsJxa1oBjj0RDC3ZKekHSCpAdijMuLWj7a2hhJDx/b/+lESd+LMf4ihDBf0g9DCHdIekvSLT04RrSgEMJ/SpolaWQIYaOkr0v6hirPq8clfULHmi0cknR74QNGS8qYZ7NCCDN0LKb2pqT/KUkxxuUhhB9KWqFjHSrvijH+tifGjZZyuaTPSloWQljScdv/Fvsz1FfWPPs0+zNUEmIkig0AAAAAKE6R0VwAAAAAADgQBQAAAAAUiwNRAAAAAEChOBAFAAAAABSKA1EAAAAAQKE4EAUAAAAAFIoDUQAAAABAoTgQBQAAAAAU6v8H8WSQxMWisyAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "filtered_image = np.zeros_like(img[..., 0])\n",
    "# here we call the compiled stencil function\n",
    "compiled_kernel(img=img, dst=filtered_image, w_2=0.5)\n",
    "plt.imshow(filtered_image, cmap='gray');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Digging into *pystencils*\n",
    "\n",
    "On our way we have created an ``ast``-object. We can inspect this, to see what *pystencils* actually does."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.38.0 (20140413.2041)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"684pt\" height=\"228pt\"\n",
       " viewBox=\"0.00 0.00 684.00 227.95\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(0.478879 0.478879) rotate(0) translate(4 472)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-472 1424.34,-472 1424.34,4 -4,4\"/>\n",
       "<!-- 140128518376696 -->\n",
       "<g id=\"node1\" class=\"node\"><title>140128518376696</title>\n",
       "<ellipse fill=\"#a056db\" stroke=\"black\" cx=\"548.645\" cy=\"-450\" rx=\"107.781\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"548.645\" y=\"-446.3\" font-family=\"Times,serif\" font-size=\"14.00\">Func: kernel (dst,img,w_2)</text>\n",
       "</g>\n",
       "<!-- 140128518374232 -->\n",
       "<g id=\"node18\" class=\"node\"><title>140128518374232</title>\n",
       "<ellipse fill=\"#dbc256\" stroke=\"black\" cx=\"548.645\" cy=\"-378\" rx=\"31.6951\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"548.645\" y=\"-374.3\" font-family=\"Times,serif\" font-size=\"14.00\">Block</text>\n",
       "</g>\n",
       "<!-- 140128518376696&#45;&gt;140128518374232 -->\n",
       "<g id=\"edge17\" class=\"edge\"><title>140128518376696&#45;&gt;140128518374232</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M548.645,-431.697C548.645,-423.983 548.645,-414.712 548.645,-406.112\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"552.145,-406.104 548.645,-396.104 545.145,-406.104 552.145,-406.104\"/>\n",
       "</g>\n",
       "<!-- 140128521227960 -->\n",
       "<g id=\"node2\" class=\"node\"><title>140128521227960</title>\n",
       "<ellipse fill=\"#56db7f\" stroke=\"black\" cx=\"52.6453\" cy=\"-306\" rx=\"52.7911\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"52.6453\" y=\"-302.3\" font-family=\"Times,serif\" font-size=\"14.00\">fshape_dst0</text>\n",
       "</g>\n",
       "<!-- 140128518531728 -->\n",