Skip to content

Instantly share code, notes, and snippets.

@jessegrabowski

jessegrabowski/sparse_jac.ipynb

Last active June 26, 2025 10:53
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save jessegrabowski/6b3dc83898ed355b850a7f698476ef71 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save jessegrabowski/6b3dc83898ed355b850a7f698476ef71 to your computer and use it in GitHub Desktop.
Download ZIP
Sparse Jacobians in Pytensor
Display the source blob
Display the rendered blob
Raw
sparse_jac.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "54289e6e",
"metadata": {},
"outputs": [],
"source": [
"import pytensor.tensor as pt\n",
"import pytensor\n",
"from pytensor.graph.basic import explicit_graph_inputs\n",
"\n",
"import numpy as np\n",
"from scipy import sparse\n",
"import networkx as nx"
]
},
{
"cell_type": "markdown",
"id": "02bed835",
"metadata": {},
"source": [
"Test system:\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
"f_1(x) &= x_1^2 + \\sin x_2, \\\\\n",
"f_2(x) &= x_2^3, \\\\\n",
"f_3(x) &= e^{x_3} + x_1, \\\\\n",
"f_4(x) &= e^{x_4} + 3x_1, \\\\\n",
"f_5(x) &= x_5^3 + x_2 x_5,\n",
"\\end{aligned}\n",
"\\qquad\n",
"x = \\begin{bmatrix} x_1 \\\\ x_2 \\\\ x_3 \\\\ x_4 \\\\ x_5 \\end{bmatrix}.\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "39f02063",
"metadata": {},
"outputs": [],
"source": [
"def get_jacobian_connectivity(variables, equations):\n",
" connectivity_pattern = []\n",
"\n",
" n_eqs = len(equations)\n",
" n_vars = len(variables)\n",
" \n",
" for i, eq in enumerate(equations):\n",
" for j, var in enumerate(variables):\n",
" if var in explicit_graph_inputs(eq):\n",
" connectivity_pattern.append((i, j))\n",
" \n",
" \n",
" return list(zip(*connectivity_pattern))"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "42dbdf04",
"metadata": {},
"outputs": [],
"source": [
"x0, x1, x2, x3, x4 = variables = pt.dscalars('x0 x1 x2 x3 x4'.split())\n",
"\n",
"equations = [\n",
" x0 ** 2 + pt.sin(x1),\n",
" x1 ** 3,\n",
" pt.exp(x2) + x0,\n",
" pt.exp(x3) + 3 * x0,\n",
" x4 ** 3 + x1 * x4\n",
"]\n",
"\n",
"n_vars = len(variables)\n",
"n_eqs = len(equations)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "362aa6f2",
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"[(0, 0, 1, 2, 2, 3, 3, 4, 4), (0, 1, 1, 0, 2, 0, 3, 1, 4)]"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"test_values = dict.fromkeys(variables, 1.0)\n",
"get_jacobian_connectivity(variables, equations)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "6613e45f",
"metadata": {},
"outputs": [],
"source": [
"# TODO: Make an Op?\n",
"def _greedy_color(\n",
" connectivity: sparse.spmatrix,\n",
" strategy: str = 'largest_first',\n",
") -> tuple[np.ndarray, int]:\n",
" assert connectivity.ndim == 2\n",
" assert connectivity.shape[0] == connectivity.shape[1]\n",
" G = nx.convert_matrix.from_scipy_sparse_array(connectivity)\n",
" coloring_dict = nx.algorithms.coloring.greedy_color(G, strategy)\n",
" \n",
" indices, colors = list(zip(*coloring_dict.items()))\n",
" coloring = np.asarray(colors)[np.argsort(indices)]\n",
" \n",
" n_colors = np.unique(coloring).size\n",
" \n",
" return G, coloring, n_colors"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "913fd5a3",
"metadata": {},
"outputs": [],
"source": [
"coords = get_jacobian_connectivity(variables, equations)\n",
"data = np.ones(len(coords[0]), dtype='bool')\n",
"\n",
"sparsity = sparse.coo_array((data, coords), (5, 5))\n",
"output_connectivity = sparsity @ sparsity.T\n",
"input_connectivity = sparsity.T @ sparsity"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "0d7c816c",
"metadata": {},
"outputs": [],
"source": [
"G, output_coloring, n_colors = _greedy_color(output_connectivity)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "f3225e84",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0, 1, 1, 2, 2])"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"output_coloring"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "5f90d65c",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAGFCAYAAABg2vAPAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAYolJREFUeJzt3Xd8TffjBvDn3kxkSYhaEUlQtRJEFiF2CYkgCULsUHsUrQ6qrbbUKKWUGLFrJKhNRNZNhCSU1mgWNUKQve/5/dEvP4rIPnc879frvr6vb+49n/skjZznns/5nCMRBEEAERERqS2p2AGIiIhIXCwDREREao5lgIiISM2xDBAREak5lgEiIiI1xzJARESk5lgGiIiI1JxmaV4kl8tx79496OvrQyKRVHUmIiIiqgSCICAzMxMNGjSAVPr2z/+lKgP37t1D48aNKy0cERERVZ87d+6gUaNGb32+VGVAX1//xWAGBgaVk4yIiIiqVEZGBho3bvxiP/42pSoDz6cGDAwMWAaIiIiUzLum+HkCIRERkZpjGSAiIlJzLANERERqjmWAiIhIzbEMEBERqTmWASIiIjXHMkBERKTmWAaIiIjUHMsAERGRmmMZICIiUnMsA0RERGqOZYCIiEjNsQwQERGpuVLdtZDeThAEZGZmIjU1FY8ePUJqaipSU1ORmZmJoqIiFBUV4e7du7h06RIGDhwIDQ0NaGpqQltbG3Xq1IGpqemLR506daCpyf8kRERUvbjnKaX09HRcuXIFcXFxiIuLwx9//IH79+8jNTUV+fn5r7xWKpWiVq1a0NTUhIaGBh4/fgwASEpKQnFxMYqKipCfn4+8vLxXtpNIJDA2NoapqSksLS1hbW2Ndu3awdraGhYWFpBKeSCHiIgqH8vAWyQlJSEoKAghISGIi4tDYmIiAEBbWxutW7dGmzZt0KtXr1c+2derVw+mpqYwNjaGhobGO98jOzv7xZGE54+HDx/i4cOHuHnzJjZt2oQHDx4AAPT09NCuXTt07NgRrq6u6Nq1K7S0tKr0Z0BEROpBIgiC8K4XZWRkwNDQEOnp6TAwMKiOXNVOEARcuXIFgYGBCAwMRFxcHLS1teHk5AQbGxtYW1vD2toa77//frXuhB88eID4+HjExcUhPj4eYWFhuHPnDoyMjNC/f3+4u7ujb9++0NPTq7ZMRESKTBAEPHv27LXp25encIuLi3Hr1i1cuHABY8aMeXEk901TuKampqhbty60tbXF/tbKrLT7b7UvAwUFBQgICMAPP/yAmzdvwtDQEP3794ebmxv69u2rcN+vIAiIi4t7UVquXLkCHR0d+Pj4YMGCBbCyshI7IhFRtSguLsbt27dffGB6PoX74MEDFBYWvvJaDQ0N1K1bF/r6+tDU1ISmpiauXr0KAGjXrt2Lc7zy8/Px+PFjZGVlvfZ+RkZGsLKyejF9a21tjbZt2yrcfuJlLAPvkJOTg02bNmHZsmW4e/cuPDw84Ofnh27duilV+0tMTMTevXuxevVqpKamwsvLC59++ilat24tdjQiokoll8sRHR2NwMBAhISE4MqVK8jJyQEANGzYENbW1mjTpg0aNWr02if72rVrl+m8q5ycnNeOKjyfwn1eOp4XDgsLixdTuP3794exsXGVfP/lwTJQgsOHD2PixIl4/Pgxhg8fjgULFuCDDz4QO1aF5OXlwd/fH99//z1SUlLg6+uLNWvWQF9fX+xoRETllp+fj+DgYAQGBiIoKAgPHjxAnTp10Lt37xdTuO3atUPdunWrNVdBQQH++uuvV6Zwo6OjoaGhga5du8Ld3R1ubm4wMzOr1lz/Ver9t1AK6enpAgAhPT29NC9XWHl5ecL06dMFAMLAgQOFv//+W+xIla6goEDYsGGDoKenJ1hZWQkxMTFiRyIiKrOnT58KX331lWBiYiIAECwsLIQ5c+YIoaGhQlFRkdjx3uiff/4R1q9fL/Tp00fQ0tISAAj9+vUTwsPDRctU2v232pSB5ORkwcbGRtDW1hbWrl0ryOVysSNVqVu3bgkdOnQQtLS0hHXr1okdh4ioVB4+fCgsWLBA0NfXF3R0dIQpU6YIV65cUbq/2c+ePRO2bNkifPDBBwIAoVu3bsLp06er/ftgGXhJRkaG0Lp1a8Hc3FyIjY0VO061yc/PF6ZOnSoAEPbs2SN2HCKit5LL5cKGDRuEmjVrCnp6esK8efOE+/fvix2rwoqLi4WDBw8KHTp0EAAIvXv3rtbvi2Xgf4qKigRXV1fBwMBAuH79uthxqp1cLhdGjBgh6OrqCtHR0WLHISJ6zdOnT4WhQ4cKAAQ/Pz8hLS1N7EiVTi6XC4cPHxbee+89wdTUVDh58mS1vG9p998qf0m7lStX4tixY9i7dy9atmwpdpxqJ5FIsGnTJrRr1w6DBg167WqJRERiunbtGmxsbHDq1Cns378fv/zyi0KdjV9ZJBIJBgwYgLi4ONjY2KBPnz5YtGiR2LFeUOkyUFxcjDVr1mD06NHo27ev2HFEo6uri61bt+Kff/7BoUOHxI5DRAQASE1NRf/+/aGvr4+4uDgMHjxY7EhVrl69ejh27BiWLFmCxYsXY926dWJHAqDiZeD06dNISUmBn5+f2FFE9/7778PZ2Rm//vqr2FGIiJCfnw8PDw/k5eXh999/h7m5udiRqo1UKsVnn32G6dOnY/r06Th79qzYkVS7DBw9ehRWVlawtbUVO4pC8PHxwblz5954ZS0iour0zTffICYmBoGBgWjcuLHYcUTx448/omfPnvDy8kJ2draoWVS6DBQWFqJ27dqQSCRiR1EItWvXBgAUFRWJnISI1FlBQQE2bNiAiRMnwt7eXuw4otHU1MQvv/yCJ0+eYN++faJmUekyQEREiufIkSNITU3FhAkTxI4iOnNzc/Tu3RsbN24UNYdKl4H33nsPN2/eRG5urmgZsvOLcO1eOmJTnuLavXRk54v3qTw+Ph76+vqoWbOmaBmIiM6cOYOWLVuiTZs2YkdRCN7e3pDJZMjMzBQtg6Zo71wNRo0aha+++gr79+/HyJEjq+19bz3MxM6oFATfSEXKkxy8fPMHCQAz45pwaWGKEXZmaFaveu4dUFRUBH9/f/j4+CjVjZiISPUUFxfztusvef6zkMvlomVQ6TJgaWmJHj16YPXq1fD29oaWllaVvt+dJzn49NBVhN5+DA2pBMXy1+8BJQBIfpKDgKhkbI1MQherOvh2UBs0Nq7aT+t79+7FvXv3eFiOiBSC8O575KkNRfhZqPQ0AQAsXrwY8fHxmD59epX+wPdcTEHPlSGISEgDgDcWgZc9fz4iIQ09V4Zgz8WUKssWHx8PPz8/DBkyBDY2NlX2PkREpWFubo7r168jIyNDtAyKNIUrk8lgYmKCWrVqiZZBLW5hvHnzZowfPx5r1qzB1KlTK338tcG3sPzUzQqPM7d3c0x1aVYJif7fw4cP0alTJ5iYmCA0NFTUXzYiIgC4e/cumjRpgnXr1lXrdWAUcQo3Pz8fDRs2hK+vL3788cdKH7+0+2+1KAMAMHv2bKxcuRJz587FN998U2nz5nsupmDBwauVMhYAfO/RBl62lXP/68jISAwbNgz5+fmIjo5W27W8RKR4BgwYgL///hsxMTFVflJzaaZwn3v+fHVN4f7888+YOnUq/vzzT7z//vuVPn5p998qP03w3PLly/Hjjz9i9erV6Ny5MxISEio85p0nOfjy8LVKSPf/vjh8DXee5FRoDLlcjqVLl6JLly5o2LAhZDIZiwARKZSvv/4aycnJGD16dJWeOKfIU7gXLlzArFmz4OfnVyVFoCzUpgxIpVLMnj0b4eHhSEtLg7W1Nb788kukpaWVe8xPD11FUQm/WPL8HDwN9sfDPZ/jzurhSP7OFc9Cd5Y4ZpFcwKeHynekQRAEnDhxAp07d8bChQuxYMEChISEoEmTJuUaj4ioqrRr1w47duzAb7/9hkWLFlXJOV1rg29hwcGryC+Sv7ME/FexXEB+kRwLDl7F2uBblZ7t77//xuDBg+Hk5IQ1a9ZU+vhlpTZl4DlbW1vExsZi3LhxWLZsGZo0aYKPP/4YDx48KNM4tx5mIvT24xJ/weS5mciMOwmhuBA1m5fuKlvFcgGhtx/jdmrp15vK5XIcPHgQtra2+PDDDyGXy3Hu3Dl8/fXX0NRU6QUjRKTEBg0ahKVLl2LJkiUYNWpUpa6z33MxpVLO5QKA5aduYm8lHiEICgqCra0tjIyMsH///ipf6VYaalcGAMDAwAArV65EUlISpk+fjo0bN8Lc3Bze3t7Ys2cP0tPT3znGzqgUaEhLvsyxhqEpGs/cg/dGfAejrr6lzqchlWCHrORfPEEQEBMTg88//xwtW7bE4MGDYWBggDNnziAyMhLdunUr9fsREYllwYIF2LFjBwIDA9G+fXtcunSpwmOWdQo3M/4kkr9zRcqPQ976msqYws3Ly8P06dPh7u4OZ2dnREVFwcTEpEJjVha1LAPPmZqa4ttvv0VycjIWL16MW7duYdiwYahbty769OmD9evXIzEx8Y2Hr4JvpL7zsJNEIinXfRGK5QKCb6a+9vWsrCycPn0aU6dOhZmZGWxtbfHzzz/Dzs4OEREROHfuHHr06MF7MRCRUhkxYgRiY2NhaGgIBwcH+Pn5Vei8rndN4b6sKPMxnp7zh4aeccmvq8AUbmFhIbZs2YK2bdtiw4YNWLt2LQ4dOgRj45LfszrxGDIAIyMjzJ8/H/Pnz0dKSgqCgoIQGBiIadOmobi4GEZGRrC2tn7xaP5BG6RUsCG+S0paDg4e+R03/riCuLg4xMXF4datWxAEAWZmZvDw8IC7uzs6d+6sEIeYiIgqwsrKChEREVi1ahV+/PFHbN68GcOGDcMnn3yCDz74oNTjPJ/CLa0nJ36GbuNWkOrqI+dG+Ftf9/IUrpVp6ZYd5ubmYvPmzVi2bBlSUlLg7u6OAwcOKORlmNVmaWF5PHnyBBEREYiPj3+xQ759+za0TJuiwdiynfBRnJOOuz+NgKHTMBh1GVGqbe75T4NOziO8//77aNmyJT744AO0bdsWLVu2hIaGBqRS6Yv/LcuDRw6ISJHl5uZi06ZNWLZsGe7cuQNnZ2e4u7vDzc0NFhYWJW676PA1BEQll+qEwaw/gvHk1Do0GL8ezy4EIOdGOMzm7H/r6zWkEoy0a4JFA1uVmP3MmTMIDAxEUFAQnj59Cm9vb3zyySdo3br1OzNVttLuv3lkoATGxsZwdXWFq6vri69lZmbi0IVYfBFW9TeUkGhoISsrCzExMYiJianUsctTIp4/yrttdW/HrNJyT1URialGjRqYNm0a/Pz8sGfPHuzbtw+ffPIJZs+ejTZt2sDd3R19+/ZFu3btXruQWmmmcAGgOPsZnp79FbW7jYamQZ1S5Xo+hbsI/18G5HI5kpKSEBYWhqCgIJw4cQI5OTlo3rw5xo0bhwkTJsDKyqpsPwARsAyUkb6+PjrYtAPCwqr8vbZt2QwzfSnkcvkbH8XFxW99rqRHdW9X2m2LiooULqsqkEgkSlFcFO09lTGrqhU/bW1tjBo16sVKg5MnTyIwMBBr1qzBkiVLIJFI0Lx58xdTuC1atUXKk9KdK/Dk1DpoGTeEnk2/MmVKTsvBL5u24PqVWMTFxSE+Pv7FZZXt7Ozw+eefw93dXfTrBpQVy0A5mJvUggRAVd5aQgJgQDd71NLhfyIxCYKg9CVLEbZ707ZvKn9i/3xUgbIUl/Ju27BhQ0yYMAGPHz/Go0ePkJqaisjISAQGBkJu2KBUU7jZf4Uj53Y06o/5qVwFasZnX8OitjbatWuHfv36wdraGu3bt4epqWl5/pMpBO5pyqGWjibMjGsiuQpPIjQzqckioAAkEgk0NDSgoaEhdhSqBs/LnyKUJUXdrrzbFhUVVWlWIyMjyGu/+3C/vCAXT06vh0GHAdDUM4Y8L+vf//byf29UJM/LAqSakGrrvnWM8xfC4NC8fuX80ikI7m3KyaWFaalOUsn9OwbywjwIBbkAgMK0O8j+698phhqWHSHVev0XTkMqgUtz5W2YRMqK5U+5XbuXjv5rSp7CledkQJ79DBnRh5ARfei15++s8kaNZvYwHfzZW8cw0Kva+xWIgWWgnEbYmWFrZNI7X5d2ch2KM/7/mgE5f4Uh539loOGkzZAavV4GiuUCfOwr52ZFRETqojRTuBp6tVFv2LevfT1dth/5d/6A6dBFkNYs4e5+/3sfVcMyUE7N6umji1UdRCSklXh0oNFH/mUaV0MqgaOFSanXsRIR0b9KM4Ur0dSGbpO2r3096+pZQCJ943MvU9UpXLW+AmFFfTuoDTTfcUnishGgKZXg20GKd0EKIiJl4NLC9J2Xii8vVZ7CZRmogMbGNbG4hItPlJ0EzdIvo76BdiWOSUSkPkbYmZX5DoUAUMd1VokXHAJUewqXZaCCvG3NMLd380oZq7tJFk6uX4QBAwZU6t27iIjUxfMp3MreuWlIJehiVUdlp3BZBirBVJdm+M6jDXQ0pWU+PKUhlUBHU4rvPdrAf64Xjh8/joiICHTp0gX//PNPFSUmIlJNRUVFMEk4haLCfODdV9svNVWfwmUZqCTetmY4M6srHC3+vR3lO29v/L/nHS1McGZWV3jZ/nvoqWfPnggPD8eTJ09gZ2eH+Pj4qg1ORKQi/vnnH3Tv3h1rv18MF4PHQCVekfGrga3Q2Fj1lhQ+xxsVVYFbDzOxMyoFwTdTkZKW88oyFwn+PRvVpbkpfOzN3nrI6f79+3B1dcXNmzexf/9+9OnTp1qyExEpo5MnT8LHxwc6OjrYs2cPOnfujLXBt7D81M0Kj/1x7xaY4qL49xd4k9Luv1kGqlh2fhGS0rJRUCSHtqYU5ia1Sr0sJSsrC97e3jhx4gTWr1+PCRMmVHFaIiLlUlRUhC+++AJLly5Fv379sG3bNtSp8/9XItxzMQVfHr6GIrlQphMLNaQSaEol+GpgqxdHbpUR71qoIGrpaKJVA8Nybaunp4fAwEDMmDEDEydORGJiIr7++mtIpZzdISK6e/cuhg0bhsjISHz//feYO3fua38fvW3N4GRZB58euorQ24+hIZWUWAqeP+9oYYJvB7VR6amBl7EMKDhNTU2sXbsWlpaWmDt3LhITE7Flyxbo6r79utlERKru+PHjGDlyJGrUqIGQkBA4OTm99bWNjWsiYJxdpUzhqipOEyiRAwcOwMfHBx07dkRgYCBMTEzEjkREVK2Kiorw+eef47vvvkP//v2xbdu2cv0trMgUrjLhOQMqSiaTYeDAgTAyMsKxY8dgZaWcJ7UQEZXVnTt3MGzYMMhkMnz33XeYPXs2p03fobT7b/4UlYy9vT0iIyMhkUjg4OCAyMhIsSMREVW5Y8eOwcbGBikpKbhw4cIbzw+g8uNPUglZWloiMjISLVu2hIuLC/bvL/kSmkREyqqwsBDz589H//794eDggNjYWDg6OoodS+WwDCgpY2NjnD59Gh4eHhg6dCiWL1+OUsz4EBEpjTt37qBbt25YsWIFli9fjsOHD/NcqSqiemdLqBEdHR3s2LEDFhYW+Pjjj5GQkICffvoJmpr8z0pEyu3o0aPw9fWFnp4eQkNDYW9vL3YklcYjA0pOKpXi66+/xq+//oqNGzfCzc0NWVlZYsciIiqXwsJCzJs3DwMGDICTkxNiY2NZBKoBy4CKGD9+PI4dO4bQ0FA4Ozvj3r17YkciIiqTlJQUdO3aFStXrsSPP/6IoKAgGBsbix1LLbAMqJDevXsjPDwcjx49gp2dHa5evSp2JCKiUjly5Aisra3xzz//IDQ0FLNnz4akEm80RCVjGVAxbdq0QVRUFOrUqQMnJyecPn1a7EhERG9VWFiIuXPnYuDAgejSpQunBUTCMqCCGjRogAsXLqBLly7o168f/P39xY5ERPSa5ORkODs7Y/Xq1VixYgUCAwM5LSASlgEVpa+vj6CgIIwfPx7jxo3DZ599xqWHRKQwDh8+DBsbG9y/fx9hYWGYNWsWpwVExDKgwjQ1NbFu3Tr88MMP+Oabb+Dj44P8/HyxYxGRGisoKMCcOXPg5uYGZ2dnxMbGws7OTuxYao8L0lWcRCLBxx9/DHNzc4wcORJ3797FoUOHeCiOiKpdUlISvLy8EBsbi1WrVmH69Ok8GqAgeGRATQwdOhTnzp3DtWvX4OjoiISEBLEjEZEaCQoKgo2NDVJTUxEeHo4ZM2awCCgQlgE14ujoCJlMBrlcDnt7e8hkMrEjEZGKKygowKxZs+Du7g4XFxfExsbC1tZW7Fj0HywDasbKygqRkZFo3rw5XFxccPDgQbEjEZGKSkpKQufOnfHzzz9j9erVOHDgAIyMjMSORW/AMqCGTExMcObMGbi5uWHIkCFYuXIlVxoQUaUKDAyEjY0NHj9+jPDwcJ4foOBYBtSUrq4udu3ahfnz52P27NmYPn06iouLxY5FREquoKAAM2fOxKBBg9C9e3dcvnyZ0wJKgKsJ1JhUKsXSpUvRtGlTfPTRR0hOTsbu3btRq1YtsaMRkRJKTEyEl5cX4uPjsWbNGkyZMoVHA5QEjwwQJk6ciKNHjyI4OBhdu3bF/fv3xY5ERErm0KFDsLGxQVpaGiIiIjB16lQWASXCMkAAgL59+yIsLAwPHjyAvb09rl27JnYkIlIC+fn5mDFjBjw8PNCrVy9cvnwZHTp0EDsWlRHLAL3Qrl07yGQyGBkZwdHREWfPnhU7EhEpsISEBDg5OeGXX37B2rVrsW/fPhgaGoodi8qBZYBe0ahRI4SGhsLBwQF9+/bF1q1bxY5ERArowIEDsLGxwbNnzxAZGcnzA5QcywC9xsDAAEeOHMGYMWMwZswYfPnll1x6SEQA/p0WmDZtGoYMGYI+ffrg0qVLaN++vdixqIK4moDeSEtLCxs2bIClpSUWLFiAxMREbNq0Cdra2mJHIyKR/P333/Dy8sLVq1exbt06TJo0iUcDVATLAL2VRCLB/PnzYW5uDl9fX9y5cwcHDx5E7dq1xY5GRNVs//79GDduHExNTSGTyWBjYyN2JKpEnCagd/Ly8sKZM2dw5coVODo6IjExUexIRFRN8vLyMHXqVAwdOhQffvghLl26xCKgglgGqFQ6d+4MmUyGwsJC2NvbIzo6WuxIRFTFbt++DUdHR2zatAnr16/H7t27YWBgIHYsqgIsA1RqzZo1Q2RkJCwtLdGtWzcEBgaKHYmIqsi+ffvQvn17ZGVlQSaT8fwAFccyQGVSt25dnD17Fq6urvDw8MDq1avFjkRElSgvLw8fffQRvLy80L9/f8TExMDa2lrsWFTFWAaozGrUqIE9e/Zg7ty5mDlzJmbMmMGbHBGpgFu3bsHBwQH+/v7YsGEDdu3axWkBNcHVBFQuUqkUP/zwAywsLDBlyhQkJydj586dvMkRkZLau3cvJkyYgPfeew9RUVFo166d2JGoGvHIAFXIpEmTcOTIEZw5cwbdunXDgwcPxI5ERGWQl5eHyZMnw9vbG66urrh06RKLgBpiGaAK69evH0JDQ3Hv3j3Y29vj+vXrYkciolK4efMm7O3tsWXLFmzcuBE7d+6Evr6+2LFIBCwDVClsbGwgk8lgYGAAR0dHBAcHix2JiEqwe/dudOjQAbm5uYiKisKECRO4WkCNsQxQpWncuDHCwsLQqVMn9OnTB9u3bxc7EhH9R25uLvz8/DB8+HAMHDgQMTExnBYglgGqXAYGBvj9998xatQo+Pr64quvvuJNjogUxI0bN2Bvb4/t27fj119/xY4dOzgtQAC4moCqgJaWFn799VdYWFhg4cKFSEhIwMaNG3mTIyIR7dq1C35+fmjYsCGioqLQtm1bsSORAuGRAaoSEokEn376KXbt2oXdu3ejb9++ePbsmdixiNRObm4uJk6ciBEjRsDd3R0xMTEsAvQalgGqUsOGDcOZM2cQFxcHJycnJCUliR2JSG389ddfsLOzQ0BAADZv3ozt27dDT09P7FikgFgGqMp16dIFkZGRyMvLg729PWJiYsSORKTydu7ciY4dO6KwsBDR0dEYO3YsVwvQW7EMULVo0aIFIiMjYW5ujq5du+LIkSNiRyJSSbm5uZgwYQJ8fHzg4eGBixcvok2bNmLHIgXHMkDVxtTUFMHBwejbty/c3d2xdu1asSMRqZTn0wI7d+6Ev78/tm3bxmkBKhWWAapWNWrUwG+//YZZs2Zh2rRpmDVrFm9yRFQJduzYgY4dO6KoqAgXL17EmDFjOC1ApcYyQNVOKpVi+fLl+Pnnn/HTTz9hyJAhyMnJETsWkVLKycnB+PHjMXLkSAwZMgQXL15Eq1atxI5FSoZlgETz0UcfISgoCKdOnYKLiwsePnwodiQipfLnn3/Czs4Ou3btwpYtW7B161beOZTKhWWAROXq6ooLFy4gJSUFDg4O+Ouvv8SORKQUtm/fjo4dO0Iul+PixYsYPXq02JFIibEMkOg6dOiAqKgo1KxZEw4ODggJCRE7EpHCysnJwdixY+Hr6wsvLy9ER0dzWoAqjGWAFIKZmRnCw8PRoUMH9OrVCzt37hQ7EpHCuX79Ojp16oS9e/di69at8Pf357QAVQqWAVIYhoaGOHbsGEaMGAEfHx98/fXXvMkR0f9s27YNtra2AICLFy/C19dX5ESkSnijIlIo2tra8Pf3h6WlJT7//HMkJCRgw4YN0NLSEjsakSiys7MxdepUbN26FWPHjsWaNWtQs2ZNsWORimEZIIUjkUjw2WefwdzcHGPHjsWdO3ewf/9+GBoaih2NqFpdu3YNnp6eSEpKwvbt2zFy5EixI5GK4jQBKSwfHx+cPn0aMTExcHJyQkpKitiRiKrN1q1bYWtrC6lUipiYGBYBqlIsA6TQunbtioiICGRnZ8POzg6XL18WOxJRlcrOzoavry/GjBmDESNGICoqCi1bthQ7Fqk4lgFSeC1btoRMJoOZmRmcnZ1x9OhRsSMRVYk//vgDtra2OHDgAAICAvDrr7/y/ACqFiwDpBTq1auH4OBg9OrVC25ubli3bp3YkYgqjSAI8Pf3R6dOnaCpqYmYmBj4+PiIHYvUCMsAKY2aNWti//79mD59OqZMmYK5c+dCLpeLHYuoQrKysuDr64tx48bBx8cHUVFReP/998WORWqGqwlIqWhoaGDlypWwsLDAzJkzkZSUhICAANSoUUPsaERldvXqVXh6euLOnTvYsWMHRowYIXYkUlM8MkBKadq0aTh06BCOHz+O7t27IzU1VexIRKUmCAI2b96MTp06QUtLC5cuXWIRIFGxDJDSGjhwIEJCQpCYmAgHBwfcuHFD7EhE75SVlYWRI0di/PjxGDVqFKKiotCiRQuxY5GaYxkgpdaxY0fIZDLo6OjAwcEBoaGhYkcieqsrV66gY8eOCAoKwq5du7BhwwZOcZFCYBkgpWdubo6IiAhYW1ujZ8+e2L17t9iRiF4hCAJ+/fVX2NnZQVdXF5cuXcKwYcPEjkX0AssAqQQjIyOcOHEC3t7eGD58OL799lve5IgUQmZmJnx8fDBx4kT4+voiMjISzZs3FzsW0Su4moBUhra2NrZu3QpLS0ssXLgQCQkJWL9+PW9yRKKJj4+Hp6cn7t27h927d8Pb21vsSERvxCMDpFIkEgm++OILbNu2Ddu3b0f//v2RkZEhdixSM4IgYOPGjbCzs0ONGjVw+fJlFgFSaCwDpJJGjRqFEydOIDo6Gp07d8adO3fEjkRqIiMjA8OHD4efnx/Gjh0LmUyGZs2aiR2LqEQsA6SyunfvjoiICGRkZMDe3h5xcXFiRyIVFxcXh44dO+L333/H3r17sW7dOujq6oodi+idWAZIpX3wwQeQyWRo0KABunTpguPHj4sdiVSQIAjYsGED7O3toaenh8uXL8PT01PsWESlxjJAKu+9997D+fPn0b17dwwYMAC//PKL2JFIhWRkZGDYsGGYNGkSxo0bh4iICFhZWYkdi6hMuJqA1EKtWrVw8OBBzJ49G5MnT0ZCQgK+++47SKXsw1R+sbGx8PT0xMOHD7Fv3z4MHTpU7EhE5cK/hKQ2NDQ0sHr1aqxatQrLly+Ht7c3cnNzxY5FSkgQBKxfvx4ODg4wMDDA5cuXWQRIqbEMkNqZMWMGDh48iKNHj6Jnz554/Pix2JFIiWRkZMDb2xsfffQRJkyYwGkBUgksA6SW3N3dcf78edy+fRsODg64deuW2JFICcTGxqJ9+/Y4ceIEfvvtN6xZswY6OjpixyKqMJYBUludOnWCTCaDpqYmHBwcEB4eLnYkUlCCIGDdunWwt7eHkZERLl++jCFDhogdi6jSsAyQWmvatCkiIiLQpk0b9OjRA3v37hU7EimY9PR0eHp6YsqUKfDz80N4eDgsLS3FjkVUqbiagNRe7dq1ceLECYwfPx7e3t5ISkrCvHnzIJFIxI5GIrt06RI8PT2RlpaGAwcOwMPDQ+xIRFWCRwaIAOjo6GD79u344osvsGDBAkyaNAlFRUVixyKRCIKAtWvXwtHREcbGxrh8+TKLAKk0lgGi/5FIJFi8eDG2bNkCf39/DBgwAJmZmWLHomqWnp6OoUOHYtq0aZg8eTLCwsJgYWEhdiyiKsUyQPQfo0ePxokTJxAREYEuXbrg7t27YkeiahITE4P27dvjzJkzOHjwIFatWsXVAqQWWAaI3qBHjx4IDw/H06dPYW9vj/j4eLEjURUSBAFr1qyBo6MjTExMEBsbi0GDBokdi6jasAwQvUXr1q0hk8lQr149dO7cGSdOnBA7ElWBZ8+eYciQIZg+fTqmTJmCsLAwNG3aVOxYRNWKZYCoBPXr10dISAi6desGV1dXbNy4UexIVIkuXryI9u3b49y5czh06BBWrlwJbW1tsWMRVTuWAaJ30NPTQ2BgICZNmgQ/Pz988sknkMvlYseiChAEAatXr4aTkxPq1q2L2NhYuLu7ix2LSDS8zgBRKWhoaGDNmjWwtLTEnDlzkJiYiK1bt0JXV1fsaFRGT58+xdixYxEYGIhZs2bhu+++49EAUnssA0SlJJFIMGvWLDRp0gQjRoxAz549ERQUBBMTE7GjUSlFR0fDy8sLz549Q2BgINzc3MSORKQQOE1AVEYeHh4IDg7GzZs34eDggNu3b4sdid5BEASsWrUKnTt3Rr169RAbG8siQPQSlgGicrC3t0dkZCQkEgkcHBwQGRkpdiR6i6dPn2LQoEGYNWsWpk+fjgsXLsDc3FzsWEQKhWWAqJwsLS0RGRmJli1bwsXFBb/99pvYkeg/oqKiYGNjgwsXLiAoKAjLly/n+QFEb8AyQFQBxsbGOH36NAYPHgxPT08sW7YMgiCIHUvtCYKAFStWoHPnzqhfvz5iY2MxcOBAsWMRKSyeQEhUQTo6OtixYwcsLCwwb948JCQkYM2aNdDU5D8vMTx58gRjxozB4cOHMXfuXHz77bfQ0tISOxaRQuNfK6JKIJFIsGTJEpibm8PPzw8pKSnYu3cv9PT0xI6mVmQyGby8vJCVlYUjR47A1dVV7EhESoHTBESVaNy4cTh27BhCQ0Ph7OyMe/fuiR1JLQiCgB9//BFdunRBw4YNERsbyyJAVAYsA0SVrHfv3ggPD8ejR49gZ2eHq1evih1JpT158gRubm6YO3cuZs+ejZCQEJiZmYkdi0ipsAwQVYE2bdogKioKdevWhZOTE06dOiV2JJUUGRkJGxsbhIeH4+jRo/j+++95fgBRObAMEFWRBg0a4MKFC+jSpQv69euHzZs3ix1JZcjlcixfvhzOzs5o1KgR4uLi0L9/f7FjESktlgGiKqSnp4egoCBMmDAB48ePx2effcalhxWUlpYGNzc3fPzxx5gzZw7Onz+Pxo0bix2LSKlxNQFRFdPU1MS6detgaWmJjz/+GImJifD394eOjo7Y0ZROREQEvL29kZOTg99//x39+vUTOxKRSuCRAaJqIJFIMHfuXOzbtw8HDhxAr1698OTJE7FjKQ25XI5ly5bB2dkZZmZmiIuLYxEgqkQsA0TVaOjQoTh37hz+/PNPODo6IiEhQexICu/x48cYMGAA5s2bh3nz5uH8+fNo1KiR2LGIVAqnCYiqmaOjIyIjI9GvXz/Y29vj8OHDsLe3r9T3EAQBz549Q2pq6hsf6enpKCoqQnFxMZKSkhATE4NBgwZBU1MTGhoa0NbWhomJCUxNTd/4qFmzZqXmfZvw8HB4e3sjLy8Px48fR9++favlfYnUDcsAkQisrKwQGRkJd3d3uLi4YMeOHRg8eHC5xsrMzMSVK1cQHx+PuLg4xMXF4erVq8jLy3vldZqamqhbty5MTU1haGgILS0taGhovHhdbm4uiouLUVxcjPz8fERHRyM1NfWN0xl16tSBtbX1K48WLVpU2iWYn08LLFy4EI6Ojti9ezcaNmxYKWMT0eskQilObc7IyIChoSHS09NhYGBQHbmI1EJeXh7GjBmDvXv3Yvny5Zg1axYkEkmJ2+Tn5yM4OBiBgYE4e/Ysbt++DeDfnX2rVq3Qrl07tGvXDmZmZq98mjcyMoJUWvaZwcLCQjx+/PjFUYWHDx/i9u3bL4pHcnIygH/v0dCuXTv069cP7u7uaNu27Tu/lzd5/PgxRo0ahePHj+PTTz/F4sWLeZ8HonIq7f6b/8KIRKSrq4udO3eiadOmmDNnDhISErBq1arXdn4ZGRk4duwYAgMDcezYMWRmZsLCwgL9+vVDhw4dYG1tjZYtW1bJCgUtLS3Ur18f9evXf+PzT58+xZUrVxAXF4fIyEisWLECixYtQtOmTeHm5gZ3d3d07twZGhoa73yvsLAweHt7Iz8/HydOnECfPn0q+9shojcRSiE9PV0AIKSnp5fm5URUDhs3bhQ0NDQEV1dXITMzUxAEQXjw4IEwf/58QV9fXwAgdOjQQViyZIlw9epVQS6Xi5z4zfLz84WTJ08KkydPFho0aCAAEJo3by74+/sLBQUFb9ymuLhYWLp0qaChoSE4OzsLd+/erebURKqptPtvlgEiBXLixAlBT09PaN26tTB27FhBV1dX0NPTE+bPny8kJyeLHa/MiouLhbCwMMHNzU0AIJiZmQlr164VcnJyXrwmNTVV6Nu3ryCRSITPPvtMKCwsFDExkWop7f6b5wwQKRBBELBw4UIsXboUUqkUkyZNwtdff43atWuLHa3Crl69iqVLl2Lv3r1o2rQp9uzZg9zcXHh7e6OwsBA7duxA7969xY5JpFJKu//mdQaIFMSTJ08wePBgLF26FKNHj8YHH3yAHTt24NKlS2JHqxRt2rTBrl27cP36dRgbG8Pe3h5du3aFlZUV4uLiWASIRMQyQKQA/vrrL1hbW+P8+fMIDAzEli1bEBERAUdHR3z44YfYsmWL2BErTe3atWFoaIji4mIIggADAwPUqVNH7FhEao1lgEhkaWlpcHV1hb6+PuLj4+Hm5gYA0NfXx5EjRzB27FiMHTsWX3zxhdLf5CgkJATW1taIj4/HqVOncPToUZw8eRJTpkxR+u+NSJlxaSGRiAoKCjBkyBA8e/YM0dHRr919T1NTE7/88gssLS0xf/58JCYmYtOmTUp3kyO5XI6lS5fiiy++gLOzM3bu3IkGDRoAADZu3IgxY8agVatWmDlzprhBidQUywCRiNauXYuwsDCcPXsWFhYWb3yNRCLBvHnz0KRJE/j6+uLu3bs4ePCg0pxUmJqaCh8fH5w5cwaff/45vvjii1euOTB69Ghcu3YNc+bMwYABA2BpaSliWiL1xGkCIpEIgoBffvkFXl5ecHZ2fufrvby8cPbsWVy9ehWOjo5ITEyshpQVc/78+VemBRYvXvzGiw999dVXMDAwwKZNm0RISUQsA0QiCQkJwa1btzBhwoRSb+Pk5ITIyEgUFRXB3t4e0dHRVZiw/IqLi7FkyRL06NED77//PuLi4tCzZ8+3vr5GjRoYOXIk/P39UVBQUI1JiQhgGSASTVhYGIyNjUt1VOBlzZo1Q2RkJKysrNCtWzccOnSoihKWz8OHD9G3b198+eWX+Pzzz3H69Om3Xsr4ZR4eHkhNTcWtW7eqISURvYxlgEgkgiBAR0enXDfzqVOnDs6ePQtXV1cMHjwYq1atqvyA5RAcHAxra2tcvXoVp0+fxqJFi0p1TwLg3/s0AOCqAiIRsAwQiUgul5d7W11dXezZswcff/wxZs2ahRkzZqC4uLgS05VecXExvvrqK/Ts2RMffPAB4uLi0KNHjzKNUZGfBRFVDFcTEImkRYsWePjwIf7880+0bNmyXGNIpVJ8//33aNq0KaZMmYKkpCTs2rULtWrVKtX22flFSErLRkGRHNqaUpib1EItnbL9WXjw4AFGjBiB4OBgLFq0CAsXLiz10YCXBQcHo0aNGjAzMyvztkRUMbw3AZFI8vPz0ahRI4waNQo//vhjhcc7duwYPD090bJlSxw5cgTvvffeG19362EmdkalIPhGKlKe5ODlPwASAGbGNeHSwhQj7MzQrJ5+ie957tw5DB8+HBKJBLt27YKLi0u5ssvlclhaWsLFxQX+/v7lGoOIXsd7ExApOB0dHfj6+sLf3x/JyckVHq9fv34IDQ3FvXv3YG9vj+vXr7/y/J0nORi5OQq9Vl1AQFQykv9TBABAAJD8JAcBUcnoteoCRm6Owp0nOa+9V3FxMRYvXoyePXuiTZs2iIuLK3cRAIDt27cjKSmpTCsriKjysAwQiWj+/PkwMjLCwIEDkZWVVeHxbGxsEBUVBUNDQzg6OuLcuXMAgD0XU9BzZQgiEtIAAMXykg8IPn8+IiENPVeGYM/FlBfPPXjwAL1798ZXX32FxYsX48SJE6hXr165M0dFRWHSpEkYNWoU7O3tyz0OEZUfpwmIRPbHH3/AwcEBPXr0wL59+6CtrV3hMTMyMjB06FAEBwdj5Hc7cfZRzQqPObd3c7SUp2DEiBGQSCTYvXs3unXrVqExExIS4OjoCCsrK5w9e1bpLrNMpOhKu/9mGSBSAL///jsGDRqEdu3aYffu3bCysqrwmIWFhRg463v8qW9TCQn/lXb8JziYCggICKjQ0QAAOHDgAMaPHw8TExNERETA1NS0klIS0XM8Z4BIifTv3x+RkZF4+vQp2rdvj927d1d4zAeZhUio3QF47cyAf8kLcvHkzEbcXTsKycsG4Z7/NGRfD3nreIIgwPTDqfh118EKFYHc3Fx89NFHGDJkCHr06IGLFy+yCBCJjGWASEF06NABly9fxsCBAzF8+HB069YNp0+fLvdFeD49dBVFcgH/rhF43aOD3yL76lkYOg1DPc/F0KnfDI8PL0P2tfNvfL1EIgGkGvgs6Fq58uTm5mLNmjVo0aIFtmzZgl9++QW//fab0txwiUiVsQwQKRADAwMEBAQgKCgI2dnZ6N27N+zs7BAUFFSmi/LcepiJ0NuP33qiYO7fF5GXFAvjPh9B3+ZD6DZpC5MPp0PX3AZPg/0hyN988aJiuYDQ249xOzWz1FkyMjLw/fffw9zcHLNmzULXrl0RFxcHPz+/cl19kYgqH8sAkYKRSCQYOHAgoqOjcfLkSdSoUQPu7u5o0qQJpk2bhrNnz6KwsLDEMXZGpUBD+vYdbc7NSEi0a6Dm+51f+bpe254oznqC/Hs337qthlSCHbKUtz4PAE+ePEFAQAAGDx6MBg0a4IsvvoC7uztu3ryJgIAAtGjRosTtiah6sQwQKSiJRILevXsjJCQEERER8PDwwOHDh9GzZ0+YmprCx8cHu3fvxo0bN167DHHwjdQSlw8WPEqGlkkjSKSvXilQq645AKDw8duve1AsFxB8M/WVr+Xk5CA6OhqrV69Gjx49YGpqilGjRuGff/7BwoUL8ffff2PDhg2wsLAo40+BiKoDL0dMpAQcHBzg4OCAVatWIS4uDoGBgQgMDMTOnTsBADVr1kTbtm1hbW2N99tYI/lJoxLHk+dmQtPo9SsUSmvo/+/5jBK3T07LwdffLcO1+MuIi4vDzZs3IZfLoaWlhe7du2Pt2rUYOHAgGjRoUM7vmIiqE8sAkRKRSCSwsbGBjY0NFi9ejEePHiE+Ph7x8fGIi4tDWFgYtgadQT3fVaUZrKQn37n58g3b0LqhEXr06IE5c+bA2toarVq1Qo0aNUr9/RCRYmAZIFJidevWRc+ePdGzZ88XX4u6/RBem2NK3E5aQ/+Nn/7luZkvnn+Xc+cvoH0T4zImJiJFxHMGiFSMXk3dd75Gu645CtPuvrZqoPBREgBAq06Td46ho1X2OxMSkWJiGSBSMeYmtd55kL9mcwcIBbnIuRH+ytez/jgHDT1j6DRoXuL2kv+9DxGpBk4TEKmYWjqaMDOuieQ33G3wuRqWHaFrboMnJ9dBnp8DrdoNkH09BHkJl2AyYM5rqwz+y8ykJmrp8M8Hkargv2YiFeTSwhQBUcklLi+s6/EpnoVsR3roThTnZULLuBHqDPwYtT7oWuLYGlIJXJrz8sFEqoRlgEgFjbAzw9bIpBJfI9WuAeNefjDu5VemsYvlAnzszSqQjogUDc8ZIFJBzerpo4tVnUr/B64hlaCLVR1Ymb57tQERKQ+WASIVVFRUBKPbx1FUmA+U80ZHrxEEaEol+HZQm8oZj4gUBssAkYq5c+cOunXrhvXLvkYPw7R3XFyoDCQSIGYfhKzHlTMeESkMlgEiFXLs2DHY2NggJSUFFy5cwJbPJ2Bu75KXCZaWr01tFPx1HnZ2drh8+XKljElEioFlgEgFFBYWYv78+ejfvz/s7e0RGxsLR0dHAMBUl2b4zqMNdDSlJd7J8E00pBLoaErxvUcbLPZ0hEwmg5mZGZydnXH06NGq+FaISAQsA0RK7vm0wIoVK7Bs2TIcPnwYJiYmr7zG29YMZ2Z1haPFv19/Vyl4/ryjhQnOzOoKL9t/Vw/Uq1cPwcHB6NWrF9zc3LBu3boq+I6IqLpxaSGREvv9998xatQo1KpVCxcuXICDg8NbX9vYuCYCxtnh1sNM7IxKQfDNVKSk5eDl0wsl+PeCQi7NTeFjb/bGVQM1a9bE/v37MXfuXEyZMgWJiYn4/vvvIZXyswWRspIIwrtPNc7IyIChoSHS09NhYGBQHbmIqASFhYVYuHAhli1bhgEDBmDr1q0wNi77TYOy84uQlJaNgiI5tDWlMDepVaYrC/7000+YOXMmPDw8EBAQwDsWEimY0u6/eWSASMmkpKTA29sbFy9exI8//ohZs2ZBUs4VA7V0NNGqgWG5s0yfPh1NmjTBsGHD0L17dxw+fBh169Yt93hEJA4e1yNSIkePHoWNjQ3++ecfhIaGYvbs2eUuApXFzc0NISEhSExMhL29PW7evClqHiIqO5YBIiVQWFiIjz/+GAMGDICTkxNiY2Nhb28vdqwXbG1tIZPJoKOjAwcHB4SGhoodiYjKgGWASMElJyfD2dkZq1atwooVKxAUFFSu8wOqmrm5OcLDw9G2bVv07NkTe/bsETsSEZUSywCRAjty5AhsbGxw//59hIWFVej8gOpQu3ZtnDx5El5eXhg2bBi+++47lOIcZSISGcsAkQIqLCzE3LlzMXDgQDg7OyM2NhZ2dnZixyoVbW1tbNu2DV9++SU++eQT+Pn5obCwUOxYRFQCriYgUjDJycnw9vbGpUuXsHLlSsyYMUOhjwa8iUQiwaJFi9C0aVOMHz8eKSkp2LdvH5cmEykoHhkgUiCHDx+GjY0NHjx4gLCwMMycOVPpisDLfH19ceLECchkMnTp0gV3794VOxIRvQHLAJECKCgowJw5c+Dm5oauXbvi8uXL6NSpk9ixKkWPHj0QHh6OZ8+ewc7ODnFxcWJHIqL/YBkgEllSUhK6dOmCNWvWYNWqVTh48CBq164tdqxK1apVK0RFRaF+/fro0qULjh8/LnYkInoJywCRiIKCgmBjY4PU1FSEh4cr5fkBpfXee+8hJCQELi4uGDBgADZu3Ch2JCL6H5YBIhEUFBRg1qxZcHd3h4uLC2JjY2Frayt2rCpXq1YtHDp0CJMnT4afnx8WLFgAuVwudiwitcfVBETVLCkpCZ6enoiLi8Pq1asxbdo0lT0a8CYaGhr46aefYGlpidmzZyMxMRHbtm2Drq6u2NGI1BaPDBBVo8DAQNjY2ODx48cIDw/H9OnT1aoIPCeRSDBz5kzs378fhw8fRs+ePfH48WOxYxGpLZYBompQUFCAmTNnYtCgQejevTsuX76sFtMC7+Lh4YHz58/j5s2bcHR0xO3bt8WORKSWWAaIqlhiYiI6d+6M9evXY82aNdi/fz+MjIzEjqUw7OzsIJPJIJVKYW9vj4iICLEjEakdlgGiKnTo0KFXpgWmTp2qltMC72JhYYGIiAi0atUK3bt3x2+//SZ2JCK1wjJAVAXy8/MxY8YMeHh4oGfPnrh8+TI6duwodiyFZmxsjFOnTmHw4MHw9PTEsmXLeJMjomrC1QRElSwhIQGenp64evUq1q5di48++ohHA0pJR0cHO3bsgIWFBebNm4eEhASsWbMGmpr8U0VUlfgvjKgSHThwAGPHjkWdOnUQERGBDh06iB1J6UgkEixZsgTm5ubw8/NDSkoK9u7dCz09PbGjEaksThMQVYL8/HxMnz4dQ4YMQZ8+fXD58mUWgQoaN24cjh07htDQUDg7O+PevXtiRyJSWSwDRBWUkJAAJycnbNiwAT///DP27t0LQ0NDsWOphN69eyM8PByPHj2CnZ0drl69KnYkIpXEMkBUAfv374eNjQ2ePXuGyMhInh9QBdq0aYOoqCjUqVMHTk5OOH36tNiRiFQOywBROeTn52PatGkYOnToi2mB9u3bix1LZTVo0AAXLlxAly5d0K9fP2zevFnsSEQqhWWAqIz+/vtvODo6YuPGjVi3bh327t0LAwMDsWOpPH19fQQFBWH8+PEYP348PvvsMy49JKokXE1AVAa//fYbxo8fD1NTU8hkMtjY2IgdSa1oampi3bp1L5YeJiYmwt/fHzo6OmJHI1JqPDJAVAp5eXmYMmUKPD098eGHH+LSpUssAiKRSCT4+OOPsW/fPhw4cAC9evXCkydPxI5FpNRYBoje4fbt23B0dMTmzZuxfv167N69m9MCCmDo0KE4d+4crl+/DkdHRyQkJIgdiUhpsQwQlWDfvn1o3749srKyIJPJMGnSJK4WUCCOjo6QyWSQy+Wwt7eHTCYTOxKRUmIZIHqDvLw8fPTRR/Dy8kL//v0RExMDa2trsWPRG1hZWSEyMhLNmzeHi4sLDh48KHYkIqXDMkD0H7du3YKDgwP8/f2xYcMG7Nq1i9MCCs7ExARnzpyBm5sbhgwZgpUrV3KlAVEZcDUB0Uv27t2LCRMm4L333oNMJuPRACWiq6uLXbt2oWnTppg9ezYSEhKwatUqaGhoiB2NSOHxyAAR/p0WmDx5Mry9veHq6opLly6xCCghqVSKpUuXYsOGDVi/fj0GDRqE7OxssWMRKTyWAVJ7N2/ehL29PbZs2YKNGzdi586d0NfXFzsWVcDEiRNx9OhRBAcHo2vXrrh//77YkYgUGssAqbXdu3ejQ4cOyM3NRVRUFCZMmMDVAiqib9++CAsLw4MHD2Bvb49r166JHYlIYbEMkFrKzc3FpEmTMHz4cAwcOBAxMTFo166d2LGokrVr1w4ymQxGRkZwcnLCuXPnxI5EpJBYBkjt3Lx5Ew4ODti2bRt+/fVX7Nixg9MCKqxRo0YIDQ2Fvb09+vTpg23btokdiUjhsAyQWvnvtMD48eM5LaAGDAwMcOTIEYwePRqjR4/GokWLuPSQ6CUsA6QWcnNz4efnh+HDh8PNzQ0xMTFo27at2LGoGmlpaWHjxo1YunQpFi9eDF9fXxQUFIgdi0gh8DoDpPJu3LgBT09P3Lx5E5s2bcLYsWN5NEBNSSQSLFiwAObm5vD19cXdu3dx4MAB1K5dW+xoRKLikQFSaTt37kSHDh1QUFCA6OhojBs3jkWA4O3tjTNnziA+Ph5OTk5ISkoSOxKRqFgGSCXl5uZiwoQJ8PHxwaBBg3Dx4kW0adNG7FikQLp06YKIiAjk5+fDzs4OFy9eFDsSkWhYBkjl/PXXX7Czs8POnTuxefNmbN++HXp6emLHIgXUokULyGQyWFhYoGvXrggKChI7EpEoWAZIpezcuRMdO3ZEYWEhoqOjeX4AvVPdunVx7tw59OvXD4MGDcJPP/0kdiSiascyQCohJycH48ePh4+PDzw8PHDx4kW0bt1a7FikJGrUqIF9+/Zhzpw5mDFjBmbOnIni4mKxYxFVG64mIKX3119/YejQofj777/h7++P0aNH82gAlZlUKsWyZcvQtGlTTJs2DcnJydi5cydq1qwpdjSiKscjA6TUAgIC0LFjRxQXF+PixYsYM2YMiwBVyEcffYSgoCCcPn0a3bp1w8OHD8WORFTlWAZIKeXk5GDcuHEYNWoUhgwZgosXL6JVq1ZixyIV4erqigsXLuDu3buwt7fHn3/+KXYkoirFMkBK5/r16+jUqRN2796NLVu2YOvWrahVq5bYsUjFtG/fHjKZDHp6enB0dMT58+fFjkRUZVgGSKls374dtra2EAQBFy9exOjRo8WORCrMzMwMYWFh6NixI3r37o0dO3aIHYmoSrAMkFLIycnB2LFj4evrCy8vL0RHR3NagKqFoaEhjh07Bh8fH4wcORJLlizhTY5I5XA1ASm869evw9PTE4mJidi6dSt8fX3FjkRqRktLC5s3b4aFhQU+//xzJCYmYsOGDdDS0hI7GlGl4JEBUmjbtm17ZVqARYDEIpFI8NlnnyEgIAA7duxAv379kJ6eLnYsokrBMkAKKTs7G2PGjMHo0aNfTAt88MEHYscigo+PD06fPo2YmBg4OTkhJSVF7EhEFcYyQArn2rVr6NSpE/bt24dt27bB39+fqwVIoXTt2hURERHIzs6GnZ0dLl++LHYkogphGSCFsnXrVtja2kIikSAmJgajRo0SOxLRG7Vs2RIymQxmZmZwdnbG0aNHxY5EVG4sA6QQsrOz4evrizFjxmD48OGIjo5Gy5YtxY5FVKJ69eohODgYvXr1gpubG9atWyd2JKJyYRkg0f3xxx+wtbXFgQMHEBAQgE2bNvF68KQ0atasif3792P69OmYMmUK5s6dC7lcLnYsojLh0kISjSAI2Lp1K6ZMmQJLS0vExMTg/fffFzsWUZlpaGhg5cqVaNq0KWbOnImkpCQEBASgRo0aYkcjKhUeGSBRZGVlwdfXF2PHjsWIESMQFRXFIkBKb/r06Th06BCOHTuG7t27IzU1VexIRKXCMkDV7vm0wMGDBxEQEIBff/2V0wKkMtzc3BASEoLExEQ4ODjgxo0bYkcieieWAao2giBg8+bN6NSpE7S0tBATEwMfHx+xYxFVOltbW8hkMujo6MDBwQGhoaFiRyIqEcsAVYusrCyMGjUK48ePh4+PD6cFSOWZm5sjPDwc7dq1Q8+ePbF7926xIxG9FcsAVbmrV6+iY8eOCAwMxM6dO7Fx40aeWEVqoXbt2jh58iS8vLwwfPhwLF26lDc5IoXE1QRUZZ5PC0ybNg3NmzdHTEwMWrRoIXYsomqlra2Nbdu2wcLCAp9++ikSEhKwbt063uSIFAqPDFCVyMrKwsiRIzFhwgSMGjUKMpmMRYDUlkQiwaJFi7B161Zs3boVrq6uyMjIEDsW0QssA1Tprly5go4dOyIoKAi7du3Chg0bOC1ABMDX1xcnTpxAVFQUOnfujLt374odiQgAywBVIkEQ8Ouvv8LOzg66urq4dOkShg0bJnYsIoXSo0cPhIeHIz09HXZ2doiLixM7EhHLAFWOzMxM+Pj4YOLEifD19UVkZCSaN28udiwihdSqVStERUWhfv366NKlC44fPy52JFJzLANUYfHx8ejYsSMOHz6M3bt345dffuG0ANE7vPfeewgJCYGLiwsGDBiAjRs3ih2J1BjLAJWbIAjYuHEj7OzsUKNGDVy+fBne3t5ixyJSGrVq1cKhQ4cwefJk+Pn5YcGCBbzJEYmCSwupXDIyMuDn54c9e/Zg8uTJWLFiBXR1dcWORaR0NDQ08NNPP8HS0hKzZ89GYmIitm3bxn9PVK1YBqjM4uLi4OnpiQcPHmDPnj3w8vISOxKRUpNIJJg5cybMzMwwYsQI9OzZE4GBgahTp47Y0UhNcJqASk0QBGzYsAH29vaoVasWLl26xCJAVIk8PDxw/vx53Lx5E46Ojrh9+7bYkUhNsAxQqWRkZGDYsGGYNGkSxo0bh8jISDRr1kzsWEQqx87ODjKZDFKpFPb29oiIiBA7EqkBlgF6p9jYWHTo0AHHjh3D3r178fPPP3M+k6gKWVhYICIiAq1atUL37t2xb98+sSORimMZoLcSBAHr16+Hg4MD9PX1cfnyZXh6eoodi0gtGBsb49SpUxg8eDC8vLzwww8/8CZHVGV4AiG9UUZGBiZMmIB9+/ZhypQpWL58OY8GEFUzHR0d7NixAxYWFpg/fz4SExOxZs0aaGryTzdVLv5G0WtiY2MxdOhQPHr0CL/99huGDBkidiQitSWRSLBkyRKYm5vDz88PycnJ2Lt3L/T19cWORiqE0wT0giAIWLduHezt7WFoaIjLly+zCBApiHHjxuHYsWMICwuDs7Mz/vnnH7EjkQphGSAAQHp6Ory8vDBlyhRMnDgRERERsLS0FDsWEb2kd+/eCA8Px+PHj2Fvb48rV66IHYlUBMsA4fLly+jQoQNOnjyJ3377DWvWrIGOjo7YsYjoDdq0aYOoqCjUqVMHnTt3xqlTp8SORCqAZUCNCYKAn3/+GQ4ODjAyMuK0AJGSaNCgAS5cuIAuXbqgX79+2Lx5s9iRSMmxDKip9PR0eHp6YurUqfDz80N4eDinBYiUiL6+PoKCgjBhwgSMHz8en332GZceUrlxNYEaunTpEjw9PZGWloYDBw7Aw8ND7EhEVA6amppYt24dLCwsMG/ePCQmJsLf35/TfFRmPDKgRgRBwNq1a+Ho6AhjY2NcvnyZRYBIyUkkEnz88cfYt28fDhw4gF69euHJkydixyIlwzKgJp49e4YhQ4Zg2rRpmDx5MsLCwmBhYSF2LCKqJEOHDsW5c+fw559/wtHREQkJCWJHIiXCMqAGYmJi0L59e5w9exYHDhzAqlWreBiRSAU5OjoiMjIScrkc9vb2kMlkYkciJcEyoMIEQcCaNWvg6OiIOnXqIDY2ltMCRCrOysoKkZGRaN68OVxcXHDw4EGxI5ESYBlQUc+nBaZPn44pU6YgLCwMTZs2FTsWEVUDExMTnDlzBm5ubhgyZAhWrlzJlQZUIq4mUEEXL16El5cXnj59ikOHDsHd3V3sSERUzXR1dbFr1y40bdoUs2fPRkJCAlatWgUNDQ2xo5EC4pEBFSIIAlavXg0nJyfUrVsXsbGxLAJEakwqlWLp0qXYsGED1q9fD3d3d2RlZYkdixQQy4CKePr0KTw8PDBz5kxMnToVoaGhMDc3FzsWESmAiRMn4ujRozh//jy6du2K+/fvix2JFAzLgAqIjo5G+/btcf78eQQGBmLFihXQ1tYWOxYRKZC+ffsiLCwMDx8+hL29Pa5duyZ2JFIgLANKTBAErFq1Cp07d4apqSliY2Ph5uYmdiwiUlDt2rWDTCaDkZERHB0dcfbsWbEjkYJgGVBST58+xaBBgzBr1ixMmzaN0wJEVCqNGjVCaGgoHBwc0LdvX2zdulXsSKQAuJpACUVHR8PT0xMZGRkICgrCwIEDxY5ERErEwMAAR44cwUcffYQxY8YgMTERixYtgkQiETsaiYRlQIk8nxaYN28eOnTogJCQEDRp0kTsWESkhLS0tLBx40ZYWlrik08+QWJiIjZt2sTzjdQUpwmUxJMnT+Du7o7Zs2djxowZuHDhAosAEVWIRCLBggULsHv3buzduxd9+vTB06dPxY5FIlC7IwNFRUVIS0vD48ePUVBQgKKiIhQXF+P+/fuoU6cOtLW1oampCU1NTejr68PU1BS1atUS9fCZTCaDl5cXMjMzcfjwYQwYMEC0LESkery9vdGwYUO4u7vDyckJx44d4zlIakYilOIalRkZGTA0NER6ejoMDAyqI1e5PX36FPHx8YiLi8Pff/+N1NTUVx5paWllvixnjRo1YGpqClNTU9StWxempqaoX78+WrduDWtrazRv3hyampXfqwRBwMqVKzF//nx07NgRe/fuhZmZWaW/DxERANy4cQP9+vVDVlYWjh49CltbW7EjUQWVdv+t1GUgOzsb586dQ0xMDOLi4hAXF4eUlBQA/16K08rKCvXq1XuxI3/5UadOHejo6EBTUxMaGhrIy8tDjRo1UFxcjOLiYhQVFSEzM/O1MvH8cffuXdy5c+fFez0vBtbW1ujatStatWpVoaMJT548wejRo3HkyBHMnTsX3377LbS0tCrl50ZE9DaPHj3CwIEDER8fj927d3O5spJT2TLw6NEjHDlyBEFBQTh16hTy8vJgamoKGxsbtGvX7sUOuVmzZlXyaf1lT58+xZUrV14Ukbi4OFy7dg2FhYWwtLSEu7s73N3d4eDgUKbrgUdGRsLb2xtZWVnYtm0bXF1dq/C7ICJ6VW5uLkaOHImDBw9i1apVmD59epW/pyAIePbs2YsPXLm5uS+mcR88eICCggI0btz4xTSurq7uiw93RkZGkEp5CtybqFQZKC4uxr59+7B+/XqEh4dDEAQ4OTnB3d0dbm5usLKyqvZMb5Ofn4/g4GAEBgYiKCgIDx48QN26dTF48GDMmTOnxKyCIGDFihVYsGABbG1tsWfPHk4LEJEo5HI55s+fj+XLl2PGjBn48ccfK3STI0EQcP/+/RfTuH/++ScePnz4Yuf/6NEjFBYWlmtsTU1N1K1b98U0rqmpKd5///0XHw4bNWqktssmVaIMFBQUICAgAN999x1u376Nnj17wtvbG66urqhXr1615SgvuVyO6OhoHDp0CNu2bcOjR4/g7e2NTz75BK1bt37ltWlpaRg9ejSOHj2KefPm4euvv+a0ABGJbt26dZg2bRoGDhyInTt3ombNmqXaLjs7G6dOnUJkZOSLI6ePHj0C8O91Dlq3bo369eu/cRq3bt260NPTg4aGxoupXLlcDkEQXhwtyMnJwaNHj944jfvgwQNcu3YNT548AQDUrl37RTGws7ND3759YWhoWGU/M0VS6v23UArp6ekCACE9Pb00L68wuVwubNu2TWjcuLEAQPDw8BBiYmKq5b2rSk5OjrB27VrBzMxMACAMGjRISElJEQRBECIiIoTGjRsLxsbGwtGjR0VOSkT0qiNHjgi1atUSbG1thQcPHrz1dampqcLmzZuFAQMGCLq6ugIAoUmTJoKbm5vw5ZdfCocOHRISExMFuVxe5Znlcrlw584d4ciRI8KSJUuEwYMHC5aWlgIAQUtLS+jTp4+wfv164Z9//qnyLGIq7f5b4cpAenq6MGzYMAGA4OnpKVy7dq3K37M65efnC/7+/kKjRo2E2rVrC76+voKmpqbg6Oj4ohwQESmaS5cuCfXr1xfMzc2F69evv/i6XC4XgoKChG7duglSqVSQSCRC586dheXLlwu3bt0SMfGbJScnCz/99JPQvXt3QUNDQwAg2NvbC7t37xaKiorEjlfplLIMxMTECJaWloK+vr6wa9euKn0vsd26dUuoV6+eAEBo3769kJmZKXYkIqISJScnC61btxaMjIyEM2fOCLt37xbatGkjABCcnJyETZs2CQ8fPhQ7ZqmlpaUJ27dvF3r16iUAEJo1ayZs3rxZyM/PFztapVG6MnD16lVBT09P6Nixo3D79u0qex9FEB4eLjRq1EgwNjYWJk+eLGhrawve3t7VcuiMiKginj17JtjY2AgABABCnz59hAsXLogdq8Kio6OFQYMGCQCExo0bC/v37xc7UqUo7f5bIdZiPHr0CAMGDIClpSWCg4NhaWkpdqQqIZfL8cMPP8DZ2RlNmjRBXFwc1q1bh507d2LPnj34+uuvxY5IRPRWhYWF+PbbbxEbG4v69esDABwdHdG5c2eRk1Wcra0tDh48iD/++APt27fHkCFDMHnyZOTm5oodrXpUZrMoD7lcLnTt2lUwNTUVkpOTK318RfHo0SOhX79+AgBhwYIFQkFBwSvPf/XVVwIA4fDhwyIlJCJ6u6SkJMHOzk7Q1NQUfvjhB6GoqEhYsmSJAEAYPXq0Sh1al8vlwoYNGwRdXV2hdevWwp9//il2pHJTmmmCkJAQAYBw7NixSh9bUYSFhQmNGjUSTExM3vp9yuVyoXv37oKdnV01pyMiKllaWprQrFkzoUmTJoJMJnvluYCAAEFLS0vo0aOH8PTpU3ECVpGrV68KLVu2FBo0aCDcvXtX7DjlojTTBBs3boSVlRX69u0rdpRKJ5fL8f3336Nr164wNzdHXFwcPvzwwze+ViKRYMqUKYiKisLVq1erOSkR0ZsVFhbC09MTaWlpOHv2LOzs7F553sfHB6dPn8alS5fQuXNnJCcni5S08rVu3Rrnzp2DVCqFm5sbcnJyxI5UZUQtA/n5+di/fz/Gjx+vcleHevz4MVxdXbFgwQLMmzcPwcHBaNSoUYnbDBgwAKamptixY0c1pSQiKtnChQsREhKCAwcOvPV8rq5duyIiIgLZ2dmwt7fHpUuXqjll1Xnvvfdw+PBh/Pnnn5gyZYrYcaqMqGUgLy8P+fn5KnfCYFhYGKytrXHx4kUcP34c3377banuk6ClpYUmTZrwfuJEpBAyMjLw888/45NPPkG3bt1KfG3Lli0hk8lgZmYGZ2dnHD16tHpCVgMbGxt89913CAgIwD///CN2nCoh+jSBKpHL5fjuu+/QrVs3WFhYIC4uTiWnP4hIPezatQt5eXmYOHFiqV5fr149BAcHo3fv3nBzc8O6deuqOGH1GTVqFHR0dLBlyxaxo1QJUcuAtrY2pFIp7t+/L2aMSvHo0SP0798fn3zyCebPn49z586hYcOGZRpDLpfjwYMHpb72NxFRVdq3bx/69OnzzinOl9WsWRP79+/H9OnTMWXKFMydOxdyubwKU1YPQ0NDDBkyBHv37hU7SpUQtQzUqFED/fv3x9atW8WMgez8Ily7l47YlKe4di8d2flFZdo+NDQU1tbWiImJwYkTJ/DNN9+U6/bJ586dw507dzB06NAyb0tEVNmysrLKVASe09DQwMqVK7F69WqsWLECnp6eKrFev3HjxsjKyhI7RpUo+x6rkk2cOBEDBgxAdHQ0OnXqVG3ve+thJnZGpSD4RipSnuTg5Vs3SgCYGdeESwtTjLAzQ7N6+m8c4/lqgc8++wydO3fGrl27ynw04GW//PILWrZsCUdHx3KPQUSkKKZPn44mTZpg2LBh6N69O4KCgmBqaip2LHoD0c8Z6Nu3Lz744AOMHDmyWk6cu/MkByM3R6HXqgsIiEpG8n+KAPDvNTaTn+QgICoZvVZdwMjNUbjz5NUlJY8ePUK/fv2wcOFCfPrppzh79myFisDWrVtx4MABzJ07V+VWVhCRctLX16/wUkE3NzeEhIQgMTERDg4OuHHjRiWlq35JSUkl3wZYiYleBjQ1NREUFITHjx9j6NChKCwsrLL32nMxBT1XhiAiIQ0AUCz/bw141fPnIxLS0HNlCPZcTAEAXLhwAdbW1rh8+TJOnDiBJUuWlGta4LmwsDBMnDgREyZMwJgxY8o9DhFRZRo2bBhOnz5d4UJga2sLmUwGHR0dODg4IDQ0tNTbVnQat7I8ffoUBw4cwLBhw0R5/6omEQSh5D0i/l1eYmhoiPT09CprRc/PQO3Xrx+2bt2K2rVrV+r4a4NvYfmpmxUep6P2fRz62g9dunTBrl270KBBgwqN9/vvv2PkyJFo27YtTp06BW1t7QpnJCKqDFlZWahfvz6mTp2KpUuXVni8p0+fwsPDAxEREdi6detbd6yVMY1b2VasWIH58+fjzp07eO+996rlPStDafffoh8ZeM7FxQUHDx58cTJeeHh4pY2952JKpRQBAIgpqA/PT1bhzJkzFSoCBQUFmDNnDlxdXeHg4ICDBw+yCBCRQtHT08OMGTOwbNkynDp1qsLj1a5dGydPnoSXlxeGDx+OpUuX4uXPo5U1jVvZoqKi8Omnn2Ls2LFKVQTKQmGODDyXkpKC4cOHQyaTYd68eZgzZw5MTEzKPd6dJznouTIE+UWvL20peJiAZxe2o+BRMuQ56ZBoakPTuCH027tCr7XLW8fU0ZTizKyuaGxcviWAFy5cwOzZs3HlyhV8//33mDlzJs8TICKFVFxcjAEDBiAiIgIymQzvv/9+hccUBAGLFy/G4sWLMX78eKxbtw4H4u7jy8PXUCQX3jmF+zINqQSaUgkWD2wFb1uzCmf7rzt37qBTp05o2rQpzp07B11d3Up/j6pU2v23wpUBACgqKsKSJUuwbNkySKVSTJ48GXPmzClXIxu5OQoRCWlv/OXKS76C7D8vQKdRK2jqm0BemIfsa+eR8+cFGHbxgZGT9xvH1JBK4GhhgoBxdm98/k0EQcDJkyfxzTffICwsDG3btoW/vz86dOhQ5u+JiKg6paenw8HBAenp6di9ezecnZ0rZdxt27Zh/PjxaDd8Hh43qPgqqrm9m2OqS7NKSPavS5cuwcvLCwUFBbh48SLq1atXaWNXF6WbJniZpqYmFi9ejKSkJEyfPh0bN26Eubk5Jk2ahODgYBQVle4EklsPMxF6+/FbW6Zuk7Yw6TsVeq1doNukLWpadUJdt3nQbtACWfEn3zpusVxA6O3HuJ2a+c4MDx48wK+//gpbW1t8+OGHKCwsxOHDhxEXF8ciQERKwdDQEGfOnEGzZs3g4uKCxYsXo7i4uMLj+vr6YuHm3yulCADA8lM3sfd/J3pXhCAIWLVqFRwcHGBkZITg4GClLAJloZBl4DlTU1N8++23SE5OxmeffYajR4+ie/fuqFevHnx9fXHo0CFkZ2e/dfudUSnQkJb98LtGDQNIJCX/aDSkEuyQvfmX7ubNm/jhhx/g6OiIBg0aYNKkSTA2NsaZM2cQGRmJAQMGcFqAiJRKgwYNcPbsWXzxxRf46quv0LVr1zKtCniTO09ysPvWm0tFblI8Hv++Cv9snISUHwfj7tpRSN2/BPkPbpc45heHr1XoHIJLly7hww8/xKxZszB16lRERESo3P1z3kQhpwneRhAEXLp0CUFBQQgMDMQff/wBXV1d2NjYwNraGu3atYO1tTXatGmDmjVrouuyYCSX4pdCEOSAIECel4Wcv8Lw5MxGGPeaBH2bN99u+LkmJjWxw8sKcXFxLx6xsbFITExEjRo10KdPH7i7u6N///6oU6dOZf0YiIhEFRISgqlTp+KPP/5Aly5dsHDhQvTu3bvMH3JKmsZ9dGgpinMzUev9ztCq0xjFOenIiD6Egge3Yer5FWqYt3vjmOWZxgX+vZLsN998g5MnT8LKygqrVq1C//79yzSGIlLqcwZK6/bt2zh69ChiYmIQFxeHv/76C8XFxZBKpWj2QRvk9f8aKMUvZ9qJtciKO/Hv/9HQhHGPidBv3++d2wmCgDsrhkIozEPt2rVfFJJu3bqhV69evMcAEaksuVyOI0eO4JtvvsHFixfRoUMH+Pr6ws3NDWZm7z6R79bDTPRadeGtzxdnP4NGLaNX37MgF/9smADtOk1Qb9g3JY5/ZpYzrExLXnb44MEDHD58GNu3b0d4eDhat26NTz/9FEOHDq3QtWMUiVqUgf/Ky8vDtWvXEBcXh5Arf+NCTadSbVeUnorinHTIc54h53Y0suJOwqjbaBjaebxz28zfPoVm5gNoampCQ0MDUqn0tcfbvv6uR3m3E+M9lXU7iUTCKRuiChAEAWfOnMHKlStx5swZFBYWon379nB3d4e7uztat279xn9jiw5fQ0BUcplWDgDAg12fojgrDQ0nbnjrazSkEoy0a4JFA1u99tyNGzcQGBiIoKAgyGQySKVSODs7Y+bMmXB1dYVUqtCz52WmlmXgZbEpTzFofUS5tk07+TOy4k+h0dTt0KhpWOJr3WrcRB1JFuRy+RsfxcXFb32upIeybKcKdyOTSCQsZmqw3fMHVZ309HQcP34cgYGBOHbsGDIzM2FqavpiCvf5o3nz5uixMrRU07gvk+dl4+76sdBt0hamHgtLfG0T45r41a3hK9O48fHxuHfvHmrUqIG+ffu+mMatyPJ1RVfa/bdqHAd5A23N8v+j16nfHFmxx1H07ME7y8DE8WPRqkHJr1FlgiBAEASlKC7Ktl1hYaHCZS3FZweFpwzFRVm2e9O2tWvXxrhx4+Dr64urV6/i+vXr+PvvvxEQEIBly5YBAHT1DFFvyo5STeO+7Mnp9RAK82Do6PXO1yalZeP91u0gFOahYcOGsLa2xpgxY2BnZ4cePXpwGvc/VLYMmJvUggR47epVpZGXfAWQSKFpVPJ1DST/ex919vwwOz9xqYf/lj9lK1jVtV15tn1e/hTpe6wqxTWNy1wEnl0IQPa186jdyw8671m98/USiQSb9x3GgM42PIG7FFS2DNTS0YSZcc0SD0OlHV8DqU5NaNdvDo1aRijOyUDOjTDk/BkKAzuPdx4VMDOpiVo6KvsjJHoNy596KemoX0XKx5+PcjHv9KNS53gWtgvpEXth5DwKBh0GlHo76/YdUadO5d7nRlWp9J7MpYVpiSeo6DR8H1lXziDr6lnI87Mh1dKFlmlTmLjOKfFyxMC/J6i4NOd9uYlIdUkkEmhoaEBDQ6NSx5XeSwdKWQaehe1CetguGHYeDkNHzzK9T0Wmi9WNSpeBEXZm2BqZ9Nbn9dr2gl7bXuUau1guwMe+8q+DTUSk6ko7jfssfPe/RcDRC0adh5fpPTiNWzYqXZua1dNHF6s65boKYUk0pBJ0sarzzjWsRET0uufTuCXJiDqI9NCd0LXogBqWtsj/569XHu/CadyyUfmf1LeD2qDnypAyr2UtiaZUgm8Htam08YiI1M27pnFzbkcDAPISLuFBwqXXnm+y4Ohbx+Y0btmpfBlobFwTiwe2woKDVyttzK8Gtir37YuJiOjd07jvjfiu3GNzGrfsVHqa4DlvWzPM7d28Usb6uHcLeFXBPbOJiNQJp3EVi1qUAQCY6tIM33m0gY6mtMy/fBpSCXQ0pfjeow2muLx7fSsREb3bt4PaQLOSywCncctHbcoA8O8RgjOzusLR4t9LT76rFDx/3tHCBGdmdeURASKiSvR8GrcycRq3fFT+nIH/amxcEwHj7HDrYSZ2RqUg+GYqUtJyXlniIsG/Z6K6NDeFj70ZDzcREVURb1szPM7Kx/JTNys8Fqdxy09lb1RUFtn5RUhKy0ZBkRzamlKYm9TikhQiomq052IKvjx8DUVyoUyrvzSkEmhKJfhqYCsWgTdQ+xsVlUUtHU21vtkQEZHYvG3N4GRZB58euorQ24+hIZWUWAqeP+9oYYJvB7Xh1EAFsQwQEZFC4DSueDhNQERECovTuBXDaQIiIlJ6nMatHmq1tJCIiIhexzJARESk5lgGiIiI1BzLABERkZpjGSAiIlJzLANERERqjmWAiIhIzbEMEBERqTmWASIiIjXHMkBERKTmWAaIiIjUHMsAERGRmmMZICIiUnMsA0RERGqOZYCIiEjNaZbmRYIgAAAyMjKqNAwRERFVnuf77ef78bcpVRnIzMwEADRu3LiCsYiIiKi6ZWZmwtDQ8K3PS4R31QUAcrkc9+7dg76+PiQSSaUGJCIioqohCAIyMzPRoEEDSKVvPzOgVGWAiIiIVBdPICQiIlJzLANERERqjmWAiIhIzbEMEBERqTmWASIiIjXHMkBERKTmWAaIiIjU3P8BEf7aPBZUedoAAAAASUVORK5CYII=",
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"nx.draw_networkx(G)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "668b188d",
"metadata": {},
"outputs": [],
"source": [
"assert output_coloring.size == sparsity.shape[0]\n",
"\n",
"projection_matrix = (\n",
" np.arange(n_colors)[:, None] == output_coloring[None, :]\n",
").astype(float)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "4d2355a9",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[1., 0., 0., 0., 0.],\n",
" [0., 1., 1., 0., 0.],\n",
" [0., 0., 0., 1., 1.]])"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"projection_matrix"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "94ccad32",
"metadata": {},
"outputs": [],
"source": [
"projected_eqs = projection_matrix @ pt.stack(equations)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "fb144e67",
"metadata": {},
"outputs": [],
"source": [
"def coo_to_csc(rows, cols, data, shape):\n",
" rows = pt.as_tensor_variable(rows)\n",
" cols = pt.as_tensor_variable(cols)\n",
" n_rows, n_cols = shape\n",
" \n",
" order = pt.argsort(cols)\n",
" csc_data = data[order]\n",
" csc_indices = rows[order]\n",
" sorted_cols = cols[order]\n",
" \n",
" counts = pt.bincount(sorted_cols, minlength=n_cols)\n",
" csc_indptr = pt.concatenate(([0], pt.cumsum(counts).astype(int)))\n",
" \n",
" return csc_data, csc_indices, csc_indptr, shape "
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "3a4dc7a3",
"metadata": {},
"outputs": [],
"source": [
"# TODO: Is this really the fastest way to do this?\n",
"# TODO: Try with the input mapping (I guess this uses R_op instead? Or L_op?)\n",
"compressed_jac = pt.stack(pt.jacobian(projected_eqs, variables), axis=-1)\n",
"\n",
"\n",
"# TODO: Make COO constructor in pytensor directly?\n",
"rows, cols = sparsity.coords\n",
"compressed_index = (output_coloring[rows], cols)\n",
"data = compressed_jac[compressed_index]\n",
"\n",
"\n",
"# Alternative scheme if we had done an input_mapping\n",
"# row, col = sparsity.coords\n",
"# compressed_index = (row, coloring[col])\n",
"# data = compressed_jac_val[compressed_index]\n",
"from pytensor import sparse as pts\n",
"sparse_jac = pts.CSC(*coo_to_csc(rows, cols, data, (n_eqs, n_vars)))"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "c7404821",
"metadata": {},
"outputs": [],
"source": [
"f_sparse_jac = pytensor.function(variables, sparse_jac)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "5f827e79",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"matrix([[2. , 0.54030231, 0. , 0. , 0. ],\n",
" [0. , 3. , 0. , 0. , 0. ],\n",
" [1. , 0. , 2.71828183, 0. , 0. ],\n",
" [3. , 0. , 0. , 2.71828183, 0. ],\n",
" [0. , 1. , 0. , 0. , 4. ]])"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f_sparse_jac(*np.ones(5,)).todense()"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "6b3ab6e6",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"40.9 μs ± 211 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)\n"
]
}
],
"source": [
"%timeit f_sparse_jac(*np.ones(5,))"
]
},
{
"cell_type": "markdown",
"id": "fda1891a",
"metadata": {},
"source": [
"## Sanity Check"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "238e8bc4",
"metadata": {},
"outputs": [],
"source": [
"dense_jac = pt.stack(pt.jacobian(pt.stack(equations), variables), axis=-1)\n",
"dense_fn = pytensor.function(variables, dense_jac)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "60556abe",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[2. , 0.54030231, 0. , 0. , 0. ],\n",
" [0. , 3. , 0. , 0. , 0. ],\n",
" [1. , 0. , 2.71828183, 0. , 0. ],\n",
" [3. , 0. , 0. , 2.71828183, 0. ],\n",
" [0. , 1. , 0. , 0. , 4. ]])"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dense_fn(*np.ones(5,))"
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "f78bdac9",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"22.1 μs ± 140 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)\n"
]
}
],
"source": [
"%timeit dense_fn(1., 1., 1., 1., 1.)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.12.9"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment