From ba0866592a881dac95232e42dfa757ab6c9cc602 Mon Sep 17 00:00:00 2001 From: Felix Koehler Date: Wed, 10 Apr 2024 10:51:38 +0200 Subject: [PATCH] Add a notebook on performance hints --- examples/performance_hints.ipynb | 593 +++++++++++++++++++++++++++++++ 1 file changed, 593 insertions(+) create mode 100644 examples/performance_hints.ipynb diff --git a/examples/performance_hints.ipynb b/examples/performance_hints.ipynb new file mode 100644 index 0000000..88e3e1c --- /dev/null +++ b/examples/performance_hints.ipynb @@ -0,0 +1,593 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Performance Hints\n", + "\n", + "How to run `Exponax` even faster than it already is 😉.\n", + "\n", + "This is beyond some general insights:\n", + "\n", + "* Whenever the `exponax.rollout` or `exponax.repeat` function transformations\n", + " are used, they internally perform a `jax.jit` over the timestepper. Hence,\n", + " there is no need to wrap the resulting function in a `jax.jit` again.\n", + " However, when using the timesteppers directly, it can be advantageous to\n", + " Just-In-Time compile them.\n", + "* The number of total degrees of freedom scale exponentially with the number of\n", + " dimensions; so does the cost of the spatial FFT and hence the cost of\n", + " simulation. As a good guideline on a modern GPU:\n", + " * 1d: Highest still nice `num_points` is between 10'000 to 100'000. For most\n", + " problems, 50-500 points are likely sufficient.\n", + " * 2d: Highest still nice `num_points` is around 500 (-> 25k total DoF per\n", + " channel). For most problems, 50-256 points are likely sufficient.\n", + " * 3d: Highest still nice `num_points` is around 48 (-> 110k total DoF per\n", + " channel). In general, 3d sims will be tough.\n", + "* The produced trajectory array is as large as the number of time steps\n", + " performed. Hence, if the underlying discretization already has a lot of\n", + " total DoF, the trajectory array can become quite large. If you are only\n", + " interested in every n-th step, consider wrapping the time stepper in a\n", + " `RepeatedStepper`\n", + "* Some usages of `jax.vmap` only work efficiently on GPUs & TPUs, on CPUs JAX\n", + " resorts to sequential looping." + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [], + "source": [ + "import jax\n", + "import jax.numpy as jnp\n", + "import exponax as ex\n", + "import equinox as eqx" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Temporal Rollout in Fourier space\n", + "\n", + "The methods of `Exponax` advance a state to the next time step in Fourier space.\n", + "If the stepper is called with a state in physical space, it is first transformed\n", + "to Fourier space, then advanced, and finally transformed back to physical space.\n", + "This is done for each time step. We can also integrate it directly in Fourier\n", + "space and then backtransform the entire trajectory.\n", + "\n", + "This especially saves compute for lower orders of EDTRK integrators (in the\n", + "greatest sense for linear PDEs) that perform fewer FFTs per time step." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAksAAAG2CAYAAABvWcJYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/SrBM8AAAACXBIWXMAAA9hAAAPYQGoP6dpAACCKElEQVR4nO3deXxU1fk/8M9smez7MgkkJGELyA4Sg7gSFrEK6s+KYlGq+JVKq2JdsC5FrYpfta2Wamtdvy6orVpURCIIqISwhjWELQvZ9z2ZzHJ/f8zcSSJZJsm9c2cmn/frlZdm5s69Zw6TmWfOec5zVIIgCCAiIiKibqmVbgARERGRO2OwRERERNQLBktEREREvWCwRERERNQLBktEREREvWCwRERERNQLBktEREREvWCwRERERNQLBktEREREvWCwRERERNQLjwqWdu7ciWuuuQZxcXFQqVT44osv+nzM9u3bMW3aNOj1eowaNQrvvPPOecesX78eiYmJ8PX1RWpqKvbs2SN944mIiMgjeVSw1NzcjMmTJ2P9+vVOHZ+Xl4err74aV1xxBbKzs3HffffhzjvvxLfffus45uOPP8bq1avx5JNP4sCBA5g8eTLmz5+PiooKuZ4GEREReRCVp26kq1Kp8Pnnn2Px4sU9HvPwww/j66+/xtGjRx23LVmyBHV1ddi8eTMAIDU1FRdeeCH+9re/AQCsVivi4+Px29/+Fo888oisz4GIiIjcn1bpBsgpMzMT6enpXW6bP38+7rvvPgBAe3s79u/fjzVr1jjuV6vVSE9PR2ZmZo/nNRqNMBqNjt+tVitqamoQEREBlUol7ZMgIiIiWQiCgMbGRsTFxUGt7nmyzauDpbKyMsTExHS5LSYmBg0NDWhtbUVtbS0sFku3x5w4caLH8z733HNYu3atLG0mIiIi1zp37hyGDx/e4/1eHSzJZc2aNVi9erXj9/r6eiQkJCAvLw9BQUGSXcdkMuH777/HFVdcAZ1OJ9l5qSul+lkQBDz7TS7+faAEADDWEIgXrr8ACeEBLmuDq/E17RrsZ9dgP7uGnP3c2NiIpKSkPj+7vTpYMhgMKC8v73JbeXk5goOD4efnB41GA41G0+0xBoOhx/Pq9Xro9frzbg8PD0dwcLA0jYftBeLv74+IiAj+IcpIqX5+88c8fHasDhpffwT76nCq1oRb/+84nr1+IhZNGeaydrgSX9OuwX52Dfaza8jZz+L5+kqh8ajVcP2VlpaGrVu3drktIyMDaWlpAAAfHx9Mnz69yzFWqxVbt251HEMkh++Ol+OZr48DAP6wcBy+ve9SzEwKR3O7BfduyMbD/z6MNpNF4VYSERHgYcFSU1MTsrOzkZ2dDcBWGiA7OxuFhYUAbNNjy5Ytcxx/99134+zZs3jooYdw4sQJ/P3vf8cnn3yC+++/33HM6tWr8cYbb+Ddd99FTk4OVq5ciebmZixfvtylz42GjqPF9fjdhoMQBOCW1ATcMTsJhhBffHhnKn43ZzRUKuDjfefw+o4zSjeViIjgYdNw+/btwxVXXOH4Xcwbuu222/DOO++gtLTUETgBQFJSEr7++mvcf//9+Otf/4rhw4fjX//6F+bPn+845qabbkJlZSWeeOIJlJWVYcqUKdi8efN5Sd9EUiirb8Md7+5FS7sFl4yOxNprL3AM/2o1aqyeOwbDQn3x8H+OYMOec1h1xShoNR71nYaIyOt4VLB0+eWXo7eyUN1V57788stx8ODBXs+7atUqrFq1arDNI+rTi1tyUd5gxOjoQPztlmnQdRMILZ46DM9/cwJlDW3YeaoSV6YwcCciUhK/shK5iNFswbfHygAAzyyegBC/7hMV9VoNbphmW8L60Z5zLmsfERF1j8ESkYvsOl2NxjYzooP0uDAxvNdjb7owHgCw7UQFKhraXNE8IiLqAYMlIhfZdKQUALBgggFqde/LVEfHBGH6iDBYrAI+3V/kiuYREVEPGCwRuYDJYsWW47Z6Xgsnxjr1mCX20aVP9p2D1eqRWzgSEXkFBktELpB5phr1rSZEBvr0OQUnunpSLIL0WhRUt2D32WqZW0hERD1hsETkAuIU3PwLDND0MQUn8vfR4topcQCADXuZ6E1EpBQGS0QyM1usjlVwzk7BiZZcmAAA2Hy0DLXN7ZK3jYiI+sZgiUhmWXk1qG0xITzAB6lJzk3BiSYOD8EFccFot1jx+cFimVpIRES9YbBEJLOOKbiYAVXjFhO9N+wt7ONIIiKSA4MlIhlZrIJjCu6qCf2bghNdO2UYNGoVTpY3oai2RcrmERGRExgsEcloT14NqpraEeKnQ9rIiAGdI8RPhwnDQhznIyIi12KwRCSjb47apuDmjY/pdh84Z11kz3XKOstgiYjI1RgsEcnEahWw+ah9FdykgU3BiVKTbcHSnnwGS0RErsZgiUgmJ8oaUdFoRICPBhePjBzUuaaPCIdKBeRVNXOvOCIiF2OwRCST7HN1AIApCaHw0Q7uTy3ET4fxscEAgN3MWyIicikGS0QyOWQPliYPD5XkfKlJtgTxPXnc+oSIyJUYLBHJxDGyFB8qyflmMsmbiEgRDJaIZNBkNONkRSMA6YOlUxVNqG4ySnJOIiLqG4MlIhkcKaqHIABxIb6IDvaV5JzhAT4YExMIANjLVXFERC7DYIlIBoeK6gAAkyUaVRKJeUtZTPImInIZBktEMsgurAMg3RScSKy3xLwlIiLXYbBEJAO5RpbEvKWcsgbUt5gkPTcREXWPwRKRxMob2lBa3wa1Cpho39NNKtFBvkiODIAgAPsKOLpEROQKDJaIJCaWDBgTE4QAvVby8ztKCDBviYjIJRgsEUlM6vpKP+fIW2KwRETkEgyWiCQmV3K3SFwRd7S4Hk1GsyzXICKiDgyWiCRksQo4UlwPQPrkblFcqB+Gh/nBYhWwv6BWlmsQEVEHBktEEjpT2YQmoxn+PhqMiQmS7Tri6NI+FqckIpIdgyUiCYn5ShOGhUCjVsl2ncnxtlV2x0oaZLsGERHZMFgikpAYLE2VaQpOdEFcMADgWEm9rNchIiIGS0SSOmQPluTKVxKlGIKhUgHlDUZUcVNdIiJZMVgikkhruwUnyhoByLcSThSg1yIpIgAAcJxTcUREsmKwRCSRYyX1sFgFRAXpERviK/v1xjum4hgsERHJicESkUQ6F6NUqeRL7haNZ94SEZFLMFgikoijvtJwafeD68kFcbbrHC/lyBIRkZw8Llhav349EhMT4evri9TUVOzZs6fHYy+//HKoVKrzfq6++mrHMbfffvt59y9YsMAVT4W8zMnyJgC25GtXGB9ru05eVTOaWcmbiEg2HhUsffzxx1i9ejWefPJJHDhwAJMnT8b8+fNRUVHR7fGfffYZSktLHT9Hjx6FRqPBjTfe2OW4BQsWdDnuo48+csXTIS9isQo4U2kLluQsRtlZVJAe0UF6CAJwooyjS0REcvGoYOnll1/GihUrsHz5cowfPx6vv/46/P398dZbb3V7fHh4OAwGg+MnIyMD/v7+5wVLer2+y3FhYWGueDrkRQqqm9FutsJXp8bwMD+XXVest8QVcURE8tEq3QBntbe3Y//+/VizZo3jNrVajfT0dGRmZjp1jjfffBNLlixBQEBAl9u3b9+O6OhohIWF4corr8QzzzyDiIiIHs9jNBphNHbUtmlosH1QmUwmmEym/jytXonnkvKcdD4p+jmnpA4AMDIqABaLGRaLFC3rW4ohEN/nVuJIUZ1HvE74mnYN9rNrsJ9dQ85+dvacKkEQBMmvLoOSkhIMGzYMu3btQlpamuP2hx56CDt27EBWVlavj9+zZw9SU1ORlZWFmTNnOm7fsGED/P39kZSUhDNnzuDRRx9FYGAgMjMzodFouj3XH//4R6xdu/a82z/88EP4+/sP8BmSJ9tSpMLX5zS4MNKKW0dbXXbd7GoV3j6pQXyAgN9PclGERkTkJVpaWnDLLbegvr4ewcE955t6zMjSYL355puYOHFil0AJAJYsWeL4/4kTJ2LSpEkYOXIktm/fjjlz5nR7rjVr1mD16tWO3xsaGhAfH4958+b12tn9ZTKZkJGRgblz50Kn00l2XupKin7O+OQwcK4Ml00di4WXJkncwp5dUNOCt0/+iLI2NebOnwudxr1n1vmadg32s2uwn11Dzn4WZ4b64jHBUmRkJDQaDcrLy7vcXl5eDoPB0Otjm5ubsWHDBjz11FN9Xic5ORmRkZE4ffp0j8GSXq+HXq8/73adTifLH4xc56WuBtPPpyubAQApsSEu/bdKjgpGkF6LRqMZhXVGl63EGyy+pl2D/ewa7GfXkKOfnT2fe38N7cTHxwfTp0/H1q1bHbdZrVZs3bq1y7Rcdz799FMYjUbceuutfV6nqKgI1dXViI2NHXSbaWgwW6w4aw+WXLUSTqRWqzBOLE5ZzCRvIiI5eEywBACrV6/GG2+8gXfffRc5OTlYuXIlmpubsXz5cgDAsmXLuiSAi958800sXrz4vKTtpqYmPPjgg9i9ezfy8/OxdetWLFq0CKNGjcL8+fNd8pzI8xXWtKDdYoWfTuPSlXAisd4Stz0hIpKHx0zDAcBNN92EyspKPPHEEygrK8OUKVOwefNmxMTEAAAKCwuhVneN/3Jzc/Hjjz9iy5Yt551Po9Hg8OHDePfdd1FXV4e4uDjMmzcPTz/9dLfTbETdEYtRjooOhFot/zYnP+coH1DKbU+IiOTgUcESAKxatQqrVq3q9r7t27efd9vYsWPR04I/Pz8/fPvtt1I2j4agU+WNAIDR0YGKXN+x7UlJAwRBcMm+dEREQ4lHTcMRuaOTFbaRpdEuzlcSjYoOhE6jQkObGUW1rYq0gYjImzFYIhokcWRpTIwyI0s+WrUjsfxYCafiiIikxmCJaBA6r4QbHa3MyBLAbU+IiOTEYIloEAoUXgkn4oo4IiL5MFgiGoRTCq+EE10wzJbkzWCJiEh6DJaIBsGxEk6hfCXROPvIUllDG2qb2xVtCxGRt2GwRDQIjpVwCuYrAUCgXothobZpwFP2NhERkTQYLBENgtIr4ToT23DS3iYiIpIGgyWiAVJyT7juiG04xWCJiEhSDJaIBqjzSjhxCkxJYlFMcfsVInfQ2m7BK1tPIeN4udJNIRowj9vuhMhdiCM4Sq+EE4nTcKcqOLJE7qG+xYQ739uLvfm1AICHF6Tg7suSuSUPeRyOLBENkDiCo/RKONEo+950VU3tqOGKOFJYeUMbfvmPTOzNr4WP1vZRs27zCTzzdQ6s1u736yRyVwyWiAZIXHXmDvlKAODvo0V8uG06kEnepKSzlU24/u+7kFveiOggPf57z8X4w8JxAIA3f8zD6k+y0W62KtxKIucxWCIaIHdaCScaE80kb1LWsZJ63Ph6JorrWpEY4Y//rJyFcbHBWHFpMv5802Ro1Sp8kV2CO9/bhzaTRenmEjmFwRLRALjLnnA/xyRvUtqjnx9FdXM7JgwLxr9XzkJ8uL/jvuumDscbt82An06DnScr8f7uAgVbSuQ8BktEA1DoZivhRKy1RErKq2rGoXN10KhVeOv2CxEZqD/vmCvGRuMPV9um5D7MKoQgMH+J3B+DJaIByKuyjSolRga4xUo4kaPWEqt4kwI2ZpcAAC4eFYnoIN8ej1s8dRgCfDQ4W9WMzLPVrmoe0YAxWCIaADFYSo4MULglXY2MCoRKBdQ0t6Oqyah0c2gIEQQB/z1UDABYPCWu12MD9VosmjoMgG10icjdMVgiGoCOkSX/Po50LT8fDRLsOSKciiNXOlrcgLOVzdBr1Zh3gaHP42+ZmQAA+PZYGQN7cnsMlogGQAyWkiLdZyWcaLRjRRyn4sh1/pttG1VKHx+DQH3f9Y4nDAvB5PhQmCwCPt1XJHfziAaFwRLRAOQ7giX3GlkCmORNrmexCvjysC1fadHk3qfgOluaahtd+mhPIQtVkltjsETUT63tFpTUtwFwz5Gljg11ObJErpGVV43yBiOCfbW4bGyU04+7ZlIcgny1KKxpwY+nq2RsIdHgMFgi6qeCGtuoUoifDmH+OoVbcz5x+5WTFY1clk0uIa6CWzgxFnqtxunH+flocMO04QCY6E3ujcESUT/lVXaUDXDHDUFHRgVCrQLqWkyoZOIsycxotmDTkVIAwLV9rILrzi32qbiMnHKUN7RJ2jYiqTBYIuqnvGr3LBsg8tVpMCLC1jZOxZHcduRWoqHNjJhgPVKTIvr9+DExQbgwMQwWq4BP9p6ToYVEg8dgiaifHCNLEe4ZLAHA6GgmeZNr/PeQbQru2slx0AywQOvS1BEAgE/3F3HqmNwSgyWifsq3jywlRblvsDSGe8SRCzQZzfjueDkAYNGUYQM+z7wLYqDTqFBY04KC6hapmkckGQZLRP3kqLHkziNL9iTvUxxZIhllnqmG0WxFYoQ/LogLHvB5/H20mD4iDADww6lKqZpHJBkGS0T90NBmQlVTOwD3q97dWcfIElfEkXyyz9UCAC5MDB/0YodLRttKDuw8xRIC5H4YLBH1g1iMMjJQjyBf9ysbIEqOCoBGrUJDmxkVjVwRR/I4dK4eADAlIXTQ57rUHixlnqmGyWId9PmIpMRgiagf3HUD3Z/TazUYEcE94kg+VquAQ0V1AIDJw0MHfb4L4oIR5q9Dk9GM7HN1gz4fkZQYLBH1g7tuoNudMdFM8ib5nK1qRmObGb46NcYaggZ9PrVahdn20aUfTjJvidwLgyWifsh34w10f07cI+50BUeWSHqH7KM/E+JCoNNI81FyyehIAMxbIvfDYImoH/LceAPdnxsZLa6I48gSSc8xBRcfKtk5xWDpcFEd6lraJTsv0WAxWCJykiAIOOtBI0uj7dNwpyqauCKOJCfmFU2RMFiKDfHDqOhAWAVg15lqyc5LNFgeFyytX78eiYmJ8PX1RWpqKvbs2dPjse+88w5UKlWXH19f3y7HCIKAJ554ArGxsfDz80N6ejpOnTol99MgD1TT3I7GNjMAOJKn3VlyVABUKqC+taPcAZEU2kwW5JQ2AJA2WAI6Rpd+4FQcuRGPCpY+/vhjrF69Gk8++SQOHDiAyZMnY/78+aioqOjxMcHBwSgtLXX8FBQUdLn/hRdewCuvvILXX38dWVlZCAgIwPz589HWxg0dqStxCm5YqB98dc7vrK4UX50GCeG2oO50BafiSDo5pQ0wWQREBPhgeJifpOcWSwjsPFnJEVFyGx4VLL388stYsWIFli9fjvHjx+P111+Hv78/3nrrrR4fo1KpYDAYHD8xMTGO+wRBwF/+8hc89thjWLRoESZNmoT33nsPJSUl+OKLL1zwjMiTeNJKONGoKCZ5k/TEKbjJ8aGDLkb5c6nJ4dBpVCiua0U+tz4hN+ExwVJ7ezv279+P9PR0x21qtRrp6enIzMzs8XFNTU0YMWIE4uPjsWjRIhw7dsxxX15eHsrKyrqcMyQkBKmpqb2ek4amjuRu966x1NkocdsTjiyRhMSVcFLUV/o5fx8tZowIB8CtT8h9aJVugLOqqqpgsVi6jAwBQExMDE6cONHtY8aOHYu33noLkyZNQn19PV588UXMmjULx44dw/Dhw1FWVuY4x8/PKd7XHaPRCKOxoypyQ4Nt7t5kMsFkMg3o+XVHPJeU56TzOdvPZyttAUdCmJ/H/JskhdumSE6VN7pFm/madg25+1kcWZoYFyjLNS4eGY7Ms9XYkVuBm2cMfINeufH17Bpy9rOz5/SYYGkg0tLSkJaW5vh91qxZGDduHP7xj3/g6aefHvB5n3vuOaxdu/a827ds2QJ/f+mnaDIyMiQ/J52vr34+nKcBoEJV3nFsqjvW67Huwjb7psWxc9XYtGmT0s1x4GvaNeTo52YTkF9t++goO74Hm+RYD9MEAFr8eLICX361CRKVcZINX8+uIUc/t7Q4N9XrMcFSZGQkNBoNysvLu9xeXl4Og8Hg1Dl0Oh2mTp2K06dPA4DjceXl5YiNje1yzilTpvR4njVr1mD16tWO3xsaGhAfH4958+YhOHjgO2//nMlkQkZGBubOnQudzn33IfN0zvSz1SrgkX1bAVhx44JLkRjhGVNxTUYzXj66DQ0mFS6+Yi5C/JR9HfE17Rpy9vMPp6qAfQcwItwfNy6aLem5RVargDfPbEdtiwkxEy7CzMRwWa4zWHw9u4ac/SzODPXFY4IlHx8fTJ8+HVu3bsXixYsBAFarFVu3bsWqVaucOofFYsGRI0ewcOFCAEBSUhIMBgO2bt3qCI4aGhqQlZWFlStX9ngevV4PvV5/3u06nU6WPxi5zktd9dbPpfWtaDVZoVGrkBgVLFnFYrmF6XSIDfFFaX0bCmrbMD3YPZLT+Zp2DTn6+WipbTp6akKorP+Gs0ZF4uvDpThQ2ICLR8f0/QAF8fXsGnL0s7Pn84x3fLvVq1fjjTfewLvvvoucnBysXLkSzc3NWL58OQBg2bJlWLNmjeP4p556Clu2bMHZs2dx4MAB3HrrrSgoKMCdd94JwLZS7r777sMzzzyDjRs34siRI1i2bBni4uIcARkR0JHcnRDu7zGBkmgUK3mThDqvhJPThSPCAAD7CmplvQ6RMzxmZAkAbrrpJlRWVuKJJ55AWVkZpkyZgs2bNzsStAsLC6FWd3yQ1dbWYsWKFSgrK0NYWBimT5+OXbt2Yfz48Y5jHnroITQ3N+Ouu+5CXV0dZs+ejc2bN59XvJKGNkfZAA8oRvlzo6ID8cOpKtZaokETBMGxEk7qYpQ/N92+Iu5AYS2sVgFqtbQlCoj6w6OCJQBYtWpVj9Nu27dv7/L7n//8Z/z5z3/u9XwqlQpPPfUUnnrqKamaSF7IkzbQ/bnO254QDUZRbSuqm9uh06gwLla6/MzupMQGwU+nQWObGacrmzAmJkjW6xH1xrPmE4gUkldlWzHhSQUpReI0HEeWaLDEKbhxscGyV7HXadSO0at9+ZyKI2UxWCJyQkG1OA3nGavgOhttD5aK61rRbDQr3BryZK6aghNNt+ct7WfeEimMwRJRH6xWAQU19pElDwyWwgJ8EBHgAwA4U8nRJRq4w8X1AIBJMlTu7o4YLB0oZLBEymKwRNSHsoY2tJut0KpViAv1zMR/TsWRFM7YXz8pBtfkD01LsAVLeVXNqGoy9nE0kXwYLBH1Id8+BRcf7g+th5UNEI3mHnE0SLXN7ahubgcAJEe5ZoQ1xF/nmEY+wKk4UpBnvvMTuVCBfefzER5YNkA0KoojSzQ4Z6tsr51hoX7w93HdQmpH3hKn4khBDJaI+pDvwcndotH2ZdcMlmigxNeOq0aVRI5giSviSEEMloj6UFDlBSNL9qmMgupmtJksCreGPNGZStuXhpFRrq01JgZLh4vrYTTztUvKYLBE1AdvGFmKDtIjyFcLq9DxfIj6Q0zuHhnt2mApKTIA4QE+aDdbcazEuU1PiaTGYImoF4IgeEXOkkql4h5xNChi2YlRLh5ZUqlUjlVxnIojpTBYIupFZaMRrSYL1CpgeJjnBktAR3FK5i1Rf7WZLCi01xobGe36EVYWpySlMVgi6kW+fVRpWJgffLSe/efCWks0UAXVLbAKQJCvFlGBepdfv/OKOEEQXH59Is9+9yeSmSdvc/JzHRvqNircEvI04hTcyKhAqFQql19/0vAQ6DQqVDYaca6m1eXXJ2KwRNQLb8hXEokjS3lVzTBbrAq3hjyJOBo5ysXJ3SJfnQYXxIUAAPYX1ijSBhraGCwR9cIbVsKJhoX6wU+ngckiOKYXiZzReWRJKcxbIiUxWCLqRcfIkucHS2q1CmPs256cLOdUHDmvI1hS7u9ghj1Y2scVcaQABktEPRAEodPIkudPwwEdlbwZLJGzrFYBZyrsBSkVmoYDgGn2YOlkeSOajWbF2kFDE4Mloh7UtpjQ2GaGSmXbRNcbjGWwRP1U2tCGVpMFOo0KCQr+HcQE+yI2xBdWAThcVK9YO2hoYrBE1ANxVMkQ7AtfnUbh1khjtH0aLreMwRI5R6zcPSIiADqNsh8ZUxNCAQAHz3EqjlyLwRJRD8SyAd6wEk401mAbWcqvbuE+W+QUd8hXEk2JDwUAZBfWKdoOGnoYLBH1IN++ga43rIQTGYJ9EaTXwmIVcLaSe8RR39xhJZxoqn3bk4Pn6licklyKwRJRDzpGlrwnWFKpVBhjYN4SOU/pGkudTYgLgUZtK05ZUt+mdHNoCGGwRD2qbDTCNISLF4q1iLxlJZxoDJO8qR/O2Ecg3WFkyc9HgxR7sM+pOHIlBkt0HkEQ8PftpzHz2e9w1V9/QOEQLWDojSNLADrVWuIecdS7+lYTKhuNAIBkN8hZAjqSvLOZ5E0uxGCJumgzWXDvhmy8sDkXgmAbgl/895+wN39obTFQ32JCbYsJgHcleAMsH0DOE/OVYoL1CPLVKdwamynx9rwljiyRCzFYIofS+lbc+HomNh4qgVatwiNXpWDCsGDUNLdj6RtZ+OxAkdJNdJmCGtuoUlSQHgF6rcKtkZZYmLKwpgWt7VwRRz0740b5SiJxRdyR4vohnSZArsVgiQAAh4vqcM2rP+FIcT3C/HV4/85U3H3ZSHzyP2lYcIEB7RYrVn9yCC9+mzskVqF4a74SAEQG+iA8wMcxckjUE3fKVxIlRwYg2FcLo9nKemHkMgyWCIIg4L4N2ahqMiLFEISNq2bjouQIAIC/jxZ/XzoNv7l8JADgb9+fxjdHy5RsrksUVHlnvhJgWxE3Opp7xFHf3KlsgEitVmGyfXTpYCHzlsg1GCwRMs9W42xVMwL1Wnxyd9p5W3uo1So8tCAF/3NZMgDgrR/zlGimS3nzyBLQUZySwRL1RpyGc6dgCQCmisHSuTpF20FDB4Mlwkd7zgEAFk2JQ3AvSZx3XJwEnUaFfQW1OOLlezN560o4kVg+IJfBEvWg3WxFQY3tS4M75SwBHcUpWT6AXIXB0hBX3WTE5qOlAICbZyb0emx0sC+unhgLAHh7l3ePLnWMLHl3sHSK5QOoB4U1zbBYBQT4aBATrFe6OV2I03Bnq5pRb1+1SiQnBktD3H8OFMFkETBpeAgmDAvp8/jbL04CAHx1qNRRf8XbNBnNqGqyPbcEL52GE2stFde1orGNHzZ0PnE7nKSoAKhUKoVb01V4gI9jijy7qE7ZxtCQwGBpCBMEwTEFd0sfo0qiKfGhmBIfinaLFR9mFcrZPMWIU3Bh/jqE+LlHbRmphfr7IDrINlpwiiviqBvnalsBACPC3XN0dQqTvMmFGCwNYbvP1iCvqhkBPhpcMznO6cctvzgRAPB+VgHazd5X5yTPvhIuKdI9PySk4kjy5vJr6sY5e77S8HA/hVvSPTFYymaSN7kAg6Uh7KM9tpGhRVOH9avw4lUTYhEdpEdloxHf2POdvEm+PVhK9PJgaXS0uCKOI0t0PjFYSgh3z6loR5L3ubohUfuNlMVgaYiqaW7HZnu9JGen4EQ+WjVuvWgEAOCtn/Klbpri8qpsHxLJXh4sjTWw1hL1rNAeLMWHuWewNC42GD5aNepaTI4FGURy8bhgaf369UhMTISvry9SU1OxZ8+eHo994403cMkllyAsLAxhYWFIT08/7/jbb78dKpWqy8+CBQvkfhqK+8/+IrRbrE4ndv/czTMT4KNR49C5Oq/LGcirso20ePvI0hjuEUc9EAQB52rde2TJR6vGBXHBAJi3RPLzqGDp448/xurVq/Hkk0/iwIEDmDx5MubPn4+Kiopuj9++fTtuvvlmfP/998jMzER8fDzmzZuH4uLiLsctWLAApaWljp+PPvrIFU9HMbbEbtsUXF/lAnoSFaR35Dm9sytfqqa5BW8vGyAS94iraDSirqVd4daQO6lsMqLNZIVaBcSFumfOEgBMs0/FHWCwRDLzqGDp5ZdfxooVK7B8+XKMHz8er7/+Ovz9/fHWW291e/wHH3yA3/zmN5gyZQpSUlLwr3/9C1arFVu3bu1ynF6vh8FgcPyEhYW54uko5kBhLc4OILH75269yBZofXe83GsSvetbTKhptgUO3p7gHajXYpj9g5B5S9TZuRrbSrjYED/4aN33Y2LGCNt79b58BkskL4/ZTr29vR379+/HmjVrHLep1Wqkp6cjMzPTqXO0tLTAZDIhPDy8y+3bt29HdHQ0wsLCcOWVV+KZZ55BREREj+cxGo0wGjtqDDU0NAAATCYTTCbpataI55LynACw44RtJO7yMVHQq4UBn398TAAiA31Q1dSO3WcqkJbcc5+5s879fKrcNqoUHaSHzyD6xlOMig5AcV0rjpfUYerwINmvJ9drmroabD/nVdqmZoeH+br1v9WkYR2V6GsaWxHk69qPNL6eXUPOfnb2nB4TLFVVVcFisSAmJqbL7TExMThx4oRT53j44YcRFxeH9PR0x20LFizA9ddfj6SkJJw5cwaPPvoorrrqKmRmZkKj0XR7nueeew5r16497/YtW7bA31/6+f2MjAxJz/fNMTUANQJairFpU9GgzpXsp0ZVkxpvb96L2kTPHl3KyMjAvkoVAA2CVW3YtGmT0k2SnbbJ9lrI2HMMYVVHXHZdqV/T1L2B9vO2ItvfAZqq3f7vIEKvQbVRhTc+y0BKqDKr4vh6dg05+rmlxbnFAR4TLA3W888/jw0bNmD79u3w9fV13L5kyRLH/0+cOBGTJk3CyJEjsX37dsyZM6fbc61ZswarV692/N7Q0ODIhwoODpaszSaTCRkZGZg7dy50OmmKIxrNVjy0dxsAK379i0sxMmpwU02qo2XY8/FhFJqDsHDhxZK00dU69/OpnQXA6bOYNno4Fi68QOmmyc56uBRbPz2CJl0YFi5Mlf16crym6XyD7eednx8FzpUgbdIYLLw8WYYWSmdr8xFsPFwKrWE0Fl45yqXX5uvZNeTsZ3FmqC8eEyxFRkZCo9GgvLy8y+3l5eUwGAy9PvbFF1/E888/j++++w6TJk3q9djk5GRERkbi9OnTPQZLer0eev35eyXpdDpZ/mCkPO/BohoYzVZEBuoxNjZk0NsYXJZigEZ9BGcqm1HWaEK8m66ccYZOp0NhbRsAIDk6aEi8+U1OsE1J55Y3QqXWQKtxTX6KXH8r1NVA+7nI/neQGBno9v9OFyZHYOPhUhw816BYW/l6dg05+tnZ87lv5t7P+Pj4YPr06V2Ss8Vk7bS0tB4f98ILL+Dpp5/G5s2bMWPGjD6vU1RUhOrqasTGxkrSbneTdbYaAJCaFC7Jfk8hfjpMt69I2Z7b/apET5JfPTSqd4sSIwIQ4KNBm8mKs/ZinERF9q1OPOHLj/j+c7CwFhardxSntFgFnKtpYbFNN+IxwRIArF69Gm+88Qbeffdd5OTkYOXKlWhubsby5csBAMuWLeuSAL5u3To8/vjjeOutt5CYmIiysjKUlZWhqcm28qepqQkPPvggdu/ejfz8fGzduhWLFi3CqFGjMH/+fEWeo9x259mCpYuSw/s40nmXp0QBAL7PrZTsnEoQBAF5lUMrWFKrVRhvr1VzrKRe4daQO2g3W1FSLwZL7ls2QDTWEIRAvRbN7RacKHNuSsUdWa0C9hfU4I8bjyH12a245IXv8fB/DsNs8excUG/hUcHSTTfdhBdffBFPPPEEpkyZguzsbGzevNmR9F1YWIjS0o7tN1577TW0t7fj//2//4fY2FjHz4svvggA0Gg0OHz4MK699lqMGTMGd9xxB6ZPn44ffvih22k2T9dutmJ/gW2JbaqEK9euGBsNANh1pgptJotk53W1muZ2NBrNUKnctxCfHC6IsxUlPVrsuR80JJ2SulYIAuCrUyMq0P3fBzVqFaYmhAIADhR4XgkBQRDwjx1ncMkL3+OG1zLxzq58VDXZVlt/sq8Id79/wKPfV72Fx+QsiVatWoVVq1Z1e9/27du7/J6fn9/rufz8/PDtt99K1DL3d6S4Dm0mK8IDfDA6OlCy86YYgmAI9kVZQxuy8mpw2Zgoyc7tSmIxyrgQP/jqul8J6Y0u4MgSddJ5mxMppupdYfqIMPxwqgr7Cmrxq7REpZvTLxv2nsNz39hWdAfqtZg3PgbXTI5Dq8mC+z7Oxnc55Vj25h68cdsMhPgxL0opHjWyRIOz+2wNAGBmojT5SiKVSoUrxKm4E56bt5RnD5aGyhScSBxZOlbcAKuX5HzQwLn7NifdmW4vTrnfw0aWTpU3Yu2XxwAAv71yFPY9lo6Xb5qCK1KisXBiLP7v1zMRpNdiT34NbvpHJioa2hRu8dDFYGkI2X1W+nwl0WVjbFNxnpzkXSBucxLpOR8SUhgdEwgfjRqNRrPjg5KGLsfIkgcFS1PiQ6FW2RLTyz0koGgzWfDbjw6izWTFpWOicH/6mPNGtFOTI/Dx/6QhMlCPE2WNuPmN3TAxh0kRDJaGCJNFnnwl0cWjIqDTqJBf3YI8D11Vle8YWZJuitIT6DRqjDXYKiEfK2He0lBXVOM5K+FEQb46jDXYppM9ZXTp2U05OFHWiMhAH7x042So1d2P9o+PC8ZnK2chIsAHZyqbselIabfHkbwYLA0RR4rr0dJuQai/DmNjpN/WIshXhwsTbSNWnjoVl18lroTznA8JqUwYZvugOVrMvKWhThxZ8qRpOMCz9onbcqwM72UWAABe+uUURAX1nkifEOGP22YlAgDe/DGPJQUUwGBpiMiy5ytdmBje4zeYwRJXxW0/6XklBAQBKLB/SCRGDK2cJaBT3hJHloY8cSrWE8oGdNaRt1SjcEt6V1rfiof+cxgAcNelyU4viFmamgAfrRqHi+qxz0NGz7wJg6UhIstRX0m+zW7FJO/dZ6vR0m6W7TpyqG8HWk1WaNQqj5p+kErnFXH81jp0NbSZUNdi21g0Psyz/g7EYOlYSQNa2913qf2zm06grsWEicNC8Pt5Y51+XESgHtdPHQYAePOHPLmaRz1gsDQEmC1W7M2zfdtKTZI+uVs0MioQw0L90G62OpLJPUVlm220LT7MDzoXbfnhTsbFBkOjVqGqqR0VjUalm0MKOWcfXY0I8EGA3rMqywwP80NMsB5mq4BDRXVKN6dbZyub8NXhEgDA8zdMhI+2f+81v56dBADYcrzM8W9FrjH0PhWGoGMlDWhutyDIV4txsdJt9PtzKpUKs0dFAuiY9vMUlfYFNEOtbIDIV6dxbKrMvKWh65wHroQTqVQqty8h8Nr2MxAEIH1ctGPquz/GxAThktGRsArA2z/lS99A6hGDpSFAnIKbmRgOjUz5SqKZ9pGrPfkeFiy12volcYgGSwAwgXlLQ54nlg3obPoI2/uPOwZLRbUt+PxgMQDgnitGDfg8d9hHlz7Zdw6NbSZJ2kZ9Y7A0BIjFKOXMVxKJwdKRonqPyluqGOIjSwAce8RxZGnoOmcvG5DgYcndoumOFXE1brep7j92nIXZKmD2qEhMtW/+OxCXjYnCqOhANBnN+HjvOQlbSL1hsOTlBEHAoXN1AIDpiQP/A3XW8DA/xIb4wmwVkF1YJ/v1pCLmLA3lYGnCMI4sDXWdtzrxRBPighHkq0VDmxlH3CjoL29ow8f7bIHNYEaVANt0468vto0uvbMr3+2CQm/FYMnLVTYaUd3cDrUKGGeQL19JpFKpHPWWPGUqzmIVUGUfWRqKZQNE4shScV0rapvbFW4NKcETtzrpTKtRI80+gv7jKfcpYfLGzrNoN1sxY0SYJDsoXD9tGML8dSiqbcWWY2UStJD6wmDJy4mjBMlRgfDzcc3msBeKeUt5nhEsldS3wiKooNOoEBfqmdMPUgj21WFEhO1DkqNLQ4/VKnhk9e6fu2S0bZHJD6eqFG6JTU1zOz7IKgQArLpylCT7cvrqNLjpwgQAwNes6O0SDJa83PFS24feeBlXwf2cWJ7gQGEt2s3uv4+RuM3JiHB/2RPg3V1Hkrf7TGGQa1Q0GtFusdUaiw3xVbo5AzZ7tK3e24HCWjQblc+bfOvHPLSaLJg4LMTpApTOmDveVgT4h1NVnIpzAQZLXu64fYRALDroCqOiAhHqr0ObyYqjHvChm18lVu723G/TUnEkeXNkacgR85WGhfpB68G1xhIj/DEs1A8mi6D46HZjmwnv7soHYMtVkmJUSTR5eCiCfbWobzW5bV0pb+K5fxHkFHGEYLwLgyW1uiNvaa8HTMU5RpYYLHVK8nb/IJek1VFjybOnolUqldtMxX12oBiNRjOSowIwb3yMpOfWatS4xD6KtiPXffKzvBWDJS/WZDQ7AgFXTsMBtppOALDXA5K88xwb6A7d5G6ROAKZV9XMGi5DjKduoNud2fZg6cfTygURgiDg3cx8AMDtsxJl2ZNTnNbb4YH7cXoaBkte7IQ9X8kQ7IuIwN53tZaaWG9pb34trG4+n3660hYsjY4OVLglyosM1CMh3B+C4J6F/Ug+4kq44R5aNqCzi0dGQqUCTpY3obyhTZE2/HS6GmcrmxGo1+L6acNlucal9mDpUFEdV7DKjMGSFxNXNLlyCk50QVww/H00qG814WRFo8uv76zGNhNK621vpuJ2H0PdTA9bzUjS8OStTn4uLMAHE+1Tyj8qNBUnjirdMG0YAmXaZ88Q4osUQxAEAfjhtHus/vNWDJa8mBLJ3SKtRu2opuvOH7qnK5oAAME6ASF+OoVb4x46RgXd99+NpFdSZ/vSMDzMs3OWROI+lT8qEEScq2nB1pxyAMCv0hJlvZY4FbeTU3GyYrDkxY6V2pO7XZyvJHIUp3TjYOmUPVgy+Lv3VKEriaUfDp2rR5vJonBryBXMFivK7NNVw7yk1lhH3lIVBMG1f9/vZxXAKthqPo2SeXr/0k55S65+nkMJgyUvZbJYcbLMFggoMQ0HdARLe/Nr3PaPWBxZMnjH54MkEsL9ER2kR7vFimz7Vjnk3cobjbBYBeg0KkS5OL9RLtNHhMFPp0FloxG55a5LBWgzWRx7ti2TeVQJAGYkdjzPnFL3TXnwdAMKlsxmM7777jv84x//QGOj7R+npKQETU1NkjaOBu5MZRPaLVYE6bWK7fM0NSEUOo0K5Q1Gx0obd3PS/ibKkaUOKpWKeUtDTEmdrXJ3bIifLKu2lKDXahyv4x9Oum4qbuOhEtS1mDAs1A9XpkTLfj29VoNZI21bvHBVnHz6HSwVFBRg4sSJWLRoEe655x5UVtr+cdatW4ff//73kjeQBuZYsS1faVxssGJvfr46DSYNDwXgvh+6p8rFkSUGS52lMlgaUoprbcFSXKjnVu7ujqPekovylgRBcBSh/FXaCJftCHDZWHEqrsIl1xuK+h0s3XvvvZgxYwZqa2vh59cxd3Hddddh69atkjaOBs6xzYlCU3Aidx6haDaaUWz/Rs1puK5mJtm+qe4vqIXJ4v5b1tDgiH8Hw0I9fyVcZ2Le0p68apfk3x0orMWxkgbotWrcNCNe9uuJxCTvffm1aHKDLV68Ub+DpR9++AGPPfYYfHx8utyemJiI4uJiyRpGg3NcwbIBnblzccozlbZRpchAHwRwIVwXo6NtW9a0miw4Wsxq3t6uI1jyrpGlsTFBiArSo81kxQEX1A17+6d8AMCiKXEIC/Dp/WAJjYgIQGKEP8xWAbtYQkAW/Q6WrFYrLJbzI/SioiIEBQVJ0igaHEEQOrY5UWglnGiavXxAfnULKhuNirbl58QpuFGsr3SeLlvWuGGgS9ISc5aGeUnZAJFKpcKl9i1BvjlaJuu1impbHNe4fVaSrNfqDqt5y6vfwdK8efPwl7/8xfG7SqVCU1MTnnzySSxcuFDKttEAFde1oqHNDJ1GhTExygawIX46jImxLZ09UOheFaHFYplyL+31VDM9oPQDSUMMluK8pGxAZ4umxAEAvjxcAqNZvqm4t3/Kh8UqYPaoSEVG9DvyllhCQA79DpZeeukl/PTTTxg/fjza2tpwyy23OKbg1q1bJ0cbqZ/EKbhR0UHw0SpfHWL6CNuHrrttn3GaI0u96pxv5u5b1tDACYLQKcHb+4Kli0dFIjpIj7oWE74/Ic+oS2ObyVEu4M5LXD+qBAAXJUdAq1ahqLYVJfXKbPHizfr9STp8+HAcOnQIjz76KO6//35MnToVzz//PA4ePIjoaPmXSVLfjilYubs7M+xTcfvcbDpHLEjJkaXuiVvWNLSZXVqnhlyrodWM5nbbiIu3FKTsTKNW4bqpwwAAnx0okuUaH+8rRpPRjNHRgY7pMFfz99EiJdY2k5BdWKdIG7zZgDas0Wq1uPXWW6VuC0nEsRJO4Xwl0YxEW7B0tLgBbSYLfHUahVsEtLZbHBuHjooORHWOwg1yQ+KWNT+cqsKevBqMc5PXE0lLTO6OCPBxi79NOVw/bTj+sfMsvs+tQE1zO8IlTL62WIF3MwsA2EaVVCrl6lRNiQ/F0eIGZJ+rxdWTYhVrhzfqd7D03nvv9Xr/smXLBtwYkoa7rIQTJYT7IzJQj6omI44U1zsSh5V0prIJggCEB/ggwoWrVjxNalK4I1i6bVai0s0hGRR7cb6SaKwhCBOGBeNocQO+PFQi6Wv5YLUKZQ1GRAb6YNGUYZKddyCmxIfh/d2FrLwvg34HS/fee2+X300mE1paWuDj4wN/f38GSwqra2l3vPm5S7CkUqkwfUQovj1Wjv0FtW4RLJ1icrdTxH+rrDzbljVKfmsmeThWwnlxsAQA108djqPFx/HZgSLJgiVBEPB9qS2bZVlaouIjc1PiQwEAR4rrYbJYodMon7PqLfrdk7W1tV1+mpqakJubi9mzZ+Ojjz6So43UD+IUXHy4H4J93ad40Ax7kve+fPdI8hbLBogr9ah7k+ND4aNRo6rJiPxq99yyhgZnKIwsAcC1U+KgUatwqKgepyukycHbk1+LomYVfHVq3HrRCEnOORjJkQEI8tWizWRFbhnzDKUkSdg5evRoPP/88+eNOpHriUHA2Bj3GFUSTbfnLR0orHWLZa1icvfoaNYG642vTuP4tronr1rZxpAsir20xtLPRQbqcbk9+fqzA9IUUH7rJ1uu0nVT4iTNgxootVqFyfYtpjgVJy3Jxui0Wi1KSkqkOl2P1q9fj8TERPj6+iI1NRV79uzp9fhPP/0UKSkp8PX1xcSJE7Fp06Yu9wuCgCeeeAKxsbHw8/NDeno6Tp06JedTkJU4veRuIyYT4kLgo1WjprkdZ6ualW4OTtlXd43mNFyfUpNto4JyLbsmZZV4afXu7lw/bTgA4PODxYMuh3G4qA7bciuhgoDls5QfVRKJX24YLEmr38HSxo0bu/z897//xeuvv45bb70VF198sRxtdPj444+xevVqPPnkkzhw4AAmT56M+fPno6Ki+80Dd+3ahZtvvhl33HEHDh48iMWLF2Px4sU4evSo45gXXngBr7zyCl5//XVkZWUhICAA8+fPR1ubZ9apEEeWRrtZsOSjVWPy8BAAwH6Fp+LaTBYU1thXwrlZP7mjBRMMAIDvcyvQ2GZSuDUkNW+usfRzc8ZFI9hXi9L6NmSeHfhIqdUq4In/HgMATIsUkBTpPrXaxGDpEIMlSfU7WBIDDvHn+uuvxx//+EdMmjQJb731lhxtdHj55ZexYsUKLF++HOPHj8frr78Of3//Hq/717/+FQsWLMCDDz6IcePG4emnn8a0adPwt7/9DYBtVOkvf/kLHnvsMSxatAiTJk3Ce++9h5KSEnzxxReyPhe5nBZrB0W53/SSuxSnPFvZDKsAhPrrEBWoV7QtnmB8bDCSowJgNFvxXU65om1pM1mwYU8h7npvH17bfoabhg6S0WxBhX0bIm9P8AZs08q/mGyr6P2fQdRc+s+BImSfq0OAjwaLRrjXRtNTEkIBAKcrm/jlRkL9Xg1ntSrzwmhvb8f+/fuxZs0ax21qtRrp6enIzMzs9jGZmZlYvXp1l9vmz5/vCITy8vJQVlaG9PR0x/0hISFITU1FZmYmlixZ0u15jUYjjMaOfc4aGmxJ1SaTCSaTdC9O8VzOnrOmuR3Vze0AgIQwH0nbIoUpw20B3N78GkXbdqK0DoCtcrfZbO53Pw9FV0+Iwavfn8V/DxbjFxNiBnyegfZ1bUs7PtpThP/LKkRVk+01vuV4Of6x4wxuS0vAsosSEOznPgsalOZsPxfZR1j1WjWCfFRD4m9g0SQDPswqxJeHSnDX7ESM7GcF/4ZWE57/5gQAYOWliQhpznWrfgvRqzE81BdFdW04kF+NWSMjlG7SoMn5Hu3sOQdUlFIJVVVVsFgsiInp+kYdExODEydOdPuYsrKybo8vKytz3C/e1tMx3Xnuueewdu3a827fsmUL/P39+34y/ZSRkeHUcWcaAECLcL2A7d9tkbwdg9VsAgAtzlY149P/bkKAQp9tmwvVANTwaa3pksPmbD8PRYEtAKDFzlOVkvzbOdvXggBsLlJhW4ka7VZb2YJQHwEzIgUcrlGhotWEv247g3/sOI30YVakxwlgdYMOffXzqXoVAA1CtBZ88803rmmUwgQBGBeqRk6dGivf/hG/u8ACdT9eM5/lqVHdrEaMn4DYxlxA7X7vHVEaNYqgxidb96AuV/kFNVKRo59bWpxb5etUsPTz0ZnevPzyy04f66nWrFnTpU8aGhoQHx+PefPmIThYulVoJpMJGRkZmDt3LnS6vj+dPtp7DjiWg4kjorBw4TTJ2iGlf+X/hLNVzQgfOwNzUpTZHuerD7OB4gpcMWMcFqaN6Hc/D1Wfl+3CifImCMMmYeGM4QM6R3/7+uXvTmFzUR4AIMUQhDtnJ2LhhBjoNGpYrAK+OVqGv+84i1MVzfiqUIMrZk7EtZNZudjZfv7sYDFw/BhGD4vEwoXTXdhCZU29uBULX92FvEYLaiImYNlFCU49LresET9m7QYg4PlfzkDqiGC3fO8oDy3AwW9y0RZgwMKFU5VuzqDJ+R4tzgz1xalg6eDBg06dTM6CdZGRkdBoNCgv75ozUV5eDoPB0O1jDAZDr8eL/y0vL0dsbGyXY6ZMmdJjW/R6PfT683NddDqdLH8wzp73bJUtUXOsIdit/nA7m5EYhrNVzcguasSCicpUuz1jX403Lja0Sz/J9e/nLa6ZMgwnvs3FpmPlWJo2uM1CnenrN3aexWs7bIHSH68Zj9tmJXZ5j9EBuG56AhZNjcf/bsnFa9vP4PlvT2LuhFi3qjGmpL76uazBNgUxPMx/SL32R0Tp8MjCcXj8i6N4KeMU5l0Qi/jw3mcFBEHA05tyYbEKuGqCAVeMMzimcNztvWO6vZjsoaIGaLVarykmK0c/O3s+pxK8v//+e6d+tm3bNqhG98bHxwfTp0/H1q1bHbdZrVZs3boVaWlp3T4mLS2ty/GAbRhPPD4pKQkGg6HLMQ0NDcjKyurxnO6sI7nbfVd4zXAkeSuzqa7RbEGBvbiiu60YdHfXTLIlxmaeqUZlo7GPowfn472F+NMm24Z9D84fi9sv7nnPLbVahfvSRyMpMgCVjUb8OeOkrG3zJiVDpCBld5bOTMDMxHC0tFvw6OdH+qz/9kFWIbLyauCrU+MPV49zUSsH5oK4EGjVKlQ1GR11tGhwPKoW+urVq/HGG2/g3XffRU5ODlauXInm5mYsX74cgG1fus4J4Pfeey82b96Ml156CSdOnMAf//hH7Nu3D6tWrQJgGwm777778Mwzz2Djxo04cuQIli1bhri4OCxevFiJpzgoji083DgIEItTHiqqR7vZ9YsFzlQ0w2IVEOSrRXQQV8L1R0KEPybHh8IqAN8cLZXtOpuOlGLNZ0cAAP9zaTJ+c/nIPh+j12qw9toLAADv7sp37I9IvSupHxoFKbujVqvw/A0T4aNV44dTVfj3/u5Xx1msAp7blIPHvrCVnLnn8lEYHiZ9bqqUfHUapMTaFtSw3pI0BpTgvW/fPnzyyScoLCxEe3t7l/s+++wzSRrWnZtuugmVlZV44oknUFZWhilTpmDz5s2OBO3CwkKo1R3x36xZs/Dhhx/isccew6OPPorRo0fjiy++wIQJExzHPPTQQ2hubsZdd92Furo6zJ49G5s3b4avr2cVaGtoM6G8wfZt3533O0uODECYvw61LSYcLanHtIQwl17/SHEdAFuRTG8ZmnalaybF4tC5Onx5qATL0hIlP3/mmWrcu+EgrAJw88x4PHJVitP/TpeOicLCiQZsOlKGx/97FJ/+TxrU/cncHYI6aix51vudVJKjAnF/+his23wCT391HGH+Ppg9OtKxx1tdSzt++9FB/HCqCgBw92Uj8ZsrRinZZKdNiQ/F0eIGZBfW4Rf2UWEauH6PLG3YsAGzZs1CTk4OPv/8c5hMJhw7dgzbtm1DSEiIHG3sYtWqVSgoKIDRaERWVhZSU1Md923fvh3vvPNOl+NvvPFG5Obmwmg04ujRo1i4cGGX+1UqFZ566imUlZWhra0N3333HcaMGSP785CaOAVnCPZ163wN26a6tqm4PXmun4o7VFQPAJgUL/9r1Rv9YlIcVCpgb36tYwpHKqX1rVj14QGYLAKunhiLZxZP7HdA+/gvxsPfR4P9BbX49yDq6AwFgiA4pmiGh7r3SImcVlyShAnDgtHQZsad7+3D9KczcM8HB/BhViGu/dtP+OFUFfx0Gvztlql45KoUaDwkAJ8SL47i1ynbEC/R72Dp2WefxZ///Gd8+eWX8PHxwV//+lecOHECv/zlL5GQ4NyKApLeaTet3N2di+zbZ+weRAXdgTpsf+MQ90+i/jGE+OJCe/Lo14elm4prN1txzwcHUN3cjnGxwXjpl5MH9KEUG+KHe+eMBgA8/80J1LW09/GIoaumuR1GsxUqFRATMnSnpLUaNd66/ULcljYChmBfNLdb8PWRUjz6+REU1rQgPtwPn/1mlseNzoiVvI8U18Nkca/CmZ6o38HSmTNncPXVVwOwJV03NzdDpVLh/vvvxz//+U/JG0jOEfOVRrpxcrfoomRbkbS9eTUwu/CPuM1kwYlSWz9NGs6RpYG6xl4B+cvD0u0F+eymHBworEOQrxav3zrNMQ0yEL+enYTR0YGoaW7H37adlqyN3kYcVYoK1EOvHXh/e4PoIF+sXTQBux65El/cczFWXj4SY2OCMHd8DDbeMxvjYt1rY3JnJEcGIMhXizaTFblljUo3x+P1O1gKCwtDY6Ot44cNG+bYZ62urs7p4k4kvVMVnjOyNC42GMG+WjS3W3DUhYm4x0sbYLYKiAjwGRJbO8jlqgkGaNQqHC6qx67TVYM+33+zi/HOrnwAwJ9/OQUjIga3z5ZOo8YjV6UAAL7ILnZpQO5JHBvoDsHk7p6o1SpMiQ/FwwtS8O39l+KNZTMQFuCjdLMGRHwuAJO8peB0sCQGRZdeeqmjiuaNN96Ie++9FytWrMDNN9+MOXPmyNNK6pNjA91o99sT7uc0ahVmJtlGl1w5FXfY/oYxaTiTuwcjMlCPpam2KffH/nt0UKsaT5Y34pH/2Fa+3XPFSKSPH/hWKp1dOiYKYf46VDW1Y/dZZcpUuLuiIbSB7lAlphscsedq0sA5HSxNmjQJqampmDhxIm688UYAwB/+8AesXr0a5eXluOGGG/Dmm2/K1lDqWUu72TGkPtqNV8J1pkTe0mExuZv5SoP2wLyxiAz0wdnKZrzxw9kBnaO6yYj/+b/9aDVZcPGoCKyeO1ay9uk0alw10VZoduOhYsnO601K6toADI0NdIeq8XG26cOcMpbSGCyng6UdO3bgggsuwHPPPYdx48bhtttuw08//YRHHnkEGzduxEsvvYSwMNcuAyebMxW2itSRgT4eM2SsRN6SuCpkMlfCDVqIn85RmO/VbadwrqZ/U/CNbWbc/vZe5FU1Y1ioH15ZMlXyVUbX2nOrNh8tg9FskfTc3sAxDcdgyWuJuVa5ZY2cjh4kp4OlSy65BG+99RZKS0vx6quvIj8/H5dddhnGjBmDdevW9brxLMnLk5K7Ra7OW2psM+GsfZsTjixJY/GUYUhNCkebyYq1Xx53+nEmK/CbDw/iSHE9wgN88N4dMxERKP1qrAsTwxETrEdDmxk7Tw4+t8rbFA/h6t1DxYhwf/j7aGA0W5Ff3ax0czxavxO8AwICsHz5cuzYsQMnT57EjTfeiPXr1yMhIQHXXnutHG2kPnhScrfI1XlLR4rrIQi2b9GRMnwwD0UqlQrPLJ4ArVqF73LKkXG8vM/HmC1WvHtSjd15tQjw0eDd5TNlC/I1apVjufeXh6RbuectOrY6GZoFKYcCtVqFsQZbHuvxUq6IG4xBbXcyatQoPProo3jssccQFBSEr7/+Wqp2UT94UnJ3Z67MWxLzlTgFJ63RMUG485JkAMAfNx5Da3vP012CIOCJL3NwpFYNnUaFN26bgYkyl3AQyxxkHC9HS7tZ1mt5ktZ2C6qbbTWohnJByqFAnIrLKWXe0mAMOFjauXMnbr/9dhgMBjz44IO4/vrr8dNPP0nZNnLSmUoxWPKckSUASBvpurwlsRglp+Ck97s5oxAX4oviulbc8NoufHmoBBZrx6akVquAb46U4pq//YhP9xdDBQF/+eUkzBoZKXvbJg8PQUK4P1pNFmzNqZD9ep5C3BMuwEeDYL8B7XpFHoLBkjT6FSyVlJTg2WefxZgxY3D55Zfj9OnTeOWVV1BSUoI33ngDF110kVztpB60mSwosM9Fu/OecN0ZZwhGiJ/OJXlLh86JK+E4siQ1fx8t1v2/SfDTaXC8tAG//eggrnxpO97fXYD/7C/C3D/vwMoPDuBocQP8dGrcMsqKeRKVCOiLSqXCNZNtq+I4FdehpFO+EstoeLfx9mCJm0sPjtNfKa666ip89913iIyMxLJly/DrX/8aY8dKt9SXBiavqhlWAQj21SIqyLNycdRqFWYmhSPjeDl2n612FFCTWlWTEcV1rVCpgInDGCzJ4ZLRUdj1yJV4NzMf7+7KR0F1i2OXdgAI8tXi9lmJWDpzOLJ2fOfStl0zOQ7rvz+D7bmVaGgzufXeia5SwuTuISPFEASVCqhoNKK6ySjLYoqhwOmRJZ1Oh3//+98oKirCunXrGCi5iY7k7iCP/IYolhCQM29JnIKzlf/nB6VcwgJ8cF/6GPz0yJX44zXjMSLCH3Ehvnh4QQp2PXIlHpg3FhEKlLYYGxOE0dGBaLdY8e1RrtoFOmosMVjyfgF6LUaE2/LScpjkPWBOjyxt3LhRznbQAJ0ut734PS1fSSQmeYt5S1rNoNYcdEucguPmua7h76PF7Rcn4faLk5RuCgDbVNy1k+PwUsZJfHm4FDfOiFe6SYpzjCyFcCXcUDAuNhj51S3IKW3A7NHy5wp6I+k/mcilTtuTuz0tX0nkiryljuRuTsENVb+wr4r76XQVauyrwIYyMcGbI0tDA5O8B4/BkocTywZ4arAk5i0B8kzFCYLQsc2JTDlR5P6SIgMwNiYIFqvg0i123FWpfRouljWWhgQxWDrOYGnAGCx5MLOloyqrpwZLQEfeUuYZ6T/EiutaUd3cDq1a5VgVQkOTOOWbNcSDJUEQHNW7udXJ0DAu1laD73RFE7f+GSAGSx6sqLYVJosAX50acSGe+6Y3a2RHkrfUhQPFUaWxhiD46jSSnps8S6o9KM/Kq1G4JcqqbTHBaLbVNTMwZ2lIGBbqh2BfLcxWAafti4KofxgseTCxGGVSZCDUEm9C6kophiAMD/OD0WzFzpOVkp77EItRkt2FibaRpRNljahrGbp5S2Jyd2SgHnotv0AMBSqVqlPeElfEDQSDJQ8mBksjowIUbsngqFQqLLjAAMC2Q7yUdtun9qYyX2nIiwrSO/5W9gzh0aUSxxQcR5WGEiZ5Dw6DJQ92ttKWr5Qs00akrjR/gi1Y2nqiAu1mabY+KW9owyH7NNzlY6MkOSd5Nk7FdQRLsR48dU/9N57B0qAwWPJg3jKyBADTEsIQGahHY5sZmRIl4H6XUw4AmBIfiuhgfosmINW+8jIrb+gmeZfWsyDlUNR5ZEkQhD6Opp9jsOTBzthHlkZ6wciSRq3CvAts+4VJNRWXcdwWLM110T5k5P5Sk2wjS8dLGtDQZlK4Ncoodmx1wi8QQ8nomEBo1CrUtphQ3mBUujkeh8GSh6ptbncU10v2gpElAJhvz1vKOF7eZdf6gWg2mrHrtG30gMESiQwhvhgR4Q+rAOzPr1W6OYrgyNLQ5KvTIDnS9llxvLRe4dZ4HgZLHupslW0KLi7EF/4+Tu9a49bSkiMQ5KtFVZMRBwoH90G282Ql2i1WjIjw99itYEge4lTc7iE6FdeRs8SRpaFmfBxXxA0UgyUPdabCPgXnRYGAj1aN9HG2UaDBbniaYc9XmjsuxiM3GCb5iFNxWWeHXpK32WJFeYNtZIkFKYceVvIeOAZLHupMlZjc7T3BEgDMF/OWjpUNOAnRbLFi24kKAEA6p+DoZ1LtlbyPFNej2ShtEVR3V95ohFUAdBoVIgP1SjeHXIzlAwaOwZKHEkeWvCVfSXTpmCj46tQoqm3FsQFurLuvoBZ1LSaE+uswY0SYxC0kTzc8zB/DQv1gsQrYXzC08pbEKThDiK9HF7KlgRG3PcmrakZrO7c96Q8GSx7qbKV3jiz5+2hx6WhbTaQtxwY2FfedfRXclSnR0Gr4EqfziXlLQ604JWssDW3RQb6ICPCBIIDbnvQTP0k8kMliRWFNCwDvC5YAYIG9QOXmAQRLgiA48pXmcQqOeiBOxQ21eksldcxXGurGxNhGl06UcSquPxgseaCC6haYrQICfDSICfa+vIM5KTHQqlU4Wd7kGEFz1qmKJhRUt8BHq8Ylo1m1m7onJnkfOlePNtPQmY4orWeNpaFurMEWLJ0s54q4/mCw5IHEyt3JUYFeudIrxF+HWaMiAQBv/pjXr8eKhSgvHhmBAL13lFQg6Y2I8EdMsB7tFuugy1R4Ek7DkRgs5ZZzGq4/GCx5oLOOyt3eldzd2T2XjwQAbNh7zhEcOqOjardBlnaRd1CpVJg5BEsIcBqOxGm4XE7D9QuDJQ90xkuTuztLTY7AnJRoWKwCXvw216nHFFa3IPtcHQAgfVy0jK0jbzDTnuQ9pEaW7NNwsZyGG7LGxNg+N8objKhraVe4NZ6DwZIH6jwN580eWpACtQr45mhZnx9ogiBgzeeHAQCXjI7kxrnUp8nDQwDY6i0NhY1FW9rNqGux7YfHrU6GriBfnWNk8SSn4pzmMcFSTU0Nli5diuDgYISGhuKOO+5AU1PP/9A1NTX47W9/i7Fjx8LPzw8JCQn43e9+h/r6rnviqFSq8342bNgg99MZMEEQOqbhor13Gg6wza3fMG04AOD5TSd6/UD7aM85/HS6Gr46NZ5eNMFVTSQPNtYQBK1ahboWk2NzWW8mTsEF6rUI9tUp3BpSUkfeEpO8neUxwdLSpUtx7NgxZGRk4KuvvsLOnTtx11139Xh8SUkJSkpK8OKLL+Lo0aN45513sHnzZtxxxx3nHfv222+jtLTU8bN48WIZn8ngVDe3o77VBJUKSIzw7mAJAFbPGwO9Vo09+TWOqtw/V1Tbgj99fRwA8OD8FCRGen+/0ODptRpH/sbRYu/fWJQr4UjEvKX+84jlQjk5Odi8eTP27t2LGTNmAABeffVVLFy4EC+++CLi4uLOe8yECRPwn//8x/H7yJEj8ac//Qm33norzGYztNqOpx4aGgqDwTMSgs/YC4kND/ODr06jcGvkFxvih+UXJ+H1HWewbvMJXD42GppOlYcFQcCaz46gud2CGSPCcPusROUaSx5n4rAQHC9twJHieiyYEKt0c2TFlXAkGmuwpXCcLOM0nLM8IljKzMxEaGioI1ACgPT0dKjVamRlZeG6665z6jz19fUIDg7uEigBwD333IM777wTycnJuPvuu7F8+fJel+QbjUYYjUbH7w0NtujcZDLBZDL156n1SjxX53OeKrddKzkiQNJrubM7L07AR3sKcLK8CW/sPI1bUxOg19oGRT/ZV4QfTlVBr1Xj2cXjYbWYYe1n2Zzu+pnk4W59PS7W9qFx+Fyd27RJCt3187lq2/S9IVjvVc9VSe72enbWyAh/AEBueQPa29vdvgSNnP3s7Dk9IlgqKytDdHTX1U1arRbh4eEoK3OuynNVVRWefvrp86bunnrqKVx55ZXw9/fHli1b8Jvf/AZNTU343e9+1+O5nnvuOaxdu/a827ds2QJ/f3+n2tMfGRkZjv//Ll8NQA1VUwU2bdok+bXc1eXRKvy3QIPnN5/ES9/mIjFIQHIQsKNMBUCFq4aZkLNnB3IGcY3O/Uzycpe+bmgEAC0O5lfh6683wc0/M/qtcz/vPW1772gsL8CmTfmKtckbucvr2VkmK6CGBvWtZmz47zcI8VG6Rc6Ro59bWlqcOk7RYOmRRx7BunXrej0mJ2cwH382DQ0NuPrqqzF+/Hj88Y9/7HLf448/7vj/qVOnorm5Gf/7v//ba7C0Zs0arF69usv54+PjMW/ePAQHBw+6vSKTyYSMjAzMnTsXOp0tIfOz/zsAlFbhygsvwMIL4yW7lrubY7YiaHMuNh0tQ02zCacaVDhln26fGh+C5389s8v0XH90188kD3fr6zaTBX89vg1NZmDa7CsRG+Id+Tzd9fPHb+8DKmtw+YWTsXDq+akL1H/u9nruj7+f+Qlnq5ox/IKZuGR0pNLN6ZWc/SzODPVF0WDpgQcewO23397rMcnJyTAYDKio6JrcazabUVNT02euUWNjIxYsWICgoCB8/vnnfXZ0amoqnn76aRiNRuj13W8lotfru71Pp9PJ8gfT+bz51bYoeIwhxOP+OAdDpwOeuW4Snl48EWcqm7D7bA2y8mpQ2diG56+fBF/94L8ayfXvR+dzl77W6XQYHR2IE2WNyClvRkJkkNJNklTnfi5rsKUODI8IcIu+9ybu8nruj3GxwThb1YwzVa24crxntF2Ofnb2fIoGS1FRUYiK6nv/rrS0NNTV1WH//v2YPn06AGDbtm2wWq1ITU3t8XENDQ2YP38+9Ho9Nm7cCF/fvr81ZmdnIywsrMdASUltJgvOefEGus5QqVQYFR2EUdFBuPWiEUo3h7zAxGEhOFHWiKPF9Zh/gWcs9OgvQRAc5RFYvZsA24q4r4+UsnyAkzyidMC4ceOwYMECrFixAnv27MFPP/2EVatWYcmSJY6VcMXFxUhJScGePXsA2AKlefPmobm5GW+++SYaGhpQVlaGsrIyWCy2DOAvv/wS//rXv3D06FGcPn0ar732Gp599ln89re/Vey59qagugVWAQj21SIy0EMmmYnc3ER7cUpvLh9Q22KC0WwFABi8ZKqRBkdcEZdbxmDJGR6R4A0AH3zwAVatWoU5c+ZArVbjhhtuwCuvvOK432QyITc315GsdeDAAWRlZQEARo0a1eVceXl5SExMhE6nw/r163H//fdDEASMGjUKL7/8MlasWOG6J9YP3r6BLpESJgwTK3k3QBAEr/zbEssGRAbqodd6f8kR6ptYa+lURSMsVmHAOZ9DhccES+Hh4fjwww97vD8xMbFLhefLL7+8zy0MFixYgAULFkjWRrmdttdYGhU9NKfgiOQwzhAMtQqoajKivMHolSMvJY4pOO97bjQwIyICoNeq0Way4lxNC4v59sEjpuHIRgyWhmq+EpEc/Hw0GB3t3ZW8WZCSfk6jVmG0fVNd5i31jcGSB+HIEpE8OqbivDNYKq237QvHDXSps45tTxgs9YXBkoewWgWcrWKwRCSHicNs9dG8dWRJXAnHfeGosxRuqOs0BkseoriuFW0mK3w0asSH8dshkZS8fWSJ03DUHXFk6SRHlvrEYMlDiFNwSZEB0Gr4z0YkpfFxtiTvikYjKhralG6O5ErqbM9pGL9oUSdj7SNLeVXNMJr7uanmEMNPXQ/BfCUi+fj7aB0LJ46WeNfoUrvZivJGe7DEnCXqxBDsiyBfLcxWAWcrm5VujltjsOQhHCvhGCwRyWKiOBVX5NxeUZ6irL4NggD4aNUsZktdqFQqR97SSeYt9YrBkoc4XcmRJSI5XeCleUtFdbZCvcNC/byy4CYNjpi3dIJ5S71isOQBBEHomIZjjSUiWYgjS8e8bBquuJZ7wlHPHCviGCz1isGSB6hubkd9qwkqFZAcxSqrRHK4IC4YKpWtJlFlo1Hp5kjGkdzNYIm6MdZgK5vBYKl3DJY8wBl74l18mD98ddzXiUgOAXotkiJsX0ZySr0nb6lYnIbjSjjqxlj7NFxxXSsa2kwKt8Z9MVjyAKftwRLzlYjkNdYLpySK6zgNRz0L8dchzr4foje97qXGYMkDnGGwROQSKfYpCW9KdhVzlrjVCfVE/JLgTa97qTFY8gBnKpncTeQKjpGlcu+YhrNaBZTY94Ubzmk46kFH3pJ3vO7lwGDJA4gjS6yxRCQvcWXQqfImWKyCwq0ZvOrmdrSbrVCrAEMI94Wj7o2LtY8slXJkqScMltxcmxkob7CtzOE0HJG8EsL94afTwGi2Ir/a8ysai/lKMcG+0HGbJOrB2E4b6gqC539JkAP/etxcue29DlFBeoT46ZRtDJGXU6tVGBNj+1LiDcmuLBtAzkiODIRWrUJjm9kxbUtdMVhyc+Wttoq7zFcicg1vSnYtrmdyN/XNR6t2zFyc8KKyGVJisOTmysRgiVNwRC4hJrt6w4dGqTiyxORu6oM3fUmQA4MlNydOwzFYInKNlE75G56umNNw5KQUVvLuFYMlN1fOkSUilxK/YRfWtKCl3axwawanRCxIyZEl6kOKY2TJ80dU5cBgyY0ZzVZU2XPtGCwRuUZkoB6RgT4QBOBkeZPSzRmUYrHGEkeWqA/il4Szlc0wmi0Kt8b9MFhyYwXVzRCgQqBei+ggvdLNIRoyUrygSF+rGWhss42MMcGb+hIb4otgXy3MVgFnKjy/bIbUGCy5MUcxyqgAqFQqhVtDNHR4Q7Jrra08G0L9dQjQa5VtDLk9lUrV8SXBSyrYS4nBkhs73SlYIiLX8YYNdWvabV+wmNxNzvKGLwlyYbDkxsShUAZLRK6V0ulDw1MrGosjSwyWyFkpbrrtyQ+nqnCuyZbHqxQGS25M3ECXwRKRa42ODoJKBdQ0t6Oyyah0cwakxmgfWeJKOHJSipuOqD7y+TG8eESLHAVrnzFYcmOxob4I0Qms3k3kYn4+GiRG2L6kuNsHh7M4skT9NSbGFiyVNbShrqVd4dbYVDcZUdFohAqCYysiJTBYcmP/vHUanpphwYgIf6WbQjTkjI1xz2/Zzqo1MmeJ+ifIV4fh9pFId8lbEtsR6Qv4+yi3UIHBEhFRNzw92bVGHFniNBz1g7tNxYlTb3H+yuYOMlgiIuqGu31o9IfRbEWDiSNL1H/u9iUhx55szmCJiMgNiR8aJ8sbYbF61oq4Mnvlbl+dGuEBPgq3hjyJWGvJXbY9EUeWhim8zonBEhFRN0ZEBMBXp4bRbEVBtWdVNC627wkXF+LHgrbUL+KI6smyRlgV/pJgslhxusK2KpwjS0REbkijVjlWB7nLlISzSuwjS3Ghvgq3hDxNUmQAfLRqNLdbUFjTomhbzlY2o91iRYBeg3CFd/xisERE1IOxnhos2UeWhjFYon7SatQYZx9dOlJcr2hbxKnAlBhb3TMleUywVFNTg6VLlyI4OBihoaG444470NTU+47gl19+OVQqVZefu+++u8sxhYWFuPrqq+Hv74/o6Gg8+OCDMJvNcj4VIvIQHdueuEf+hrOK6+wjSyFM7qb+mzAsBABwVOFgSUzuFqcGleQxuysuXboUpaWlyMjIgMlkwvLly3HXXXfhww8/7PVxK1aswFNPPeX43d+/o2aRxWLB1VdfDYPBgF27dqG0tBTLli2DTqfDs88+K9tzISLPIAZLp8p7/2LmbjiyRIMx0R4sKT2yJCZ3jzUEApWKNsUzRpZycnKwefNm/Otf/0Jqaipmz56NV199FRs2bEBJSUmvj/X394fBYHD8BAcHO+7bsmULjh8/jvfffx9TpkzBVVddhaeffhrr169He7t7VC8lIuWI03D51c1oM1kUbo3zHCNLLBtAA9B5ZEnJvRE7T8MpzSNGljIzMxEaGooZM2Y4bktPT4darUZWVhauu+66Hh/7wQcf4P3334fBYMA111yDxx9/3DG6lJmZiYkTJyImJsZx/Pz587Fy5UocO3YMU6dO7facRqMRRmPHflENDbZ/UJPJBJPJNKjn2pl4LinPSedjP7uOp/V1qK8aoX461LWacKKkDhfEBff9IIVZrQJK7QneUQFaj+lrT+Rpr2dnJYX7QqdRoaHNjDMVDRgR7vpdJGqa21HeYPucTYrQoxTy9LOz5/SIYKmsrAzR0dFdbtNqtQgPD0dZWVmPj7vlllswYsQIxMXF4fDhw3j44YeRm5uLzz77zHHezoESAMfvvZ33ueeew9q1a8+7fcuWLV2m+aSSkZEh+TnpfOxn1/Gkvo7QalAHFf6d8RMKoty/3lJ9O2C2aqGGgCN7fsBxVg6QnSe9np1l8NXgXLMK73+1A1MjXf+6P1mvAqBBpF7Arh3fA5Cnn1tanFvxp2iw9Mgjj2DdunW9HpOTkzPg8991112O/584cSJiY2MxZ84cnDlzBiNHjhzwedesWYPVq1c7fm9oaEB8fDzmzZvXZZpvsEwmEzIyMjB37lzodDrJzktdsZ9dxxP7OstyHGf2FME/diQWzhujdHP6dLCwDti/ByE+wIJ5ntPPnsgTX8/OyjQfx4a9RdDFjMTC+a5/3ZfvKgCO52LayBjMnXuBbP0szgz1RdFg6YEHHsDtt9/e6zHJyckwGAyoqKjocrvZbEZNTQ0MBoPT10tNTQUAnD59GiNHjoTBYMCePXu6HFNeXg4AvZ5Xr9dDrz+/6INOp5PlD0au81JX7GfX8aS+HhcbAqAIpytbPKLNpY22fMtwvWf1syfzxn6eHB+GDXuLcLysUZHndrLCVgh2fFyI4/py9LOz51M0WIqKikJUVFSfx6WlpaGurg779+/H9OnTAQDbtm2D1Wp1BEDOyM7OBgDExsY6zvunP/0JFRUVjmm+jIwMBAcHY/z48f18NkTkjcTClCfLPaPWUmG1bVohwtf9pwzJfU10JHk3QBAEl1eCF1fCjYt1jzxBj1gNN27cOCxYsAArVqzAnj178NNPP2HVqlVYsmQJ4uLiAADFxcVISUlxjBSdOXMGTz/9NPbv34/8/Hxs3LgRy5Ytw6WXXopJkyYBAObNm4fx48fjV7/6FQ4dOoRvv/0Wjz32GO65555uR46IaOgRg6Wi2lY0Gd2/BptYdTmSwRINwpiYIPho1KhvNeFcTatLr222WB3lOsYZGCz1ywcffICUlBTMmTMHCxcuxOzZs/HPf/7Tcb/JZEJubq4jWcvHxwffffcd5s2bh5SUFDzwwAO44YYb8OWXXzoeo9Fo8NVXX0Gj0SAtLQ233norli1b1qUuExENbWEBPogKsn15OuUBo0sF9mApgt/3aBB8tGpHnTFX11vKq7JtcxKo12J4mHuUv/CI1XAAEB4e3msBysTExC71IOLj47Fjx44+zztixAhs2rRJkjYSkXcaGxOEykYjTpU3YWpCmNLN6ZU4DceRJRqsCcNCcKS4HkeK63H1pFiXXfe4oxhlENRqFSxuUOLMY0aWiIiUIk7F5br5yFKbyYKyBluNpUgW76ZBmqjQtifiXozusM2JiMESEVEfxsQEAnD/JO+iWtuoUoBegwCPmTcgd9V52xNXVvJ2t+RugMESEVGfxjg21HXvYKnAPgWXEOav+C7t5PnGGAKh06hQ32pCUa3rkrxP2DfQHRfLkSUiIo8xOto2slTRaERdi/vuGymuhIsPd4+kWPJseq3G5Unetc3tjqnksW6yEg5gsERE1KcgXx2G2TelPWlf0uyOHCNLCuzlRd6p81ScK+TYN89NCPdHoN595pIZLBEROUHMW3LnJO9z4siSmyy3Js83wcVJ3sdLbMGSOyV3AwyWiIicIuYtuXOtJbHGEkeWSCquTvLeX1ALAJiSECr7tfqDwRIRkRPGxrh3krfVKjhylhKYs0QSGWsIgk6jQl2L/EnegiBgnz1YmjEiXNZr9ReDJSIiJ3TeI86Vy6idVdFoRLvZCo1ahdgQFlkiaei1GsdrX+6puHM1rahsNEKnUWHS8BBZr9VfDJaIiJwwKjoQKhVQ22JCVZP7rYgrqLbt0j4s1A86Dd/aSTriVNxhmYOlfQU1juv56jSyXqu/+BdFROQEX50GiREBANyzOGUh85VIJtPsW/zszauR9Tp78+1TcInuNQUHMFgiInKaWG/JHfOWHMFSBIMlktZFyREAgENFdWhpN8t2nf32kaXpI9xv/0UGS0REThIL9J2qcL9gSayxNIIjSySx+HA/DAv1g8kiYJ999Edq9S0mRw0zBktERB5sjBuviOM0HMlFpVI5Rpd2namW5RoHCm1BWFJkACID9bJcYzAYLBEROUkcWTpZ3uR2K+I4DUdyShtpC5Yyz8oTLO1z4yk4gMESEZHTEiMCoFWr0GQ0o6S+TenmODS2mVDTbFuhx5ElkoMYLB0trkdjm0ny84vTexcmMlgiIvJoPlo1RkaJSd4NCremgziqFB7ggyBfncKtIW80LNQPIyL8YbEK2Jsv7aq4drMV2efqAADT3awYpYjBEhFRP6TE2qbickrdJ2+p0J7cHc9RJZJRmj1vKVPivKVjJfUwmq0I89dhZFSApOeWCoMlIqJ+SDEEAwBOuFGStziyxJVwJCe58pbE/eCmjwiDSqWS9NxSYbBERNQP4sjSiVL3mYYTN9AdweRukpE4snSspAH1LdLlLYn5Su46BQcwWCIi6pdx9pGls1XNaDNZFG6NzbkaTsOR/KKDfZEcFQBBALLypBld6rJ5rpsmdwMMloiI+iUmWI8wfx0sVgGnK5qUbg4AFqQk13HkLUk0FVdQ3YKqJiN8NGrHHnTuiMESEVE/qFQqR95SjhtMxZksVhTXtQJgjSWSnyNvSaIkb3FUaeJw99s8tzMGS0RE/eTIW3KDJO/SujZYrAJ8tGrEBPkq3RzycmIl7xNljY7aXoMh7gc3w02LUYoYLBER9dM4x4o45UeWCmqaAdiKUarV7rmSiLxHZKAeY2JstcayJJiK60juZrBERORVOtdaUnrbE+4JR64mVd7S2comnKpogkatwoWJ7rsSDmCwRETUb2NigqBWATXN7ahsMiraFrEgJYMlchUxb2mwm+p+kV0CALh0dCTCAnwG3S45MVgiIuonX50GSZG2SsMnFK7kXcBgiVwsNSkCKhVwuqLJUbaivwRBwBcHiwEAi6cOk7J5smCwREQ0ACmx7pG3lF9ty1liQUpylbAAH8dU3Md7zw3oHAcK61BY04IAHw3mjTdI2TxZMFgiIhqAcQbl94izWAWcrbIFS6OiAxVrBw09t6QmAAA+3ncOJou1348XR5XmTzDAz8d9SwaIGCwREQ2AO9RaKqxpQbvZCr1WjeFhHFki15k33oDIQB9UNhqxNae8X49tN1vx1WFbvtJ1HjAFBzBYIiIaEHFF3JnKJrSb+//NWgony22jWqOiA6Fh2QByIR+tGjfOiAcAfJBV2K/H7jxZidoWE6KC9Jg1MlKO5kmOwRIR0QAMC/VDkF4Lk0XA2Spltj0Rt1sZExOkyPVpaLv5QttU3A+nqhyrMp3xebZtCm7R5DiPCfIZLBERDYBKpeqo5K1Q3lLnkSUiV0uI8Mclo20jQx/tdW50qaHNhO+O26btPGEVnIjBEhHRAI2zr4jLUWhF3MlyjiyRspbaE70/3XfOqenozUfLYDRbMTo6EBfEBcvdPMl4TLBUU1ODpUuXIjg4GKGhobjjjjvQ1NTz0Hd+fj5UKlW3P59++qnjuO7u37BhgyueEhF5ODHJW4mRJYtVwJlKMVjiyBIpY864GEQF6VHV1I6M430neneuraRSecYUHOBBwdLSpUtx7NgxZGRk4KuvvsLOnTtx11139Xh8fHw8SktLu/ysXbsWgYGBuOqqq7oc+/bbb3c5bvHixTI/GyLyBh3bnrh+ZIkr4cgd6DRq3GRP9P5wT0Gvx5bWtzq2SFk0JU72tklJq3QDnJGTk4PNmzdj7969mDFjBgDg1VdfxcKFC/Hiiy8iLu78TtdoNDAYuha6+vzzz/HLX/4SgYFdv4WFhoaedywRUV/G2qe/KhqNqG4yIiJQ77JrcyUcuYslM+Oxfvtp/HS6GnlVzY7q9p2ZLVY89eVxCAIwMzHc4wJ8jxhZyszMRGhoqCNQAoD09HSo1WpkZWU5dY79+/cjOzsbd9xxx3n33XPPPYiMjMTMmTPx1ltvKb4xJhF5hgC91lE5O7fMtVNxXAlH7mJ4mD8uHxMFAFjz2WHUNLd3ud9iFfD7Tw/hm6Nl8NGocd/c0Uo0c1A8YmSprKwM0dHRXW7TarUIDw9HWVmZU+d48803MW7cOMyaNavL7U899RSuvPJK+Pv7Y8uWLfjNb36DpqYm/O53v+vxXEajEUZjx+aZDQ22IXiTyQSTyeTs0+qTeC4pz0nnYz+7jjf29ZjoQBRUt+BocR0uHBHisuueKK0HACRH+J3Xn97Yz+6I/dxh5WVJ2H22GrvP1uDqV37Aq0smY/LwEFitAh7beBxfZJdAq1bhrzdNwoUJIf3qMzn72dlzqgQFh1EeeeQRrFu3rtdjcnJy8Nlnn+Hdd99Fbm5ul/uio6Oxdu1arFy5stdztLa2IjY2Fo8//jgeeOCBXo994okn8Pbbb+PcuZ73u/njH/+ItWvXnnf7hx9+CH9/zxpaJKLB+eacCpuLNJgZZcXSUa4rTvnCIQ2KW1S4c6wFE8M5Gk7KK2kB3s7VoKJNBY1KwPWJVpS2qPBjuRoqCLhtjBVTI9zrtdrS0oJbbrkF9fX1CA7ueXWeosFSZWUlqqurez0mOTkZ77//Ph544AHU1tY6bjebzfD19cWnn36K6667rtdz/N///R/uuOMOFBcXIyoqqtdjv/76a/ziF79AW1sb9Pru8w+6G1mKj49HVVVVr53dXyaTCRkZGZg7dy50Op1k56Wu2M+u4419/V1OBVZ+mI2UmEB8uWpW3w+QgMUqYNLTW9FutuK7+2djRHjXL2ne2M/uiP18vsY2Mx75/Ci2HK9w3KZSAf97/YQBJ3XL2c8NDQ2IjIzsM1hSdBouKiqqz+AFANLS0lBXV4f9+/dj+vTpAIBt27bBarUiNTW1z8e/+eabuPbaa526VnZ2NsLCwnoMlABAr9d3e79Op5PlD0au81JX7GfX8aa+np5o2339ZEUT2q0qBOjlf1stqmp2rIRLigruMcHbm/rZnbGfO4TrdPjHr2bgjR/OYt3mXFisAp6/fiL+n73a92DI0c/Ons8jcpbGjRuHBQsWYMWKFXj99ddhMpmwatUqLFmyxLESrri4GHPmzMF7772HmTNnOh57+vRp7Ny5E5s2bTrvvF9++SXKy8tx0UUXwdfXFxkZGXj22Wfx+9//3mXPjYg8W3SwL2JDfFFa34ajxfVITY6Q/ZpcCUfuTKVS4a5LR+KyMdFobjdjWkKY0k0aNI8IlgDggw8+wKpVqzBnzhyo1WrccMMNeOWVVxz3m0wm5ObmoqWl6/40b731FoYPH4558+add06dTof169fj/vvvhyAIGDVqFF5++WWsWLFC9udDRN5j8vBQlNaX4VBRnUuCJa6EI08w1uA9r0+PCZbCw8Px4Ycf9nh/YmJit0v+n332WTz77LPdPmbBggVYsGCBZG0koqFpcnwoNh8rw6Fz9S65HveEI3Itj6izRETkzibH20oGZJ+rc8n1TnFPOCKXYrBERDRIk4aHQqUCiutaUdHYJuu1uCcckesxWCIiGqRAvRaj7VNih2WeiiusaYGRe8IRuRSDJSIiCUweHgoAOFRUJ+t1TnElHJHLMVgiIpLA5PhQAPLnLZ3iSjgil2OwREQkgSn2YOnQuTpYrfJtjMCVcESux2CJiEgCYw1B0GvVaGgzI7+6WbbrcCUckesxWCIikoBOo8aEYbYSAnLlLXElHJEyGCwREUlETPLOLqyT5fxcCUekDAZLREQScRSnLJKnfEBuWQMAroQjcjUGS0REEhGTvHNKGmA0WyQ//778WgAdK++IyDUYLBERSSQh3B9h/jq0W6w4Udoo+fn3FtiCpZmJ4ZKfm4h6xmCJiEgiKpVKtnpLLe1mHCu2Te/NSAyT9NxE1DsGS0REEnJU8pY4WMourIPZKiAuxJfJ3UQuxmCJiEhCYt5StsTlA/ba85VmcAqOyOUYLBERSWjScNuKuLOVzahvNUl23r35NQCAC5MYLBG5GoMlIiIJRQTqMSLCNk22N69GknOaLVYcKLSNLF3IfCUil2OwREQksUtGRwIAvs+tkOR8x0sb0NJuQbCvFmOiuc0JkasxWCIiktiVKdEAgO9PVEAQBr+pbud8JTWLURK5HIMlIiKJpSVHQq9Vo6S+DSftG98Ohjidx5IBRMpgsEREJDE/Hw3SRkYAALadGNxUnCAI2FdgC5ZYjJJIGQyWiIhk0HkqbjDyqppR1dQOH60aE+0r7YjItRgsERHJ4IqxtmBpf2Et6lsGXkJA3A9uyvBQ6LUaSdpGRP3DYImISAbx4f4YFR0Ii1XAzlOVAz6PWF+J+UpEymGwREQkEymm4liMkkh5DJaIiGQiTsVtP1kJi7X/JQQqGtuQX90ClQqYlsCRJSKlMFgiIpLJjMQwBOm1qGlux+EB7BUn5iulGIIR4qeTuHVE5CwGS0REMtFp1LhkjL2a9wCm4hxTcMxXIlIUgyUiIhmJU3Hb+rn1iSAI+OFUFQBb5W4iUg6DJSIiGV02NgoAcLS4ARUNbU4/7qfT1Thd0QR/Hw0uGx0lV/OIyAkMloiIZBQd5ItJ9mKS23OdLyHwrx/PAgB+OSMeIf7MVyJSEoMlIiKZiVNxXx4ucer4U+WN2J5bCZUKWH5xoowtIyJnMFgiIpLZoilx0KpV+OFUFX605yH15s0f8wAA88bHYEREgNzNI6I+MFgiIpJZclQgbr1oBADg6a+Ow2yx9nhsVZMRnx0sBgDceUmyS9pHRL1jsERE5AL3pY9GqL8OueWN+GjvuR6Pe393AdrNVkyOD8WMESwZQOQOGCwREblAqL8P7k8fAwB4eUsu6lvP31y3zWTB/2UWAADunJ0ElUrl0jYSUfc8Jlj605/+hFmzZsHf3x+hoaFOPUYQBDzxxBOIjY2Fn58f0tPTcerUqS7H1NTUYOnSpQgODkZoaCjuuOMONDU1yfAMiGioW5qagNHRgahtMeHVrafOu/+/2cWobm7HsFA/XDXBoEALiag7HhMstbe348Ybb8TKlSudfswLL7yAV155Ba+//jqysrIQEBCA+fPno62to9bJ0qVLcezYMWRkZOCrr77Czp07cdddd8nxFIhoiNNq1HjsF+MBAO/sysfZyo4vZmaLFf/6wZbYffusRGg1HvP2TOT1tEo3wFlr164FALzzzjtOHS8IAv7yl7/gsccew6JFiwAA7733HmJiYvDFF19gyZIlyMnJwebNm7F3717MmDEDAPDqq69i4cKFePHFFxEXFyfLcyGioeuyMVG4MiUa205U4MF/H8ZYQxCOlTTgRGkDjGYrAnw0uGlmvNLNJKJOPCZY6q+8vDyUlZUhPT3dcVtISAhSU1ORmZmJJUuWIDMzE6GhoY5ACQDS09OhVquRlZWF6667rttzG41GGI1Gx+/19fUAbFN6JtP5eQgDZTKZ0NLSgurqauh0LEonF/az67Cvbe5Ji8H2IwXYe7IFe0923O6vV2PlrJEwNTegunng52c/uwb72TXk7OfGxkYAtgGW3nhtsFRWVgYAiImJ6XJ7TEyM476ysjJER0d3uV+r1SI8PNxxTHeee+45x0hXZ0lJSYNtNhENcb9dB/xW6UYQDTGNjY0ICQnp8X5Fg6VHHnkE69at6/WYnJwcpKSkuKhFzlmzZg1Wr17t+N1qtaKmpgYRERGSrl5paGhAfHw8zp07h+DgYMnOS12xn12Hfe0a7GfXYD+7hpz9LAgCGhsb+0y7UTRYeuCBB3D77bf3ekxy8sCKshkMtpUk5eXliI2NddxeXl6OKVOmOI6pqOi6E7jZbEZNTY3j8d3R6/XQ6/VdbnN2hd5ABAcH8w/RBdjPrsO+dg32s2uwn11Drn7ubURJpGiwFBUVhagoeXbTTkpKgsFgwNatWx3BUUNDA7Kyshwr6tLS0lBXV4f9+/dj+vTpAIBt27bBarUiNTVVlnYRERGRZ/GYtamFhYXIzs5GYWEhLBYLsrOzkZ2d3aUmUkpKCj7//HMAgEqlwn333YdnnnkGGzduxJEjR7Bs2TLExcVh8eLFAIBx48ZhwYIFWLFiBfbs2YOffvoJq1atwpIlS7gSjoiIiAB4UIL3E088gXfffdfx+9SpUwEA33//PS6//HIAQG5urmNlGgA89NBDaG5uxl133YW6ujrMnj0bmzdvhq+vr+OYDz74AKtWrcKcOXOgVqtxww034JVXXnHNk+qDXq/Hk08+ed6UH0mL/ew67GvXYD+7BvvZNdyhn1VCX+vliIiIiIYwj5mGIyIiIlICgyUiIiKiXjBYIiIiIuoFgyUiIiKiXjBYcmPr169HYmIifH19kZqaij179ijdJI/23HPP4cILL0RQUBCio6OxePFi5Obmdjmmra0N99xzDyIiIhAYGIgbbrgB5eXlCrXY8z3//POOMh4i9rF0iouLceuttyIiIgJ+fn6YOHEi9u3b57hfEAQ88cQTiI2NhZ+fH9LT03Hq1CkFW+x5LBYLHn/8cSQlJcHPzw8jR47E008/3WUvMfZz/+3cuRPXXHMN4uLioFKp8MUXX3S535k+rampwdKlSxEcHIzQ0FDccccdXcoJSYnBkpv6+OOPsXr1ajz55JM4cOAAJk+ejPnz559XcZyct2PHDtxzzz3YvXs3MjIyYDKZMG/ePDQ3d+xYev/99+PLL7/Ep59+ih07dqCkpATXX3+9gq32XHv37sU//vEPTJo0qcvt7GNp1NbW4uKLL4ZOp8M333yD48eP46WXXkJYWJjjmBdeeAGvvPIKXn/9dWRlZSEgIADz589HW1ubgi33LOvWrcNrr72Gv/3tb8jJycG6devwwgsv4NVXX3Ucw37uv+bmZkyePBnr16/v9n5n+nTp0qU4duwYMjIy8NVXX2Hnzp2466675GmwQG5p5syZwj333OP43WKxCHFxccJzzz2nYKu8S0VFhQBA2LFjhyAIglBXVyfodDrh008/dRyTk5MjABAyMzOVaqZHamxsFEaPHi1kZGQIl112mXDvvfcKgsA+ltLDDz8szJ49u8f7rVarYDAYhP/93/913FZXVyfo9Xrho48+ckUTvcLVV18t/PrXv+5y2/XXXy8sXbpUEAT2sxQACJ9//rnjd2f69Pjx4wIAYe/evY5jvvnmG0GlUgnFxcWSt5EjS26ovb0d+/fvR3p6uuM2tVqN9PR0ZGZmKtgy7yIWMA0PDwcA7N+/HyaTqUu/p6SkICEhgf3eT/fccw+uvvrqLn0JsI+ltHHjRsyYMQM33ngjoqOjMXXqVLzxxhuO+/Py8lBWVtalr0NCQpCamsq+7odZs2Zh69atOHnyJADg0KFD+PHHH3HVVVcBYD/LwZk+zczMRGhoKGbMmOE4Jj09HWq1GllZWZK3yWMqeA8lVVVVsFgsiImJ6XJ7TEwMTpw4oVCrvIvVasV9992Hiy++GBMmTAAAlJWVwcfH57xNkWNiYlBWVqZAKz3Thg0bcODAAezdu/e8+9jH0jl79ixee+01rF69Go8++ij27t2L3/3ud/Dx8cFtt93m6M/u3kfY18575JFH0NDQgJSUFGg0GlgsFvzpT3/C0qVLAYD9LANn+rSsrAzR0dFd7tdqtQgPD5el3xks0ZB0zz334OjRo/jxxx+VbopXOXfuHO69915kZGR02VaIpGe1WjFjxgw8++yzAGxbQB09ehSvv/46brvtNoVb5z0++eQTfPDBB/jwww9xwQUXIDs7G/fddx/i4uLYz0MIp+HcUGRkJDQazXkrhMrLy2EwGBRqlfdYtWoVvvrqK3z//fcYPny443aDwYD29nbU1dV1OZ797rz9+/ejoqIC06ZNg1arhVarxY4dO/DKK69Aq9UiJiaGfSyR2NhYjB8/vstt48aNQ2FhIQA4+pPvI4Pz4IMP4pFHHsGSJUswceJE/OpXv8L999+P5557DgD7WQ7O9KnBYDhvwZPZbEZNTY0s/c5gyQ35+Phg+vTp2Lp1q+M2q9WKrVu3Ii0tTcGWeTZBELBq1Sp8/vnn2LZtG5KSkrrcP336dOh0ui79npubi8LCQva7k+bMmYMjR44gOzvb8TNjxgwsXbrU8f/sY2lcfPHF55W+OHnyJEaMGAEASEpKgsFg6NLXDQ0NyMrKYl/3Q0tLC9Tqrh+VGo0GVqsVAPtZDs70aVpaGurq6rB//37HMdu2bYPVakVqaqr0jZI8ZZwksWHDBkGv1wvvvPOOcPz4ceGuu+4SQkNDhbKyMqWb5rFWrlwphISECNu3bxdKS0sdPy0tLY5j7r77biEhIUHYtm2bsG/fPiEtLU1IS0tTsNWer/NqOEFgH0tlz549glarFf70pz8Jp06dEj744APB399feP/99x3HPP/880JoaKjw3//+Vzh8+LCwaNEiISkpSWhtbVWw5Z7ltttuE4YNGyZ89dVXQl5envDZZ58JkZGRwkMPPeQ4hv3cf42NjcLBgweFgwcPCgCEl19+WTh48KBQUFAgCIJzfbpgwQJh6tSpQlZWlvDjjz8Ko0ePFm6++WZZ2stgyY29+uqrQkJCguDj4yPMnDlT2L17t9JN8mgAuv15++23Hce0trYKv/nNb4SwsDDB399fuO6664TS0lLlGu0Ffh4ssY+l8+WXXwoTJkwQ9Hq9kJKSIvzzn//scr/VahUef/xxISYmRtDr9cKcOXOE3NxchVrrmRoaGoR7771XSEhIEHx9fYXk5GThD3/4g2A0Gh3HsJ/77/vvv+/2/fi2224TBMG5Pq2urhZuvvlmITAwUAgODhaWL18uNDY2ytJelSB0KkNKRERERF0wZ4mIiIioFwyWiIiIiHrBYImIiIioFwyWiIiIiHrBYImIiIioFwyWiIiIiHrBYImIiIioFwyWiIiIiHrBYImIvE5lZSVWrlyJhIQE6PV6GAwGzJ8/Hz/99JPSTSMiD6RVugFERFK74YYb0N7ejnfffRfJyckoLy/H1q1bUV1drXTTiMgDcWSJiLxKXV0dfvjhB6xbtw5XXHEFRowYgZkzZ2LNmjW49tprAQAqlQqvvfYarrrqKvj5+SE5ORn//ve/u5zn4YcfxpgxY+Dv74/k5GQ8/vjjMJlMXY758ssvceGFF8LX1xeRkZG47rrrHPcZjUb8/ve/x7BhwxAQEIDU1FRs375d9udPRNJjsEREXiUwMBCBgYH44osvYDQaezzu8ccfxw033IBDhw5h6dKlWLJkCXJychz3BwUF4Z133sHx48fx17/+FW+88Qb+/Oc/O+7/+uuvcd1112HhwoU4ePAgtm7dipkzZzruX7VqFTIzM7FhwwYcPnwYN954IxYsWIBTp07J88SJSDbcSJeIvM5//vMfrFixAq2trZg2bRouu+wyLFmyBJMmTQJgG1m6++678dprrzkec9FFF2HatGn4+9//3u05X3zxRWzYsAH79u0DAMyaNQvJycl4//33zzu2sLAQycnJKCwsRFxcnOP29PR0zJw5E88++6yUT5eIZMaRJSLyOjfccANKSkqwceNGLFiwANu3b8e0adPwzjvvOI5JS0vr8pi0tLQuI0sff/wxLr74YhgMBgQGBuKxxx5DYWGh4/7s7GzMmTOn2+sfOXIEFosFY8aMcYx0BQYGYseOHThz5oy0T5aIZMcEbyLySr6+vpg7dy7mzp2Lxx9/HHfeeSeefPJJ3H777X0+NjMzE0uXLsXatWsxf/58hISEYMOGDXjppZccx/j5+fX4+KamJmg0Guzfvx8ajabLfYGBgQN+TkSkDI4sEdGQMH78eDQ3Nzt+3717d5f7d+/ejXHjxgEAdu3ahREjRuAPf/gDZsyYgdGjR6OgoKDL8ZMmTcLWrVu7vdbUqVNhsVhQUVGBUaNGdfkxGAwSPzMikhtHlojIq1RXV+PGG2/Er3/9a0yaNAlBQUHYt28fXnjhBSxatMhx3KeffooZM2Zg9uzZ+OCDD7Bnzx68+eabAIDRo0ejsLAQGzZswIUXXoivv/4an3/+eZfrPPnkk5gzZw5GjhyJJUuWwGw2Y9OmTY5VdEuXLsWyZcvw0ksvYerUqaisrMTWrVsxadIkXH311S7tEyIaJIGIyIu0tbUJjzzyiDBt2jQhJCRE8Pf3F8aOHSs89thjQktLiyAIggBAWL9+vTB37lxBr9cLiYmJwscff9zlPA8++KAQEREhBAYGCjfddJPw5z//WQgJCelyzH/+8x9hypQpgo+PjxAZGSlcf/31jvva29uFJ554QkhMTBR0Op0QGxsrXHfddcLhw4dl7wMikhZXwxHRkKNSqfD5559j8eLFSjeFiDwAc5aIiIiIesFgiYiIiKgXTPAmoiGH2QdE1B8cWSIiIiLqBYMlIiIiol4wWCIiIiLqBYMlIiIiol4wWCIiIiLqBYMlIiIiol4wWCIiIiLqBYMlIiIiol4wWCIiIiLqxf8HWuQF3iMmmKEAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "NUM_SPATIAL_DIMS = 1\n", + "DOMAIN_EXTENT = 3.0\n", + "NUM_POINTS = 100\n", + "DT = 0.1\n", + "\n", + "burgers_stepper = ex.stepper.Burgers(NUM_SPATIAL_DIMS, DOMAIN_EXTENT, NUM_POINTS, DT)\n", + "\n", + "u_0 = ex.ic.RandomTruncatedFourierSeries(\n", + " NUM_SPATIAL_DIMS,\n", + " cutoff=5,\n", + " max_one=True,\n", + ")(NUM_POINTS, key=jax.random.PRNGKey(0))\n", + "\n", + "ex.viz.plot_state_1d(u_0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The state of the initial condition is a 1x100 tensor with **real** floating\n", + "point values" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "((1, 100), dtype('float32'))" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "u_0.shape, u_0.dtype" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For an integration in Fourier space, we have to transform it to Fourier space.\n", + "Important, whenever we are dealing with FFTs in `Exponax`, we need to use\n", + "`rfft`." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "u_0_hat = jnp.fft.rfft(u_0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Its shape is `(1, 51)`, and has complex values. (JAX' `complex64` type is composed of two `float32` values)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "((1, 51), dtype('complex64'))" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "u_0_hat.shape, u_0_hat.dtype" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can use the familiar `ex.rollout` function transformation but need to\n", + "transform the step function in Fourier space." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "trj_hat = ex.rollout(burgers_stepper.step_fourier, 100, include_init=True)(u_0_hat)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Using a `jnp.fft.irfft` will batch over all time steps. (We need to inform the\n", + "number of points because we used the real-valued FFT)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "trj = jnp.fft.irfft(trj_hat, n=NUM_POINTS)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAG2CAYAAACeUpnVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/SrBM8AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA7dklEQVR4nO3df3hU1Z3H8c/MZDKTkB+AkQQwFNZi1RUBQdJou9o2K2tdutauS9EKi7Y+dSkrZrsCLT9qbY3a1WVb2LKldds+Wwu1P2wrSpdGxdqiaAB/VbFaWyglAUQyECA/Zu7+0e2EzD3IZWYyd2bO+/U893k6h3NnztzbJF/P93vuCTiO4wgAAMBiQb8HAAAA4DcCIgAAYD0CIgAAYD0CIgAAYD0CIgAAYD0CIgAAYD0CIgAAYD0CIgAAYD0CIgAAYD0CIgAAYD1fA6InnnhCM2bM0KhRoxQIBPTggw+e9JzHH39cF1xwgSKRiN75znfqm9/85qCPEwAAFDdfA6Kuri5NnDhRq1at8tT/jTfe0BVXXKH3ve992r59uxYsWKCPf/zj+tnPfjbIIwUAAMUskC+buwYCAf3oRz/SlVdeecI+Cxcu1Pr16/Xiiy8m2z760Y/q4MGD2rBhQw5GCQAAilGJ3wM4FZs3b1ZTU9OAtunTp2vBggUnPKe7u1vd3d3J14lEQgcOHNBpp52mQCAwWEMFAABZ5DiODh06pFGjRikYzH6Cq6ACovb2dtXW1g5oq62tVSwW09GjR1VWVuY6p6WlRbfddluuhggAAAbRrl27dMYZZ2T9fQsqIErH4sWL1dzcnHzd2dmpMWPG6F+/u0mR8opk+8XvGOY6tyrqvjzhoHtWqcQQqIYM/UIpM1IhwwRVah9JMk1kmc41zXiZzvUaV5vez/C1vL1XeqedULrjAAAUpkOHDmn8+PGqrKwclPcvqICorq5OHR0dA9o6OjpUVVVlnB2SpEgkokgk4m4vr1B0SH9ANMRwgSuiYVebKSAKpxsQmc7LICAKEhABAIrcYJW7FNRziBobG9Xa2jqgbePGjWpsbPRpRAAAoBj4OkN0+PBhvfbaa8nXb7zxhrZv367hw4drzJgxWrx4sXbv3q1vf/vbkqRPfvKTWrlypW699VZdf/31evTRR/W9731P69evP+XP3vq7txQu60m+Pq/OPUMUNhRtDSkNudqMMzOGtXuBlMaAYzhP7hNNs0ampYEJw4LBoGFuJpHBrFHC8MFeZmtM480kxjeNwytmlwAAqXydIXr22Wc1efJkTZ48WZLU3NysyZMna9myZZKkPXv2aOfOncn+48aN0/r167Vx40ZNnDhR99xzj77+9a9r+vTpvowfAAAUB19niC699FK93WOQTE+hvvTSS7Vt27ZBHBUAALBNQRVVZ1N7e0yhaDz5+q2jva4+1Yai6kjCPalmLI42xHlOSorMFAoasmgyxoxZTvskDG1BwwebitmymUYzyXaGy2u6jdQaANijoIqqAQAABgMBEQAAsB4BEQAAsJ61NURHDvUo1NP/9Q/3xF19uvvcbb2GpzCGHfdSfFNNTupSfNPDFb3yWrfjdSm+1+dcmYrgvdQVZVKPk4MSKiOW9gOAPZghAgAA1iMgAgAA1iMgAgAA1rO2hqivNyEn0F/pc7i7z9XnWJ+7EihhKA4yPSfIy7ODEobqmEy288g247OJ0n2vNJ9VdCJ+PcPIq0zqj1JRjwQAg48ZIgAAYD0CIgAAYD1rU2bBUEDBkv5cxBHDsntT2iNuyIWZd57PZHSFxetS/FTZTqOZ+LVkP5tY/g8Ag48ZIgAAYD0CIgAAYD0CIgAAYD1ra4iGVEUUikTfto+pXsi0FYZX7lqb9As8vCzrz/T9MlnZn/pdvdQUSf7VFZkUQ/lNtmvZqEkCUKyYIQIAANYjIAIAANYjIAIAANaztoZo2GlDVFI2JPk6UuKODbO9PYbXOprBZqqDCnqt8TG0DXZUnYu6IhObao28SrcmidojAPmOGSIAAGA9AiIAAGA9a1Nm579jqCLlFcnXFVH3pTBN85tSS6ZsgJcUQdBwptfUQp5k3zxLd3uPE/ErjWbiJYtUYLcr67ym2kitAfALM0QAAMB6BEQAAMB6BEQAAMB61tYQTRhZpfKKyuTr3oR7QXk45G0pvqkUxtSWbn2E1yXxfkl3KX4u6opM/KhTYQm/N5lsNUL9EYBMMEMEAACsR0AEAACsR0AEAACsZ20N0TmnD1FFZf9ziHbFjrn6lIdDrrZwyNuzg0yRZmp5jLH2yHBetuV7TdJgy6dnGKWi1ih9+VxDBiD/MUMEAACsR0AEAACsZ23KrG5IiSor+r9+T7zU1SduWBZe6nEpfsgwL5+6VYfnrUE8LuEvhkxAtpfie5XPaTQTU3Yoj4ebV0itATBhhggAAFiPgAgAAFiPgAgAAFjP2hqiiniXKuP98WB1pMzV51jcXWxgqhcqMRQbmMpeUsuPTCUKlq+IN8qnuiKTfKk1Ycl+dhVaXRmAzDBDBAAArEdABAAArEdABAAArGdtDVGo602Fgt3J1xWVY92deuKuJlPtSthQWGB4XJGrdsP4rCKPNQpeSxn82qYjkTqOLL+/X3VFJoVWa8IzjNLnpa4sn+89gBNjhggAAFiPgAgAAFjP2pRZsPuwgqX989/lw92xoWnrDtOUeSjNLTgy2X4j26mwYljun89ptHxPo7BkP3sKLYUK4E+YIQIAANYjIAIAANYjIAIAANaztoYo0XNUie7+eLA0kLpQXIoYioOyWWuRiyXxxVAblIl8qSsqlroSluynp1juP1DMmCECAADWIyACAADWIyACAADWs7aGKBAMKRAM9b/uPebqEw2Xu9riXp7dfwLuKiXYrFjqSqgrSk+x3H+gWDBDBAAArEdABAAArGdtykylUSlS1v860efqUuK420Ih9yUzTX2b0gihlNdes2+mpeNekaZzy5el+CbFkkZJ/RoF+BV8USz3HyhEzBABAADrERABAADrERABAADrWVtDlIhUKhGpSL4OJOLuToal+AFDDVEgaLiMhjYvdRWmaqGExwoMU21MJhFvNuuPTO+VT9E4dUXIV9x/IDfy6W8SAACALwiIAACA9QiIAACA9XwPiFatWqWxY8cqGo2qoaFBW7Zsedv+K1as0Lve9S6VlZWpvr5et9xyi44dc9f6nIwTrZQTrUoeSvS5jkBft/voNRx9Pa5DcfcRcBIDjmBAriMg92HqZxIIBFxHJoKGwyaO47gOpMcxHACQT3z9G7du3To1Nzdr+fLl2rp1qyZOnKjp06dr7969xv7333+/Fi1apOXLl+vll1/WN77xDa1bt06f+cxncjxyAABQTHwNiO6991594hOf0Ny5c3Xuuedq9erVKi8v13333Wfs/6tf/UoXX3yxrrnmGo0dO1aXXXaZZs2addJZJQAAgLfjW0DU09OjtrY2NTU19Q8mGFRTU5M2b95sPOeiiy5SW1tbMgD67W9/q4cfflgf/OAHT/g53d3disViAw5JSkQqlIhWJg85CdcR6DvmOoypsHiv4ehzHa6UXEoKjTQavEo47qPQkEZLXzHcfyDf+PYcov379ysej6u2tnZAe21trV555RXjOddcc43279+v97znPXIcR319ffrkJz/5timzlpYW3XbbbVkdOwAAKC4F9R/4jz/+uO644w7953/+p7Zu3aof/vCHWr9+vW6//fYTnrN48WJ1dnYmj127duVwxAAAoBD4NkNUU1OjUCikjo6OAe0dHR2qq6sznrN06VJdd911+vjHPy5JmjBhgrq6unTjjTfqs5/9rIJBd3wXiUQUiUSy/wUAAEDR8G2GqLS0VFOmTFFra2uyLZFIqLW1VY2NjcZzjhw54gp6QqGQJPPWC2+nRyUDjkCiz3309bqP3mOuw1RXZFzGn1IvZKxbkrd6oUzqirIt3bqihOHIdyzFB4Di5OteZs3NzZozZ46mTp2qadOmacWKFerq6tLcuXMlSbNnz9bo0aPV0tIiSZoxY4buvfdeTZ48WQ0NDXrttde0dOlSzZgxIxkYAQAAnCpfA6KZM2dq3759WrZsmdrb2zVp0iRt2LAhWWi9c+fOATNCS5YsUSAQ0JIlS7R7926dfvrpmjFjhr74xS/69RUAAEARCDiWzfnHYjFVV1dr1x/bVVVVlWwvO7Tb1TfQ1+tqc4LumSgnHHW3lRjqlkpKU97LEI8a2kw3yLTM1nO/HNzydNNfBVXl//8yfZRBthTDDuhF8BV8Uwz3H3g7sVhMdXV16uzsHPD3O1t8nSHyU9xxFD8uMDAFP4E+w5YgJe7gR4k+Q5shsImn9AsY/vw7hlDC1C8Dpj/glsXFRckU/PJHEgC8KcT/IAcAAMgqAiIAAGA9a1NmruXh8R53J1ObKX1lqCtSyJD6SkmHBQzpMcfQFshyyiyfmWqP7Pn2MCVuyfp5Q8oUyAx/awAAgPUIiAAAgPUIiAAAgPWsrSGKhAKKhPoT7MYl9n3u5fSBoLvNVPdjXD6fRabaAK81BKZ+2ZYaaRfCthzpMj2yIF+eTQQA8IYZIgAAYD0CIgAAYD0CIgAAYD1ra4iCR2MKhvtrP5wjhz2dFwiXutsMNSTG5wml9HMShsoan0LUbG/nUcw1QwCA4sMMEQAAsB4BEQAAsJ69KbPDexUMHEm+Thx6y91nSFUuh4QTYDsPAMBg4+8KAACwHgERAACwHgERAACwnrU1RImO3ylxuDz52jl2xNXHKRviaivWDRkyWWIPAEChY4YIAABYj4AIAABYj4AIAABYz9oaot4/vKbe8mjydbDc/cyh4LDTPb2XY9j2QgF3rOnqF8xuPEoVEGCvYLEWOAI5wgwRAACwHgERAACwnrUps0Ov/0GK9u9cX/2XZ3k70ZAK89yWRQmP+TGv/TJh+872AVPKFGnhSgLwCzNEAADAegREAADAegREAADAetbWEB3evVcqDSdfV5/zTlefQEmpq80JGi6Zafm8h7oix2PtEcvp4QXLrgEgfcwQAQAA6xEQAQAA6xEQAQAA61lbQ9S197AC4f6vHygJu/oESqOuNifkra7ICYbcHxpK6ZfB84tMdUVenznkOOlXJfnxzCGi9uJEyVP6qBcDso+/NQAAwHoERAAAwHrWpsz6jvWpL97/OlA2xN0pbEiPGVJrMqTRzLvdB0/ex/1OGW2/UWjpsXyXL9t0kDIBgOxihggAAFiPgAgAAFiPgAgAAFjP2hqikmiJSo5bdh+sGOrq45R4XHYfctcVOYZtP5SyPN+0dYepXshrXRH1QshXlDylj3oxIDeYIQIAANYjIAIAANYjIAIAANaztoaooq5KFaX9tT+h6tNcfRKlZa42pyRiaDt5vZDkrhmiXsgtnyJ0njmUvgIccl4oxHsNFIt8+vsDAADgCwIiAABgPWtTZsPOqldlWX/6yxkyzNXHKa1wt4XdS/FNW3fEjWmuga8zSYXlc9qrEJEeS18BDjkvFOK9BooZM0QAAMB6BEQAAMB6BEQAAMB61tYQDTlvoiqGlCdfJ4YYlt1HhrjaHMNyelO9kJf6oGzXAXldiZ8n5TLWK8QakgIccl4oxHsN2IYZIgAAYD0CIgAAYD0CIgAAYD1ra4iC77xAwcrK5Ou+cvdziHoNRT5xQ8FQBrto+CKT8RZD/ZEfzxwqxBqSAhxyXijEew2AGSIAAAACIgAAAGtTZvtKa3UsUpV8HTgWd/UJekytpJuByWRqPZNZea/fa7DlIhr3a0uOQkubFNhwfVNo9xWAd8wQAQAA6xEQAQAA6xEQAQAA61lbQ9S257DKD/UXBLyrptzVJxJyx4uGJoUMdSqmWoPUJsdwnteSF6+1DKZuCcO6e691RaYl+3lSkkS9kEEeDy3v5fN9BZB9vs8QrVq1SmPHjlU0GlVDQ4O2bNnytv0PHjyoefPmaeTIkYpEIjrrrLP08MMP52i0AACgGPk6Q7Ru3To1Nzdr9erVamho0IoVKzR9+nTt2LFDI0aMcPXv6enRX//1X2vEiBH6/ve/r9GjR+v3v/+9hg4dmvvBAwCAouFrQHTvvffqE5/4hObOnStJWr16tdavX6/77rtPixYtcvW/7777dODAAf3qV79SOByWJI0dOzaXQwYAAEXIt4Cop6dHbW1tWrx4cbItGAyqqalJmzdvNp7zk5/8RI2NjZo3b55+/OMf6/TTT9c111yjhQsXKhQKGc/p7u5Wd3d38nUsFpMkbXy5Q6XlXcn26gvOcJ17+pBSV1ukxF1YEDLUGnipNXLkLsgJGGp0TPU9hl1FjDUPpl06iqE0gnohtzweWl7J53sIwD++1RDt379f8XhctbW1A9pra2vV3t5uPOe3v/2tvv/97ysej+vhhx/W0qVLdc899+gLX/jCCT+npaVF1dXVyaO+vj6r3wMAABQ+34uqT0UikdCIESP0ta99TVOmTNHMmTP12c9+VqtXrz7hOYsXL1ZnZ2fy2LVrVw5HDAAACoFvKbOamhqFQiF1dHQMaO/o6FBdXZ3xnJEjRyocDg9Ij51zzjlqb29XT0+PSksNKa5IRJFIxNX+3Cv7FIocSb6+5KzTXX3ChlzY0GjY1VZqyNYZd5QPpfYxpN8MIapxmbwhQeI1jWaSyVJ8L/I98s7nNEoeDy2v5PM9BJD/fPs7VVpaqilTpqi1tTXZlkgk1NraqsbGRuM5F198sV577TUlEv1/+l999VWNHDnSGAwBAAB44et/uDc3N2vNmjX61re+pZdfflk33XSTurq6kqvOZs+ePaDo+qabbtKBAwd0880369VXX9X69et1xx13aN68eX59BQAAUAR8XXY/c+ZM7du3T8uWLVN7e7smTZqkDRs2JAutd+7cqWCwP2arr6/Xz372M91yyy06//zzNXr0aN18881auHChX18BAAAUgYDjGKtdilYsFlN1dbXqrv6yguGyZPuCGy919T3n9ApXW025OzU3xFBEVGqoP0ptKjH28bYNiKmfqeQnky0+vNYQeemWyVRktpfY50utSZ4MI6/ky70BkH9isZjq6urU2dmpqqqqrL9/vte6AgAADDoCIgAAYD0CIgAAYD1fi6r9dGj3awqU9D+f6A8Hprn6jKqKutrKw+56IdPzigIBQ6yZUiARMD44yNtmGwHDphzZfjaRiddynnQj7WKpF6IUZiBqgwDkO2aIAACA9QiIAACA9axNmXUfelOBUP8S+j0Hj7n6HO7pc7Ud6XVfskiJO64MG/bgSKRkuUxPPDBt52HKvxTDsxIKMT1me+aH1BeAYsUMEQAAsB4BEQAAsB4BEQAAsJ61NUSpDnT1uNpix9w1RMOicVebaSl+tMS94D2YshQ/ZKhIcb+TZNxcxfJajmzXshTr5aTmBwC8SXuG6PXXX9eSJUs0a9Ys7d27V5L0yCOP6KWXXsra4AAAAHIhrYBo06ZNmjBhgp5++mn98Ic/1OHDhyVJzz33nJYvX57VAQIAAAy2tAKiRYsW6Qtf+II2btyo0tL+pevvf//79dRTT2VtcAAAALmQVg3RCy+8oPvvv9/VPmLECO3fvz/jQeVCuKxywNYd3Ud7XX0OGWqIuuPu2qBjhrbyuLvwpzSlQMhUG5T6rCJJMuwMkle8RNVenznEs4TcqAMCgMGX1gzR0KFDtWfPHlf7tm3bNHr06IwHBQAAkEtpBUQf/ehHtXDhQrW3tysQCCiRSOiXv/ylPv3pT2v27NnZHiMAAMCgSitldscdd2jevHmqr69XPB7Xueeeq3g8rmuuuUZLlizJ9hgHxZARYxUM9+9m39fjXk5/1LB1x7E+d3qs15Ay602425yURfWm9JhpOw/TjvWmJftemc4MGlJameys4SVFVsxL50lzAUBhSSsgKi0t1Zo1a7Rs2TK98MILOnz4sCZPnqzx48dne3wAAACDLqMHM9bX16u+vj5bYwEAAPBFWjVEH/nIR3TXXXe52u+++25dffXVGQ8KAAAgl9KaIXriiSf0uc99ztV++eWX65577sl0TDkxdPQYhUrL37ZPj7FeyFDjY6gFMnRTPKVjuAALTfzY/M6vq1SAtwcAkKa0/r4dPnx4wAMZ/ywcDisWi2U8KAAAgFxKKyCaMGGC1q1b52pfu3atzj333IwHBQAAkEtppcyWLl2qq666Sq+//rre//73S5JaW1v13e9+Vw888EBWBwgAADDY0gqIZsyYoQcffFB33HGHvv/976usrEznn3++fv7zn+uSSy7J9hgHRe2YapVEhyRfxw3PEkqt+ZGkhOE5QXHTHhxpMj2/J5hBFU1OtsLI4jOHclG2Q20QACBV2svur7jiCl1xxRXZHAsAAIAv/Fg0BAAAkFfSmiGKx+P693//d33ve9/Tzp071dPTM+DfDxw4kJXBDaaJ44YpUl6RfP3SH9yr40J5nFsxjczrVhtet+kwRcvZ3LU+k6ubx7emaHCJAeSTwf6dlNYM0W233aZ7771XM2fOVGdnp5qbm3XVVVcpGAwan08EAACQz9IKiL7zne9ozZo1+pd/+ReVlJRo1qxZ+vrXv65ly5bpqaeeyvYYAQAABlVaAVF7e7smTJggSaqoqFBnZ6ck6W//9m+1fv367I0OAAAgB9KqITrjjDO0Z88ejRkzRmeeeab+93//VxdccIGeeeYZRSKRbI9xUEw6Y6jKKyqTrw8e6XX1KS1xx4um+puQsc39mak1SaY6GFOb99ogw7mGfvlcL5TvtUF5PjwAQJrSmiH68Ic/rNbWVknS/PnztXTpUo0fP16zZ8/W9ddfn9UBAgAADLa0ZojuvPPO5P+eOXOmxowZo82bN2v8+PGaMWNG1gYHAACQC2k/mPF4jY2NamxszMZbAQAA5FzaAdGOHTv0la98RS+//LIk6ZxzztH8+fP1rne9K2uDG0xnDi9TRWV58vWO4eWuPqbnEEUNdUVhQz/js35SX3usFzLXGnmt78lu1YvXOqV03ysT1PcAANKVVg3RD37wA5133nlqa2vTxIkTNXHiRG3dulXnnXeefvCDH2R7jAAAAIMqrRmiW2+9VYsXL9bnP//5Ae3Lly/Xrbfeqo985CNZGRwAAEAupBUQ7dmzR7Nnz3a1f+xjH9OXvvSljAeVC6cPCatySDj5un5YmavPkd64qy1iSpkZ1tibUmuhlCbzcv2Tp9pM73XCczNYYp9JeizddBhpLwCAH9JKmV166aX6xS9+4Wp/8skn9d73vjfjQQEAAORSWjNEH/rQh7Rw4UK1tbXp3e9+tyTpqaee0gMPPKDbbrtNP/nJTwb0BQAAyGcBx3GcUz0pGPQ2sRQIBBSPu9NOforFYqqurtbzb/xBlZVVyfaNrx9w9TWlzGrKSw1tYU/9ysIDr1vEkPcqNaTfTE+9NqXpSJkBAIpVLBZTbV2dOjs7VVVVdfITTlFaM0SJRCLb48i5ynBIVaWh5OuRFe4tR9465t7Oozwc8tQWNgQ7JSlRgqkOyBRImJb/Z3tLDq/BD4EOAKAYnVIN0ebNm/XQQw8NaPv2t7+tcePGacSIEbrxxhvV3d2d1QECAAAMtlMKiD7/+c/rpZdeSr5+4YUXdMMNN6ipqUmLFi3ST3/6U7W0tGR9kAAAAIPplFJm27dv1+233558vXbtWjU0NGjNmjWSpPr6ei1fvlyf+9znsjrIwRAJBRQ5rg7HVPNj4jllZsgtpbaZ+qSm1SRzDZEpjWaKbk39TDJJj2U1HeYUfjo2rwTSWkgKANY5pd+Wb731lmpra5OvN23apMsvvzz5+sILL9SuXbuyNzoAAIAcOKWAqLa2Vm+88YYkqaenR1u3bk0uu5ekQ4cOKRx2r7gCAADIZ6cUEH3wgx/UokWL9Itf/EKLFy9WeXn5gAcxPv/88zrzzDOzPkgAAIDBdEo1RLfffruuuuoqXXLJJaqoqNC3vvUtlZb2197cd999uuyyy7I+yMEQCgyszamOersUpq07TNt0REtOXh+USb2Q6ZlD2V4677k2iLqf/JXte0NNEoAidUoBUU1NjZ544gl1dnaqoqJCodDAYuIHHnhAFRUVWR0gAADAYEvrwYzV1dXG9uHDh2c0GAAAAD8w/w0AAKyX1gxRMXD+//izsrCpJsd9eUzbY5j2Hys11f14eA6RqYbItNWGqV9e1wYVYp0R9TJu2byPXF8AeYTfSAAAwHoERAAAwHrWpsx64o564v1JszLD0vnSoONqC5rSV4aw0rQsPjXNlUkqLOtpr0JMaQ020kODy+v15doByAF+0wAAAOsREAEAAOsREAEAAOtZW0PUm3DUk+ivEYoaindCYW/xYrrbYxi7mOoq3KVMmdW32FQvlC/1J5lc83z5Dn7h2gHIAX5bAAAA6+VFQLRq1SqNHTtW0WhUDQ0N2rJli6fz1q5dq0AgoCuvvHJwBwgAAIqa7wHRunXr1NzcrOXLl2vr1q2aOHGipk+frr17977teb/73e/06U9/Wu9973tzNFIAAFCsAo7jmCpUcqahoUEXXnihVq5cKUlKJBKqr6/X/PnztWjRIuM58Xhcf/VXf6Xrr79ev/jFL3Tw4EE9+OCDnj4vFoupurpar/xutyqrqpLtwyPuip5AzxFvX8JUpxD0PdbMe06+13fk+/hSFdp48wnXDsh7sVhMtXV16uzsVNVxf7+zxdffAj09PWpra1NTU1OyLRgMqqmpSZs3bz7heZ///Oc1YsQI3XDDDSf9jO7ubsVisQEHAADA8XwNiPbv3694PK7a2toB7bW1tWpvbzee8+STT+ob3/iG1qxZ4+kzWlpaVF1dnTzq6+szHjcAACguBbXs/tChQ7ruuuu0Zs0a1dTUeDpn8eLFam5uTr6OxWKqr69X3JGO27lDwWOHXOcGjhlmkwypMCdouIymttRpecM0vWPYzsMokyn+XKQHPKQMAxksp85Juq3QtpYotPHmE64dYD1fA6KamhqFQiF1dHQMaO/o6FBdXZ2r/+uvv67f/e53mjFjRrItkfjTL7KSkhLt2LFDZ5555oBzIpGIIpHIIIweAAAUC1//c6e0tFRTpkxRa2trsi2RSKi1tVWNjY2u/meffbZeeOEFbd++PXl86EMf0vve9z5t376ddBgAAEiL7ymz5uZmzZkzR1OnTtW0adO0YsUKdXV1ae7cuZKk2bNna/To0WppaVE0GtV555034PyhQ4dKkqsdAADAK98DopkzZ2rfvn1atmyZ2tvbNWnSJG3YsCFZaL1z504FB2EJeyjwp+PPAkfd9ULBo2+5TzTUBjmhUndbiSFNFxp4rqn2KGCqUTC2eXtagrEmyVQvke2apESa9UEe77XX+qO8qTXKp9oT6mXSl+2fHQB5w/fnEOXaiZ5DdPrRP7r6+hEQmYMfj20GBVekneXgN2+edZQv4zgVhThmP3CdgJwo6ucQAQAA5AMCIgAAYD3fa4j8EikJKlrSHw+a0mOJAx2utkBp1NUWHOKeukuYag0SKam1YJ+7T8iQkssgtWYsNTI9/0gZ1JVks4bGa+1RodUaFWLtSSGO2Q+p14lrBBQkfnIBAID1CIgAAID1rE2ZlQfjKg/Gk6/j+9yrzBKxN11twXJvle3BIe42J3Vq3bQ6zfBepnVixhRPBqvRAh5j47RTa9le6m1KrWWwQs2UWiONZlCIY841HmsAFCR+IgEAgPUIiAAAgPUIiAAAgPWsrSEKxToUcrqSr3s7drr6OD3HvL1ZSdhTWyD1qdGmZfKGtzevnDcsnTctz88yU62Rp+1BvNZLZFKjQl2RPwpxzPmA6wbkFX76AACA9QiIAACA9QiIAACA9aytIdLeN6Qj5cmXx/5o2O0+7L48pq07TG1OxNCvZOBzhwKBHvd5hqEGjFttGPrJvRVIRnVFHrfp8FRXlEm9BHVFhoHk+X/LsJ1FegrxXgNFgp80AABgPQIiAABgPWtTZj2vPa+e8v601uHd+1x9yk8f5moLV7uX4jt9ve62bkO/1JRZiSElZ0rJJNypMJmWuptm2wOG93N3887rViDOwE/xtDT/FN6fNFqBpVYKbbwArMNvJAAAYD0CIgAAYD0CIgAAYD1ra4j2tb2iY5H+7TW6Dx529SmtHOI+0VAvZGxLxE/e1meqDTLEqKYtPgzv75jCW+MyeUM3Q1tGUsacWlMkDUJdUR5JrSsa9JoiqWiuHVJwX4Gc4KcKAABYj4AIAABYj4AIAABYz9oaor3P/UFdxz0HKDrMvdVGVa9hK4y4oXbH0GZ61o2TUvdjqgMy1vyYnnPjPtPM9MydkA/1LPleQ5PlZxPBA2pj0se1A7KOnyAAAGA9AiIAAGA9a1Nmu1/er/JgKPl67NSRWX3/1PSYJAVS0jLGPt4/wH2uYd29Y1xjn+Xp9jTfz6al+L5s7yEVxbUDgFzgNyMAALAeAREAALAeAREAALCetTVEbxzpUzTQX18x1tAnYFh2HQiFDD2zx1hXY6oD8cp4ruF7ZbK0v1hrUliKn3vUPAHwCb9pAACA9QiIAACA9QiIAACA9aytIeo41qfIcbUJJVH3pQiWGi5P0F1DZKorChj6ZZWp1sJYfzHI48h31KTAFvx/HcgIPy0AAMB6BEQAAMB61qbMuh1HznFL3KPVEVefkmipqy1QEna/manNtDzbw5Jt49YVAABgUDFDBAAArEdABAAArEdABAAArGdtDVF5MKDIcfU6kaEVrj7hIVFXW6DU0GasITr5UnzPS/NNS2e9tgEAgJPiLygAALAeAREAALAeAREAALCetTVEo8pKFD1uW4voadWuPqHycldbIGKqITI8r8iwnYerrshUQ5TtOiDqigAAOCn+WgIAAOsREAEAAOtZmzI7a3iZyoPHp8yqXH2C5ZWuNtOye+PWHR628zCl1RxXi+T4tcS+GNJtxfAdAACDjr8WAADAegREAADAegREAADAetbWEI28oE4V4f6vH6kZ7uoTKBvibjNt3RE2LLs3LcVPaXOO2zqkv1P6Marx/UwMn2GsU/J4btHW6QSL9HsBAFz4jQ8AAKxHQAQAAKxHQAQAAKxnbw1Rw1mqjEaSr0tOG+nqE4iatu4oc7eZ6ooMzyFKfe6QE3RfflOb57odU1se18FkUvNUaDzXaAHp4v9jQEb4CQIAANYjIAIAANazNmVWOWmKqob0p8SClUNdfQKm3egNbcb0mGEpvpN6bgapsKwvsc/2cno/pu9JGbgV2jUptPECKBr89gEAANYjIAIAANbLi4Bo1apVGjt2rKLRqBoaGrRly5YT9l2zZo3e+973atiwYRo2bJiampretj8AAMDJ+F5DtG7dOjU3N2v16tVqaGjQihUrNH36dO3YsUMjRoxw9X/88cc1a9YsXXTRRYpGo7rrrrt02WWX6aWXXtLo0aM9f27wLyYrWFnR//pYp7tTb4+3NysxLJ83LoEf2M/rEntX7dEJ+uVkiX2atUZ5v8Q+jx9PALhQawVkXcBxHMfPATQ0NOjCCy/UypUrJUmJREL19fWaP3++Fi1adNLz4/G4hg0bppUrV2r27Nkn7R+LxVRdXa39r7SpahADInPA4kNAlO2i6nwOiDI5d5ADIt+eQ1RofzgLbbx+4TrBQrFYTLV1ders7FRVVVXW39/Xn6qenh61tbWpqakp2RYMBtXU1KTNmzd7eo8jR46ot7dXw4e7N2eVpO7ubsVisQEHAADA8XwNiPbv3694PK7a2toB7bW1tWpvb/f0HgsXLtSoUaMGBFXHa2lpUXV1dfKor6/PeNwAAKC4+F5DlIk777xTa9eu1eOPP65o1L19hiQtXrxYzc3NydexWEz19fU6WjVa4eOm3MqchOvcYKDL/YaGfp6fHZSaIivE9JhHnlNk6crj9JjkU4qs0NIohTZeAEXN14CopqZGoVBIHR0dA9o7OjpUV1f3tuf+27/9m+688079/Oc/1/nnn3/CfpFIRJFI5IT/DgAA4Ot/opWWlmrKlClqbW1NtiUSCbW2tqqxsfGE59199926/fbbtWHDBk2dOjUXQwUAAEXM95RZc3Oz5syZo6lTp2ratGlasWKFurq6NHfuXEnS7NmzNXr0aLW0tEiS7rrrLi1btkz333+/xo4dm6w1qqioUEVFxQk/J9W+I706FupNvh4TcZ/rJOLuExN93j7AyyozU1op39Nj2Uxz+PGZg4D0mAeFNt58wrUDcsL3gGjmzJnat2+fli1bpvb2dk2aNEkbNmxIFlrv3LlTweMCgq9+9avq6enR3//93w94n+XLl+tzn/tcLocOAACKhO/PIcq1Pz+HaOtru1RZ2V9UPSboXo4fPGp4NhEzRJ66eSqq9muGKMtF1cwQeVBo480nXDtAUpE/hwgAACAf+J4y88tbx/rUG+6f7amvLnP1cXqOuNoCmcyQpJ6bycyPSQ5mg9JeTs9sUPoKcYagEMecD7hugG/46QMAANYjIAIAANYjIAIAANaztoaoqychdfc/Z8gpcW/94ZSUuk/s6/H2AV5qV7K9UszjZ2R1pdipjCVb550I9UIoJNxrIK/wEwkAAKxHQAQAAKxnbcosGAgoGDwuJZS6E/0J2pxsXrF8T4V5/dzBPE9id/p8UohjzgdcNyDv8VMKAACsR0AEAACsR0AEAACsZ20NUVUkpIrISb6+qa4om/zYUDXTsWTz3GKtDTLJl3GcikIccz7gugEFiZ9cAABgPQIiAABgPQIiAABgPWtriKqjIVVGQ/0NTsLVx1STE8hiXZFvNT+5eL9Brg/Km9ogqThqRorhO/iFawcUBX6SAQCA9QiIAACA9axNmQ2LhFR1fMos0evuZNyNfhAHdYLPzIkcLIH3glRYlhXDd/AL1w6wCj/xAADAegREAADAegREAADAetbWEJX1dKqsu3+pvRMKuzuFrL08J5RXNT7p4jsgFdcTsB6/BQAAgPUIiAAAgPUIiAAAgPWsLZIJxdoVcg4nX8erR7v6JEpK3Sfmc62BYfuRgpPP1zcXbP/+2cb1BOARvy0AAID1CIgAAID1rE2ZJd7co0R3eX9D1Uh3J9N0exZ3uzfKJO1lGm8xpNEKEama7OFaAsgBftMAAADrERABAADrERABAADrWVtDFI8dUDx+NPnaGBka6oXijruboUkBD2MImjp5rJfw8v55xa9aJupP8gP3AUCe47cUAACwHgERAACwHgERAACwnrU1ROrrlfqO+/qGGgdTvZCxhshxNwYCJ6/ySbf2SDpB/ZHp/bJdu5FuLRA1JMWB+wigSPHbDQAAWI+ACAAAWM/elFlJ+E/H/3NCYVcXU3qsL2FKdLkFjAmxgQyZNmOqzZQeM43NlEXzmlozMZ6azZQJ24oMPlJcAOAJvy0BAID1CIgAAID1CIgAAID1rK0hCg09XaGK8uTreDjq6tNnKNSJmwp/DAKGbsGU+qCEoUgnaHj/hMfF+KZ6IY8lT8ZzPZ6a/jYi1LcAAPIEf5EAAID1CIgAAID1CIgAAID1rK0hCtScoUBlRfJ1j+OODXsT7ufkxA1FOam1QZJk2rnDSanKMdUZOab3MlTzmCJZU62R1+cQ5aLWKFUGj0gCACCrmCECAADWIyACAADWszZl1jfsDPVVVSVfH+l1p8d6PeaRUlNhkrd0WMgQjiYMy+6DhuRSJkv2vaa9TCmtTFJrXj7ThNQaAGCwMUMEAACsR0AEAACsR0AEAACsZ20NUUdPUEe6++PBuGGJvYmpNibkcdm9q2rGUAhUiHVFJqm1Rl6X/5uku6z/RKhJAgCkYoYIAABYj4AIAABYj4AIAABYz9oaot+8eVRDuvu//qjKqKtPiaHwxVTjYworQ4ZKFSelyVQbFEjtpMzqb0yVUdmuKzJJPTWbzy/KVLZrktJFLRMA5A9miAAAgPUIiAAAgPXyIiBatWqVxo4dq2g0qoaGBm3ZsuVt+z/wwAM6++yzFY1GNWHCBD388MOn/Jnbdsf07B86k8fBY72uozfhuI54QobDcR+O4Ug5z5HhcNxHwng4rsN0rleO47gO0+d6fr+Uwyvzd01/HPnMeP85ODg4OE54DCbfA6J169apublZy5cv19atWzVx4kRNnz5de/fuNfb/1a9+pVmzZumGG27Qtm3bdOWVV+rKK6/Uiy++mOORAwCAYhFwnFOZR8i+hoYGXXjhhVq5cqUkKZFIqL6+XvPnz9eiRYtc/WfOnKmuri499NBDybZ3v/vdmjRpklavXn3Sz4vFYqqurtaSH7cpOqQi2f5X44a7+lZHw662sKHqN2wqqjYVZHvY3NXrQx5DpocwejzXaxQcMJycbtFztguIc1F8DQDIH7FYTHV1ders7FTVcZuzZ4uvq8x6enrU1tamxYsXJ9uCwaCampq0efNm4zmbN29Wc3PzgLbp06frwQcfNPbv7u5Wd3d38nVnZ+ef2o8cHtCv65A7+An1ui+PKSAqSTcgyuCp16ZzTQEMAREAoBgcOnRIkjRY8zi+BkT79+9XPB5XbW3tgPba2lq98sorxnPa29uN/dvb2439W1padNttt7navzTrkjRHDQAA/PLmm2+quro66+9b9M8hWrx48YAZpYMHD+od73iHdu7cOSgXFN7FYjHV19dr165dgzL9iVPD/cgf3Iv8wb3IH52dnRozZoyGD3eXuGSDrwFRTU2NQqGQOjo6BrR3dHSorq7OeE5dXd0p9Y9EIopEIq726upq/s+dJ6qqqrgXeYT7kT+4F/mDe5E/gsHBWQ/m6yqz0tJSTZkyRa2trcm2RCKh1tZWNTY2Gs9pbGwc0F+SNm7ceML+AAAAJ+N7yqy5uVlz5szR1KlTNW3aNK1YsUJdXV2aO3euJGn27NkaPXq0WlpaJEk333yzLrnkEt1zzz264oortHbtWj377LP62te+5ufXAAAABcz3gGjmzJnat2+fli1bpvb2dk2aNEkbNmxIFk7v3LlzwPTYRRddpPvvv19LlizRZz7zGY0fP14PPvigzjvvPE+fF4lEtHz5cmMaDbnFvcgv3I/8wb3IH9yL/DHY98L35xABAAD4zfcnVQMAAPiNgAgAAFiPgAgAAFiPgAgAAFjPuoBo1apVGjt2rKLRqBoaGrRlyxa/h1T0WlpadOGFF6qyslIjRozQlVdeqR07dgzoc+zYMc2bN0+nnXaaKioq9JGPfMT1AE5k35133qlAIKAFCxYk27gXubN792597GMf02mnnaaysjJNmDBBzz77bPLfHcfRsmXLNHLkSJWVlampqUm/+c1vfBxxcYrH41q6dKnGjRunsrIynXnmmbr99tsH7JnFvRg8TzzxhGbMmKFRo0YpEAi49ib1cu0PHDiga6+9VlVVVRo6dKhuuOEGHT48cM/Sk7EqIFq3bp2am5u1fPlybd26VRMnTtT06dO1d+9ev4dW1DZt2qR58+bpqaee0saNG9Xb26vLLrtMXV1dyT633HKLfvrTn+qBBx7Qpk2b9Mc//lFXXXWVj6Mufs8884z+67/+S+eff/6Adu5Fbrz11lu6+OKLFQ6H9cgjj+jXv/617rnnHg0bNizZ5+6779aXv/xlrV69Wk8//bSGDBmi6dOn69ixYz6OvPjcdddd+upXv6qVK1fq5Zdf1l133aW7775bX/nKV5J9uBeDp6urSxMnTtSqVauM/+7l2l977bV66aWXtHHjRj300EN64okndOONN57aQByLTJs2zZk3b17ydTwed0aNGuW0tLT4OCr77N2715HkbNq0yXEcxzl48KATDoedBx54INnn5ZdfdiQ5mzdv9muYRe3QoUPO+PHjnY0bNzqXXHKJc/PNNzuOw73IpYULFzrvec97TvjviUTCqaurc770pS8l2w4ePOhEIhHnu9/9bi6GaI0rrrjCuf766we0XXXVVc61117rOA73IpckOT/60Y+Sr71c+1//+teOJOeZZ55J9nnkkUecQCDg7N692/NnWzND1NPTo7a2NjU1NSXbgsGgmpqatHnzZh9HZp/Ozk5JSm7Q19bWpt7e3gH35uyzz9aYMWO4N4Nk3rx5uuKKKwZcc4l7kUs/+clPNHXqVF199dUaMWKEJk+erDVr1iT//Y033lB7e/uAe1FdXa2GhgbuRZZddNFFam1t1auvvipJeu655/Tkk0/q8ssvl8S98JOXa79582YNHTpUU6dOTfZpampSMBjU008/7fmzfH9Sda7s379f8Xg8+QTsP6utrdUrr7zi06jsk0gktGDBAl188cXJp4u3t7ertLRUQ4cOHdC3trZW7e3tPoyyuK1du1Zbt27VM8884/o37kXu/Pa3v9VXv/pVNTc36zOf+YyeeeYZ/fM//7NKS0s1Z86c5PU2/c7iXmTXokWLFIvFdPbZZysUCikej+uLX/yirr32WkniXvjIy7Vvb2/XiBEjBvx7SUmJhg8ffkr3x5qACPlh3rx5evHFF/Xkk0/6PRQr7dq1SzfffLM2btyoaDTq93CslkgkNHXqVN1xxx2SpMmTJ+vFF1/U6tWrNWfOHJ9HZ5fvfe97+s53vqP7779ff/mXf6nt27drwYIFGjVqFPfCItakzGpqahQKhVyrZTo6OlRXV+fTqOzyqU99Sg899JAee+wxnXHGGcn2uro69fT06ODBgwP6c2+yr62tTXv37tUFF1ygkpISlZSUaNOmTfryl7+skpIS1dbWci9yZOTIkTr33HMHtJ1zzjnauXOnJCWvN7+zBt+//uu/atGiRfroRz+qCRMm6LrrrtMtt9yS3FSce+EfL9e+rq7OtTiqr69PBw4cOKX7Y01AVFpaqilTpqi1tTXZlkgk1NraqsbGRh9HVvwcx9GnPvUp/ehHP9Kjjz6qcePGDfj3KVOmKBwOD7g3O3bs0M6dO7k3WfaBD3xAL7zwgrZv3548pk6dqmuvvTb5v7kXuXHxxRe7Hj/x6quv6h3veIckady4caqrqxtwL2KxmJ5++mnuRZYdOXJkwCbikhQKhZRIJCRxL/zk5do3Njbq4MGDamtrS/Z59NFHlUgk1NDQ4P3DMi4JLyBr1651IpGI881vftP59a9/7dx4443O0KFDnfb2dr+HVtRuuukmp7q62nn88cedPXv2JI8jR44k+3zyk590xowZ4zz66KPOs88+6zQ2NjqNjY0+jtoex68ycxzuRa5s2bLFKSkpcb74xS86v/nNb5zvfOc7Tnl5ufM///M/yT533nmnM3ToUOfHP/6x8/zzzzt/93d/54wbN845evSojyMvPnPmzHFGjx7tPPTQQ84bb7zh/PCHP3RqamqcW2+9NdmHezF4Dh065Gzbts3Ztm2bI8m59957nW3btjm///3vHcfxdu3/5m/+xpk8ebLz9NNPO08++aQzfvx4Z9asWac0DqsCIsdxnK985SvOmDFjnNLSUmfatGnOU0895feQip4k4/Hf//3fyT5Hjx51/umf/skZNmyYU15e7nz4wx929uzZ49+gLZIaEHEvcuenP/2pc9555zmRSMQ5++yzna997WsD/j2RSDhLly51amtrnUgk4nzgAx9wduzY4dNoi1csFnNuvvlmZ8yYMU40GnX+4i/+wvnsZz/rdHd3J/twLwbPY489ZvwbMWfOHMdxvF37N99805k1a5ZTUVHhVFVVOXPnznUOHTp0SuMIOM5xj+IEAACwkDU1RAAAACdCQAQAAKxHQAQAAKxHQAQAAKxHQAQAAKxHQAQAAKxHQAQAAKxHQASgoPzjP/6jrrzySr+HAaDIsNs9gLwRCATe9t+XL1+u//iP/xDPkwWQbQREAPLGnj17kv973bp1WrZs2YANUCsqKlRRUeHH0AAUOVJmAPJGXV1d8qiurlYgEBjQVlFR4UqZXXrppZo/f74WLFigYcOGqba2VmvWrFFXV5fmzp2ryspKvfOd79Qjjzwy4LNefPFFXX755aqoqFBtba2uu+467d+/P8ffGEC+ICACUPC+9a1vqaamRlu2bNH8+fN100036eqrr9ZFF12krVu36rLLLtN1112nI0eOSJIOHjyo97///Zo8ebKeffZZbdiwQR0dHfqHf/gHn78JAL8QEAEoeBMnTtSSJUs0fvx4LV68WNFoVDU1NfrEJz6h8ePHa9myZXrzzTf1/PPPS5JWrlypyZMn64477tDZZ5+tyZMn67777tNjjz2mV1991edvA8AP1BABKHjnn39+8n+HQiGddtppmjBhQrKttrZWkrR3715J0nPPPafHHnvMWI/0+uuv66yzzhrkEQPINwREAApeOBwe8DoQCAxo+/PqtUQiIUk6fPiwZsyYobvuusv1XiNHjhzEkQLIVwREAKxzwQUX6Ac/+IHGjh2rkhJ+DQKghgiAhebNm6cDBw5o1qxZeuaZZ/T666/rZz/7mebOnat4PO738AD4gIAIgHVGjRqlX/7yl4rH47rssss0YcIELViwQEOHDlUwyK9FwEYBh0e+AgAAy/GfQgAAwHoERAAAwHoERAAAwHoERAAAwHoERAAAwHoERAAAwHoERAAAwHoERAAAwHoERAAAwHoERAAAwHoERAAAwHoERAAAwHr/B+3+cz6hQ2bxAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "ex.viz.plot_spatio_temporal(trj)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The trajectory is identical to the one obtained by simulation in physical space" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Array(True, dtype=bool)" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "jnp.allclose(\n", + " trj,\n", + " ex.rollout(burgers_stepper, 100, include_init=True)(u_0),\n", + " atol=1e-5,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Ensemble simulation\n", + "\n", + "One particular feature of `Exponax` that is highly relevant for the integration\n", + "with deep learning is the batched execution.\n", + "\n", + "Rather straightforward, we can `jax.vmap` a timestepper to operate in muliple\n", + "states at once." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "NUM_SPATIAL_DIMS = 1\n", + "DOMAIN_EXTENT = 3.0\n", + "NUM_POINTS = 100\n", + "DT = 0.1\n", + "\n", + "burgers_stepper = ex.stepper.Burgers(NUM_SPATIAL_DIMS, DOMAIN_EXTENT, NUM_POINTS, DT)\n", + "\n", + "ic_gen = ex.ic.RandomTruncatedFourierSeries(\n", + " NUM_SPATIAL_DIMS,\n", + " cutoff=5,\n", + " max_one=True,\n", + ")\n", + "\n", + "one_u_0 = ic_gen(NUM_POINTS, key=jax.random.PRNGKey(0))\n", + "\n", + "multiple_u_0 = ex.build_ic_set(\n", + " ic_gen, num_points=NUM_POINTS, num_samples=10, key=jax.random.PRNGKey(0)\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "((1, 100), (10, 1, 100))" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "one_u_0.shape, multiple_u_0.shape" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [], + "source": [ + "one_u_1 = burgers_stepper(one_u_0)\n", + "# burgers_stepper(mutliple_u_0) # This will fail because the vanilla timestepper is single-batch only\n", + "multiple_u_1 = jax.vmap(burgers_stepper)(multiple_u_0)" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "((1, 100), (10, 1, 100))" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "one_u_1.shape, multiple_u_1.shape" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "When using `jax.vmap`, we essentially share the same dynamics across all initial\n", + "states. What if we wanted to simulate the same state but with three different\n", + "dynamics? We could create a list of timesteppers and then loop over time (for\n", + "example, with a list comprehension). However, sequential looping is slow. There\n", + "is an easy way to also use JAX' automatic vectorization for that. For this we\n", + "create an ensemble of three different Burgers steppers (this will only work if\n", + "the parameter we vmap over does not change the shape the timesteppers attribute\n", + "arrays)." + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [], + "source": [ + "DIFFUSIVITIES = jnp.array([0.1, 0.3, 0.7])\n", + "\n", + "burgers_stepper_ensemble = eqx.filter_vmap(\n", + " lambda nu: ex.stepper.Burgers(\n", + " NUM_SPATIAL_DIMS, DOMAIN_EXTENT, NUM_POINTS, DT, diffusivity=nu\n", + " )\n", + ")(DIFFUSIVITIES)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If we inspect the single timestepper PyTree structure next to the ensemble\n", + "timestepper PyTree structure, we see an additional batch axis in the internal\n", + "arrays." + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(Burgers(\n", + " num_spatial_dims=1,\n", + " domain_extent=3.0,\n", + " num_points=100,\n", + " num_channels=1,\n", + " dt=0.1,\n", + " dx=0.03,\n", + " _integrator=ETDRK2(\n", + " dt=0.1,\n", + " _exp_term=c64[1,51],\n", + " _nonlinear_fun=ConvectionNonlinearFun(\n", + " num_spatial_dims=1,\n", + " num_points=100,\n", + " dealiasing_mask=bool[1,51],\n", + " derivative_operator=c64[1,51],\n", + " scale=1.0,\n", + " single_channel=False\n", + " ),\n", + " _coef_1=f32[1,51],\n", + " _coef_2=f32[1,51]\n", + " ),\n", + " diffusivity=0.1,\n", + " convection_scale=1.0,\n", + " dealiasing_fraction=0.6666666666666666,\n", + " single_channel=False\n", + " ),\n", + " Burgers(\n", + " num_spatial_dims=1,\n", + " domain_extent=3.0,\n", + " num_points=100,\n", + " num_channels=1,\n", + " dt=0.1,\n", + " dx=0.03,\n", + " _integrator=ETDRK2(\n", + " dt=0.1,\n", + " _exp_term=c64[3,1,51],\n", + " _nonlinear_fun=ConvectionNonlinearFun(\n", + " num_spatial_dims=1,\n", + " num_points=100,\n", + " dealiasing_mask=bool[3,1,51],\n", + " derivative_operator=c64[3,1,51],\n", + " scale=1.0,\n", + " single_channel=False\n", + " ),\n", + " _coef_1=f32[3,1,51],\n", + " _coef_2=f32[3,1,51]\n", + " ),\n", + " diffusivity=f32[3],\n", + " convection_scale=1.0,\n", + " dealiasing_fraction=0.6666666666666666,\n", + " single_channel=False\n", + " ))" + ] + }, + "execution_count": 20, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "burgers_stepper, burgers_stepper_ensemble" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First task is to make three different predictions from the single initial\n", + "condition." + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [], + "source": [ + "ensembled_u_1 = eqx.filter_vmap(lambda stepper: stepper(one_u_0))(\n", + " burgers_stepper_ensemble\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This adds a three-dimensional batch axis to the state" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(3, 1, 100)" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ensembled_u_1.shape" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can also use both vmapping over the ensemble and the multiple initial states" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [], + "source": [ + "ensembled_multiple_u_1 = eqx.filter_vmap(\n", + " lambda stepper: jax.vmap(stepper)(multiple_u_0)\n", + ")(burgers_stepper_ensemble)" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(3, 10, 1, 100)" + ] + }, + "execution_count": 24, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ensembled_multiple_u_1.shape" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Be mindful about the order of nested vmappings as they affect the order of axes\n", + "in the returned arrays." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "jax_fresh", + "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.10.13" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}