Skip to content

GitLab

Commit 1e02cdc7 authored by Martin Bauer's avatar Martin Bauer
Browse files

Separated modules into subfolders with own setup.py 

This restructuring allows for easier separation of modules into
separate repositories later. Also, now pip install with repo url can be
used.

The setup.py files have also been updated to correctly reference each
other. Module versions are not extracted from git state
parent b91e6990
Hide whitespace changes
Inline Side-by-side
doc/notebooks/demo_benchmark.ipynb 0 → 100644
View file @ 1e02cdc7
%% Cell type:code id: tags:
``` python
from pystencils.session import *
import timeit
%load_ext Cython
```
%% Cell type:markdown id: tags:
# Demo: Benchmark numpy, cython, pystencils
In this benchmark we compare different ways of implementing a simple stencil kernel in Python.
The benchmark kernel computes the average of the four neighbors in 2D and stores in a second array. To prevent out-of-bounds accesses, we skip the cells at the border and compute values only in the range `[1:-1, 1:-1]`
%% Cell type:markdown id: tags:
## Implementations
The first implementation is a pure Python implementation:
%% Cell type:code id: tags:
``` python
def avg_pure_python(input_arr, output_arr):
for x in range(1, input_arr.shape[0] - 1):
for y in range(1, input_arr.shape[1] - 1):
output_arr[x, y] = (input_arr[x + 1, y] + input_arr[x - 1, y] +
input_arr[x, y + 1] + input_arr[x, y - 1]) / 4
```
%% Cell type:markdown id: tags:
Obviously, this will be a rather slow version, since the loops are written directly in Python.
Next, we use *numpy* functions to delegate the looping to numpy. The first version uses the `roll` function to shift the array by one element in each direction. This version has to allocate a new array for each accessed neighbor.
%% Cell type:code id: tags:
``` python
def avg_numpy_roll(input_arr, output_arr):
neighbors = [np.roll(input_arr, axis=a, shift=s) for a in (0, 1) for s in (-1, 1)]
np.divide(sum(neighbors), 4, out=output_arr)
```
%% Cell type:markdown id: tags:
Using views, we can get rid of the additional copies:
%% Cell type:code id: tags:
``` python
def avg_numpy_slice(input_arr, output_arr):
output_arr[1:-1, 1:-1] = input_arr[2:, 1:-1] + input_arr[:-2, 1:-1] + \
input_arr[1:-1, 2:] + input_arr[1:-1, :-2]
```
%% Cell type:markdown id: tags:
To further optimize the kernel we switch to Cython, to get a compiled C version.
%% Cell type:code id: tags:
``` python
%%cython
import cython
@cython.boundscheck(False)
@cython.wraparound(False)
def avg_cython(object[double, ndim=2] input_arr, object[double, ndim=2] output_arr):
cdef int xs, ys, x, y
xs, ys = input_arr.shape
for x in range(1, xs - 1):
for y in range(1, ys - 1):
output_arr[x, y] = (input_arr[x + 1, y] + input_arr[x - 1, y] +
input_arr[x, y + 1] + input_arr[x, y - 1]) / 4
```
%% Cell type:markdown id: tags:
And finally we also create a *pystencils* version of the same stencil code:
%% Cell type:code id: tags:
``` python
src, dst = ps.fields("src, dst: [2D]")
update = ps.Assignment(dst[0,0],
(src[1, 0] + src[-1, 0] + src[0, 1] + src[0, -1]) / 4)
kernel = ps.create_kernel(update).compile()
def avg_pystencils(input_arr, output_arr):
kernel(src=input_arr, dst=output_arr)
```
%% Cell type:code id: tags:
``` python
all_implementations = {
'pure Python': avg_pure_python,
'numpy roll': avg_numpy_roll,
'numpy slice': avg_numpy_slice,
'Cython': None,
'pystencils': avg_pystencils,
}
if 'avg_cython' in globals():
all_implementations['Cython'] = avg_cython
else:
del all_implementations['Cython']
```
%% Cell type:markdown id: tags:
## Benchmark functions
We implement a short function to get in- and output arrays of a given shape and to measure the runtime.
%% Cell type:code id: tags:
``` python
def get_arrays(shape):
in_arr = np.random.rand(*shape)
out_arr = np.empty_like(in_arr)
return in_arr, out_arr
def do_benchmark(func, shape):
in_arr, out_arr = get_arrays(shape)
timer = timeit.Timer('f(a, b)', globals={'f': func, 'a': in_arr, 'b': out_arr})
calls, time_taken = timer.autorange()
return time_taken / calls
```
%% Cell type:markdown id: tags:
## Comparison
%% Cell type:code id: tags:
``` python
def bar_plot(*shape):
names = tuple(all_implementations.keys())
runtimes = tuple(do_benchmark(all_implementations[name], shape) for name in names)
for runtime, name in zip(runtimes, names):
assert runtime >= runtimes[names.index('pystencils')], runtimes
speedups = tuple(runtime / min(runtimes) for runtime in runtimes)
y_pos = np.arange(len(names))
labels = tuple(f"{name} ({round(speedup, 1)} x)" for name, speedup in zip(names, speedups))
plt.text(0.5, 0.5, f"Size {shape}", horizontalalignment='center', fontsize=16,
verticalalignment='center', transform=plt.gca().transAxes)
plt.barh(y_pos, runtimes, log=True)
plt.yticks(y_pos, labels);
plt.xlabel('Runtime of single iteration')
plt.figure(figsize=(8, 8))
plt.subplot(3, 1, 1)
bar_plot(16, 16)
plt.subplot(3, 1, 2)
bar_plot(128, 128)
plt.subplot(3, 1, 3)
bar_plot(1024, 1024)
```
%%%% Output: display_data
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAkQAAAHnCAYAAABZgEKHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XmclWX9//HXG1A2BRRIEcVxARVBEUcwDUJFMTWBoq+BPwUVy1yCEvcVxcK0chcVZUzQyhIDQUNj0wJ12CFcExMxRVlcWBT4/P449+CZmTMbzsJw3s/H4zy4z3Vfy+cej/CZ67ru+ygiMDMzM8tmdWo6ADMzM7Oa5oTIzMzMsp4TIjMzM8t6TojMzMws6zkhMjMzs6znhMjMzMyynhMiMzMzy3pOiMzMzCzrOSEyMzOzrFevpgOw6tOiRYvIycmp6TDMzMyqzZw5cz6OiJZl1XNClEVycnLIz8+v6TDMzMyqjaR3y1PPS2ZmZmaW9ZwQmZmZWdZzQmRmZmZZzwmRmZmZZT1vqs4ii95fS86Vkyq1z2UjT63U/szMzGqCZ4jMzMws6zkhMjMzs6znhMjMzMyynhMiMzMzy3pOiMzMzCzr1fqESNIgSXtVQb+jJbVPjpdJalGBtndI6p4cXyzpLUlRWh+SBkp6M3kNrGCsL0jarSJtzMzM7Gu1PiECBgGVnhBFxOCI+HdF20naHTg6ImYmRf8EegIlfpdK0uYGoCvQBbihggnOY8CFFY3VzMzMUmosIZKUI+k1SY9KWijpL5IaSTpB0vi0eidKekpSXUl5khZLWiTpF5L6AbnAOEnzJTWUdKSkGZLmSPq7pFZJP9Ml3SrpFUlvSOqWlNeVdHvS50JJl6TVzy0Sc2NJkyQtSOI4I8Ol9QOeK3gTEfMiYlkZP45ewPMRsSoiVgPPAycXGbuppNclHZS8f0LS+cnpCUD/MsYwMzOzEtT0DNFBwIMRcRjwKalZjqnAIZJaJnXOAcYAnYDWEdEhIjoCYyLiL0A+cGZEdAI2AXcD/SLiSOAR4Ja08epFRBdgKKkZGYCfAPsBRyRxjCsl3pOBFRFxeER0IC3xSXMsMKdCPwVoDbyX9n55UrZVRKwFLgbyJP0Y2C0iHkrOrQbqS2pewXHNzMyMmk+I3ouIfybHY4HvRESQWgL6f5KaAd8GngX+A+wv6W5JJ5NKoIo6COgAPC9pPnAtsHfa+aeSP+cAOclxT2BURGwCiIhVpcS7COiZzDR1S5KUoloBK0u76AyUoSyKFUQ8n8RwLzC4yOmPyLB0KOknkvIl5W9elylcMzMzq+mEqOg/+gXvxwD/j9Qy0JMRsSmZBTkcmA5cBIzO0J+AJRHRKXl1jIiT0s5vTP7czNdfW6IMcWQONuIN4EhSScmvJV2fodp6oEF5+kuzHNgn7f3ewIqilSTVAQ5Jxti9yOkGSXnRmB+MiNyIyK3bqGkFwzIzM8sONZ0QtZH07eS4P/ASQESsIJUQXAvkASR3aNWJiL8C1wGdk3afAbsmx68DLQv6lLSTpEPLiGEKcIGkekmboonGVsndbOsiYixwe1oM6ZYCB5YxZlF/B06StFuymfqkpKyoXyT99wcekbRTEpeAPYFlFRzXzMzMqPmEaCkwUNJCUjMe96edG0dqSa3gTq/WwPRkKSwPuCopzwNGJeV1SW1qvlXSAmA+cEwZMYwG/gssTNoMKKVuR+CVZKxrgBEZ6kwCehS8kfRzSctJzfoslDQ6Kc8tOE6W6W4GXk1eNxVdupPUjtQy2aUR8SIwk1TCCKlZq9kFy35mZmZWMUpt2amBgaUc4Jlkc3Km8/cA8yLi4eqMqzJIegk4LSLWVNN4dwITIuIfpdWr36pttBp4R6WO7W+7NzOz7ZmkORGRW1a9mp4hykjSHOAwUhuta6NLgTbVON7ispIhMzMzK1m9sqtUjeTZPBlnh5Jb5mutiHi5msd7qDrHMzMz29FslzNEZmZmZtXJCZGZmZllvRpbMrPq17F1U/K9CdrMzKwYzxCZmZlZ1nNCZGZmZlnPCZGZmZllPSdEZmZmlvW8qTqLLHp/LTlXTqqx8f1UazMz2155hsjMzMyynhMiMzMzy3pOiMzMzCzrOSEyMzOzrOeEyMzMzLJerU2IJO0p6Y+S3pb0b0mTJbUrpX4nSaekvb9R0rAqim2opLOT4x9JWiJpi6TcUto8IukjSYu3YbzbJR3/TWI2MzPLZrUyIZIkYDwwPSIOiIj2wNXAHqU06wScUsr5yoqtHnAu8HhStBj4ATCzjKZ5wMnbOOzdwJXb2NbMzCzr1cqECDgO+CoiRhUURMT8iHhR0mOSeheUSxon6XTgJuAMSfMlnZGcbi9puqT/SPp5WptfSlqcvIYmZTmSlkp6KJnxmSKpYYbYjgfmRsSmJK6lEfF6WRcUETOBVaXVkfS3tJmnn0oal7R9F2guac+yxjEzM7PiamtC1AGYU8K50cA5AJKaAscAk4HrgT9FRKeI+FNS92CgF9AFuEHSTpKOTNp3BY4Gzpd0RFK/LXBvRBwKrAF+mGH8Y0uJ7Zv6CXC9pG7ApcAlaefmJmObmZlZBdXWhKhEETEDOFDSt4D+wF8LZmsymBQRGyPiY+AjUktu3wHGR8QXEfE58BTQLan/TkTMT47nADkZ+mwFrKycqyksIj4kldhNAy6NiPQZpY+AvYq2kfQTSfmS8jevW1sVYZmZmdV6tTUhWgIcWcr5x4AzSc30jCml3sa0482kvspEFaxf1HqgQSl9fFMdgU8onvw0SMYuJCIejIjciMit26hpFYZlZmZWe9XWhGgqUF/S+QUFko6S9N3kbR4wFCAiliRlnwG7lqPvmUAfSY0kNQb6Ai9WILalwIEVqF9ukroA3wOOAIZJ2i/tdDtSG7jNzMysgmplQhQRQSpROTG57X4JcCOwIjn/IanEJH12aBqpTdTpm6oz9T2XVEL1CvAyMDoi5lUgvGeB7gVvJPWVtBz4NjBJ0t+T8r0kTU6r9wQwCzhI0nJJ56V3Kqk+8BBwbkSsILWH6BGl7EQqCcuvQJxmZmaWUCq32LFIagQsAjpHRLVvnJE0Hrg8It6spvH6krrW60qrV79V22g18I7qCCkjf9u9mZlVN0lzIqLE5wAWqJUzRKWR1BN4Dbi7JpKhxJWkNldXl3rAb6txPDMzsx1Kpk3BtVpEvAC0qeEYXgfKfPZQJY73ZHWNZWZmtiPa4WaIzMzMzCrKCZGZmZllPSdEZmZmlvV2uD1EVrKOrZuS7zu9zMzMivEMkZmZmWU9J0RmZmaW9ZwQmZmZWdbzHqIssuj9teRcOammwzDzU8vNbLvjGSIzMzPLek6IzMzMLOs5ITIzM7Os54TIzLZ7Tz/9NN27d+db3/oWDRs2ZN9996VPnz4899xzW+vk5eUhiWXLllVrbJdccgnf//73C5VdffXVnHTSSTRv3hxJ5OXlldh+9erVDB06lDZt2lC/fn323ntvBg0atE2xLF68mJ/+9KcceeSR7Lzzzkgqtf7s2bM5+eSTadasGY0bN6Zjx4788Y9/3Hp+3rx5NGrUiP/+97/bFI9ZbeKEyMy2a3fddRd9+/albdu2PPzww0yaNIlrr70WgKlTp26td+qppzJr1ixatWpVbbG9/fbbPPDAA9xwww2Fyu+++27Wr1/PaaedVmr71atX853vfIcXXniBESNG8Pzzz3P77bez6667blM8c+bMYfLkybRp04bc3NxS606aNInu3buz55578vjjj/O3v/2N888/nw0bNmytc8QRR3DiiSdy3XXXbVM8ZrWJIqKmY7BqUr9V22g18I6aDsOsQneZtWnThiOPPJLx48cXO7dlyxbq1Km53+suueQSZs+ezauvvlqovCCut956i7Zt2zJmzJiMsz4XXHABzz77LIsWLaJJkybfOJ70n8e1117LLbfcQqa/4z/77DMOOOAABgwYwB13lP53wuTJk+nduzfvvvsue+211zeO0ay6SZoTEaX/hkCWzxBJypG0ODnOlXRXJfV7h6TuyXGepHckzU9enUpo00bSFElLJf1bUk4FxvujpLaVEbvZ9mbVqlXsueeeGc+lJ0NFl8wGDRqEpIyv6dOnb223YMECTj/9dHbbbTcaNmzIsccey4svvlhmXBs3bmTs2LEMGDCg1LhK8sUXX/CHP/yBwYMHV0oyVN5xAZ588klWrlzJpZdeWmbdk046iSZNmpS67Ge2I8jqhChdRORHxM+/aT+SdgeOjoiZacWXRUSn5DW/hKZ/AG6LiEOALsBHFRj2fuDybYvYbPvWpUsXHn30UW677TbeeOONcre77rrrmDVrVqHXscceS6NGjWjTpg0Ac+fO5ZhjjmHVqlU89NBD/PWvf6V58+b07NmTOXPmlNr/7NmzWbNmDd26ddum65ozZw7r169njz32oF+/fjRs2JBddtmFPn368M4772xTn+X10ksvsfvuu7No0SI6duxIvXr12GeffRg+fDibN28uVLdevXp8+9vfLrRfy2xHVOkJUTLrslTSQ5KWJLMeDZNz0yXlJsctJC1LjgdJelrSxGQ25WJJv5Q0T9LsJMkoaH+HpH9JWiypi6Q6kt6U1DKpU0fSW5JaFInru2mzNPMk7VrkfA9JzyTHu0gaI2mRpIWSfpiUnyRplqS5kp6UtEuGH0E/oEJ/c0hqD9SLiOcBIuLziFhXpE49Sa9K6pG8/7WkW5LTLwI9JflBm7bDGTVqFAceeCCXX345Bx10EC1atKB///5MmTKl1HYHHHAARx999NbXSy+9xKxZsxg3bhz7778/AJdddhlt2rRh6tSp9OvXj1NOOYXx48ez//77c/PNN5fa/+zZs5HEYYcdtk3XtWLFCgCGDRtG3bp1mTBhAg8++CDz5s2jR48efPbZZ9vUb3nHXrduHQMGDGDQoEG88MILDBw4kJtvvplhw4YVq3/EEUfwyiuvsGXLliqLyaymVdUMUVvg3og4FFgD/LAcbToAA0jNjtwCrIuII4BZwNlp9RpHxDHAhcAjEbEFGAucmZzvCSyIiI+L9D8MuCgiOgHdgPWlxHIdsDYiOkbEYcDUJMG6FugZEZ2BfOCXGdoeCxT91fKWJLH6vaT6Gdq0A9ZIeipJ1m6TVDe9QkRsAgYB90s6ETgZGJ6c2wK8BRxetGNJP5GULyl/87q1pVyy2fapXbt2zJs3jxkzZnDNNdfQqVMnxo8fT69evRgxYkS5+pg4cSJXXHEFt956K3369AFg/fr1zJgxgx/96EfUqVOHTZs2sWnTJiKCnj17MnPmzFL7XLFiBU2aNGHnnXfepusqSC72228//vjHP3LiiScyYMAA/vznP/Pf//6XsWPHblO/5R17w4YNXH/99Vx66aX06NGDESNGcP7553Pvvfeydm3hvytatmzJxo0bWbVqVZXFZFbTqioheidtaWgOkFOONtMi4rOIWAmsBSYm5YuKtH8CIFmSaiKpGfAIXydN5wJjMvT/T+B3kn4ONEsSjJL0BO4teBMRq4GjgfbAPyXNBwYC+2Zo2wpYmfb+KuBg4Chgd+CKDG3qkUrShiX19ieV/BQSEUuAx0j9bM6NiC/TTn8EFNvxGBEPRkRuROTWbdS0hMs1277VrVuX7t27M2LECF544QX+85//0LFjR4YPH87q1atLbbtgwQIGDBjAeeedV2j2Y9WqVWzevJmbb76ZnXbaqdDrnnvuYfXq1aXOiGzYsIH69TP9flM+zZs3B6Bnz56Fbo/v2rUrTZo0Yd68edvcd3nHPvHEEwuVn3TSSXz11VcsWbKkUHnDhg2BVBJptqOqqiWWjWnHm4GGyfEmvk7CGpTSZkva+y0UjrPoLRMREe9J+lDS8UBXvp4tSq80UtIk4BRgtqSewIai9RLKMI6A5yOifwltCqwn7doi4oPkcKOkMaSSnqKWA/Mi4j8Akp4mlYA9nKFuR1KzbnsUKW9A6bNeZjuMvfbai8GDBzNkyBDefPNNunTpkrHehx9+yOmnn87RRx/NfffdV+hcs2bNqFOnDhdddBFnn312xvalbVJu3rx5mclYaQ499FCAEp8VVJV3z5U0dsEdaUXHLpgZatGi0E4Esx1KdW+qXgYcmRz328Y+zgCQ9B1Sy1oFc7ujSS2d/TkiNhdtJOmAiFgUEbeSWu46uJQxpgAXp7XdDZgNHCvpwKSskaR2GdouBQ5Ma9sq+VNAH2BxhjavArsV7IMCjgf+neEafgA0B7oDdyWzYwXaAUuKtjGr7d57772M5a+99hpAiXegbdiwgd69e9O4cWOefPJJ6tUr/Ptf48aN6datGwsWLKBz587k5uYWe5Xm4IMP5quvvmL58uXbcFWw9957k5uby5QpUwrdGj9r1iw+/fRTjjrqqG3qtzwKlg2LbpT++9//ToMGDejQoUOh8nfeeYd99tln60yR2Y6oujfh3g78WdJZwNSyKpdgtaR/AU1ILY8VmEBqqSzTchnAUEnHkZqx+jfwLKnlrUxGAPcmt+RvBoZHxFOSBgFPpO0DuhYoetvLJOCnpBI0gHFJoiNgPnABpG7zBy6IiMERsVnSMOAfSeI0B3govdNkD9NI4IRkRuwe4E5goKQ9gPVps1FmO4wOHTpw3HHH0bdvX/bbbz8+/fRTJk+ezKhRo/i///u/rXeMFTV06FDmzp1LXl7e1uSpQPv27WnSpAm/+93v6N69O7169eK8886jVatWfPzxx8ydO5fNmzczcuTIEuPq3r07AK+88gp77713oXMzZsxg5cqV/O9//wMgPz+fXXZJ3YPRr9/XvwuOHDmSXr160a9fPwYPHszKlSu55pprOPjggwvdzn/jjTcyfPhw3nnnHXJyckqMad26dUyePBn4OmH8y1/+AkBOTs7WJK9Dhw4MGjSI66+/ni1bttC5c2deeOEFRo8ezXXXXbc11gIvv/zy1us121HVqgczSpoODIuI/AzncoHfR8S23QNbiSS9BJwWEWuqabxfAJ9GRKYltq38YEbbXlTkwYyjRo1i8uTJLFiwgA8//JC6devSrl07+vfvz9ChQ7duas7Ly+Occ87ZmjT06NGDGTNmZOxz2rRp9OjRA4ClS5cyfPhwpk6dytq1a2nZsiWdO3fmggsu4JRTTik1tq5du9K+fXvGjCn8e1hpYxf9O/fZZ5/l+uuvZ9GiRTRu3JhTTz2V2267jT32+HpV/LLLLuPuu+/mf//7H82aNSva5VbLli1jv/32y3hu4MCBhZ4l9OWXX3LTTTfx6KOP8uGHH5KTk8NFF13EkCFDCrV777332HfffZkwYUKZT9422x6pnA9m3CESIklXAj8DzoyIl2oitiLxdCU1Y7OwmsY7B3isjI3iTohsu1GRhGh7lpeXx5AhQ/jggw9o1KhRlY1zzDHH0KlTp2L7oKrDrbfeyv3338/bb79N3bp1y25gtp0pb0JUqx7MGBE9Ms0ORcTIiNh3e0iGACLi5epKhpLxxpSVDJlZ5TvrrLNo3bp1lSYq69atY8GCBVxxRaYbVKvWhg0buPPOO7npppucDNkOr1YlRGZm25O6devyyCOPVOnsUKNGjfjiiy/Yd99MT/moWsuWLWPIkCGcddZZ1T62WXXzk43NzL6Bgidh74gOPvhgDj64tBtyzXYcToiySMfWTcnfQfZumJmZVSYvmZmZmVnWc0JkZmZmWc8JkZmZmWU9J0RmZmaW9bypOossen8tOVdOqukwzMxq3I7ycFCrPJ4hMjMzs6znhMjMzMyynhMiMzMzy3pOiMzMzCzrOSEyMzOzrOeEKANJPSQ9kxwPknRPCfX6SLq+SFk/SSEpN63sKklvSXpdUq8S+tpP0suS3pT0J0k7VyDeiyWdU976ZmZmVljWJkSSKuORA5cD96X1uSvwc+DltLL2wI+BQ4GTgfsk1c3Q163A7yOiLbAaOK8CcTySjGtmZmbboEoSIkk5kpZKekjSEklTJDVMzk0vmD2R1ELSsuR4kKSnJU2U9E4y6/FLSfMkzZa0e1r7OyT9S9JiSV0k1UlmVlomdeokMzItisR1o6QHJU0B/iCpgaQxkhYl4xxXgWtsB2yMiI/Tim8GfgNsSCvrDfwxIjZGxDvAW0CXIn0JOB74S1L0KNAnw5h3FcxISeolaaakOhGxDlgmqUvRNmZmZla2qpwhagvcGxGHAmuAH5ajTQdgAKmE4RZgXUQcAcwCzk6r1zgijgEuBB6JiC3AWODM5HxPYEGRZKXAkUDviBgAXAQQER2B/sCjkhqU8/qOBeYWvJF0BLBPRDxTpF5r4L2098uTsnTNgTURsamUOgBXAmckidtdwDnJtQPkA92KNpD0E0n5kvI3r1tbviszMzPLMlWZEL0TEfOT4zlATjnaTIuIzyJiJbAWmJiULyrS/gmAiJgJNJHUjNSyUUHSdC4wpoQxJkTE+uT4O8BjSV+vAe8C7coRJ0ArYCWkZqSA3wOXZqinDGWxDXVIZoLOB54H7omIt9NOfwTslaHNgxGRGxG5dRs1zXQdZmZmWa8qE6KNaceb+fprQjaljVt0Nia9zZa091so/DUjRZOFiIj3gA8lHQ90BZ4tIa4v0o4zJSLltZ6v49+V1OzW9GQJ8GhgQrI0uBzYJ63d3sCKIn19DDRL29eUqU6BjsAnFE9+GiQxmZmZWQXVxKbqZaSWrQD6bWMfZwBI+g6wNiIK1oJGk1o6+3NEbC5HPzNJltmSPUFtgNfLGcNS4ECAiFgbES0iIicicoDZwOkRkQ9MAH4sqb6k/UgtJb6S3lFEBDCNr38eA4G/FR1Q0r6kZqGOAL4nqWva6XbA4nLGbmZmZmlqIiG6HfiZpH8BLcqqXILVSftRFL4bawKwCyUvlxV1H1BX0iLgT8CgiNhYRpsCM4Ejkg3RJYqIJcCfgX8DzwEXFSRrkiZLKpjpuQL4paS3SO0peji9n2Sch4FhEbGC1HWPTtvzdCzwQjljNzMzszRKTU7UHpKmk0oK8jOcyyV163qxzcVVFMudwMSIqNFEJNnQ/cuIOKu0evVbtY1WA++opqjMzLZf/rb77CFpTkTkllVvh3kOkaQrgb8CV1XjsL8CGlXjeCVpAVxX00GYmZnVVpXxcMJqFRE9SigfCYys5lg+JLVMV6Mi4vmajsHMzKw222FmiMzMzMy2lRMiMzMzy3q1bsnMtl3H1k3J90ZCMzOzYjxDZGZmZlnPCZGZmZllPSdEZmZmlvWcEJmZmVnW86bqLLLo/bXkXDmppsMwMzPLqCafIO4ZIjMzM8t6TojMzMws6zkhMjMzs6znhMjMzMyynhMiMzMzy3rbVUIkabOk+ZIWS3pSUqNS6uZIGpD2fpCke6oorj6Srk+Ou0uaK2mTpH5pdfaVNCeJf4mkC9LO9Ze0SNJCSc9JapFhjN7J+fmS8iV9pwLxtZT03De9TjMzs2xV5QmRpIrc2r8+IjpFRAfgS+CCUurmAANKOV+ZLgfuS47/CwwCHi9S5wPgmIjoBHQFrpS0V3L9dwLHRcRhwELg4gxj/AM4PGl/LjC6vMFFxErgA0nHlv+SzMzMrECZCVEyE/OapEeTGYy/FMzcSFpWMNshKVfS9OT4RkkPSpoC/EFSXUm3SXo16eOn5YjtReBASTdLGpIWzy2Sfg6MBLolMyq/SE7vlczAvCnpN2ltCmZoFku6Na3886S/BZJmS9ojw/W3AzZGxMcAEbEsIhYCW9LrRcSXEbExeVufr3+2Sl6NJQloAqwoOk5EfB4RkbxtDETROpKOSn5+DSQ1TmaiOiSnnwbOzPiTNDMzs1KVd4boIODBZIbjU+DCcrQ5EugdEQOA84C1EXEUcBRwvqT9SmqYzKp8D1gEPAwMTMrrAD8GxgFXAi8mM0q/T5p2As4AOgJnSNpH0l7ArcDxyfmjJPVJ6jcGZkfE4cBM4PwM4RwLzC3H9ZKMtxB4D7g1IlZExFfAz5JrWQG0T64pU/u+kl4DJpGaJSokIl4FJgAjgN8AYyNicXI6H+hWnjjNzMyssPImRO9FxD+T47FAefa3TIiI9cnxScDZkuYDLwPNgbYZ2jRM6uSTWpp6OCKWAZ9IOiLpZ15EfFLCmP+IiLURsQH4N7AvqQRsekSsjIhNpJKp7kn9L4FnkuM5pJbhimoFrCzH9RIR7yVJ44HAQEl7SNqJVEJ0BLAXqSWzq0poPz4iDgb6ADeXMMxNwIlALqmkqMBHSf+FSPpJsicpf/O6teW5DDMzs6xT3v09RZdvCt5v4uukqkGROl+kHQu4JCL+XsY465M9NEWNJrVvZ0/gkVLab0w73kzq+lRK/a/SlqkK6heLCWhaSh/FRMQKSUtIzdi8m5S9DSDpz6Rmt0prP1PSAZJaFCzVpdkd2AXYidTPvODn3CCJtWhfDwIPAtRv1bbYMpyZmZmVf4aojaRvJ8f9gZeS42WklsYAflhK+78DP0tmS5DUTlLjCsQ5HjiZ1GxPQVL1GbBrOdq+DHxXUgtJdZP4Z1Rg7KWkZnxKJWlvSQ2T491ILbW9DrwPtJfUMql6YtJn0fYHJnuMkNQZ2BnINBP2IHAdqZmuW9PK2wGLM9Q3MzOzMpR3hmgpqSWgB4A3gfuT8uHAw5KuJpV4lGQ0qeWouck/+itJLQuVS0R8KWkasCYiNifFC4FNkhYAecDqEtp+IOkqYBqp2aLJEfG38o5Nam/RbyUpIkLSUaQStN2A70saHhGHAock9SIZ5/aIWAQgaTgwU9JXpGaMBiXlFyQxjiKVUJ6d1FkPnJE2e0VS/2xgU0Q8niR3/5J0fERMBY4jtffIzMzMKkhF/s0tXkHKAZ5JboWvEclm6rnAjyLizRoY/05gYkS8UN1jl5ekmaQ2sWdMDCG1ZNZq4B3VGJWZmVn5VcW33UuaExG5ZdXbrh7MmImk9sBbpDZMV3sylPgVUOJDImtashz3u9KSITMzMytZmUtmyV1eNTY7FBH/BvavqfGTGD4kdbv7dil5MOPTNR2HmZlZbbXdzxCZmZmZVTUnRGZmZpb1KvI9Y1bLdWzdlPwq2LBmZmZW23mGyMzMzLKeEyIzMzPLek6IzMzMLOs5ITIzM7Os54TIzMzMsl6ZX91hOw5JK0l9l9q2aAqsrcRwtscYKqv/b9LPtrTQxMIBAAAgAElEQVStSJvy1i1PvRbAx+Uct7by5756+tkePvf+zKfsiJ/5fSOiZZm1IsIvv8p8AQ/u6DFUVv/fpJ9taVuRNuWtW556QH5Nfyaq+uXPffX0sz187v2Zr9zPQ22MwUtmVl4TazoAqj6Gyur/m/SzLW0r0qa8dbeH/97bg+3h5+DP/TdvU56628N/6+3B9vBzqJEYvGRmZttEUn6U4xukzXYU/szv2DxDZGbb6sGaDsCsmvkzvwPzDJGZmZllPc8QmZmZWdZzQmRmZmZZzwmRmZmZZT0nRGZWqST1kPSipFGSetR0PGbVRVJjSXMknVbTsVjFOSEys60kPSLpI0mLi5SfLOl1SW9JurKMbgL4HGgALK+qWM0qSyV97gGuAP5cNVFaVfNdZma2laTupJKZP0REh6SsLvAGcCKpBOdVoD9QF/h1kS7OBT6OiC2S9gB+FxFnVlf8Ztuikj73h5H6ao8GpP4feKZ6orfKUq+mAzCz7UdEzJSUU6S4C/BWRPwHQNIfgd4R8WugtKWB1UD9qojTrDJVxude0nFAY6A9sF7S5IjYUqWBW6VyQmRmZWkNvJf2fjnQtaTKkn4A9AKaAfdUbWhmVaZCn/uIuAZA0iCSWdIqjc4qnRMiMyuLMpSVuNYeEU8BT1VdOGbVokKf+60VIvIqPxSrDt5UbWZlWQ7sk/Z+b2BFDcViVl38uc8yTojMrCyvAm0l7SdpZ+DHwIQajsmsqvlzn2WcEJnZVpKeAGYBB0laLum8iNgEXAz8HVgK/DkiltRknGaVyZ97A992b2ZmZuYZIjMzMzMnRGZmZpb1nBCZmZlZ1nNCZGZmZlnPCZGZmZllPSdEZmZmlvWcEJlZlZK0WdJ8SYslTZTU7Bv01UPSMWnvL5B0duVEWq7xu0laklxPw3LUHy2p/TaOlSNpcQXb/Cut7YBtGbeUvq/ONJbZjsLPIcoiLVq0iJycnJoOw8zMrNrMmTPn44hoWVY9f7lrFsnJySE/P7+mwzAzM6s2kt4tTz0vmZmZmVnWc0JkZmZmWc8JkZmZmWU97yHKIoveX0vOlZMqtc9lI0+t1P7MzMxqgmeIzMzMLOs5ITIzM7Os54TIzMzMsp4TIjMzM8t6TojMzMws6zkhMjMzs6xX6xMiSYMk7VUF/W79UkZJyyS1qEDbOyR1T44vlvSWpCitD0kDJb2ZvAZWMNYXJO1WkTZmZmb2tVqfEAGDgEpPiCJicET8u6LtJO0OHB0RM5OifwI9gRK/SyVpcwPQFegC3FDBBOcx4MKKxmpmZmYpNZYQScqR9JqkRyUtlPQXSY0knSBpfFq9EyU9JamupDxJiyUtkvQLSf2AXGCcpPmSGko6UtIMSXMk/V1Sq6Sf6ZJulfSKpDckdUvK60q6PelzoaRL0urnFom5saRJkhYkcZyR4dL6Ac8VvImIeRGxrIwfRy/g+YhYFRGrgeeBk4uM3VTS65IOSt4/Ien85PQEoH8ZY5iZmVkJanqG6CDgwYg4DPiU1CzHVOAQSS2TOucAY4BOQOuI6BARHYExEfEXIB84MyI6AZuAu4F+EXEk8AhwS9p49SKiCzCU1IwMwE+A/YAjkjjGlRLvycCKiDg8IjqQlvikORaYU6GfArQG3kt7vzwp2yoi1gIXA3mSfgzsFhEPJedWA/UlNS/asaSfSMqXlL953doKhmVmZpYdajohei8i/pkcjwW+ExFBagno/0lqBnwbeBb4D7C/pLslnUwqgSrqIKAD8Lyk+cC1wN5p559K/pwD5CTHPYFREbEJICJWlRLvIqBnMtPULUlSimoFrCztojNQhrIoVhDxfBLDvcDgIqc/IsPSYUQ8GBG5EZFbt1HTCoZlZmaWHWo6ISr6j37B+zHA/yO1DPRkRGxKZkEOB6YDFwGjM/QnYElEdEpeHSPipLTzG5M/N/P197gpQxyZg414AziSVFLya0nXZ6i2HmhQnv7SLAf2SXu/N7CiaCVJdYBDkjF2L3K6QVJuZmZmFVTTCVEbSd9OjvsDLwFExApSCcG1QB5AcodWnYj4K3Ad0Dlp9xmwa3L8OtCyoE9JO0k6tIwYpgAXSKqXtCmaaGyV3M22LiLGArenxZBuKXBgGWMW9XfgJEm7JZupT0rKivpF0n9/4BFJOyVxCdgTWFbBcc3MzIyaT4iWAgMlLSQ143F/2rlxpJbUCu70ag1MT5bC8oCrkvI8YFRSXpfUpuZbJS0A5gPHlBHDaOC/wMKkzYBS6nYEXknGugYYkaHOJKBHwRtJP5e0nNSsz0JJo5Py3ILjZJnuZuDV5HVT0aU7Se1ILZNdGhEvAjNJJYyQmrWaXbDsZ2ZmZhWj1JadGhhYygGeSTYnZzp/DzAvIh6uzrgqg6SXgNMiYk01jXcnMCEi/lFavfqt2kargXdU6tjLRp5aqf2ZmZlVJklzIiK3rHo1PUOUkaQ5wGGkNlrXRpcCbapxvMVlJUNmZmZWsnplV6kaybN5Ms4OJbfM11oR8XI1j/dQdY5nZma2o9kuZ4jMzMzMqpMTIjMzM8t6NbZkZtWvY+um5HsTtJmZWTGeITIzM7Os54TIzMzMsp4TIjMzM8t6TojMzMws63lTdRZZ9P5acq6cVNNhbPf89G0zs+zjGSIzMzPLek6IzMzMLOs5ITIzM7Os54TIzMzMsp4TIjMzM8t6tTYhkrSnpD9KelvSvyVNltSulPqdJJ2S9v5GScOqKLahks5Ojm+T9JqkhZLGS2qWoX4DSa9IWiBpiaThFRzvj5LaVlb8ZmZm2aZWJkSSBIwHpkfEARHRHrga2KOUZp2AU0o5X1mx1QPOBR5Pip4HOkTEYcAbwFUZmm0Ejo+Iw5M4T5Z0dAWGvR+4fNujNjMzy261MiECjgO+iohRBQURMT8iXpT0mKTeBeWSxkk6HbgJOEPSfElnJKfbS5ou6T+Sfp7W5peSFievoUlZjqSlkh5KZnGmSGqYIbbjgbkRsSmJa0rBMTAb2Ltog0j5PHm7U/KK9DqS6kl6VVKP5P2vJd2SnH4R6JkkY2ZmZlZBtTUh6gDMKeHcaOAcAElNgWOAycD1wJ8iolNE/CmpezDQC+gC3CBpJ0lHJu27AkcD50s6IqnfFrg3Ig4F1gA/zDD+saXEdi7wbKYTkupKmg98BDwfES+nn0+SqkHA/ZJOBE4GhifntgBvAYdn6PcnkvIl5W9et7aEsMzMzLJbbU2IShQRM4ADJX0L6A/8NW2GpqhJEbExIj4mlYjsAXwHGB8RXySzNk8B3ZL670TE/OR4DpCToc9WwMqihZKuATYB40qIe3NEdCI1g9RFUocMdZYAjwETgXMj4su00x8Be2Vo82BE5EZEbt1GTTMNbWZmlvVqa0K0BDiylPOPAWeSmukZU0q9jWnHm0l9lYkqWL+o9UCD9AJJA4HTgDMjIjK02Soi1gDTSc0AZdKR1OxU0f1SDZKxzczMrIJqa0I0Fagv6fyCAklHSfpu8jYPGApbZ1UAPgN2LUffM4E+khpJagz0JbVHp7yWAgemxXUycAVwekSsy9RAUsuCu8+SfUk9gdcy1PsB0BzoDtxV5I61dqQSRTMzM6ugWpkQJbMsfYETk9vulwA3AiuS8x+SSkzSZ4emkdpEnb6pOlPfc0klVK8ALwOjI2JeBcJ7llTCUuAeUonY88nYowAk7SVpclKnFTBN0kLgVVJ7iJ5J71RSC2AkcF5EvJH0e2dybg9gfUR8UIE4zczMLKEyVnBqJUmNgEVA54io9p3EksYDl0fEm9U03i+ATyPi4dLq1W/VNloNvKM6QqrV/G33ZmY7DklzIiK3rHq1coaoNJIKlpvurolkKHElqVmf6rIGeLQaxzMzM9uh7HDPrYmIF4A2NRzD68Dr1TheaRvHzczMrAw73AyRmZmZWUU5ITIzM7Ost8MtmVnJOrZuSr43DJuZmRXjGSIzMzPLek6IzMzMLOs5ITIzM7Os54TIzMzMsp43VWeRRe+vJefKSTUdRq3hJ1abmWUPzxCZmZlZ1nNCZGZmZlnPCZGZ1Zinn36a7t27861vfYuGDRuy77770qdPH5577rmtdfLy8pDEsmXLqjW2Sy65hO9///uFyq6++mpOOukkmjdvjiTy8vKKtfvggw+46qqryM3NpWnTprRs2ZITTjiBmTNnFqu7efNmfv/739OhQwcaN25Mq1at6Nu3LwsXLtymmD/77DOGDRtGjx49aNKkCZKYPn16sXpvvPEGQ4YM4bDDDmOXXXahVatWnH766SxYsKBY3XXr1nHDDTfQrl07GjZsyD777MPZZ59d7L9H7969ueiii7YpbrPtgRMiM6sRd911F3379qVt27Y8/PDDTJo0iWuvvRaAqVOnbq136qmnMmvWLFq1qr7vS3777bd54IEHuOGGGwqV33333axfv57TTjutxLZz5szhT3/6E7179+Yvf/kLeXl5NGjQgB49evDMM88UqnvdddcxbNgw+vTpw8SJE7nzzjt5++23Oe6441i+fHmF4/7kk0945JFHqFevHieeeGKJ9aZMmcK0adMYOHAgEydO5L777mPlypV07dqVOXPmFKo7ePBgbrvtNs4//3wmT57MiBEjmDlzJieccAKff/751no33ngjDz30EG+88UaF4zbbHigiajoGqyb1W7WNVgPvqOkwag1vqq5abdq04cgjj2T8+PHFzm3ZsoU6dWru97VLLrmE2bNn8+qrrxYqL4jrrbfeom3btowZM4ZBgwYVqrNmzRp22WUX6tX7+p6VTZs2ceihh7LHHnsUminaa6+96NGjB48//vjWstdee41DDjmEUaNG8dOf/rRCcUcEkgB44YUXOPHEE5k2bRo9evQoVO/jjz/eOstVYO3ateTk5PD973+fP/zhDwCsX7+eXXfdlcsvv5xf/epXW+s+99xzfO973+O5556jV69eW8u7dOlCbm4u9913X4XiNqtKkuZERG5Z9bJ6hkhSjqTFyXGupLsqqd87JHVPjsdJel3SYkmPSNqplHZNJL0v6Z4KjveCpN2+adxm1WnVqlXsueeeGc+lJ0NFl8wGDRqEpIyv9OWhBQsWcPrpp7PbbrvRsGFDjj32WF588cUy49q4cSNjx45lwIABpcZVkmbNmhVKhgDq1atHp06deP/99wuVf/nllzRp0qRYe0glXxWVnuCUpkWLFsXqNm3alHbt2hWKcdOmTWzevLncMf74xz9m3LhxrF+/vsKxm9W0rE6I0kVEfkT8/Jv2I2l34OiIKPg1cBxwMNARaAgMLqX5zcCMbRj2MeDCbWhnVmO6dOnCo48+ym233VahZZbrrruOWbNmFXode+yxNGrUiDZt2gAwd+5cjjnmGFatWsVDDz3EX//6V5o3b07Pnj2LLQkVNXv2bNasWUO3bt2+0fWl+/LLL5k1axaHHHJIofILL7yQsWPH8re//Y1PP/2U//znP1x44YXsvffenHHGGZU2fnmsWrWKxYsXF4px11135ayzzuKuu+5i2rRpfP755yxZsoTLLruMww8/nBNOOKFQH927d+fTTz9l1qxZ1Rq7WWWo9IQomXVZKukhSUskTZHUMDk3XVJuctxC0rLkeJCkpyVNlPSOpIsl/VLSPEmzkySjoP0dkv6VzLh0kVRH0puSWiZ16kh6S1KLInF9V9L85DVP0q5FzveQ9ExyvIukMZIWSVoo6YdJ+UmSZkmaK+lJSbtk+BH0A7buCI2IyZEAXgH2LuHndiSwBzClhPNNk5mmg5L3T0g6Pzk9Aehfwn8Ss+3SqFGjOPDAA7n88ss56KCDaNGiBf3792fKlIz/C2x1wAEHcPTRR299vfTSS8yaNYtx48ax//77A3DZZZfRpk0bpk6dSr9+/TjllFMYP348+++/PzfffHOp/c+ePRtJHHbYYZV2rTfeeCPLly/niiuuKFR+0003cdVVV/GDH/yApk2bcsABB7BkyRKmT5/O7rvvXmnjl8cll1xCRDB06NBC5WPGjKFv374cf/zx7LrrrnTo0IGvvvqK559/np133rlQ3cMPP5w6deowe/bs6gzdrFJU1QxRW+DeiDgUWAP8sBxtOgADgC7ALcC6iDgCmAWcnVavcUQcQ2pG5JGI2AKMBc5MzvcEFkTEx0X6HwZcFBGdgG5AaXO61wFrI6JjRBwGTE0SrGuBnhHRGcgHfpmh7bFAsV9Bk6Wys0hLltLO1QF+C1xWUkARsRa4GMiT9GNgt4h4KDm3GqgvqXkp12S2XWnXrh3z5s1jxowZXHPNNXTq1Inx48fTq1cvRowYUa4+Jk6cyBVXXMGtt95Knz59gNS+lxkzZvCjH/2IOnXqsGnTJjZt2kRE0LNnz4x3e6VbsWIFTZo0KfaP/bZ6/PHHGTlyJNddd12xWaf777+fESNGcO211zJt2jSefPJJdt11V0466SRWrFhRKeOXx69//Wsef/xx7rnnHg488MBC56699lrGjh3L7bffzowZM3jsscf45JNP+N73vscXX3xRqO5OO+1E06ZNqzV2s8pSVU+qfici5ifHc4CccrSZFhGfAZ9JWgtMTMoXAem/qj0BEBEzkz03zYBHgL8BdwDnAmMy9P9P4HeSxgFPRcTyUtbbewI/LngTEaslnQa0B/6ZtNuZVLJWVCtgZYby+4CZEZFpE8OFwOSIeK+0PQAR8bykHwH3AocXOf0RsBfwSXqhpJ8APwGo26RliX2b1YS6devSvXt3unfvDqSSkZNPPpnhw4dz0UUXsdtuJW+NW7BgAQMGDOC8885j2LBhW8tXrVrF5s2bufnmm0ucDSpt0/aGDRuoX7/+N7iqr02cOJFBgwZx3nnnMXz48ELnVq1axS9+8Qsuu+yyQueOP/54cnJyuO222/j9739fKXGUZtSoUVx99dWMGDGCc889t9C5JUuWMHLkSEaPHs155523tbxr1660a9eO0aNHM2TIkEJtGjZs6D1EVitVVUK0Me14M6m9MwCb+HpWqkEpbbakvd9C4TiL3hYXSSLxoaTjga58PVuUXmmkpEnAKcBsST2BDSXErwzjCHg+IspamlpPkWuTdAPQEijplpFvA90kXQjsAuws6fOIuLJIP3WAQ5IxdgfS78ttQIZZr4h4EHgQUneZlRG7WY3aa6+9GDx4MEOGDOHNN9+kS5cuGet9+OGHnH766Rx99NHF7mhq1qwZderU4aKLLuLss8/O2L60zdHNmzdn9erV234RiX/84x/86Ec/om/fvjzwwAPFzr/xxhts3LiRo446qlD57rvvzgEHHMDSpUu/cQxleeyxx7jwwgu59NJLueaaa4qdX7RoEUCxGNu2bUuzZs0yxrhq1SpatGhRrNxse1fdm6qXAUcmx/22sY8zACR9h9Sy1tqkfDSppbM/R8Tmoo0kHRARiyLiVlLLXQeXMsYUUstTBW13A2YDx0o6MClrJKldhrZLgQPT2g4GegH9k+W9YiLizIhoExE5pJb2/lA0GUr8Ium/P7D1jjWlppX2JPXzNasV3nvvvYzlr732GkCJd6Bt2LCB3r1707hxY5588slid3Q1btyYbt26sWDBAjp37kxubm6xV2kOPvhgvvrqq216DlCBWbNm0bt3b0444QTGjh2bMQEruL5XXnmlUPmqVat46623aN269TaPXx7jx4/nnHPOYfDgwdx+++0Z65QU4xtvvMGaNWuKxfi///2PDRs2cNBBB1VN0GZVqLq/3PV24M+SzgKmllW5BKsl/QtoQmp5rMAEUktlmZbLAIZKOo7UjNW/gWdJLW9lMgK4N7klfzMwPCKekjQIeEJSwXz6tUDR22MmkZoJGp28HwW8C8xKlsOeioibks3lF0REaXedbZUkX4OBLhHxmaSZyfg3kEoyZ0fEpvL0ZbY96NChA8cddxx9+/Zlv/3249NPP2Xy5MmMGjWK//u//9t6x1hRQ4cOZe7cueTl5W1Nngq0b9+eJk2a8Lvf/Y7u3bvTq1cvzjvvPFq1asXHH3/M3Llz2bx5MyNHjiwxroLlu1deeYW99y58D8SMGTNYuXIl//vf/wDIz89nl11S91b065f6He+1117j1FNPpUWLFlx22WXF7mo7+uijAcjJyeG0007jtttuo06dOnz3u9/lk08+4Te/+Q0bN27kZz/72dY206dP57jjjsv43KOinn32Wb744outszszZszg448/pnHjxnzve98DYObMmfTv35/DDjuMQYMGFdoEXb9+fY444ggAunXrxuGHH86ll17K6tWryc3N5b///S8jRoygadOmDBw4sNDYL7/8cqGfoVltUqsezChpOjAsIvIznMsFfh8RlXev7DaS9BJwWkSsqabx7gQmRMQ/SqvnBzNWjB/MWLVGjRrF5MmTWbBgAR9++CF169alXbt29O/fn6FDh27d1JyXl8c555zDO++8Q05ODj169GDGjMxPp0h/COHSpUsZPnw4U6dOZe3atbRs2ZLOnTtzwQUXcMopp5QaW9euXWnfvj1jxhT+/aq0sQv+Li2ItyTpf+euW7eO3/72tzzxxBO8++67NGnShM6dO3PDDTcUWi6cNGkSp512Gs8++ywnn3xyqbHn5OTw7rvvFivfd999tz7L6cYbbyy2pylTPUg9/fpXv/oVEyZMYPny5bRo0YJjjjmGm266qdhM0Pnnn8+8efPIzy/2V7RZjSnvgxl3iIRI0pXAz4AzI+KlmoitSDxdgfURsW1fSFTx8c4vuOOsNE6IKsYJUfbKy8tjyJAhfPDBBzRq1Kimw+Hqq69mwoQJLFq0qNwPX6xuGzZsoFWrVtx+++2FNmCb1bTyJkS16sGMEdEj0+xQRIyMiH23h2QIICJerq5kKBmvzGTIzMrvrLPOonXr1tvNV1DMmDGDq6++ertNhgAeeOABvvWtbxVbRjOrLap7D5GZ2Xavbt26PPLII8ydO7emQwHgn//8Z02HUKb69euTl5dXbJO7WW1Rq5bM7JvxklnFeMnMzKz2K++SmVP5LNKxdVPy/Y+8mZlZMbVqD5GZmZlZVXBCZGZmZlnPCZGZmZllPSdEZmZmlvW8qTqLLHp/LTlXTqrpMGw74bvozMy+5hkiMzMzy3pOiMzMzCzrOSEyMzOzrOeEyMzMzLKeEyIzMzPLek6IKkBSD0nPJMeDJN1TQr0+kq5Pq7dS0vzkNbiENv0lLZK0UNJzklpUIK7bJR2/LddkZmZmToiKkVQZjyK4HLgv7f2fIqJT8hpdwph3AsdFxGHAQuDiCox3N3DlNwnYzMwsm1VrQiQpR9JSSQ9JWiJpiqSGybnpknKT4xaSliXHgyQ9LWmipHckXSzpl5LmSZotafe09ndI+pekxZK6SKoj6U1JLZM6dSS9VXT2RdKNkh6UNAX4g6QGksYkMzbzJB1XgWtsB2yMiI8r8qNJXo0lCWgCrMjQ998knZ0c/1TSOICIeBdoLmnPCoxpZmZmiZqYIWoL3BsRhwJrgB+Wo00HYADQBbgFWBcRR/z/9u49Tus5///446nohKJYyWFYhZRCck5WTrt+YvElFq1YZ6tdS+yyrPWVZb/rsEiR7BbLsmw5HxO2MKWD5BCycozIKVJevz+u99Q111xzzUzNzNV0Pe+329z6XJ/36fWZ+dxmXr3f7891AROBY7PqtYmIXYFTgZER8T0wGjg6lfcDplWTrOwA9I+Io4DTACKiOzAAuFVSy1pe327AlJxzh6alsLskbZzbICK+A04BZpBJhLoCN+fp+xfAhZL2AH4NnJFVNiWNbWZmZnVUjITorYiYmo4nA2W1aPNkRHwREfOABcC4dH5GTvvbASJiArC2pHbASJYlTccDt1QzxtiIWJiOdwf+nvp6BXgb6FKLOAE6AvOyXo8DytJS2GPArbkNJK1OJiHaDtiQzJLZebn1IuJD4ELgSeDXETE/q/ij1Da3719IKpdUvuTrBbW8BDMzs9JSjITo26zjJSz7+JDFLIsndzYmu833Wa+/p/LHj0ROu4iId4AP06bjnYAHq4nrq6xjVRt9zRaSFX9EfBIRFfGOIDMTlatnqvtGRARwJ7BrNf13Bz6havLTMo1dSUQMj4heEdGrWeu2dboQMzOzUrEybaqew7Jk4bDl7OMIAEm7AwsiomJK5CYyS2d3RsSSWvQzgbTMlvYEbQK8WssYZgFbVLyQ1DGr7KBUnutdoGvFXidgn3z1JPUGDiAzk3S2pM2yirsAL9UyRjMzM8uyMiVEVwKnSPoPUOtHznN8mtoPAwZlnR8LrEn1y2W5rgeaSZoB3AEMzJrlqckEYLu0ORrgzLSBfBpwJjCwoqKkqQAR8R5wMTBB0nQyM0b/m92ppBZkZpiOT/V/DYxUxupkkrDyWsZoZmZmWZRZoWn6JI0Hzo6IKklBenrtLxGxRyPFcjUwLiIea6TxDgG2j4gLCtVr0bFzdDzuqsYIyZoAf9q9mZUCSZMjoldN9VamGaIGIWkIcDd5Nik3oP8FWjfieM2BPzfieGZmZquU+ngTwpVCRPSt5vxQYGgjx/IhmWW6xhrvn401lpmZ2apolZ8hMjMzM6uJEyIzMzMreavMkpnVrHuntpR7I62ZmVkVniEyMzOzkueEyMzMzEqeEyIzMzMreU6IzMzMrOR5U3UJmfHuAsqG3F/sMMzMzPIq5jvoe4bIzMzMSp4TIjMzMyt5TojMzMys5DkhMjMzs5LnhMjMzMxKnhMiMzMzK3krTUIkaYmkqZJekvRPSa0L1C2TdFTW64GS/tpAcR0s6cKscealOKdKOiGr3uUp9pckHZGnn2slfVngehZm9TusjjE+Jmmdul6bmZmZZTRoQiSpLu9ztDAiekZEN2ARcHKBumXAUQXK69M5wPVZr+9IcfaMiJsAJP0E2B7oCewE/EbS2hUNJPUC2tUwzhtZ/Ra69nz+DpxaxzZmZmaWFEyI0szFK5JulTRd0l0VMzeS5kjqkI57SRqfji+SNFzSI8DfJDWTdOd0l3UAACAASURBVIWkF1IfJ9UirqeBLSRdIumXWfFcKulMYCiwR5pNGZyKN5T0kKTXJf0pq80ASTPSzM3lWee/TP1NkzRJ0g/yXH8X4NuI+LiGeLsCT0XE4oj4CpgG7J/6aAZcQSaxWm6S2kp6VdKW6fXtkk5MxWOBASvSv5mZWSmrzQzRlsDwiNgW+JzazUTsAPSPiKOAQcCCiNgR2BE4UdJm1TVMs0oHADOAm4Hj0vnVgCOBMcAQ4Ok0m/KX1LQncATQHThC0saSNgQuB36UyneUdHCq3waYFBE9gAlARXKRbTdgSs65Q7OSw43TuWnAAZJapyRxL6Ci7HRgbES8X+gbBmwm6UVJT0naI7cwIhakvkZJOhJYJyJGpLJPgRaS2ue2k/QLSeWSypd8vaCGEMzMzEpTbRKidyLi2XQ8Gti9Fm3GRsTCdLwvcKykqcBzQHugc542rVKdcuC/wM0RMQf4RNJ2qZ8XI+KTasZ8PCIWRMQ3wMvApmQSsPERMS8iFpNJpvqk+ouA+9LxZDLLcLk6AvOyXo8DylJy+BhwK0BEPAI8APwHuB2YCCxOCdnhwLXVxFzhfWCTiNgO+BVwW/aSW4WIeJRMongdcEJO8UfAhnnaDI+IXhHRq1nrtjWEYWZmVppqs8cnqnm9mGUJVcucOl9lHQs4IyIermGchRHRM8/5m4CBwAbAyALtv806XkLm2lSg/ncRETn1q8QELM0icpKxEWRmnyrKLgUuBZB0G/A6sB2wBTBbEkBrSbMjYovsQSLi24r4I2KypDeALmSSw6XSLNnWKa51gblZxS3TeTMzM6uj2swQbSJpl3Q8AHgmHc8hszQGcGiB9g8Dp0haHTL7ciS1qUOM95DZj7Nj6gvgC2CtWrR9DthTUoe0l2cA8FQdxp5FJqEBQFLHrLKDUjlpn1T7dLwtsC3wSETcHxEbRERZRJQBX+cmQ6nNeik+JG1OZgbtzTzxDE5jDgBGZn1PRSZhnFOHazMzM7OkNjNEs4DjJN1IZtbjhnT+YuBmSeeTSTyqcxOZ5agp6Q/3PODgAvUriYhFkp4EPouIJen0dDJLUtOAUcCn1bR9X9J5wJNkZoseiIh/13ZsMnuL/ixJaTbpTEkHkZkdm09m5gpgdeDpNAv0OfCztERXrdRPr4i4kMwy3h8kLSYzW3VyRMzPqd+FzDJZ74j4QtIE4HfA78kkppNqGtPMzMzy07JVozyFUhlwX3oUvijSMtEU4PCIeL0I418NjIuIxxp77NpKMY6NiMcL1WvRsXN0PO6qRorKzMysbuYM/Um99ylpckT0qqneSvPGjPlI6grMJrNhutGToeR/gWrfJHIl8VJNyZCZmZlVr+CSWXrKq2izQxHxMrB5scZPMXxI5n1+VloVj9+bmZnZ8lmpZ4jMzMzMGoMTIjMzMyt5dfmsMWviundqS3kDbFgzMzNr6jxDZGZmZiXPCZGZmZmVPCdEZmZmVvKcEJmZmVnJc0JkZmZmJa/gR3fYqkXSPODtFeymLbCgEdrVtn5N9QqVFyrrAHxci/FXFsv7cynWOKvSfVSovKndR9A491J9jtEY91J93Uc11fHvpIYZZ9OIWK/GniLCX/6q9RcwvDHa1bZ+TfUKlddQVl7s73Vj/FyKNc6qdB8VKm9q91F9/owba4zGuJfq6z6qqY5/JxV3HC+ZWV2Na6R2ta1fU71C5ct7LSujxrqW+hpnVbqP6jJOU9AY11KfYzTGvVRf91FNdXwfFXEcL5mZ5SGpPGrx6chmhfg+svrie6nheYbILL/hxQ7AVgm+j6y++F5qYJ4hMjMzs5LnGSIzMzMreU6IzMzMrOQ5ITIzM7OS54TIbDlIaiNpsqQDix2LNU2StpY0TNJdkk4pdjzWdEk6WNIISf+WtG+x42mqnBBZSZE0UtJHkl7KOb+/pFclzZY0pBZdnQvc2TBR2squPu6jiJgVEScD/wP4ceoSVU/30r0RcSIwEDiiAcNdpfkpMyspkvoAXwJ/i4hu6Vwz4DVgH2Au8AIwAGgGXJbTxfHAtmTeRr8l8HFE3Nc40dvKoj7uo4j4SNJBwBDgrxFxW2PFbyuP+rqXUrs/A2MiYkojhb9KaV7sAMwaU0RMkFSWc7o3MDsi3gSQ9A+gf0RcBlRZEpO0F9AG6AoslPRARHzfoIHbSqU+7qPUz1hgrKT7ASdEJaieficJGAo86GRo+TkhMoNOwDtZr+cCO1VXOSJ+CyBpIJkZIidDBnW8jyT1BX4KtAAeaNDIrKmp070EnAH0A9pK2iIihjVkcKsqJ0RmoDznalxLjohR9R+KNWF1uo8iYjwwvqGCsSatrvfSNcA1DRdOafCmarPM/742znq9EfBekWKxpsv3kdUX30tF4ITILLNhsbOkzSStARwJjC1yTNb0+D6y+uJ7qQicEFlJkXQ7MBHYUtJcSYMiYjFwOvAwMAu4MyJmFjNOW7n5PrL64ntp5eHH7s3MzKzkeYbIzMzMSp4TIjMzMyt5TojMzMys5DkhMjMzs5LnhMjMzMxKnhMiMzMzK3n+6I4S0qFDhygrKyt2GGZmZo1m8uTJH0fEejXVc0JUQsrKyigvLy92GGZmZo1G0tu1qeclMzMzMyt5TojMzMys5DkhMjMzs5LnhMjMzMxKnjdVl5AZ7y6gbMj99drnnKE/qdf+zMzMisEzRGZmZlbynBCZmZlZyXNCZGZmZiXPCZGZmZmVPCdEZmZmVvKafEIkaaCkDRug35skdU3HcyR1qEPbqyT1ScenS5otKQr1Iek4Sa+nr+PqGOtjktapSxszMzNbpsknRMBAoN4Toog4ISJerms7SesCO0fEhHTqWaAfUO1nqaQ2vwd2AnoDv69jgvN34NS6xmpmZmYZRUuIJJVJekXSrZKmS7pLUmtJe0u6J6vePpL+JamZpFGSXpI0Q9JgSYcBvYAxkqZKaiVpB0lPSZos6WFJHVM/4yVdLul5Sa9J2iOdbybpytTndElnZNXvlRNzG0n3S5qW4jgiz6UdBjxU8SIiXoyIOTV8O/YDHo2I+RHxKfAosH/O2G0lvSppy/T6dkknpuKxwIAaxjAzM7NqFHuGaEtgeERsC3xOZpbjCWBrSeulOj8HbgF6Ap0ioltEdAduiYi7gHLg6IjoCSwGrgUOi4gdgJHApVnjNY+I3sBZZGZkAH4BbAZsl+IYUyDe/YH3IqJHRHQjK/HJshswuU7fBegEvJP1em46t1RELABOB0ZJOhJYJyJGpLJPgRaS2tdxXDMzM6P4CdE7EfFsOh4N7B4RQWYJ6GeS2gG7AA8CbwKbS7pW0v5kEqhcWwLdgEclTQV+B2yUVf6v9O9koCwd9wOGRcRigIiYXyDeGUC/NNO0R0pScnUE5hW66DyU51xUORHxaIrhOuCEnOKPyLN0KOkXksollS/5Ol+4ZmZmVuyEKPePfsXrW4CfkVkG+mdELE6zID2A8cBpwE15+hMwMyJ6pq/uEbFvVvm36d8lLPvYEuWJI3+wEa8BO5BJSi6TdGGeaguBlrXpL8tcYOOs1xsB7+VWkrQasHUaY92c4pbpfG7MwyOiV0T0ata6bR3DMjMzKw3FTog2kbRLOh4APAMQEe+RSQh+B4wCSE9orRYRdwMXANundl8Aa6XjV4H1KvqUtLqkbWqI4RHgZEnNU5vcRGOp9DTb1xExGrgyK4Zss4Atahgz18PAvpLWSZup903ncg1O/Q8ARkpaPcUlYANgTh3HNTMzM4qfEM0CjpM0ncyMxw1ZZWPILKlVPOnVCRiflsJGAeel86OAYel8MzKbmi+XNA2YCuxaQww3Af8Fpqc2RxWo2x14Po31W+CPeercD/SteCHpTElzycz6TJd0Uzrfq+I4LdNdAryQvv6Qu3QnqQuZZbJfR8TTwAQyCSNkZq0mVSz7mZmZWd0os2WnCANLZcB9aXNyvvK/Ai9GxM2NGVd9kPQMcGBEfNZI410NjI2IxwvVa9Gxc3Q87qp6Hdufdm9mZiszSZMjoldN9Yo9Q5SXpMnAtmQ2WjdFvwY2acTxXqopGTIzM7PqNa+5SsNI782Td3YoPTLfZEXEc4083ojGHM/MzGxVs1LOEJmZmZk1JidEZmZmVvKKtmRmja97p7aUexO0mZlZFZ4hMjMzs5LnhMjMzMxKnhMiMzMzK3lOiMzMzKzkeVN1CZnx7gLKhtxf7DCKyu+sbWZm+XiGyMzMzEqeEyIzMzMreU6IzMzMrOQ5ITIzM7OS54TIzMzMSl6TTYgkbSDpH5LekPSypAckdSlQv6ekH2e9vkjS2Q0U21mSjk3HV0h6RdJ0SfdIaldNm/0lvSpptqQhdRzvH5I610fsZmZmpahJJkSSBNwDjI+IH0ZEV+B84AcFmvUEflygvL5iaw4cD9yWTj0KdIuIbYHXgPPytGkGXAccAHQFBkjqWodhbwDOWZG4zczMSlmTTIiAvYDvImJYxYmImBoRT0v6u6T+FecljZF0EPAH4AhJUyUdkYq7Shov6U1JZ2a1+ZWkl9LXWelcmaRZkkZIminpEUmt8sT2I2BKRCxOcT1ScQxMAjbK06Y3MDsi3oyIRcA/gP7ZFSQ1l/SCpL7p9WWSLk3FTwP9UjJmZmZmddRUE6JuwORqym4Cfg4gqS2wK/AAcCFwR0T0jIg7Ut2tgP3IJCS/l7S6pB1S+52AnYETJW2X6ncGrouIbYDPgEPzjL9bgdiOBx7Mc74T8E7W67np3FIpqRoI3CBpH2B/4OJU9j0wG+hRzbhmZmZWQFNNiKoVEU8BW0haHxgA3J01Q5Pr/oj4NiI+Bj4is+S2O3BPRHwVEV8C/wL2SPXfioip6XgyUJanz47AvNyTkn4LLAbG5GmjfJeS59pmAn8HxgHHp9mkCh8BG+YZ9xeSyiWVL/l6QZ5hzMzMrKkmRDOBHQqU/x04msxMzy0F6n2bdbyEzEeZ5EtOCtXPtRBomX1C0nHAgcDREVEl0SEzI7Rx1uuNgPeqiaE7mdmp3P1SLdPYlUTE8IjoFRG9mrVuW02XZmZmpa2pJkRPAC0knVhxQtKOkvZML0cBZ8HSWRWAL4C1atH3BOBgSa0ltQEOIbNHp7ZmAVtkxbU/cC5wUER8XU2bF4DOkjaTtAZwJDA2t5KknwLtgT7ANTlPrHUhkyiamZlZHTXJhCjNshwC7JMeu58JXESaVYmID8kkJtmzQ0+S2USdvak6X99TyCRUzwPPATdFxIt1CO9BMglLhb+SScQeTWMPA5C0oaQH0piLgdOBh1Pcd2YlcqT6HYChwKCIeC31e3Uq+wGwMCLer0OcZmZmlij/Ck7TJqk1MAPYPiIafeOMpHuAcyLi9UYabzDweUTcXKhei46do+NxVzVGSCstf9q9mVlpkTQ5InrVVK9JzhAVIqkf8ApwbTGSoWQImc3VjeUz4NZGHM/MzGyVssq9b01EPAZsUuQYXgVebcTxCm0cNzMzsxqscjNEZmZmZnXlhMjMzMxKnhMiMzMzK3mr3B4iq173Tm0p91NWZmZmVXiGyMzMzEqeEyIzMzMreU6IzMzMrOR5D1EJmfHuAsqG3F/sMFZJfgdsM7OmzTNEZmZmVvKcEJmZmVnJc0JkZmZmJc8JkZnVu3vvvZc+ffqw/vrr06pVKzbddFMOPvhgHnrooaV1Ro0ahSTmzJnTqLGdccYZ/L//9/8qnTv//PPZd999ad++PZIYNWpUte1HjBjBVlttRYsWLdhyyy0ZNmxYpfLPP/+cP/zhD+y66660b9+edu3aseuuu3LvvfcWjOuzzz5jgw02QBKPPfbYcl3bF198wdlnn03fvn1Ze+21kcT48ePz1v3++++57LLLKCsro2XLlvTo0YO77767Up3333+f8847j169etG2bVvWW2899t57byZMmFAwjjfffJPWrVsjidmzZ1cq69+/P6eddtpyXZ9ZQ3JCZGb16pprruGQQw6hc+fO3Hzzzdx///387ne/A+CJJ55YWu8nP/kJEydOpGPHjo0W2xtvvMGNN97I73//+0rnr732WhYuXMiBBx5YsP2IESM46aSTOPTQQ3nooYc4/PDDOfXUU7nhhhuW1vnvf//L9ddfz5577sno0aO544476NKlC4cccgjXXXddtX2fe+65SFqh6/vkk08YOXIkzZs3Z5999ilY94ILLuCiiy7i9NNP58EHH2TnnXfm8MMP54EHHlhaZ/Lkydxxxx3079+fu+66i1GjRtGyZUv69u3LfffdV23fp556Km3bts1bdtFFFzFixAhee+215btIswaiiCh2DNZIWnTsHB2Pu6rYYayS/JTZMptssgk77LAD99xzT5Wy77//ntVWK97/w8444wwmTZrECy+8UOl8RVyzZ8+mc+fO3HLLLQwcOLBSncWLF7PhhhtywAEHcOutty49f/zxxzN27Fjef/99Vl99db766isk0bp160rt9957b15//XX++9//Vonr2WefZd999+Xaa69l0KBBPProo/Tr16/O1xcRS5Oqxx57jH322Ycnn3ySvn37Vqr30UcfsfHGGzNkyBAuvvjiSjHOmzeP6dOnA5lZqzXXXJPmzZc9kLx48WK22WYbfvCDH+SdKbrtttsYPHgw5513HoMHD+b1119niy22qFSnd+/e9OrVi+uvv77O12hWV5ImR0SvmuqV9AyRpDJJL6XjXpKuqad+r5LUJx2PkfSqpJckjZS0ejVt/iRppqRZkq5RHf6rKOkxSevUR+xmK2r+/PlssMEGecuyk6HcJbOBAwciKe9X9rLPtGnTOOigg1hnnXVo1aoVu+22G08//XSNcX377beMHj2ao446qmBc1Zk4cSLz5s3jZz/7WaXzxxxzDJ988gnPPPMMAG3atKmSDAH06tWL9957r8r57777jpNOOokhQ4aw+eab1xhHIbX9tfHwww+zaNGiKtfys5/9jBkzZvDWW28B0K5du0rJEEDz5s3p2bMn7777bpV+P/30U371q19x5ZVX0q5du2rHP/LIIxkzZgwLFy6sVbxmjaGkE6JsEVEeEWeuaD+S1gV2joiK/zqNAbYCugOtgBPytNkV2A3YFugG7AjsWYdh/w6cugJhm9Wb3r17c+utt3LFFVfUaVnkggsuYOLEiZW+dtttN1q3bs0mm2wCwJQpU9h1112ZP38+I0aM4O6776Z9+/b069ePyZMnF+x/0qRJfPbZZ+yxxx7LdV0zZ84EoFu3bpXOb7PNNgC8/PLLBdtPmDCBrbbaqsr5P/3pTyxatIhzzjlnueJaHjNnzqRFixZVZm5qcy2LFi1i4sSJbL311lXKzjnnHLbaaiuOOeaYguP36dOHzz//nIkTJy5H9GYNo94TojTrMkvSiDTj8YikVqlsvKRe6biDpDnpeKCkeyWNk/SWpNMl/UrSi5ImpSSjov1Vkv6TZlx6S1pN0uuS1kt1VpM0W1KHnLj2lDQ1fb0oaa2c8r6S7kvHa0q6RdIMSdMlHZrO7ytpoqQpkv4pac0834LDgKU7RyPigUiA54GN8rQJoCWwBtACWB34MCe+tmmmacv0+nZJJ6biscCAgj8Ys0YybNgwtthiC8455xy23HJLOnTowIABA3jkkUcKtvvhD3/IzjvvvPTrmWeeYeLEiYwZM2bpzMlvfvMbNtlkE5544gkOO+wwfvzjH3PPPfew+eabc8kllxTsf9KkSUhi2223Xa7rmj9/PgDrrFN5MnbdddetVJ7P8OHDmTRpEuedd16l87Nnz+aPf/wj1113HS1atFiuuJbH/PnzadeuXZUZpdpcy0UXXcTcuXM599xzK51/5pln+Nvf/larZbAePXqw2mqrMWnSpOWI3qxhNNQMUWfguojYBvgMOLQWbboBRwG9gUuBryNiO2AicGxWvTYRsSuZGZGREfE9MBo4OpX3A6ZFxMc5/Z8NnBYRPYE9gEJztRcACyKie0RsCzyREqzfAf0iYnugHPhVnra7AVX+q5qWyo4hK1mqEBETgSeB99PXwxExK6fOAuB0YJSkI4F1ImJEKvsUaCGpfZ5xfyGpXFL5kq8XFLhks/rRpUsXXnzxRZ566il++9vf0rNnT+655x72228//vjHP9aqj3HjxnHuuedy+eWXc/DBBwOwcOFCnnrqKQ4//HBWW201Fi9ezOLFi4kI+vXrV+OTT++99x5rr702a6yxxnJdV8V+y7pufB4/fjxnnnkmxxxzDEcffXSlslNOOYX+/fvXuAG6vmXvNco9X8htt93G0KFDueCCCyrNtC1atIiTTjqJwYMH07Vr1xrHX3311Wnbtm3eJUSzYmmoj+54KyKmpuPJQFkt2jwZEV8AX0haAIxL52eQWUqqcDtAREyQtLakdsBI4N/AVcDxwC15+n8W+D9JY4B/RcTcAr/Y+gFHVryIiE8lHQh0BZ5N7dYgk6zl6gjMy3P+emBCRFTZ7CBpC2Brls0ePSqpT9ayW0Ucj0o6HLgO6JHTzUfAhsAnOW2GA8Mhs6k679Wa1bNmzZrRp08f+vTpA2SSkf3335+LL76Y0047rcosS7Zp06Zx1FFHMWjQIM4+++yl5+fPn8+SJUu45JJLqp0NKrRp+5tvvlmhWZjs2ZPsJ+MqZlMqyrO98MILHHTQQfzoRz/i5ptvrlR255138uyzz1JeXs5nn30GwJdffgnAV199xYIFC6p9UmtFrbvuunz66adVEqNPP/202msZN24cAwcOZNCgQZU2YgNcddVVzJ8/nzPPPHPptXz99ddA5q0AvvjiC9Zaq9KkPK1atfIeIlupNFRC9G3W8RIye2cAFrNsVqplgTbfZ73+nspx5v5Rj4h4R9KHkn4E7MSy2aLsSkMl3Q/8GJgkqR/wTTXxK884Ah6NiJqWphaSc22Sfg+sB5xUTZtDgEkR8WWq/yCwM1ApIZK0GpnEaSGwLjA3q7glhWe9zIpmww035IQTTuCXv/wlr7/+Or17985b78MPP+Sggw5i5513rrL00q5dO1ZbbTVOO+00jj322LztC22Obt++/dI/+MujYn/NzJkzKyVEFfttcmdGZsyYwX777UfPnj25++67WX31ys9TvPzyyyxcuHBpv9kOPvhg2rZtuzS5qG/bbLMN3377LW+88UalfUTVXcvjjz/O4YcfziGHHMKNN95Ypb+XX36ZDz74gE6dOlUp23777enRowdTp06tdH7+/Pl06NChSn2zYmnsD3edA+xAZi/NYcvZxxHAk5J2J7OsVbEOdBOZpbO/R8SS3EaSfhgRM4AZknYhs9F5am695BEyy1NnpbbrAJOA6yRtERGzJbUGNoqI3F2js4AtgPGp7QnAfsDeaXkvn/8CJ0q6jEzitSeZ2a5cg1P/5wMjJe0SEd+lJ9I2IPP9NSuqd955h4033rjK+VdeeQWg2ifQvvnmG/r370+bNm345z//WeXppjZt2rDHHnswbdo0tt9++zo/vr/VVlvx3XffMXfuXDbaKN9WvsJ22WUXOnTowJgxYyo9Ej969GjWXXdddtttt6XnXn/9dfbZZx8233xz7rvvPlq1alWlv4EDB1Z5HH7q1KkMHjyYK6+8kp122qnOMdbW/vvvzxprrMGYMWMqvSfT6NGj6datG5ttttnScxMnTqR///7svffejB49Ou/3fciQIVXepuChhx7i8ssvZ/To0Wy55ZaVyj744AO++eabKufNiqmxE6IrgTslHQM8UVPlanwq6T/A2mSWxyqMJbNUlm+5DOAsSXuRmbF6GXiQzPJWPn8kk/y8lOpfHBH/kjQQuF1Sxbz774DchOh+MjNBN6XXw4C3gYlpavpfEfGHtLn85Ig4AbgL+BGZ5cEAHoqIcdmdSupC5gm13hHxhaQJafzfk0kyJ0XE4mqux6zRdOvWjb322otDDjmEzTbbjM8//5wHHniAYcOG8T//8z9LnxjLddZZZzFlyhRGjRq1NHmq0LVrV9Zee23+7//+jz59+rDffvsxaNAgOnbsyMcff8yUKVNYsmQJQ4cOrTauiuW7559/vkpC9NRTTzFv3jw++OADAMrLy1lzzcwzE4cdlvm/2+qrr84ll1zCqaeeSqdOnejXrx9PPPEEI0eO5Nprr126N+mjjz5in332YdGiRVx88cVVntjabrvtaNGiBWVlZZSVleWNtUePHuy+++5LX48fP5699tor7/sj5XrwwQf56quvmDFjxtJr+/jjj2nTpg0HHHAAAOuvvz6DBw/msssuY6211mL77bfnjjvu4IknnuDf//730r5eeeUVfvKTn9ChQwd+85vfVHmSb+eddwYyyWbuE3QVb6ew0047VXma7bnnngOW/UzMVgZN6o0ZJY0Hzo6I8jxlvYC/RMTyPVNbjyQ9AxwYEQ0z3111vKuBsRHxeKF6fmPGhuM3Zlxm2LBhPPDAA0ybNo0PP/yQZs2a0aVLFwYMGMBZZ521NHEYNWoUP//5z3nrrbcoKyujb9++PPXUU3n7zH5zwVmzZnHxxRfzxBNPsGDBAtZbbz223357Tj75ZH784x8XjG2nnXaia9eu3HJL5f83FRo793fkjTfeyJ///GfefvttNtlkEwYPHsyppy5714uK5KU6FdebT0Xb3DdmvP/++znwwAN58MEH2X///QteY1lZGW+//XaV85tuummlj0lZsmQJl112GSNGjOCDDz5gyy235MILL1yaAMKyn1F1Cv39qGib740ZTzzxRF588UXKy6v8Kjerd6rlGzOuEgmRpCHAKcDREfFMMWLLiWcnYGFETG+k8U6seOKsECdEDccJUdMwatQofvnLX/L+++/nffPEldX555/P2LFjmTFjxgp/vEexffPNN3Ts2JErr7ySQYMGFTscKwG1TYia1BszRkTffLNDETE0IjZdGZIhgIh4rrGSoTRejcmQmWXeVbpTp05N7iMjnnrqKc4///wmnwxBZoZt/fXX57jjjit2KGaVNPYeIjOzomnWrBkjR45kypQpV6y2UgAAEQ9JREFUxQ6lTp599tlih1BvWrRowahRo6psmjcrtia1ZGYrxktmDcdLZmZmK6faLpk5RS8h3Tu1pdx/uM3MzKpoUnuIzMzMzBqCEyIzMzMreU6IzMzMrOQ5ITIzM7OS503VJWTGuwsoG3J/scOwEucn8sxsZeQZIjMzMyt5TojMzMys5DkhMjMzs5LnhMjMzMxKnhMiMzMzK3lOiOpAUl9J96XjgZL+Wk29gyVdmI43lfS4pOmSxkvaqJo2l0p6R9KXyxHX6ZJ+Xtd2ZmZmluGEKIek+ngrgnOA69PxlcDfImJb4A/AZdW0GQf0Xs7xRgJnLmdbMzOzkteoCZGkMkmzJI2QNFPSI5JapbLxknql4w6S5qTjgZLulTRO0ltpNuRXkl6UNEnSulntr5L0H0kvSeotaTVJr0taL9VZTdJsSR1y4rpI0nBJjwB/k9RS0i2SZqRx9qrDNXYBvo2Ij9OprsDj6fhJoH++dhExKSLer6Hva7JmnvaTNEHSahHxNTBH0vImVGZmZiWtGDNEnYHrImIb4DPg0Fq06QYcRWYG5VLg64jYDpgIHJtVr01E7AqcCoyMiO+B0cDRqbwfMC0rWcm2A9A/Io4CTgOIiO7AAOBWSS1reX27AVOyXk9j2TUeAqwlqX0t+8o1BDgiJWjXAD9P1whQDuyR20DSLySVSypf8vWC5RzWzMxs1VaMhOitiJiajicDZbVo82REfBER84AFZJaXAGbktL8dICImAGtLakdmOakiaToeuKWaMcZGxMJ0vDvw99TXK8DbQJdaxAnQEZiX9fpsYE9JLwJ7Au8Ci2vZVyVpJuhE4FHgrxHxRlbxR8CGedoMj4heEdGrWeu2yzOsmZnZKq8YH93xbdbxEqBVOl7MsgQtdzYmu833Wa+/p/I1RE67iIh3JH0o6UfATiybLcr1Vdaxqg+/RguBpZlHRLwH/BRA0prAoRGxIlM13YFPqJr8tExjm5mZWR2tTJuq55BZtgI4bDn7OAJA0u7AgqzE4yYyS2d3RsSSWvQzgZQ4pT1BmwCv1jKGWcAWFS/SfqiK7/N5ZGasloukTYFfA9sBB0jaKau4C/DS8vZtZmZWylamhOhK4BRJ/wE61FS5Gp+m9sOAQVnnxwJrUv1yWa7rgWaSZgB3AAMj4tsa2lSYAGwnqWKWqS/wqqTXgB+Q2QMFgKSpWcd/kjQXaC1prqSLsjtN/d0MnJ1mnQYBN2XtbdoNeKyWMZqZmVkWReSuMjVNksaTSRbK85T1Av4SEVU2HTdQLFcD4yKiURIUSdsBv4qIYwrVa9Gxc3Q87qrGCMmsWv60ezNrTJImR0SvmuqtTDNEDULSEOBuMstVjeV/gdaNOF4H4IJGHM/MzGyVUoxN1Q0iIvpWc34oMLSRY/mQzDJdY433aGONZWZmtipa5WeIzMzMzGrihMjMzMxK3iqzZGY1696pLeXe0GpmZlaFZ4jMzMys5DkhMjMzs5LnhMjMzMxKnhMiMzMzK3neVF1CZry7gLIh9xc7DDMzs7yK+U72niEyMzOzkueEyMzMzEqeEyIzMzMreU6IzMzMrOQ5ITIzM7OSt9IkRJKWSJoq6SVJ/5TUukDdMklHZb0eKOmvDRTXwZIuzBpnXopzqqQT8sQ/VdLYrPOnS5otKSR1KDDOJpIekTRL0suSyuoQ4z8kdV6+KzQzM7MGTYgk1eWx/oUR0TMiugGLgJML1C0DjipQXp/OAa7Pen1HirNnRNyUdX5h1vmDss4/C/QD3q5hnL8BV0TE1kBv4KM6xHhDitPMzMyWQ8GEKM3EvCLpVknTJd1VMXMjaU7FjIekXpLGp+OLJA2X9AjwN0nNJF0h6YXUx0m1iOtpYAtJl0j6ZVY8l0o6ExgK7JFmYwan4g0lPSTpdUl/ymozQNKMNPN0edb5L1N/0yRNkvSDPNffBfg2Ij6uRcx5RcSLETGnUB1JXYHmEfFoavNlRHydU6d5+h72Ta8vk3RpKn4a6FfHBNTMzMyS2swQbQkMj4htgc+BU2vRZgegf0QcBQwCFkTEjsCOwImSNquuYfqjfgAwA7gZOC6dXw04EhgDDAGeTrMxf0lNewJHAN2BIyRtLGlD4HLgR6l8R0kHp/ptgEkR0QOYAJyYJ5zdgCk55w7NSg43zjrfUlJ5Sq4Opm66AJ9J+pekF1MC2Sy7QkQsBgYCN0jaB9gfuDiVfQ/MBnrkdizpFymu8iVfL6hjWGZmZqWhNgnROxHxbDoeDexeizZjI2JhOt4XOFbSVOA5oD2Qb79Lq1SnHPgvcHOaWflE0napnxcj4pNqxnw8IhZExDfAy8CmZBKw8RExLyUUY4A+qf4i4L50PJnMMlyujsC8rNfjgLKUHD4G3JpVtklE9CKzlHeVpB9WE2c+zYE9gLNTzJuTSX4qiYiZwN9THMdHxKKs4o+ADfO0GR4RvSKiV7PWbesQkpmZWemozRJLVPN6McsSqpY5db7KOhZwRkQ8XMM4CyOiZ57zN5FJDjYARhZo/23W8RIy16YC9b+LiMipXyUmYGkWkZOMjSAz+1RR9l769820fLgd8EaB8bPNJZPsvQkg6V5gZzIzZLm6A58BuUt8LVO8ZmZmVke1mSHaRNIu6XgA8Ew6nkNmaQzg0ALtHwZOkbQ6ZPblSGpThxjvIbM8tGPqC+ALYK1atH0O2FNSh7QENQB4qg5jzwK2qHghqWNW2UGpHEnrSGqRjjuQWWp7uQ7jvACsI2m99PpH+dpL+imZGbY+wDWS2mUVdwFm1mFMMzMzS2qTEM0CjpM0HViXzBNNkNm/crWkp8nMsFTnJjJ/3KdIegm4kTp8qGxaFnoSuDMiKsaZDixOG6IHF2j7PnBeaj8NmBIR/67t2GT2Fm0nqWKm6UxJMyVNA85k2bLW1kB5Ov8kMDQiXgaQdKakucBGwHRJN6XzvSqO03WdDTwuaQaZma0R2YGkRGsoMCgiXgP+Clydyn5AZobt/Tpcm5mZmSVatmqUpzDzXjj3pUfhiyJtpp4CHB4Rrxdh/KuBcRHxWGOPXVspKfw8IvItsS3VomPn6HjcVY0UlZmZWd00xKfdS5qc9vgWtNK8MWM+6XH02WQ2TDd6MpT8L1Dtm0SuJD6j8gZvMzMzq4OCS1fpKa+izQ6lZafNizV+iuFDYGyNFYsoIm4pdgxmZmZN2Uo9Q2RmZmbWGJwQmZmZWcnzRz2UkO6d2lLeABvWzMzMmjrPEJmZmVnJc0JkZmZmJc8JkZmZmZU8J0RmZmZW8pwQmZmZWckr+NEdtmqRNA94ux66agssqId+VqSvurarbf2a6i1veQfg41qMvzKpz59zY41TKvcTNL17yvfTitf3/bRMXb7Hm0bEejXWigh/+atOX8DwYvdV13a1rV9TveUtB8qL/XMr5s+5scYplfsplTWpe8r304rX9/3UMD/nii8vmdnyGLcS9FXXdrWtX1O9FS1vShrrWnw/rfg4TYHvpxWv7/tpmXq/Fi+ZmTUCSeVRi09bNqst31NWn3w/eVO1WWMZXuwAbJXje8rqU8nfT54hMjMzs5LnGSIzMzMreU6IzMzMrOQ5ITIzM7OS54TIrMgkHSxphKR/S9q32PFY0yZpc0k3S7qr2LFY0ySpjaRb0++lo4sdT2NxQmS2AiSNlPSRpJdyzu8v6VVJsyUNKdRHRNwbEScCA4EjGjBcW8nV0/30ZkQMathIramp4731U+Cu9HvpoEYPtkicEJmtmFHA/tknJDUDrgMOALoCAyR1ldRd0n05X+tnNf1damelaxT1dz+ZZRtFLe8tYCPgnVRtSSPGWFTNix2AWVMWERMkleWc7g3Mjog3AST9A+gfEZcBB+b2IUnAUODBiJjSsBHbyqw+7iezfOpybwFzySRFUymhiZOSuVCzRtSJZf+7gswvl04F6p8B9AMOk3RyQwZmTVKd7idJ7SUNA7aTdF5DB2dNWnX31r+AQyXdwKr1cR8FeYbIrP4pz7lq3wE1Iq4Brmm4cKyJq+v99AngxNpqI++9FRFfAT9v7GCKzTNEZvVvLrBx1uuNgPeKFIs1fb6frKH43srihMis/r0AdJa0maQ1gCOBsUWOyZou30/WUHxvZXFCZLYCJN0OTAS2lDRX0qCIWAycDjwMzALujIiZxYzTmgbfT9ZQfG/VzB/uamZmZiXPM0RmZmZW8pwQmZmZWclzQmRmZmYlzwmRmZmZlTwnRGZmZlbynBCZmZlZyXNCZGYNStISSVMlvSRpnKR2K9BXX0m7Zr0+WdKx9RNprcbfQ9LMdD2talH/pvTp4cszVpmkl+rY5j9ZbY9annEL9H1+vrHMVhV+HyIza1CSvoyINdPxrcBrEXHpcvZ1EfBlRFxZjyHWZfxhwHMRcUsjjFUG3BcR3ZajbV/g7Ig4sA5tmkXEkgLlS3+OZqsizxCZWWOaSPqk9jTbc19FgaS/ShqYjudIuljSFEkzJG2VEoSTgcFphmYPSRdJOju1GS/pL5ImSJolaUdJ/5L0uqQ/Zo3zM0nPpz5ulNQsN0hJe0t6MY09UlILSScA/wNcKGlMTv02ku6XNC3NhB2RFVOvdPylpEtTnUmSfpDO/zC9fkHSHyR9mSeeZpKuSHWmSzop3zc3q+1QYI90jYOra59+Bk9Kug2Ykc7dK2lymgn7RTo3FGiV+huTPZYyrkjXPSPr2vum679L0iuSxkjK92GiZisFJ0Rm1ihS4rE3tf+spI8jYnvgBjKzHXOAYcBfIqJnRDydp82iiOiT6v0bOA3oBgyU1F7S1sARwG4R0RNYAhydE2dLYBRwRER0B5oDp0TETSn230REpTbA/sB7EdEjzeg8lCe2NsCkiOgBTABOTOevBq6OiB2p/oM1BwELUp0dgRMlbVZNXYAhwNPp+/SXGtr3Bn4bERVLe8dHxA5AL+BMSe0jYgiwMPWXe+0/BXoCPYB+wBWSOqay7YCzgK7A5sBuBWI2KyonRGbW0FpJmgp8AqwLPFrLdv9K/04GymrZpiLZmgHMjIj3I+Jb4E0yn+q9N7AD8EKKaW8yf6izbQm8FRGvpde3An1qGHcG0E/S5ZL2iIgFeeosAipmxLKvaRfgn+n4tmr63xc4NsX8HNAe6FxDTLVt/3xEvJVV90xJ04BJZL5nNY2zO3B7RCyJiA+Bp8gkXRV9z42I74Gp1P7naNbomhc7ADNb5S2MiJ6S2pJJCE4DrgEWU/k/ZS1z2n2b/l1C7X9XVbT5Puu44nVzQMCtEXFegT7qvKwTEa9J2gH4MXCZpEci4g851b6LZZs263JNFTGdEREP1zW2Qu3TXqOvcl73A3aJiK8ljafqzyVf39XJ/hnU9ZrNGpVniMysUaRZkzOBsyWtDrwNdE37c9qSma2pyRfAWisQxuPAYZLWB5C0rqRNc+q8ApRJ2iK9PobMrEe1JG0IfB0Ro4Erge3rENMk4NB0fGQ1dR4GTknfNyR1kdSmQJ+536fatm8LfJqSoa2AnbPKvqton2MCcETap7Qemdm05wvEZrZSckJkZo0mIl4EpgFHRsQ7wJ3AdGAM8GItuhgHHFKxqXo5xn8Z+B3wiKTpZJbvOubU+Qb4OfBPSTPIzC4Nq6Hr7sDzaUnqt8Afa6if7SzgV5KeT7HkW267CXgZmKLMo/g3Uni2ZTqwOG3gHlyH9g8BzdP35hIyyVqF4cD03A3lwD1pvGnAE8A5EfFBoQs2Wxn5sXszsyKS1JrMsmJIOhIYEBH9ix2XWanxeq6ZWXHtAPw1PZL+GXB8keMxK0meITIzM7OS5z1EZmZmVvKcEJmZmVnJc0JkZmZmJc8JkZmZmZU8J0RmZmZW8pwQmZmZWcn7/4EI+iOdz5+1AAAAAElFTkSuQmCC)
%% Cell type:markdown id: tags:
All runtimes are plotted logarithmically. Next number next to the labels shows how much slower the version is than the fastest one. For small arrays Cython produces faster code than *pystencils*. The larger the arrays, the better pystencils gets.
doc/notebooks/demo_wave_equation.ipynb 0 → 100644
View file @ 1e02cdc7
This source diff could not be displayed because it is too large. You can view the blob instead.
doc/sphinx/api.rst 0 → 100644
View file @ 1e02cdc7
API Reference
=============
.. toctree::
:maxdepth: 3
kernel_compile_and_call.rst
simplifications.rst
datahandling.rst
configuration.rst
field.rst
finite_differences.rst
plot.rst
ast.rst
doc/sphinx/ast.rst 0 → 100644
View file @ 1e02cdc7
*********************************************
For developers: AST Nodes and Transformations
*********************************************
AST Nodes
=========
.. automodule:: pystencils.astnodes
:members:
Transformations
===============
.. automodule:: pystencils.transformations
:members:
doc/sphinx/configuration.rst 0 → 100644
View file @ 1e02cdc7
*************
Configuration
*************
.. automodule:: pystencils.cpu.cpujit
\ No newline at end of file
doc/sphinx/datahandling.rst 0 → 100644
View file @ 1e02cdc7
************
DataHandling
************
.. autoclass:: pystencils.datahandling.DataHandling
:members:
\ No newline at end of file
doc/sphinx/field.rst 0 → 100644
View file @ 1e02cdc7
*****
Field
*****
.. automodule:: pystencils.field
:members:
\ No newline at end of file
doc/sphinx/finite_differences.rst 0 → 100644
View file @ 1e02cdc7
******************
Finite Differences
******************
.. automodule:: pystencils.fd
:members:
doc/sphinx/kernel_compile_and_call.rst 0 → 100644
View file @ 1e02cdc7
*****************************************
Creating and calling kernels from Python
*****************************************
Creating kernels
----------------
.. autofunction:: pystencils.create_kernel
.. autofunction:: pystencils.create_indexed_kernel
.. autofunction:: pystencils.create_staggered_kernel
Code printing
-------------
.. autofunction:: pystencils.show_code
GPU Indexing
-------------
.. autoclass:: pystencils.gpucuda.AbstractIndexing
:members:
.. autoclass:: pystencils.gpucuda.BlockIndexing
:members:
.. autoclass:: pystencils.gpucuda.LineIndexing
:members:
doc/sphinx/plot.rst 0 → 100644
View file @ 1e02cdc7
**********************
Plotting and Animation
**********************
.. automodule:: pystencils.plot2d
:members:
doc/sphinx/simplifications.rst 0 → 100644
View file @ 1e02cdc7
***************************************
Assignment Collection & Simplifications
***************************************
AssignmentCollection
====================
.. autoclass:: pystencils.AssignmentCollection
:members:
Simplifications
===============
.. automodule:: pystencils.simp
:members:
doc/sphinx/tutorials.rst 0 → 100644
View file @ 1e02cdc7
Tutorials
=========
These tutorials are a good place to start if you are new to pystencils.
All tutorials and demos listed here are based on Jupyter notebooks that you can find in the pystencils repository.
It is a good idea to download them and run them directly to be able to play around with the code.
.. toctree::
:maxdepth: 1
/notebooks/01_tutorial_getting_started.ipynb
/notebooks/02_tutorial_basic_kernels.ipynb
/notebooks/03_tutorial_advection_diffusion.ipynb
/notebooks/04_tutorial_phasefield_spinodal_decomposition.ipynb
/notebooks/05_tutorial_phasefield_dentritic_growth.ipynb
/notebooks/06_tutorial_datahandling.ipynb
/notebooks/demo_assignment_collection.ipynb
/notebooks/demo_benchmark.ipynb
/notebooks/demo_wave_equation.ipynb
pacxx/benchmark.py deleted 100644 → 0
View file @ b91e6990
import os
from time import perf_counter
import subprocess
from tempfile import TemporaryDirectory
from pystencils import create_data_handling
from pystencils.backends.cbackend import CBackend
from jinja2 import Environment, FileSystemLoader
from pystencils.backends.cbackend import generate_c
script_path = os.path.dirname(os.path.realpath(__file__))
PAXX_ROOT = '/local/bauer/code/pacxx/install'
DEFAULT_PAXX_COMPILE_OPTIONS = ('-Ofast', '-march=native')
def generate_benchmark_code(target_file, kernel_ast, target):
assert target in ('cpu', 'gpu')
assert hasattr(kernel_ast, 'indexing'), "AST has to be a CUDA kernel in order to create a PACXX kernel from it"
backend = CBackend()
function_body = kernel_ast.body
f_sizes = {f.shape[-1] for f in kernel_ast.fields_accessed}
assert len(f_sizes) == 1
env = Environment(loader=FileSystemLoader(script_path))
result = env.get_template("benchmark_template.cpp").render(f_size=f_sizes.pop(),
code=backend(function_body),
target=target)
with open(target_file, 'w') as f:
f.write(result)
def pacxx_compile(source, executable, options=DEFAULT_PAXX_COMPILE_OPTIONS):
command = ['pacxx++', *options, source, '-o', executable, ]
env = os.environ.copy()
env['PATH'] = "{}:{}".format(env.get('PATH', ''), os.path.join(PAXX_ROOT, 'bin'))
env['LD_LIBRARY_PATH'] = "{}:{}".format(env.get('LD_LIBRARY_PATH', ''), os.path.join(PAXX_ROOT, 'lib'))
try:
subprocess.check_output(command, env=env, stderr=subprocess.STDOUT)
except subprocess.CalledProcessError as e:
print(" ".join(command))
print(e.output.decode('utf8'))
raise e
def run_paxx_benchmark(executable, domain_size, iterations):
assert len(domain_size) == 3
arguments = [executable, *domain_size, iterations]
arguments = [str(e) for e in arguments]
output = subprocess.check_output(arguments)
return float(output) / iterations
def paxx_benchmark(ast, domain_size, iterations, target='cpu', compile_options=DEFAULT_PAXX_COMPILE_OPTIONS):
"""Generates, compiles and runs the kernel with PAXX
Args:
ast: pystencils AST object (has to be generated for CUDA, even when run on CPU with pacxx)
domain_size: x, y, z extent of spatial domain
iterations: number of outer iterations
target: either 'cpu' or 'gpu' to specify where pacxx should run the kernel
compile_options: compile options for pacxx
Returns:
seconds for one outer iteration
"""
with TemporaryDirectory() as base_dir:
code = os.path.join(base_dir, 'code.cpp')
executable = os.path.join(base_dir, 'bench')
generate_benchmark_code(code, ast, target)
pacxx_compile(code, executable, compile_options)
time_per_iteration = run_paxx_benchmark(executable, domain_size, iterations)
return time_per_iteration
def lbm_performance_compare(domain_size, iterations, **lb_params):
"""Runs benchmark with pacxx and with normal pystencils backends.
Args:
domain_size: 3-tuple with size of spatial domain
iterations: number of outer iterations
**lb_params: parameters passed to lbmpy to choose lattice Boltzmann algorithm & optimization options
Returns:
dictionary with measurements of time per iteration for different backends
"""
import pycuda.driver as drv
from lbmpy.creationfunctions import create_lb_ast
if 'optimization' not in lb_params:
lb_params['optimization'] = {}
lb_params['optimization']['target'] = 'cpu'
cpu_ast = create_lb_ast(**lb_params)
lb_params['optimization']['target'] = 'gpu'
gpu_ast = create_lb_ast(**lb_params)
# print kernel code of CPU or GPU version - just for comparison, files are not used
with open("pystencils_cpu_code.c", 'w') as f:
print(generate_c(cpu_ast), file=f)
with open("pystencils_gpu_code.cu", 'w') as f:
print(generate_c(gpu_ast), file=f)
cpu_kernel = cpu_ast.compile()
gpu_kernel = gpu_ast.compile()
f_sizes = {f.shape[-1] for f in cpu_ast.fields_accessed}
assert len(f_sizes) == 1
f_size = f_sizes.pop()
dh = create_data_handling(domain_size, default_target='gpu', default_layout='fzyx')
dh.add_array('src', values_per_cell=f_size)
dh.add_array('dst', values_per_cell=f_size)
dh.fill('src', 0)
dh.fill('dst', 0)
# to keep it simple we run outer loop directly from Python
# make domain size large enough, otherwise we measure the python call overhead
def run_benchmark(kernel):
dh.all_to_gpu()
for i in range(10): # warmup
dh.run_kernel(kernel)
drv.Context.synchronize()
start = perf_counter()
for i in range(iterations):
dh.run_kernel(kernel)
drv.Context.synchronize()
return (perf_counter() - start) / iterations
return {
'pystencils_cpu': run_benchmark(cpu_kernel),
'pystencils_gpu': run_benchmark(gpu_kernel),
'pacxx_cpu': paxx_benchmark(gpu_ast, domain_size, iterations, target='cpu'),
'pacxx_gpu': paxx_benchmark(gpu_ast, domain_size, iterations, target='gpu'),
}
if __name__ == '__main__':
no_opt = {
'openmp': 8, # number of threads - pacxx uses also HT cores
'split': False,
'vectorization': False,
'gpu_indexing_params': {'block_size': (64, 8, 1)},
}
only_vectorization = {
'openmp': 4,
'split': False,
'gpu_indexing_params': {'block_size': (64, 8, 1)},
'vectorization': {'instruction_set': 'avx',
'assume_inner_stride_one': True,
'nontemporal': False},
}
best = {
'openmp': 4,
'split': True,
'gpu_indexing_params': {'block_size': (64, 8, 1)},
'vectorization': {'instruction_set': 'avx',
'assume_inner_stride_one': True,
'nontemporal': True}
}
res = lbm_performance_compare(stencil='D3Q19', relaxation_rate=1.8, compressible=False,
domain_size=(512, 128, 32), iterations=500,
optimization=only_vectorization)
cpu_speedup = ((res['pacxx_cpu'] / res['pystencils_cpu']) - 1) * 100
gpu_speedup = ((res['pacxx_gpu'] / res['pystencils_gpu']) - 1) * 100
print("Time for one kernel call [s]")
for config_name, time in res.items():
print(" {0: <16}: {1}".format(config_name, time))
print("CPU {:.02f}% GPU {:.02f}%".format(cpu_speedup, gpu_speedup))
pacxx/benchmark_template.cpp deleted 100644 → 0
View file @ b91e6990
#include <PACXX.h>
#include <vector>
#include <sstream>
#include <iostream>
#include <chrono>
using namespace pacxx::v2;
size_t division_round_up(size_t a, size_t b)
{
if( a % b == 0)
return a / b;
else
return (a / b) + 1;
}
int main(int argc, char** argv)
{
{% if target == 'cpu' %}
Executor::Create<NativeRuntime>(0);
{% elif target == 'gpu' %}
Executor::Create<CUDARuntime>(0);
{% endif %}
if( argc != 5 ) {
std::cout << "Usage: ./benchmark xSize ySize zSize iterations" << std::endl;
return 1;
}
Dimension3 domainSize;
int64_t iterations;
auto &exec = Executor::get(0);
std::stringstream( argv[1] ) >> domainSize.x;
std::stringstream( argv[2] ) >> domainSize.y;
std::stringstream( argv[3] ) >> domainSize.z;
std::stringstream( argv[4] ) >> iterations;
// add ghost layers to be comparable to pystencils native backend
domainSize.x += 2;
domainSize.y += 2;
domainSize.z += 2;
int64_t totalSize = domainSize.x * domainSize.y * domainSize.z * {{f_size}};
std::vector<double> src( totalSize, 0.0 );
std::vector<double> dst( totalSize, 0.0 );
auto & dsrc = exec.allocate<double>(src.size());
auto & ddst = exec.allocate<double>(dst.size());
dsrc.upload(src.data(), src.size());
ddst.upload(dst.data(), dst.size());
double * _data_src = dsrc.get();
double * _data_dst = ddst.get();
const int64_t _size_src_0 = domainSize.x;
const int64_t _size_src_1 = domainSize.y;
const int64_t _size_src_2 = domainSize.z;
// fzyx layout
const int64_t _stride_src_0 = 1;
const int64_t _stride_src_1 = domainSize.x;
const int64_t _stride_src_2 = domainSize.x * domainSize.y;
const int64_t _stride_src_3 = domainSize.x * domainSize.y * domainSize.z;
auto pacxxKernel = [=]( range & config ) {
struct Vec3D {int x; int y; int z; };
const Vec3D blockDim = { config.get_block_size(0), config.get_block_size(1), config.get_block_size(2) };
const Vec3D blockIdx = { config.get_block(0), config.get_block(1), config.get_block(2) };
const Vec3D threadIdx = { config.get_local(0), config.get_local(1), config.get_local(2) };
{{ code|indent(8) }}
};
size_t blockSize[] = {64, 8, 1};