diff --git a/notebook/diagonal.ipynb b/notebook/diagonal.ipynb new file mode 100644 index 0000000..842674c --- /dev/null +++ b/notebook/diagonal.ipynb @@ -0,0 +1,110 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "from jacobi import propagate, jacobi\n", + "import numpy as np\n", + "from matplotlib import pyplot as plt\n", + "import timeit" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "def fn(x):\n", + " return np.log(3 * x + 1)\n", + "\n", + "n = np.geomspace(1, 1e2, 10).astype(int)\n", + "\n", + "t = {False: [], True: [], \"noopt\": []}\n", + "for ni in n:\n", + " x = np.linspace(1, 10, ni)\n", + " xcov = np.eye(ni)\n", + " for key in t:\n", + " if key in (False, True):\n", + " timer = timeit.Timer(f\"propagate(fn, x, xcov, diagonal={key})\",\n", + " globals=locals())\n", + " k, ti = timer.autorange()\n", + " t[key].append(ti / k)\n", + " else:\n", + " import jacobi._propagate\n", + "\n", + " orig = jacobi._propagate._try_reduce_jacobian\n", + " jacobi._propagate._try_reduce_jacobian = lambda x: x\n", + "\n", + " timer = timeit.Timer(f\"propagate(fn, x, xcov, diagonal=False)\",\n", + " globals=locals())\n", + " k, ti = timer.autorange()\n", + " t[key].append(ti / k)\n", + "\n", + " jacobi._propagate._try_reduce_jacobian = orig\n" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAG1CAYAAAAV2Js8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/av/WaAAAACXBIWXMAAA9hAAAPYQGoP6dpAABzgElEQVR4nO3dd3yNd//H8dfJ3okgEgRBxBa7RmuUGi0td0s3pVTtqpZuOmhr1EqrQ2l7d/t129QMtam9R2wRsuc51++PtLmFICJxJTnv5+NxHnKu+TmHk/N2fcdlMQzDQERERMQOOZhdgIiIiIhZFIRERETEbikIiYiIiN1SEBIRERG7pSAkIiIidktBSEREROyWgpCIiIjYLQUhERERsVtOZhdQ2NlsNk6dOoW3tzcWi8XsckRERCQXDMMgPj6esmXL4uBw7es+CkI3cOrUKYKDg80uQ0RERPIgKiqK8uXLX3O9gtANeHt7A5lvpI+Pj8nViIiISG7ExcURHByc9T1+LQpCN/Bvc5iPj4+CkIiISBFzo24t6iwtIiIidktBSEREROyWmsbygc1mIy0tzewyRG4rZ2dnHB0dzS5DROSWKAjdorS0NI4cOYLNZjO7FJHbzs/Pj8DAQE0tISJFloLQLTAMg9OnT+Po6EhwcPB15ykQKU4MwyApKYlz584BEBQUZHJFIiJ5oyB0CzIyMkhKSqJs2bJ4eHiYXY7IbeXu7g7AuXPnCAgIUDOZiBRJuoRxC6xWKwAuLi4mVyJijn//A5Cenm5yJSIieaMglA/UP0Lslf7ti0hRpyB0DREREdSsWZPGjRubXYqIiIgUEAWhaxg0aBC7d+9m48aNZpciIiIiBURBSLK0bt2a4cOHA1CpUiWmTJliaj151bt3bx544AFTzn35eygiIoWfRo1JjjZu3Iinp6fZZZgip34vLVq0YM2aNSZUIyJSfGXYMnByMDeKKAhJjkqXLm12CaaaPXs2HTt2zHqukYEiIvnv+RXP4+3izfCGwynlXsqUGtQ0lo8MwyApLcOUh2EYN1VrYmIiTz75JF5eXgQFBTFp0qRs669sGps8eTJ16tTB09OT4OBgBg4cSEJCQrZ9Pv30U4KDg/Hw8KBbt25MnjwZPz+/bNt89NFHVKlSBRcXF8LCwvjqq6+yrbdYLHz22Wd069YNDw8PQkND+e2337LWW61W+vbtS0hICO7u7oSFhTF16tSbeu258e+Myf8+/P39uXDhAo888gjlypXDw8ODOnXq8O233173OB9++CGhoaG4ublRpkwZHnzwwax1NpuN8ePHZ72WevXqMXfu3Hx/LSIihdGak2v4M+pP5h2eR2xqrGl16IpQPkpOt1Lz9UWmnHv3mx3wcMn9X+cLL7zAypUr+fXXXwkICODll19my5YthIeH57i9g4MD06ZNIyQkhMOHDzNw4EBefPFFPvzwQwAiIyMZMGAA7733Hl27dmXp0qW89tpr2Y7x888/M2zYMKZMmUK7du34448/eOqppyhfvjxt2rTJ2m7s2LG8//77TJgwgenTp/PYY49x7Ngx/P39sdlslC9fnh9//JGSJUuydu1a+vfvT1BQED169Mix9nHjxjFu3Ljrv3+7d1OhQoXrbpOSkkLDhg0ZNWoUPj4+zJs3jyeeeIIqVarQpEmTq7bftGkTQ4cO5auvvqJ58+bExMSwevXqrPXjx4/nv//9LzNnziQ0NJRVq1bx+OOPU7p0aVq1anXdWkREirJ0azrvbXgPgEdqPEIVvyqm1WIxbvZSgp2Ji4vD19eX2NhYfHx8sq1LSUnhyJEjhISE4ObmRlJaRpEIQgkJCZQsWZL//ve/PPTQQwDExMRQvnx5+vfvz5QpU6hUqRLDhw+/ZsffuXPnMmDAAKKjowF4+OGHSUhI4I8//sja5vHHH+ePP/7g0qVLQGY/m1q1avHJJ59kbdOjRw8SExOZN28ekHlF6NVXX+Wtt94CMq9ceXl5sWDBgmxNVZcbPHgwZ86cybqa0rt3by5dusQvv/yS9dpiYmKu+55UqlQJJyenrBrc3NyyzZT83//+N8cO2Pfddx/Vq1dn4sSJQGZn6fDwcKZMmcJPP/3EU089xYkTJ/D29s62X2pqKv7+/ixdupRmzZplLX/66adJSkrim2++uW69hcWVnwERkdz4fOfnfLD5A0q6leT3br/j7eJ9451u0vW+vy+nK0L5yN3Zkd1vdjDt3Ll16NAh0tLSaNq0adYyf39/wsLCrrnP0qVLGT9+PHv37iUuLo6MjAxSUlJISkrCw8ODffv20a1bt2z7NGnSJFsw2rNnD/3798+2TYsWLa5q2qpbt27Wz56envj4+GTd0woy53j6/PPPOX78OMnJyaSlpV3zSta/r83f3/+a63PywQcf0K5du6znQUFBWK1Wxo0bxw8//MDJkydJS0sjNTX1mrdXad++PRUrVqRy5cp07NiRjh07ZjX5HTx4kKSkJNq3b59tn7S0NOrXr39TtYqIFCXnks7x8faPARjecHiBhKCboSCUjywWy001TxUVR48e5b777uPZZ5/lnXfewd/fnzVr1tC3b1/S0tLy/T5rzs7O2Z5bLBZsNhsA3333HSNHjmTSpEk0a9YMb29vJkyYwPr16695vLw0jQUGBlK1atVs27z77rtMnTqVKVOmZPWXGj58OGlpaTke09vbmy1btrBixQoWL17M66+/zpgxY9i4cWNW/6p58+ZRrly5bPu5urpet1YRkaJs8ubJJGUkUbd0XbpW6Wp2OQpC9qhKlSo4Ozuzfv36rC//ixcvsn///hz7pmzevBmbzcakSZNwcMjsX//DDz9k2yYsLOyqySevfF6jRg0iIyPp1atX1rLIyEhq1qyZ69ojIyNp3rw5AwcOzFp26NCh6+4zYMCAa/Yf+lfZsmVzde7777+fxx9/HMjs7Lx///7r1u/k5ES7du1o164db7zxBn5+fvz555+0b98eV1dXjh8/rv5AImI3Np/dzLzD87Bg4eWmL+NgMX/MloKQHfLy8qJv37688MILlCxZkoCAAF555ZWskHOlqlWrkp6ezvTp0+nSpQuRkZHMnDkz2zZDhgzhrrvuYvLkyXTp0oU///yTBQsWZJuT54UXXqBHjx7Ur1+fdu3a8fvvv/PTTz+xdOnSXNceGhrKl19+yaJFiwgJCeGrr75i48aNhISEXHOfvDSNXevcc+fOZe3atZQoUYLJkydz9uzZawahP/74g8OHD3PXXXdRokQJ5s+fj81mIywsDG9vb0aOHMlzzz2HzWajZcuWxMbGEhkZiY+PT7awKCJSHFhtVsavHw9A99Du1CpZy+SKMpkfxcQUEyZM4M4776RLly60a9eOli1b0rBhwxy3rVevHpMnT+a9996jdu3afP3114wfPz7bNi1atGDmzJlMnjyZevXqsXDhQp577rlsHWgfeOABpk6dysSJE6lVqxYff/wxs2fPpnXr1rmu+5lnnqF79+707NmTpk2bcuHChWxXhwrSq6++SoMGDejQoQOtW7cmMDDwujNY+/n58dNPP9G2bVtq1KjBzJkz+fbbb6lVK/PD/9Zbb/Haa68xfvx4atSoQceOHZk3b951Q52ISFH14/4f2XdxH94u3gxtMNTscrJo1NgN3MyoMcmuX79+7N27N9uQcSle9BkQkdy4mHKR+36+j7i0OF5q8hKP1ni0wM+pUWNy202cOJH27dvj6enJggUL+OKLL7LmGRIREfs1fet04tLiCC0RSo+w6/fZvN0UhCTfbNiwgffff5/4+HgqV67MtGnTePrpp80uS0RETLT7wm7m7s+c5+3lJi+bfm+xKxWuaqRIu3IkmYiI2DebYWPc+nEYGHQK6USjwEZml3QVdZYWERGRAvHH4T/Yfn477k7uPN/webPLyZGCkIiIiOS7hLQEJm+aDED/uv0p41nG5IpypiAkIiIi+W7m9plcSLlARZ+KPFnzSbPLuSYFIREREclXhy8d5us9XwMwqvEoXBxdTK7o2hSEREREJN8YhsH4DePJMDJoXb41d5a/0+ySrktBSLK0bt2a4cOHA1CpUiWmTJliaj151bt37+vO+CwiIgVn2fFl/HX6L1wcXHix8Ytml3NDCkKSo40bN9K/f3+zy7it5syZg8Viue7j6NGjZpcpIlJoJWckM2HjBAB61epFsE+wyRXdmIKQ5Kh06dJ4eHiYXcZt1bNnT06fPp31aNasGf369cu2LDj4fx/qtLQ0E6sVESl8Zu+czanEUwR6BvJ0naIxoa6CkJ1KTEzkySefxMvLi6CgICZNmpRt/ZVNY5MnT6ZOnTp4enoSHBzMwIEDSUhIyLbPp59+SnBwMB4eHnTr1o3Jkyfj5+eXbZuPPvqIKlWq4OLiQlhYGF999VW29RaLhc8++4xu3brh4eFBaGgov/32W9Z6q9VK3759CQkJwd3dnbCwMKZOnZov74m7uzuBgYFZDxcXFzw8PLKejx49mv/85z+88847lC1blrCwsKyaf/nll2zH8vPzY86cOVnPo6Ki6NGjB35+fvj7+3P//ffr6pKIFCsnE07y+c7PARjZaCQezkXjP9MKQvnJMCAt0ZzHTd4794UXXmDlypX8+uuvLF68mBUrVrBly5Zrbu/g4MC0adPYtWsXX3zxBX/++Scvvvi/tt/IyEgGDBjAsGHD2LZtG+3bt+edd97Jdoyff/6ZYcOG8fzzz7Nz506eeeYZnnrqKZYvX55tu7Fjx9KjRw/+/vtvOnfuzGOPPUZMTAwANpuN8uXL8+OPP7J7925ef/11Xn755evOaj1u3Di8vLyu+zh+/Hiu3rdly5axb98+lixZwh9//JGrfdLT0+nQoQPe3t6sXr2ayMhIvLy86Nixo64qiUixMWHjBFKtqTQJbMI9Fe8xu5xc0y028lN6Eowra865Xz4FLp652jQhIYFZs2bx3//+l7vvvhuAL774gvLly19zn387UUPm1aK3336bAQMGZN1Udfr06XTq1ImRI0cCUK1aNdauXZstLEycOJHevXszcOBAAEaMGMFff/3FxIkTadOmTdZ2vXv35pFHHgEyQ8y0adPYsGEDHTt2xNnZmbFjx2ZtGxISwrp16/jhhx/o0SPnG/kNGDDgmuv+VbZs7v7ePD09+eyzz3Bxyf1Q0O+//x6bzcZnn32GxWIBYPbs2fj5+bFixQruuafo/MIQEcnJ2pNrWXZ8GY4WR15q8lLW77qiQEHIDh06dIi0tDSaNm2atczf3z+rqScnS5cuZfz48ezdu5e4uDgyMjJISUkhKSkJDw8P9u3bR7du3bLt06RJk2xBaM+ePVd1wG7RosVVTVt169bN+tnT0xMfHx/OnTuXtSwiIoLPP/+c48ePk5ycTFpaGuHh4des3d/fH39//2uuvxl16tS5qRAEsH37dg4ePIi3t3e25SkpKRw6dChf6hIRMUu6NZ3xG8YD8Ej1R6haoqrJFd0cBaH85OyReWXGrHMXkKNHj3Lffffx7LPP8s477+Dv78+aNWvo27cvaWlp+d6p2tnZOdtzi8WCzWYD4LvvvmPkyJFMmjSJZs2a4e3tzYQJE1i/fv01jzdu3DjGjRt33XPu3r2bChUq3LA2T8+rr7pZLBaMK5om09PTs35OSEigYcOGfP3111ftW7p06RueU0SkMPt6z9ccjTuKv5s/z4Y/a3Y5N80uglC3bt1YsWIFd999N3Pnzi24E1ksuW6eMlOVKlVwdnZm/fr1WV/+Fy9eZP/+/bRq1eqq7Tdv3ozNZmPSpEk4OGR2K7uyT05YWBgbN27MtuzK5zVq1CAyMpJevXplLYuMjKRmzZq5rj0yMpLmzZtnNa8BN7yqkp9NYzkpXbo0p0+fznp+4MABkpKSsp43aNCA77//noCAAHx8fPJ8HhGRwuZ80nk+2v4RAMMbDMfHpej9jrOLIDRs2DD69OnDF198YXYphYKXlxd9+/blhRdeoGTJkgQEBPDKK69khZwrVa1alfT0dKZPn06XLl2IjIxk5syZ2bYZMmQId911F5MnT6ZLly78+eefLFiwIFs78QsvvECPHj2oX78+7dq14/fff+enn35i6dKlua49NDSUL7/8kkWLFhESEsJXX33Fxo0bCQkJueY++dk0lpO2bdsyY8YMmjVrhtVqZdSoUdmuaj322GNMmDCB+++/nzfffJPy5ctz7NgxfvrpJ1588cXr9s0SESnMPtj8AUkZSdQtVZf7q95vdjl5Yhejxlq3bn1V/wx7N2HCBO688066dOlCu3btaNmyJQ0bNsxx23r16jF58mTee+89ateuzddff8348eOzbdOiRQtmzpzJ5MmTqVevHgsXLuS5557Dzc0ta5sHHniAqVOnMnHiRGrVqsXHH3/M7Nmzad26da7rfuaZZ+jevTs9e/akadOmXLhwIdvVITNMmjSJ4OBg7rzzTh599FFGjhyZrbnQw8ODVatWUaFCBbp3706NGjXo27cvKSkpukIkIkXW1nNb+f3w71iw8FLTl3CwFM1IYTGu7Nxwm61atYoJEyawefNmTp8+zc8//3zV7REiIiKYMGECZ86coV69ekyfPp0mTZrc1HlWrFjBjBkzbrppLC4uDl9fX2JjY6/60kpJSeHIkSOEhIRk+8KXTP369WPv3r2sXr3a7FKkgOgzIGKfrDYrD897mL0xe+ke2p2xzcfeeKfb7Hrf35czvWksMTGRevXq0adPH7p3737V+u+//54RI0Ywc+ZMmjZtypQpU+jQoQP79u0jICAAgPDwcDIyMq7ad/HixbfU90NuzsSJE2nfvj2enp4sWLCAL774Imt4vYiIFB//d+D/2BuzF29nb4bWH2p2ObfE9CDUqVMnOnXqdM31kydPpl+/fjz11FMAzJw5k3nz5vH5558zevRoALZt25Zv9aSmppKampr1PC4uLt+OXdxt2LCB999/n/j4eCpXrsy0adN4+umiMcW6iIjkzqWUS0zbOg2AQfUHUdK9pMkV3RrTg9D1pKWlsXnzZl566aWsZQ4ODrRr145169YVyDnHjx+fbcI+yb3rze4sIiLFw4xtM4hNjaWqX1V6hvU0u5xbVqh7NkVHR2O1WilTpky25WXKlOHMmTO5Pk67du146KGHmD9/PuXLl79uiHrppZeIjY3NekRFReW5fhERkeJkb8xeftz/IwAvN30ZJ4dCfT0lV4r+K8iFmxme7erqiqurawFWIyIiUvQYhsG49eOwGTY6VupI48DGZpeULwr1FaFSpUrh6OjI2bNnsy0/e/YsgYGBJlUlIiJif/44/Adbz23F3cmd5xs9b3Y5+aZQByEXFxcaNmzIsmXLspbZbDaWLVtGs2bNTKxMRETEfiSmJ/LB5g8A6FenH4GexedihOlNYwkJCRw8eDDr+ZEjR9i2bRv+/v5UqFCBESNG0KtXLxo1akSTJk2YMmUKiYmJWaPIREREpGB9vP1jziefp4J3BXrV6nXjHYoQ04PQpk2baNOmTdbzESNGANCrVy/mzJlDz549OX/+PK+//jpnzpwhPDychQsXXtWBOr9FREQQERGB1Wot0POIiIgUZkdij/DVnq8AGNVkFC6OLiZXlL9Mbxpr3bo1hmFc9ZgzZ07WNoMHD+bYsWOkpqayfv16mjZtWuB1DRo0iN27d19141B7M2fOHPz8/G64ncVi4Zdffinwesy2YsUKLBYLly5duqXjVKpUiSlTpuRLTdcyZswYwsPDC/QcIlK8GYbBexveI8OWwV3l7+Ku8neZXVK+Mz0ISeHWs2dP9u/fn/Xc7C/X2xm4WrduzfDhw7Mta968OadPn8bX1/eWjr1x40b69+9/S8e4XE7vy8iRI7P1rxMRuVnLo5YTeSoSZwdnRjUeZXY5BcL0pjEp3Nzd3XF3dze7jELDxcUlX0Ysli5dOh+quT4vLy+8vLwK/DwiUjylZKTw/sb3AehVqxcVfCqYXFHB0BUhO/PHH3/g5+eX1fdp27ZtWCyWrNuVADz99NM8/vjjQPamsTlz5jB27Fi2b9+OxWLBYrFka8KMjo6mW7dueHh4EBoaym+//Zbt3CtXrqRJkya4uroSFBTE6NGjs90jLqfmovDwcMaMGZO1HqBbt25YLJas5znZsWMHbdu2xd3dnZIlS9K/f38SEhKy1vfu3ZsHHniAsWPHUrp0aXx8fBgwYABpaWlZ61euXMnUqVOzXuvRo0evahr79/35448/CAsLw8PDgwcffJCkpCS++OILKlWqRIkSJRg6dGi2/maXv9Y5c+ZknePyx7+ve+PGjbRv355SpUrh6+tLq1at2LJlS7Zj5fS+XHn1zmaz8eabb1K+fHlcXV2z+tv96+jRo1gsFn766SfatGmDh4cH9erVK7BZ3EWkcJu9azYnE05SxqMM/er0M7ucAqMglI8MwyApPcmUh2EYuarxzjvvJD4+nq1btwKZ4aRUqVKsWLEia5uVK1fSunXrq/bt2bMnzz//PLVq1eL06dOcPn2anj3/N7362LFj6dGjB3///TedO3fmscceIyYmBoCTJ0/SuXNnGjduzPbt2/noo4+YNWsWb7/9dq7f33/7a82ePZvTp09fs/9WYmIiHTp0oESJEmzcuJEff/yRpUuXMnjw4GzbLVu2jD179rBixQq+/fZbfvrpp6zbq0ydOpVmzZrRr1+/rNcaHByc4/mSkpKYNm0a3333HQsXLmTFihV069aN+fPnM3/+fL766is+/vhj5s6dm+P+PXv2zDrH6dOn+fbbb3FycqJFixYAxMfH06tXL9asWcNff/1FaGgonTt3Jj4+/qbel6lTpzJp0iQmTpzI33//TYcOHejatSsHDhzItt0rr7zCyJEj2bZtG9WqVeORRx7J8abGIlJ8nUo4xawdswAY2WgkHs4eJldUcNQ0dg15GTWWnJFM028KviN3TtY/uj5X/1B9fX0JDw9nxYoVNGrUiBUrVvDcc88xduxYEhISiI2N5eDBg7Rq1eqqfd3d3fHy8sLJySnH5qHevXvzyCOPADBu3DimTZvGhg0b6NixIx9++CHBwcHMmDEDi8VC9erVOXXqFKNGjeL111/HweHGmfzf5iQ/P7/rNk998803pKSk8OWXX+Lp6QnAjBkz6NKlC++9917WiEMXFxc+//xzPDw8qFWrFm+++SYvvPACb731Fr6+vri4uODh4XHDprD09HQ++ugjqlSpAsCDDz7IV199xdmzZ/Hy8qJmzZq0adOG5cuXZwuOl7+v/zY/Hjp0iEGDBjFu3Djat28PQNu2bbNt/8knn+Dn58fKlSu57777cv2+TJw4kVGjRvHwww8D8N5777F8+XKmTJlCRERE1nYjR47k3nvvBTLDba1atTh48CDVq1e/7vsgIsXHxE0TSbWm0jiwMR0qdTC7nAKlK0LXUJxHjbVq1YoVK1ZgGAarV6+me/fu1KhRgzVr1rBy5UrKli1LaGjoTR+3bt26WT97enri4+PDuXPnANizZw/NmjXDYrFkbdOiRQsSEhI4ceLErb+oy+zZs4d69eplhaB/z2Wz2di3b1/Wsnr16uHh8b/w2KxZMxISEm76/nIeHh5ZIQgy74VXqVKlbP1zypQpk/VeXEtsbCz33Xcf9957Ly+88ELW8rNnz9KvXz9CQ0Px9fXFx8eHhIQEjh8/nusa4+LiOHXqVNZVpn+1aNGCPXv2ZFt2+d9jUFAQwA1rF5HiY92pdSw5tgRHiyOjm4zO9nu7ONIVoXzk7uTO+kfXm3bu3GrdujWff/4527dvx9nZmerVq9O6dWtWrFjBxYsXc7walBvOzs7ZnlssFmw2W673d3BwuKqJLz09PU+13E45ve6bfS+sVis9e/bEx8eHTz75JNu6Xr16ceHCBaZOnUrFihVxdXWlWbNmWf2Z8tvltf/7C/Bm/h5FpOhKt6Xz7oZ3AegZ1pNqJaqZXFHBUxDKRxaLpUi0o/7bT+iDDz7ICj2tW7fm3Xff5eLFizz//LXvIePi4pKnSSZr1KjB//3f/2EYRtaXa2RkJN7e3pQvXx7IbPo6ffp01j5xcXEcOXIk23GcnZ1veP4aNWowZ84cEhMTs64KRUZG4uDgQFhYWNZ227dvJzk5OatZ6q+//sLLyyurL1BeX2tePPfcc+zYsYNNmzbh5uaWbV1kZCQffvghnTt3BiAqKoro6Ohs29zoffHx8aFs2bJERkZmC7qRkZE0adIkH1+JiBRl3+z5hsOxh/F382dQ/UFml3NbqGnMDpUoUYK6devy9ddfZ3WKvuuuu9iyZQv79++/7hWhSpUqZd0GJTo6mtTU1Fydc+DAgURFRTFkyBD27t3Lr7/+yhtvvMGIESOy+ge1bduWr776itWrV7Njxw569eqFo6PjVedftmwZZ86c4eLFizme67HHHsPNzY1evXqxc+dOli9fzpAhQ3jiiSeyzUielpZG37592b17N/Pnz+eNN95g8ODBWfVUqlSJ9evXc/ToUaKjowvsqsjs2bP58MMPmTlzJhaLhTNnznDmzJmsUW6hoaF89dVX7Nmzh/Xr1/PYY49dNaVBbt6XF154gffee4/vv/+effv2MXr0aLZt28awYcMK5HWJSNESnRzNR9s/AmBYg2H4uPiYXNHtoSBkp1q1aoXVas0KQv7+/tSsWZPAwMBsV02u9J///IeOHTvSpk0bSpcuzbfffpur85UrV4758+ezYcMG6tWrx4ABA+jbty+vvvpq1jYvvfQSrVq1yuon88ADD2TrewMwadIklixZQnBwMPXr18/xXB4eHixatIiYmBgaN27Mgw8+yN13382MGTOybXf33XcTGhrKXXfdRc+ePenatWvWkHXI7DTs6OhIzZo1KV269E31ybkZK1euxGq10rVrV4KCgrIeEydOBGDWrFlcvHiRBg0a8MQTTzB06FACAgKyHSM378vQoUMZMWIEzz//PHXq1GHhwoX89ttveeoPJiLFzwebPyAxPZHaJWvzQNUHzC7ntrEYuR13bafi4uLw9fUlNjYWH5/s6TglJYUjR44QEhJyVXOGFG69e/fm0qVLdnFbkIKkz4BI8bDt3DaeWPAEAF93/pq6peveYI/C73rf35fTFaFriIiIoGbNmjRu3NjsUkRERAqM1WZl/IbxADxQ9YFiEYJuhoLQNRTn4fMiIiL/+ungT+y+sBsvZy+GNbC/PoMaNSZ26fJbg4iI2KvY1FimbZkGwKDwQZRyL2VyRbefrgiJiIjYqRlbZ3Ap9RJV/arSs/rVM9/bAwUhERERO7QvZh8/7P8BgJeavISzg/MN9iieFITygQbeib3SjNMiRZNhGIxbPw6bYeOeivfQJMh+J1ZVH6Fb4OzsjMVi4fz585QuXbrY349F5F+GYZCWlsb58+dxcHDAxcXF7JJE5CYsOLKALee24OboxshGI80ux1QKQrfA0dGR8uXLc+LECY4ePWp2OSK3nYeHBxUqVMiajVtECr+k9CQmbZoEQL+6/QjyCjK5InMpCF1DREQEERERN7zXlJeXF6GhoUXi5qAi+cnR0REnJyddCRUpYj7++2POJZ+jvFd5etXqZXY5ptPM0jeQ25kpRURECrujsUfp9ls3MmwZTG87ndbBrc0uqcBoZmkRERHJYhgG7218jwxbBi3LtaRV+WvfYNueKAiJiIjYgZUnVrLm5BqcHJwY1XiUmrX/oSAkIiJSzKVaU3lvw3sA9KrZi0q+lcwtqBBREBIRESnm5uycw4mEEwR4BNC/bn+zyylUFIRERESKsdMJp/lsx2cAPN/weTycPUyuqHBREBIRESnGJm6aSIo1hYZlGtIppJPZ5RQ6CkIiIiLF1PrT61l8bDEOFgdeavKSOkjnQEFIRESkGEq3pfPuhncB6BnWkzD/MJMrKpwUhK4hIiKCmjVr0rhxY7NLERERuWnf7f2Og5cOUsK1BIPCB5ldTqGlmaVvQDNLi4hIUROdHE2Xn7uQkJ7AG83e4MFqD5pd0m2nmaVFRETs1NQtU0lIT6BmyZp0q9rN7HIKNQUhERGRYuTv83/zy8FfAHi56cs4OjiaW1AhpyAkIiJSTNgMG+PXjwega5Wu1Ctdz+SKCj8FIRERkWLil4O/sPPCTrycvXiu4XNml1MkKAiJiIgUA7GpsUzZPAWAZ+s9Syn3UuYWVEQoCImIiBQDH277kIupF6niW4VHajxidjlFhoKQiIhIEbf/4n6+3/c9AKObjsbZwdnkiooOBSEREZEizDAMxq8fj9Ww0r5ie+4IusPskooUBSEREZEibNHRRWw6uwk3RzdGNhppdjlFjoLQNegWGyIiUtglpScxcdNEAPrW6UtZr7ImV1T06BYbN6BbbIiISGE1dctUPtvxGeW8yvHrA7/i6uhqdkmFhm6xISIiUowdjzvOF7u+AODFxi8qBOWRgpCIiEgR9N7G90i3pdOibAvaBLcxu5wiS0FIRESkiFl1YhWrTqzCycGJUU1GYbFYzC6pyFIQEhERKUJSram8u+FdAJ6o+QQhviEmV1S0KQiJiIgUIV/u+pKo+ChKu5fmmbrPmF1OkacgJCIiUkScSTzDpzs+BWBEoxF4OnuaXNGtiU9JZ/Oxi6bW4GTq2UVERCRXEtMTGbVqFMkZyTQIaMC9IfeaXVKeGYbBz1tPMm7+XjJsNpY/35oSni6m1KIgJCIiUsglpCUwcNlAtp7bipezF6/e8WqR7SC961Qsb/y6i03/XAkKKeXJmbgUBSERERG5WnxaPAOWDuDv83/j7eLNJ+0/IbREqNll3bTYpHQmLdnHf/86hs0ADxdHhrQNpU/LSrg6OZpWl4KQiIhIIRWXFseAJQPYEb0DHxcfPrnnE2qVrGV2WTfFZjP4cXMU7y3cR0xiGgD31Q3ilXtrEOTrbnJ1CkIiIiKFUmxqLP2X9Gf3hd34ufrx6T2fUt2/utll3ZTtUZd4/bddbI+6BEBogBdju9aiedVS5hZ2GQUhERGRQuZSyiX6L+nPnpg9lHAtwaf3fEqYf5jZZeVaTGIaExbt5buNURgGeLk6MbxdKL2aV8LZsXANWFcQEhERKURiUmLot7gf+y/ux9/Nn8/u+azI9Amy2gy+2XCciYv2EZucDkD3+uUY3ak6AT5uJleXMwUhERGRQuJC8gWeXvw0By8dpKRbSWZ1mEUVvypml5Urm4/F8Novu9h9Og6A6oHevPVAbRpX8je5sutTELqGiIgIIiIisFqtZpciIiJ2IDo5mqcXPc2h2EOUdi/NZx0+o7JvZbPLuqHz8am8u2Av/7flBAA+bk6M7BDGo00q4FTImsFyYjEMwzC7iMIsLi4OX19fYmNj8fHxMbscEREphs4nnafv4r4ciT1CgEcAn3f4nIo+Fc0u67oyrDa+XHeMD5bsJz41A4CejYJ5sWMYJb1cTa4u99/fuiIkIiJiorOJZ3l68dMcjTtKoGcgn9/zOcE+wWaXdV1/Hb7AG7/uYt/ZeADqlvflzftrEx7sZ25heaAgJCIiYpIziWfos6gPUfFRBHkGMavDLIK9C28IOhObwrj5e/ht+ykASng482LH6vRoFIyjQ9Gc6VpBSERExASnE07TZ1EfTiScoJxXOWZ1mEU5r3Jml5WjtAwbsyOPMG3ZARLTrFgs8FjTCoy8Jww/D3NujZFfFIRERERus5MJJ+m7qC8nE05S3qs8n3f4nCCvILPLytHqA+d547ddHD6fCECDCn68eX9tapfzNbmy/KEgJCIichtFxUfRd1FfTieepoJ3BWZ1mEWgZ6DZZV3l5KVk3v5jNwt2ngGglJcLozvVoHv9cjgU0WawnCgIiYiI3CbH447TZ1EfziadpZJPJWZ1mEWAR4DZZWWTkm7ls9WHmbH8ICnpNhwdLDzZrCLPta+Gj5uz2eXlOwUhERGR2+Bo7FH6LurLueRzhPiGMOueWZT2KG12Wdn8ufcsY3/fzbELSQA0CfHnzftrUT2w+E4foyAkIiJSwA7HHqbvor5EJ0dTxbcKn3X4jFLuhefGo8cuJPLWH7tZuuccAGV8XHm5cw261iuLxVJ8msFyoiAkIiJSgA5dOkTfRX25kHKB0BKhfNr+U0q6lzS7LACS06x8tOIgM1cdJi3DhpODhb53hjCkbShervYREezjVYqIiJjgwMUDPL34aWJSYggrEcan93xKCbcSZpeFYRgs2nWWt/7YzclLyQC0rFqKMV1rUTXAy+Tqbi8FIRERkQKwL2Yf/Rb342LqRWr41+DTez7F19X8IeeHzicw5rddrD4QDUA5P3deu68GHWoFFvtmsJwoCImIiOSzPRf20G9JP2JTY6lVshYft//Y9BCUmJrB9D8PMmvNYdKtBi6ODjzTqjIDW1fF3cXR1NrMpCAkIiKSj3Zd2EW/xf2IT4unTqk6zGw/Ex8X80ZdGYbBH3+f5p15ezgTlwJAm7DSvNGlFpVKeZpWV2GhICQiIpJPdpzfwTNLniE+PZ56pevxUbuP8HbxNq2e/WfjeePXXaw7fAGAYH933rivFu1qljGtpsJGQUhERCQfbD+/nQFLBpCQnkCDgAZ82O5DPJ3NueISn5LOlKUHmLP2KFabgauTA4PaVKX/XZVxc7bfZrCcKAiJiIjcoq3ntvLs0mdJTE+kYZmGfHj3h3g4e9z2OgzD4OetJxk3fy/RCakAdKhVhlfvrUmw/+2vpyhQEBIREbkFm85sYuCygSRnJNMksAnT2043JQTtOhXLG7/uYtOxiwCElPJkTNdatKpWuGavLmwUhERERPJo45mNDFo2iOSMZO4IuoNpbafh7uR+W2uITUpn0pJ9/PevY9gM8HBxZEjbUPq0rISrk5rBbkRBSEREJA/+Ov0XQ5YNIcWaQvOyzZnaZipuTm637fw2m8GPm6N4b+E+YhLTALivbhCv3FuDIN/bG8aKMgWha4iIiCAiIgKr1Wp2KSIiUsisPbmWocuHkmpNpWW5lkxpMwVXR9fbdv7tUZd4/bddbI+6BEBogBdj769F8yqF5/5lRYXFMAzD7CIKs7i4OHx9fYmNjcXHp/jefVdERHJnzck1DPtzGGm2NFqVb8Xk1pNxcXS5Lee+kJDKxMX7+G5jFIYBXq5ODG8XSq/mlXB2dLgtNRQVuf3+1hUhERGRXFp1YhXDlw8n3ZZOm+A2TGo1CWdH5wI/b3KalVlrDjNz5WESUjMA6F6/HKM7VyfA+/Y1xxVHCkIiIiK5sPz4ckasHEGGLYN2Fdrxfqv3cXYo2BCUYbUxd/MJPli6n7NxmcPha5fz4Y0utWhcyb9Az20vFIRERERuYNmxZYxcNZIMWwb3VLyHd+96t0BDkGEYLNtzjvcW7uXAuQQAypdw54UOYXSpWxYHB/u7OWpBURASERG5jsVHFzNq1SgyjAw6VerEuDvH4eRQcF+fW49fZPz8vWw4GgOAn4czQ9qG8vgdFTQcvgAoCImIiFzDwiMLGb16NFbDyn2V7+OtFm8VWAg6Ep3IhEV7mb/jDACuTg70aRnCgFZV8HUv+H5I9kpBSEREJAfzDs/j5TUvYzNsdK3SlTebv4mjQ/5fkYlOSGX6sgN8vf44GTYDiwUebFCeEfdU03xAt4GCkIiIyBV+P/Q7r0a+is2w0a1qN8Y0H4ODJX+HpyelZTBr9RFmrjxEYlrmnHVtwkozqlN1qgdqupbbRUFIRETkMr8c/IXXI1/HwOA/of/h9Wav52sIyrDa+GHTCaYs3c+5+MyRYHXL+zK6U3VNiGgCBSEREZF//N/+/2PsurEYGPQM68nLTV/OtxBkGAZLdp/lvYV7OXQ+EYBgf3de7FCde+sEaSSYSRSEREREgB/2/cBbf70FwKPVH2V0k9FYLPkTTrYcv8j4+XvYeDTzzvAlPJwZencojzWtiIuTZoQ2k4KQiIjYvW/3fsu49eMAeLzG47zY+MV8CUGHzycwYdE+Fuz830iwp+8M4ZlWVfBx00iwwkBBSERE7NrXe77m3Q3vAtC7Vm9GNBxxyyHofHwq05Yd4JsNx7HaDBws8GDD8jzXXiPBChsFIRERsVtf7vqSCZsmANCndh+GNxh+SyEoMTWDz1Yf4ZNV/xsJ1rZ6AKM6Vics0Dtfapb8pSAkIiJ2afbO2UzePBmAfnX6MaT+kDyHoAyrje83RfHBkgNEJ2SOBKtX3peXOtfgjsol861myX8KQiIiYnc+2/EZU7dMBeDZes/ybL1n8xSCDMNg8T8jwQ7/MxKsYkkPXugQxr11gvKts7UUHAUhERGxKzO3zyRiWwQAg8IHMaDegDwdZ/OxGMbN38vmY5kjwfw9XRjatiqPaiRYkaIgJCIidsEwDD7a/hEfbf8IgGENhvF0nadv+jiHzifw/sK9LNp1FgA3ZweeblmZZ1pVxlsjwYocBSERESn2DMNgxrYZfPL3JwCMaDiCp2o/dVPHOBefwtSlB/huY1TWSLAejYIZ3q4agb5uBVG23AZ5CkKxsbFYrVb8/f2zLY+JicHJyQkfH90jRURECgfDMJi6ZSqzds4CYGSjkfSq1SvX+yekZvDpqsN8uvowSf+MBGtXI3MkWGgZjQQr6vIUhB5++GG6dOnCwIEDsy3/4Ycf+O2335g/f36+FCciInIrktKTmLBpAnP3zwVgdJPRPFbjsVztm2618d3GKKYu3U90QhoA9YL9eLlTdZpqJFixYTEMw7jZnfz9/YmMjKRGjRrZlu/du5cWLVpw4cKFfCvQbHFxcfj6+hIbG6srXSIiRcjms5t5LfI1ouKjAHi56cs8Uv2RG+5nGAaLdp3h/YX7OBydORKsUkkPXuhQnc51AjUSrIjI7fd3nq4IpaamkpGRcdXy9PR0kpOT83JIERGRfJGckcy0LdP4es/XGBgEegYyttlYmpdrfsN9Nx2NYdz8PWw5fgmAkp4uDL07lEeaVNBIsGIqT0GoSZMmfPLJJ0yfPj3b8pkzZ9KwYcN8KUxERORmbTu3jVcjX+VY3DEAuod2Z2SjkXi7XL8vz8FzCby3cC9LdmeOBHN3duTpO0Pof5dGghV3eQpCb7/9Nu3atWP79u3cfffdACxbtoyNGzeyePHifC1QRETkRlIyUojYFsEXu77AwCDAPYAxzcdwZ/k7r7vfubgUPlh6gB82/W8kWM/GFXiuXSgBPhoJZg/yFIRatGjBunXreP/99/nhhx9wd3enbt26zJo1i9DQ0PyuUURE5Jr+Pv83r6x5haNxRwHoWqUrLzZ+EV9X32vuk5CawScrD/Hp6iMkp2eOBGtfswyjOoZRNUAjwexJnjpL2xN1lhYRKZxSral8uO1D5uyag82wUcq9FGOajaFVcKtr7pNutfHthuNMXXqAC4mZI8HCg/14uXMNmoT4X3M/KXoKtLM0wKFDh5g9ezaHDx9mypQpBAQEsGDBAipUqECtWrXyelgREZEb2hm9k1fXvMqh2EMA3Ff5PkY3GX3Nq0CGYbBg5xkmLNrHkX9GgoWU8uTFDmF0rK2RYPYsT0Fo5cqVdOrUiRYtWrBq1SrefvttAgIC2L59O7NmzWLu3Ln5XaeIiAhp1jRmbp/J5zs/x2pY8Xfz5/Vmr3N3hbuvuc+mozG8M38PW/8ZCVbKy4Vhd4fycJMKODtqJJi9y1MQGj16NG+//TYjRozA2/t/balt27ZlxowZ+VaciIjIv3Zf2M0ra17h4KWDAHSq1ImXmr5ECbcSOW5/OjaZ8fP38tv2U0DmSLB+d1Wm/12V8XLVHaYkU57+JezYsYNvvvnmquUBAQFER0ffclH5KSoqiieeeIJz587h5OTEa6+9xkMPPWR2WSIikkvp1nQ+2fEJn/79adZVoFfveJX2FdvnuH1KupXPVh8mYvkhktOtWCzQs1EwI9pX00gwuUqegpCfnx+nT58mJCQk2/KtW7dSrly5fCksvzg5OTFlyhTCw8M5c+YMDRs2pHPnznh6eppdmoiI3MC+mH28suYV9l3cB8A9Fe/hlTtewd/t6o7NmTNCn+Wd+buJismc3LdhxRKM6VKLOuWvPYJM7Fue7zU2atQofvzxRywWCzabjcjISEaOHMmTTz6Z3zXekqCgIIKCggAIDAykVKlSxMTEKAiJiBRi6bZ0Zu2YxcfbPybDyMDP1Y9X7niFjpU65rj9/rPxjP19F5EHM2/xFOjjxkudq9O1Xll1hJbrylMvsXHjxlG9enWCg4NJSEigZs2a3HXXXTRv3pxXX331po61atUqunTpQtmymf9Yf/nll6u2iYiIoFKlSri5udG0aVM2bNiQl7LZvHkzVquV4ODgPO0vIiIF78DFAzw27zEitkWQYWRwd4W7+fn+n3MMQbFJ6Yz5bRedpq4m8uAFXJwcGNymKsueb8X94eUUguSG8nRFyMXFhU8//ZTXX3+dHTt2kJCQQP369fM0mWJiYiL16tWjT58+dO/e/ar133//PSNGjGDmzJk0bdqUKVOm0KFDB/bt20dAQAAA4eHhOd77bPHixZQtWxaAmJgYnnzyST799NPr1pOamkpqamrW87i4uJt+TSIicvMybBnM3jmbD7d/SIYtAx8XH15p+gqdQjpdFWisNoPvNh5n4qJ9XExKB+CemmV49d6aVCjpYUb5UkTly4SKVquVHTt2ULFiRUqUyLn3fq6KsVj4+eefeeCBB7KWNW3alMaNG2eNRrPZbAQHBzNkyBBGjx6dq+OmpqbSvn17+vXrxxNPPHHdbceMGcPYsWOvWq4JFUVECs6hS4d4Zc0r7LqwC4DW5VvzerPXKe1R+qptNxyJYcxvu9h9OvM/qqEBXrzRpRYtQ0vd1pqlcMvthIp5ahobPnw4s2bNAjJDUKtWrWjQoAHBwcGsWLEiTwXnJC0tjc2bN9OuXbusZQ4ODrRr145169bl6hiGYdC7d2/atm17wxAE8NJLLxEbG5v1iIqKynP9IiJyfRm2DGbtmMVDvz/Ergu78HbxZlzLcUxrO+2qEHTqUjKDv9lCj4/Xsft0HD5uTrzRpSbzh92pECR5lqemsblz5/L4448D8Pvvv3P48GH27t3LV199xSuvvEJkZGS+FBcdHY3VaqVMmTLZlpcpU4a9e/fm6hiRkZF8//331K1bN6v/0VdffUWdOnVy3N7V1RVXV9dbqltERG7scOxhXlvzGn9H/w3AneXuZEzzMQR4BGTbLiXdyierDvPhioOkpNuwWOCRJhV4vn01Snrp97XcmjwFoejoaAIDAwGYP38+PXr0oFq1avTp04epU6fma4G3qmXLlthsNrPLEBGRf1htVv67579M2zKNNFsaXs5ejGoyivur3J+tL5BhGCzceYa35+3h5KXM4fCNK5XgjS61qF1Ow+Elf+QpCJUpU4bdu3cTFBTEwoUL+eijjwBISkrC0dEx34orVaoUjo6OnD17Ntvys2fPZgUxEREpOo7GHuW1yNfYdn4bAC3KtmBM8zEEemb/nb73TBxjf9vNusOZw+GDfN14qXMNutQN0kgwyVd5CkJPPfUUPXr0ICgo8x/kv3141q9fT/Xq1fOtOBcXFxo2bMiyZcuyOlDbbDaWLVvG4MGD8+08IiJSsGyGjW/2fMPULVNJsabg6ezJC41eoHto92zB5lJSGpOX7Oe/fx3DZoCLkwMD7qrMgNZV8HDRbTEk/+XpX9WYMWOoXbs2UVFRPPTQQ1l9ahwdHXM9kutfCQkJHDx4MOv5kSNH2LZtG/7+/lSoUIERI0bQq1cvGjVqRJMmTZgyZQqJiYk89dRTeSk91yIiIoiIiMBqtRboeUREirvjccd5LfI1tpzbAsAdQXfwZvM3CfIKytrGajP4ZsNxJi3ex6V/hsN3qh3Iy51rEOyv4fBScPJl+PytWLFiBW3atLlqea9evZgzZw4AM2bMYMKECZw5c4bw8HCmTZtG06ZNb0t9uR1+JyIi2dkMG9/t/Y4pW6aQnJGMh5MHzzd6noeqPZTtKtBfhy8w5rdd7D0TD0C1Ml6M6VKL5lU1EkzyLrff36YHocJOQUhE5OadiD/B62tfZ+OZjQA0CWzCmy3epJzX/+5HeeJiEuPn72XejtMA+Lo7M6J9NR5rWgEnxzzN7iKSJbff32pwFRGRfGMzbPy470cmbZ5EckYy7k7ujGg4gh5hPXCwZIab5DQrH686xEcrDpGaYcPBAo82rcCI9mH4e7qY/ArE3igIiYhIvjiVcIrX177O+tPrAWhYpiFvtXiLYO/M+zsahsH8HWcYN/9/w+GbhPgzpkstapbVFXcxx00FoT///JNWrVrl6xB5EREp2gzD4P8O/B8TN00kMT0RN0c3hjccziPVH8m6CrTndBxjftvF+iMxAJT1dePle2twbx0Nhxdz3VQQevrpp7l06RIdO3bk/vvvp1OnTsW234xGjYmI3NiZxDO8sfYN1p5aC0D9gPq81eItKvpUBOBiYhqTluzjm/XHsRng6uTAgFZVGNCqCu4u+k+1mO+mO0v//fff/Pbbb/z222/s2LGDli1b0rVrV+6//34qVKhQUHWaRp2lRUSuZhgGvxz8hfc3vk9CegKujq4MrT+Ux2o8hqODIxlWG1+vP87kJfuJTc4cDn9vnSBe6lyd8iU0HF4K3m0ZNXbq1KmsULR8+XLCwsLo2rUrXbt2pVGjRnk9bKGiICQikt3ZxLOMWTeGNSfXAFC3dF3ebvE2Ib4hAKw9FM3Y33az72zmcPjqgd680aUWzaqUNK1msT+3ffh8YmIiCxcu5Ndff2X+/PmMGDGCl19+OT8ObSoFIRGRTIZh8Pvh33l3/bvEp8fj4uDCkPpDeKLmEzg6OBIVk8S4+XtYsPMMAH4ezjzfvhqPNNFweLn9TJ1HyGq1EhMTQ+nSpfP70LedgpCICJxPOs+b695kxYkVANQpVYe3W7xNZb/KJKdZ+WjFQT5edThrOPzjd1RkRPtq+HloOLyYw9R5hBwdHYtFCBIRsXeGYTDvyDzGrx9PXFoczg7ODAwfSO9avXG0OPLb9lOMn7+H07EpANxR2Z83utSiRpD+4yhFg+YREhGRHEUnR/PWurf4M+pPAGqWrMnbLd4mtEQou07FMva33Ww4mjkcvpyfO6/cW4NOtQM1HF6KFAWha9DweRGxV0npSXy5+0tm75xNUkYSTg5OPFvvWZ6q/RTxyQYv/7yD7zZkDod3c3bg2VZVeaZVZdycNRxeih7da+wG1EdIROxFhi2DXw/+SsS2CM4nnwegdsnajGk+hsq+ofz3r2N8sGQ/cSkZANxXN4iXOtegnJ+7mWWL5KjA+wgdOnSI2bNnc+jQIaZOnUpAQAALFiygQoUK1KpVK6+HFRGR28wwDFadWMUHmz/gUOwhAMp5lWN4g+F0qNSByIMXGDxnNQfOJQBQI8iHMV1q0rSyhsNL0ZenILRy5Uo6depEixYtWLVqFe+88w4BAQFs376dWbNmMXfu3PyuU0RECsDO6J1M2jSJTWc3AeDr6suAugPoEdaDUxfTGfDfzSzadRaAEh7OPH9PGI80qYCjg/oBSfGQpyA0evRo3n77bUaMGIG3t3fW8rZt2zJjxox8K05ERApGVHwU07dMZ8HRBQC4OLjweM3H6VunL6djLLzw4y5+334KmwGODhaeuKMiw9uFaji8FDt5CkI7duzgm2++uWp5QEAA0dHRt1yUiIgUjEspl/j474/5bt93ZNgysGChS5UuDA4fzLmL7oz8bj+Ld5/N2r51WGle6lSDsEDv6xxVpOjKUxDy8/Pj9OnThISEZFu+detWypUrly+FiYhI/knJSOHrPV8za8cs4tMzb33RomwLhjcYzqXY0rz4/UFWH8j8j6zFAh1rBTKoTVVql/M1s2yRApenIPTwww8zatQofvzxRywWCzabjcjISEaOHMmTTz6Z3zWKiEgeWW1W5h2Zx/St0zmTmHnri7ASYYxoOILU+Kq8+sNBNh3L7CDt6GDh/vCyDGxdhaoBugIk9iFPw+fT0tIYNGgQc+bMwWq14uTkhNVq5dFHH2XOnDk4Ohb9uSQun0do//79Gj4vIkXO2pNrmbx5Mvsu7gMg0DOQweFDcExqyEcrDrPrVBwALk4O9GhUnmfuqkKwv+4ML8XDbbnX2PHjx9m5cycJCQnUr1+f0NDQvB6q0NI8QiJS1OyN2cvkTZNZd3odAN7O3jxVuy/eqW34ZNVxDp9PBMDDxZHHmlag352VCfBxM7NkkXx3W+41VqFCBSpUqHArhxARkXxyOuE0M7bN4PdDv2Ng4OTgRI/QhymZ0Yk586M5eWkPAD5uTvRuEcJTzStRwlOjwMS+5SkIGYbB3LlzWb58OefOncNms2Vb/9NPP+VLcSIicmNxaXF8tuMzvt79NWm2NADaV+hIkK0bPyxL4nx8FAClvFzo27Iyj99RAW83ZzNLFik08hSEhg8fzscff0ybNm0oU6aMbrAnImKCNGsa3+39jk92fEJsaiwA9Us3JJge/LHKidjkiwCU9XXjmVZV6Nk4WPcDE7lCnoLQV199xU8//UTnzp3zux4REbkBm2Fj0dFFTN0ylZMJJwGo6B1CJYeeLP/Ln1VpNiCdkFKePNuqCg/UL4eLk4O5RYsUUnkKQr6+vlSuXDm/axERkRvYeGYjkzZNYteFXQD4u5aiokN31m8JYWeGBbBRPdCbQW2q0rlOkG6FIXIDeQpCY8aMYezYsXz++ee4u+uuwyIiBe3gxYN8sOUDVp1YBYCbowflHTqzY0c9jlkz+/vUr+DH4DZVaVs9QF0WRHIpT0GoR48efPvttwQEBFCpUiWcnbN3utuyZUu+FCciYu/OJZ3jw20f8vPBn7EZNhwtjpShNQf2NuN8hhcALaqWZFDrqjSrUlIBSOQm5SkI9erVi82bN/P444+rs7SISAFISEtg9q7ZfLnrS1KsKQCUMBoSdag1l9JKA9CuRgAD21SlQYUSZpYqUqTlKQjNmzePRYsW0bJly/yup9C4fGZpEZHbJd2Wztz9c5m5fSYxKTEAuNuqEH38HuKTK+JggfvqZd4Go0aQJnkVuVV5mlm6evXq/PDDD9StW7cgaipUNLO0iNwOhmGw7PgypmyZwrG4YwA428oQd6o9GfG1cHJwoHuDcjzbuiohpTxNrlak8CvQmaUnTZrEiy++yMyZM6lUqVJeaxQREWDrua1M2jSJ7ee3A+Bg8ybp7N3EX2qMq5MzjzevQL+7KlPOT4NTRPJbnoLQ448/TlJSElWqVMHDw+OqztIxMTH5UpyISHF2JPYIU7dMZdnxZZkLDBdSo1uSFtMKL2dPnm5dkb4tQyjl5WpuoSLFWJ6C0JQpU/K5DBER+xGdHM3M7TOZu38uVsMKhoW0S41Ji26Hr0tJBt8dQq9mlfD10G0wRApankeNiYjIzUlKT+LL3V/y+c7ZJGckAZARX53Uc50o5VqB/h0q80iTCni63tL9sEXkJuT60xYXF5fV2SguLu6626pTsYjI/2TYMvj14K9M3zqDCynRAFiTy5N6rhNl3Woz4N4q/KdBed0HTMQEuQ5CJUqU4PTp0wQEBODn55fj3EGGYWCxWDTkXESEzN+Jq06sYsLGSRyLPwKALc2f1HMdqOjejEH3hdK1XlmcHHUfMBGz5DoI/fnnn/j7+wOwfPnyAitIRKQ42Bm9k3HrJrAjJnOmfSPDg9TotlTzaM+Q+2twT81AHHQfMBHT5ToItWrVKuvnkJAQgoODr7oqZBgGUVFR+VediEgRExUfxTtrJxN5ZikAhs2JtJgW1PHsxtBudbkrtJRm4xcpRPLUIy8kJCSrmexyMTExhISEqGlMROzOpZRLvLN2GouifsLAimFYyIitT33vh3mue1OahPibXaKI5CBPQejfvkBXSkhIwM3N7ZaLEhEpKlKtqYxf8wk/H/kKmyUZgIyEUBp6P86o/9xNnfK+JlcoItdzU0FoxIgRAFgsFl577TU8PDyy1lmtVtavX094eHi+FmgW3WtMRK7HMAzmbP+ZGdumkWa5ABawpQTRyPsJXn3wAULLeJtdoojkwk3da6xNmzYArFy5kmbNmuHi4pK1zsXFhUqVKjFy5EhCQ0Pzv1KT6F5jInKltSe28PKKd7hg3Q+Ake5Lfe9HePPuJwkppQAkUhgUyL3G/h0t9tRTTzF16lQFAxGxKyfiTjFi6Xj2xK8AwLA5E+xwHx90GUb1MiVNrU1E8iZPfYRmz56d33WIiBRaSelJjF0Vwfyob8GSDoB7alNeb/E899WqYXJ1InIrNI+7iMg12Awbn2+fy4fbp5POJbAAKSE8UW0II+66WxMhihQDCkIiIjlYFbWe11aNIybjMJA5I3Rz/yeZ8OAT+Hm63GBvESkqFIRERC5zPDaK55e9w974SAAMqyvlLF344L5B1AwqZXJ1IpLfFIRERICEtATGrp7GwqgfwZKBYVhwT2nOqy2Gc3+d6maXJyIFREFIROya1WZl1vbvmPn3h6QTl9kPKDmUJ0KHMKJ1K5zVD0ikWFMQEhG7tfzYGl5fM55LGccBsKWW4g6/3rzXvSelvDVLvog9UBASEbtz+OJhXlj+DvvjNwBgWN0pa3Rl0r0DqVNO9wQTsScKQiJiN2JTYxm7egpLTvwEFhuG4YBb8p280nwoD9QN1V3hReyQgpCIFHvptnQ+2/Y1n+yYSQaJYAEjsQaPhg7i+TYtcXVyNLtEETGJgpCIFFuGYbD06ArGrn2P2IyTAFhTynCHX2/e7f4gAeoHJGL3FIREpFjaH7OfUSve4WD8FgBsGZ4E2R5gYqd+hFfQfcFEJJOC0DVEREQQERGB1Wo1uxQRuQkxKTG8ueYDlp34FSwGhs0Rl6RWvNRsMA/Wr6p+QCKSjcUwDMPsIgqzuLg4fH19iY2NxcfHx+xyROQa0qxpfLztCz7f+SkZJANgS6jDw1WeZWTb5ri7qB+QiD3J7fe3rgiJSJFmGAYLjizmnXUTiMs4C4A1uRyNfXrx7uMPEOTrbnKFIlKYKQiJSJG1K3oXL698h8MJOwCwpXtTxtqN9zv1oVFF9QMSkRtTEBKRIudc0jneipzEipML/ukH5IRzQlteaDqAhxtVxcFB/YBEJHcUhESkyEjJSGHmts+Zs+tzrKSCBaxx4fSo8gwvPHoHnq76lSYiN0e/NUSk0DMMg98PzePd9ZOIz4gGwJpUgQbeTzL+sS4E+3uYXKGIFFUKQiJSqG0/v51XVr3DsYQ9ANjSfSmd9h/Gd3iCZlVKmVydiBR1CkIiUiidTjjNW2snsPr0EgAMmwtOcXfzQtN+PNakKo7qByQi+UBBSEQKlaT0JD7c+ilf7fkCG+kYhgVrXEMeDOnHC480xsfN2ewSRaQYURASkULBZtj4+cAvTNgwlURrDAAZiSGEez7JuEc7EVLK0+QKRaQ4UhASEdNtPLORN9aMJyrxAAC2tJKUTO3O2/c8zF3VAkyuTkSKMwUhETFNVFwUb617n3VnVgBgWF1xiGvP84378OQdVXBydDC3QBEp9hSEROS2i0+LZ8aWj/hu37fYyMAwLGRcakq3Sn0Y9XBD/DxczC5RROyEgpCI3DZWm5Uf98/lg03TSbLGApCREEpt9ycY90g7Qst4m1yhiNgbBSERuS22ntvKK6veJCrxIADW1NKUSPkPb7b/D22rl8Fi0XB4Ebn9FIREpEBFJ0czNvL9zPuCAYbVDculjgxr9ARPNa+Ki5P6AYmIeRSERKRApNvS+Wz7f/nk74/IIDlzPqDYRnSt8DQv9mhISS9Xs0sUEVEQEpH8t+bEX7y88i0uZhwHwJpcnnruT/H2o52pUtrL5OpERP5HQUhE8s2phDOMWPo2u2JXAmDL8CAgoxvj73maOyrrvmAiUvgoCInILUu3pvPmqpn8euwLDEsqhmHBNbk5LzYZTo8GYeoILSKFloLQNURERBAREYHVajW7FJFC7fudy5iw6V1SLWfAAqRU5NEqw3i+dVtcnRzNLk9E5LoshmEYZhdRmMXFxeHr60tsbCw+Pj5mlyNSaPx9+igjlr3JWetGAIwMT5r6PsmETn3w93QzuToRsXe5/f7WFSERuSkxSUkMXzCFLXFzsThk3h2+vGM7Jt87ipqBZcwuT0TkpigIiUiuZFhtvLP8Z/7v6AwM52gsDuBurcord7zC/TUbmV2eiEieKAiJyHUZhsH327YzYeP7pLnuAGewWH14pOpARrV8BAcHTYgoIkWXgpCIXNOmY2d5cclUzjkuwOKaAYYDjfy7Mrn9C5RwV585ESn6FIRE5CrHLyQyeuF3bE/6EgeXGCxAoEttJrYdQ70yYWaXJyKSbxSERCRLbHI67y1dwy/HP8TRay8OLuBKCZ5vNJKHa3bRfEAiUuwoCIkIaRk2Pl+7nw+3fYzNezmOXlYshiNdQh7hleZD8HD2MLtEEZECoSAkYscMw2De36d5Z+X3xHn8hIPvJSxAdd9GvNf6NSr7VTa7RBGRAqUgJGKnNh+L4Y35yzlk+xon3wM4AD5OpXmj+Uu0r9ROzWAiYhcUhETszJHoRMYt2MrKc9/hUnINThYrjjjTq1ZvBoT3w93J3ewSRURuGwUhETsRk5jG1KX7+W73bziVnodrqTgAmpZpyRvNXybYJ9jkCkVEbj8FIZFiLiXdypy1R4lYvZoM/59xKXsYgDLuZXmt2cu0Cm5lcoUiIuZREBIppmw2g1+3n2TC4u1ccP4D5/JrcbLYcHZwoX/dfjxV+ylcHV3NLlNExFQKQiLF0NpD0bwzfzf74lfiWmY+Lk4JALQNvpsXm7xAOa9yJlcoIlI4KAiJFCMHzsbz7oK9LD+yDdfAX3EvdwyACt4VeLnpy7Qo18LkCkVEChcFIZFi4Fx8Ch8sOcD3m/fiXGoJHiF/YbEYuDm68Uy9Z3iy5pO4OLqYXaaISKGjICRShCWlZfDpqiN8vOoAae4bcK+8EAenRAA6VOrAyEYjCfQMNLlKEZHCS0FIpAiy2gzmbo5i0uL9RKcfwq3sr7i7RwFQ2bcyLzV9iTuC7jC5ShGRwk9BSKQIMQyDlfvP8+6Cvew9fwbX0ovx9NsAFgMPJw8Ghg/k0eqP4uzobHapIiJFgoKQSBGx61Qs4+fvZc3Bczj7bcS7yiJwTAKgc0hnnm/0PAEeASZXKSJStCgIiRRyp2OTmbR4P/+35QQW1+N4hvyKg9tJAEJLhPJyk5dpFNjI5CpFRIomBSGRQio+JZ2ZKw8xa80RUm1xuAQuxMVvEwBezl4Mrj+YnmE9cXLQx1hEJK/0G1SkkEnLsPH9xuNMWXqAC4nJOJdYj0+Zpdgsmc1g91e5n+ENh1PKvZTJlYqIFH0KQiKFhM1m8Nv2U0xesp/jMUk4uh/FL/R3rE4nsQE1/GvwctOXCQ8IN7tUEZFiQ0FIxGSGYbBszzkmLt7H3jPxWJxj8K2wDJvnZqyAj4sPQ+sP5cFqD+Lo4Gh2uSIixYqCkIiJ/jp8gfcX7mXL8UtYHOPxLrsCB9/12MjAgoXuod0Z1mAYJdxKmF2qiEixpCAkYoKdJ2N5f9E+Vu0/Dw5JeJRZjYv/WqykYgOaBjXluQbPUatULbNLFREp1hSERG6jg+cSmLxkH/N3nAFLGm6l1uJRejXpJGIF6pSqw9AGQzUrtIjIbaIgJHIbnLyUzNSl+5m7+QQ2IwOXEhvxDlxBGrGkA1X9qjKk/hDaBLfBYrGYXa6IiN0o9kHo0qVLtGvXjoyMDDIyMhg2bBj9+vUzuyyxExcSUolYfoj//nWMNGsGTr5b8S+7nFSiSQPKeZVjUPggOod0VkdoERETFPsg5O3tzapVq/Dw8CAxMZHatWvTvXt3SpYsaXZpUozFp6Tz6eojzFp9mMS0DJy8d1Gq7DJSHU6TCpRyL8WAugPoHtpd9wUTETFRsQ9Cjo6OeHh4AJCamophGBiGYXJVUlylpFv5at0xPlxxkItJ6Th6HKBkpaWkOR4jlcyh8H1q9+HRGo/i7uRudrkiInbPwewCVq1aRZcuXShbtiwWi4Vffvnlqm0iIiKoVKkSbm5uNG3alA0bNtzUOS5dukS9evUoX748L7zwAqVKaUZeyV8ZVhvfbjhO6wkreGf+HmJth/Cv8jkeFWeR5ngMdyd3+tXpx4L/LKBvnb4KQSIihYTpV4QSExOpV68effr0oXv37let//777xkxYgQzZ86kadOmTJkyhQ4dOrBv3z4CAjLvtB0eHk5GRsZV+y5evJiyZcvi5+fH9u3bOXv2LN27d+fBBx+kTJkyBf7apPiz2Qzm7TjN5CX7ORKdiIPrGUqELCPDbQfpgLODMz3CevB0nad1SwwRkULIYhSidiKLxcLPP//MAw88kLWsadOmNG7cmBkzZgBgs9kIDg5myJAhjB49+qbPMXDgQNq2bcuDDz6Y4/rU1FRSU1OznsfFxREcHExsbCw+Pj43fT4pngzDYMX+80xYuI/dp+OwOF/AJ+hPbJ5bAAMHiwNdq3Tl2XrPUtarrNnliojYnbi4OHx9fW/4/W36FaHrSUtLY/Pmzbz00ktZyxwcHGjXrh3r1q3L1THOnj2Lh4cH3t7exMbGsmrVKp599tlrbj9+/HjGjh17y7VL8bXpaAzvL9zHhqMxWJzi8C63HIvPBmxYAWhfsT2DwwdT2a+yyZWKiMiNFOogFB0djdVqvaoZq0yZMuzduzdXxzh27Bj9+/fP6iQ9ZMgQ6tSpc83tX3rpJUaMGJH1/N8rQiK7T8UxcfE+/tx7DhwTL5sNOg0DaFG2BUMaDKFWSc0GLSJSVBTqIJQfmjRpwrZt23K9vaurK66urgVXkBQ5R6MTmbxkP79tPwWWVFxLReJRejUZJGMFwkuHM7TBUBoHNja7VBERuUmFOgiVKlUKR0dHzp49m2352bNnCQwMNKkqsRdnYlOYuuwAP2yKwmqk4VxiPd6BK0knngygWolqDGswjDvL3anZoEVEiqhCHYRcXFxo2LAhy5Yty+pAbbPZWLZsGYMHDza3OCm2Liam8dHKQ3yx9iipGek4+W6hZNBy0iwxpAMVvCswKHwQHUM64mAxfQYKERG5BaYHoYSEBA4ePJj1/MiRI2zbtg1/f38qVKjAiBEj6NWrF40aNaJJkyZMmTKFxMREnnrqqQKtKyIigoiICKxWa4GeRwqPxNQMZq05wqerDhOfmoaT905Kll1GmsNZ0oAAjwAG1BvAA1UfwNlBs0GLiBQHpg+fX7FiBW3atLlqea9evZgzZw4AM2bMYMKECZw5c4bw8HCmTZtG06ZNb0t9uR1+J0VXaoaVr/86TsTyg1xITMXRcz++ZZeS7hQFgJ+rH0/XeZqeYT1xc3IzuVoREcmN3H5/mx6ECjsFoeIrw2rjp60nmbr0ACcvJePofhSfskvIcDkEgIeTB71q9eLJmk/i5eJlcrUiInIzisU8QiIFwTAMFu48w8TF+zh0PhEH11P4hSzF6rabDMDFwYWHqz9M3zp98XfzN7tcEREpQApCYjcMw2DNwWgmLNrH3ydisThH41NhKYbnNqyAo8WRB6o+wIB6Awj01KhEERF7oCAkdmHL8YtMWLiPdYcvYHGKxavcnzj4bMTABkCnSp0YGD6QSr6VzC1URERuKwWha9CoseJh35l4Ji7ex5LdZ7E4JuJeZgUu/n9hIx0DuKv8XQypP4Tq/tXNLlVEREygztI3oM7SRVNUTBIfLNnPz9tOYlhScPVfg3vpNVhJAaBBQAOGNRhGgzINTK5UREQKgjpLi106F5fCjOUH+XbDcdJtaTiX+AuvMivJIAErUMO/BkMbDKVF2RaaDVpERBSEpHiITUpn5qpDzI48Qkp6Os5+m/APXEG65SIZQCWfSgyuP5j2FdtrNmgREcmiICRFWlJaBrMjj/LxykPEpaTh5PM3/pWXke5wnnQg0DOQgfUG0qVKF5wc9M9dRESy0zeDFElpGTa+23ic6X8e5Hx8Co5e+ygRuoQMp5OkA/5u/vSr04+Hwh7C1dHV7HJFRKSQUhCSIsVqM/h120k+WLqfqJhkHN0P41dlCVaXI2QAXs5e9K7Vm8drPo6ns6fZ5YqISCGnIHQNGj5fuBiGwZLdZ5m4eB/7zybg4HYS35Al2Nz2YgVcHV15tPqj9KndBz83P7PLFRGRIkLD529Aw+fNt/ZQ5mzQW49fwsHlHJ6BS8HzbwCcLE78p9p/6F+3PwEeASZXKiIihYWGz0uR9/eJS0xYtI/VB6KxOF3Cs9wyHH02YWBgwcK9le9lYL2BBPsEm12qiIgUUQpCUugcPBfPpMX7WbDzDBbHBNwDl+NSYj02MjCANsFtGFx/MNVKVDO7VBERKeIUhKTQOHExialLD/B/W05gs6TgWnoV7qUisZKKDWgS2IShDYZSr3Q9s0sVEZFiQkFITHcuPoUPlx/im/XHSbOl4FJiHR5lVmElEStQu2RthjYYyh1Bd2g2aBERyVcKQmKaCwmpfLzqMF+uO0pKehrOfpsoEbicDEssVqCKbxWG1B9C2wptFYBERKRAKAjJbXcxMY1PVx9mztqjJKWl4+TzNyUqLyXDIZoMoJxXOQaGD+TekHtxdHA0u1wRESnGFISuQfMI5b/Y5HRmrTnC52uOkJCajqPnfkpWWkKa4wkygJJuJelftz8PVnsQF0cXs8sVERE7oHmEbkDzCN26+JR0Zkce5dPVh4lPycDB7Tj+5ZeQ6nwAyJwN+qnaT/F4jcfxcPYwuVoRESkONI+QmC4xNYMv1h3lk1WHuZSUjoPLOUpXWUaKy3ZSAWcHZx6p/gj96vTTbNAiImIKBSHJd8lpVr5ef4yPVhziQmIaFqdYSoesINVtPSnYcLA40KVyFwaFDyLIK8jsckVExI4pCEm+SUm38u2G43y44hDn41PBIYnSFSLJ8FpNipEGZE6GOLT+UKqWqGpytSIiIgpCkg/SMmx8vymKiD8PciYuBSxplCq/AYvvClJsCWBAg4AGPNfwOcIDws0uV0REJIuCkORZutXG/20+wfQ/D3LyUjJgpVTQdpxLLSUhIwZsEFoilOENhnNnuTs1F5CIiBQ6CkJy0zKsNn7Zdoppyw5wPCYJMCgZsBfPwKVcTD9JagaU9SzL4PqD6RzSWXMBiYhIoaUgJLlmtRn88fcppi49wOHoRAD8Sx6lRPBSzqUeJC0dSriWoH/d/vQI66G5gEREpNBTEJIbstkMFuw8w5Sl+zlwLgEAP7+zlA35k6iU7ZxLBQ8nD3rV6kWvWr3wdPY0uWIREZHcURC6Bs0sDYZhsHj3WT5Ysp+9Z+IB8Pa6RJVqqzmUHElUCjg5ONGjWg/61+1PSfeSJlcsIiJyczSz9A3Y48zShmGwfN85Ji/Zz86TcQB4eyRRvcZ6DiQvw2pkYMFC58qdGRQ+iGDvYJMrFhERyU4zS8tNMwyD1QeimbxkP9uiLgHg4ZZGeO1tHEiZz96kFABalmvJ8AbDCfMPM7FaERGRW6cgJACsPRTNB0v2s/HoRQDcXKw0qbuXQ+m/sSMxFoC6pesyvMFwGgc2NrNUEREp7AwD0pMgJQ5SYiH1nz//fWQ9/+fPBz4EJ1dTSlUQsnMbj8YwefF+1h2+AICLE7QIP8px289sTTwLQGXfygxtMJS2wW01F5CIiD2wZvwTVi7lEGauF24uW2fLyP35Oo4Hr4ACeznXoyBkp7Yev8jkJftZfSAaABdHC3fVP8cZh5/YFH8EgDIeZRgUPoguVbrg5KB/KiIiRYJhQFriNYLLpVwEmzhIT8yfWiwO4OYLrj7g5gNufv/87PvP83/WObnlz/nyQN9udmbHiVg+WLqfP/eeA8DJwULb8ERi3X9hfczfAPi4+NCvTj8erv4wbib+4xQRkX8YBsRGwcnNcH4fJF+6/lUZI59GPDt7Zg8s2X6+PMz4Xh1u3HzBxRMKeUuCgpCd2HM6jg+W7Gfx7szmLkcHC+3DraR5z2PduUhIBjdHNx6v+ThP1X4KHxf7GCEnIlIoJcXAqS1wcktm+Dm5GRLP39wxHJxyvvri5nuNMHP5cz9w9QZH5wJ5eYWJglAxt/9sPFOXHmDejtMAOFigQ11nnEovYeXJRRjJBo4WR/4T+h8G1BtAaY/SJlcsImJn0pPh9N//BJ9/Qk/M4au3c3CCMrUhsA54lro62FwZZpw9Cv3VmMJAQaiYOnQ+gWnLDvDb9lMYRuZnoX1tT/zLrWJR1M9knMzsxNahUgeG1B9CRZ+KJlcsImIHbNbMpq1/A8/JzXBud84di/2rQLmG/3sE1gFndVfIbwpCxcyxC4lMW3aQn7eewPbPVJntavpSsfIm/jj2LUnHkgC4I+gOhjccTq2StUysVkSkGDMMiD1xWejZAqe25twR2TPgstDTAMrWBw//21+zHVIQKiZOXExixp8H+XHzCaz/JKC2NfypWW03vx9/n/WHYgCoWbImwxsMp1nZZmaWKyJS/CTFZAadbP16zl29nYtXZtAp1wDKNsgMP77l1YxlEgWhaygq9xo7HZtMxPKDfL8xinRrZgBqFVaSJrWP8fvx6WzcfxKAij4VGVJ/CO0rtsfB4mBmySLFi2FkjtRx8QYHfbbsRnoynNmRvYnrmv16amVv4ipVDRwcb3/NkiPda+wGCuu9xs7FpfDhikN8s/44aVYbAC2qluTuBjHMOzGL/Rf3A1DavTQD6g2gW2g3nB2Kf+9/kQJhs0HcycwvumyPI3DxSOYMuhbHzAnhvAIymzm8yvzvudc/zz3/+dnNV//7L0psVojenz30nN11jX49lXPo1+N++2sW3WusuIpOSGXmikN89dcxUjMyA1CTEH/ub5rOkjMz+WDHZgC8nb3pU6cPj9V4DHcnfQhFbsiakTlPy+Uh59+fLx4Fa+r19zesEH8683Ejjq5XBKTSVwSnMv8LVK5e+fLyJJcu79fz7/D1U1shLeHqbT1LQ7lG6tdTxCkIFREXE9P4ZPVhvlh7lKS0zOa6BhX8eLSlG6ujv+Td7X8C4OLgwqM1HqVv7b74ufmZWLFIIZSRBpeO5Rx2Lh27/i0BHJyhRKXM//Fne4SAT9nMCe4SzkLCucx+If/+nHAWEs7/73lqbGaoio3KfNyIs+e1ryx5lfnnUTpzmUYU3bzki//06blBvx5nz//16/n3ao/69RQLCkKFXGxSOrPWHObzyKMkpGb+kq5X3pfed/mxJe573tr2GzbDhoPFgfur3M/A8IEEegaaXLWIidKTM6/gXNWMdTjzf/qG7dr7OrlBiZD/BZzLA49v+ev363B2B5+gXNSX8k9QOndZUPrnz8QrlqUnZY4wuvhPE9yNuPneuFnOq0zmHDR2MFHeVdJTcujXc+jq7f7t11P2stBTOkz9eoopBaFCKj4lnc/XHOWzNYeJT8kMQDWDfHimTRn2p/7K29u/I82WBsDdFe5maP2hVParbGbJIrdPakJmMLiyv07M4cy+PNfj7Jlz0PGvDN5BBd/h2dkN/CpkPm4kNeGfgHT+iitMlwWmf9dZ0/53m4ULB258bI+SN26W8yqTGa4MW+bVMsOa2V/GlnHZnxn/W5/jssuW57TMZv3nuJcvv3LZZftlLc9p2WX1XVlr4jn165EcKQgVMompGcxZe5RPVx/mUlI6AGFlvBnUtjxnWMp7O0cRnx4PQKMyjRjecDj1Stczs2SRgpF86Yqwc9nPCWevv6+rTw5NWP88vAKKTnOGq1fmo2SV629nGJk307zyKtOVV5gSzmUGJ8MKSRcyH/bGs/QV8/U0UL8eO6cgVEgkp1n56q+jzFx5mJjEzCs9VUp7MuTuyqS4reODv98gOjnzTvFhJcIY1mAYLcu1xFJUfqGLXMkwMvtn5NSEFXP4xl/S7v7XDjse/kUn7OQHiwXcS2Q+Soddf1ubNXO+m8Qcri5deZXpen8HFsfMJiSHy/7MWuaUeWXNwSnnZVnLrzhGTsuuOq7jFdvnYpmD0z9z94SDb7B9/duQG1IQMllKupVvNxznwxWHOB+fOSqlUkkPhtxdBXe/3Xy4fSDH4o4BUM6rHIPrD6ZzSGfNBVSc2GyZTQaG9Z8/bZc1DVy2LmuZ9bJtrlxny75NTuuyHf+K51ed25rDuhz2z2n7rHVXbG9Lz+yrE3M4sxnnejwDru6Y/O+f7iVuz99PcePgmNm52qt0Zj+Y67GmQ2r8FSHECSwOChNSbCgImSQ1w8oPm04Q8edBzsSlAFC+hDtD7w6lbGAU07e9wO4duwHwd/PnmbrP8FC1h3AuKh0cbbbMkTHpyZCRChkplz2uXJ4KGcmZv3QN4+ov3azntv99sV61jS2H51fud419zD6XvfMue40+OyGZd78W8zg6q9lIij0FIZMM/KwjcelJZCQ9RlnfMAa3DaVmpVgitr/JX7v+AsDDyYPetXvzZM0n8XT2vPmTGMb/Qsa/oSP9sjCSm3By1X457J+ew/Gsafn8jtk5i8P/mg2y/nTIfGRb5pjZ/GBxvGKdwxX7XfbzVce8cvvLt8vhfBaHq7e/bn2OmZ2S/StnDkd38TD73RURO6YgZJIDLqe46O7Ac+77uauhNx9Ffcb4vTsBcMKBh/1q0s+nJv7R0fDnuKuDyVXB5YqAkp5y4wngbheLY+YIDCfXzOHJWY9/njv/8zyrP8CVX+SW/z2//AvZ8s+6K/fJab+rtslhv8u/vK96fvkX/rX2uUY9Dped81r1ODjkcK7LQoyIiBQIBSGTWMhsX9+S/gfTti/FarFgMQzuS0hi0KVLlDtyFJifr2fMDCOXhZDchBMnV3Byv8Zyt+yPa61z1D8zEREpnPQNZRZnd7ClsNIjc66KVjZXhjiUIszfFwJuMpzkJrg4OKlzo4iIyBUUhEzi71OemEsHCS8dzvCGw2lYpqHZJYmIiNgdBaFriIiIICIiAqvVWiDHn3H3DM4mnqV+QH3NBSQiImISi2EYhtlFFGZxcXH4+voSGxuLj4+P2eWIiIhILuT2+1vDUURERMRuKQiJiIiI3VIQEhEREbulICQiIiJ2S0FIRERE7JaCkIiIiNgtBSERERGxWwpCIiIiYrcUhERERMRuKQiJiIiI3VIQEhEREbulICQiIiJ2S0FIRERE7JaT2QUUdoZhAJl3sRUREZGi4d/v7X+/x69FQegG4uPjAQgODja5EhEREblZ8fHx+Pr6XnO9xbhRVLJzNpuNatWqsXnzZiwWS672ady4MRs3brzuNnFxcQQHBxMVFYWPj09+lFrk5eZ9M9Ptrq+gzpdfx72V4+Rl35vZJ7fb6nOYnT6Dt+d89vAZzO32BfkZNAyD+Ph4ypYti4PDtXsC6YrQDTg4OODi4nLdNHklR0fHXP+F+vj46BfwP27mfTPD7a6voM6XX8e9lePkZd+b2edmj6/PYSZ9Bm/P+ezhM3iz2xfUZzA3393qLJ0LgwYNKtDtJVNhf99ud30Fdb78Ou6tHCcv+97MPoX931JhVdjfN30G8+84Bf0ZzOs5zKCmMZPExcXh6+tLbGxsof4fmEhxps+hiLkKw2dQV4RM4urqyhtvvIGrq6vZpYjYLX0ORcxVGD6DuiIkIiIidktXhERERMRuKQiJiIiI3VIQEhEREbulICQiIiJ2S0FIRERE7JaCUCH0xx9/EBYWRmhoKJ999pnZ5YjYpW7dulGiRAkefPBBs0sRsUtRUVG0bt2amjVrUrduXX788ccCOY+GzxcyGRkZ1KxZk+XLl+Pr60vDhg1Zu3YtJUuWNLs0EbuyYsUK4uPj+eKLL5g7d67Z5YjYndOnT3P27FnCw8M5c+YMDRs2ZP/+/Xh6eubreXRFqJDZsGEDtWrVoly5cnh5edGpUycWL15sdlkidqd169Z4e3ubXYaI3QoKCiI8PByAwMBASpUqRUxMTL6fR0Eon61atYouXbpQtmxZLBYLv/zyy1XbREREUKlSJdzc3GjatCkbNmzIWnfq1CnKlSuX9bxcuXKcPHnydpQuUmzc6udQRG5dfn4ON2/ejNVqJTg4ON/rVBDKZ4mJidSrV4+IiIgc13///feMGDGCN954gy1btlCvXj06dOjAuXPnbnOlIsWXPoci5suvz2FMTAxPPvkkn3zyScEUakiBAYyff/4527ImTZoYgwYNynputVqNsmXLGuPHjzcMwzAiIyONBx54IGv9sGHDjK+//vq21CtSHOXlc/iv5cuXG//5z39uR5kixVpeP4cpKSnGnXfeaXz55ZcFVpuuCN1GaWlpbN68mXbt2mUtc3BwoF27dqxbtw6AJk2asHPnTk6ePElCQgILFiygQ4cOZpUsUuzk5nMoIgUrN59DwzDo3bs3bdu25YknniiwWhSEbqPo6GisVitlypTJtrxMmTKcOXMGACcnJyZNmkSbNm0IDw/n+eef14gxkXyUm88hQLt27XjooYeYP38+5cuXV0gSyUe5+RxGRkby/fff88svvxAeHk54eDg7duzI91qc8v2Icsu6du1K165dzS5DxK4tXbrU7BJE7FrLli2x2WwFfh5dEbqNSpUqhaOjI2fPns22/OzZswQGBppUlYh90edQxHyF6XOoIHQbubi40LBhQ5YtW5a1zGazsWzZMpo1a2ZiZSL2Q59DEfMVps+hmsbyWUJCAgcPHsx6fuTIEbZt24a/vz8VKlRgxIgR9OrVi0aNGtGkSROmTJlCYmIiTz31lIlVixQv+hyKmK/IfA4LbDyanVq+fLkBXPXo1atX1jbTp083KlSoYLi4uBhNmjQx/vrrL/MKFimG9DkUMV9R+RzqXmMiIiJit9RHSEREROyWgpCIiIjYLQUhERERsVsKQiIiImK3FIRERETEbikIiYiIiN1SEBIRERG7pSAkIiIidktBSERuu9atWzN8+HCzy8hiGAb9+/fH398fi8XCtm3bzC4pS6VKlZgyZYrZZYgUW7rXmIjYvYULFzJnzhxWrFhB5cqVKVWqlNklichtoiAkIsWC1WrFYrHg4HDzF7oPHTpEUFAQzZs3L4DKRKQwU9OYiJ1q3bo1Q4cO5cUXX8Tf35/AwEDGjBmTtf7o0aNXNRNdunQJi8XCihUrAFixYgUWi4VFixZRv3593N3dadu2LefOnWPBggXUqFEDHx8fHn30UZKSkrKdPyMjg8GDB+Pr60upUqV47bXXuPzWh6mpqYwcOZJy5crh6elJ06ZNs84LMGfOHPz8/Pjtt9+oWbMmrq6uHD9+PMfXunLlSpo0aYKrqytBQUGMHj2ajIwMAHr37s2QIUM4fvw4FouFSpUq5XiMPn36ULduXVJTUwFIS0ujfv36PPnkkzlu/8knn1C2bFlsNlu25ffffz99+vQBMgPY/fffT5kyZfDy8qJx48YsXbo0x+NB7v5OAHbu3EmnTp3w8vKiTJkyPPHEE0RHR2etnzt3LnXq1MHd3Z2SJUvSrl07EhMTr3lekeJMQUjEjn3xxRd4enqyfv163n//fd58802WLFly08cZM2YMM2bMYO3atURFRdGjRw+mTJnCN998w7x581i8eDHTp0+/6txOTk5s2LCBqVOnMnnyZD777LOs9YMHD2bdunV89913/P333zz00EN07NiRAwcOZG2TlJTEe++9x2effcauXbsICAi4qraTJ0/SuXNnGjduzPbt2/noo4+YNWsWb7/9NgBTp07lzTffpHz58pw+fZqNGzfm+BqnTZtGYmIio0ePBuCVV17h0qVLzJgxI8ftH3roIS5cuMDy5cuzlsXExLBw4UIee+wxABISEujcuTPLli1j69atdOzYkS5dulwz0OXGpUuXaNu2LfXr12fTpk0sXLiQs2fP0qNHDwBOnz7NI488Qp8+fdizZw8rVqyge/fu6P7bYrdu+/3uRaRQaNWqldGyZctsyxo3bmyMGjXKMAzDOHLkiAEYW7duzVp/8eJFAzCWL19uGIZhLF++3ACMpUuXZm0zfvx4AzAOHTqUteyZZ54xOnTokO3cNWrUMGw2W9ayUaNGGTVq1DAMwzCOHTtmODo6GidPnsxW391332289NJLhmEYxuzZsw3A2LZt23Vf58svv2yEhYVlO1dERITh5eVlWK1WwzAM44MPPjAqVqx43eMYhmGsXbvWcHZ2Nl577TXDycnJWL169XW3v//++40+ffpkPf/444+NsmXLZp03J7Vq1TKmT5+e9bxixYrGBx98YBhG7v5O3nrrLeOee+7JdsyoqCgDMPbt22ds3rzZAIyjR4/e8PWK2ANdERKxY3Xr1s32PCgoiHPnzt3SccqUKYOHhweVK1fOtuzK495xxx1YLJas582aNePAgQNYrVZ27NiB1WqlWrVqeHl5ZT1WrlzJoUOHsvZxcXG56jVcac+ePTRr1izbuVq0aEFCQgInTpy4qdfZrFkzRo4cyVtvvcXzzz9Py5Ytr7v9Y489xv/93/9lNad9/fXXPPzww1n9mBISEhg5ciQ1atTAz88PLy8v9uzZc0tXhLZv387y5cuzvW/Vq1cHMpvi6tWrx913302dOnV46KGH+PTTT7l48WKezydS1KmztIgdc3Z2zvbcYrFk9Wn598vauKzJJD09/YbHsVgs1z1ubiQkJODo6MjmzZtxdHTMts7LyyvrZ3d392wBp6DZbDYiIyNxdHTk4MGDN9y+S5cuGIbBvHnzaNy4MatXr+aDDz7IWj9y5EiWLFnCxIkTqVq1Ku7u7jz44IOkpaXleLzc/J0kJCTQpUsX3nvvvav2DwoKwtHRkSVLlrB27dqsJstXXnmF9evXExISkqv3QaQ40RUhEclR6dKlgcw+Jf/Kz/l11q9fn+35X3/9RWhoKI6OjtSvXx+r1cq5c+eoWrVqtkdgYOBNnadGjRqsW7cuW3iIjIzE29ub8uXL39SxJkyYwN69e1m5ciULFy5k9uzZ193ezc2N7t278/XXX/Ptt98SFhZGgwYNstXRu3dvunXrRp06dQgMDOTo0aPXPF5u/k4aNGjArl27qFSp0lXvnaenJ5AZTFu0aMHYsWPZunUrLi4u/Pzzzzf1XogUFwpCIpIjd3d37rjjDt5991327NnDypUrefXVV/Pt+MePH2fEiBHs27ePb7/9lunTpzNs2DAAqlWrxmOPPcaTTz7JTz/9xJEjR9iwYQPjx49n3rx5N3WegQMHEhUVxZAhQ9i7dy+//vorb7zxBiNGjLipofZbt27l9ddf57PPPqNFixZMnjyZYcOGcfjw4evu99hjjzFv3jw+//zzrE7S/woNDeWnn35i27ZtbN++nUcfffS6V85y83cyaNAgYmJieOSRR9i4cSOHDh1i0aJFPPXUU1itVtavX8+4cePYtGkTx48f56effuL8+fPUqFEj1++FSHGiICQi1/T555+TkZFBw4YNGT58eNZIq/zw5JNPkpycTJMmTRg0aBDDhg2jf//+Wetnz57Nk08+yfPPP09YWBgPPPAAGzdupEKFCjd1nnLlyjF//nw2bNhAvXr1GDBgAH379r2pUJeSksLjjz9O79696dKlCwD9+/enTZs2PPHEE1it1mvu27ZtW/z9/dm3bx+PPvpotnWTJ0+mRIkSNG/enC5dutChQ4dsV4xycqO/k7JlyxIZGYnVauWee+6hTp06DB8+HD8/PxwcHPDx8WHVqlV07tyZatWq8eqrrzJp0iQ6deqU6/dDpDixGIbGTIqIiIh90hUhERERsVsKQiIiImK3FIRERETEbikIiYiIiN1SEBIRERG7pSAkIiIidktBSEREROyWgpCIiIjYLQUhERERsVsKQiIiImK3FIRERETEbikIiYiIiN36f8fa+QXpBTwNAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "for key, val in t.items():\n", + " plt.plot(n, val, label=f\"diagonal={key}\" if key in (False, True) else \"without optimization\")\n", + "plt.loglog()\n", + "plt.legend()\n", + "plt.xlabel(\"number of x values\")\n", + "plt.ylabel(\"time / sec\");" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3.10.7 ('venv': venv)", + "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.7" + }, + "orig_nbformat": 4, + "vscode": { + "interpreter": { + "hash": "d2f38bb2cd6c72c971581fd65aa33d771e91734dcc3715df26c859a426a25f42" + } + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/src/jacobi/_jacobi.py b/src/jacobi/_jacobi.py index 75f465d..0c7a28c 100644 --- a/src/jacobi/_jacobi.py +++ b/src/jacobi/_jacobi.py @@ -154,7 +154,7 @@ def jacobi( # more like student's t rei = c[-1, -1] ** 0.5 - # update estimates that have significantly smaller error + # update estimates that have smaller estimated error sub_todo = rei < re[todo] todo1 = todo.copy() todo[todo1] = sub_todo diff --git a/src/jacobi/_propagate.py b/src/jacobi/_propagate.py index 070c388..33d81c9 100644 --- a/src/jacobi/_propagate.py +++ b/src/jacobi/_propagate.py @@ -47,6 +47,15 @@ def propagate( If ycov is a matrix, unless y is a number. In that case, ycov is also reduced to a number. + Notes + ----- + For callables `fn` which perform only element-wise computation, the jacobian is + a diagonal matrix. This special case is detected and the computation optimised, + although can further speed up the computation by passing the argumet `diagonal=True`. + + In this special case, error propagation works correctly even if the output of `fn` + is NaN for some inputs. + Examples -------- General error propagation maps input vectors to output vectors:: @@ -135,7 +144,7 @@ def _propagate_full(fn, y: np.ndarray, x: np.ndarray, xcov: np.ndarray, **kwargs _check_x_xcov_compatibility(x_a, xcov) - jac = np.asarray(jacobi(fn, x_a, **kwargs)[0]) + jac = jacobi(fn, x_a, **kwargs)[0] y_len = len(y) if y.ndim == 1 else 1 @@ -143,6 +152,11 @@ def _propagate_full(fn, y: np.ndarray, x: np.ndarray, xcov: np.ndarray, **kwargs jac = jac.reshape((y_len, len(x_a))) assert np.ndim(jac) == 2 + # Check if jacobian is diagonal, count NaN as zero. + # This is important to speed up the product below and + # to get the right answer for covariance matrices that + # contain NaN values. + jac = _try_reduce_jacobian(jac) ycov = _jac_cov_product(jac, xcov) if y.ndim == 0: @@ -156,7 +170,7 @@ def _propagate_diagonal(fn, y: np.ndarray, x: np.ndarray, xcov: np.ndarray, **kw _check_x_xcov_compatibility(x_a, xcov) - jac = np.asarray(jacobi(fn, x_a, **kwargs)[0]) + jac = jacobi(fn, x_a, **kwargs)[0] assert jac.ndim <= 1 ycov = _jac_cov_product(jac, xcov) @@ -187,7 +201,7 @@ def wrapped(x, *rest): xcov = xcov_parts[i] _check_x_xcov_compatibility(x_a, xcov) - jac = np.asarray(jacobi(wrapped, x_a, *rest, **kwargs)[0]) + jac = jacobi(wrapped, x_a, *rest, **kwargs)[0] ycov += _jac_cov_product(jac, xcov) if y.ndim == 0: @@ -198,12 +212,16 @@ def wrapped(x, *rest): def _jac_cov_product(jac: np.ndarray, xcov: np.ndarray): if xcov.ndim == 2: - return np.einsum( - "i,j,ij -> ij" if jac.ndim == 1 else "ij,kl,jl", jac, jac, xcov - ) - elif jac.ndim == 2: + if jac.ndim == 2: + return np.einsum("ij,kl,jl", jac, jac, xcov) + if jac.ndim == 1: + return np.einsum("i,j,ij -> ij", jac, jac, xcov) + return jac**2 * xcov + assert xcov.ndim < 2 + if jac.ndim == 2: if xcov.ndim == 1: return np.einsum("ij,kj,j", jac, jac, xcov) + assert xcov.ndim == 0 # xcov.ndim == 2 is already covered above return np.einsum("ij,kj", jac, jac) * xcov assert jac.ndim < 2 and xcov.ndim < 2 return xcov * jac**2 @@ -213,3 +231,23 @@ def _check_x_xcov_compatibility(x: np.ndarray, xcov: np.ndarray): if xcov.ndim > 0 and len(xcov) != (len(x) if x.ndim == 1 else 1): # this works for 1D and 2D xcov raise ValueError("x and cov have incompatible shapes") + + +def _try_reduce_jacobian(jac: np.ndarray): + if jac.ndim != 2 or jac.shape[0] != jac.shape[1]: + return jac + # if jacobian contains only off-diagonal elements + # that are zero or NaN, we reduce it to diagonal form + ndv = _nodiag_view(jac) + m = np.isnan(ndv) + ndv[m] = 0 + if np.count_nonzero(ndv) == 0: + return np.diag(jac) + return jac + + +def _nodiag_view(a: np.ndarray): + # https://stackoverflow.com/a/43761941/ @Divakar + m = a.shape[0] + p, q = a.strides + return np.lib.stride_tricks.as_strided(a[:, 1:], (m - 1, m), (p + q, q)) diff --git a/test/test_propagate.py b/test/test_propagate.py index c30ba28..4ca395f 100644 --- a/test/test_propagate.py +++ b/test/test_propagate.py @@ -78,20 +78,23 @@ def fn(x): @pytest.mark.parametrize("ndim", (1, 2)) -def test_cov_1d_2d(ndim): +@pytest.mark.parametrize("diagonal", (False, True)) +@pytest.mark.parametrize("len", (1, 2)) +def test_cov_1d_2d(ndim, diagonal, len): def fn(x): - return x + return 2 * x - x = [1, 2] - xcov_1d = [3, 4] - xcov_2d = np.diag(xcov_1d) + x = [1, 2][:len] + xcov = [3, 4][:len] + if ndim == 2: + xcov = np.diag(xcov) - y, ycov = propagate(fn, x, xcov_1d if ndim == 1 else xcov_2d) + y, ycov = propagate(fn, x, xcov, diagonal=diagonal) - assert np.ndim(ycov) == 2 + assert np.ndim(ycov) == ndim - assert_allclose(y, x) - assert_allclose(ycov, xcov_2d) + assert_allclose(y, np.multiply(x, 2)) + assert_allclose(ycov, np.multiply(xcov, 4)) def test_two_arguments_1(): @@ -224,11 +227,6 @@ def fn(x): y, ycov = propagate(fn, x, xcov, diagonal=True) - # Beware: this produces a matrix with all NaNs - # y_ref, ycov_ref = propagate(fn, x, xcov) - # The derivative cannot figure out that the off-diagonal elements - # of the jacobian are zero. - y_ref = [2, np.nan, 5] assert_allclose(y, y_ref) @@ -239,3 +237,10 @@ def fn(x): # ycov_ref = jac @ np.array(xcov) @ jac.T ycov_ref = [[12, np.nan, 8], [np.nan, np.nan, np.nan], [8, np.nan, 80]] assert_allclose(ycov, ycov_ref) + + # propagate now detects the special case where jac is effectively diagonal + # and does the equivalent of propagate(fn, x, xcov, diagonal=True), which + # is nevertheless faster + y2, ycov2 = propagate(fn, x, xcov) + assert_allclose(y2, y_ref) + assert_allclose(ycov2, ycov_ref)