Skip to content
Snippets Groups Projects
01_tutorial_getting_started.ipynb 190 KiB
Newer Older
{
 "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",
    "*pystencils* is a package that can speed up computations on *numpy* arrays. All computations are carried out fully parallel on CPUs (single node with OpenMP, multiple nodes with MPI) or on GPUs.\n",
    "It is suited for applications that run the same operation on *numpy* arrays multiple times. It can be used to accelerate computations on images or voxel fields. Its main application, however, are numerical simulations using finite differences, finite volumes, or lattice Boltzmann methods. \n",
    "There already exist a variety of packages to speed up numeric Python code. One could use pure numpy or solutions that compile your code, like *Cython* and *numba*. See [this page](demo_benchmark.ipynb) for a comparison of these tools.\n",
    "![Stencil](../img/pystencils_stencil_four_points_with_arrows.svg)\n",
    "As the name suggests, *pystencils* was developed for **stencil codes**, i.e. operations that update array elements using only a local neighborhood. \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",
    "But lets not dive too deep into the concepts of *pystencils* here, they are covered in detail in the following tutorials. Let's instead look at a simple example, that computes the average neighbor values of a *numpy* array. Therefor we first create two rather large arrays for input 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": [
      "3.84 ms ± 36.2 µ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": "iVBORw0KGgoAAAANSUhEUgAAAdkAAAAoCAYAAACy/oAKAAAACXBIWXMAAA7EAAAOxAGVKw4bAAANAUlEQVR4Ae2di3HdNhaGJY0L0Ho78HagdTpwOnA2FSTuwBlXkHE68KaC2OnA7mAdd+B0kFgdeP8PAigQF3w/gEsdzOASjwPg8P8BHgIgeS+/fv16Ya4+BC4vL19Lq7+8Zrc6fhRXn4gr75kON/LfyP8g/9THdbi4kNwvLhD9JPWR8yHUF4lZcCYCCb6L+ErqQiPjaiYvuWIJvou4ov6kPpKML1DIuASrRdgnddFanbhjZM3XhYE6yx/yN4EXhd/Lv4nirwn7dGSfR/EvQc6nPVH+5yATyf0Ry4Ww5J6EcA1H6dPgUIM+OR2k4yp8gb28cbXhNWktrugHNfElXWycaBx2jM/NrmljcDcDu+GAzhE+lOYHrsTuDT9Eyl+T5sPBqHJBfhdkFXYz3BD38sg4oxzSKSPv6ghpXrYlF+dNCatu9H0/toxkX8ozc3fHuJzSnss/i9NqCks3DKNUWs6X6pnNlcp2YhjrtkZYbWX5VfqD4Qoc5Xbny7BfhrvnLdt/u8aGMO8cW8ob7PNmZKOLYxfIe6aLtGt5etIX+TfyLQOjeDC2Qa6VH+vqy1OXKxPnpWHJYHgX3fGpPJ0XnfGf0zZyccn9KN8Yd4XptM2NA2XkqG+Rbrm210iTXoGHRXz5c5zFlcoOYrjSuQ7y68/j0FxFfXI3voSrYX9/LZiL+yCG6TgZM7aG+rwZ2cqMrB/AzEhZIqYz4U9mhUrDGEn8fgaVhpXPhf+kbEaOztcYujR/alx1of9YI8tsoHWjoHhrECmOIWsZ3qk6bSmP/vKL+FL52Vyp7CCGa56/P98sv8o7PFdgKVeEr4eO/RLcwxjowzDIhKNkB8eWZHr7/JUEzFWCgDbyn8izjMsG/rfyl1LtJ/ln5CVqfqu4exAqSXdRyUM8vlMmKvdK4d+i+C5BryPn9XfS4K3iGC7nhANxsEkx8BJlDugjv5gvj8MsrnzZQQz3QujoXIFjrXwdHfsluM/p/2PH1hDuV3MatzKbIcCSLXuTjROB7klhHf9sEu8C/9HhQ5LWRD3xxLkTO3HqQMyEg+OhiTHGOMivdQxGEyMaO4zu4zhBYW4CGsOb5JWKrsLXQq6mYLgXToflCgAr5+uw2C/EfU7fnzK2OnE3IzsH+u3KYFyYoTZOxpBNd2azjVMa5DPzgdg+919lpvVdqzx7nM5o+7r66tgjL53J0ibnFzv0bZ1LnFkovCZfS7kag+FeMB2dK3Csla+jY78U9zljYMzY6sT90ZwWrcxmCLxQzc9l+GiA2R2G5k/dwf1OQuQwsqR/itJOgsp/obpe45UZZrQ8BEU7wbm6QqTQkVkr5xtcOoslPeARZGo4rsbXClyNwXAvzA7NFSBWzNehsV8B9zljYMzY6sTdjOwcyDcqow7E3dDJhyTS5iTHMvG/0vRcXLKtWXBORml0kJaTYWZpdkzZ79TGSflWZd0Rzhd3fXdofomnNxDcTeaMb1No78DafM3kagqGDUQb83t4rgByT74a4oYDh8d+Ju7DyJ1KTBlbnbj3GlkNxBu1yyzoqfxbnVw8AzpVyVLOEQE6R2rkuIBgyDv3fNc4UbVxqz5GR84Zz49JG8iETp9kPZjoCVcTMWyA2phf4+oO6dX4aogbDhj2dw9SnlzThqFrS0wcW5249+7JqpFP8uyDoTCvKMxyupB+lmdZ0lx9CGC41uYm28HpA/Lv5eN89oe/C7AojweyfqeDhzR/rGFZO1Fp92gXV70YduC+RPmYv1w9xtUdKlvwZdjnelw7rQv3IJXFsGOc9I6tUKGOnX3+KhLKBtUwy4a4WbMaf9Hk5XRO3FxlCARjJp6yHW+Kur6T0ilZ/cCgvpPnwa3gaINVkdCnmDGzPM5NWNg75tWlxuiGgjr+W372jV5UjwuqPfa+Y91SkeriXVyNwPAE9zknJ7zgtI/fUO2D5wog1uSrBPZq8+zGSB/uIzA8GScjxtZgn79UJUEoe5RiXDB5F3DUHmBaicozKJ+qPAPvbJ3Og6/q8LDRrJuNmk9c54axyT1gtYnaao/+lD7M1duWyvBd0tX6EDqoQW7+BvfAexXbOXMJV3Nwn3N6xtU9anvztSb29BedydmNEdCvCffBmaz0ZdaxxLAsLQ9mNTj2o9N9whr0WqyDNzS77Lf7gftpitK+DDdrD97N5WoO7nPANq7aqO3Jl2F/j31NuLeMrEhyy0E6snT3Uv5GauOzy3SRHLI/enn3V2wKv5GnHGvVNz5+Vstz0ts56c6NAn81d+uTjniAc2brW7vHwnH01oF0YgmHJWTejzN3h8AcribhPgdo46oTtc35Muyz2NeBuy5e7B2gIYaE73E2H/dWGCPZ+o5sJN/7916+TpYbFOz+vq6Xa9ockp2ar/YX/QWUynOR51wHP7I/Vbfa5HWOYXnI9Yka9JNO3Jitjr0/15c1nOMcHYyr/mvKHEy3LLM1X6p/9XFy7mMEPmvA/RFa+LsgPhH3kxSLZxmEecK4NYOTvJudKj1e9uOdyrisou4LPbEMaS2nuvgwffM+puJ0ln/K/8UxzmsVjCIDZZidcxc/eclb5ZjF/yr/g8q3MIiaP0xQ5zhpn3SPE5dOZ7VnugcmtGFc7YX0Ou1szZeNkzxPNeDujKzU4+EmZmxvE1WZ3eYuvLz/hXFm5ksZ/iElZ8S6ylOc8s6wu8hdnOXKxrAqn6fbqDv3tKkrpvzeMoAsGZauebAnvQkITZ8cJQ8mGHyM66+Kn8iMTPgZHUbKmtjKCHge6Yepo78/Vv73aYbi9JXOPpeRt6QVEDCuVgBxRhWG+wzQphTRxQTxk79tUhoXITJbf0OGvC/DhSssJyPX+ks1xYfKM0ts/b2a4oN/KxTaD8cxZbwuk/8qTeUwsui02XJ2dB5gaH4iBgG/qUdhzdL47OVi42peX53KE/LG1Tys4z5aAnfP3YO/pj3SXQzGEJ8u67q7fwHVmqFKPiwVMztzeUrDGLHJjDEKs8WnSsMgN+VpS/Fb0uVeyf/sQvohTwfqdrPkkK4j8tkZ8dgytIlu+Ei/qIl8ULK/qAz6M+Ne7fWRXGtqa/ZUOVefpW2LgPG1Lb5r1m5crYnmtLoM+4uLqwiy/0VhgiyjOcMrQ8OyLQYQh8FhKbVxAtLtm+kYDCx5fCmqMdzUoXiog3weSGryo7xghJHBYXQfu9DpT6hvTJnfVDy3bHhaa5TidWS5uXXOkYgFDQFDwBAwBAyBLAJXMiIYKGZrwWAxq2Sfk+XcYDS/iQwoRq/1l2OSZybbPLykMO5a3pVXPmFmkc6oKt60hWDiqD91lO9zY8qgS0vvvgrjPOnN6yPcJJgzBAwBQ8AQMARGI/DIS/KQBw/3YCxxGF0MEkvApLH3GhwfLWBmSxwD7YypDFH6cA8zP2aA7t1L5cdPiWJkgwFXsOWYtVJvcF2z2JDPcUyZoGtcbkqYB6jYn26Wv6cUNllDwBAwBAyBh4eAM7IyHBig3NOUJ2mSxTjGBjOLmpfrmznGhpQ6gtG9Tiok7mbASTrRKWWY7Y4x2Jlm3N5yOlPPylmiIWAIGAKGgCEQELgKgZ2PGLyWMfWGHqOZM4Qfc/pNLEO9wSjnqrO0HgQ0i3df7eoRObcsViQO+RUp46q+rnimnBxijJTGvpSRxdjl9mX5Pm0zexY47IM2f3umOE8I21+llbmG8MBb7gaojDYLW+UGzd+kLaypyuLGVX20nB0nBxojRbEPe7K7dknIk7F0r+3EFzqFeWWG7yCHJ3l55acxulKS2S+vBvGUsNsDHlFGos7xCg5gm5uIAJyoyGEM7MTTPytx46o+uoyTcpzUgH0RI+shZ9baGMtAA0YzhNOj8tib/YeAaz3p21cmqoO/29vln2aiNs8+KKxZcbiVt6X2ytk0ruojyDgpx0kt2JdaLuYjFRjTyUZPwGFgux6EyjLqy9hfpWXRGUx8Ia4OuXc5eObnJ2Bc1ceZcVKOkyqwL2ZkPe51/BVRuU5Qdcu6OeH1K7s5qZqlO+WMq/pIMk7KcVIT9kWNrGZIPL32twBhSXKUmzGrwlDY6zej0L0XEifsf7MnbsvE97BUGTKu6qPFOCnHSW3Yl9yTdSzoIp5+xGJVdlR/5x7vqg0dr7JXws5uTs6DV+OqPp6Mk3KcVIV90ZlsOQ6s5T4EdCfIvjffejZXOQLGVX0EGSflOKkRezOy5fpDlS2rk7JMzLeqJz1cVuXJHFwp46o+go2TcpzUiv2lLqblULGWq0NAHZVZ7PfyfJUrduxt8yrPW/nPtgwfQ1MmbFyVwb2vVeOkD51t82rF3ozstrwfpnZ14C86mQ8yrvHHQQ5zfkc6EeOqPjaNk3KclMbelovLcX+OLbOUbO48EDCu6uPJOCnHSTHszciWI/0sWtZdIO8y81eHdNJnCr+Td39feBYn8ICUNK7qI9s4KcdJLdjbcnG5PmAtGwKGgCFgCBwcgf8DpsUufKsWZaMAAAAASUVORK5CYII=\n",
       "$\\displaystyle {{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+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAClpJREFUeJzt3W9sVWcdwPHv0z+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": {
      "needs_background": "light"
     },
     "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 than we needed 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": [
      "639 µs ± 35 µ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": "iVBORw0KGgoAAAANSUhEUgAAAKoAAAAZCAYAAAChKLVZAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAFdklEQVR4Ae2bgVHcMBBFOYYCSFJBoAMgHUAHUEKgA5hUwEAHJB2EdABUEKADSAWE64D8p0ge2ZaNvfb57sCaMbJW0mr3a7Urycfk5eVlxZImk8mZ+q3r2dDzqOdEvKbKxzQiYEagyq4mFkP1zC7UFwNdUflS2YbK22YJx47vHoE6u1o1orNf6Heq8pYGwruOaUTAikClXVkNFUFiowwhP6ZZhR37vW8EYhvK7GrNgolC/Gah35Yv3xboY3FEoDECdXbVxaPGAnCwanWYYpug5zBmMr73g4Bw3X8j27DMrkyHqRhOAeJO/1oNRzG97l19uC24VJ+9unZjXRkBYXcn6oWen8JvqjLRDOyvVP6l3CXROeDiPNyB15OXJpP8ObvqZKhihkfcFhiNjRSk1A+wD5YVRHSYVxJ2zxqbhR6nc2F5EhPUjjY3oi/dTUzKrsyhX8x2BcRmMFLKesJeNcYs9+6FuB3KSDUeYfA4J8QcCz3Ig4c814P3xDiZg5yRirYiGgeR66669yAv4jROGi9pV6bDlJhhkLjmUxTxUuBVD/x7XQaoY8ivQ6i+7jFlmBVdTkX/owfDXvhUZ1cmQ5XGN3rcPjPW3q/imJR7lyDu6kHtlnLflFNmCQrMhzB/1LOr9+slELnSrjJDlTIYXjiFf9E7ng/Dwnt+0vMkZd3KVP5BZUvC+1YC1kYGy+Bt+0TycB33EPSHj+rAhr32Z9EJs4MljR3miTlDNr4S3lcIAN5EukrcK/r1Qo4whJ/drqQgDFD4jHdfBogHPce+zPtdqLfm4nGl5zDVX/SZyCC+LA6nR2rcOpr6OUyUs296jtvCU49I/zFrmquPWR7GUALD9TCe3lkwVOA13fzFuejMpXnu1Ncsr/r2NqfBo34TU/YzIbFKAeC7J7D/7CNcw/Ov51nMhpKhOG6y7D3mb1+JRyrqzz57cC8lI8zt71UmtCMHV1bFDzEiObwxmHmk3ubUXU/hnqXwNGiiMquW1ZkDJdRbc/HlaoVrqdIEd5VB/Tnc4fmKiUn6qKdoaLTjYIIRllIsj5f7VG2zQ4loeDHuKTNazKRveWLexXeNhZHiObkByOmpOrZueNRJsV9cnoW84tmfXUmBVLhgEkzhMsUv0MQTQ02GqNAm5GrXiwziYw5dyIK8XpY43DL5VG4FeZvm6mOWR325xC+FcdFYpMhTwlY0J2tT+Yrt1N8sb4KXeU5XJUgu+RUIreT1cg1tBcL+q2FoxjK0lRyPey/Qs4ijMsY7Fa3qANN2jKbtUxGDvgHT1G8tqqJJ0zF7add1Tp2higmX4kFZB0Y8Cao7juq7CE5YArhSGlCG0tivENhX58KpynPZn2pcPpumvjQxZ8XFJJJLyB8vMk+efdbnnK6KGUoSUsJqze1LVe8ULXgUq5Z4oBLQA8vQVvackXpZwcrt49sy69j+TOOzH82SyuxNWfzJvbbo4D2L6Ci21anvOV3TUIQLPsfxayauXFD40APCPSF72HD6p9glMbk5oD2zIWVoKz/3yT88NvQNJ+vBJ1/zwAk/GCteEifCdqruLndH9V/1DJ16ndNOP0qxaC6guZPdA3RL/7Z9NB6HAf5NJnk6N/Bjoe2IXykyNOHVtzx1Y2ostnMcvsLiqmuerBtS3qQAnjgPQyVUtf7FVZ0SdXV+sogKrfdp6suWKPtfMM+LmwsWmsmjdpGnTs9UncbiNuBJspoX6ZDypnQItMENlYGlPFuAIwE4iFcNyrbNJSfen69Tbuvj5eZ3n+aJbyuDtb03sKX8mV9K53kZKiFp4X84rcnm0MQ9JInwWfdN3TValD9+US28M2iK11wMFeEEJAcBLqj7Oqg11fnNtxO2bK+uFz1itZmIf3yMaza5icDpAAAAAElFTkSuQmCC\n",
       "$\\displaystyle 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": "iVBORw0KGgoAAAANSUhEUgAAAIYAAAAYCAYAAAA/FYWiAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAEqElEQVRoBeWagVUUMRCGOZ4FoFYgdIBagdIBagVoB/isgIcdoBUodqBWINABdIBcB/h9YbMvt+zevd3N3cXHvDdkk1xm/kwmmckuk9vb240hNJlMjhl3De/AT+ADZE0pi6cK+xZAt+Er+OP/gj23cTttoWP0ZcCdwO/juKr+M9ZLLsGqQ29HjDyfwuex/pDKebbYpHMoeVJEOufhdawUXu438B1R32XneHo8NOq0xaMhlmBXfWiMe079V6Ot5GoMIWKM4S9tKxl7bmzpvGtbDHKMFBk7zVj9Fn6Vtpf6jFOnJ50wdyusZ6ViXhauebbYHKMUpzhk/G/4CCUXY2Stcaw5x4NNPht2r20xYUEbff2rOIjJqAndXv/R6xsBbg2xBe5maFwfqAWaK8zeBqWn8Lccm/KeLXSMsQw445Qetj9W1qrGg/U9fLIqfWP1gNWQHZL8KEv8tsX60BIZ92zRO5TgWVvwDRxjM3Jr8n1Gb0LWPmxYWgmhyxvUDoYMJ4V1uG0+2fAgf+wcvVafgTlN8t2Qf8eA7LJF7+QTYFOECSZmsOLS0Na/WymZKgcwhBy5WBVWHeRNqbirxQvOnGJkLUaF7nm26O0YFTCN+AHBxjrjnJ77XKehLJ1Mlj2W3YE1FY79I0CnYPQtbU7qtEXtGCyyxjLWSC9hwbjgHrEu/jXAPlOak3gDKeYW0hP7Y+cwjxJ5Xm0v47wdQ582MdY/o31VG8HT4gLdrsU72A0ptlMwpKFFfH3WsdsWCEZ+2EHHPld1HeQSPqzqPo9OcqL8Zolsj/Sgq9m3qM44DZEVe5RH6YLcpBiomwvR1C9pZ8ygOTLO+blIOmP9GUL9kOtSJ/w8Z7NFdIR4bUvravYqRxEMVH9f6GuURb9H/iCjVdiyYgeLJ0IwNuW9rJ+2n/KiOTX7lQn3dn7GRMdA5Kwz0ie+G39jH5TNFuE9hscPgqdKlqg7eZ1kVHKjjJSQK3B3YZOcvDeathh6BY7OxDA39lQezxrdl3chhAqaNlfAF2J1m+2R6F/GHNWpHWbe2qIrXDPp26PvF/V864jAcEqkJYoE0tu7Uxl9ntE1aDe16ciFHTk6sHYIu1FdkDHeh9023fPaGDN4joxtDeW06xjimQkxEUfVN2gdNxk8Q3idk5dmkpq7prL/ZsbuKXWBkeuTlLrOMqVt1Ym3+gxxXXTvpB1ri+AYCPHli8e5FI76dPL0HSb9d78q5O8SsbsQTYMbWtexYY7Q64vFuEbR+n7V1lEDppy22ESYjuCdPjgE5UxeQb8GUnm6c2haPy0Z+4xTVLq0Uci/Vjn7apP+QOenqBc8OolftQ9sy22LCTJV8BX+A0tfYGOXiY5XJHMQ25ZGTMr4662nNaHrUlwZZynYW2Rrj2AXcM44TRe+tH3oHBsywq2DNjepG7b+qt2Cd9Q6Zvm6moIf8pzDaEP09hkDRq+GL3AKj+/e9D/MMZ1UKY7hqeXJ5E5YO7GIhlZPsOAE1MXn1TVcC4cArGQUM8dFc3i06Aer6C/FIZK5ejPz2I6ko/juYnDiWeAc49xayyJOjFZka2xkd5tkxmu7uYX/t7HqK+oaLbCx8Q9xlEVrw3lWOQAAAABJRU5ErkJggg==\n",
       "$\\displaystyle 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": "iVBORw0KGgoAAAANSUhEUgAAAHoAAAAZCAYAAAD+OToQAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAE6klEQVRoBe2ai3EUORBAvS4C8HERHM7AHBlABvZdBEAGpojABRkAEVxBBvgiADsDTASAMzDvaaUtjXZm16OZXW/BdJWsX3er1a1utcY7u7m52auB2Wz2CroDygPKFeUFvK6pJ9hBDcxqDB2N/AZaDbxH/z3VA/oPd3CPk0hoYL9SC8cF3Rn9Iwyud0+wgxqoNbRbyY2aQnY+toPb/X1FulezdUL0YUF3FPufi/GpuyMaGOLR+RZMzHolY4Z5yrOcydQeRwPo9bi8RquSsVwcGIbsGy9/no+vakNjtv4emier8Ka5bg1EvX+LGH9S/4c+LxMF8ybIOl9ImPdoVBcY6ZFm3714QHNBMUvvRTfhh6ewTqL+Hid90H7jWOpbAwEvjVWHbk7MY5gdwih4sn1KuqtdqBXA8XB8hm5+0lqxxhtkPcPY6Xgch3EaQR49Vf2dZ5KYBH/P+jqQCfJ52ntVMgaxBjVknyl4XECDn8T2quoFk1PIXqWhjjl0rXMFB8tRMGqXPn32fqW8rjI0hP9Twj1LvYB4ihb9soGg4fkF3la8uVz/F+jrJNe31Z/2QOdXHpCFoelouJQFP6ItUw2j93rZf4PwNbVh4Q/rCtD785DTYNFHhgbhhjqZPD4nv6T9uxxz6sa78i8V6tgWQG++ZG1t8i/FZEzZTGy79Or4SUiGaGjkVyCnvgb/Qjl1LLYbl33C7VPD5yPlWRsN4xuRAb4errCPtnVXjUEXdEKtgn/kuPKkMNQvoYSmSh7o1I/G8HA1dEhfWx23ySIu5SJ59Es6xvMEnhJP7Ns44P07RriVZyNpiPyttiVDtmR3M3rsp4hh7lHu33uxy4u6GQ+fOcKgyS6Jm3K8Q+Zz5sroor4PwjsahIMcgb6e53HtuuSd7g3w/QHRCXyXFDRUBuhNDvW8EvSE+5TSUOJdIUtrApnLE+U+AzdcXRIypnf5Tl2MOZ5gbHnkG9dU5saXScb1Wp9YT5hr6JY5w/xFCNVMNmom3ERVuCt55X14aujF+y+fK9tjyQCfqlCZ5FHeKIvOEPREX+WpI72robt1fWiq5YHWEL10hTKmoZWnEdKVBQiy7tvKIZ4AhxonI8cZ0A5hZB39hmVYt3w5r8dforQ8JGp8s9/LEnnDfdfz+uuCtqgVolkwNIr1o4IhTgjhL98Ec6fZ/Byr7q+CuPASbFGGpbXXDKjYUoF3dT+bRx202MLfAXjw2pxT+a/3IdKwfm0JBqZu3MvMB8TiRINWBZ7IpR8nbFmGvoI3jBxlVVchj+nLbAh+dL4P8DBxDYA8Oug/lKfzkaW/6vt8xh8R31FShmlGZ8z3wjeV9w4qszyHe0NUkt/Gy2RiYzKwpnei39Vbk6Z1m4iKzPWj7EE/8GwcgnW8nB8qT+QR/pFE+5qiI5ootl4jrKcNnw7+7xVMegELm1CYHfZWUq+FIvIYis3XhZ/Z7d/IvxSZcryu9tjydK3jOGvpQCZvh3dhaL3hIYv7Nt84xM0alTz9vQBar7TFb+EiL18OS8+Y2zIeIs9t10h4rKXnz79oooBez4Mx8Fnc+23n/02JjEafxZMlyj36s3MMnZY8kDV5c7Dv1j0aAVJI2fkfHuARJl2+QwXvZvOL1rswYOzQH2TXmZ4jb7gi78TQ6gNBTCL8eDJKoifPCeYaQLdej34OXeRBPwGnY432F5GebwAAAABJRU5ErkJggg==\n",
       "$\\displaystyle 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": "iVBORw0KGgoAAAANSUhEUgAAANoAAAAZCAYAAABXYTDBAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAGi0lEQVR4Ae2bjXHVOBSF8zIUwEIH0AGwFSx0EJYKCB2EoQImdAB0AHRAqGCTdEC2grDpIHs+ja7jH8nPT5btZ2zNaKzfq6Mj3asrvWRze3t7kBI2m82p+t1XfKR4pfhWsm70XcPKwGIZiOnFJkXRvLCP6ouCHSj/VZ9Hyj9dLMPrxBfPQJteHCayc1Tr9175JxqI020NKwNLZSCqF6mKBpFlpTKXsVy2VLLXeS+bgbIOFHpxL4UTuYiPa/2e+Px5rXzNrgwshoE2vehzopUJ5GFkp8cQ3EzF47KQNR1noC9f6n+EjPgI+1UzN7wR9u70Qlp40CdqAITxMNJZjtrzWvl9lz5LbpuLL8mxR6vOazUE78JxoYiRvY98BTyij4pH5fGU3wu8ZUxd08Je0YukV0cJcUFWB7KeavA3vqjTR/0g+qX6uVfLTp0W3CgXX5KDgfsh3id9HRaO/4QDLOXwQbjelgv2BW8ZU5e0cDf0Itl1lLDnGvSxKRl5RburRfF4EOdjKZnGw2U6iQLa84qcfIlzLudnffnIwCkG9oPiN0WUi31UUTKVceruC17gdAriJqgXSY8hEoZCcTS+h3SPgFPtpU+3fSD0RVuDta7CQG6++CnmX0U2+lThKqRYETD7gDcCrVrcphdJiibxPxQ5+vGhi+AtUJGvJwTEXcbVbnUZ6+QE8kPwxRpJ7pXic6XPAsPuVdHM8Eb1olA0EY/i2Cvgn0pjSVEMTq+HiteatLOC+v6hfErg9Isu7i4YUgZP6eMxcXrjxhDgiYcc3B4XfJt3ylz7IviiTTFXtYFH3AqMzANF5LxSm7b70lB8gQvvo8Cn9KhBfNhegwd+LuJB7TICYjK8fm0Na7peaHLMjcmekvZ5BP9UPPF50hdWn/qVjO+Kx6H+Kh8Eg+SyWd08QuO2lakfhoaLe/EapjQKQxl/cqaPy8OPy5s85Zmr8cfcGq+sKvtp7UNfLyM7X5LL+iavp/omc+o5gxv34ujz8AyZnLKNF1GVT4IXjIpZ9MJONKwxvrAFLAyT/+QLuH/lcPeQ+cvLrH/GwlAfty2Pa8x9oji9lGcOBBaBQJtvalPnB4/gQhYRa8wJ9gzrqHY3SlvgSbstDMUXa2D428YfpE4cVO7ocOd5go/6H0OAYSq82fake96vbwDlsThYlgohlPUJkstJwLN+w2Xpi0H9ce9wzeqBDcVGrysC7VCi4AMOeFQP3k9qg6FpBLXhdLOfKsrK6NqqHivtnq2Vdqee8swdfkPKqeK7oD6D8GW4Na/N3WjNlNpl5bQ5wl2JxkLJOLl4gays1VR4NW7FMCqfrheaVOioZoMkuVsheVYmmWycoHtgbeyrdlkwSE6Sm6N+KFErBi+bNsE5+f7OZVQaxWXjopj0Ibb+0K/6QfiSXDc343rXr/oncco4CngADbdVZXAT5FLlk+Etc+PxJenFoTpXgrcelDVOnUrDtEwnF2BgDF2Rm1UNuTImw9q0uWG4Rc7d1KLxZ2r8wM9Jwil5rLrQKWzyh+LrgQYw7DbWWN/YfI3D8wCQKfE6OH33pFM0CeFHXZuoI0Kb4dImrLqTUr0Vp3xZXEhrhBExNMYOFWj+NyqHg2ehevB6jmjXcLFVjxUmYMFRNNyiIqgv91/cTaeERUU1MRRfjAnuKcIXzT300sq+u/S813FNgjfnnjyUMCbIZjBLU9k0qneTjBBQJ2Rbno3bIHlkDNswluu5vzX++Fl4cXPMEP2l9N8qM8Wy/p+V4H5mnsE7tTFjZm3IW72Vlb9D8cUatI1bxpA7fSoeKo9AymOEMMDB+7LKR8crTFn1AheGxWZT/KNIwNK6S6m+3Ce4w9nrI9nk4MFzL6m4YyofDINkc5/g6T3pLyFK2H5JDg8aDxWZQ+F6qQ3GiFfGG0UC8/mqNm4zq55Fw4DRnzrawUHltzblK8H3y86X5LKur4XPjEVl3G0Z9e/LaZkv0nAb/e+PKfBqTNYpm170+qPibQsSqtcE2Gwvyhs11C5XWd9NkQtHqpzcfPkNxGNExdjtgm9MTueGN8bjFIrGabnzX/zHJrCt3C8UpzKnyOyC8GflS/Jwe6/FR9IJD4Fjcjo3vLENNrqi+YXi94g3Y51qscnPpVybLQtfXkEm/zeZrrzPDW/bvA7bKges49JbuRAPONbvIDoXXzx6xR4c9pGnueGNcjjJiQYaWSsuwfzQm+WhJTrD36SiL1/qjwt6NhcvYm54t22z/wEpCTBaXuAzLwAAAABJRU5ErkJggg==\n",
       "$\\displaystyle 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": "iVBORw0KGgoAAAANSUhEUgAAANAAAAAZCAYAAABAQ6AIAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAGKUlEQVR4Ae2cj3HUOBSHswwFwF0FQAeB6yB0ACUAHcBQQSZ0AHQAdEBSwRE6SKjgIB3kvs+RPbbX9tqyrfUmqxkhS5aefnp/9J7kDavr6+uDmLRarU4Y94D8mHxJfgetK8p92nMgmgO7plerGAMKi/zIWA3ngPpXisfUn0Zzbj/wznNgF/XqXqTUXtTGHVM/hAF6o33acyCWAzunV7EGJIPKxpKHbuW2WCbux91tDpR1aPF6dT9GVoRqT2rjDkP9R619X91zoDcHdlGvxnigMmO8UBh0iWC4R35dJrJ/noYD8PXFLQmnB+vVNBzsTyXqEqFMHkFlt3HsHm/K7V3PjPH27itjnnf1279b5wC8O6f1I/kL/LuirveX99+pf6PMEu1e7LipZRc9oXlnCvAP1qttLG6UB2KRepAHCKm38YRFnlEOHbMN/ixxTs8IGtAf+O83CA3qqmw81E2vyBrRzqURejXrWsGlZ5fvOoAsRZ2BHAmRI4onufGE+m/qPzPKLf/QT6P7Qb8kOyPzebPjFfuHFkhJmyfAI99OXRP5X/K3Jl7Spnc6Jb8ds/YJ8AKxf2K+KL3qP8OwnuDRWD6T5bs8L4yH54MoA4KoYYMu9jgwWFp6lJc+bEjveL8P3TYwqeP1JQYhD/ukYzr9Ii9i89gEeKRebSIf9d6NiIGZXoPPzb9y1R5lQBAxBNMSKyFCmIzm5gQALfiAfkm8TzOKu9OqPOD5JfmIZ73W0lOUXm1zUYUBwWQNIr8V+4dndzkVXm/zN/k/hJDtZJQPqcckrbdVkEMwxEw+dEwJj9f2F/n6pcM7eeP54xHt7lLJEnPnclJmYvNXIW2hs/x2B23lO+9mSyUeOsdcejUb/k2EMwMKi3yPELLQIAjoO4MVzAfqFzyrJGNDAUO3iteinqWEGPIp+5QZT8B2RGdxl9fvZuAFSlLjcU5ydgNHmRvyBRifg6XJSJTdVi5sUsqUuTxSKKchybP4KN7kHug9sxov58ldzR32U2hwkinCLmn+DjTrRSoM9Xkb6whErB7STe7g9fW7GTQprP1nSwi8cn6kbogmDm/m6h+4xSG/NbptpGQyhQ99z4WT8iH7DuROAYBiJ6Wu9/GsUhHW2Jmh+wcaL6G7pnhjMTC+bQdSef4i1w3A5Xgg1zjWUhlPwH1M38ID0eYVst9ZirYykanxlGnXn5lL4zGs81a0sk7eGYKf076qjyvX58ALzSR6VV7HnM+sRx7L64fw88ZeeNBQKpkOKofXn5X2sXVoakAeaDfSnQoDdAy1otci3oBFZchwU1cp5dFh3ta3ZEw0HsYaRmoMFf7R5uYhnjXe0pZhrY/pW2d8NN76HAFjtCzq9FLXwa8ByedCF/IQjrabFHYsK2teInQZU/QKJ2bGMBS/Huonwio8NHWN6oq2toP70Dn69nfeiocJA/WypqbfIrZ535sRif6dW6bQb4tAulbYGoF0DSq/ywyIyd1lToOSKCR3uEI5eP+Wpk81JbLb0KTwFehaSohhbe4NDZ6F6kq7lfMPOLw8aDr0KrO6kefLEn/Z+PP22cuUMoUvWzkD3WORMt/QIDMcysq5h/eZACYwHgWmUa790V1iDOIYkirGE7DKq+ycOITQBH1PmN8YvEjUDSvclBrPcrTL7zmiCci2p4XLtB1495vc0xe97vOk2/dHiP46Wk+jIF4HQfmdQ2+U38ZZHZOyq/EGAikxNEzf2eTO9jnwxo75TVdypUQO3rjlRqRXcXMzLO76FvWM96/IqdOSZTqIF8EW3KRyJ3NGm+s7H/1r7EFI6MzEfpfwm0VlZx9Kp29/5jM8ney3cIGZz8C/5kn7YJoaT9eczOWO6aVDbvRd3RvfpcTbCGDhjdswIEOOpwi1KZafnF1BifSig88BjDW0Lf6vh0DLm8S2j5Yb8Y/Bs5F4rQNzebAufkFSe92rmhJvL0AL65TcgFw/QjGUe4NSJ/FCsTwHp97yBJxZCBtw+3c3jd9+YueZY1xQ/DOwRnnKOTDdRprbMiBDi8X/QR1KaMx7GARvGNT1m7PQbRlFMPbFb1LL4FY8iq0YkHARsAdgP/xNdUERz4VbNhLeGib7WWLRHv42sP1//R3DyYnQa4YAAAAASUVORK5CYII=\n",
       "$\\displaystyle 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": "iVBORw0KGgoAAAANSUhEUgAAAH8AAAAyCAYAAAB1ewShAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAE60lEQVR4Ae2di1EbMRCGbYYCHNKB6YCQCkI6gKQCoAMYKmBIB4QKMqQDoAOgg9AB4A6c/zskj7g522efdEa+3RlFOj1X+0ur1QOn3+v1XuUGciH9Go/Hp2GEhfOVQL/fvxD3J6UejPqKAPxDuccg8UXgj4JvC2YsAYHP5N4KurCj8NWmi3gS2E9BogXXSAJuIk8msxsMvY016uPad0Wg7cu9evCadtjP/Kb1WPlEEnBAX6l6NPNQrmyfLd2ygb+06Nop6FT2Aa1pIBzJ24/Vsqn9WJLMsB4DP0PQYrFs4MeSZIb1GPgZghaLZQM/liQzrMfAzxC0WCwb+LEkmWE9Bn6GoMVi2cCPJcl26ol2uge7dsLXDmiNWtHJ3qUq4FZuz1V0p7h7hR90AvjbxS3sdRp8CZB77mcntc/y/0iY4dX2wgJNUUA8Haeot5NqX6AP5B4k0BsJlocrvxRGpXKB0hnqJPhC91ruXqDfBkgPFX4JvmsHNZC4ai2/lKldflUZO6f2BRLrJm47FLoGwvfwuwvhLs583iaOBHbnXy4ln/maaayl3ENDX+UQPiqWd2QYWc9uzVWwFWLWP4ov2v8ph8GHFrguLQOKWm9KCr4D/kxCLV4C65tBcCN3CeD6/qfwSA6DKzk5fnw7u54vIuBF7lRxf32GdfdTq/0zCfA8ECIzjFnv96ZsYYpXKkGeNoI7Atnz4NvD+LsqDRCftpb+zJnvBHG3YM8PJVi/Vz5XmJntCVV76+Pkh9a2z/POj8DDpD7aVX18V633bP3QTLtylXypLOcCLBtlYmnbUjrLSJl4GT13gKvsuFwwxrfaLjpcVddM8BGWCn2pKlgnzpUPsyK4YgkII2eFm/JQUTfAhwOynAXNVEnipZJ3Ace7uqHSl16+VHYqSJXMRIhMrfYnLEpAzHqocla9JbXyL1ppKsBKq9IKrTDWdiNJwWdGyKESoUJdaoT7JQEj6yRIf8uV/l9sEE74PF++RTQcW8BVD07PTy1f/biQu5S7cX65X1PrSQa+GAFsTtIK0OW/O0RROrMPYc9SwcoSl9zgw6LHGC1IvCCwH3L82Vo2JL6xQdg5Hcsh3y252jbaZsKecuuEkIdikqNPjJ4jhbmhwrjqieGyxU10clK7BwjO8cLgYyB+U/xEKyVnIk4D2BrI0xNa7UH9wv6Yu3wlA1+NI9Sylbu0QeR7F8sXf5XGW6z6W6yHgeuB9lo0jJvKSjLwp7a4ngkrsRM0gN/dT0i03qhG684lA3+uiOZncFpufsb0ObABOKX0GmBmixszUy0xGwlgw4hZDtBqL60GfjbwTmdUwHMyORDwHJfXJlP7tUXVTkYBOVBLtW9BlX9P+bc98O6bX1aptXPhZ1m46GDrZW6FMhAOAH/hcVCYQcDN5wlxLsyjzQInfWPgsW1my+cdt6ZogalYKp1yrzbzJYUPRHVuQf22DrY50GHAcJg2IQFfy+Az8Cci+xCBhW5BBfKnJlybwddEepHLVsxY1nPUeBIy8JOItXmlMtz8gU2yAyQDvzlO0WoQ4K3eghr40aBrVpHborV6C2oGXzPMYpZu/RbUwI8JX4O6nLHX6i2oqf0GgOVe1MDPHcEG/Bv4DYSXe1G/5vPUKjwStJ9czx3ZgH9hyxHwVhA1JOx/b5/EkOw/WwilkXlY4HPXX/4T8tF/Lpu0ZsotEIsAAAAASUVORK5CYII=\n",
       "$\\displaystyle \\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": "iVBORw0KGgoAAAANSUhEUgAAAKoAAAAZCAYAAAChKLVZAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAFdklEQVR4Ae2bgVHcMBBFOYYCSFJBoAMgHUAHUEKgA5hUwEAHJB2EdABUEKADSAWE64D8p0ge2ZaNvfb57sCaMbJW0mr3a7Urycfk5eVlxZImk8mZ+q3r2dDzqOdEvKbKxzQiYEagyq4mFkP1zC7UFwNdUflS2YbK22YJx47vHoE6u1o1orNf6Heq8pYGwruOaUTAikClXVkNFUFiowwhP6ZZhR37vW8EYhvK7GrNgolC/Gah35Yv3xboY3FEoDECdXbVxaPGAnCwanWYYpug5zBmMr73g4Bw3X8j27DMrkyHqRhOAeJO/1oNRzG97l19uC24VJ+9unZjXRkBYXcn6oWen8JvqjLRDOyvVP6l3CXROeDiPNyB15OXJpP8ObvqZKhihkfcFhiNjRSk1A+wD5YVRHSYVxJ2zxqbhR6nc2F5EhPUjjY3oi/dTUzKrsyhX8x2BcRmMFLKesJeNcYs9+6FuB3KSDUeYfA4J8QcCz3Ig4c814P3xDiZg5yRirYiGgeR66669yAv4jROGi9pV6bDlJhhkLjmUxTxUuBVD/x7XQaoY8ivQ6i+7jFlmBVdTkX/owfDXvhUZ1cmQ5XGN3rcPjPW3q/imJR7lyDu6kHtlnLflFNmCQrMhzB/1LOr9+slELnSrjJDlTIYXjiFf9E7ng/Dwnt+0vMkZd3KVP5BZUvC+1YC1kYGy+Bt+0TycB33EPSHj+rAhr32Z9EJs4MljR3miTlDNr4S3lcIAN5EukrcK/r1Qo4whJ/drqQgDFD4jHdfBogHPce+zPtdqLfm4nGl5zDVX/SZyCC+LA6nR2rcOpr6OUyUs296jtvCU49I/zFrmquPWR7GUALD9TCe3lkwVOA13fzFuejMpXnu1Ncsr/r2NqfBo34TU/YzIbFKAeC7J7D/7CNcw/Ov51nMhpKhOG6y7D3mb1+JRyrqzz57cC8lI8zt71UmtCMHV1bFDzEiObwxmHmk3ubUXU/hnqXwNGiiMquW1ZkDJdRbc/HlaoVrqdIEd5VB/Tnc4fmKiUn6qKdoaLTjYIIRllIsj5f7VG2zQ4loeDHuKTNazKRveWLexXeNhZHiObkByOmpOrZueNRJsV9cnoW84tmfXUmBVLhgEkzhMsUv0MQTQ02GqNAm5GrXiwziYw5dyIK8XpY43DL5VG4FeZvm6mOWR325xC+FcdFYpMhTwlY0J2tT+Yrt1N8sb4KXeU5XJUgu+RUIreT1cg1tBcL+q2FoxjK0lRyPey/Qs4ijMsY7Fa3qANN2jKbtUxGDvgHT1G8tqqJJ0zF7add1Tp2higmX4kFZB0Y8Cao7juq7CE5YArhSGlCG0tivENhX58KpynPZn2pcPpumvjQxZ8XFJJJLyB8vMk+efdbnnK6KGUoSUsJqze1LVe8ULXgUq5Z4oBLQA8vQVvackXpZwcrt49sy69j+TOOzH82SyuxNWfzJvbbo4D2L6Ci21anvOV3TUIQLPsfxayauXFD40APCPSF72HD6p9glMbk5oD2zIWVoKz/3yT88NvQNJ+vBJ1/zwAk/GCteEifCdqruLndH9V/1DJ16ndNOP0qxaC6guZPdA3RL/7Z9NB6HAf5NJnk6N/Bjoe2IXykyNOHVtzx1Y2ostnMcvsLiqmuerBtS3qQAnjgPQyVUtf7FVZ0SdXV+sogKrfdp6suWKPtfMM+LmwsWmsmjdpGnTs9UncbiNuBJspoX6ZDypnQItMENlYGlPFuAIwE4iFcNyrbNJSfen69Tbuvj5eZ3n+aJbyuDtb03sKX8mV9K53kZKiFp4X84rcnm0MQ9JInwWfdN3TValD9+US28M2iK11wMFeEEJAcBLqj7Oqg11fnNtxO2bK+uFz1itZmIf3yMaza5icDpAAAAAElFTkSuQmCC\n",
       "$\\displaystyle 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.40.1 (20161225.0304)\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=\"#ffffff\" stroke=\"transparent\" 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\">\n",
       "<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=\"#000000\" 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\" fill=\"#000000\">Add</text>\n",
       "</g>\n",
       "<!-- Pow(Symbol(x), Integer(2))_(0,) -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>Pow(Symbol(x), Integer(2))_(0,)</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" 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\" fill=\"#000000\">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\">\n",
       "<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=\"#000000\" d=\"M126.2854,-216.5708C122.0403,-208.0807 116.8464,-197.6929 112.1337,-188.2674\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"115.237,-186.6477 107.6343,-179.2687 108.976,-189.7782 115.237,-186.6477\"/>\n",
       "</g>\n",
       "<!-- Mul(Pow(Symbol(x), Integer(2)), Add(Integer(5), Symbol(x), Symbol(y)))_(1,) -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>Mul(Pow(Symbol(x), Integer(2)), Add(Integer(5), Symbol(x), Symbol(y)))_(1,)</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" 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\" fill=\"#000000\">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\">\n",
       "<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=\"#000000\" d=\"M143.7146,-216.5708C147.9597,-208.0807 153.1536,-197.6929 157.8663,-188.2674\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"161.024,-189.7782 162.3657,-179.2687 154.763,-186.6477 161.024,-189.7782\"/>\n",
       "</g>\n",
       "<!-- Symbol(x)_(0, 0) -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>Symbol(x)_(0, 0)</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" 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\" fill=\"#000000\">x</text>\n",
       "</g>\n",
       "<!-- Pow(Symbol(x), Integer(2))_(0,)&#45;&gt;Symbol(x)_(0, 0) -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>Pow(Symbol(x), Integer(2))_(0,)&#45;&gt;Symbol(x)_(0, 0)</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M83.7307,-146.7307C73.803,-136.803 60.6847,-123.6847 49.5637,-112.5637\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"51.7933,-109.8436 42.2473,-105.2473 46.8436,-114.7933 51.7933,-109.8436\"/>\n",
       "</g>\n",
       "<!-- Integer(2)_(0, 1) -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>Integer(2)_(0, 1)</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" 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\" fill=\"#000000\">2</text>\n",
       "</g>\n",
       "<!-- Pow(Symbol(x), Integer(2))_(0,)&#45;&gt;Integer(2)_(0, 1) -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>Pow(Symbol(x), Integer(2))_(0,)&#45;&gt;Integer(2)_(0, 1)</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M99,-143.8314C99,-136.131 99,-126.9743 99,-118.4166\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"102.5001,-118.4132 99,-108.4133 95.5001,-118.4133 102.5001,-118.4132\"/>\n",
       "</g>\n",
       "<!-- Pow(Symbol(x), Integer(2))_(1, 0) -->\n",
       "<g id=\"node6\" class=\"node\">\n",
       "<title>Pow(Symbol(x), Integer(2))_(1, 0)</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" 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\" fill=\"#000000\">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\">\n",
       "<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=\"#000000\" d=\"M171,-143.8314C171,-136.131 171,-126.9743 171,-118.4166\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"174.5001,-118.4132 171,-108.4133 167.5001,-118.4133 174.5001,-118.4132\"/>\n",
       "</g>\n",
       "<!-- Add(Integer(5), Symbol(x), Symbol(y))_(1, 1) -->\n",
       "<g id=\"node9\" class=\"node\">\n",
       "<title>Add(Integer(5), Symbol(x), Symbol(y))_(1, 1)</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" 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\" fill=\"#000000\">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\">\n",
       "<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=\"#000000\" d=\"M190.3082,-149.1278C207.3555,-137.763 232.4019,-121.0654 251.5344,-108.3104\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"253.4799,-111.2199 259.8589,-102.7607 249.5969,-105.3956 253.4799,-111.2199\"/>\n",
       "</g>\n",
       "<!-- Symbol(x)_(1, 0, 0) -->\n",
       "<g id=\"node7\" class=\"node\">\n",
       "<title>Symbol(x)_(1, 0, 0)</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" 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\" fill=\"#000000\">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\">\n",
       "<title>Pow(Symbol(x), Integer(2))_(1, 0)&#45;&gt;Symbol(x)_(1, 0, 0)</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M155.7307,-74.7307C145.803,-64.803 132.6847,-51.6847 121.5637,-40.5637\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"123.7933,-37.8436 114.2473,-33.2473 118.8436,-42.7933 123.7933,-37.8436\"/>\n",
       "</g>\n",
       "<!-- Integer(2)_(1, 0, 1) -->\n",
       "<g id=\"node8\" class=\"node\">\n",
       "<title>Integer(2)_(1, 0, 1)</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" 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\" fill=\"#000000\">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\">\n",
       "<title>Pow(Symbol(x), Integer(2))_(1, 0)&#45;&gt;Integer(2)_(1, 0, 1)</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M171,-71.8314C171,-64.131 171,-54.9743 171,-46.4166\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"174.5001,-46.4132 171,-36.4133 167.5001,-46.4133 174.5001,-46.4132\"/>\n",
       "</g>\n",
       "<!-- Integer(5)_(1, 1, 0) -->\n",
       "<g id=\"node10\" class=\"node\">\n",
       "<title>Integer(5)_(1, 1, 0)</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" 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\" fill=\"#000000\">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\">\n",
       "<title>Add(Integer(5), Symbol(x), Symbol(y))_(1, 1)&#45;&gt;Integer(5)_(1, 1, 0)</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M270.2854,-72.5708C266.0403,-64.0807 260.8464,-53.6929 256.1337,-44.2674\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"259.237,-42.6477 251.6343,-35.2687 252.976,-45.7782 259.237,-42.6477\"/>\n",
       "</g>\n",
       "<!-- Symbol(x)_(1, 1, 1) -->\n",
       "<g id=\"node11\" class=\"node\">\n",
       "<title>Symbol(x)_(1, 1, 1)</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" 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\" fill=\"#000000\">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\">\n",
       "<title>Add(Integer(5), Symbol(x), Symbol(y))_(1, 1)&#45;&gt;Symbol(x)_(1, 1, 1)</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M287.7146,-72.5708C291.9597,-64.0807 297.1536,-53.6929 301.8663,-44.2674\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"305.024,-45.7782 306.3657,-35.2687 298.763,-42.6477 305.024,-45.7782\"/>\n",
       "</g>\n",
       "<!-- Symbol(y)_(1, 1, 2) -->\n",
       "<g id=\"node12\" class=\"node\">\n",
       "<title>Symbol(y)_(1, 1, 2)</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" 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\" fill=\"#000000\">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\">\n",
       "<title>Add(Integer(5), Symbol(x), Symbol(y))_(1, 1)&#45;&gt;Symbol(y)_(1, 1, 2)</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M298.3082,-77.1278C315.3555,-65.763 340.4019,-49.0654 359.5344,-36.3104\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"361.4799,-39.2199 367.8589,-30.7607 357.5969,-33.3956 361.4799,-39.2199\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.files.Source at 0x7fc7dc51b2e8>"
      ]
     },
     "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": "iVBORw0KGgoAAAANSUhEUgAAALAAAAAcCAYAAADBaTXLAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAG2klEQVR4Ae2bj3HVOBDG85gUELgKLnQQoIIjHYTr4KCDZKggEzrgqIAjHQAVHKQDchUE0kHu+ymSkfRk2ZZl+5G8nVFs/dv9tFqtVvLL6ubmZqc2rVarM/HcU9pXulQ6kZxrPRejTcS0mDLukOBVbQO2hvJWfDHcHeU/6LGv/JOl9LaJmJbSxV2T+yAekCb7QOl5XD4gfxS1PVUennjjpWgTMS2li19SruznOAU8MGCMV41ey1t+SjUeUOYbqwsd/LIBrKo19eVvCqZqg7sHjC5kn+zmAe0GuZ2dd8r/EZUNysr4H0cdWBTQl9vH/H83EdP8Wvi1JeJUZcCHSi/1/rcbTeOBVfFWhe9V6byTazP2yYFu8UNcNIhBmKSbfRQX8dhmK2lAuj1Cx13sZJsnanOitntNWxXyTmde9sjXSuKHoXCgq8ZzLK+hmNCJ0sexcu9rf+nuqxKL39iWnuzIOMsjXyfKu8N+1lbUjli4sSnTWAXVDU08Ad0I8sEu9V6CSX2YAG5Rsord1qf1I939UMI5+uks1pfqcRRf4/I4rzbG2bpyF0JgbB+VqpBcPLcYjyXkFQzJK7lYuIqMoUxKMKkPevmicZgrwaEyS9pLJttp8sRdwm9snwp40N0bpXMlQgDsgmdAKiN0Jc7Njt3OxSW4YLCrFywa6x978wA/jBVDxaOfOiF6x5BfKC1CIzCh6MNFQN8doZcpg20ZHleu/ylh8DnCVrGpc24h8JbXdgXodTR9FgcWRHDlUZF/CcDBmOzCJmyYzfuWDOwu9cFGpHe863O95xzqN437T8aOAfOFLDtJYohBsp1Cz5TwTHhuvO1vSlcSaFaNng+Vn5xmwMQWlVMiu01vvUyuEAnw8HCV+c3NCbJVx3wRz/+u8mvK5iLJdraDvsDG2eiiRT46Z7fO6Z6+e3a8JvZtPWXTUKkJuvUOGFbAsUDoYd47g2/a1kqSOTkmyeBMwJ1jEvdUGMSXhWN02ya7rVz9zDzpya76w28HTyUVpcfTVq4+xXjgKUKPze2W3llIVOBl13Srcuwra0+qx3HC4wAP/Egp54Ffq57YxBErCBDuMplYJNff9av5nAMTY/yeAT0Hhoz4sMp62H9tKR4snhNi+ZxXCxlWyslIgzOE8oQI4OAqLf7ohVR0joPK0bWtfLTSC970kxhjiGuEm1ad68BWZG4rYmBrHScsmAOTZHD980LjTE76WAzqz0EXTxkTk9fmVDgQJQ/DPh6L/VRtm8OQyvBYfFBqynzBtfH4vON3ycJ48bTcSAQLTXV4VzwwtpkktUFHzM8hHrgxzlRrMYrrUfraNUiq71Rlm4BpLAb1T+pQk8OWzb1z0tDadOrwqD/zwwS7HRKng1FAycVIRW088JRc93Ei/iWisyl2ucCA6deDWOCGHugvLrspuC1O/+2jiHTP6UonxNRnKzMDmxBDieLw0BfOoC0DjJqbposShiP6IDdFLDAo9fsYbLHLqF3/SwyY1eAK9BqSJoeLdVdvAPmKUN2xVx92nig3EyaU2LqwZ8JQosGUV1sk/hX4f2QrsfdlTNhRvMgoh8CPTebIzct3DJjA/2mqtSYJQWwDbiUFAbnqjbBotadYBWXqxxXIjRLXOoNIfSbBlACBt0opn61xLgwJWJ1FgfeyWMFb7UtrJ4KfDc4kn3i3IeWJfTHAZCyvcnTeGuqoDmrsblcZJsrcqSUMERfPJ0B+jcU1DEJfWlDG+NSnibVU14uQIx4ouk1ujs8kmBICmfBA+V6buTB4Inu/Elu/s/NFJ3fS7zKK3gL6NtQ8c+PgjBiviuERmuXuonGmfynlCCNnDkzwzoMTXfJeTiDW7upqlUkmB5bmjrAW31p8hI0bmll/yGN1UnQPnBq3+LEIs/eqqX6urDYexzf1lCxCVT7AZG2O8SgZHZn/ibMelSA/eTJW40mI1Tm3zCEDET62uyfCmLxiHMKrb1vJNOcNyeyKA9dYqm9w6re8zHWT+BV54DF41gB2FEgWV4tXwtp6A+ON6SE6IgaG3L3cbW6GvxbI1QyiikVIQYRHhE9sfbMQk0IqFHagfn7Yg0Fz91tkvGAYiQcWvcjaA1FAq/FaRvwG4tzpqPmvZDFgsO9VQcw7OUle8K8hkwssFGAV+0F6CQ6whewm7SasHNYwYqjrNwe3rTbkr7Bz5nglPQeH0Bie2hHW8YHJXAn6BszW9VkVyZN3zOg+5aU0PDDeYfCB9T7pqXSs0i+hGl+Du4zXfb1rQt3GgBEuRhyqDsVotpgPuVvaaqBLA9aJ8Cu2YCcMDBgm1oiJw4rjpi4w2/qtBoZqQHZJbL/2z8FrBjyU8bb9VgNLauB//cXFCFN3t9QAAAAASUVORK5CYII=\n",
       "$\\displaystyle \\left( x^{2}, \\  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.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": "iVBORw0KGgoAAAANSUhEUgAAACkAAAAeCAYAAACxHzfjAAAACXBIWXMAAA7EAAAOxAGVKw4bAAACvklEQVRYCe2Y7XETMRCGzxkX4EkJpgMDHTgdJJQQOvAMJaSEQAVM6CCUQOggLgHSgXkecfLodLpzdEDuDzuz0Wm1H69Wq1WS5nA4NFO4aZoNfD/FttZmSaAqWiwWgnvfGq2rjCcqL9zVFALsFrtb7F9Nsa+xOatRnkv3P8i/lfnBi0PNrQhyAz/C1t0N9bdnfHEaO+5PoHmCBXYNe6tnoWIm2yxegugC/uZIFr/OgpCgRZDIbS82+QgsjoojWQ4vQkPH/Zboxfojy2v4lnXr1e87ePcv0XaaOcGsvdfwO/gHHDJIRuMLg2gGAkDv7QaGz9CutDaHrHfcZHPd5qpUhzOksWl6IEERWg0Z+z6GKNnMmNqkNXx32l0J5OCliRFx0mvsOobvo86pEd2dfuKY6XshQ4cJ8rzGEBroLpfHuWvwOpm7a2+7/BjlYyN6XlA3Gu4D3/bkTkzm+gtxSpfmJ4vFS4NcQEfnMYgj5M6fC9KndpvZ62QVZX7DAXjnuEmxC/JQPX5g7TM8mdoYXk5bXEpPTI5HDFjnHvu6AxLhG61QGLrZG9aGNqDpc0iAkiBSEvR5KuDbhGxzkL7VRYDuKHPwp9M8k/rzFFPy1bs4I/gl7GWQTHf8DoLkhyCLT2WiU/OZZy2f68tsr8ykb/AKoGEXHOdH5kOUH9GQ3pg8bjTPmvO8lEIJLFkIINvxasS7BrnjEfXykheChAi0lDl/LUxJnf3yROZSAx2P1WVxAwDSxp53JcDWYZjzHeofHfvkl2S9VQvx9r0+ieKgDEv727GXqQtFEK4psKaPfZbvDWzvvUx9qwN7irJ/GvfiKoe3vYWScpRhoONOsLh2apxih82DfqtABoMJ/1oRIHx8Sk9tqI2jzbXf3u5a8pcC394aOidYvNUn7fBvfft31e9OI9JaxkF1Zmpi4N+yOtb+L70YupAa3+UnAAAAAElFTkSuQmCC\n",
       "$\\displaystyle {{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.fields(\"img(4): [2D]\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAywAAAAsCAYAAACHZjwAAAAACXBIWXMAAA7EAAAOxAGVKw4bAAARw0lEQVR4Ae1di7HVOBbkUS8AYCIAMoCZCGbIgE8EO2TAFBFMMRkAEWxBBjARPCAD2Ah2eBmw3bo6Xln+Stfyt1Xla1ufI6lPW9K5kuyLHz9+3JATAkJACAgBISAEhIAQEAJCQAgsicDFxcUr5H8Lxz0c33D8AVvl+kIGC6CQEwJCQAgIASEgBISAEBACQmAxBLyx8hq2CQ2VG7h/h9M93D+UwbKYWpSxEBACQkAICAEhIASEgBAQAkQABspXnB4FBssD3H/Gcf8mfuSEgBAQAkJACAgBISAEhIAQEAJLI8ClYOau/cW9S/PRWQgIASEgBISAEBACQkAICAEhsAQCmFm5H+XLGRa6T5phOQGhXyEgBISAEBACQkAICAEhIATWgwA34LtN970GC9aS/c5jPeVWSYSAEBACQkAICAEhIASEgBBYKwKwHW7h4Ib5bIf0NFY+YtblLwrpXBKGiC8Q/ow78xlRTggIASEgBISAEBACQkAICAEh0IcAbIdr2BGvcXzOsSOQjpMlt5D2ueXT+pYwRPwNEWgZ3WWmFllnISAEhIAQEAJCQAgIASEgBITAEAKwJzj58QtsiSdDcS3c2yB8U9gf9PP3/zQMFgTwYy3/wfErIn9hZDkhIASEgBAQAkJACAgBISAEhEAKArAr+FpiflvlzVA6xOUm+7c4/gzicpblSZvB4tacpVhDgVBdCgEhIASEgBAQAkJACAgBISAEOEPC1xTTaBlctYW43xGPEyc1B5vkomaweMvGfaAFge4rk7UUO7xBnbmph+AQUNbZvY0AZzkhkISAuJQE12EiixeHUXVSRcWLJLgOE1m82L+qj6hj1PnsyZDYYKGx8u0osyueNJymcsaZB/Qe7vWigf23GZPWUFyaFM7dCBMvdqPKSSsiXkwK526EiRe7UWVnRY6qY9Tbvlh/G2PsrL3x1WuNIYwzDBToNrl0or2vgMdRdbhm7oHHIgrSrRDoRUBc6oXnsIHixWFV31tx8aIXnsMGihf7V/0hdQwjhXvieWR/KqUyWCCEhsoXm23A9VEcDTVzZvWFfhamsxAYQiDkjbg0hNZxwsWL4+g6pabiRQpax4krXuxf10fV8Wuo9mWuesPvsDyFkHBXfq7MzaSDcXY/KixnmOg+nU76FQLjEBCXxuF0tFjixdE0Pq6+4sU4nI4WS7zYv8YPruOP0DA/KPkbcOB1krvJ2H4JFDeeJwtIys1HRn7vcHzISVs4DTfga9N9YZCnFr9SPolLUys6UZ54kQjYQaKLFwdRdGI1xYtEwDYYXTpeVmkwUrhfnKtPHuWUxBksSMgPRd6AsDm/u+LeGJBT6BJpQGQOMD8Cg79KyJfM4gishk/iUnFdp2QgXqSgdZy44sVxdJ1SU/EiBa1txpWOl9UbJ0aczZFaDPeWMAywuK6MUzTxEqlUeZuMj/pzE9BD1J8fp5ETAtkIiEvZ0O06oXhxvnqBITer8i2Ou/lTS
Loading full blame...