{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# The need for speed without bothering too much: An introduction to `numba`\n", "\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Do you really need the speed?\n", "\n", "* Write your Python program.\n", "* Ensure it executes correctly and does what it is supposed to.\n", "* Is it fast enough?\n", "* If yes: Ignore the rest of the presentation.\n", "* If no:\n", " 1. Get it right.\n", " 2. Test it's right.\n", " 3. Profile if slow.\n", " 4. Optimise (C, C++/Cython/`numba` and save yourself the pain).\n", " 5. Repeat from 2.\n", " \n", "> We *should forget* about small efficiencies, say about 97% of the time: **premature optimization is the root of all evil**.\n", "\n", "> Yet we should not pass up our opportunities in that critical 3%. A good programmer will not be lulled into complacency by such reasoning, he will be wise to look carefully at the critical code; but only **after** that code has been identified.\n", "\n", "

Donald Knuth

" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# The need for speed\n", "\n", "* For many programs, the most important resource is **developer time**.\n", "* The best code is:\n", " * Easy to understand.\n", " * Easy to modify.\n", "* But sometimes execution speed matters. Then what do you do?\n", "* Go find a compiler!" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# A Python compiler?\n", "\n", "* Takes advantage of a simple fact:\n", " * Most functions in your program only use a small number of types.\n", "* → Generate machine code to manipulate only the types you use!\n", "* LLVM (Low Level Virtual Machine) library already implements a compiler backend.\n", " * It is used to construct, optimize and produce intermediate and/or binary machine code.\n", " * A compiler framework, where you provide the \"front end\" (parser and lexer) and the \"back end\" (code that converts LLVM's representation to actual machine code).\n", " * Multi platform.\n", " * LLVM optimizations (inlining, loop unrolling, SIMD vectorization etc)." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# How can `numba` help?\n", "\n", "* If you have big `numpy` arrays with your data (remember `pandas` uses `numpy` under-the-covers), `numba` makes it easy to write simple functions that are fast that work with that data.\n", "* `numba` is an open source JIT (Just-In-Time) compiler for Python functions.\n", "* From the types of the function arguments, `numba` can often generate a specialized, fast, machine code implementation at\n", "runtime.\n", "* Designed to work best with numerical code and `numpy` arrays.\n", "* Uses the LLVM library as the compiler backend.\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# `numba` features\n", "\n", "* Numba generates optimized machine code from pure Python code using the LLVM compiler infrastructure. With a few simple annotations, array-oriented and math-heavy Python code can be JIT compiled to performance similar as C, C++ and Fortran, without having to switch languages or Python interpreters.\n", "* `numba` supports:\n", " * Windows, OSX and Linux.\n", " * 32 and 64 bit CPUs and NVIDIA GPUs (CUDA).\n", " * `numpy`\n", "* Does *not* require any C/C++ compiler.\n", "* Does *not* replace the standard Python interpreter (all Python libraries are still available).\n", "* Easy to install (`conda` not longer required): `pip install numba` (wheels for Windows/Linux/OSX are available, no need to compile anything)." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# How `numba` works\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# `numba` modes of compilation\n", "\n", "* **object mode**: Compiled code operates on Python objects. Supports nearly all of Python, but generally cannot speed up code by a large factor. Only significant improvement is the compilation of loops that can be compiled in *nopython* mode.\n", " * In object mode, Numba will attempt perform *loop lifting*, i.e. extract loops and compile them in *nopython* mode.\n", " * Works great for functions that are bookended by uncompilable code, but have a compilable core loop.\n", " * All happens automatically.\n", "* **nopython mode**: Compiled code operates on native machine data. Supports a subset of Python, but runs close to C/C++/FORTRAN speed." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# **nopython** mode features\n", "* Standard control and looping structures: `if`, `else`, `while`, `for`, `range`.\n", "* `numpy` arrays, int, float, complex, booleans, and tuples.\n", "* Almost all arithmetic, logical, and bitwise operators as well as functions from the math and numpy modules.\n", "* Nearly all `numpy` dtypes: `int`, `float`, `complex`, `datetime64`, `timedelta64`.\n", "* Array element access (read and write).\n", "* Array reduction functions: `sum`, `prod`, `max`, `min`, etc.\n", "* Calling other `nopython` mode compiled functions.\n", "* Calling `ctypes` or `cffi` wrapped external functions." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Supported Python features\n", "\n", "* Built-in types: `int`, `bool`, `float`, `complex`, `bytes`.\n", "* Container types: generators, lists, tuples and sets (with some restrictions).\n", "* Built-in functions: Most builtin functions are supported, with some restrictions.\n", "* Standard lib modules: `array` (limited support), `cmath`, `collections`, `ctypes`, `enum`, `math`, `operator`, `functools`, `random`.\n", "* Third-party modules: `cffi` - Similarly to `ctypes`, `numba` is able to call into `cffi` declared external functions." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Supported `numpy` features\n", "\n", "* `numba` integrates seamlessly with `numpy`, whose arrays provide an efficient storage method for homogeneous sets of data and `numpy` *dtypes* provide type information useful when compiling.\n", "* `numba` understands calls to `numpy` *ufuncs* and is able to generate equivalent native code for many of them.\n", "* `numpy` arrays are directly supported in `numba`. Access to `numpy` arrays is very efficient, as indexing is lowered to direct memory accesses when possible.\n", "* `numba` is able to generate *ufuncs* and *gufuncs*. This means that it is possible to implement ufuncs and gufuncs within Python, getting speeds comparable to that of ufuncs/gufuncs implemented in C extension modules using the `numpy` C API." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# `numpy` arrays and dtype objects\n", "\n", "* `numpy`'s main object is the homogeneous multidimensional array. It is a table of elements, all of the same type, indexed by a tuple of positive integers. In `numpy` dimensions are called *axes*. The number of axes is *rank*.\n", "* A data type object (an instance of `numpy.dtype` class) describes how the bytes in the fixed-size block of memory corresponding to an array item should be interpreted. It describes the type of the data (integer, float, Python object, etc.), size of the data (how many bytes is in e.g. the integer), byte order and some other parameters that further describe the data if it is e.g. a sub-array or an aggregate of other data types.\n", "\n", "```python\n", ">>> import numpy as np\n", ">>> a = np.arange(15).reshape(3, 5)\n", ">>> a\n", "array([[ 0, 1, 2, 3, 4],\n", " [ 5, 6, 7, 8, 9],\n", " [10, 11, 12, 13, 14]])\n", ">>> a.shape\n", "(3, 5)\n", ">>> a.ndim\n", "2\n", ">>> a.dtype.name\n", "'int64'\n", ">>> np.ones( (2,3,4), dtype=np.int16 ) # dtype can also be specified\n", "array([[[ 1, 1, 1, 1],\n", " [ 1, 1, 1, 1],\n", " [ 1, 1, 1, 1]],\n", " [[ 1, 1, 1, 1],\n", " [ 1, 1, 1, 1],\n", " [ 1, 1, 1, 1]]], dtype=int16)\n", "```" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Basic usage\n", "\n", "* **Lazy compilation**: Compilation will be deferred until the first function execution. `numba` will infer the argument types at call time, and generate optimized code based on this information. `numba` will also be able to compile separate specializations depending on the input types. \n", "```python\n", "from numba import jit\n", "#\n", "@jit\n", "def f(x, y):\n", " return x + y\n", "```\n", "* **Eager compilation**: `int32(int32, int32)` is the function’s signature. In this case, the corresponding specialization will be compiled by the `@jit` decorator, and no other specialization will be allowed. This is useful if you want fine-grained control over types chosen by the compiler (for example, to use single-precision floats).\n", "```python\n", "from numba import jit, int32\n", "#\n", "@jit(int32(int32, int32))\n", "def f(x, y):\n", " return x + y\n", "```" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Compilation options\n", "\n", "* **`nopython`**: `numba` has two compilation modes: `nopython` mode and `object` mode. The former produces much faster code, but has limitations that can force `numba` to fall back to the latter. To prevent `numba` from falling back, and instead raise an error, pass `nopython=True`.\n", "* **`nogil`**: Whenever `numba` optimizes Python code to native code that only works on native types and variables (rather than Python objects), it is not necessary anymore to hold Python’s global interpreter lock (GIL). `numba` will release the GIL when entering such a compiled function if you passed `nogil=True`.\n", "* **`parallel`**: Enables an experimental feature that automatically parallelizes (and performs other optimizations for) those operations in the function known to have parallel semantics. This feature is enabled by passing `parallel=True` and must be used in conjunction with `nopython=True`:\n", "* **`cache`**: To avoid compilation times each time you invoke a Python program, you can instruct `numba` to write the result of function compilation into a file-based cache. This is done by passing `cache=True`:" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Creating `numpy` universal functions with `@vectorize`\n", "\n", "* `numba`'s vectorize allows Python functions taking scalar input arguments to be used as `numpy` ufuncs. Creating a traditional `numpy` ufunc is not not the most straightforward process and involves writing some C code. `numba` makes this easy. Using `@vectorize()`, `numba` can compile a pure Python function into a ufunc that operates over `numpy` arrays as fast as traditional ufuncs written in C.\n", "* Using `@vectorize`, you write your function as operating over input scalars, rather than arrays. `numba` will generate the surrounding loop (or kernel) allowing efficient iteration over the actual inputs." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 0. 2. 4. 6. 8. 10.]\n", "[0. 0.4 0.8 1.2 1.6 2. ]\n" ] } ], "source": [ "from numba import vectorize, float64\n", "import numpy as np\n", "\n", "@vectorize([float64(float64, float64)])\n", "def f(x, y):\n", " return x + y\n", "\n", "a = np.arange(6)\n", "print(f(a,a))\n", "a = np.linspace(0, 1, 6)\n", "print(f(a,a))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Compiling python classes with `@jitclass`\n", "\n", "`numba` supports code generation for classes via the `@jitclass` decorator. A class can be marked for optimization using this decorator along with a specification of the types of each field. We call the resulting class object a **jitclass**.\n", "* *All methods* of a jitclass are compiled into `nopython` functions. The data of a jitclass instance is allocated on the heap as a C-compatible structure so that any compiled functions can have direct access to the underlying data, bypassing the interpreter." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "import numpy as np\n", "from numba import jitclass # import the decorator\n", "from numba import int32, float32 # import the types\n", "\n", "spec = [\n", " ('value', int32), # a simple scalar field\n", " ('array', float32[:]), # an array field\n", "]\n", "\n", "@jitclass(spec)\n", "class Bag(object):\n", " def __init__(self, value):\n", " self.value = value\n", " self.array = np.zeros(value, dtype=np.float32)\n", "\n", " @property\n", " def size(self):\n", " return self.array.size\n", "\n", " def increment(self, val):\n", " for i in range(self.size):\n", " self.array[i] = val\n", " return self.array" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# More features\n", "\n", "* **Flexible specializations with `@generated_jit`**: Sometimes you want to write a function that has different implementations depending on its input types. The `@generated_jit` decorator allows the user to control the selection of a specialization at compile-time, while retaining runtime execution speed of a JIT function.\n", "* **Creating C callbacks with `@cfunc`**: Interfacing with some native libraries (for example written in C or C++) can necessitate writing native callbacks to provide business logic to the library. The `@cfunc` decorator creates a compiled function callable from foreign C code, using the signature of your choice.\n", "* **Automatic parallelization with `@jit`**: Setting the `parallel` option for `@jit` enables an experimental Numba feature that attempts to automatically parallelize and perform other optimizations on (part of) a function.\n", "* **CUDA support**: CUDA has an execution model unlike the traditional sequential model used for programming CPUs. In CUDA, the code you write will be executed by multiple threads at once (often hundreds or thousands). Your solution will be modeled by defining a thread hierarchy of grid, blocks and threads.`numba`'s CUDA support exposes facilities to declare and manage this hierarchy of threads. The facilities are largely similar to those exposed by NVidia’s CUDA C language." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Example 1: Summation" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "457 µs ± 19.4 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n", "46.4 ms ± 1.21 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n", "168 ns ± 5.99 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)\n", "165 ns ± 0.0791 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)\n", "Speed-up factor: 2715.0x\n" ] } ], "source": [ "from numba import jit\n", "\n", "def psum(x):\n", " res = 0\n", " for i in range(x):\n", " res += i ** 2 + 1 + i\n", " return res\n", "\n", "nsum = jit(psum)\n", "\n", "t1 = %timeit -c -o psum(1000)\n", "%timeit -c psum(100000)\n", "t2 = %timeit -c -o nsum(1000)\n", "%timeit -c nsum(100000000)\n", "print(\"Speed-up factor: {:.1f}x\".format(int(t1.average / t2.average)))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Example 2: Fibonacci series" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "71.6 µs ± 1.85 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n", "918 ns ± 11.8 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n", "Speed-up factor: 78.0x\n" ] } ], "source": [ "from numba import jit\n", "\n", "def fib(n):\n", " a, b = 0, 1\n", " for i in range(n):\n", " a, b = b, a + b\n", "\n", " return a\n", "\n", "nfib = jit(fib, nopython=True)\n", "\n", "t1 = %timeit -c -o fib(1000)\n", "t2 = %timeit -c -o nfib(1000)\n", "print(\"Speed-up factor: {:.1f}x\".format(t1.average / t2.average))" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "%load_ext cython" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "%%cython\n", "\n", "def cfib(int n):\n", " cdef int i, a, b\n", " a, b = 0, 1\n", " for i in range(n):\n", " a, b = b, a + b\n", " return a" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "645 ns ± 1.52 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n", "Speed-up factor (cython vs numba): 1.4x\n" ] } ], "source": [ "t3 = %timeit -c -o cfib(1000)\n", "print(\"Speed-up factor (cython vs numba): {:.1f}x\".format(t2.average / t3.average))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Example 3: Mandelbrot fractal" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "from timeit import default_timer as timer\n", "from matplotlib.pylab import imshow, jet, show, ion\n", "import numpy as np\n", "\n", "from numba import jit\n", "\n", "def mandel(x, y, max_iters):\n", " \"\"\"\n", " Given the real and imaginary parts of a complex number,\n", " determine if it is a candidate for membership in the Mandelbrot\n", " set given a fixed number of iterations.\n", " \"\"\"\n", " i = 0\n", " c = complex(x,y)\n", " z = 0.0j\n", " for i in range(max_iters):\n", " z = z * z + c\n", " if (z.real * z.real + z.imag * z.imag) >= 4:\n", " return i\n", "\n", " return 255\n", "\n", "def create_fractal(min_x, max_x, min_y, max_y, image, iters):\n", " height = image.shape[0]\n", " width = image.shape[1]\n", "\n", " pixel_size_x = (max_x - min_x) / width\n", " pixel_size_y = (max_y - min_y) / height\n", " for x in range(width):\n", " real = min_x + x * pixel_size_x\n", " for y in range(height):\n", " imag = min_y + y * pixel_size_y\n", " color = mandel(real, imag, iters)\n", " image[y, x] = color\n", "\n", " return image" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "5.82 s ± 88.8 ms per loop (mean ± std. dev. of 3 runs, 1 loop each)\n", "47.9 ms ± 403 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)\n", "Speed-up factor: 121.5x\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAD8CAYAAAB9y7/cAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztvX+wHNd13/k93T0zD6AlUZLLDkNyl1LMIgMsgMWDYNFWRBhAADz0uEynStqSy1VmvNpi1f5M4qRiKq6Ka3/8Ee+mYju1W3JYVrJySmvLYZyIpTfEDwMwyFKVuSCABSBAZMjIWgkWbdkliSIBvDcz3Sd/3Hu7b/d0z3T3dM9099xP1dSb193TffvXt0+fe+45xMwwGAwGQ/uwlt0Ag8FgMFSDEXiDwWBoKUbgDQaDoaUYgTcYDIaWYgTeYDAYWooReIPBYGgpCxd4ItogoteJ6E0ienbR2zcYDIZVgRYZB09ENoD/AOAYgNsALgH4OWa+tbBGGAwGw4qwaAv+xwG8ycxfZ+YhgN8D8NSC22AwGAwrgbPg7T0I4Fva/7cBfFRfgIieAfAMANiwD+zEexfXOkN7oGU3ACirEZG1ECXMJDy6+x0AwBs33wPE38pJ/o7ld/XXD5dLf5Ov+A3fDKTPzRbuYMjbmS6uRQt8UqMip5iZnwPwHAC8lz7AH6Wji2hXSPwGqmQbpm+7SsiqgboXPccJbSd1TVqWPjGYRo4D+sseNi+fgrv3KDAcAQhFm2wLg1sX0V8/gc0rpwEA7p4j4K1t8TDw/ajA+zHVZb/YvmSE49srd+XVrXtJvOL/YeZlFy3wtwE8rP3/EIBvL7gN6RhxbzS1EPZ5mCbu02BfiPv+4wBxsB4iW3yX6+DRCP31E+CtLbA3CsQ9sR266JJVqVCSRdWJvLrfWij0WVi02lwC8CgRfYiIugA+BeCFBbchmarFnSwj7hVSK3Ev6TxnEncA8BnurkPAeIzBtbNApwvYNuA4oLU1kOPA3X0YGI3Bd++BR2Mh7EnintqYaq/dys/fit57C7XgmXlMRP8DgNMAbAD/kplvLrINS2FFL65FsdLiruMz3L1HMbh+Npjk7j0KeD7gecINo0S9hllkK7XkxQZWzpJfaJhkXhbmg6/SejfiXhm1EnZgvnMd25dEgU/zwavvtnDJUKcDQLhk4HniuyeFTRf4mAU/oQVpYttkn7zYQLXrr5hX/D/ED/i7texkXS2MuFdG7cS9RGaK+7RpQKQjFUC6UNcUY8mXh1Ggqqx3I+6VUUtxL8l6L+SaUeuwCIMb5wEA/QMbgSU/wbxv7QvoT1qIT34F7tH27+EyWIELZ1m0WdzTl4mtnyiYRkQyUsYKxB0ANi+fCt0yVbXNXOe1Z7XPUNnW+4pYBcuALGq9uGdyzRAFy5JtAd2OaINFcPcdQ//ABgCgf9AVbgg5L8u6C709VHi9L+R8t/yeNT74smjxRbJsainsQKnnfEJckzpU1b8dJxDuwfVzcPcdE9+vnomuw7alO8YGHAJ5Png4lKNaOdxOnnDJxMZX59Ou3B8fbqiVfvnVVaVFDGoyzE1rxX2a392yxPWpPvp0S1icg5sXApcMdTugmDW+eWkgBN+yMbhxHuQ4QsBmdNYW7gNouiUvNrSY7SwQY8HPSwsvirqwCuIenS796o4jrG8Vu65vmghgX8S3kwVybGxeHqRuanBNxMQzs2y3H1rt+rqVyPs+iCjcbnxU6zQqHDVqLPlirKY6Geu91tTW3w6ULu4TeWakuA5uXgBsO2JRR6xr5sxC1D/oioeG40TXob8llHlPVGT0GEs+P+3Zk0XT8s6ZZVFbYQfm71DNIO5q4JK775gQ+V4PZGuDmWTEjGrP5uVTs7fteeCtLWA81tpjTXayqjeHpLbnpen3RtPbL2nHXuShDEulJSe/brRa3OOr00ei7lgT1rptg7od2TkqrHNybLFt2wY6DtDpCuG37OzCG/Hji4ia4NNxQN2u6Ljt9URb9NGx81CBEbTQa6QF97nxweelBSe9jrRS3FP2KSLujgNYNsjhyEAlBY+9cD1EgU9d+N4zWO+AcOdYtjDnxmPA8oVfXcXP37wQpBI++djHI+0s5Iuf2OFyffML88eLjTXaJ79aamV877Wj9v72IuKe4I5JXr9cxveCqBh33zER9ggZy24RIP3nsG0xDcgu7hCpC6jbweDaWVCvG74R2LZ4GwDA47HISKlCJlPSIMxFmWGlxpLPRHNbvgwafKLrSG2FHSjdag9WG+soZc8DRmO4e4+KlL6+Jz4IQx3JcYK/m5fSI2amoX7HzIBtRfLE99dPaO23IoYQxdw7c2PuoYWyOtkk57XezYVZKq0S94z7EhVLrWMVAKnOTdWJ2u1MDlwqgf6BDcCywMOR8MfHHhru40+Kh4DnBdZ85iyTeSnB9bEwV024wcVuLwGTTbJsjLiXSm3FvSJhB6YMIGIOQiOJpJ+7kzAqtSSSXDv9gy54OBIunPf8EPDOu5E6mhFfPDCfP16nBN/8Qv3xDcQo1yyMuJdKK8Q9q49drTrj22Mgoj4HfvhFwMMR4Htwdx8G37kj3EYxJkfb1j9uvhKa1FasisDXcPj1KtJ4cc8p7ECBof9q/Qt0BQyunRUuGc/Ll4GyJiK/8OuqQbrQnJYumgadxLpT20iZrFEyBYQdmCLu0yJUyAr88JHOz4roH9gQhUF8jpb0S2paYkbKkkW+KfddQ9rZjFbOgwmNXCq1FHagUmEHiiXtIqIgXJIcO3lUaYn0D7oizl7FxMdJeBClivySrfmlXGcNEHnTyZpEA05cE6iluGcV9nk2MU9VJkmeOPfC+L7Ypm0BHkDEmKu7UrW/rA7YnG4q0+E6SbsFvmYFDFaJ2ol7XYR9SgGPwa2Lc20/L/GHiLvrEOBhMstkzG0zEVUTpyyhrzA7ZWnUfKRruwU+L0bcS6Fx4l5CewuJewx392FhTZMVpCRYKEqQ9YIgQDGRByaP6zypDjKK6FKs+BqLfHsVLa/1bsS9FGol7lk67Uqw2osXyYhlk2Qf8HwRsrj/+FztykPRztxCEUJF/fV1vz9r2j5jwRtKoXbCPos52ltI0KeU4AvwWSYCsyob6BTH3XcMGA/h7joEZj/RWgcQKQiik8mST2La8U+zwDNaykvzxdfQkq/nY2fR1PTp2xRqI+55wh71nxHl+uQinnd9Wu3VlPZVSeAKigtijhDPUlILR7YxxdI392ou2mnB57ngzAUzF7UQ9zwDldRPZl0jScKbpzh1SmfqzG0sGHffMYCH0YnKio/74vX5QOR4qONZem6rpA7bDJayseIF7RT4rBhxn4uli3ve9ALIMfiICNb974P/g3fC+qVFBLnIYCef4e45Aup0sHnldP5tZsTddyzIXAmkuFvi7U+p4zp1HRmJlCdc64nNjcZihK06XEFHcL2ENEKN2mYUzlCIpYp73hGPurgrl4ltA51OYuk6VTpv89IA1n07RabH3J32lG61x0vl6eX4NCrP9Mp+KJjxB2DaA0iv4aqOS2yf8rqyJpa3LCHsPoO63WiFKf0YzbgGlm6A1ID2CXzWC8tY74VZ2o1TRNhj4m69/35QtwNa64U1SJVYqY9lgWwL/QMbwoLu9cI86WnXlyZ85DiiFJ4qw9fR/o/lWhft0vzN+v75XFk0jXDNhFWdAESOFSUIt1jGksfHFg+/bjcaDZQg9EX7NJgZtHOHKGe4Y0ezRL4m+lKPViyamhz8prG0nDJFcpSktNP/wTuA4wSl8ZSQERFefP1lUeBaxqIH7oiR9FHrOdzjn2C7IpfM4NZFoNsBOg4Gty6K/ycElMJ0BPo+WtF1qipOZRGERVKCUMqHzeC1l8Ji30q45f4HbScSlaZiRULm6luI/Zbv3gP1uuDxODqvCdZ5DXRm+S0wNILGCDuQGCUTIP3p7p4jolxdrwfrve8RBa0BIf63LmJw8wI2r5wO0vbmHcjk7j+OwfVzwYOkf2AjFHkp7Go7g1sXxQAn25KCKQttywLZuTp3Z9A/6Eq3irTcpUhTrxu+wTjyWPR6QLcTtbhtC4PXXgpG3ZJj48XXX4Z1//uiuXNiD4VMRATcCnLzgBkvvv6y2F7SeairFV8D2tXJmuUmrMFTtWk0qv7lrLYyi/znkCJ87UzwHQCo0wkWdfcfB/xx2DQicFI8eMwtAccR7h+NIC2AEsHYfHIcgChYrn/QLVyebxrxdfYPbITbXD8B9nyR6AwIYvHdPUdEOmEAICv4jZ7qgLe3xfFJirxJiaGfmK9QbzjdDty9RwFvjJOPfgyIZ8rRC4/UqGMzwpLb1a6SfbME3oh7bhYm7mWcm5S2TnQcKreMtJip15sUPil2okYqhzcpWSKqI+G+IRKujTT6BzZE9kZgOakIMqALvo5KKSzeNOyJgVjuvmPAaCiPmTxWebVFO0/U60bOy8lHPybWFy8jGA+FnCKmS01EVqLI5ynZZxTPkMhC/O3KBbMIcQcS47bBLPy7kz8M/cwWgdbWgF5PuAw6TuCfjnQWWiTyySDZb755+ZQoi2dbpfvVyyKxpN/6CcCyQb2ucDsljLIdXDs7mSwtqY8iaZ6+jIqY8fzgvLi7DoW/0/sBDDNpjwVvrPfSWIiwl8mU9k6WmhN+5xff+Arc3YeFYM/A3XcM1O0E1mT/wAZ4NALG48lc6iR86YPr5wrtSl3pr58AbDvRbSTyyo9BlgW+c3d64ZC4O0YtF3dzyX4KWlsDb23JsoJ++Fbg+7mt+KWnEi7Jis9jwRcWeCJ6GMDvAPgrAHwAzzHzbxLRBwB8EcAjAL4B4L9i5u+RuNN+E4AL4C6Av83MV6Ztwwj84qlM3Ks6/jPamyTw5DiiuLTnYfPKaWFNe17EB56FwG2hu28sAiy7ti6YqhC+cj/VfQVEXVjurkMY3LoorHNVTQqxMQGdLqjbAd+9F67XF/lyqNuB/+4dseyKuWkW5aIZA/j7zPzXATwB4L8nol0AngVwjpkfBXBO/g8AJwE8Kj/PAPjsHNuOYsS9FCoR9yrLsBUQd8T85P2DLng4yleLFLIDVg91dJww/HKFozaCyJ8U3L3CYAvcOb0e0OuFIZlyHSpl8ualQRB5JMJYbdEhPRqnb6eu9/sS2lV4i8z8lrLAmfkdAF8D8CCApwB8Xi72eQA/K78/BeB3WPDHAO4nogcKt9xQGqX728v0radRtL3MOPnYx8F37ohX/7v3REcq+7k6BQdXz4SDpIhAnU5gtS8qE2RdEKGXWpinbaf2UYA5CD0F5LFiX44sdsTHtkV/hlpOhq4CEA/nOSzhpYdMLljkSwmTJKJHAOwH8AqAH2XmtwDxECCiH5GLPQjgW9rPbstpb8XW9QyEhY817CyhcTV9mteE0i74RR7nedvs+2AoX69Izwsrf2qAwdUzQXjlQkrs1ZTNSwMZfRO6pdw9RyDKQ2nYtuis7YahqP2DLsjxwJ4fPCD7BzbAvh9Y6IOrZ2RaY5mXxhvPfhjXNWxywcwt8ET0QwD+LYC/y8w/mPJ6ljRj4iwx83MAngOED37e9hmSaaSwA5nFfWaUhe+DLQvkA7AhhKdAZAZZ1kqLuyJ+DAY3zgcRRQBAvW5i4rSkTtvE4ynPexlBIUuv3brAh89cAk9EHQhx/wIz/4Gc/OdE9IC03h8A8B05/TaAh7WfPwTg2/NsXzZiyjxjvScxt7gv67jO0+5pIyotKtwpasQ9GXf/ccD2hKVNlLuPI1iPGkncFX0ceOddYRWWOLq3zRS+U2VUzOcAfI2Z/5k26wUAT8vvTwP4kjb9F0jwBIC3lSunEoy4JzKXuFftV59GVeJuqITB1TOiM1SmXSjyAO0f2JAZLz3w9raw9rsdvPjGVwBo1nzatVFnDVhQ2+YJk/wbAF4GcAMiTBIA/hGEH/73AfxnAL4J4JPM/F35QPg/AWxAhEn+IjO/Om0bmcIkm9aTviTmFvZlkrPtidEz4czIdLKtyQE6htJQic2K5rUP0iQol4ptg7e3o3HwQHqZv7qGTIoGFPrZQuLgF8FMgTfiPpPCwl6XY1ig/fHc4rGZkXkqc+Tg5oXK8r8Y8tNfPyFGssrBZEEMfNJAJyBd4IH6DnwSjcj9kzwC365kY4YIhcS9LsIOlC/ukXWrBFii0LW79yhAo9zbM1QDe7GwVd3nHhf3giy9s3UBtE/g6yRQSyK3sNfxmFUh7rGkY5HlPR+w2n2zNwmyrTDEjli4aWrsbShMxRE17XPR1FGsFsiqijuQ0TUTz3mitqelGSDbqrQWqiE/qdkkFQVcNGJWDfQvp8CvhovGiHuEXMJe5+NUtrgnCHtkeV3c1TSTqbA2uHuOiERj7E1fUM8Nn4NauGkqtOJrfKcbspArzcAywxxnodclzUmiuKsUtJaFF9/4CshxwpJ8cXG3LQxunAd1OiDHNrHtS6a/fgLu/uMiln40Dv3vRWPf63rNL4DmWvArTmssdmCuGPcJcY9Hyajh7jcviLhq2hbpCQCArMgISyPsNYEI2B5OzUxZ6uZabMW3R+DrLmIl0SphB8oXdxX66HOQ5bG/fgKbV06LakoWidwzctvG114fRIlEBkb3oqGRCRDRpPgXdNO0mWYK/Ar6SI2wR0nKG/PiG18Jik9gPBZZCPXMjhYBsAGHghqj7v7jK5f9sW6o5GLwPBHNtKpUYMU3U+DjNEHQClBKRMysdSza4qlA2JX1fvKxj+PF19MHKpHjRAYyqVGWhuWyeflUkJVzKpY1UXaxLBdOLdw0FdBOZWw4hfKz6+KuOiyzrENftspc2XOuP8gpPrHe6CXs7joUFJWIMzFK1baNi6YmkOOIMofaNVJq3dWWGoGzaGYcvH7iW3LiSkkpEKumU8q5nceqKeGBMfUmTwiHpG43U51VQ30J0gw7DnhrK5qiAJjwy2eOiZ/h/qiNBT+jne2Og2+R/73UBGBS3GnHDvBQRCAocZxL6JdQASeT5ZaWhsDz4O4+DNq5w+SVaSD9AxtA1w/6Rdxdh8L0wEH2yGpcNW1007TD/G0Yc5fISxB3BQ+H4ottixtBCxWsM5HSbvPieaKjFVIwDI1h8/KpSKd3aqbPtIFrwGrXxI3RPAtepyHumUqrJ8l1qwucVE1Ln2V5M0/UyJQhZ3Vwyc0t4lPSEDCzWP9wBHfPEVBn+ftrmA8iAiurXb9+9etAlvirw/U9NyVG0zRb4GtO5WXxEtbPwxFg+yIWfMcaaHsbIEtY9pYF0l5tq7gZKn9byFq8g0Wt1aKVhAz1QPjjtTQFCdcXdRxxnj0vCH8tStvcNEbgS6SSiu1F3lJUUentbaDTFR1VCs1/2QTXTYQilZnYF+GQtm188g2if9AVxgrGon/Jh7DigWgnqyUGtdHaGvju3dCKb/qgp5Ks+Gb4OBQ1jJ5R/vTSxX1W3pjY9iLWuBoBSFY090rw2zmPnfTtw7aBTkd81DT1KZu0dcZTANsycZhtRwo1G3FvFpuXBmGZPxLXGnU74lrWrjMica7JtoJSfoaQeqhkw6hM1MMNlLOebkfEee9Ym5xXVITjnVuOI0TVtstZf4Ztps0nIgxuXRSfG+fDeqBmpGpzCR7YMiHczh2B0aLEfXDzApg5CK9s3JtphRiBz0jloi42kk3cZ7VBWTmeJwb93NuaTMQF5LO248uqbbCPwa2LsN77nmSRn0foZ/1e2x9h2VEwOrV/YEMIgWOn/drQAMhxMLhxXgyCgqjtOrh1Ub49OoAjvcz3toJ+prm3WZconBIMvWb64BfonlnYyS5rn1TCLSLhqrHsyfkTSZoSfJvxebF1kLrBJCosMff605h1o2qumaTC2SYzZDtIGmns7j8O2LIztS5iXFOaY8Ev8LVrIdZ6dIPZl53RJpX7HL0eqNMRHTVW1E8d5ErXP2r+NF+6lmMdti3yp3e7cPccAba3J5fL+8aQ1YeviztRaMUZVoLB1TMYXD8nrPurZ4TgdxwRImzcMxGaI/ALYKGiLjZY+tvIyUc/JkLFRkPwaCRebXs9rchFioDatngw2HbyAyB+47AvlmUOohWCDrDIPk4R+iIds0k3sAmFXEmUdT+4ekb0uXS60QXmyX1UlzeDOfWheQJfgXtm4cIuNlru6pTwMYM9L4gLdvceFYJvhSNFSblYYtEIg9deklEoM6xsIAg7ZM8XHWEUdnSmFrxOe1jM3rnJ3+nbYB/uniPoH3TzrdfQGtzdh0U4MHM7BjuVRPMEvkSWIuxiw8V+l7WtWqGEwfVzoF5PuDGU39y2hVXvOFLshcvF3X9chFamibQWrQLPE7nUr50V29i5A+j1xJvAzh0TD5DgQWLbYt60N4VZD4OktvkM3t42KYBXCHffMfF392Fxfft+8bJ+LWUlBX5pwi42vpjtMANkiWRN29vCVaMiEm5eEDHGjoPBay9BFZ0eXD0TZGJMs/IjIWjjMdx9x+DuPx7GLXc7Iiwx/sZAJFxA3Y74AJXEyzOzseRXAHf/cXH97TkiwiTv3p17FGsbaU664KDi/XyisFTf2rzintL2xKLTQBC2qEIIg+07DqjbwealQVBNR8WK9w+64Lv3xPKWDYyGkyMC9XVJ6z8p1tzddSjSlsGN82Ib6yfAd++Fr9JZra60BFMWBVE0/YOuGdS0AkSuIZWjRl5HEU3Tr90CI0Nrk7ZAa3t70wXPIZBL7zRZxshb3w8TcPkALDG6Fb4HHovjEQ8n5PEYsC1QpyPKqJElfqcudF3cZ8bja/ts2Vo5vVE4/Jw5kj4hlTRrX41alRhxbyeRsn4s6rbOFHdDQwR+ztCnRlvts1Y/69jERR4eYInwMuXK0EVxcPWMKFIthd/dc0QUqbYxIejUEa6W1KpIjgP4XnD+9O2oJFJBzvpYZkAAU1046q0kNZ2soVVsXj4lAwaEyEfEPY26WN9LpNU++KX62kUDyltXwn5kHpIdt449D+6eIyKZU4LfMiLYjiPdMHKoeKcT6fycWvJOuXmkj1/PzT64eSHMF2NRNGfOjLBJfb/dPUem7LihTZDjgNZ6y27GciioJc0R+Jw7uJIuGSC7P9tngP2ZNUnD2HixP5tXTouoGSX0M347uHY2FO+4aDtOkEgKnW6QNCpe/CPeWQsg8sBz9x83HasrwOaV0+DRKDqx4oFNS9eROWlGJytRLsFc+kkpW9zzWu/6iFVt2kRna7cDsqyZw/r7BzbEq3DJKXf76yfAno/BtbPCP7+1Jd4oZr1aq+Mh3wxg2WHmQUOrcfccAUbjaAd9lg5WoHD63bp1tLa3kzUDrRP3qrCtTOIOVJjXxbaDB9XmpYHwscY7deMknF8j7quBu/84Ivnhs3bQrzANUaNsrIq4Z/K9Jy1j26JTstsBOc7SE3JtXhpE2jC4fi4cFav55yMfhYrdl6GXhvZDRMI1KK8Nsq0wD3wVNQhaQGuOSmvFPR65MkvcU0agiuySPtx9x4Rfva7I/N5BJ2zScaXQBWWKaq8Om1dOiz6gG+dF/43jrFYUVQGNmVuViMgmoqtE9GX5/4eI6BUieoOIvkhEXTm9J/9/U85/JPtGpjezteIe30zWDqW0gU8AMB6Dt7bA47F85a0X+oCrAJWUTU/Opqo1mdfz1URmMwViA+oUdfGbL5kylOnvAPia9v+vAfh1Zn4UwPcAfFpO/zSA7zHzjwH4dbncbGZFarRZ3LV9y12lRgs1nPhtAy7+zUsDce6Tzq/WyWqq96wmg6tnRGz8vmPVjzVZtsbMwVxHhogeAtAH8NvyfwJwBMDzcpHPA/hZ+f0p+T/k/KM0593Z5AM/kyLirqo2aVZ74m+1SJo6uzjIcaTP1Zr0w8skabPCPA3thmwrTDRmmGDeR99vAPiHANTR/SCA7zOzLO+D2wAelN8fBPAtAJDz35bLRyCiZ4joVSJ6dYTt+OxwuWWLewW53BM3k0fcAVC3K0vVOWEBajWcP+i4tEGOXYuO1mkon6s41vEsk5ZJLLai6BlDeTwG37kjvtc45HtZFO5tI6KfBvAdZr5MRD+lJicsyhnmhROYnwPwHCDi4Iu2r9HkeXjF66R6HgavvSReXVUZPUBYORYJ94xFtRb2OGRbYA9hHLOWe97knlkd3P3HxTiJ8Rbcx5+U6Qq2lt2sWjOPCfoxAD9DRN8A8HsQrpnfAHA/EakHx0MAvi2/3wbwMADI+e8D8N0iG66F9V4FmgsiOmozpQJSvASfzDkDiNhw6nVFvHuvC3JExsfBtbOJmR/rzOaV06J4th5VY6Xkije0lsHVM8Jo8Tkx0VhAA/qYCpNTeworFTN/hpkfYuZHAHwKwHlm/nkAFwB8Qi72NIAvye8vyP8h55/nAu9UrRZ3tQnlR89Tn1TD3XVIVLjxRFm9zSunG2WxJ7F5+VQ4Elfmmjexz6uFKvARoMmHcc8kU8Ud8ssAfomI3oTwsX9OTv8cgA/K6b8E4NlMa6vTyNAq2hIbwENEoPt2RhJ6JbdlynyZZwa+16oiCJtXTodvNsY9s3IMrp0NBsGRLe9F07k6lVJGvDDzHwH4I/n96wB+PGGZLQCfnGc7rcgMOatoh2WBt7ZBaz3wcAQiAmcR6aSQSJ8Bx075QUOx7SAFsmEFUX1J4zGA0eT8NrtnClD/ZGN2OBinkWX2MrQ54msXE0BrPQxunBcl9zxppUyzVqxYTLgKhZxScclgaCJBNafxGPD96QnGdAomGwt/Xg+tfMU7kznZWGPMoNqLe1LelBltnqhxquN5IkY9UhUpxScfF3fVbq3D1mBoC5tXTou3WuOemUmNk5LUgGniPscDJ3nw0eS2eDTC4OYFUf1IummCDHr6ulT4o2oXicEf5NgAsxkMZGg8orKY6Fvi0Xh6JSdDQCMs+KVY70nintEyn7rajNY0e34Yx25HrfHgY1sYvPaSlhddy/XOvhF3Qyvor58Q4u55E+JeZxdzHWiEwC+cuLjPKerADHdMGj7D3X1YptG1J1PoqnaqFLvxRFxmpKehBWxeOS06V3eshRONeyYTRuDjRHze8ws7kDPFb3xZ9uHuORKmSO10QzFXdU5VZ2og/Jrv3dwIhhYwuHYW2N5Or94EmAiaBIzAK+IWcEnCnstqT1rWZ2HJ7z8uRqFeOxvkwlbLD66dFda944BdjB+EAAAgAElEQVS63dDat+3GD3AyrDYqpbW7+7BwzxhyUftO1lT/e97QxaQQqbi1XiJ5k4RNG9REvW7El07dDjYvDdA/sBEIuB4GGWSINNa7oen4LPLOeMPktASGqdQ+Dv6JzonoxIoKWitBLuN45LLaVYrf+DpsO2yfbZvC0oaV5eRjHw9HZKe5Z4DZLpoVjIOvvQUfkFfYM8Sgw5K1P7WY2iRxzir6uYVd/Ch5nlbUAj4D5KO/fsJExRhWApV3hre2QqvdWO+5aYYPvshgo7RVxfzitLaWumz8N2kCPjFvWvZH9SFlmWuDl+T3oIiBFiVDjiOWNxhWgMG1s6G419jLUHfqbcHncYtnsdgnJ4pYc2XNZ7AQJtaTJyeKlpKAZGgjq1j32GhUchww+eK10rZNci3DyuDuPhykIQjEXbs36+xWrhv1t+CzWO9FxF3heal+8KkWuG2D1nqzsz6KBkwuI8McybZiud8pEH5ybAyunxOFLUw0jGFVMGkISqPeFnwWE76ouAfuEQKRHVjKkcyNSZaCvj7bBvksfqOm679JeWgQEcA+BjcuiiHYxBPVivR1Gb+7YRVwdx0Ki3kAtXPNkEW16WjNSs0FHpMCruVcSXOXsOclF4RIcJHQWi+wjt3Hnwx/4/uTAh3PASPzwBDs4KIk5Vv3vMirZCRnjLTS++snAMsH4AHkCLfMeCza5ZgYdsNq4O45Iq12rVJTXhomvIui/gKPWDpdG6EQKnzpp2YWqXYdB6Rb48yRPOJEBDiOSAGgY9sg9sXF0nGCC46IRJ3Tx5+U7RAiPrh+Dv0DGyJ/uxaCNbh5Ae6uQyDbCS5cNRAJ47H4i9Ayd/cfF24YY6kbVoj+gQ3waGTcMRVSax98YD/HLXG9HqeymFOXtUDdTvhdWdJJESlW1P0CQFj6FsHde1TMV+uS8Fg+RNRIWLLE6DuZ0RGAeFiodespBdR+qnBNg2GFYCPslVPvgU7WB/mJ3kkA2kAkWfwiKNmlo7tXkqJb4i4a2xY5XiTu3qPg7WG4Dn2dE6uSYq8XAVbTez1gezs6TX94SD879bpiGfmmEbyZ2DbIcUzUjGElSHXRRAp5zIiiyeKimXOgk1jF8vUyz0Cn+gt8ZyN9gZgfPmlfIu4d7SKhrhBXUnldxmPhN4+Xx4uvU1ufdd9O8NbWZI4M6SpKQh+hOrh1McxzLacF3x0HZFvGbWNYGSKdrPH49yxhkhWPZBWrWL5etqii04yD6TPYEyW70k46J8TRAgAPpaXuecLaHo+F9ex54SdpncxBToyIuKvpegRAwu9Z63x19x2Lbm80Di9C7WI0KX8Nq8Dg1sWoq7JmlcjqIO55qbnAl0P6E98HOo6w2pNib30//EyuFDwcpQq5vlza/MG1s5NWh6xaozqA3b1HweNxmEDMYGgxg5sXZL/XSkhT5dT/KGZ5rZIpdaeuJsnKZwYsaTHM6vCRQs/SjaN/pj4I9G2pj+eBPR8nH/t40C79E/xkJB8gPmPz8iljyRtWgsGti3jx9ZejAxDjIcqGTNQ7TFJpHfvZRrTqvuy0Vcqwx3BCtlwX0/oqIp2peSIDUi5UhvRHwpP77cHddwxke4nLGwxtQ1zvQyEB8X4xQ2bqLfA6WUUemG3NI/ST8zvvzpVBMv6bzNaFHukT/40+yMryAYiRtptXTLpgw2qgp8Y++ejH5Nt2tnxRhpDmCDyQPJx/HuSDgMGlFfwoJPIJUTfBenwGrXVyRdP0D7rBuk2opaHJ9A9sANYdk1WyIM0SeEUJ4U4RVM51oBShj7hsZjFF5IP1bQ/h7juGwbWzYvTf2AN1O+DhKLB01GhYEIH9sdyt+nexGAzT4NEIL77+sqzq5AVWPBGZrJIZqL3As8/pZfvK24j4qws9MLfY57LmxQ8S8t+IQVKDa2fFCFnfA1RIp0Xor58IBktNXO6mM8rQcFQ0Ge3cAb5z1/jjc2JMPJ34m0EJca/TYvTD7SS/kajfDW6cF6kShiMxIItlnnifhevGV9E5vvjI3zGzCa80tAJmFpE1SUkEFVUbgg3ECHwcJZ4KFYI5p9hnEnqxYPCVbAuDWxeFT10fhKXaow+K8rzwfz90EZmMlIam4+4/DozHYbI/iQmXnE0jBH4pI8iS/PwlCX0WgrTDAPjuveTt+gx392Eh/sEGQnfT5pXTJnbe0HgGV8+ILxaFsfGmfykT9c5FQx/gj9rHAaB6P3waJVSUSl11Wr55IlDHAa2ticFOo/HkgyGWDC1aE5aCMn8AQL2eiaYxtAZ3zxHw1nbgj4/cG9MMsDmDM+qSqiBPLprad7IqFtLZmrzh2SKfdOIztHVqtI3apl44JMlXL6exXtPVZwSDpGwrmjvfYGg6yoIvWhxkhTDvOVko8uTXffczXDuJAu44oqp8rzc7DUJ8PTHaFC7ZP+gGH8Pq0V8/IeoU79wRTFuEL74u1nteGnXnL/UglxF7P0Xwg05Y3xfVqO7eBY/G4HfejeaxSUteNuUB0KqOVt83RZlXGB6N4O4+HPrl45hImgiNEvilU/YAqzShHw6zZalMnMxBzVfY7aoU1T+wAfZFemgwGyt+xXD3Hg0iyU4+9nExsUVvp1Uw19EhovuJ6Hkieo2IvkZEP0FEHyCis0T0hvz7frksEdE/J6I3ieg6Ea1n2khMVJv6qjSVmNAHBQ9mWalpDwHHAfW6oE5ncl5DUeKuKFyc2dBYJmooa64ZEzKZzLydrL8J4BQzf4KIugB2AvhHAM4x8z8homcBPAvglwGcBPCo/HwUwGfl39wsrcNVbFz8LSsfjo4f5sQJc9FMETFlvagRsNK9o9rGw2EYSuk4IMcORDL1FbeG9NdPiGHq8WNPLXzYG1Jx9x8H7DEAD+TPLAfUTnJ6EQqrFBG9F8CTAD4HAMw8ZObvA3gKwOflYp8H8LPy+1MAfocFfwzgfiJ6oOj2W0vckp+6rGblq2UtAtmWCK9Uoig/PBwGA6Wa5N4QLhntwpbf6xziayifwdUzGNw4D9qxA4PXXsKLb3xFlNyMu2mMHz5gHjP0wwD+AsC/IqKrRPTbRHQfgB9l5rcAQP79Ebn8gwC+pf3+tpwWgYieIaJXiejVEbZTN750V03Z/nidPCIPRESeR2P4d+6ChyNRoFy5f4ZS8D0f8D3wcAR3//FapzIIHkKqfGLkY9wzq0okaGDHWtDPZNw0k8wj8A6AdQCfZeb9AO5AuGPSSDr6E+rFzM8x80eY+SMd9KY2wIi8vrxmyWepQK/aXuNOKh6ORMcaJ0QeyWlNehMxlA85DqjbqfV1vEzm8cHfBnCbmV+R/z8PIfB/TkQPMPNb0gXzHW35h7XfPwTg2zO3wpz8aFCzl+mPFw2oxh8PJPvkpy6vFRHRRrpO/FY7XnUd4do/sAHwSLxx6Ohpndk3g7hWlP76iTBCbFTtNbB0Q3IOCisTM/8ZgG8R0WNy0lEAtwC8AOBpOe1pAF+S318A8AsymuYJAG8rV87sjdX8dXxB7cvsc44M3Y61zSIMbpwHbLuWr7TKZbR5+VR0P5KSwBlWFh6PAd/H5qUBBq+9NLmA8cMDmD+K5n8E8AUZQfN1AL8I8dD4fSL6NIBvAvikXHYAwAXwJoC7ctlSUE/YVlrymhUPZLDkdStem8aWBYIQ9/5BF4Or9bTceTSCu+eIdCMl5P6ORdI0KRrIMB/q4c9jD/DG4PE9uLsPg9Z6AH6w3MYtggKGZP2TjdFREQKYUTyXKvKiAeWvM7ZPMy3vpFqvKleNloiMOp2lj3LtH9gI2iAGsmiW+qwcP3Iw10R8tKGVqLTBKhJMGTuq4A3UIDhF/Pop+KZdGxeNbP8r/h9mTjbWup6JpZ+MBbhr5n4oqwIhGakq0qa/fiKIy+8fdKPRMWnnceKmZeGPNbSewdUzQWe7ugc4qJFQczfukmidwAMtFPmE/Sks8iqCxmcRK58F5sSQyqzC3z+wkRjtwp4vCjnsPTrZxmmo+Soe3vNrHe5pKI/BzQvhP0bUZ9KYdMF5fdytjq6ZRtZwMXls+usnRNm/FNz9xwFPRClsXj0jR5WKG4vs6UIc1IuFGJClu2MAiPqyKqe3iobxp1e+iqRElvugCo8b2o94W7sXTqjYxbx0Y3FOWmnBK5Z+cip212S24uOiT5YYEeg4icnIItawGhzl+XD3HpXRC54YLDWj5it7fqTMoC7u7p4jkfKD7PniM2Of4qUPBzfOz9h5Q1tw9x4Fbw+X3YzlUFBLmmHBq1wrRX667Aibsiz5WDRNZuIVn2wbsAju/uMYXJ203PvrJ8Bj4bpx9x0T1ru6uMaxXDDjMZjUuqLRLP0DG8BYjkS2CLDF9pM6ysJ9zJB3R8LMIB8iz47paG09/QMbAHmAxQBTmIuGKN2Kt2jlw2lbbcHrND6XfNqqM1jxkegZMSFSBMTdfzy+Urh7jgjXh97pGS/2PeWYBta6Kgiut0X+NrO4T5vP/uRgKEPr2Lx8CoNrZ0Fraxjcugh0ZA4adR0nla40NCRMEggt+BKs4aX65udp/5R2Bxe2JtwiEZMQdhX/ro9cdfcfF2FmoxFAFFjB7q5DgXjr10dq3VfbDix4Zc27uw9HhB0WAY4TZrrc2hLT83aUxSz54OFFlsiYKXPg13WErqE83N2Hw4ABOfAJQHqoZE5Da+kuXkWs3XnCJJvhoimZpbptFtX5qoSv1xOCu/cowFthxApZwHgEjg0ccncfDvzrAGakPLAxuHY26odXnaV6LVkA8Bjkj8RDIXgDKFIKMWEgl4J9MBMGRtxbj0hVMMbg+gUx2AmIiLxBsJIWfJylCH3R/Uhpa8SCl0WJybZCq3kkO6f0XC66hWIRaG0NvLUls1Am3Cj6a7BtY3DzgoyykdkeVacpMhbkyHLtpb1y622RFjz1ulMjggztpH9gAzwcgmXGVNbfHFfcgm+eD74Cfzb7vPiTWfJ+BA9qJe5KhK+fAzl2GKkio1B4NAbLm0F8/OniDkRfgfX2q2gZ6VfntIEns+rKJu9Y8vL6+snC4OYFI+4ryublU4BlC5ckjB9ep3kCXyELF/oKHlbU7YLW1mQnlC1Cy7a2Z3do+n5gASUW+Y67bDA50CnxbXBWofA8r9Rp67HIpA1eQdT15+4/Ho6F0MV9jnu5rtZ7Xprjg58jVDL3phbpo89bAnBGuCQPh8B4LNwzgHC7vPOu9vtYBajIutPdMpHf+D7Q7cjOWWE9i5DKrfzrT5o+a7CWVqKQ5duK6VRdHUQo7xjACP31ExhcPQ338SfBPK584FPTaI7A6yyoo1J/ilcu9mXtk/R9MwCy7TBaBRBCmkd49Xkxoeet7cDP7+4/DrJJxCXnXX/WbcbRRB7DEdzHnwQ6DqjbxeblU8Kil+tZdkI1Q3GSRlpvXjktggEsCkZVAzAdrAkYF01GFuK+Kctlo/njBzfOC3dNEXGftpzeieoz/Dt3o2GR6jfz3HQZO2mD0a2jcSDmm5cG4LEXJDMzNBMej+HuOSLeECXurkMT4yvQ7YiC82063yXoQXOiaICoi2YZeV5iVGrVz9q/hG2T7GCV/0h/fA98567o+MwzsGjqtrVtyMFTpa4/bXtJxOL/ybbEsVPHhwjkOKYDtqEEtQEskd6at7Ymfeu2La7z4RB85y4AFcVVLF1w3f3v7Y2DX6AfPguVunBmuWxmpS6wLBEXPhxOdn7OK75pNV+rsp6muWx0V41lyePiAx5kvniRLzw+yMtQb/oHXdHp78tUGSOGP5TZT/XrzLJAKmGdnF9no3XRNEvgdZaVrTGF+FO/FMHP0QEbCQ2T5fgGty7C3XUIZFuhqyJBhDNlb0xjka/E0wY5SVSOmuDhp42yNTQH9TAWA++0tBbx6833AdsG37sXHVxXFyu8KCW5a5sr8DWn1EictIdZkhVv27B+6D7AsmRkgS9EWnW+5rRuyraGFhqjvOzqXob56XaAe9LXnhKhFaaanl8Ua+OeKYn6mMBFqHsxboSds3N30saLTivi6/R9+O/egf/9t8WF73nCTRNPHbAkgoFWsU9mkjpuNetOPUAGty6C1tYiSdUMzWNw9Uw4KjptwJv81OH6rhvGgl8wc/vtk6x5n8GWsI6p44jQMS3KpQkXfmpSsyLIjl/jc28e/QMbAHPQKe7uOgQgeWCeztRarE2jRMPVCPwSKSz2Sb55KfIYjtI7QYtS9IYp6CJRbZ4q9Gn++I6TmDnT0BBYlJJ0dx8G9bpg7+5kiO80cZ9n001/MCTQPIGPR9LUrLO1KIV89nGh9xnsj+f3PZd1oceSmeVlplUfF3nLAjwP7r5jGFwz4t5EWC8Gc/fe5AIViXttKNnt3DyBbzmFrPoEoQeQXVQXYbnkbVOMiVTFwXpjIq/HwM9A5TIxI12XT1DrV9b/BTJmJI2Tdi03oL+uCppv+gKtPXm5O2bjHbF6FaZpn0UyxzandsgygxwHg5sXAIj6nfEEZO7+40H1KnfvUZFLp20WYAPpH9iYzFCaRIXWexvdM0DTRrLqxK25FrhpZpG7U7bux2QOV9KENS/z4JNti9C6oSwu0u2ALAubl0+J4e6sOqA5qDJl4uTrgbvvGOB70Xq9emhkFoEvaMHXQuAzGqrtHcm64uT20+fNVLlo5nDbTHTESlcNe56Im7YskVXT87D52kvhdsbaG0QdbmpDQJDyFyJ6JiisnUAucV9hanrnF6ClbpokcsfVp8XQ14V58nbHoytUrHQsNtrdfVhY6p4XflhmopRlDPsHNiZy3BsWT3/9BNDrATvWglxHVdIk6z0v7XHRAPW1VCumNa6bMlw28dBJPQFbfFm1TdUxS7L0X7djQixrgrv3KHh7GGaOnDa2YwXcM0DbS/YpEqsH1dhKrZC5O2PrwpwdsGIdsZGuQa3Y6NiASM4S9uVfOc3zTIWomjC4fg7U685OMlgHkS5Khfei8cG3iNb46GdlyiyC3kk3LX2BRSa1cA3or58AADkq25uxdHFqYb1XSM3u7BKoo2W6YFph0Re48TKnLk56zZcuGhNRUw/Y80UY62gIjGRsfIK7rfFUfN+17GgZdHJbJ0ro6yL281pXaUnJ4vPIwuDGeVCnM9/2DKWhagoHTOtHMaTS3E5WRdqJrpvbYckUTltch+OYs+2JMfLhzMh0fXCUoXzmHS3s7j4c9JGo0cw8Hk9mjyzQwbp090xBQ2o1OlkNuSh8MdfBqs/Z9qlGS3we+6IsnKESeOyBx8V86P2DLqBKMNo2Xnz9ZQxeewnkmK7DrMwl8ET094joJhF9lYh+l4jWiOhDRPQKEb1BRF8koq5ctif/f1POf6SMHUilLm6GGlFZTvpFME+7p/rj6/sG22Tc/cfF+AJfjDdQKSLysHlpIMXdArod9A+6sohNfJBbA633BVFY4InoQQD/E4CPMPN/AcAG8CkAvwbg15n5UQDfA/Bp+ZNPA/geM/8YgF+Xy81PjV1MdWXui3tZVn1VN6XPcPcfDyI38mAGRiUjBpXJ8FPPL9x5Td0OqNcDWRZ4eztMX9D0ztYF3TvzHiUHwA4icgDsBPAWgCMAnpfzPw/gZ+X3p+T/kPOPUtU9JcaKT2Vuaz5c0WKPc8Y2zyw0blnCAlT+/YIVr9j30V8/sfJx8xOJ3fYejRgC7u7DgRWf9FDULXx9/ualATYvDUREjVfOdbZ0632B90thgWfmPwXwTwF8E0LY3wZwGcD3mVnl/LwN4EH5/UEA35K/HcvlPxhfLxE9Q0SvEtGrI2wXbZ4hI6Vd7Iu06udtc1zcgSAbZR766yfEg8HzV3rka/+gC5YpH9x9x4S4j8fRjKWeB2xvw91zZMInHywvv/NwGKwLEMdZDXh68fWXRfqCWbahMe4AzOeieT+EVf4hAH8VwH0ATiYsqu7GpDMycacy83PM/BFm/kgHvaLN01ZoTvQsSrdoFiH0RdtMhBdffxnYsQZ0uhjcuih8vFa+nCfu/uMimoMZ8EWREWXFFnH1NJnNSwNxvj1fxK0PR0HUy8SoYaJIUrHAWvc8UZ5PFfzw/OABrCo8MbOIquk4oE5DO1oXrEfzHKW/CeBPmPkvAICI/gDATwK4n4gcaaU/BODbcvnbAB4GcFu6dN4H4LtzbD8kXuVpYn47qj5VCftcPJQyfaXh9yqO/4wRrxNFQqTv1n38SQxeC33C1OmAk4p5T2Fw9YyIvgkyU47BviesUTJGRdzdFZwLeR24e46A1tbAW/fCPDM6NqJZPz0PfHccnEPrvp3grZQ3fNO5GjDPXfdNAE8Q0U7pSz8K4BaACwA+IZd5GsCX5PcX5P+Q889znYPwV5DS/PLJK6/Gsp9jxKsKj9y8fAqDq2dSY7VTI0BUHhv18TyZzbJ9At8/6Kb2M/QPumHCNj2Pe+zD0k1z8rGPg4cj+D94BzweRyz9QBLYFxE4AGBb4TplbiH/+283r1zfEq6LwhY8M79CRM8DuAJgDOAqgOcAbAL4PSL63+S0z8mffA7AvyaiNyEs90/N0/D8DTZWfFYqseajGxB/yzofeXPX+D7cx58ELGFxk+NM5J/pH9gAmBNzobh7jsgEZjHLUwocus7Uot/9AxuNKhMoyumNJkeXSvT9dB9/cuqbECdZ6zoqOoY0t5nnRwtvT3vTqqv1vqSH/lyOLGb+VQC/Gpv8dQA/nrDsFoBPzrO9GY2Z3fFiyEzlIi82Iv6WIfRZRF4LrWNmEEiIeGJmUinuvqgu5O46JKb3eoC3HVqqEiIl7tbUkEDRcVhd8qx5SXowbV45LTpH/fQHorvvGDAeI4yvQL4QZq1wC1sWXrx5IfDPU68LHg6D+WLVpmhLFlbLpG3hq3OVLMziKct1k9DeRPFWQm/bGNw4nyjIiRklfQbubYlEWLEScoF7YTxO7WTtr59I3de6dMzyeJw8Q5bRg+enu0b0tMxpy8RdN8GGOfzINyweicgc3tqKtjG1Lq+x3uOslsADRuRzUqlffnJj84t9hraSbYFskWBMCavyx7v7jgWW4+DaWeEm0HKhRPzLwIRYseeDt7YSY703r5wW7gZPxIXr8HgcVJYCwgLh/fUTpQ6m6h90ReWq9RNiG/uOBeGI7p4jQbiijrvrUBjqyOGgJdWu/voJDK6dlYWzeVLckwR92jztYUmdjijCspYQUaefa3NfJ9L8ZGNxsrhpjC++EJW7bJI3Wux3sbaSSjVr26BeF4Pr54SgqVS0FoHuuw98954QC9sGra1h89IA7q5Dwl0DzI60seTD49bFQDjJtoSbQyXOmrVvWnUpWBQUDS8L4VNXD9NY4RPZrsHNC6GwA5PuKNsG7dwB3t4GhqNoJ6giR1RSQOCDpyAZHLodcV48L+rDzyjwbbPe8yQbW02BB4zIF2QpIi82nP83WluVwFvvfY9wQ5AFjIahcEMIMXt+OAiq1wM5DvjuvSCDIYBJCzV2zdFaT1j+o2HYN0BWcjhgrJ3Bflok1mFVk6Pe3Xcs7CTWxV0Tw4n8+lofBhFh8NpLIj5dHZssnaBZsKxA4K33vgf8zrtAxwHfuZvse69rSb6K3ipWO5tk1geWeaUrxFJvlrzum7hY+b4Izbu3Bb57V46y9IMi3DwcheF8ks1LAzEQSi/oPdE2jnx4axt8927gq2fPTwwHDF0+nOpaqqoASTDYKEHco4OTkt1Raj4Ph2F4aIIbJm2fJ8IiEyAi8DvvimNnxL0QDR0OZlgmuUsDlt8AyAbMXlaLrmFmkN4xqoRcT1ylhGs0BnW6orPP28oXEZLFko1H9KgoHMuPWvFVYtkAPBHUk+CKmVUV6+SP/SSA5H3N4xmYzA7pA50OYNuTud+N3z0X7bPg82AukLlY+qjArBa9ZiEnWo1JkR3M4Dt3wnQEacsmfYL2pXQ6qvVMgTqdyssHBp3ISaTtRyzaJb7Ps6zyaUSs+tEIfO8eeDgsLO6rbr0DbRX4XNZWfU5GE1m6yItGZBd69ZNpQpRFvKduZ4rYx5dLYoGFvwfXzgZvDBx/85j2cEpxxZTFhPuqKZZ7zdpmXDSGuVnIoKhsDRF/p7lulFBobps85M5wHc9dHh+Ql5TbfEEPTZUFEhiL4+HpIjrl4RejkkCNAkU8oovVwPCoAe204PNSs6duE1lovPwsslj0BduatYMwuq1YW9IsYyAIjSxSASkvm5cGwoJ3HNCOHbN/UIW4+5z8SaLu92kN29degc974dXw5BjmJIvIz/lQKizySVgUxMBX7X/XGVw7C/ZleGhaO8sQ96xCnkSO+3MphkZN9aO9Al+Emp6kJlErSx6o1JoPNpHVmk+y5NV0izC4dRGDG+dBvRLqIOQkUuwkQ19B9gdbATGf2JgR96K0b6BTnCIJyMwgqFKohV8+TpZzO2e7Z/rp4z53NWpTjoBdFP0DGyIPPhAU2mA1sndKqOdMzShTZHOK5yoI/GoPdDLUhlpZ8oqc0TaFNlHUaPI5yI0zLf96qfgcVlEq43yVXQIy10/aL+55ab/AF7nZan7SmkRtRT6rf36Oztgiv+HRCO6+YyK6ZVbu9DnZvHwK5NjBfs4sVI4p+1VCf0a4kfz3nxH3ZNov8EVpwMlrCrUUeSD7OV6gyMPzhEUtE55VzeblU6BeN3RLxd1HGon7U6awA8257xrSTiPw02jISWwCtet8VeQdDZtzH5JFcVoKALk8WakVocqGPZEegWwrX5x/2cJe8H5b+HXVIF1YDYGfpyO5QSezCdRS5IF85zmn0Ge15ANxtShMBrYABtfOArYF2rEDg9deAiW8OUy6b1bQam8gqyHw82IuwFKptchXJPSpIi9FXRd3vbO1bFTHbf/AhigoIv8fXD8HWFZYmnAaS+xInVyFsd6nsTqpCuJDxHP/3jfhkyVSm/QGSeQ917H0B6mrVVkjY5BtAx1H+N/lSNaqctEElaN4CBBh8+qZsADI6I6sQzslM+QSO1InV2HEfRZGsQxLo7aWPFDsZs5g0UcSeqkqUxaJARlXOFgAAAqmSURBVE47dojc8zLDY/+gK0r2aZ+iBJb61TOibKAfpkvYvHIaYfEPP+LSLF3cy6q/u2ia2GaswkCnOPNY8YCx4iuittb8POc7ZZ9Cd4wF6nZEeb89R0SKguvnIstGKi8RgTodbF4+hf5BN3MnbH/9hHh7sKywQIdCFh4HxEPA/8E7keLZpQp8iSK5UOOgZuJuBjpVSc1OtqFiKigArotmMHJ0NAbG46AIdlDPVYtT15OU8Xic2aLfvHIaGI9DcdcjgjxPbG/XIfD2tgiZTGjn3KkGjLgvhdXxwSvm9cUDxh9fAUuvEjWNec53in8+qC5lWaIyEsnKSqomrMxHw6NRKDIewJ4ninfbHjhjnHz/wAaA4UThE7FOBkZ3xdsExpEataVQskAacc+HUamitODk15Ha+uXntUIT9kvViZX/QFVJYmYMrp0VQj4eT2ZgZOkn9zwp3hlwnMkc+LF6qzwaB+6Zua33pvraFU1uu4YR+HloyUVQN2or8kAlLpt41AoRCReN54WFu5MyVmZ8qxApCcTLekTc08rwzUtF98XCrosW3der56IBynHTBOsy7poqaFUYpY5WBBzQQieVsFqWqAOblIfGssTyoKBjVJHU6do/6ILHYwyunonmtolFyoRtK6Fj1Yh7rVhdZSozeqhlF0VdWBVLPiqoftS6jvwuLEno7j0q4tkBUf0p4YHAY9lxu/domIo4yULXXDWFo+oqdMkYcS/O6oVJ6pRlxQfrW93nZdXU1povMYwyMQ9MWu54EtWfaG1NdMQCgGUHKQ76BzbEdC8U3mAQU/yBolEoJUGFwmjEfRITJpmVsh9uTe9YqjG1tearsuSDZRKqQMlOUPZ88NZWEPpIjh0MaNq8fEpko2S//IyPkfa04Hpvwz6ksNoCb2gUtRX5eSgUoRJ2hnIwKtUHj0aBH75/0I3Udc3y8MhtvVcsjAs53y0Wd2DVXTSKsl01wXrN87MqaueyWZKrJrK8bYuskL2esPK1gU2R6BkgMWIml8A3XdwbLOzGRVMXGnwR1Z3aWfNVu2oSf5ewzbJCHafRdHFfIYzAA+X74iPrNiJfFbUTgirPdVbRlpb95uVTIpTStkVn7Iy31MzWe9PFfcX6yYzAL4IVuqAWTe0qRRU910WteIUlomoG188FKYH7B11QpyPSHnQcGXkzxy3fBnFfMWaebSL6l0T0HSL6qjbtA0R0lojekH/fL6cTEf1zInqTiK4T0br2m6fl8m8Q0dPV7M4cVN0XsYIX1yJphcjHV5P1mrRIpiKwRQKy4UgkF7t7L8wtY9mAbYN6XVDHkZZ9DrE34t5Ispzh/xtAPOHFswDOMfOjAM7J/wHgJIBH5ecZAJ8FxAMBwK8C+CiAHwfwq+qhsFKs2OvhoqmdNZ+XOdpOnQ4G186KzJGqY3U8Bra3ge1tkZb45gUMrp8D7dgBWluLdNRObYcR98YyM1UBM79ERI/EJj8F4Kfk988D+CMAvyyn/w4L0+OPieh+InpALnuWmb8LAER0FuKh8btz70GZlJnCYOp2THqDKqlFmoOSznGkCpTMPqlfo0QkCnRfPgUAIu2vJmgiayVEWuBdh4S7ZrgtZlrWbN9+k8V9hYVdUTQXzY8y81sAwMxvEdGPyOkPAviWttxtOS1t+gRE9AyE9Q8A23/Iz381abnKmO96+2EAf1lOQyqnSW0F8rY3IZXLAln4sbUfUN/eTF7g7ak/b9K10KS2AtW09z/PumDZycaSzCaeMn1yIvNzAJ4DACJ6lZk/Ul7zqqVJ7W1SW4FmtbdJbQWa1d4mtRVYfnuLvkP+uXS9QP79jpx+G8DD2nIPAfj2lOkGg8FgqIiiAv8CABUJ8zSAL2nTf0FG0zwB4G3pyjkN4DgRvV92rh6X0wwGg8FQETNdNET0uxCdpD9MRLchomH+CYDfJ6JPA/gmgE/KxQcAXAhH4F0AvwgAzPxdIvpfAVySy/0vqsN1Bs9l35Va0KT2NqmtQLPa26S2As1qb5PaCiy5vbXORWMwGAyG4phYPYPBYGgpRuANBoOhpdRW4Ilog4hel2kPnp39i8rb8zARXSCirxHRTSL6O3J67rQNC2yzTURXiejL8v8PEdErsq1fJKKunN6T/78p5z+yhLbeT0TPE9Fr8hj/RM2P7d+T18FXieh3iWitLse3aelFUtr7f8hr4ToR/Tsiul+b9xnZ3teJ6IQ2vXLNSGqrNu8fEBET0Q/L/5d+bCMV2+vyAWAD+I8APgygC+AagF1LbtMDANbl9/cA+A8AdgH43wE8K6c/C+DX5HcXwIsQYwCeAPDKEtr8SwD+HwBflv//PoBPye+/BeC/ld//OwC/Jb9/CsAXl9DWzwP4b+T3LoD763psIQbp/QmAHdpx/dt1Ob4AngSwDuCr2rRcxxLABwB8Xf59v/z+/gW29zgAR37/Na29u6Qe9AB8SOqEvSjNSGqrnP4wRGTg/w/gh2tzbBdxQxQ4iD8B4LT2/2cAfGbZ7Yq18UsAjgF4HcADctoDAF6X3/8FgJ/Tlg+WW1D7HoLIE3QEwJflRfaX2k0THGN5Yf6E/O7I5WiBbX2vFEyKTa/rsVUjsz8gj9eXAZyo0/EF8EhMMHMdSwA/B+BfaNMjy1Xd3ti8vwXgC/J7RAvUsV2kZiS1FcDzAPYB+AZCgV/6sa2riyZzaoNlIF+x9wN4BbG0DQBmpW1YFL8B4B8CUAk5Pgjg+8w8TmhP0FY5/225/KL4MIC/APCvpEvpt4noPtT02DLznwL4pxAhwm9BHK/LqO/xBfIfy2Vfvzr/NYQlDNSwvUT0MwD+lJmvxWYtva11FfjMqQ0WDRH9EIB/C+DvMvMPpi2aMG0h+0BEPw3gO8x8OWN7ln28HYjX3s8y834AdxBmKE1iqe2V/uunIFwEfxXAfRCZVNPatOzjO42504tUCRH9CoAxgC+oSQmLLa29RLQTwK8A+MdJsxOmLbStdRX4WqY2IKIOhLh/gZn/QE7Om7ZhEXwMwM8Q0TcA/B6Em+Y3ANxPRGpwm96eoK1y/vsAZBmIVha3Adxm5lfk/89DCH4djy0A/E0Af8LMf8HMIwB/AOAnUd/jCzQwvYjsfPxpAD/P0pcxpV3Lau9fg3jQX5P320MArhDRX6lDW+sq8JcAPCqjEroQHVMvLLNBREQAPgfga8z8z7RZedM2VA4zf4aZH2LmRyCO3Xlm/nkAFwB8IqWtah8+IZdfmLXGzH8G4FtE9JicdBTALdTw2Eq+CeAJItoprwvV3loe34Q21D69CBFtQKQg/xlmvqvNegHAp2Rk0ocgak/8v1iSZjDzDWb+EWZ+RN5vtyGCMf4MdTi2VXWalNCR4UJEqvxHAL9Sg/b8DYjXqOsA/j/5cSF8qecAvCH/fkAuTwD+L9n+GwA+sqR2/xTCKJoPQ9wMbwL4NwB6cvqa/P9NOf/DS2jnfwngVXl8/z1EdEFtjy2A/xnAawC+CuBfQ0R11OL4QtRZeAvACEJwPl3kWEL4vt+Un19ccHvfhPBTq3vtt7Tlf0W293UAJ7XplWtGUltj87+BsJN16cfWpCowGAyGllJXF43BYDAY5sQIvMFgMLQUI/AGg8HQUozAGwwGQ0sxAm8wGAwtxQi8wWAwtBQj8AaDwdBS/hMF0G/BSb9nEQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "image = np.zeros((500 * 2, 750 * 2), dtype=np.uint8)\n", "t1 = %timeit -c -o -r 3 create_fractal(-2.0, 1.0, -1.0, 1.0, image, 20) # without numba\n", "image = np.zeros((500 * 2, 750 * 2), dtype=np.uint8)\n", "create_fractal, mandel = jit(create_fractal), jit(mandel)\n", "t2 = %timeit -c -o create_fractal(-2.0, 1.0, -1.0, 1.0, image, 20) # with numba\n", "print(\"Speed-up factor: {:.1f}x\".format(t1.average / t2.average))\n", "imshow(image)\n", "show()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Example 4: Diffusion controlled reaction" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "7.37 s ± 109 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n" ] } ], "source": [ "import numpy as np\n", "\n", "RHO_0 = 0.25 # initial particle density\n", "MCS = 10000 # Monte Carlo steps\n", "RUNS = 10\n", "DIM_0, DIM_1 = 81, 81\n", "\n", "def init_lattice(rho=RHO_0, dim0=DIM_0, dim1=DIM_1):\n", " \"\"\"Place particles on lattice with density rho.\"\"\"\n", " lat = np.zeros(shape=(dim0, dim1), dtype=np.int32)\n", "\n", " for i in range(lat.shape[0]):\n", " for j in range(lat.shape[1]):\n", " if np.random.random() <= rho:\n", " lat[i, j] = 1\n", "\n", " return lat\n", "\n", "def execute_mcs(lat):\n", " \"\"\"Move each particle on the lattice once in a random direction.\"\"\"\n", " lat_t0 = np.copy(lat) # lattice at time t0\n", "\n", " for i in range(lat_t0.shape[0]):\n", " for j in range(lat_t0.shape[1]):\n", " if lat_t0[i, j] == 1: # we found a particle\n", " r = np.random.randint(4)\n", " i1, j1 = i, j\n", "\n", " # move particle up/down/left right\n", " if r == 0: i1 += 1\n", " if r == 1: i1 -= 1\n", " if r == 2: j1 += 1\n", " if r == 3: j1 -= 1\n", "\n", " # boundary conditions\n", " if i1 < 0: i1 = 0\n", " if i1 > lat.shape[0] - 1: i1 = lat.shape[0] - 1\n", " if j1 < 0: j1 = 0\n", " if j1 > lat.shape[1] - 1: j1 = lat.shape[1] - 1\n", "\n", " # check trap\n", " if i1 == lat.shape[0] // 2 and j1 == lat.shape[1] // 2:\n", " # we hit the center (trap), remove the particle\n", " lat[i, j] = 0\n", " elif not lat[i1, j1]:\n", " # new position is empty, move particle\n", " lat[i, j] = 0\n", " lat[i1, j1] = 1\n", "\n", "def calc_rho(lat):\n", " return np.sum(lat) / lat.size\n", "\n", "def execute_simulation(runs=RUNS, mcs=MCS, rho_0=RHO_0):\n", " rho_t = np.zeros(mcs, dtype=np.float)\n", "\n", " for r in range(runs):\n", " lat = init_lattice()\n", " for m in range(mcs):\n", " rho_t[m] += calc_rho(lat)\n", " execute_mcs(lat)\n", "\n", " for i in range(rho_t.shape[0]):\n", " rho_t[i] /= (runs * rho_0)\n", "\n", " return rho_t\n", "\n", "t1 = %timeit -c -o execute_simulation(runs=1, mcs=1000)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "59.4 ms ± 1.43 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n", "Speed-up factor: 124.1x\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAD8CAYAAAC2PJlnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd8VGW+x/HPL50kGAgkEEqAUAVCFxCWKiKoKCgoWFddu+tadl253rvu3rWsrmvvbS2rqOtiwwKKWCiiFCF0QguBQCBAKAECyXP/yHFvNlKSIcnJTL7v12temTlzZub3zAn5cs7zPOeYcw4REZFAhPldgIiIBC+FiIiIBEwhIiIiAVOIiIhIwBQiIiISMIWIiIgETCEiIiIBU4iIiEjAFCIiIhKwCL8LqGoNGzZ0LVu29LsMEZGgMn/+/O3OuaTjrRfyIdKyZUvmzZvndxkiIkHFzDaUZz0dzhIRkYApREREJGAKERERCZhCREREAqYQERGRgClEREQkYAoREREJmELkKN6cm8X7Czf5XYaISI0W8pMNAzV5QTaFRcWM7t7U71JERGos7YkcxbCOjVicnc/GHQV+lyIiUmMpRI5ieMdGAExdusXnSkREai6FyFGkJcXTs0V9Xv9uA0XFzu9yRERqJIXIMVzeryUb8gqY9H2W36WIiNRICpFjOLNzY/q3acD/frSMb1Zt87scEZEaRyFyDBHhYTw5oQcp9WL4nw+WUHi42O+SRERqFIXIcdSPi+LuUR11WEtE5AgUIuUwpH0yfVol8sBnK5izJs/vckREagyFSDmYGY+N707dmAjunLyYYo3WEhEBFCLl1jghhjtHdmBDXgEzVub6XY6ISI2gEKmAM9NTaJ0Ux90fLqWg8LDf5YiI+E4hUgHREeHcf14Xsnfu58HPVvpdjoiI7xQiFdS7VSIX9mrOK7PX82lGjt/liIj4SiESgD+e04kuzRL4/b8Wk7v7gN/liIj4RiESgDpR4Tx6YTcOHC7mjx8t9bscERHfKEQClJYUz29Oa8snGVv4fNlWv8sREfGFQuQEXDMwjZYNYnnyy9V+lyIi4guFyAmIDA/jiv6tWJSdz+w12/0uR0Sk2ilETtDYns1ISYhh4uQMzR0RkVpHIXKC4qIj+MPZJSdovOPdxX6XIyJSrRQilWBkegrXDExjyuIcFmbt9LscEZFqoxCpJDef1paG8VE88NkKnNMJGkWkdlCIVJL46Ah+PbQt363dwYeLNvtdjohItVCIVKIJvVNpWq8OD362kn0H1ckuIqFPIVKJoiLCeGx8Nzbn7+d/P1rmdzkiIlVOIVLJerVM5NqBrXl73ka+XrXN73JERKqUQqQK3Hp6W1onxXH7Oz+SvbPA73JERKqMQqQKREeE8/iE7uwvLOL6fyxQ/4iIhCyFSBXp1CSBv5zfhYxN+fzq1XkU6brsIhKCFCJVaFTXJtwzujNz1uZx81sLFSQiEnIi/C4g1F3cJ5XcPQd5fPpq4qLCeXBsV79LEhGpNAqRKmZm3HZ6Ow4cKuL5b9bSJjmeawa29rssEZFKocNZ1eT3IzowoG1D7vtkBTNW5vpdjohIpVCIVJPwMOOx8d1p0SCW299ZxMYdGvorIsFPIVKNEuOieOnyUzhcVMwtb/+ojnYRCXoKkWrWJjmeP57TifkbdvLYF6v8LkdE5IQoRHwwpntTRndrwuNfZjJ3bZ7f5YiIBEwh4gMz44/ndCIuKpz/nbKMQ0XFfpckIhIQhYhP6sVG8dC4rizdvJuHpq30uxwRkYAoRHw0Mj2FC3s157mv1/LanPV+lyMiUmEKEZ/dM6YzfdMS+fOUZSzZlO93OSIiFaIQ8VlkeBhPXdSDBnHRXP3aPHLy9/tdkohIuSlEaoAG8dG8eHkvdhUc4spX5lF4WB3tIhIcFCI1ROemCTxyYVeW5+zmzsmLcU4TEUWk5lOI1CAjOqdw45DWTF6wicenZ/pdjojIceksvjXMb4e3Z0NeAY98sYpi57hlWFvMzO+yRESOSCFSw5gZf7ugK9v3HuSx6auJjQrn2kE6dbyI1Ew6nFUDRUeEM+nqvpzRqRH3f7qCp7/KVB+JiNRI2hOpocyMB8/vSkHhAh78bCUN4qK48JRUv8sSEfkP2hOpwRJiI3ntyt60b1SXP320jMzcvX6XJCLyHxQiNZyZ8dylPYmJDOe6f8xnf2GR3yWJiPybQiQItGwYxyMXdiMzdy8PTl3hdzkiIv+mEAkSg9olcUnfVP4+az2PfbFaHe0iUiOoYz2I/OHsTuzYV8gjX6xiUfYunpjQnbhobUIR8Y/2RIJIVEQYT0zowe9HdOCrlbmMfmoWG3cU+F2WiNRiCpEgEx5mXD+4Nc9e0pOc/ANc8coPbMjb53dZIlJLKUSC1PBOjXnu0p6s3baX295ZpDP/iogvFCJBrH+bhtw9qhPzN+zkyld+IG/vQb9LEpFaRiES5C7v15J7Rndmzto8Bj/0FbMzt/tdkojUIgqREHBJ3xZ8cGN/6sVGcufkDPYePOx3SSJSSyhEQkTnpgncP6YL2TsLuOyluRw4pJntIlL1FCIh5BdtG/KX87uwIGsXv560kMNF6mwXkaqlEAkxF/RqzsSRHfh82Vb+8OFSios1s11Eqo5CJARdO6g15/doxptzs+j/wJd8kpHjd0kiEqIUIiHqoXFd+PPozkRHhHHDGwv41as/sKug0O+yRCTEKERClJlxad8WTLt1EFf0b8mXK3IZ9eRMVm/d43dpIhJCFCIhLioijLtHdeLNq/uy58Bhznt6NrM0l0REKolCpJbom9aAydf3o2HdaC5+cS6PfbGaQxq9JSInKChDxMzSzOwlM3vX71qCSVpSPJ/cPIBhJzfikS9WMeqJmSzM2ul3WSISxMoVImb2GzNbYmZLzeyWQD/MzF42s1wzW3KE50aY2UozyzSzO4/1Ps65tc65qwKtozarExXO85f25K9ju7B19wHOe2Y2V/z9e51SXkQCctwQMbPOwNVAb6ArcLaZtS2zTrKZ1S2zrM0R3u4VYMQRPiMceAoYCXQEJphZRzNLN7MpZW7J5WybHEVYmDGuV3Om3z6YawamMWdtHqf97Wse+XyVruEuIhVSnj2Rk4HvnHMFzrnDwNfAmDLrDAI+MLMYADO7Gni87Bs5574BdhzhM3oDmd4eRiHwFnCucy7DOXd2mVtu+Zsnx5IYF8XEkSfz2W8GMrRDMo9NX8345+ewc5+GAotI+ZQnRJYAA82sgZnFAmcCzUuv4Jz7J/AZ8JaZXQxcCVxQgTqaAhtLPc72lh2RV8uzQHczm3iUdUaZ2fP5+fkVKKN2atkwjmcu6cHvzmjPoux8LnpxruaUiEi5HDdEnHPLgQeAzykJikXAz04T65x7EDgAPAOc45zbW4E67EgffYya8pxz1znnWjvn7j/KOh85565JSEioQBm1l5lx45A2PHdpT1Zv3cOEF+bqbMAiclzl6lh3zr3knOvhnBtIyeGo1WXXMbMBQGfgPeDuCtaRzX/u3TQDNlfwPaQSnNGpMfedl87ynN2c//RssvLU4S4iR1fe0VnJ3s9U4DxgUpnnuwMvAOcCVwCJZnZPBer4AWhrZq3MLAoYD3xYgddLJbqgV3OeuqgH67bvY8jfvuKvU1eow11Ejqi880T+ZWbLgI+AG51zZScXxALjnHNrnHPFwOXAhrJvYmaTgDlAezPLNrOrALwO+5uAqcBy4B3n3NKAWiSV4qwuKXzym18wsG1Dnpqxhgufn8P8DTtwTmcFFpH/Z6H+R6FXr15u3rx5fpcR1N5fuIk7/rWYwsPF9Eitx0PjupKWFO93WSJShcxsvnOu13HXU4hIeeQXHOL9HzfxwGcrKDxczDndmjChdyq9WtTH7EjjIkQkmJU3RILytCdS/RJiI7m8X0um3z6I0d2b8tGizYx7dg63vv0jBYUaxSVSWylEpEJSEurw0LiuzL7zNK4f3Jr3f9zM+c/MYe22iozoFpFQoRCRgCTVjeb3Izpw35h01uTuZcRj3+rMwCK1kEJETshFfVKZdutA+rVuwCNfrKLPfdN5eNpKTVQUqSUUInLCWjaM45UrevPiZb1o36guj3+ZyflPz9aZgUVqAY3Okko3ZfFmbn37R4odnNu1CcM7NWZIhySiI8L9Lk1Eyqm8o7MiqqMYqV3O7tKE9KYJPPrFaj5buoXJCzfRMD6a205vx/hTmhMWpiHBIqFCeyJSpQ4cKmLGilye+DKTZTm76ZFajz+d05n0ZjoxpkhNpnkiUiPERIYzMj2FKb/+BX84uyNLNu9m1JMzuXnSQjbt2u93eSJygrQnItUqb+9Bnv9mLS/OXEexc4zp1pT7zksnJlL9JSI1ifZEpEZqEB/NxDNP5vNbBzK2RzMmL9zEyMe+ZcWW3X6XJiIBUIiIL9KS4nlwbBceGteVzbv2c9bjM7npzQVk79SwYJFgohAR35gZY3s249vfD2H8Kc35JCOHUU/MZGFW2SsNiEhNpRAR3yXXjeHeMem8f2N/IsLDGPfsHO54dxHfrc2j8LBOoyJSk6ljXWqU7XsP8l+TM5i2bCsADeOjuHpAGr8akEa45peIVBtdT8SjEAlOOfn7WbRxF69/t4FZmXk0T6zDMxf3pHNTzS8RqQ4anSVBLSWhDiM6p/CPq/rw0Liu7C8sYvRTs3hl1jpdolekBlGISI32U+f7P6/rR4/U+vzxo2Xc9f4S9ukswSI1gkJEgkKrhnFMuqYvl/RN5c25WQx/5Bse+GwF2/Yc9Ls0kVpNISJBIzzMuGd0Oq9d2Zu46HCe+WoNIx/7lvcXbtLFsER8ohCRoDOwXRLTbh3Ev64/lZjIMG55+0cmPP+dLoQl4gOFiAStni0SmX77IO4Y0Z55G3bS594vuOPdRezYV+h3aSK1hq4nIkEtOiKcGwa3oUdqfd76Pot35mXz/sLNDOuYzNiezRjSPhkzzS8RqSoKEQkJfdMa0DetARec0pw35mbxzaptfJKxhd6tEnn5l6cQH61fdZGqoH9ZElL6tW5Iv9YNOXCoiCe+XM1TM9Zw+sNfc+3ANM7o3JiUhDp+lygSUjRjXULatKVb+MunK1i7fR8AA9o25PcjOmjmu8hx6LQnHoWIFBc7FmTtZPaaPJ77eg37CotITYxlZOfGXDuoNYlxUX6XKFLjKEQ8ChEpLb/gEG98v4HZmXnMzNxO/dhILu7TgusGt1a/iUgpChGPQkSOJiM7n3s+XsbcdTtIS4rjd8PbM6RDsi7VK4JC5N8UInI8UxZv5v5PVrBp135io8K5blBrrh2URnSEwkRqL4WIRyEi5XGoqJhvV2/jpZnrmJWZR9N6dfjTOZ0Y1rGR36WJ+EIh4lGISEU455i+PJeJ72Wwbc9BBrZL4vweTRnaIZm6MZF+lydSbcobIupJFCnFzBjWsRGntErkjncXMWPlNr5ZtY2IMOO+89IZ2bmxwkSkFO2JiBxDcbFjzto8bp60kLx9hdSNieDx8d0Z3D5Jp1ORkKYrG4pUgrAwo3+bhsy6cyjPXNyD2KhwrnjlBya88J2uZSKCQkSkXGIiwxmZnsL02wfz2+HtWLBhF8Mf+ZrZmdv9Lk3EVwoRkQqIj47gpqFteevavsRFR3DRi3O54Y35rNq6x+/SRHyhEBEJQI/U+nz86wFcdmoLvlq5jVFPzOTZr9eQu/uA36WJVCt1rIucoOydBdz69o/8sH4nAEPaJ3HZqS0Z1C6JsDB1vktw0hBfkWrSrH4s71x7KvM37OQf321g2rKtzFi5jeS60fxqQCsu7duSOlGa/S6hSXsiIpVsf2ER05Zt4eWZ61iUnU+YwaB2SQxsl8Tobk2pr7MGSxDQjHWPQkT84lzJHJMPf9zMh4s2U1BYRL3YSG4e2paB7RrSJrmu3yWKHJVCxKMQkZrAOcfsNXnc+/FyluXsBuDSvi34rzNP1qEuqZHUJyJSg5iVTFr8+OZfsGbbXp75ai2vf7eBVVv3cO+YdNokx/tdokhANMRXpBqZGW2S6/K3C7pyz+jOzNuwk9Mf+ZqJkzPIzN1LqB8ZkNCjPRERn1zStwV90xrw2PTV/HPeRiZ9n0XLBrHcMaIDZ6an+F2eSLmoT0SkBti0az9Tl2xh0vdZrM7dy4C2Dbn19Hb0SK3vd2lSS6lj3aMQkWCy7+BhnvlqDS/NXMf+Q0U0T6zDn8/tzOD2yX6XJrWMQsSjEJFglF9wiHfmbeS5b9ayfe9Bzu6SwrhezRnULsnv0qSWUIh4FCISzHYfOMTD01bx1g9ZHDhUTIfGdbmkbwvG9Wqma8BLlVKIeBQiEgr2Fxbx99nreOeHjazPK6B1Uhw3n9aWwe2TSaijKy1K5VOIeBQiEkqcc0xZnMNfPl3Bpl37CQ8zTuuQzLWD0ujZItHv8iSEaLKhSAgyM0Z1bcKIzo2Zlbmdjxfn8N7CTUxbtpWhHZK566yTaZ2kiYtSfbQnIhLkdhUU8uK363jqq0ycg5uGtOGmoW2IiVSfiQROh7M8ChGpLTbk7eN37y7m+3U7aHxSDCM6N+aWYW2pF6uzBkvFlTdEdNoTkRDRokEcb13dlxcu60WHlLq8Mns9Zz8xk4zsfL9LkxCmEBEJIWFhxukdG/HKFb155YpT2HPgMGOensXDn6+ioPCw3+VJCFKIiISowe2T+fL2QfRv05DHp6/m1Pu/5PNlW/0uS0KMQkQkhDWIj+aVK07h7788hZSEGK59fR5/nrKM3QcO+V2ahAgN8RUJcWbGkA7J9GxZn7veW8JLM9fxaUYOF/VJ5fyezUhJqON3iRLEtCciUkucFBPJExO68+JlvTAzHpq2ilPv/5L/fj+D3D0H/C5PgpSG+IrUUss27+bprzKZsjiHqPAwzuvRlEv6tqBz0wS/S5MaQPNEPAoRkWNbsimfp7/K5IvluRQeLubOkR24blBrv8sSnylEPAoRkfLZua+Q85+dzdpt++jdKpG/ju1CiwZxfpclPtFkQxGpkPpxUXxy8wDuGNGeRRt3cdbjM7nj3UVk7yzwuzSpwbQnIiI/k5m7h0c+X83UpVsAOKNTY64f3Fr9JbWIDmd5FCIigcvM3cPTM9YweeEmAO4682RGpjemab06mJnP1UlVUoh4FCIiJy4rr4ArX/2BzNy9ANSPjeSK/q04o1Nj2jWKV6CEIIWIRyEiUjmKih0LsnayPGc3U5duYVZmHgCnd2zE/eel0zA+2ucKpTIpRDwKEZGqsX77Pl6etY7X5mwgKiKMc7o24eahbUltEOt3aVIJFCIehYhI1VqcvYtXZq1nyuIcHI4LT2nONQNaK0yCnELEoxARqR5ZeQXc/+lyPl2yhTqR4dw7pjOjuzUlLEz9JcFI80REpFqlNojlmUt68snNA2jbKJ7b3lnEwL/O4NOMHEL9P6u1mUJERCpVxyYn8d4N/blvTDq79x/i+jcWMPqpWbzzw0aKihUmoUaHs0Skyhw4VMRrc9bz2pwNZO/cT6OTorl1WDvG9mxGRLj+D1uTqU/EoxAR8V9xsePjjBz+8ukKNu3aT0SY0bJhHGN7NuOSvi2Ij9aljWoahYhHISJScxw8XMSnGVtYnrObuet28OPGXURFhHFe96bcNyZdnfA1iELEoxARqbm+X7eDV2ev5+OMHFISYpjQO5UJvVNJqquJi35TiHgUIiI1m3OOD37czMuz1rE4Ox+AtKQ4Tu/YiGsHtiYxLsrnCmsnhYhHISISPJZsymfa0i1MW7aVFVv2kJoYy92jOjK4fTLhOtRVrRQiHoWISHCavnwrv3t3MTv2FdKuUTy/Hd6eYSc3Ur9JNVGIeBQiIsFrf2ERr85Zz9MzMtl94DDN6tfh4j4tGH9Kc+rrMFeVUoh4FCIiwe/AoSI+XLSZF75Zy+rcvdSNieCaAWlM6JOqswdXEYWIRyEiEjqccyzZtJu7P1zCgqxd1IkM5/J+LfndGe3VZ1LJyhsimuEjIkHDzEhvlsDkG/rz/bodPPv1Gp79eg3rtu/l/B7NGNguiZjIcL/LrFUUIiISlHq3SuSUlvV5+PNVPDkjk6lLt9L4pBj+cn46A9omac+kmuhwlogEve17DzJjRS5/+mgZew8e5uSUk7hpSBtOOzlZeyYBCuk+ETNLA+4CEpxzY4+1rkJEpPYoKDzM5AWbeGpGJjn5B0iMi+J3Z7Tn/B7NiIrQCR8rolKvJ2Jmt5rZUjNbYmaTzCwmwKJeNrNcM1tyhOdGmNlKM8s0szuP9T7OubXOuasCqUFEQldsVASX9G3BjN8O5q9ju1A3JoKJkzMY8eg3vDs/m8NFxX6XGHKOGyJm1hS4GejlnOsMhAPjy6yTbGZ1yyxrc4S3ewUYcYTPCAeeAkYCHYEJZtbRzNLNbEqZW3I52yYitVRMZDjjejVn2q0DefTCbmDw238u4uwnZrJyyx6/ywsp5d2/iwDqmFkEEAtsLvP8IOCDn/ZQzOxq4PGyb+Kc+wbYcYT37w1kensYhcBbwLnOuQzn3NllbrnlKdjMRpnZ8/n5+eVsooiEmuiIcEZ3b8r02wZx35h01mzby1mPf8sfPljCd2vz/C4vJBw3RJxzm4CHgCwgB8h3zk0rs84/gc+At8zsYuBK4IIK1NEU2Fjqcba37IjMrIGZPQt0N7OJR6n7I+fcNQkJCRUoQ0RCkZlxUZ9Upt82mFFdm/DG3CzGP/8d//VeBgcPF/ldXlArz+Gs+sC5QCugCRBnZpeUXc859yBwAHgGOMc5t7cCdRxpLN5Re/ydc3nOueucc62dc/dX4HNEpBZLbRDLIxd2Y8H/nM5FfVJ5c24WAx+cwX2fLCcnf7/f5QWl8hzOGgasc85tc84dAiYD/cquZGYDgM7Ae8DdFawjG2he6nEzfn7ITESkUiTUieS+Mek8fXEPWjWM4/lv1nLq/V9y13sZbN6lMKmI8oRIFtDXzGLNzIDTgOWlVzCz7sALlOyxXAEkmtk9FajjB6CtmbUysyhKOu4/rMDrRUQq7Mz0FN665lTev7E/Z3dJ4Y25WfzigS957us1OsxVTuXpE5kLvAssADK81zxfZrVYYJxzbo1zrhi4HNhQ9r3MbBIwB2hvZtlmdpX3GYeBm4CplATUO865pQG3SkSkAro1r8eTF/Xgs1sG0LtVIvd/uoLe907n4c9XaVjwcQTlZMOK0GRDEamI4mLHtGVb+Pus9cxdt4O60RGc0bkxw05Opl+bhpwUE+l3idUipGesV4RCREQCUVzs+GzpFqYu3cLny7ZSUFhEZLgxvFNjBrdLYnjHxiTEhm6gKEQ8ChEROVH7C4uYlbmdjzNymL58K7sPHKZ+bCTDOzZmRHpjhrQPvTnQChGPQkREKtPhomK+W7uDV+esZ/ryrRQ76N0ykVFdUxh6ciOa1qvjd4mVQiHiUYiISFU5cKiIhz9fxbvzs9mxr5DYqHBuH96es9JTaJwQ0CkGawyFiEchIiJVrbjYsWRzPv/z/hIWZZecaqlr83r8ekgbhnVs5HN1gVGIeBQiIlJdnHN8t3YH89aXHO7K21fIFf1acceI9kF3XROFiEchIiJ+2HvwMP/9Xgbv/7iZ+rGRnNejGRNHdiAiPDiua1Kp1xMREZGKiY+O4NHx3Xn+0p6cnHISL81cx7jn5rB++z6/S6tUChERkSo0vFNj3vhVHx4c24VFG3dxxqPfcO/Hy9i4o8Dv0iqFQkREpIqZGRf0as7UWwYyrGMjXvh2HQP/OoNJ32f5XdoJU4iIiFSTto3q8tRFPZh++yBSE2OZODmDa1+fx6KNu/wuLWARfhcgIlLbtE6KZ+otA3n0i9W8Nmc9U5dupXerRC7qncrI9MZERwTPSC6NzhIR8VH+/kM8+sUqvlyRy4a8AurGRHDrsHZc1CfV12HBGuLrUYiISDA4VFTM58u28uK3a1mQtYvUxFhuGtqGc7o28SVMFCIehYiIBJPiYseUjBzu/2Q5OfkHSEuKY/wpzRnTvRlJdaOrrQ6FiEchIiLBqLjY8eGizTz6xSrW5xUQHmb8emgbrhmYRmxU1XdnK0Q8ChERCXZLNuXzt2krmbFyG/VjIzm3W1OuGZhGkyo8Y7BCxKMQEZFQ4JxjxspcXp65npmZ2zGDU9MacEGv5pzbrQlmVqmfpxDxKEREJNSs3rqHN7/P4v2Fm9hZcIi0pDi6NavHjUPb0DopvlI+QyHiUYiISKgqPFzMO/M28umSHGZl5gHQr3UDbhjchn6tGxAWFvjeiULEoxARkdpgzba9vPV9Fm/MzaKgsIi2yfG8flWfgC+OVd4Q0Yx1EZEQ0DopnrvO6sjNp7XlX/OzmZmZR73YyCr/XIWIiEgIqRsTyS/7t+KX/VtVy+fpBIwiIhIwhYiIiARMISIiIgFTiIiISMAUIiIiEjCFiIiIBEwhIiIiAVOIiIhIwEL+tCdmtg3YEODLGwLbK7GcYKA21w61rc21rb1w4m1u4ZxLOt5KIR8iJ8LM5pXn3DGhRG2uHWpbm2tbe6H62qzDWSIiEjCFiIiIBEwhcmzP+12AD9Tm2qG2tbm2tReqqc3qExERkYBpT0RERAKmEDkCMxthZivNLNPM7vS7nhNhZs3NbIaZLTezpWb2G295opl9bmarvZ/1veVmZo97bV9sZj1Kvdfl3vqrzexyv9pUXmYWbmYLzWyK97iVmc316n/bzKK85dHe40zv+Zal3mOit3ylmZ3hT0vKx8zqmdm7ZrbC296nhvp2NrNbvd/rJWY2ycxiQm07m9nLZpZrZktKLau07WpmPc0sw3vN42ZWsWvqOud0K3UDwoE1QBoQBSwCOvpd1wm0JwXo4d2vC6wCOgIPAnd6y+8EHvDunwl8ChjQF5jrLU8E1no/63v36/vdvuO0/TbgTWCK9/gdYLx3/1ngeu/+DcCz3v3xwNve/Y7e9o8GWnm/F+F+t+sY7X0V+JV3PwqoF8rbGWgKrAPqlNq+vwy17QwMBHoAS0otq7TtCnwPnOq95lNgZIXq8/sLqmk378ucWurxRGCi33VVYvs+AE4HVgIp3rIUYKV3/zlgQqn1V3rPTwCeK7X8P9araTegGTAdGApM8f6BbAciym5nYCpwqnc/wlvPym770uvVtBtwkvcH1cosD9mxWM/cAAAC00lEQVTt7IXIRu8PY4S3nc8Ixe0MtCwTIpWyXb3nVpRa/h/rleemw1k/99Mv5k+yvWVBz9t97w7MBRo553IAvJ/J3mpHa3+wfS+PAncAxd7jBsAu59xh73Hp+v/dNu/5fG/9YGpzGrAN+Lt3CO9FM4sjhLezc24T8BCQBeRQst3mE9rb+SeVtV2bevfLLi83hcjPHel4YNAPYTOzeOBfwC3Oud3HWvUIy9wxltc4ZnY2kOucm1968RFWdcd5LmjaTMn/rHsAzzjnugP7KDnMcTRB32avH+BcSg5BNQHigJFHWDWUtvPxVLSNJ9x2hcjPZQPNSz1uBmz2qZZKYWaRlATIG865yd7irWaW4j2fAuR6y4/W/mD6XvoD55jZeuAtSg5pPQrUM7MIb53S9f+7bd7zCcAOgqvN2UC2c26u9/hdSkIllLfzMGCdc26bc+4QMBnoR2hv559U1nbN9u6XXV5uCpGf+wFo643wiKKkA+5Dn2sKmDfS4iVguXPu4VJPfQj8NELjckr6Sn5afpk3yqMvkO/tLk8FhptZfe9/gMO9ZTWOc26ic66Zc64lJdvvS+fcxcAMYKy3Wtk2//RdjPXWd97y8d6onlZAW0o6IWsc59wWYKOZtfcWnQYsI4S3MyWHsfqaWaz3e/5Tm0N2O5dSKdvVe26PmfX1vsPLSr1X+fjdYVQTb5SMcFhFySiNu/yu5wTb8gtKdk8XAz96tzMpORY8HVjt/Uz01jfgKa/tGUCvUu91JZDp3a7wu23lbP9g/n90VholfxwygX8C0d7yGO9xpvd8WqnX3+V9Fyup4KgVH9raDZjnbev3KRmFE9LbGfgTsAJYArxOyQirkNrOwCRK+nwOUbLncFVlblegl/f9rQGepMzgjOPdNGNdREQCpsNZIiISMIWIiIgETCEiIiIBU4iIiEjAFCIiIhIwhYiIiARMISIiIgFTiIiISMD+D2NcYN5dsRX9AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from numba import jit\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "RHO_0 = 0.25 # initial particle density\n", "MCS = 10000 # Monte Carlo steps\n", "RUNS = 10\n", "DIM_0, DIM_1 = 81, 81\n", "\n", "@jit(nopython=True)\n", "def init_lattice(rho=RHO_0, dim0=DIM_0, dim1=DIM_1):\n", " \"\"\"Place particles on lattice with density rho.\"\"\"\n", " lat = np.zeros(shape=(dim0, dim1), dtype=np.int32)\n", "\n", " for i in range(lat.shape[0]):\n", " for j in range(lat.shape[1]):\n", " if np.random.random() <= rho:\n", " lat[i, j] = 1\n", "\n", " return lat\n", "\n", "@jit(nopython=True)\n", "def execute_mcs(lat):\n", " \"\"\"Move each particle on the lattice once in a random direction.\"\"\"\n", " lat_t0 = np.copy(lat) # lattice at time t0\n", "\n", " for i in range(lat_t0.shape[0]):\n", " for j in range(lat_t0.shape[1]):\n", " if lat_t0[i, j] == 1: # we found a particle\n", " r = np.random.randint(4)\n", " i1, j1 = i, j\n", "\n", " # move particle up/down/left right\n", " if r == 0: i1 += 1\n", " if r == 1: i1 -= 1\n", " if r == 2: j1 += 1\n", " if r == 3: j1 -= 1\n", "\n", " # boundary conditions\n", " if i1 < 0: i1 = 0\n", " if i1 > lat.shape[0] - 1: i1 = lat.shape[0] - 1\n", " if j1 < 0: j1 = 0\n", " if j1 > lat.shape[1] - 1: j1 = lat.shape[1] - 1\n", "\n", " # check trap\n", " if i1 == lat.shape[0] // 2 and j1 == lat.shape[1] // 2:\n", " # we hit the center (trap), remove the particle\n", " lat[i, j] = 0\n", " elif not lat[i1, j1]:\n", " # new position is empty, move particle\n", " lat[i, j] = 0\n", " lat[i1, j1] = 1\n", "\n", "@jit(nopython=True)\n", "def calc_rho(lat):\n", " return np.sum(lat) / lat.size\n", "\n", "@jit(nopython=True)\n", "def execute_simulation(runs=RUNS, mcs=MCS, rho_0=RHO_0):\n", " rho_t = np.zeros(mcs, dtype=np.float32)\n", "\n", " for r in range(runs):\n", " lat = init_lattice()\n", " for m in range(mcs):\n", " rho_t[m] += calc_rho(lat)\n", " execute_mcs(lat)\n", "\n", " for i in range(rho_t.shape[0]):\n", " rho_t[i] /= (runs * rho_0)\n", "\n", " return rho_t\n", "\n", "t2= %timeit -c -o execute_simulation(runs=1, mcs=1000)\n", "print(\"Speed-up factor: {:.1f}x\".format(t1.average / t2.average))\n", "rho_t = execute_simulation()\n", "plt.plot(rho_t)\n", "plt.yscale('log')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# References\n", "\n", "1. [Stanley Seibert - Accelerating Python with the Numba JIT Compiler (SciPy 2015)](https://www.youtube.com/watch?v=eYIPEDnp5C4)\n", "* Stanley Seibert - Numba: A JIT Compiler for Scientific Python (Continuum Analytics 2015)\n", "* Travis E. Oliphant - Performance Python: Introduction to Numba (PyData 2015)\n", "* [Numpy documentation](https://docs.scipy.org/doc/)\n", "* [Numba documentation](http://numba.pydata.org)\n", "* [Numba examples](https://github.com/numba/numba/blob/master/examples)" ] } ], "metadata": { "celltoolbar": "Slideshow", "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.4" } }, "nbformat": 4, "nbformat_minor": 2 }