{ "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", "\n", "![Stencil](../img/pystencils_stencil_four_points_with_arrows.svg)\n", "\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.93 ms ± 40 µ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", "text/latex": [ "$\\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": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.figure(figsize=(3,3))\n", "ps.stencil.plot_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": [ "643 µs ± 8.66 µ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", "text/latex": [ "$\\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", "text/latex": [ "$\\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", "text/latex": [ "$\\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", "text/latex": [ "$\\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", "text/latex": [ "$\\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", "text/latex": [ "$\\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", "text/latex": [ "$\\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": [ "\n", "\n", "\n", "\n", "\n", "\n", "%3\n", "\n", "\n", "\n", "Add(Pow(Symbol(x), Integer(2)), Mul(Pow(Symbol(x), Integer(2)), Add(Integer(5), Symbol(x), Symbol(y))))_()\n", "\n", "Add\n", "\n", "\n", "\n", "Pow(Symbol(x), Integer(2))_(0,)\n", "\n", "Pow\n", "\n", "\n", "\n", "Add(Pow(Symbol(x), Integer(2)), Mul(Pow(Symbol(x), Integer(2)), Add(Integer(5), Symbol(x), Symbol(y))))_()->Pow(Symbol(x), Integer(2))_(0,)\n", "\n", "\n", "\n", "\n", "\n", "Mul(Pow(Symbol(x), Integer(2)), Add(Integer(5), Symbol(x), Symbol(y)))_(1,)\n", "\n", "Mul\n", "\n", "\n", "\n", "Add(Pow(Symbol(x), Integer(2)), Mul(Pow(Symbol(x), Integer(2)), Add(Integer(5), Symbol(x), Symbol(y))))_()->Mul(Pow(Symbol(x), Integer(2)), Add(Integer(5), Symbol(x), Symbol(y)))_(1,)\n", "\n", "\n", "\n", "\n", "\n", "Symbol(x)_(0, 0)\n", "\n", "x\n", "\n", "\n", "\n", "Pow(Symbol(x), Integer(2))_(0,)->Symbol(x)_(0, 0)\n", "\n", "\n", "\n", "\n", "\n", "Integer(2)_(0, 1)\n", "\n", "2\n", "\n", "\n", "\n", "Pow(Symbol(x), Integer(2))_(0,)->Integer(2)_(0, 1)\n", "\n", "\n", "\n", "\n", "\n", "Pow(Symbol(x), Integer(2))_(1, 0)\n", "\n", "Pow\n", "\n", "\n", "\n", "Mul(Pow(Symbol(x), Integer(2)), Add(Integer(5), Symbol(x), Symbol(y)))_(1,)->Pow(Symbol(x), Integer(2))_(1, 0)\n", "\n", "\n", "\n", "\n", "\n", "Add(Integer(5), Symbol(x), Symbol(y))_(1, 1)\n", "\n", "Add\n", "\n", "\n", "\n", "Mul(Pow(Symbol(x), Integer(2)), Add(Integer(5), Symbol(x), Symbol(y)))_(1,)->Add(Integer(5), Symbol(x), Symbol(y))_(1, 1)\n", "\n", "\n", "\n", "\n", "\n", "Symbol(x)_(1, 0, 0)\n", "\n", "x\n", "\n", "\n", "\n", "Pow(Symbol(x), Integer(2))_(1, 0)->Symbol(x)_(1, 0, 0)\n", "\n", "\n", "\n", "\n", "\n", "Integer(2)_(1, 0, 1)\n", "\n", "2\n", "\n", "\n", "\n", "Pow(Symbol(x), Integer(2))_(1, 0)->Integer(2)_(1, 0, 1)\n", "\n", "\n", "\n", "\n", "\n", "Integer(5)_(1, 1, 0)\n", "\n", "5\n", "\n", "\n", "\n", "Add(Integer(5), Symbol(x), Symbol(y))_(1, 1)->Integer(5)_(1, 1, 0)\n", "\n", "\n", "\n", "\n", "\n", "Symbol(x)_(1, 1, 1)\n", "\n", "x\n", "\n", "\n", "\n", "Add(Integer(5), Symbol(x), Symbol(y))_(1, 1)->Symbol(x)_(1, 1, 1)\n", "\n", "\n", "\n", "\n", "\n", "Symbol(y)_(1, 1, 2)\n", "\n", "y\n", "\n", "\n", "\n", "Add(Integer(5), Symbol(x), Symbol(y))_(1, 1)->Symbol(y)_(1, 1, 2)\n", "\n", "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "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", "text/latex": [ "$\\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", "text/latex": [ "$\\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/lTSbwQL9oQEC/aUNmX31F1jHq/gCZfoR2/SNVoOMNyiO+uxAABPFp6981Y4T0O28sSR9e9EOhEQFzqhObQAeLFodXfWXnxohOaQweIF/tX/8F1bJ8RCV88MErpZrAw4eEMFpCGhgmXgl3h+jEPXPNtaYfDAnWWOwMBcekM8HacVLzYsXLPqJp4cQZ4O04qXuxYub5q0nE1vk6eGAjfEva1JFWgJG7q59IrLjv7Gk7rI4yzHNyEQ8PJvgNj76r+BX5XjO8V/Qz3/8VBf370sfVFAYjL9JRp9aIR8hzxw80+f8OP5aqtaUQcbgrapRvQAzF7ifq7D2d6DN/C7+7aMBmoxxJ82jSXBvAUL/LbGfFiBS3pAL/VXiTqaABPtRc7aC+kYzfGLDH2XE2fsJCObULgTmKz46JzwP4Dx+8YlHLjfZED8rlmjRmyc/ge5hOE8c1hPF5YOK6tfFz3xn02rny4ZqNYkxOEcdbkQ3BPQ4iZ08ApUr+tyAUGpocGfsQMxzurC66pK+JW4c4wOFrGFb4Wf84z8rd6iE8TcDrAM4sXTI+De+HIIZ5vzckHywv5ihcT8KEFzyxeUA7cJO0F5LAMVd9gZRxzRjrxQrxo9P3iRfd4KMAm69lnehzqE1Y89jxXx7ntO/Jlx5Dcll8iUZaVg3SjHaw4GgxXPsETnM3C4iuV2ZlZGONdA4RwU6WV7yf4h7Mp9KcxU3OQ5zo1eN62AKT7Bn/ecjB1WBfpgTNNlR48KBz8Vy8eIN5IQ8xdPFxTVxZOXS3ionqIT2dqIcIzhxeuY+JzxqJAHmcs+S+Sm6mj3xwuqod4cSboEZ45vFB7AR3srf8RL9yDtetxyQQ6Vp+w8md/Ah2f277/lNNFMVNaO4/RsDb+gTA/hNM4+Jx4PGB6pg3k1CwrC+MZB8Pif/NpgNC/kuFlclDU+Jcfft9jf9xbHWsyrEx2Rjy+gaYTgy2Esa5d5URYVX9cE6fKwmU6HMS5hgHuG7NS8KNhw2V9o7BC3GzutOVBeeaP61XyCeWq4WjlXeqM8hTjBWR/DevLvHA0uBTXHXHEi5HPUIzdVPfUVZcs6sfCcD1ne8HBTltfQ56xHG1h1cywldnOiB/WQ+3FCM4BM/ECOJE7OHixyLgE+RbrRyC7mI4hW33CiT/WF1ZtkLVLdt6qjoPyJ40HmQ6OP42xpcnsOjOhAdprsHQJSPFHXq3GB2VYWCwP/pxSbAyO4VfrQL0Mq0tteRvicjlZQ0aYF8LdkoHQL+cacliGhiGVI6svTVc+8CfGtcY1lsNwHAS9agxxTYwaS+zg18AFfskEjcswxT3KsUo+tWGWU1/ISeIS4rvXBdo5zBN+xXgB2eycKs7h2pZgVn5hWUpf+7qS341OwsLiMsC/eDuDPBrPUlyOMfeQk8QLyuxKA/9ivLC6II9J2gtf1upPFpM/9uzTb4IXKGuyjsfiEMbrysdj1fv8Is4q+hFfVvECz3moW17DqU8IxpfAY/V9AsqY/Ox3pYH/Ktt3lIvkTDZYLpHoGgfdndOp6C+XE73Hg2R5hplxyUHbhyvZKL4PI2Iqi37858P5454kpMyffbxw6Ri9KDv281GrJSy22b/yT7lAGUiy4sulhvIBDu8R5zWOb7iOl3xZlbhc5ksU3sAIMjjwXPMyutXxCZhx5m92LiFfvtCCyyZd3rjnW+/4rzN1zY6sGC8gO/5+E58Fuk+n0+y/4oWHHBzobZdK8mJ2rQ9nuHpeDOlruIrjYgzlI15UOC4yLmHbjRJU/ciQvqrSBhdIoz7h/3js8tkf4sUan2OUmWN3ujY74BTS88vEtHay/6EAKA3LPvaDfMvHzeTgvvZvI+45YxL7WZraPz2Ix39DPzMPOA6qTSatSXjXy0M/i9MSxg69lm8cJ+UesmhM9c7mpMjrituXD8KIW98SidrmeuYBx3/JazjgvpUT8J+ljl119+U1bpju47LPzifgshiXkHdtlsNjRMVW/yjxGkcxXpi+fFlauWNxSp19HavnHfeH5oXhDBw6n1mPWTFe9OVt5Rtzhhy271m88nXcDC+mwmwI1758PGbiRbMNsb6nyLgEuHf2I336inWNuOoTTmMb09cmxgopOjad96VBWNF+vy9vK194RnxbgVFbCRXG6bq+iYBrCKDL2gBzSjrq181+ID/+08sC82Fyzt8T1PjffHaytEDi2RHKMj++qthmYJyfl8ektgGYlxaf16F7iZt/hx5bv/Y6vRfiENWJOFMHziEecb6Dgw2lc/DjdduMl4+x+GmNfFqES9AVnx3q859IK3y23TNE/zl4gbJw/8FH5BW+OIPZz+XEi0Sk5+BFYpFKRBcvElEVL9w4ZYlxydn9iPqEGtkP/eyv8DnmM0UXj1dOvj2/lz6MA5tqANsTPzsIoPGNUzRWOE3JwdObQBjz5hKm2KigfxjPkvyJi0deFmdbnKNi4Mcp3Fc4X+FMI4zLoiibdWxz3Hi25oF5W5nH+NEI42C1gR/qy2/a3MdB7MxwvIvrt/DjmlfixlmiRlr4r8KhbGvk01Jcsmc35jgbBBqioSvGC3CHzzZndDj9vogTL7JhL8aL7BJNmFC8yAZTvJh/XDJFP6I+wVNez74DYk3PsY1J4vHKcCMFZTIS/3F3S6x4v6fD1611cw/C+FBPunwL8mgkTCqzTR9D+SD8MY7O6fw2mWP9hvIeK2eL8VB3PisNPsFvMS4hb86G1ZZ/EVs4GqO1JTS4L8ILyCXvq+VX/r7zLTRb1H1fmVHf1fEiLK/XR2e7hPAivPA8nKRN9GWs8Tms4xqvc3kxpK+p6jqUj8d81f2IL6N4EYzdgIn6hACPqZ6XFDnQQVafgHTJ7eVQGoSvpn1HWfjHZmO8Mgbbm0hIxxkGs8idxxZ/8A8vH9LK4f4Wbqh8bmBrc6xz18b0tvhb8qP1yvqXcKXklihrtsxEPq2BS/bPhdU5vqf/5LzwOL2C7Ctcc7M/G0duGN3ls7VBXkAVg25yXgQ5TtVecAZ+tbO+4kWg8XGX4kUTp6n7kbgPiO9Zgsmfff8sqE9oH3tOreMmi5o+k+s4yCL1OeZLeq5hoLBMSe7Sx+byqRcgmb1tK0nIGiKj7Hw4WIfbARA0VN7gPl5qFha5ARpk0Mip3tARRo6unwR5RUHDt4XzaVsONFyonhgoLx804kJ8uEeG+F4Bg6X2KyD7Mg51y+HTUlwywyBuOHgfL3ecnBfI428czKv2x8A5zwZkrdJtjBcpGE7Oi6nbizXzaSpepCjM4iLvkv2VeGFAZ5wX5IX6hAx95SRZUMcpxV3Tc8z2qm9M3lmvSx9iiZ/ifrX/YHXW4hTA6TcOqJ+CQBw80Yrj0h3bkI/bhqMS40GebfI3TBqJpvJA2ZhHqXz4T4o1WpMUGeWlvMX2J0xSifFCUvm0GJc4kAPnqZu2f8/iVwuX4MXt8bBuPuZmeJGIdAleqL3o7n9a24tEnbno6kdyUJstzSTtRWpp1SekInZW/EV0nFjiVbTvfmz+AGWv9p6n1MMZLJ7cHDhzw/omDZbMRpsdKo2cKV3DAKJwKIr5UElnzcpQlnet+VggzsyP9ZPLQCCDT0tzyXELVXUGMPjGZVlt3zwSLzL4YEk2yAsrutoLQ6LAuQAvWvWlfqSA8gqKXJgX6hMK6tZEL6xjK0Zre2GBOK+l3+fsCl3WH/U3T2ndL5dzmLDAe7+XINo1a+etvrMqyo4EBxsILiVyy6VwzzdumSOhfsZxFsYj8rH8HuKClr/cDAgszSXkz2V5X8EPviGPHHwEP/fRyKj64kUESMnbpXmh9qKkdvNld/FihL7Uj+TDvvqUU/ICstQnrFDjU+p4RHthCKyl33cfYgUGWX+mXyChqxAqzoaQH9vjDEDfMioDYBdn1JtGBV97PEudkR8/XlQ8L+TDt76RpHIzIbAFLokXM5EhyEa8CMDQZYXAObxAWvUjFZL7upibF8hPY4WZKXRUHaPetDH+lTsGrgwW6suD+OxoA13U+wPqTMuvqEM+XKbzBXllWZdjC+fzuYN8Nrm8b2w91xhvzVwSL5ZjjHixHPZrzjmHF/45Vj+yZsWeWba5eOG5pLHCmfrKSX40HaO+nBx4iXFp9p7XcEkYN5tzCpHLmbgp5kiOy2jcBy0LV5oNQ2ljhTNlXA4kY6WwMjvEr5JL4Ld40aGwmbzFi5mA3lg2ObxQP7IxJWcUtzgv1CdkaGXaJEfT8UvAN+btu50o12ZYGMsP3J9jwHuo5USo9yyzH52amCgA9aAVy1c5X08kUmISEVgjl8SLRCUWiC5eFAB1ByLFix0osUAVSvNCfUIBpSWKPIqOPdfOXr3VMFiIN4RzszaXSe3u+xqJfFJ0ISAEhIAQEAJCQAgIASEgBBIRgD3BN5R9xvEQNsVZK4xqS8KCcvDtQs+R0dGWhgUQ6FIICAEhIASEgBAQAkJACAiBVARgQ3ApOidAuNH+LGOFebfOsLiAk1XEjLgf4uyMKFNOCAgBISAEhIAQEAJCQAgIgX0j4FdrvYMNMcme6k6DhTD6qZxbyOzLvmFV7YSAEBACQkAICAEhIASEgBA4FwE/u/Iz7Iesj0S25f8/C120EDRjZpUAAAAASUVORK5CYII=\n", "text/latex": [ "$\\displaystyle \\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": "iVBORw0KGgoAAAANSUhEUgAAAzUAAAAsCAYAAABc3ZYiAAAACXBIWXMAAA7EAAAOxAGVKw4bAAATH0lEQVR4Ae2djbXUthaFgUUBF1JBoAOSVJDQAUkqSOiALCrIIh0QKniLdJCkAgIdkFfBg9sBb3++ktF4/CePPLbHW2t5LOvnSN5nW9KxZM3tT58+3bIzAkbACBgBI2AEjIARMAJGwAisHYHbt2+/UB2vdDzQ8a+OX2TPXN+2USMo7IyAETACRsAIGAEjYASMgBFYNQLBoHkp+wVj5pauX+v0QNdf2ahZtepcOSNgBIyAETACRsAIGAEjYARAQEbMe50eJ0bNI12/1fHwjn7sjIARMAJGwAgYASNgBIyAETACW0CAZWfRXQfPg7sxxGcjYASMgBEwAkbACBgBI2AEjMBaEdAMzcNG3Zipwf3jmZobIPxrBIyAETACRsAIGAEjYASMwLYQYNOAaqOAXqNG69Z+5tjWvbm2RsAIGAEjYASMgBEwAkbACKwZAdkYVzr40H+SU14Mmr80e/MbAjqXnynhM8X/yG4CJLQzAkbACBgBI2AEjIARMAJGwAiUQEA2xrXsjZc63ubaG8rDpMuV8j2NdWnd/UwJv1MCLKcvKTAm9tkIGAEjYASMgBEwAkbACBgBI1AKAdkdTKR8I5vj+zEyg53CDmi/kD5cfzgyahTBn9n8V8e3SvyOxHZGwAgYASNgBIyAETACRsAIGIE5EJD9wbbM/P/M733ylY6NAV7p+DVJx2zN921GTbW2bay1lAi01wgYASNgBIyAETACRsAIGAEjkIWAjBW2acaw6V0lpnQflYYJmAMnu+X2gVETrJ/qD2wUWf1T50EOX0xGQNjyMRNKQGlgW+3UoLPdhSFgXV+YQjd2O+bfxhS2QHXNkQVAL1Ck9VYARIuoEFgrl1SvkyZWmkYNBs2/nqUpy/pAHqbUKkMxKO2Brr0JQ1moF5dmXS+ugl1XwPzbtfpH3bw5Mgqm1SWy3lanks1WaM1cUt1YWoYtck9j5Oxv+ustnSWIGQSEVR/d6GxXDoEnDVGsA3wUMG9E+XLjCFjXG1fgxqtv/m1cgWeovjlyBpBnKMJ6mwHUnYpcLZdkyPAtP8ekv5OpjRoJwJh5F2cT5LcriwBGY3TR+kzDYpzP20cg1at1vX19bu0OzL+taez89TVHzo95iRKttxIoWgYIrJlLL1W/51PUlP5PzQ8SkO4kMEWe87QgIEPxYSOYGTHcPzcn/14KAtb1pWhym/dh/m1Tb+estTlyTrTLlWW9lcNy75I2wKW/pCP+lPM71RX/aHeHlGEZFB+xZ2UeXUqSUGW90PFMx8/BHwf4Sap2r9K/DfmoK/VmCRd/2tOcSqsEKPy1jj/bpS0ayqYB3iigsApWqm/rurCec8SJE3tvb8y/HsKYHxU45kgPR4hy3zIA0Maj3Q7cWlUbIEOG789Z5fI4l1qVUaNM/NnmLQma9X9pRBw+/nmjcn7TwT7UzAxheKTTYArqdKRjWuqj8nzSGXnXkvWHzl2u2kmhK/Lc4Tw8KvMvMDh32TspbzX6tq6XZZzw33V7Y/7182/v/AAdc6SfI41Y9y0NQC7hcu/twIrbACZZKtskh2dx+Rm7cM26hbOA46MfdvyqDRD5rxXONY3FmJ3AqCM3inHzRscfktFZb8WN+mdSyTmLCxhcqV78SZBdYQTWpG/rurByM8Xtvb3ZA/90j8zQ06dkvyDaOz94nMyR8Y2K+5bxWG0p5d7bgZW3AYzxW1dh9XEsGjVYQ53GQZ+AjDgMjLblbVSc5WgM9uNH1V1i2W56k7uz6f7A+KHqXxk04fqDrmedHesC0uHzIWBdz4dthuTdtjfm3yiW7JYfoGOOjOLI6hJZb8VVstt2YANcqmwS1ZMXV6Ptk7j8jJmP0Zkm0opBfZuLhkxXfFueTYVJKXw3xLKzN/I/4ZAf42xuzFWE3TkRsK7PiXZvWV3tyUW3N+ZfLyfSyF3yAwDMkZQG2/Fbb7PoapftwEa4FMfHo7+7hyF3E5q8T/xFvQLwKgj80CMYw2rQSVbcuxqZ7CrGn1rWsx2hLNIQ915x9dIExUFgPjyirDjjE6e3vlFY9b2P0gHijzr+p4NwymibZaKDID8yI34o4qnSpx84/a0w6ssyu9opTRxg1WH25CGwQn1vWtcDeML15+JttVQ0cP+Vwr5cE5fDPUCkPbY3m+YfSpvb7ZwfwGuOjCDZQFu4xFhiNXobwMb9xPrHpavhUs+jGI2a+z1pWqMYbPPR/c8amLBZQPFDsjEiKONFU77CMCCIe9aMa14rDTuZsUytqqP8US7bvsWwqgzF0eh8jOGc5WIccjjqMuWPODyTP5XHA3ogJ8pUOLMvfybXsT4YQcVxtMxjTNEBuMhZ3wU4l+B5xHvF8cy8jjwMmAN++rzwTLCZB2k5189rzDf3WWXG59DtTQFOzK2vqfKlZzhat+Fj5Zgfx+3oWOy2lm4qR7hPOfctHe1Hgo37iZu+ru7nhE3sfw76xcApj1M6ONXWtghLHsSsNv6uMmRbQcoz1WE4dLkvuiJiuG46nf3AaPhXbwyYQWHw9FB+yMQ3OjjWSkZLL065xzjSXSt/PYuj64jDFwpPZ2UIP6q3yqo6VMXd01G5UB/8DOjsZkbA+i4LcANPnrX6+Qkl0SDXm1zwnCgPz0qVTv7KoOE5IL2umZnkjdCYTUDIUtodPbdJAW5vEjB26jU/dqr4odtutIUeSySANbBxP+FxacKOWbyDfXVaKkZNbNj7lmowQCEdA5Qc95MGOCwN65MdjQmWek1xDKC+Cw8aH97H3dVYghaXmCGXTQbehfvAqDkwkHT9iERybDOdOtKlRk6Me0W4ZF7HAMmOMtrSx2RgmfXhU51xRR7uNeh2dK0C9lM51FbO6vW9Nl0P6C3F8+DPeMkXFNDkNs9VNH4w9HnBEB3PEv8t1ct3xZ/StsSy0rPbm4DGEPYpaOfww6MJ7QbGMgZ108Gb+5LJUuGmg5cMRtuc+ZGgYo4kYHz2pm3hKscSc+pt4DlNsXE/8ZkzqW9T49IFuZRi1uWnnR/tMGpGOXUQDN4nvXElr0CjnGjApGXGCseBURpX+5Wft74Mjpp1oF64euCktAyucPwXTuVC/fFXnaOum4Oz6o1Dku4m4036A0NH8hngUW/qlDpk05nGOqVxlV95mdJOja2jNGMCQh2Q1TTOxmQfnaannAeKu6/ymzh2yg64NPXXmX4oIsjDSFylvlWvRXStcllCydsNXhQw85jyrVNvCZ7wGH7HFwTyVs8Ns5vN57TJdV4YxDQxLg1D1oEL5RblhTCgjF23N0vxD+BVNm1kW/vUyT/ytbkGf+skKoPnnnY/nXGv47s88M38uEGnjSM9uuuCdFJ4Tzmr4Ag3pTpupm/pwbNVP0rvfsL9BBw/Gqfkcik8K8Xa/FbCjgzEqIkDj7YBwEgxo5Ix+GWg1HSx3KHBcWUwNDPrOsr8J4ljiQz/YRPvLYmqZmiYPWo65KeDOJR9MLjTNesmkfl1yNysc9esTpVc+TGC0gFmEDP+FMjG/eEYLM7ihsoRDn8ozUsd6Zv6WeoyQujq9C1cFtG1yuWtYm3I6Jrd9vgOpnprPVJvpH2ntNE4QQVH3JZc+FcvtVR6NudIHY0cLn02b0Lm/911e4POBfHZ2xqVi84726eR/JufHTez77HvSMvbbX80pLsUpFP8Q+WsiCPc5ur7liE823SlPO4nboBxP5H0ExO5NEubr7rE9rltHN9G6zqMjNkf46jhyfoQXmXwEB19cK8wlha8HZKnNK0f3yv8fZpf/ng/T5AphxVa11XXH1vCYp76w66Ql6U0Vd10ZgAXZfL2Rkk+yw3pKbBK0xKH4g/q0kyTcy1ZGFzs7lbf2xz+vnIUB271h+NzlD8kM9Shxl3XBxjr+uz6VpmL6Vpl8zw0eQw+6YeMvXpTWgyVA73qGrlNbHs/4At5etMM6XdqvMrebXuje1+Mf1FfqkNn+6S4Xv5FGUNnyaEdnsQv5dstP8BVrpMjiuvU3ZBOcuL7ylHcGjhCHTbTt/Th2dSL0rqfuHkOdtsOiANF24A+/iku+3lWHsbcPH9Zm5jdEdmvlQnHcpXZnMphKdgHWV9xOpeZEG4UUv0UCyZMxycdb2NYOL9QWLpen/zkva8jXTtdzaKoPGYSAIWHt3LhmjLrt8shikYc46A584KsGMY2zXEmpwoL8ioR8selaDF9FZ78PJf/P8n15r3CA+6wVACcl3Jr1PciupYe4Da6aH4zgJ4qjqOkEXrj+ah1Krnk5TmjEaycwvC3zXjGeF5W8M1Z1tKgKnOBH5W75/ZmEf6NVdsI/o0VNTndzvkBbubIMHsusm9xP/FZ8TtvB87WBkxs8xnP4JrjmZvQjt+7IbwanHakKRnMuvnneqjigOkbXX+rG64HR9y84ln2ki59YSDGTmfRsIn15WYP/iND6diRCYMGg4d89Xc1uqRclks1DQ/C03S6rNyv+n0cZL0MYdWgUGEsx6E+b3TGIKzqLNnUrc1lfyDbJmSFYRhqDHrb8Ju9uivV91K6js9Vk4M8JxglqevUmzD9TbxmN0E4H18KfCn/K4XFddjMErbqXGl49pgZisuQ0nLP6d9re7MU/3J028m/HCEnpt0rP4DNHBkgzwX3Le4nDnW/13bg3G1AbpsfxyzN8cyh9ppXenAJ4s3s4BIw0vpoxyBg2LVEjkak6FIxycOQKCqzTbdD5Sj+iY6DpUptci4tTPfMM3Okb4UtpmuVzewJD3S91Azc5TBMDpbp6HoWvUkuvKyXqYVrGk+3HQUxEK6r41+q46D3zvZJ8SfzL8g44HVahz37u/gBJnK9bZTiL6Zv0b3As01xRPWd9GyP1ZvSuZ8o2BavuZ2Zm0vpvQ/xT/FZbb7S83L0aDyTltnmv6NMOGZKaOjsRiCgN9E0CrXT9ZUu6AjiErQ6LnjA9mDmqZlgw9dY0dz/xbpMfa9B1/ENR9RJ85rw4noLOL2Q7Dfys0EBjRgfq18q93Vr87sN8m8MKCX4x4x762zhmApcSppMfnDba2ijxsB/8RzJ1F1pvTX7heY1OiqhgwNdh3t2P3GAyukXC3NpzA3kcomNh65luJBvtLsbUrKE6plA4Q1vloDRJV1IQmHEwwhW9xKsMGZ+13VzWVt610e4SgaG0Jgdir5PykpljvLPXE7b0qZR9dpCoon6XkrX0Xi4amDLNS8uUjeH3v5WAZR1YNyfwt20wnv0b4x/OSo6mX/mVbXlcLH+KEd5Me3a+5Y1c6TUsx11kXF2P5EB1haSLsilHHhy23zGx31j6tayo1ETM/JHSrt/89WK1OdApoZ5Y/KDiMQADmuSZUhxEwFdHjmU2RxoshwH3CP2R5lKBcxcDm93YiNZqsprkpOr78V0TQcuTqKLtjduzW2Vi+tN5d9bk+IupC6b4V8m3sX5l1n+pSTP5Qf33dpGTQHEfcsU1Oo8uborojf3EzX+l+RZhEuZAI5u88PY+pHk19+yjy2rMmoCyRlc8/G7jZoe9CY24gw0MYRKuiMjCeEiA+VAhJNmd5AVXGs5MVLn0lPiiejlvRP0vbSuK90LucpYFh9YAtb2n00XrbflmVOmBhvkX7zxXbcbEYS5zxP4QZWG2qhW3blvKavNCborqTf3E2XVuai0hbkU77213YiROueMOZilwWW/9L9zk6/6ZclIFJQE23sqAiLcNTKC9XmSODoWHTRILDvA/1oHO1JFB7G+1nGSLkeUE8tj5xDeEtgJgaV1rfLZQvm99MfOfHDkscLSLc+jnqy3iMQFnZfmn9uN9ZOpiyMjdOe+ZUH1ltSbZLmfWFCXSxddkksj2o14uzljDiZYmn8AHuX0nm/r5qoEqhgNFn9UyBv+vqVUvQId2Y6A8MXwYDvps2Cr8vgT0NnLUjnsmgdZ7QICW9C19Xa5dDX/Lle3pe7sFI4or/uWUorIlHNuvak89++ZOtpK8jVzSXXDFvlpyhi2NmpQRLjJHz1InYeWwvdPYYsFOqtTOSw5mmTl5lQslHNf9+Qliw3g1qxr662hrAu8NP8uUKmFb2kKR0Lb4b6lsC5yxJ1Lb0HX7t9zlLOxtGvkkurEBMBzjSsnfaObLj9j6QxTkixp4gMdu/IIsCSo+lPQ8qIPJNIQsf52Nqf7YGaPpU02aNpRXqWurbd2ZV1gqPl3gUotfEtTOOK+pbASJoibXW/uJyZoZZtZ1sil54JyzK7ArYgfzNSQIgy6n2qw6iVFrZCdFih8zzKLcloth3PrPrCm2cb6ejj1PlOsUdfW2364aP7tR9dT79QcmYrcsvnm1pv7iWX1e87S18SlwLuTVosdGTWAKcF8+M1SKWZu7IyAETACRsAIGAEjYASMgBEwAsURkN3B7mhvdXwl22PySqOD5WdJLdkt6akK8TK0BBR7jYARMAJGwAgYASNgBIyAESiDgGwNPmdgMoXNASYbNNSmdaamirixmiiE7yZOKgR5dkbACBgBI2AEjIARMAJGwAgYgYhAWB32WrbGyd9odxo1FBamg65U0LtYuM9GwAgYASNgBIyAETACRsAIGIFTEAizNF/Lzsj+o822cv8PY4mfAbzaDhoAAAAASUVORK5CYII=\n", "text/latex": [ "$\\displaystyle \\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": "iVBORw0KGgoAAAANSUhEUgAAA5gAAAAsCAYAAAAAexdDAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAXUklEQVR4Ae2dj9XUNtbGFw4FvLAVBDpgkwo26YAkFSR0wB4qyCEdkFSwSzpItgICHZCt4EveDvien1/JyB6PLXnssTzz6ByP9ffq6rlXsq6lke99/Pjxb3ZGwAgYASNgBIyAETACRsAIGAEjYARyEbh3794r5b3R9VjXH7r+Jdvy9p4NTEFhZwSMgBEwAkbACBgBI2AEjIARMAJZCATj8rVsSQzLvyn8RrfHCv/DBmYWhM5kBIyAETACRsAIGAEjYASMgBEwAiAgg/KDbl8lBuZThd/penJfP3ZGwAgYASNgBIyAETACRsAIGAEjYARKEGBrbHS3wfP4QYzx3QgYASNgBIyAETACRsAIGAEjYASMwBQCWrl80svDCibud69g3gHhXyNgBIyAETACRsAIGAEjYASMgBGYhwAH/jSH/NjAnAfg2Utpn/P3XGev2BUaASNgBIyAETACRsAIGAEjcJEIyL640cUBPbOdymNc/qZVzR8h4i2ys6E8X0EJ7YVq+5ZTmc5Xq2syAkbACBgBI2AEjIARMAJG4JIRkH1xK1vjta53c2wNlWMB7EZln0ecfIpsRKLSu4T2pVjjrcJnKEClbJotI2AEjIARMAJGwAgYASNgBHaKgGwOFrS+kL3xdW4Tgp3CSbL/okwI/2kDMxfBDfJJSHy49H+6/inBvd+ABVdpBIyAETACRsAIGAEjYASMwBUgINuDz4zwbcufppqrvBzq87OuH5K8rGJ+bQMzQaQ2rwTX7IcueZNQWxvMjxEwAkbACBgBI2AEjIARMAL1IyDbg8+OYGRO7pxU3r+Uj8WwjpPdcs8GZgeSegLhrUDzsVIJ6o96ONs/J8KWPyLTIehEYNuceKW73YUhYFlfmEB31hzr384EtgG71pENQD+xSsvsRABdvINAjfoknk5e4LKB2RFzPQEJF+PyD69eLiuT0JFZ+m+M9tCJHivsA5SWhXpzapb15iK4agasf1ct/qzGW0eyYKoqk2VWlTh2z0yt+iS+2PqKHfJQ8+NZ57/4MyUVqqcEy8oawm3+MFshi3tm6VmPefaNPw2Y95Ic3DkClvXOBbhz9q1/OxfgGdi3jpwB5IWrsMwWBvTKyVWpTzIqOfeFa/bnEW1g1qnZGJbv4ypbnSzumisM+Ojim5k0Lqb5vn8EUrla1vuX595aYP3bm8TOz6915PyYn1qjZXYqgi6fIlCrPr0Wky9TRkv8/g5mCVrny/uNqkpPZDpfzRdek4z2J70mslKM+/3u5t9LQcCyvhRJ7rMd1r99yu2cXFtHzon2MnVZZsvgaCp3CFSuT7+Jyxvt8PtSfOIvcvfHcrMHV9evuv7ShSVrtzICwpk3GRxAUyzMUtZU1ytdL3R9H/zR2JokpfzvQjl45bs36Aofae0v9ze0FP9G16+ThM+fgQN/fMjPwrhXKm/LemE5l5CTTlz7eGP9G1EY60cDjnVkXEdqnEdYZiMyK0nyGFDXGCCjkrNK2Pn1VYkcY95RA1PE2aYJYYyI2caBlOaDrnQJONbv+yECXxIF9odJy8VIHvx5963q+VEX37phxZTBO1dO5OOlAy8fPuoOvVvR+kX3Y645lepY4rnjxTcPht/A4Nx1X0l91cjbst5W44T/VY831r9x/bt2/QAd68i4jiSpfq4kYFyK12NAtWMAi12NXVKqa5NbZCX0SHjWiprKs6LFKZ3+1EaedDjNdFWsJBP+tItMWmNQ/lvFE2bwzjlRFR7RCQzNt7p+GZOx0r5WnmpcwOBGfPFBWLuFEahJ3pb1wsItJHft48016J/aGJ/zxS/rrl0/6E7WkbxBxc+VPJz2lstjQNVjAPP7wZ2JU3o2aWCKACuYfC6DZdI5jvKrrsbNYaq0TOgA4DDL0C6oD4N+VQNT9DH2htqBIrFlFsNrSt5gsctTbtU+MH4i/hvjMoT/VHj3eqp22SUIWNYJGNt5r3a8sf5lKd3V6gfoWEeydKSqTJbZ4uLwGFDvnLSxR6TzxQuF9zPUhMn4kDGSUbTJcmr53HrWzocxco6DYFgRXNvARCZDLhqVx9KHyuwqTp2E/5myNfat/M+45MdQXhtzVWF3TgQs63OiPVrXsfHkoscb69+oTqSJV6kfAGAdSdVgH37LbBU5eQyod04a58bZZ7REDXkQPdzVcTBumGzz4P8/XRiWEOX/eQdO+Zmokw9Hmd9ZBVI8ysIbCehxNQfA6P5B6cVbaFRuUxfaQ9vihGhtfj6sVYHachNo/zlSBzKbdKIVv48DTU5nfZ2uAoa6yENaR/YBU1a3o87J2y7DfyF/8/9Q5UP/vtWFnhFPHYMvPJQXYxGaET86xnPlJy66/8oDv2wFbt0ZZdvWeWmeCuW9a1lP4Imuv5TeNtvZg+7/rLjPatLl0AZU/RrHm13rH0Jb2125fgCvdWRCySbGQeaa555HVCWzCXyqf054DKh+DIgG5qOJrnqYrMkIB8qQQEf9SxfLoDGOg31IZMtkExfviuPQhqdJmLxM/tt8CqPcivoUN+RXnrbOofRT4kS75XEOHZXHGKGtBxjMoTdWJtQF3t+P5TslTbQx6KjjVZ+O4jDmSHvRT+uHlQd5t5jIH+lynHHUn6YOpTW6FeO5y8U06HC1dcoP5g0fuqf00Ke/UjrRr3hedvyahCM/HZ2M6b6P98k5+CADyslZ3hNjXg6+CZ4Heq80+sybSCdgDvhpf6FPcBAXebm3/TWWW/uuOmM/9HizgE6sLa+59CVndLQdw3PpWD+WH4dzsT93vhN0xM+VkbFDuEZ8dvmc8BhQ/xggGTXz8dIx44EKsnIZV3T4XEO0VknCz0mytwSiU34mDRhu6X/WhrYZ8mYpzRNJtHfRonNQtnEKv5Dn77pYsfp7mtZkGPiZKPNY6Y9EZ3DVa4BcG6VyGFw/6/pO5TsYtJmW9ZS/IZhfPzI/5sB/1AmPdFUQg/IP4QXGTGSfyI+O8J9OHKvZrV4FXGMa+W5VPl3Zjjgg/1RuxB/wLXrN5EZpD3U1LvCDn8m13coIWN7LAtzDk77W9p9QE0Z8e0AV/URl6CtNPvkb45J+QH6FWbHnzXvOAV4UWdod9NukAo83CRhX6rV+XKngx5rdGwc9j+iB1cNn788JjwE9+VYWnHxO9/ltDExFMhlBuP/pZWAS0540mqQ1252k3Kx4UoY36akhELMeK9+kh0lPalyynbI1KpXOf+SgzcAy6JQ+WkZlf1Eevs/IoTT9SdogTSKVH0wwdm91/aywbrPcD/CQWTJ2sLHtZPGFAJPFEoeRjLE/Rjsadhj3cxz4fimsMBo5NCe2Gxm1cpYfWbCVmvaSt2OsKvxUF66/NZt8Q3rGSwA+N4KsGifakcZQ/pgNLIv/uNwWrsRDW4NsszkK2M/VoaF6qpd3bbKekFuK5zcCvO0LlAsC6Ot2Osbx0oWXPdFRnm/Xjuq70umTi+pFZGDgflXjzRT2A/isGoUezRg3eC7yXO879OaRaH7bT1AYvTz2DPfzKAHMOpKAcedNx8Eq5xFryiyjj6b47PU54TEgUfuN9Snh5MDLGF/mNPBTAEOx3V4Y4iDW2XJFfLyUxkOG1SHycJWWZ5LU2Tal8Add7RavwAe0j27tyilDeV3tdrLYhqm7ymBgwtNqW3hTHlQPmNDeZ2n80v5QxwEe1JtTv/KwGvKuz5fimHzAfytD+SPNAxnGtAE6TIz5z2arb/jl0NPONiyFI2adbcXk03VAI6Wp9I7+pWkl/sBDR/9LyufmPVaP4sG4xTyX3hr5Ai+DfTam9etV/OryVh2byFr1oofxZVWHB8VPyk15GGfBsx2D5IfmwVZxxbX05e+MpQrzIqfTN/tyWCsc6r3q8UYYtLI5BWfRYbwrGmuOlVH8pP7l8hpodcbmgrLo5VXrB1jJHeiI4orlnYt7mu9YPYqvRUfgA5CqmkeIn5NlJhonPSOC7uz6ORFk6zFgAX0K+jA4bgjnWf05yKf472YPZC1jfHGxupW65i2lmO28Jce6VibefLI61aQprukgPcv7c4il5alL4Vvi5V7qSt/Kw0Oz8kVi4sgPL3E1rE2CngKTZagT3nr8tXSOeVTuR5WhjRhU59haFrF5dIynheJpE9j1Xay3I/N+JoWRx9BqcKT5e1Lmufx8IzO2LUlqVi77ekf6gbwlB+Kg3+iBwlGXGj1TfJ/nY6udytpuGUxXVZv4kh/xQCemfTj0cBU3VY+wnbVKvwqzd3hUJW/hR/89u6xV7xK7K1j5eS8Zp/3tQLdVF/rXbgdXfg7WSh26ikv75l3M+r9XPd5sqH+j41NF48ZV6wfdr68jCo/KbqkuO1VPRTpS3TxiCZmJxhLPCNRh788JjwG9ecpU3xwaA6bKzOnPohnn9UNz+CE2OnEU5s1QZ9VM4XaVijRdzRt03Tnw5mBlSHHi/dOKk8K8tW9XueSHBoZpsyolf4cGabrgo31TT1453sR3VqcSGtllRAMjeJBOpHfsTjldB2+qjuWfG686oixmvQnOrTe0Z3AFRGmtzI7RU57BNxmKR1ZtefljexrdUriDocKsSPbjYpnOipzyscrV0NadyXSkiV6J1U+6h5+4mGcgDb3p1NvPUxIWLYzfjj6XlM/NO1aP0sDt4A1gLu0l8gUeWtwV7mCs8NnlrTo3k7Xq7qwigrEcP+1bePy6jspNaRiNnXSFodvHdnTMCGVG8yyhA0M0VDfj51WON2r3ZvoXZSEejo5PShvVv0hj6i46jMOz9EvlrlY/wFXuqI4o7ajspmRSkj5Wj9I21ZFQ/26eK2NY9mWivCc/I4IO7fo5IRw8BvSe6VFXSvQpp4zoFfVn5We+Tf8rtp/ui6FbFeTtQbsCE96qMOjFt+ZfKF/0s1+685855cd4668Q0IimTLCAMRyb1SqF27qUp++g33fQGnM5ZeClw/cYwTRNfP+kMA/QVV2QBXUU/5m2hLHQnj8lh7ZNQUZ08u8iLeJ0fdTFS4XUvVIcBl/rFKbsI13pf20+J4PqY4UNmTOYNi6EkWu76hKSeKBSBp1MHbRiHJ8eiSvaTVyg1+SXn5cjuJj/LvTpl9Xzf38K7t8nPOjHzSr9hq2pUd6byFo6iG6j8/2xCTk1Oo6cMuRG/2jHS9GlLP2M8blxisPfjK0hqnNTOi/7+I9yepBWJ8+aAdXL+Hmt480m+pcrzwz9yyU1O9+V6we4WUfGtecinysLPiNAb9fPCY8B5xsDZoz5zGVw/bnMXezI74OQhlHAQTYYijgm5hhjGBLEpUYAWxU4fEe35gAcKudP/HHCTzyOSQ2H62B4MJFKJzdMmKLBSnLqmDwxCYuO8JTLKQPNCNQUvaF0jCRW1Y4ZLUNl5sTBZzuhnEMgswxbfl+qTbGuLxT+p9rXTlTlZ2sxcurISvGcGBuNzMgvytf5Bh9YKR+4RR1gohldowPkiRHhTnyaLyb/IM9XgdbrGBl4jLr6VvEY5w3PpMV8vXvx4Ra98rUGMZoxQIbwW53nSuW9laxjv+rrIP2kP6YdlZswZZs+pzKj8/EFzWfyx/EafWf1fFDmob+wYhq3civ7Ju5ax5ut9K9EyEf1r4TIiXmvVT+AzToyojwX/FxZ5BkBdBfynPAYMNIPFk4qGfPjfKU/l5lmSYqJcp71EldMgvvbvm4UxzIsg23LT4jrbJeM6UrLLgNdXe32zUijtrt4xJivns/acEv5CRge28bLoL7odlbRQ58XpZm2J/qn6lH6M12dfhXLXvJdbabPHMhbcZvJWnUz3jCetdthkYEcRmJnK6HCq8hNdNHLdittCHfG10vWi3O1TbhWp39p24Pcj45PSj9Z/wKNjl6nPFyz/5h+gInc6Bil9It5tqgt6NludES8zurXuTJTvs2fEUEH/ZxIbI41xqpjurTWGDClg0rPHvOVlwWig7lMDk73VXALx1t8jMPWiVmsY1aeorXcpskzeDBFYRnodlbi0goq8rOCyEPHLgMBrdAwSLdOYfSKATNuk23Tggds96AHfb5zwvShTr/KKbSnPIXyrkHW/fGsHwb+xeUWcHol2m/lZ8cJDxT+xnCpuq+mre92qH85oCyhf+xEGVxFz2HgUvIU6gfNrmGMyoH/onWkUG5Ly6z/TOiHkc8S+B/IObTbz4kDZOZHFOoSFS2tTznMl+gTBwbeBnsrh3ab50HrO6+HSQ6g9h3bwNiu22yblKCYFLUnUipMmSZP0tjRMsof3RZCjHWX3Nnm+UJtZeUDJbA7goAwYmAEq4cJVhiWPync33qbUjnAVTQwSvv/I07LRP/XSV0xLvu+cj1D2y+zeas940x5byXraMjd9HAlzEuk1K0ht/+qAurqvGg5RXdThq/RvzP9KxHRyfpnvfr07WzpycnPoxLhxby1P1tq1ZGl+nWUQ8F962cErPo5USCwqawzdQmyB/OUqbpOTC8Z85kbj82nj7KyiYHJQCNBcDR3x4hSPP83wmDAcMCRnh4aw4Tpc100uPnPZ0YZZW0c+7s7k60QX9stCvIbMXb1b4QnhMMWFl4cfIMu6c6bFrZK9v8PrOjW0bHI23EqA+4R+07akoGV6+HNZ3xoLcl2LbRK5b2ZrCXn+P/lobfR/R0Zi8tN9T+sRWgXxMdu9K8Q88X1r7D+S8leqh+0e3CMmgOIny1zUGvKlMptEZlt/Yyg5X5OzNaZYwVLdQk6i+jTMYaOxGeN+WFe/VQ0WMgrdpsYmIFLGG4Nxci5FD49DChGN3elvZfnoRrNymbrxsq0mWSYKt/WB10k7Ax7xSMTUwwdDq6xgTkMUxMrrOYYhcdWz0dqmkw6MFgpITkOrbhPEhvJMFhPkn8vq/QJy/neGfLeWtaXtrsiX1gXmHOH+helcNXjRgRh7fsM/YClqTFqUHZ+tiwnzRlyW1JmfkYsJ8rNKc3QJXheUp8iBoPjRkzUPXeuiI2Gm7X4cv+u7Pl/JQgMyWKDLxiXGJrZLpSZZYFnV7JsRlZao2CXpXzl1KR3t0AgnZjqgJNI8ZDXhV6x4o7/ja54EjPlqSOuuBOe5TLqiXRZpecNmp0Q2FrWYYz7IPlx4jI68pXi0h0ZUU6WW0Tigu5b65/HjfqV6ZiOZMjOz5aNxLukzESLebCfERvJsoZql9SnjHEjNjl3zsFC13vxiBFc7O6pYHGhpQoIDIwovo+ZvVKnMnzssyQ/AzEnKRYbs0u1s5SO2gjPfJSe//uNbfcsJe38QkD4YgQOfVpnFXxU37NzyFH1cPowA4ddQGAPsrbcLlddrX+XK9ulWnaKjqisny1LCaKAzrllpvr8bC+Qz96y1qpP4gs75Lu589dNDUyUQA1gu+tsCxkaYy4IjkNfbsfy1ZYW+P7WBsM6khG+vwpb3s6s6tbW78h8qOeR2pT98iWWvfR7zbK23C5d+5pnXLVjjfWvDv2bM0YF2a02d4rIWEciEt37uWRm/Lu4X2qoNn0SPyzEvNSccvZ5DpttkY1KIuY5JXbW8mukMXYX7R917cq4pD3wrRvbLvmDrd3yCLBtke/7rO0w+lbTb5hXO1jxZvuljcthaVYpa8ttWFgXGGv9u0ChLtykOTriZ8vCQigkt7rM/IwolMi+s9emTy8FZ86XFY6ivvkK5lHOnIDhgAH0XIaDtz2uoA/Cd9XV8xVYHiSpdvCmaXer9IONWSmyRllbbisJu0Ky1r8KhVIZS9aRygSSwc7aMvMzIkMIF5SlFn0KenfyDkobmJUrpwTNoS1ssTp6um7lTTB7RsAIGAEjYASMgBEwAkbACFSMgGwOTph9p+sfp+6+23yLbMU418Iap04+l9C9VbYWiZgPI2AEjIARMAJGwAgYASNwIQjIzuDvVixqcbDPyX/t8grmDhQjvFFA6PzP7mSh76DJZtEIGAEjYASMgBEwAkbACBiBMyAgWwM7443sjEXO87CBeQahLVFFMDJvJPiib4AuUbdpGAEjYASMgBEwAkbACBgBI3B5CITVy89lY/y2VOv+H84LF1TGgi7oAAAAAElFTkSuQmCC\n", "text/latex": [ "$\\displaystyle {{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+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XuQXOd53/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": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "try:\n", " 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);\n", "except ImportError:\n", " print(\"No requests installed\")\n", " img = np.random.random((82, 290, 4))" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6IAAAEfCAYAAABbM3sFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmwXNV99vtnGSSEBjRLgASSAFlICCSMwAwGZGMIIkwpm1vYlF8FE/PGBYkdXLn2+1Y5JpWkyqlKXl/fVEKCY25IxY5tYjPYxmYQyAwGrNlCsxiE5gnNQgKRdf/QYenpdu9zdp/u3qe7z/dT5eJ3Wt29V+9ee3dv97N/O8QYBQAAAABAUT7U0wMAAAAAAPQuHIgCAAAAAArFgSgAAAAAoFAciAIAAAAACsWBKAAAAACgUByIAgAAAAAKxYEoAAAAAKBQNR2IhhCuDSGsCiGsDSF8rV6DAgAAAAC0rxBj7N4DQzhO0mpJV0vaIGmepM/EGJfXb3gAAAAAgHZzfA2PvUjS2hjj65IUQviBpJskZR6IhhC6d9QLAAAAAGgFO2KMI7u6Uy3R3DGS1tvfGzpuAwAAAAD0Tuvy3KmWX0RDhdt+5xfPEMKdku6sYTkAAAAAgDZSy4HoBkmn2d9jJW0qv1OM8X5J90ul0dyZM2fWsGgAAAAAQDOYO3du1Y+pJZo7T9LEEMKEEEJfSbdKeqyG5wMAAAAA9ALd/kU0xngkhHC3pCckHSfpgRjjsrqNDAAAAADQlmqJ5irG+Likx+s0FgAAAABAL1BLNBcAAAAAgKpxIAoAAAAAKBQHogAAAACAQnEgCgAAAAAoFAeiAAAAAIBCcSAKAAAAACgUB6IAAAAAgEJxIAoAAAAAKBQHogAAAACAQnEgCgAAAAAoFAeiAAAAAIBCcSAKAAAAACgUB6IAAAAAgEJxIAoAAAAAKBQHogAAAACAQh3f0wMAAADF2bt3b8nfp512WsX77dmzp4jhAAB6KX4RBQAAAAAUigNRAAAAAEChiOYCANBiFi9enOqrr7461bfffnuqFyxYkOo1a9akes6cOSXPdf3111dcxnPPPVfzOAEAyMIvogAAAACAQnEgCgAAAAAoFNHcFrJs2bJUn3POOam+7bbbUv3222+neu3atak+fPhwyXN96EPH/j+IEELF5cUYu7xP1nNu3bo11YMHD071Pffck+qf/vSnqZ47d26Xzw8AOGr8+PGpPnLkSKq/853vpHr79u2p7tOnT6oPHjxY8lzbtm1Lte/3AQBopC5/EQ0hPBBC2BZCeNVuGxZCeCqEsKbjv0MbO0wAAAAAQLvIE839N0nXlt32NUlzYowTJc3p+BsAAAAAgC6FPDGcEMJ4ST+LMU7t+HuVpJkxxs0hhFMkzY0xTsrxPGlhM2fO7OaQ29PSpUtTfc0116T6Yx/7WKo9art58+ZUe8wqT4Q27/3yPldXPLLrEbJ333031QcOHCh5jEd4lyxZkupnn322LmMCgFb20ksvpfq///u/U/3ee++leuDAgakeN25cqpcvX17yXJdffnnFZfi+GwCAzpSdZrcgxjijq8d091NmdIxxsyR1/HdUN58HAAAAANDLNLxZUQjhTkl3Nno5AAAAAIDW0N0D0a0hhFMsmrst644xxvsl3S+VRnMhvfXWW6n2qPKQIUNS/eSTT6baL07uXWk9inXcccelul+/fqkuj1jV0hkx67Ee5fUuvR7HPeGEE1I9duzYVE+bNq3kub797W+netasWan2bsHeRRgAepNTTjkl1RMmTEi1x3Tff//9VHvX3JUrV5Y8FxFcAEBP6O6nz2OSZnfUsyU9Wp/hAAAAAADaXZ7Lt/ynpJckTQohbAgh3CHpm5KuDiGskXR1x98AAAAAAHSpy2hujPEzGf90VZ3H0it4HPeiiy5K9aJFi1K9bt26VHvsdujQY5drHTNmTKo9gus8Qlsep/X4VrWyorl54r4eFdu3b1+qf/azn5Xcr3///qlesWJFqr/4xS+m2iO/q1at6nLZANAufH/r+1XvRu78tA0AAJoBJ4YAAAAAAArFgSgAAAAAoFANv3wLpMWLF6f6j//4j1PtnW+9y6x3zR02bFiqvSvt/v37U7179+5UH3/88RXrvPLEa2u5j9/uHXRHjx5dcr8DBw6k2mO6//zP/5zqv/zLv0z1qFHHLmX7/PPPdzk+AGgXeU61qKVTOoDm5VcQmD59eqpnzJiRao/m+3evvXv3pnrPnj2p9u9XP/zhD1P9yCOPlCx77dq1qb7lllsqjs9PPQPK8YsoAAAAAKBQHIgCAAAAAApFNLcAffv2TfWCBQtS7XGqQ4cOpXrEiBEVn8cjFH5/7x578ODBVHukNauzrtT4OG4W7/RYbsCAAan2mK7HjZ966qlUe0dhAKgnP71i1qxZqZ44cWKqL7/88lT76QHPPfdcg0eXTy2d0gE0L//+c9JJJ6Xar0ywdOnSVPt3r9NOOy3V/r3y3HPPTbXHdC+88MKSZQ8fPjzVflWInTt35n8B6NX4RRQAAAAAUCgORAEAAAAAhSKaW4Bt27aleuXKlan2eIRHZz/0oWP//4BfnNxjqR6/8Giudz3zx/bp06dkTHkubt5TkV2pdN14tHnXrl2pfvTRR1P9hS98IdU333xzqss7vAGA89jttddem+rPfe5zqZ43b16qV61aleoHHngg1e+9916jhggAmfx70cMPP1zxPh6hffvtt1P9wgsvdPn89913X6r9FARJmj9/fqr9e6V/Lz3jjDO6XAZ6L34RBQAAAAAUigNRAAAAAEChiOYWIISQau8G69EF72joEdwTTzyx4nN6lNcjYR6DPXz4cNVjzYrRVhuvrfXi6f54XzfeCdhjIB7f9e7CANCZcePGpdr3md/5zndS/eyzz6Z6//79qfZuld/73vdSPX78+FT7qRZFqHXfC6C1+H4o6zujd7H91Kc+lepNmzal2r9r+SkLfjqCf9eSpEsvvbTi/XxMnV0hAeAXUQAAAABAoTgQBQAAAAAUimhuwTz6kBWh8ijvkCFDKt7fO+L6/U844YRUe6dcjwF3tuxa1NJBt7PH+r95bNm7EfsF5P0+HmEGgHLeyXzFihWp9lMevEu5X+jdLxI/bdq0VOfZz9cqz+cHgPZ36qmndnmfDRs2pHrLli2pXrduXaonTJiQ6kGDBqXa47hbt24ted6RI0emesSIEammizjy4hdRAAAAAEChOBAFAAAAABSKA1EAAAAAQKE4R7QHZZ3j45cQ8HNBPX/vlxDw85GGDRuW6kOHDqW6/JIm5eeMVqNel3LJe16o80u2+LlQfj5D3mUAwMknn5zqM844I9V+yYGsy2u5ZrlEAfs8AOX8+5LX/p3Kvyd6f5J33nkn1du3by95Xj9n9MILL0y1n2+ftc8EJH4RBQAAAAAUjANRAAAAAECh+L28YB7xyrrd47j79u1LtbfQ9iiqX7Jlx44dqfb22XljY92J0dbj/rU+L3E0AN3h+w7fZ7bS5Qf8NXjNpVyA3mv37t2p9giuX3Jl9erVqfbvnqNHj65Ye+RWkq644opU+2Wu/HIxHu0FynX5i2gI4bQQwrMhhBUhhGUhhC913D4shPBUCGFNx3+HNn64AAAAAIBWlyeae0TSV2KMkyVdLOmuEMIUSV+TNCfGOFHSnI6/AQAAAADoVJfR3BjjZkmbO+p9IYQVksZIuknSzI67PShprqSvNmSUbcrjuB/60LH/T8DjER6t8Nt37tyZau+A6xEI775bRHQ1zzKI1gJoJlmxViKuAFqZd771Kyq88MILqfbvjH7K12uvvZbqLVu2pHrPnj0ly3jssccqLvv0009PtX+/BcpVdY5oCGG8pPMlvSJpdMdBqmKMm0MIozIec6ekO2sbJgAAAACgXeQ+EA0hDJT0Y0lfjjHuzfv/EMcY75d0f8dz8LMXAAAAAPRyuQ5EQwh9dPQg9Hsxxp903Lw1hHBKx6+hp0ja1qhBtpOsiEJWV1vv3Lhr166K9+lOhKzaiGzW89YSx80rKy7nfHzEfAHUU559UDNq1XEDqJ2f2nXgwIGK9/HvmNu3b0+1X7HBv1/5qWBSaTfeAQMGpNqjwN5NFyiXp2tukPRdSStijP/H/ukxSbM76tmSHq3/8AAAAAAA7SbPL6KXSfqcpKUhhMUdt/1vSd+U9KMQwh2S3pJ0S2OGCAAAAABoJ3m65r4gKSvveVV9h9Oesi4efMIJJ6T68ssvT7XHdPN0dPTuu84jFB6fkLKjtmPGjEn1jh07Uu0d1Lx7b6Miu0VEftGali9fnuopU6ak2i+s7fPyhhtuSPXDDz9c8lwvvvhiI4ZYqFdffTXVf/qnf1rxPs8880xRw2lpWfuzrH1sntMiGrWfYv/XHpYtW5bqqVOnptq/E/gpPddff32qf/zjH6f6pZdeatQQW16edSyVbsM33nhjqn/yk5+kupU+M7xr7gUXXNCDIwGy0VMZAAAAAFAoDkQBAAAAAIWq6jqi7S4rvvHxj3881d5hzGMceS/YO3bs2FR7F7M+ffpUfN483W2zYmP+nOXPk/UYfx3VXsS96A66Pr5Fixal+vzzz0/1rFmzUp313rnOusP90R/9UaoffPDBVP/yl7/MNXYcs2LFilR/4hOfSPW0adNSvWrVqlR7tHzw4MGp3rt3b6qzIlNz5sxJtcfKy/++7bbbUr1x48ZUL168WD3F15Pvhy677LJUjx49OtWvv/56qssvPN5TfP194QtfSPUpp5ySan/vjj/+2MeSd3288MILU7158+ZU+2kD5davX59q36d7R8cTTzwx1YcOHUq1z63+/fun2k+d8G6Q/vnhF4Yvj8T5Pubtt99O9SWXXJJqf31r1qxRs1m5cmWqzzvvvFTfeuutqV6yZEmqfXvK+tz0+/zBH/xByfI2bNiQat/XN7u1a9em2mOgkyZNSrVvBzt37kz1888/n2pfT08//XSqff/l81gqnct//ud/nmo/taFVI/u+zWftC/37lX/v2rbt2AUefB2Xe+GFF1K9f//+VPv27/uqu+66K9W+jolMA53jF1EAAAAAQKE4EAUAAAAAFCoU2XkvhJAWNnPmzMKWW+6tt95K9eTJk1N99tlnp9pjQh7x8FhR3uiq32/hwoWp3rRpU6o97tEIecfqkV2P9nqUzW/Pivi6enbQ9YiSx9fGjx+fao/CnX766an2CzTnjVI7jzr5c91zzz2pfu6551L91FNPVb2MduKxNEm66qpjTbbPPPPMVM+fPz/Vvl6zeFzVL8CdNRc91nvyySeX/Jtv/x7F9I7RHrny5dUrcrV06dJU/+Ef/mHJv02cODHVjz/+eKoPHjyY6pdffjnVo0aNqrgM3w6K5tHcb3zjG6n26PvPf/7zVPv+xdf3WWedlWrfJ3icsZxHX31d+v7WO5n7fsHH590nnUehs/Yp55xzTsnffr+s1+d8LjqPJ5522mmpLo+ff8C7tJfP3SuvvDLV/n551PbSSy9N9bx581L95ptvprrWfewHPAopla6Du+++O9Ues+zJCKS/F7//+7+fao9o+3r100R8f+avM8/+zGO9/v1FKp0Hu3btSvXFF1+cav9M8+9FPSkrzuz7al9nvv37enUelc2zjqXszw3v1O77oa1bt6bat3m/z5YtWzKXB7SDuXPn+p8LYowzunoMv4gCAAAAAArFgSgAAAAAoFBtHc31C73ffPPNqfbuiR7n8c6NHqEo755ai2o70RbN54NHuTya5p3i/P7VzqW8988TzfUYjatn5Nm7Eo4bNy7VHvnzmK7PrSeffLJu42g23tXzk5/8ZKrLo4Ye53vllVdS7TE8j5N5d1yvBw4cWLHO6qh8+PDhisuSSuNUHgn199qf9/bbb6/4vB6PrZbHzD7/+c+X/NuvfvWrVHu0zyNlffv2TfW5556b6n79+qW60dH/znhUxyNrHoV+5513Uu3bsnf+9fh91hwo55E8j0n68/ry/Ll8zvk+z2V12fXTOcrjf74P83F4vNaX7Z15XbXRXJ8nHoeXSiOe3rnau5C+8cYbqfbPVo8/+2elzz+vna8Lj5sPGDCg5H7etdi3Td/f+PZYRKdhjwX7dufbmu8XfK74fcaMGZNqf91Zc87jpz5/yrtkewTVP8t9GT73/+zP/qziWBvVWdf3Z7Nnz061n8LxH//xH6n2fbXHwX3e+HY3dOjQVPt27d9lfJ8vlW6rXvv68Jizf0/0deyfP9ddd12qiemi3RHNBQAAAAA0PQ5EAQAAAACFartorsc9vvjFL6bao0f3339/qj1m4dEej3J4rChrfXW2HvN0ls1S7fJqWVb5473OE8GttTtunvtlRXM9dpYVq6xVVlTHI39eP/jgg6n2i2PPmTOnbmPqKVmxNF/fHnWXSiNlHl3yi417ZC0rEj9hwoRU+3btc8C7xK5fvz7VS5YsKXku3+Z9fOvWrau4bI9f3nvvvan2+dDZRdIr8fitd4aUSuec74c8duaRSV83Re7bO+NRnY997GOp9tfm27XH5fx98E7LPjc6e50ee/bHe9zVY7Ae3fZunFlz0d93f3+GDx+e6vL4n//t792MGccSTL59+ZxzHs31eGfWWD2G6cst/7es/apvm77N+j7P30ePQHp02KOU/jqdr1epNIrp77e/X1//+tcrLqNe3XTLx+r7PY8q+3cKXx8eLff31CO4vj/yeendY33f6d39N2zYUDI+j6l6fNqX5+/77t27U/3lL3851R6bryWm6/NKkqZPn55qj17752NWxH3kyJGp9jnqc9873fqpNP79wLcbqbTrs39ueO1z08ft69hjuh6Rztr/la8boFURzQUAAAAAND0ORAEAAAAAhWq7aK5fIDwrduuxHY/2eWQlK6KapTuR2HpFWWvpXFvrOBrxGjqTFc3N0/E0z/K609U4K6oza9asVHsEaMGCBVUvoxls2rQp1R/96EdTvWjRolT7e1J+QXtfBx738jjVXXfdlWqPHnqX3fLnrcS3R4/1le93fK6sXr264nP5hd69I6nHLz1SW969sivPPfdcqi+99NKSf/N1kxVl9Xiex0+bhUd1Lr744lRnRUh9P1zra/O4nPOYqX82eITP17e/p76v8e6s3vnT9y+//vWvS5bt751HAbM6fw8bNqzi7dVGc7PmTPk4spbtp6v4duAd6b0TqHcs9vpTn/pUqn0+PPTQQ5nj8XXun/G+7/W4pr+P1W6PLusUBKl0Hfp3Df8M8XXmz/WRj3wk1X5KgD9n1n4ua9925ZVXltzP921ZXYR9DjmPlnqn66yOzFn88+Ciiy4q+TffL3hHXH/dZ5xxRqr9dZ933nmp9jngn1F5PifK+TJ8TJMmTUq1rxu/T1a3b+9S7O/DV77ylYrL9VM1gFZDNBcAAAAA0PQ4EAUAAAAAFKryFZNbzMqVK1Pt3RC9i6FHZ7yDn0ciPK6VJyraLHHcei23Ho+p5rFFd/WsNo6bd3we1fEuf48++miqv/SlL6Xao2yPPPJIrmU0A+9O6BEjj595PG7EiBElj/fHXH/99an29f/AAw+k2mOVHmvz6KE/1t8vv927dz755JMlY7rppptS7V0tfay+v/AonHce9deaFVvOozxWmSdmiWx+GobHPX07dR6j826wWTxW6Z9DPnfLY5E+Nz3q6LF+X3ZWNLda/nlV3i3Ut23f1nw93XHHHan2uPrTTz+dau8G6x1t/XPZu9h6jMvjj37/8uf1SLx3VV62bFmqfZ904YUXpjrP9rh27dpUe4zYt30pO/bty/bnuu2221Lt63jevHkVb/fTF5xH10eNGpXq8n3bjTfemGo/dcA7/Po+z1+fd9D12y+44IJU5+n06p+N3sFaKp3XHr321+TL9s8M324WLlyYao/K+vbu25zX5d/hstbt0qVLU+3drf01+Hvnnzk+Dt+uX3755YrLAnobvs0AAAAAAArFgSgAAAAAoFBtEc31KIxH5LJit9V2xHXdiePm0Upx3EZ30O3sMY2I1+YdU577ecTIYzjemdPjQ63EL8ztsXePj3qkzmOLUmmcyp/LY08e9/Ln8guBV8vjceVjevjhh1M9derUVHvHSX/vTj311FR7DNkjnR7ZRc8aO3Zsl/fx+eHvu3fj9Ci072s8uuq17wemTZtWxYjzq+XUhvLI9/r161PtEUiP1/7jP/5jqn0f5nFG3w7OPvvsVC9evDjVHmn11+DruDzWn8W73nsc0ve9GzZsyPVclZb94osvZt7Pv1P4Y7zz7SWXXFLxuXz/5/FVj/Xn4V1vfd1L0g9+8INU33LLLak+66yzUu2dwj2m69+pPFru3Yh92/J4sfO4avm26NuR70ud73t9nfljvVutR4dr5Z9F/rnhHdx9Xfr691i1z0uPtPt69ffBTwUBegN+EQUAAAAAFKrLA9EQQr8Qwm9CCEtCCMtCCH/ZcfuEEMIrIYQ1IYQfhhD6dvVcAAAAAADk+UX0sKRPxBinSZou6doQwsWS/lbSt2KMEyXtknRHJ88BAAAAAICkHOeIxqMncXxwokifjv9FSZ+Q9NmO2x+UdK+k++o/xN/l55tI0pQpU1Lt5zbMnz8/1d5mu9pza2o9L7TR51TWehmURi+jnueqVnsuaNblPWoZQ2f38+Vt37491c8880zFxw4dOrSqMfUkP+/NL/fgfFvxc8ek0jb2fg6Nn2dz+umn1zzOcv785ZfS8HOQ/DygCRMmpNovO+Dnc/k5S/6811xzTaqvuOKKVD/33HNVj71atcz3omVtN1ma8bX5mPwcx2a8vI6fS+fnqkmllzHasmVLqv3SLM7Ph/Nzq7NMnz694u1+DqHXJ598csn9/PPbz93zc9X9Uh++r/K5knUuop8n7OdK+v6s/NxvP1/S92d+LqSfT5i1b/TlVevMM89M9apVq0r+zdfZT3/601Sfe+65qc46F97fC1+X/nqy1qXr7DPDtxE/39bXTdZnhr+nfp5wPfnz+pzz99HPtfbPcj9f3F+bbzd5z4MG2l2uT8sQwnEhhMWStkl6StJrknbHGD/YCjdIGpPx2DtDCPNDCPMr/TsAAAAAoHfJdSAaY3w/xjhd0lhJF0maXOluGY+9P8Y4I8Y4o9K/AwAAAAB6l6ou3xJj3B1CmCvpYklDQgjHd/wqOlbSpgaMr6LBgweX/O3Rh71796baYyfValQct6fu353HN/pSLnnv7/E3f689CrN79+5U+xzwKEye+F9n46j2NXksySPj/hpaib9+3z7yrhd/jF8awzU60jhx4sSSvz1e5jGwRYsWpdovI+Pv3fDhw1PtsUd/39etW1fV+MrXZdH7hZ5SxCkI1WrndVl+GSOPlvt+0uODHgOtl1GjRqXaI8Hnn39+yf3yXP7KL9mSFRvNWjd+iQ3ffn1796hruYMHD6ba92Eeu/XPKL8MVL1MmjSp5G9/T33Z/t77aQuvvvpqqn2f55fb8dMLPvvZz6b65ptvTvUjjzyS6qzPjPK/s94Xv49Hczs7BaQRsmK6kycf+03GL0/m0Wbncyvr+2mr7neA7srTNXdkCGFIR32ipE9KWiHpWUmf7rjbbEmPNmqQAAAAAID2kef/SjpF0oMhhON09MD1RzHGn4UQlkv6QQjhryUtkvTdBo4TAAAAANAm8nTN/a2k8yvc/rqOni9aOI/wSNlRHe+853GHPDGQanWnw2ojYnfdicr2VPwv7/N4VMzjL7t27Uq1R4n8Plnd+OoZBczTsddjT/fcc0+qPZ7z+OOP51peM+jOHMjqMNqTvPPtnj17Uu3RPo93+5zz9+61115L9dKlSysuy7tp1lMrRbmaPWrcSuvS5Rm3b38eQ5eyu2A3Io6bxeOtHn+USvfpeTqnV/s+eodzj7F21p3Vvy/s3Lkz1ZdffnmqfT9X3tW20fw0BN+3eezZT5Hwsfpr83i2x2N9/+evP0t3Tjtoxq7UPhd9vZbH3SupZY4C7ao5vg0CAAAAAHoNDkQBAAAAAIVqfLuxBvAOd1J2bMc7WTZLHLfI+9cz+tHoLrvl/H3092j//v2pnjHj2BWBvIOpxyQ3btxYtzFV+/qy5qJfhB09y7tGbtp0rPH3mDHHLoucFbnyeendPzds2JDqu+++O9UezfXOkrXyeeY10S9kKZ8b3pXW/63ICKR/jpeffuPbYJ7u5y5rn7x8+fJU/97v/V6qjxw5kuqFCxemunxd+N8eK/bHN8spCH6lgYceeijV3mnYO/z66/HPqwMHDqT6iSeeSLWvV3+eduTz1E8L8/c9D/bVwFHNsZcEAAAAAPQaHIgCAAAAAArVktFcj5NIpbGiRscdmj2OW09FXmS+PG7lf3scyKNO69atS/XKlStTndXZr9qOdfV8r73b6sMPP1zxsVkXuG52RUTDi+Axqw9/+MOp9ihW1kXV+/fvn2qPsvn77hczzyvPum3VToxZkclmeQ2tOq+7M+6seHdPKR+D7/frNd99GQMHDky1x4B9Wy4/jcc7svupIc8880zF+/hpAEXz70he+2eO75/8tXpHZV8f3dmfuVbdvrJU+3pa9XUC9cYvogAAAACAQnEgCgAAAAAoVEtGc48/vnTYHiPJE/GqtlNu0THOPLK6ynYnBtIMMeLy98TjQ+ecc06qf/Ob36TaL0LuPEqUJ8ZVRNdc56+1lq7NRct6/XljlXnWTU/GlTxGN2DAgFR7BG3Hjh0VH+uRtd27d6d6/vz5qfb9VvnpBR9o97hWM8bUap3XPaWWMZU/thm2x8721bXsJ7Nem39+zJs3L9Xe/dS7nZfbt29fqj067LVH+XvS1KlTU/3GG2+k2jt5Z50C4++F7xdHjBjR5XJrnT/N/pnRjPsFoNXwiygAAAAAoFAciAIAAAAACtWS0dxyHtvJuoB0M8Zxq41xemTIY4QnnHBCqj2Wmle1r7Ve0Vx/Dd6lVJJ+/etfp9ovPO7dDffv359qX38egawlCljPyK7Py6xOvuhZhw4dSvXmzZtTnXWhcn/vfM55nG/YsGFdPhZAqaK3D99m835X8M8vP33E9/Ue320Wvj/zDsHeHTzre5R/v8jq8s6+DUA1+EUUAAAAAFAoDkQBAAAAAIVqi2iua3QX0nrGcT3+4nGZgwfm+EEgAAAgAElEQVQPptojpzfffHOqV6xYUbHuTgfTau9Tr+iNR2u9S64knXHGGV0+xmOPWR106xnBrdfzZsWeWinS5GP1iFrWa+vs8c3C9x2+Dfq2mff1fSAr0p4Vxe9sTO0mzykVRat1XvcUH3d31mszbo8uz3bQ6NdQ3kHXo/zecbZV5ZkreTr0Z+msE3KeZbfSHG2Hz3igKM396QoAAAAAaDsciAIAAAAACtV20dxaNKrDqv+bx3s8ZjphwoRU7927N9Vr1qxJ9SOPPJLqkSNHprp///6pLo+4dme83b1/ngvD++v32rvhStKWLVtS7TEXXzfeLbheY631MdXGL4nqtIc8c79WeeJreaP5raidX1s9NXuMuFZZr6+WbS1PtNmj2pI0YMCAVHuU3x9f/rnW6uq5P2sH1Z464euP/RlwVHt/YgEAAAAAmg4HogAAAACAQrVkNLeeMYZanivvYz2C4XHciRMnpvq3v/1tqjuLrH7g+OOPvXU33HBDrnHUS574aZ4IzwsvvJDqjRs3lvybrye/iHaj47i1dkXuznO1oqxIUi1dFZtJVvyv0V1sO3v+rHWbVTejaudN0a+n3eZ1O4y7M/WaK93Zrv0zeOnSpan+kz/5k1SfeuqpqfZTa3pSns/vPPOm2tMRurNva1XtsK8GipL7F9EQwnEhhEUhhJ91/D0hhPBKCGFNCOGHIYS+XT0HAAAAAADVRHO/JGmF/f23kr4VY5woaZekO+o5MAAAAABAe8oVzQ0hjJX0+5L+RtI94Wj24hOSPttxlwcl3SvpvgaMse5qiVXmvd3jKe+++26qx48fn+oFCxak+uSTT0711KlTU7158+ZUr1+/PtWrVq1K9erVqyuOqRnl6aArSf369Ut1tXGWRnQE7s6ye3uHwUbHWOvNI+A7d+5MtXfLHDJkSJfPkyeyVmssq9XWbTWyXluzRNmaZd1ndXrNEy9uZXleax6+/nzbd0eOHMlcln9G+Sk0vp37fqRZ+Getd9nP6hbst+/bty/V/p3Fbd26tabxtcs8/UC7vR6g3vL+Ivr/SPq/JX2wRQ2XtDvG+MFeeoOkMXUeGwAAAACgDXV5IBpCuF7SthjjAr+5wl0r/t/VIYQ7QwjzQwjzuzlGAAAAAEAbyRPNvUzSjSGE6yT1k3SSjv5COiSEcHzHr6JjJW2q9OAY4/2S7pekEEJdslXlkUePPtRyweBa47hZY/QOsK+//nqqb7755lQvWrQo1QsXLkz1nj17Uu3xIY8CNUtkrdp1nDdyW2RH3O6syzzzw2uPRvnt5RdMbzbNMs8axd8Xj9NnvS8eX/P7e2TNI/T+PJMnT65tsKibes3rrCh+u283zSjPZ/+oUaNSvWTJklQPGDAg1WeddVaqDxw4UPJ475q7Y8eOVP/iF79ItX83GTFiRK6xN8LixYtTfeutt6b6vffeS/WvfvWrVA8ePDjVvp/zqHJWR3+P5jL3s9XzVA2glXX5i2iM8X/FGMfGGMdLulXSMzHG2yQ9K+nTHXebLenRho0SAAAAANA2qumaW+6rOtq4aK2OnjP63foMCQAAAADQznJ1zf1AjHGupLkd9euSLqr/kAAAAAAA7ayqA9Fm4edjSaXnavj5DH7+Q5ZGnR/kbdEnTJiQ6o0bN6b6iSeeSPXQoUNTvXfv3lT7eaHVXsakM/W6xETR52Y24tIs9TwvNM/9sy5NU34OEoo1aNCginWec3f9PPCsS0EcPHiwhtFly3PJIM5BagzWa3PK8774du3nPjr/blH+nH7+Z//+/VN96NCh3OMsip/36uez+jrw1+r8+9awYcNS/eSTT6ba18XEiRNrG2wvwb4DOKqWaC4AAAAAAFXjQBQAAAAAUKiWjOZ6DE6STjvttFR7NG379u0VH9+IS7mU8zjLypUrU+2RP4/vbtp07Oo3J554YtXLq6SIyGktz1PPOG6zRHazHp8Vn2z2S7a4dowSrV69OtUf+9jHUu2x/mXLlqXao//+Pvq2PGPGjFT7JWE8EpelfB1nzRu/3WOBHu/2uFzWpRaK1oxzyNerrzOv/X13fvs777xT8T4exUc+RXx2TZ06NdV++ZaLLjrW+sI/x8tj7748j776Jdo+97nPpXratGmpfuSRR6oaa638u5BfpsXnuH+Pcrt27Ur1CSeckGqf+34fZGvG/R/Q0/hFFAAAAABQKA5EAQAAAACFaslo7p49e0r+njJlSqo91upROI+gZKlnbMIjcuPHj0/15MmTU71+/fpUe0SwiDhuO3TH7anIbt7H+O0jRoxI9UsvvZRqj+Z6dKvdNUt315EjR6Z6w4YNFceR1U3S37tt27al+swzz0y1d+PsTjQ3q6tl1vjWrl2b6uuvvz7Vo0ePTvW8efO6HEc9DRkyJNUrVqxIte+Tzz777ELH5DyG7TFaj2VmdVX1qGLW+zt27Nhah4gGO/XUU1O9YMGCVA8ePDjVPo+l0u643k22T58+qd65c2eqPf5bhMWLF6f6qquuSrXPcf8s8qit73d2796d6rvvvjvVvv368wBANfhFFAAAAABQKA5EAQAAAACFaslorsfMJOmFF15ItXeH9Pt5VNYjKPXqoFseNfTYinfR27t3b6r9AvfeXbNaWePrLH5bS1fgnorN9vSya1lGs0RR6yWrI3DWfaTSOb5v376Kj2l0d9c1a9aU/O2dIn2bfeONNyrex3mnyPPPPz/VHlPzbpUe/3MeD/YO21JpR1yP/HqHVo+HnnHGGaneunVrxXEUzfd5p59+eqqLPl0ii382nHLKKan2aO5bb71V8bE+p/29yvPaitCq+5dyjX4dfrqPR6x92zz55JNLHuMxeN8G/ZSbn/3sZ6n2brwzZ85M9csvv9zNUXfOx+v7W3+tvr/170UeKfbY/LPPPptqPx2h1vh5q87TPJ+DADrHL6IAAAAAgEJxIAoAAAAAKFRLRnM96iWVdjT0DojeQdfjV0VEN7O6LGZd9DyPPBHcPPfvzv1atWtukfcv5/G8iy++uOLtPi/bQWdxcO9A+dprr6Xat9m/+Iu/SPXmzZtTvWjRom6PySN0s2bNKvm3w4cPp3rp0qWpHj58eMXn8m6SHl/zSKdv+x5fy4rmZj2/VBr39Bidz5u+ffum2jv/ekTQTwP467/+61T7Oq41IugdcT/5yU+m+l/+5V9S/fOf/zzVmzZtSnXW/jJr31ZPvj3u378/1R7Nzbq/xx/986b8fewpRay/ZlHLaz3rrLNS7dv+li1bUl3ePdu76Pq88VOCvIPuqlWrUn3DDTek2j8batkGfaxSaRz37bffrnj7pEmTUu37YZ/XfrqA79P9tdUazW23eZr1etrtdQL1wC+iAAAAAIBCcSAKAAAAAChUS0Zzp0+fXvK3dzT0KJxH8jxuk9XRsJ5x0KyIYrXRjHrFcZsxptud+7dSHHfUqFGpfuWVV1Ltsb1zzz23puU1G19/5duZR7k8ZurR0sceeyzVHtf0WJzH6MaNG5fqAQMGpNqjst4B99///d9LxuSRVX9fxowZk2rf1jxO6pE6H5/vd8r3VZVMnjw51eXRPL8QvXe19M6ePm6PBfr68Nfw0EMPpfrQoUOp9nXm69gj1VJpF2GP+Xnn4DfffDPVf/M3f1NxrP7++rg7i3c3wqBBg1L9/PPPp9qjtt6x2OPcPlaPOftr+/u///uS5XlEc+7cud0cdXbHzlo/A1pJI7qW+rbl25y/b1LpdpvVPdnjsr6f+7u/+7tU+9z66le/murf/va3qfZ56fuXK6+8MtUf//jHS8bnMXg/nemkk05SJR7l9+f1Uwr89IU8cdzO3p92nqfVbptAb8YvogAAAACAQnEgCgAAAAAoVEtGc8t5B0nvXudx3BEjRqTaI1QevcmSJ0KSN3ZS7TJquX8zxnG7E01phg66eZ/L46hZcyurG2cr6U4U0KOfHgn17THrYuseFfX7e1Q2K9LpkTrvDCmVRih9fG779u2pPu+881K9ePHiVHs37Fq6SZbf//vf/36qZ8yYkepp06aleuHChRWfyztlekTQx+fx02HDhqV6+fLlmWN86aWXUu3r2WN7vh/2+3hc3Zfn243HHIvg88bfx/Hjx6f67LPPTrW/774uhw4dmuoTTjgh1f/wD/9QsjyPbnpc8/LLL0+1dz8uj4RW0m7xv1o/T2sxceLEVPt24/sBqfT9vuiii1L9xBNPpNq7z/pz+Xbnz/utb30r1T6HJkyYkOrTTz+94mM9ci+VblNe+37VT2v61Kc+pUo8IlxLd9yefE8bpd22O6An8IsoAAAAAKBQHIgCAAAAAArVFtFc7143cuTIVL/++uup9jifR6M81psnpptXns6PeTriFt19tqdisI2K7Na6jDzP5bGnrOiSR1HzdFJtB+Wdnbdu3Zpq7xTp3VY9huhdHOfMmZNqj9NmbSu+LXuEzqNyUmlEzmOZO3fuTPWFF16Y6mXLlqXaI6R+cXuP9lXLO3ZKpd14X3zxxVR750vvcOndmX3OeZzPOwV7l/Gnnnoq1Z11+s66qL3ve70Lrr8XfoqEP9Yj2UVH3Pzzw+fHa6+9lmrvouydrl999dVU+2eMr1d/nVLpXPHlrV69OtV54vv+vnid1RXeOxaXr2PfpnxeV9vlvV583kulY8+q/XX76/H5V+1nvL8//fr1K/k3Py3gn/7pn1I9ZcqUVPu+x+3evTvVvl/0OL3PE18fa9asSbXvs3xdSKXbmq8b36/efvvtqfZ97Lp161Lt89dfWx5Zc7R8TC5rnjbbHJVK51Oe19OIOQq0ulwHoiGENyXtk/S+pCMxxhkhhGGSfihpvKQ3Jf1fMcZdjRkmAAAAAKBdVPN/vXw8xjg9xvhBx4yvSZoTY5woaU7H3wAAAAAAdKqWaO5NkmZ21A9Kmivpq1l3LopH2zzO4jEavyC7R9a8q121nU1rvVhzI+K4RXTNLSIeW2QEt7PHZnXI846GHmnySI7HrNpBVuSqs0inx+A94vXyyy+n2js6Xnfddan2LqIe4/R4UxaPQ5V3zfVY3OHDh1P9V3/1V6meO3duxefqLPLbXeWx7RUrVqTau2j++te/TvULL7yQ6gsuuCDVHjnN6rx6xhlnpNojYZ3F3bwjrscVff1dffXVqfZ9r8eLx40bV/F5PFJctKzuwj5HvdOtf5b4e+fdTD1uKZW+viVLllS83+DBg1M9aNCgimP19e3bgW9D/hnoccbymKnPX4+N5tm+GmHv3r0lf/t4fd346/PIuXc83rXrWFDL90HVevPNN0v+9vfR9wvz5s1L9WWXXZbq66+/PtWLFi1KtZ8G4NuWj9VPCfDTCTx+W/5e+fhuvPHGVG/cuDHV//qv/5pqnx8+x30/Uq2sOSpVP0+bbY5K2fvCAQMGpDprjvr+z7d9j1sDvUHeX0SjpCdDCAtCCHd23DY6xrhZkjr+Oyrz0QAAAAAAdMj7i+hlMcZNIYRRkp4KIazMu4COA9c7u7wjAAAAAKBXCNXGGEMI90raL+kLkmbGGDeHEE6RNDfGOKmLx6aFzZw5s+rBVssv1uwRCo9H+O0e+fGunlk8jlPe6cwjmr48j2+UxwS7Us/us1n389fk0VKP6nid54LO3Yn11vK8LquTXR7lz+8xxjfeeCPVHg3Kim6ed9553R5Hs/B456mnnppqf/1ZnWslafny5ameMWNGqv3xHlPzbcU7YHuc3jvu5om9e2RUkp555plUL1iwINUe//Volcdja7m4e3d4PM/jWz7n/H3xSGy168xvv/baa0v+7cknn0y1Rwx9Pfl76hFDj0n6WL0LsD/Wu89OmnTs46WIzro+3/1UDd8v+mvw6KbPDY9LS6Xz2qPAn//851PtHUxXrqz8//v6aSi+zvxzxSOCzjunSqVxSK9riWVWy+d3eQdY/9sjjc7nmceqfV56lLWefFvzSKy/Lz7H/fuBb5sXXXRRl8vy5/Rt8+mnny6538KFC1PtEVkfq++jPfJ7zjnndDmOPLLmqFS6TXlM1/k8bYY5KpWeLpHVwdkj+64n5yhQBD+VSdIC6yuUqctobghhQAhh0Ae1pGskvSrpMUmzO+42W9KjVY4XAAAAANAL5Ynmjpb0cMevCsdL+n6M8ZchhHmSfhRCuEPSW5JuadwwAQAAAADtossD0Rjj65KmVbh9p6SrGjGoevHIi8esPMrmkQi/3aOXeeJr5V0m9+3bl2qPnXg0zWNTWV0qa4mlNqpDbSMuIF2+3Kxl+Hvk8SGPy1TbGTHv6/E55NEgj+S0WxzXZcWl/fV3dmFuj2b5dueRRn/vPGLocSV/v1599dXc45ekxYsXl/ztEVzvoOtj9+03T2S/UbLicn6Be4/C+Rz1dbl06dJU59nePTIplZ5q4LXv2zx+6usvq/aYns+NPBH9RvEosPMOzt7tsn///qn2bcLXdzlff1//+tdTnSeq54/1LrO+DfnzuPJ9pEfOs+K8jebvb3nnUI+W+mez8/2LvwaPxDaKR7R9n+TfAzxmPnLkyFT7/syj7i7ru4bf3+eDVNq52WOg/hk1ceLEisurl6w5KlU/T5thjkql30H89TX7HAWaVTXXEQUAAAAAoGYciAIAAAAACsWBKAAAAACgUHmvI9ryss738fOfss6v6g4/58Fb+fs5aa5e53zWevkWl3UuRPm5KD3FL7fh52c5v73ac1vL7++X+hk0aFCqp037nVOo217WZXHyXi7Hz8H88Ic/3OX9/dIvGzduzLWMPPxcrSuuuKJuz1ukrPO8fN+Wdc5Sd/g5i1OmTKnqsWvXrk21nxeadd5+M8q6bI+fp553ffs+xtfBsGHDUp11jqj3QGgHU6dO7ekh1EWe1+H7s9WrV9dlueXz5KMf/WhdnrcW7TZHpfpd2gbAUfwiCgAAAAAoFAeiAAAAAIBCNXcGqgBFxCzmzp2bar9kQaMvzVLr5Q6KvlxCtTzC98UvfjHVo0aNSvXjjz9e6JjQGNVGQNE8ETK/ZM6nP/3pVPulDJ577rlU54nTd3bprJ4yefLkijVQjv0ZABzFL6IAAAAAgEJxIAoAAAAAKFSvj+YWzSNl1UZfi+ia6zG3vn375htYheesNTqX57nee++9VHsH3WaPFPcWeTvoor316dMn1e+8806q33///VR7F2Xffr1TtUfu3bZt2+oyTgAAUCx+EQUAAAAAFIoDUQAAAABAoYjmNqFGdMStNbJb7f2PO+64VHt324MHD6a6s+imR/WyYnt5xgGgZ+3atSvVzzzzTKq9a+7o0aNT7ZFdj9+fdNJJFZ+faC4AAK2JX0QBAAAAAIXiQBQAAAAAUCiiuQWrJVrak3HcPPfxOO7AgQNT7dE8j+N65NYvYi9lR3j9MQCa34ABA1J9+PDhVPs+JSu+75Fdj/X680ycOLF+gwUAAIXhWz0AAAAAoFAciAIAAAAACkU0twfVKxLbiFhvd563b9++qfbOlx6784vbO79wfflzvfvuu6nOiuZ6LNj5OAAUY8uWLak+88wzU+0x+x07dlR8rEdzfT/g275HcwEAQGviF1EAAAAAQKE4EAUAAAAAFIpobsE8ytqduGw1j61n19w8j/V47Lhx41Ltcbysi9IvWbKk5O99+/Z1uTwXQujyPgDqa/Xq1am+/vrrU33JJZek+vnnn0/1O++8k+phw4al2uO4J598cqp9n+KPHTt2bC3DBgAATSDXL6IhhCEhhP8KIawMIawIIVwSQhgWQngqhLCm479DGz1YAAAAAEDryxvN/bakX8YYz5Y0TdIKSV+TNCfGOFHSnI6/AQAAAADoVJfR3BDCSZKukPSHkhRjfFfSuyGEmyTN7Ljbg5LmSvpqIwbZTjxC6l0gs7q79mRH3Grv7xFcj+YOHDgw1du2bUv1gQMHUu2ddSWpf//+FZft3TKHDj32I/zChQsrjmP69Oldjhv5ZMXKvfb5jda1fPnyVF955ZUl/3bZZZelet26dal+8803U/3yyy+n2rteHzlyJNX9+vVLte8Xzj777C7H5LFeAADQmvL8InqGpO2S/r8QwqIQwr+GEAZIGh1j3CxJHf8d1cBxAgAAAADaRJ4D0eMlfUTSfTHG8yUdUBUx3BDCnSGE+SGE+d0cIwAAAACgjeTpmrtB0oYY4ysdf/+Xjh6Ibg0hnBJj3BxCOEXStkoPjjHeL+l+SQoh9Mp2psOHD0/122+/nWqPNA4YMCDVHl/L0ohoba2Pee+991Ltkb2rrroq1aNGHfvh3KN2kyZNKnkuv9j95s2bK46pPM77AY/molj16gqN+tq+fXuqPdLu3Wf99gkTJqTat0VJ+sUvfpHqlStXptojuH6qgXfB9X3bW2+9leo777wz1R7Zf/HFFyuOFQAAtL4ufxGNMW6RtD6E8MGRwlWSlkt6TNLsjttmS3q0ISMEAAAAALSVvNcR/RNJ3wsh9JX0uqTbdfQg9kchhDskvSXplsYMEQAAAADQTnIdiMYYF0uaUeGfrqpwG8pkdcT1yO7OnTtTfeKJJ3b52CxFRHCzeDTX47Tz5x87PfjCCy9MtXe09a63UmkHTr+Q/ZgxY1LtHVoPHTpU8XlRP76+vc6ao1nRaRTvuOOOS/XgwYNT7V2o165dm2qPze/Zs6fkufwxI0eOTPVJJ52Uao/gbt26NdU33HBDxdvvu+++VPu86du3b6rplAsAQHvJex1RAAAAAADqggNRAAAAAEChyM4VwKOiHnnzDrpDhgxJ9e7du1M9cODAVHtXylrjtLU8Puux3q3Wo8Yeu33ttdcqPtZfsyQdPHgw1aeffnqq169fX3F53vETjbFv375U+3ty2mmnpXrXrl2p3r9/f8njfd54dNNjox75Rf1411zftvw99fdk0KBBqR4/fnzJc/n76DFdj9DPmjUr1R7Tf/rpp1O9cePGVH/oQ8f+P9E+ffqk+pxzzqn0cgAAQBvgF1EAAAAAQKE4EAUAAAAAFIpobsGmTJmS6hUrVqTaY7rOY3Aezc3qrOtx1XqqNsrr4/AL1Htc05/To55S6WvNiuN6nPm8886ranyo3rhx41Lt89Wjl3771KlTSx7v75d3W927d2+qq+0SjXx8G5w0aVKqPYLrsWqP8npsVpJmz56d6pdffjnVS5YsSfUvf/nLVO/YsSPV3hF31KhRqZ44cWKOVwEAANoJv4gCAAAAAArFgSgAAAAAoFAciAIAAAAAChVqvQxIVQsLIS1s5syZhS23FSxdujTVfq6Wvz9eDx8+PNV+qRS/3Mt7772XuTw/17LaOZB1f7/0ht/Hbz/77LNTvXbt2lT7ZSCk0kt6+HmDfp7htGnTqhk2arRq1apU+2U/fJ75vPS5KElbt25NtZ8r6PXkyZPrM1iU2LZtW6r9/Guv+/fvn+rO9gm+PfqlYPz84L59+6baz0PlXFAAANrT3Llz/c8FMcYZXT2GX0QBAAAAAIXiQBQAAAAAUCgu39Ikzj333Iq3L1u2LNUewc263ItfHqUzHn2tJabrPIKbxV9P+WUh3EknnZRqLs3SHPyyH1nWrFmT6i1btpT8m1+mxWOgXqMxfN/h+4hDhw6l2i8J5bcfOXKk5Ll83+FR7Isvvrg+gwUAAL0Cv4gCAAAAAArFgSgAAAAAoFBEc5vcOeecU/H2xYsXp3r37t2p7tevX6q9u2V55Darq22eeG3W/b32+J7XBw4cSLV305w+fXqXy0Xzoytqc6IbMQAAaDb8IgoAAAAAKBQHogAAAACAQhHNbVG1Rlk9zusXpffao7Ye5c3qeDpkyJCaxgQAAACgd+AXUQAAAABAoTgQBQAAAAAUimhuE/KOuB7Bvfnmm3tiOHX1pS99KdXf/va3U/3II4/0xHAAAAAA9IAufxENIUwKISy2/+0NIXw5hDAshPBUCGFNx3+HFjFgAAAAAEBr6/JANMa4KsY4PcY4XdIFkg5KeljS1yTNiTFOlDSn428AAAAAADoVvBtql3cO4RpJ34gxXhZCWCVpZoxxcwjhFElzY4yTunh8WtjMmTO7OeTitXNUtt0RBQYAAAAaa+7cuf7nghjjjK4eU+05ordK+s+OenSMcbMkdRyMjqr0gBDCnZLurHI5AAAAAIA2lbtrbgihr6QbJT1UzQJijPfHGGfkOSoGAAAAALS/3NHcEMJNku6KMV7T8XfLRnObPWrrsVEf6+7du3tiOHU1ZMiQVDfjus+K8krEeQEAAIBKuhPNreY6op/RsViuJD0maXZHPVvSo1U8FwAAAACgl8p1IBpC6C/pakk/sZu/KenqEMKajn/7Zv2HBwAAAABoN1V1za15YQVHc+sVwW3nqGy7q2cUmA68AAAAwO9qdDQXAAAAAICacSAKAAAAAChUtdcRbXp54rhEbXsPf089MlAWH0iyorzlPKbriOkCAAAAXeMXUQAAAABAoTgQBQAAAAAUqu2iuc5jt0RwkUdWlFcqnUMAAAAAuo9fRAEAAAAAheJAFAAAAABQqLaO5gKNcu+991a8vbNOuwAAAACO4hdRAAAAAEChOBAFAAAAABSqraO53gF1yJAhqZ45c2aq6aYL5/Oks/mQNbcAAAAAdI1fRAEAAAAAheJAFAAAAABQKA5EAQAAAACFartzRP3yGVnnf/rtLutcvzzPieaUdc5nnvM6y++T9Xgu2QIAAABUh19EAQAAAACF4kAUAAAAAFCotovmuqzIZFa8NiuuWW2Utzt6Y/y3s0ulNOKSKHmes7P54Jf9AQAAANB9/CIKAAAAACgUB6IAAAAAgEK1dTQ3S7VdTquN8nZHEfHfZtao10mnWwAAAKD55PpFNITwZyGEZSGEV0MI/xlC6BdCmBBCeCWEsCaE8MMQQt9GDxYAAAAA0Pq6PBANIYyR9KeSZsQYp0o6TtKtkv5W0rdijBMl7ZJ0RyMHCgAAAABoD3mjucdLOjGE8J6k/pI2S/qEpM92/PuDku6VdF+9B9gMiohxFhH/bTadvbxnLfQAAAU8SURBVE6iswAAAED76vIX0RjjRkl/J+ktHT0A3SNpgaTdMcYjHXfbIGlMpceHEO4MIcwPIcyvz5ABAAAAAK0sTzR3qKSbJE2QdKqkAZJmVbhrrPT4GOP9McYZMcYZtQwUAAAAANAe8kRzPynpjRjjdkkKIfxE0qWShoQQju/4VXSspE3VLHju3LlVDrV38vhqOyt/ncwPAAAAoH3l6Zr7lqSLQwj9QwhB0lWSlkt6VtKnO+4zW9KjjRkiAAAAAKCd5DlH9BVJ/yVpoaSlHY+5X9JXJd0TQlgrabik7zZwnAAAAACANhFirHhqZ2MWFsJ2SQck7ShsoeitRoh5hsZjnqEIzDMUgXmGIjDPeodxMcaRXd2p0ANRSQohzKdxERqNeYYiMM9QBOYZisA8QxGYZ3B5zhEFAAAAAKBuOBAFAAAAABSqJw5E7++BZaL3YZ6hCMwzFIF5hiIwz1AE5hmSws8RBQAAAAD0bkRzAQAAAACFKvRANIRwbQhhVQhhbQjha0UuG+0thPBmCGFpCGFxCGF+x23DQghPhRDWdPx3aE+PE60lhPBACGFbCOFVu63ivApH/b8d+7ffhhA+0nMjRyvJmGf3hhA2duzTFocQrrN/+18d82xVCOH3embUaCUhhNNCCM+GEFaEEJaFEL7UcTv7M9RNJ/OM/RkqKuxANIRwnKR/lDRL0hRJnwkhTClq+egVPh5jnG5twb8maU6McaKkOR1/A9X4N0nXlt2WNa9mSZrY8b87Jd1X0BjR+v5NvzvPJOlbHfu06THGxyWp43PzVknndDzmnzo+X4HOHJH0lRjjZEkXS7qrYy6xP0M9Zc0zif0ZKijyF9GLJK2NMb4eY3xX0g8k3VTg8tH73CTpwY76QUk39+BY0IJijM9Jervs5qx5dZOkf49HvSxpSAjhlGJGilaWMc+y3CTpBzHGwzHGNySt1dHPVyBTjHFzjHFhR71P0gpJY8T+DHXUyTzLwv6slyvyQHSMpPX29wZ1PjmBakRJT4YQFoQQ7uy4bXSMcbN0dOcoaVSPjQ7tJGtesY9Dvd3dEYt8wE4tYJ6hJiGE8ZLOl/SK2J+hQcrmmcT+DBUUeSAaKtxGy17Uy2Uxxo/oaJzorhDCFT09IPQ67ONQT/dJOlPSdEmbJf19x+3MM3RbCGGgpB9L+nKMcW9nd61wG/MMuVSYZ+zPUFGRB6IbJJ1mf4+VtKnA5aONxRg3dfx3m6SHdTTasfWDKFHHf7f13AjRRrLmFfs41E2McWuM8f0Y439L+o6OxdWYZ+iWEEIfHT04+F6M8ScdN7M/Q11Vmmfsz5ClyAPReZImhhAmhBD66ujJyY8VuHy0qRDCgBDCoA9qSddIelVH59fsjrvNlvRoz4wQbSZrXj0m6X90dJu8WNKeDyJvQLXKzsf7Ax3dp0lH59mtIYQTQggTdLSZzG+KHh9aSwghSPqupBUxxv9j/8T+DHWTNc/YnyHL8UUtKMZ4JIRwt6QnJB0n6YEY47Kilo+2NlrSw0f3fzpe0vdjjL8MIcyT9KMQwh2S3pJ0Sw+OES0ohPCfkmZKGhFC2CDpG5K+qcrz6nFJ1+los4WDkm4vfMBoSRnzbGYIYbqOxtTelPQ/JSnGuCyE8CNJy3W0Q+VdMcb3e2LcaCmXSfqcpKUhhMUdt/1vsT9DfWXNs8+wP0MlIUai2AAAAACA4hQZzQUAAAAAgANRAAAAAECxOBAFAAAAABSKA1EAAAAAQKE4EAUAAAAAFIoDUQAAAABAoTgQBQAAAAAUigNRAAAAAECh/n+xbZDBgJeQ/gAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "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": [ "\n", "\n", "\n", "\n", "\n", "\n", "%3\n", "\n", "\n", "\n", "140704680405368\n", "\n", "Func: kernel (dst,img,w_2)\n", "\n", "\n", "\n", "140704680405256\n", "\n", "Block\n", "\n", "\n", "\n", "140704680405368->140704680405256\n", "\n", "\n", "\n", "\n", "\n", "140704680405032\n", "\n", "_data_img_22\n", "\n", "\n", "\n", "140704680404416\n", "\n", "Loop over dim 0\n", "\n", "\n", "\n", "140704680404080\n", "\n", "Block\n", "\n", "\n", "\n", "140704680404416->140704680404080\n", "\n", "\n", "\n", "\n", "\n", "140704681164528\n", "\n", "_data_dst_00\n", "\n", "\n", "\n", "140704680403520\n", "\n", "_data_img_22_01\n", "\n", "\n", "\n", "140704680403352\n", "\n", "_data_img_22_0m1\n", "\n", "\n", "\n", "140704680404024\n", "\n", "Loop over dim 1\n", "\n", "\n", "\n", "140704680404360\n", "\n", "Block\n", "\n", "\n", "\n", "140704680404024->140704680404360\n", "\n", "\n", "\n", "\n", "\n", "140704680403968\n", "\n", "_data_dst_00[_stride_dst_1*ctr_1]\n", "\n", "\n", "\n", "140704680404360->140704680403968\n", "\n", "\n", "\n", "\n", "\n", "140704680404080->140704681164528\n", "\n", "\n", "\n", "\n", "\n", "140704680404080->140704680403520\n", "\n", "\n", "\n", "\n", "\n", "140704680404080->140704680403352\n", "\n", "\n", "\n", "\n", "\n", "140704680404080->140704680404024\n", "\n", "\n", "\n", "\n", "\n", "140704680405256->140704680405032\n", "\n", "\n", "\n", "\n", "\n", "140704680405256->140704680404416\n", "\n", "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ps.to_dot(ast, graph_style={'size': \"9.5,12.5\"})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*pystencils* also builds a tree structure of the program, where each `Assignment` node internally again has a *sympy* AST which is not printed here. Out of this representation *C* code can be generated:" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
FUNC_PREFIX void kernel(double * RESTRICT _data_dst, double * RESTRICT const _data_img, int64_t const _size_dst_0, int64_t const _size_dst_1, int64_t const _stride_dst_0, int64_t const _stride_dst_1, int64_t const _stride_img_0, int64_t const _stride_img_1, int64_t const _stride_img_2, double w_2)\n",
       "{\n",
       "   double * RESTRICT const _data_img_22 = _data_img + 2*_stride_img_2;\n",
       "   for (int ctr_0 = 1; ctr_0 < _size_dst_0 - 1; ctr_0 += 1)\n",
       "   {\n",
       "      double * RESTRICT _data_dst_00 = _data_dst + _stride_dst_0*ctr_0;\n",
       "      double * RESTRICT const _data_img_22_01 = _stride_img_0*ctr_0 + _stride_img_0 + _data_img_22;\n",
       "      double * RESTRICT const _data_img_22_0m1 = _stride_img_0*ctr_0 - _stride_img_0 + _data_img_22;\n",
       "      for (int ctr_1 = 1; ctr_1 < _size_dst_1 - 1; ctr_1 += 1)\n",
       "      {\n",
       "         _data_dst_00[_stride_dst_1*ctr_1] = ((w_2*_data_img_22_01[_stride_img_1*ctr_1] - w_2*_data_img_22_0m1[_stride_img_1*ctr_1] - 0.5*_data_img_22_01[_stride_img_1*ctr_1 + _stride_img_1] - 0.5*_data_img_22_0m1[_stride_img_1*ctr_1 + _stride_img_1] - 0.5*_data_img_22_0m1[_stride_img_1*ctr_1 - _stride_img_1] + 0.5*_data_img_22_01[_stride_img_1*ctr_1 - _stride_img_1])*(w_2*_data_img_22_01[_stride_img_1*ctr_1] - w_2*_data_img_22_0m1[_stride_img_1*ctr_1] - 0.5*_data_img_22_01[_stride_img_1*ctr_1 + _stride_img_1] - 0.5*_data_img_22_0m1[_stride_img_1*ctr_1 + _stride_img_1] - 0.5*_data_img_22_0m1[_stride_img_1*ctr_1 - _stride_img_1] + 0.5*_data_img_22_01[_stride_img_1*ctr_1 - _stride_img_1]));\n",
       "      }\n",
       "   }\n",
       "}\n",
       "
\n" ], "text/plain": [ "FUNC_PREFIX void kernel(double * RESTRICT _data_dst, double * RESTRICT const _data_img, int64_t const _size_dst_0, int64_t const _size_dst_1, int64_t const _stride_dst_0, int64_t const _stride_dst_1, int64_t const _stride_img_0, int64_t const _stride_img_1, int64_t const _stride_img_2, double w_2)\n", "{\n", " double * RESTRICT const _data_img_22 = _data_img + 2*_stride_img_2;\n", " for (int ctr_0 = 1; ctr_0 < _size_dst_0 - 1; ctr_0 += 1)\n", " {\n", " double * RESTRICT _data_dst_00 = _data_dst + _stride_dst_0*ctr_0;\n", " double * RESTRICT const _data_img_22_01 = _stride_img_0*ctr_0 + _stride_img_0 + _data_img_22;\n", " double * RESTRICT const _data_img_22_0m1 = _stride_img_0*ctr_0 - _stride_img_0 + _data_img_22;\n", " for (int ctr_1 = 1; ctr_1 < _size_dst_1 - 1; ctr_1 += 1)\n", " {\n", " _data_dst_00[_stride_dst_1*ctr_1] = ((w_2*_data_img_22_01[_stride_img_1*ctr_1] - w_2*_data_img_22_0m1[_stride_img_1*ctr_1] - 0.5*_data_img_22_01[_stride_img_1*ctr_1 + _stride_img_1] - 0.5*_data_img_22_0m1[_stride_img_1*ctr_1 + _stride_img_1] - 0.5*_data_img_22_0m1[_stride_img_1*ctr_1 - _stride_img_1] + 0.5*_data_img_22_01[_stride_img_1*ctr_1 - _stride_img_1])*(w_2*_data_img_22_01[_stride_img_1*ctr_1] - w_2*_data_img_22_0m1[_stride_img_1*ctr_1] - 0.5*_data_img_22_01[_stride_img_1*ctr_1 + _stride_img_1] - 0.5*_data_img_22_0m1[_stride_img_1*ctr_1 + _stride_img_1] - 0.5*_data_img_22_0m1[_stride_img_1*ctr_1 - _stride_img_1] + 0.5*_data_img_22_01[_stride_img_1*ctr_1 - _stride_img_1]));\n", " }\n", " }\n", "}" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ps.show_code(ast)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Behind the scenes this code is compiled into a shared library and made available as a Python function. Before compiling this function we can modify the AST object, for example to parallelize it with OpenMP." ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
FUNC_PREFIX void kernel(double * RESTRICT _data_dst, double * RESTRICT const _data_img, int64_t const _size_dst_0, int64_t const _size_dst_1, int64_t const _stride_dst_0, int64_t const _stride_dst_1, int64_t const _stride_img_0, int64_t const _stride_img_1, int64_t const _stride_img_2, double w_2)\n",
       "{\n",
       "   #pragma omp parallel num_threads(2)\n",
       "   {\n",
       "      double * RESTRICT const _data_img_22 = _data_img + 2*_stride_img_2;\n",
       "      #pragma omp for schedule(static)\n",
       "      for (int ctr_0 = 1; ctr_0 < _size_dst_0 - 1; ctr_0 += 1)\n",
       "      {\n",
       "         double * RESTRICT _data_dst_00 = _data_dst + _stride_dst_0*ctr_0;\n",
       "         double * RESTRICT const _data_img_22_01 = _stride_img_0*ctr_0 + _stride_img_0 + _data_img_22;\n",
       "         double * RESTRICT const _data_img_22_0m1 = _stride_img_0*ctr_0 - _stride_img_0 + _data_img_22;\n",
       "         for (int ctr_1 = 1; ctr_1 < _size_dst_1 - 1; ctr_1 += 1)\n",
       "         {\n",
       "            _data_dst_00[_stride_dst_1*ctr_1] = ((w_2*_data_img_22_01[_stride_img_1*ctr_1] - w_2*_data_img_22_0m1[_stride_img_1*ctr_1] - 0.5*_data_img_22_01[_stride_img_1*ctr_1 + _stride_img_1] - 0.5*_data_img_22_0m1[_stride_img_1*ctr_1 + _stride_img_1] - 0.5*_data_img_22_0m1[_stride_img_1*ctr_1 - _stride_img_1] + 0.5*_data_img_22_01[_stride_img_1*ctr_1 - _stride_img_1])*(w_2*_data_img_22_01[_stride_img_1*ctr_1] - w_2*_data_img_22_0m1[_stride_img_1*ctr_1] - 0.5*_data_img_22_01[_stride_img_1*ctr_1 + _stride_img_1] - 0.5*_data_img_22_0m1[_stride_img_1*ctr_1 + _stride_img_1] - 0.5*_data_img_22_0m1[_stride_img_1*ctr_1 - _stride_img_1] + 0.5*_data_img_22_01[_stride_img_1*ctr_1 - _stride_img_1]));\n",
       "         }\n",
       "      }\n",
       "   }\n",
       "}\n",
       "
\n" ], "text/plain": [ "FUNC_PREFIX void kernel(double * RESTRICT _data_dst, double * RESTRICT const _data_img, int64_t const _size_dst_0, int64_t const _size_dst_1, int64_t const _stride_dst_0, int64_t const _stride_dst_1, int64_t const _stride_img_0, int64_t const _stride_img_1, int64_t const _stride_img_2, double w_2)\n", "{\n", " #pragma omp parallel num_threads(2)\n", " {\n", " double * RESTRICT const _data_img_22 = _data_img + 2*_stride_img_2;\n", " #pragma omp for schedule(static)\n", " for (int ctr_0 = 1; ctr_0 < _size_dst_0 - 1; ctr_0 += 1)\n", " {\n", " double * RESTRICT _data_dst_00 = _data_dst + _stride_dst_0*ctr_0;\n", " double * RESTRICT const _data_img_22_01 = _stride_img_0*ctr_0 + _stride_img_0 + _data_img_22;\n", " double * RESTRICT const _data_img_22_0m1 = _stride_img_0*ctr_0 - _stride_img_0 + _data_img_22;\n", " for (int ctr_1 = 1; ctr_1 < _size_dst_1 - 1; ctr_1 += 1)\n", " {\n", " _data_dst_00[_stride_dst_1*ctr_1] = ((w_2*_data_img_22_01[_stride_img_1*ctr_1] - w_2*_data_img_22_0m1[_stride_img_1*ctr_1] - 0.5*_data_img_22_01[_stride_img_1*ctr_1 + _stride_img_1] - 0.5*_data_img_22_0m1[_stride_img_1*ctr_1 + _stride_img_1] - 0.5*_data_img_22_0m1[_stride_img_1*ctr_1 - _stride_img_1] + 0.5*_data_img_22_01[_stride_img_1*ctr_1 - _stride_img_1])*(w_2*_data_img_22_01[_stride_img_1*ctr_1] - w_2*_data_img_22_0m1[_stride_img_1*ctr_1] - 0.5*_data_img_22_01[_stride_img_1*ctr_1 + _stride_img_1] - 0.5*_data_img_22_0m1[_stride_img_1*ctr_1 + _stride_img_1] - 0.5*_data_img_22_0m1[_stride_img_1*ctr_1 - _stride_img_1] + 0.5*_data_img_22_01[_stride_img_1*ctr_1 - _stride_img_1]));\n", " }\n", " }\n", " }\n", "}" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ast = ps.create_kernel(update_rule)\n", "ps.cpu.add_openmp(ast, num_threads=2)\n", "ps.show_code(ast)" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "loops = list(ast.atoms(ps.astnodes.LoopOverCoordinate))\n", "l1 = loops[0]\n", "l1.prefix_lines.append(\"#pragma someting\")\n", "l1.is_outermost_loop" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Fixed array sizes\n", "\n", "Since we already know the arrays to which the kernel should be applied, we can \n", "create *Field* objects with fixed size, based on a numpy array:" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
FUNC_PREFIX void kernel(double * RESTRICT const _data_I, double * RESTRICT _data_dst)\n",
       "{\n",
       "   double * RESTRICT const _data_I_21 = _data_I + 1;\n",
       "   for (int ctr_0 = 1; ctr_0 < 81; ctr_0 += 1)\n",
       "   {\n",
       "      double * RESTRICT _data_dst_00 = _data_dst + 290*ctr_0;\n",
       "      double * RESTRICT const _data_I_21_01 = _data_I_21 + 1160*ctr_0 + 1160;\n",
       "      double * RESTRICT const _data_I_21_0m1 = _data_I_21 + 1160*ctr_0 - 1160;\n",
       "      for (int ctr_1 = 1; ctr_1 < 289; ctr_1 += 1)\n",
       "      {\n",
       "         _data_dst_00[ctr_1] = -2*_data_I_21_0m1[4*ctr_1] + 2*_data_I_21_01[4*ctr_1] - _data_I_21_01[4*ctr_1 + 4] + _data_I_21_01[4*ctr_1 - 4] - _data_I_21_0m1[4*ctr_1 + 4] - _data_I_21_0m1[4*ctr_1 - 4];\n",
       "      }\n",
       "   }\n",
       "}\n",
       "
\n" ], "text/plain": [ "FUNC_PREFIX void kernel(double * RESTRICT const _data_I, double * RESTRICT _data_dst)\n", "{\n", " double * RESTRICT const _data_I_21 = _data_I + 1;\n", " for (int ctr_0 = 1; ctr_0 < 81; ctr_0 += 1)\n", " {\n", " double * RESTRICT _data_dst_00 = _data_dst + 290*ctr_0;\n", " double * RESTRICT const _data_I_21_01 = _data_I_21 + 1160*ctr_0 + 1160;\n", " double * RESTRICT const _data_I_21_0m1 = _data_I_21 + 1160*ctr_0 - 1160;\n", " for (int ctr_1 = 1; ctr_1 < 289; ctr_1 += 1)\n", " {\n", " _data_dst_00[ctr_1] = -2*_data_I_21_0m1[4*ctr_1] + 2*_data_I_21_01[4*ctr_1] - _data_I_21_01[4*ctr_1 + 4] + _data_I_21_01[4*ctr_1 - 4] - _data_I_21_0m1[4*ctr_1 + 4] - _data_I_21_0m1[4*ctr_1 - 4];\n", " }\n", " }\n", "}" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "img_field, dst_field = ps.fields(\"I(4), dst : [2D]\", I=img.astype(np.double), dst=filtered_image)\n", "\n", "sobel_x = -2 * img_field[-1,0](1) - img_field[-1,-1](1) - img_field[-1, +1](1) \\\n", " +2 * img_field[+1,0](1) + img_field[+1,-1](1) - img_field[+1, +1](1)\n", "update_rule = ps.Assignment(dst_field[0,0], sobel_x)\n", "\n", "ast = create_kernel(update_rule)\n", "ps.show_code(ast)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compare this code to the version above. In this code the loop bounds and array offsets are constants, which usually leads to faster kernels." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Running on GPU\n", "\n", "If you have a CUDA enabled graphics card and [pycuda](https://mathema.tician.de/software/pycuda/) installed, *pystencils* can run your kernel on the GPU as well. You can find more details about this in the GPU tutorial." ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [], "source": [ "gpu_ast = create_kernel(update_rule, target='gpu', gpu_indexing=ps.gpucuda.indexing.BlockIndexing, \n", " gpu_indexing_params={'blockSize': (8,8,4)})" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
FUNC_PREFIX __launch_bounds__(256) void kernel(double * RESTRICT const _data_I, double * RESTRICT _data_dst)\n",
       "{\n",
       "   if (blockDim.x*blockIdx.x + threadIdx.x + 1 < 81 && blockDim.y*blockIdx.y + threadIdx.y + 1 < 289)\n",
       "   {\n",
       "      const int64_t ctr_0 = blockDim.x*blockIdx.x + threadIdx.x + 1;\n",
       "      const int64_t ctr_1 = blockDim.y*blockIdx.y + threadIdx.y + 1;\n",
       "      double * RESTRICT _data_dst_10 = _data_dst + ctr_1;\n",
       "      double * RESTRICT const _data_I_11_21 = _data_I + 4*ctr_1 + 5;\n",
       "      double * RESTRICT const _data_I_1m1_21 = _data_I + 4*ctr_1 - 3;\n",
       "      double * RESTRICT const _data_I_10_21 = _data_I + 4*ctr_1 + 1;\n",
       "      _data_dst_10[290*ctr_0] = -2*_data_I_10_21[1160*ctr_0 - 1160] + 2*_data_I_10_21[1160*ctr_0 + 1160] - _data_I_11_21[1160*ctr_0 + 1160] - _data_I_11_21[1160*ctr_0 - 1160] + _data_I_1m1_21[1160*ctr_0 + 1160] - _data_I_1m1_21[1160*ctr_0 - 1160];\n",
       "   } \n",
       "}\n",
       "
\n" ], "text/plain": [ "FUNC_PREFIX __launch_bounds__(256) void kernel(double * RESTRICT const _data_I, double * RESTRICT _data_dst)\n", "{\n", " if (blockDim.x*blockIdx.x + threadIdx.x + 1 < 81 && blockDim.y*blockIdx.y + threadIdx.y + 1 < 289)\n", " {\n", " const int64_t ctr_0 = blockDim.x*blockIdx.x + threadIdx.x + 1;\n", " const int64_t ctr_1 = blockDim.y*blockIdx.y + threadIdx.y + 1;\n", " double * RESTRICT _data_dst_10 = _data_dst + ctr_1;\n", " double * RESTRICT const _data_I_11_21 = _data_I + 4*ctr_1 + 5;\n", " double * RESTRICT const _data_I_1m1_21 = _data_I + 4*ctr_1 - 3;\n", " double * RESTRICT const _data_I_10_21 = _data_I + 4*ctr_1 + 1;\n", " _data_dst_10[290*ctr_0] = -2*_data_I_10_21[1160*ctr_0 - 1160] + 2*_data_I_10_21[1160*ctr_0 + 1160] - _data_I_11_21[1160*ctr_0 + 1160] - _data_I_11_21[1160*ctr_0 - 1160] + _data_I_1m1_21[1160*ctr_0 + 1160] - _data_I_1m1_21[1160*ctr_0 - 1160];\n", " } \n", "}" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ps.show_code(gpu_ast)" ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.8" } }, "nbformat": 4, "nbformat_minor": 2 }