heat_equation_pytensor.ipynb

{
  "cells" : [ {
    "metadata" : {
      "ExecuteTime" : {
        "end_time" : "2025-03-17T18:08:20.197470Z",
        "start_time" : "2025-03-17T18:08:19.379418Z"
      }
    },
    "cell_type" : "code",
    "source" : [ "import pytensor\n", "import pytensor.tensor as pt\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "L = 1.0\n", "T = 0.1\n", "\n", "def build_heat_equation_fn(mode=None):\n", "    nx = 50\n", "    nt = pt.scalar(\"nt\", dtype=int)\n", "    dx = L / (nx - 1)\n", "    dt = T / nt\n", "\n", "    alpha = pt.scalar(\"alpha\")\n", "    r = alpha * dt / dx**2\n", "\n", "    x = np.linspace(0, L, nx)\n", "    u = np.sin(np.pi * x)\n", "    # Check bounds are zero\n", "    np.testing.assert_allclose(u[[0, -1]], 0, atol=1e-15)\n", "\n", "    A = pt.zeros((nx - 2, nx - 2))\n", "    row_idxs, col_idxs = np.diag_indices(nx - 2)\n", "    A = A[row_idxs, col_idxs].set(1 + 2 * r)\n", "    A = A[row_idxs[:-1], col_idxs[1:]].set(-r)\n", "    A = A[row_idxs[1:], col_idxs[:-1]].set(-r)\n", "\n", "    assume_a = \"tridiagonal\"\n", "    if mode.lower() == \"numba\":\n", "        # Next best thing\n", "        assume_a = \"symmetric\"\n", "\n", "    u_history, _ = pytensor.scan(\n", "        fn=lambda utm1, A: pt.linalg.solve(A, utm1, assume_a=assume_a),\n", "        outputs_info=[u[1:-1]],\n", "        non_sequences=[A],\n", "        n_steps=nt,\n", "        strict=True\n", "\n", "    )\n", "    # Add the fixed edges\n", "    u_history = pt.pad(u_history, [[0, 0], [1, 1]], mode=\"constant\")\n", "\n", "    fn = pytensor.function([nt, alpha], u_history, mode=mode)\n", "    return fn\n" ],
    "id" : "8b38bfaeeb055722",
    "outputs" : [ ],
    "execution_count" : 1
  }, {
    "metadata" : {
      "ExecuteTime" : {
        "end_time" : "2025-03-17T18:08:23.030771Z",
        "start_time" : "2025-03-17T18:08:20.203522Z"
      }
    },
    "cell_type" : "code",
    "source" : "fn = build_heat_equation_fn(mode=\"FAST_RUN\")",
    "id" : "c46e787114deeb3",
    "outputs" : [ ],
    "execution_count" : 2
  }, {
    "metadata" : {
      "ExecuteTime" : {
        "end_time" : "2025-03-17T18:08:23.356448Z",
        "start_time" : "2025-03-17T18:08:23.137105Z"
      }
    },
    "cell_type" : "code",
    "source" : "u_history_res = fn(nt=50, alpha=2.5)",
    "id" : "51625468a0c5107",
    "outputs" : [ ],
    "execution_count" : 3
  }, {
    "metadata" : {
      "ExecuteTime" : {
        "end_time" : "2025-03-17T18:08:23.562571Z",
        "start_time" : "2025-03-17T18:08:23.365995Z"
      }
    },
    "cell_type" : "code",
    "source" : [ "# Create a heatmap of the temperature evolution.\n", "plt.figure(figsize=(8, 5))\n", "# Use imshow: note that we set extent to map array indices to physical time and space.\n", "extent = [0, L, T, 0]  # time goes from 0 to T; flipped so time increases downward\n", "plt.imshow(u_history_res, extent=extent, aspect='auto', cmap='hot')\n", "plt.colorbar(label='Temperature')\n", "plt.xlabel('Position, x')\n", "plt.ylabel('Time, t')\n", "plt.title('Heat Equation: Temperature Evolution')\n", "plt.show()" ],
    "id" : "97373d7c55d6085",
    "outputs" : [ {
      "data" : {
        "text/plain" : [ "<Figure size 800x500 with 2 Axes>" ],
        "image/png" : "iVBORw0KGgoAAAANSUhEUgAAAqUAAAHWCAYAAABOhvDAAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAZHFJREFUeJzt3XtcVNX6P/DPgMwMouAFATUUNRVvSaESVGLFN0pTOZWimaLHsot30lJTwUpRjxodtUg7ajfTTPPr1ziYkZQmHfNWdrzfzQRvAYoKOrN/f/hjcpxB1uxZsIeZz/v14lVunll77csMz6y197N1iqIoICIiIiLSkJfWHSAiIiIiYlJKRERERJpjUkpEREREmmNSSkRERESaY1JKRERERJpjUkpEREREmmNSSkRERESaY1JKRERERJpjUkpEREREmmNSSuSBunXrhm7dumndDSLVUlNTodPppLcbFhaGwYMHS2+XiCrGpJQ0t2zZMuh0Omzfvt3u77t164b27dtXah8yMzORmpoqHN+tWzfodDq7P+Hh4ZXXUQfs3bsXqampOH78uNZdKdfx48fL3Y+3/7jydlS2GTNmYO3atVp3w66wsLByj9njjz+udffs2rp1K1JTU1FQUKB1V4joFjW07gCRK8jMzMTChQsdSkzvuusupKWl2SwPCAiQ2DP19u7di2nTpqFbt24ICwuz+t0333yjTadu06BBA3zyySdWy+bOnYvff/8d77zzjk2sp5oxYwaeeeYZJCQkaN0VuyIiIvDqq6/aLG/UqJEGvanY1q1bMW3aNAwePBh16tSx+t2BAwfg5cXxGiItMCklUikgIADPPfec1t1QRa/Xa90FAICfn5/NPlyxYgX+/PPPartvK2I2m1FaWgqj0eg2/WjcuLHbHC+DwaB1F4g8Fr8OUrX16aefIjIyEr6+vqhXrx769euHU6dOWcVs3rwZffr0QZMmTWAwGBAaGoqxY8fi6tWrlpjBgwdj4cKFAGA19SjLli1b0LlzZxiNRrRo0QIffPCBzfVwZdPYy5Yts3m9TqezGsE9ceIEXnnlFbRu3Rq+vr6oX78++vTpYzW9vWzZMvTp0wcA8PDDD1u2KScnB4D9a0rPnj2LoUOHIjg4GEajER07dsRHH31kFVPWzzlz5mDRokVo0aIFDAYDOnfujJ9//tkq9vr169i/fz/OnDnj+E67TUlJCVJSUnD33XdbjuNrr72GkpISqzidTocRI0Zg1apVaNu2LXx9fREdHY09e/YAAD744APcfffdMBqN6Natm80lAWWXiuzYsQMxMTHw9fVFs2bNkJGR4XSfPvvsM7Rr1w4GgwFZWVkAgDlz5iAmJgb169eHr68vIiMj8eWXX9q8vri4GB999JHlOJZd8zh48GCbUXDA/vWWd+rH6dOn8fe//x3BwcEwGAxo164dlixZcueD4oA5c+ZAp9PhxIkTNr+bOHEi9Ho9/vzzT8uyVatWWd7bgYGBeO6553D69Ok7rkP0PZSamorx48cDAJo1a2ZzeYi9a0qPHj2KPn36oF69eqhZsybuv/9+fP3111YxOTk50Ol0+OKLLzB9+nTcddddMBqNePTRR3H48OEK9hARARwpJRdSWFiI8+fP2yy/fv26zbLp06djypQp6Nu3L55//nmcO3cO8+fPR9euXbFr1y7LlNyqVatw5coVvPzyy6hfvz62bduG+fPn4/fff8eqVasAAC+++CL++OMPbNy40WYq+U5MJpPd/vr6+sLPzw8AsGfPHjz22GNo0KABUlNTcePGDaSkpCA4OFh4Pbf7+eefsXXrVvTr1w933XUXjh8/jvfffx/dunXD3r17UbNmTXTt2hWjRo3CP//5T0yaNAlt2rQBAMt/b3f16lV069YNhw8fxogRI9CsWTOsWrUKgwcPRkFBAUaPHm0Vv3z5cly6dAkvvvgidDodZs+ejaeeegpHjx6Fj48PgJuJTps2bZCUlGQ3URBlNpvRq1cvbNmyBcOGDUObNm2wZ88evPPOOzh48KDNtZabN2/GunXrMHz4cABAWloannzySbz22mt477338Morr+DPP//E7Nmz8fe//x3fffed1ev//PNPdO/eHX379kX//v3xxRdf4OWXX4Zer8ff//53VX367rvv8MUXX2DEiBEIDAy0JJLvvvsuevXqhQEDBqC0tBQrVqxAnz59sH79evTo0QMA8Mknn+D5559Hly5dMGzYMABAixYtVO1Le/3Iz8/H/fffb0laGzRogH//+98YOnQoioqKMGbMmArbvX79ut33gp+fH3x9fdG3b1+89tpr+OKLLywJYZkvvvgCjz32GOrWrQvg5heqIUOGoHPnzkhLS0N+fj7effdd/Pjjj1bvbbWeeuopHDx4EJ9//jneeecdBAYGAij/8pD8/HzExMTgypUrGDVqFOrXr4+PPvoIvXr1wpdffom//e1vVvEzZ86El5cXxo0bh8LCQsyePRsDBgzAf/7zH6f6TeQRFCKNLV26VAFwx5927dpZ4o8fP654e3sr06dPt2pnz549So0aNayWX7lyxWZ9aWlpik6nU06cOGFZNnz4cMWRt0NsbGy5fX3xxRctcQkJCYrRaLRa1969exVvb2+r9R07dkwBoCxdutRmXQCUlJSUO25Tbm6uAkD5+OOPLctWrVqlAFA2bdpkt/+xsbGWf6enpysAlE8//dSyrLS0VImOjlZq1aqlFBUVWfWzfv36ysWLFy2x//u//6sAUP7v//7PZpuSkpJs1n8nPXr0UJo2bWr59yeffKJ4eXkpmzdvtorLyMhQACg//vijZRkAxWAwKMeOHbMs++CDDxQASkhIiGU7FEVRJk6cqACwii07rnPnzrUsKykpUSIiIpSgoCCltLRUVZ+8vLyU//73vzbbevuxLC0tVdq3b6888sgjVsv9/Pzs7sekpCSrfVUmJSXF5nwurx9Dhw5VGjZsqJw/f95qeb9+/ZSAgAC759utmjZtWu57IS0tzRIXHR2tREZGWr1227ZtVudtaWmpEhQUpLRv3165evWqJW79+vUKAGXq1KnlbqMj76F//OMfNsf+1u25dV+PGTNGAWB1rC9duqQ0a9ZMCQsLU0wmk6IoirJp0yYFgNKmTRulpKTEEvvuu+8qAJQ9e/aUsweJqAyn78llLFy4EBs3brT5ueeee6zi1qxZA7PZjL59++L8+fOWn5CQELRs2RKbNm2yxPr6+lr+v7i4GOfPn0dMTAwURcGuXbuc6m9YWJjd/paNLJlMJmzYsAEJCQlo0qSJ5XVt2rRBfHy86vXeuk3Xr1/HhQsXcPfdd6NOnTrYuXOnqjYzMzMREhKC/v37W5b5+Phg1KhRuHz5Mr7//nur+MTERMvIFgA89NBDAG5Oc5YJCwuDoihOjZICN0e727Rpg/DwcKvj/cgjjwCA1fEGgEcffdRqSjsqKgoA8PTTT6N27do2y2/tMwDUqFEDL774ouXfer0eL774Is6ePYsdO3ao6lNsbCzatm1rs223Hss///wThYWFeOihh1Qfx4rc3g9FUbB69Wr07NkTiqJYbUt8fDwKCwuF+hIVFWX3vXDr+ZSYmIgdO3bgyJEjlmUrV66EwWBA7969AQDbt2/H2bNn8corr1hd69qjRw+Eh4fbTJlXhczMTHTp0gUPPvigZVmtWrUwbNgwHD9+HHv37rWKHzJkiNU12/beG0RkH6fvyWV06dIFnTp1sllet25dq6nBQ4cOQVEUtGzZ0m47ZdPHAHDy5ElMnToV69ats7pmDbh5uYAz/Pz8EBcXV+7vz507h6tXr9rtZ+vWrZGZmalqvVevXkVaWhqWLl2K06dPQ1EUy+/UbtOJEyfQsmVLm7uOy6b7b78W8NYkG4AlQb19H8tw6NAh7Nu3r9zp1bNnz96xb2XVEEJDQ+0uv73PjRo1slx+UaZVq1YAbl63eP/99zvcp2bNmtmNW79+Pd5++23s3r3b6lrUyqi/aa8f586dQ0FBARYtWoRFixbZfc3t22JPYGDgHd8LANCnTx8kJydj5cqVmDRpEhRFwapVq/DEE0/A398fwF/nWevWrW1eHx4eji1btlTYF9lOnDhh+QJzq1vfG7eWrKvK9waRu2FSStWO2WyGTqfDv//9b3h7e9v8vlatWgBujlT+z//8Dy5evIjXX38d4eHh8PPzw+nTpzF48GCYzeaq7nq5yktCTCaTzbKRI0di6dKlGDNmDKKjoxEQEACdTod+/fpV2TbZ2+8ArBJkWcxmMzp06IB58+bZ/f3tyWZ5fZPZZ0f7dOuIaJnNmzejV69e6Nq1K9577z00bNgQPj4+WLp0KZYvXy7UD0fOG3v9KDtfnnvuOSQlJdl9ze0zFWo1atQIDz30EL744gtMmjQJP/30E06ePIlZs2ZJad/RfVFZqvK9QeRumJRStdOiRQsoioJmzZpZRrDs2bNnDw4ePIiPPvoIgwYNsizfuHGjTWxljEw1aNAAvr6+OHTokM3vDhw4YPXvstGU24t527tb+csvv0RSUhLmzp1rWXbt2jWb1zqyTU2bNsWvv/4Ks9lsNVq6f/9+y++10qJFC/zyyy949NFHK20E8VZ//PEHiouLrUZLDx48CACWywJk9Gn16tUwGo3YsGGDVRmipUuX2sSWt466devaLQBv77yxp0GDBqhduzZMJlOFI50yJCYm4pVXXsGBAwewcuVK1KxZEz179rT8vuw8O3DggOVSiDIHDhy443noyHvI0ffG7e9XwDXeG0TuhteUUrXz1FNPwdvbG9OmTbMZfVAUBRcuXADw14jFrTGKouDdd9+1abMsAZH5hBdvb2/Ex8dj7dq1OHnypGX5vn37sGHDBqtYf39/BAYG4ocffrBa/t5779lt9/btnj9/vs2IkCPb1L17d+Tl5WHlypWWZTdu3MD8+fNRq1YtxMbGVtjG7WSVhOrbty9Onz6NxYsX2/zu6tWrKC4udqr92924cQMffPCB5d+lpaX44IMP0KBBA0RGRkrrk7e3N3Q6ndVxO378uN0nN/n5+dk9ji1atEBhYSF+/fVXy7IzZ87gq6++qnD9ZX14+umnsXr1avz22282vz937pxQO6KefvppeHt74/PPP8eqVavw5JNPWiX/nTp1QlBQEDIyMqwuZ/j3v/+Nffv2WSoS2OPIe8jR98a2bduQm5trWVZcXIxFixYhLCzM7rXCRKQOR0qp2mnRogXefvttTJw4EcePH0dCQgJq166NY8eO4auvvsKwYcMwbtw4hIeHo0WLFhg3bhxOnz4Nf39/rF692u61XWXJxqhRoxAfHw9vb2/069fvjv0oLCzEp59+avd3ZYXEp02bhqysLDz00EN45ZVXLIleu3btrBIJAHj++ecxc+ZMPP/88+jUqRN++OEHywjdrZ588kl88sknCAgIQNu2bZGbm4tvv/0W9evXt4qLiIiAt7c3Zs2ahcLCQhgMBjzyyCMICgqyaXPYsGH44IMPMHjwYOzYsQNhYWH48ssv8eOPPyI9Pd3qBiFRskpCDRw4EF988QVeeuklbNq0CQ888ABMJhP279+PL774Ahs2bLB7LbJajRo1wqxZs3D8+HG0atUKK1euxO7du7Fo0SLL9coy+tSjRw/MmzcPjz/+OJ599lmcPXsWCxcuxN13321zbkRGRuLbb7/FvHnz0KhRIzRr1gxRUVHo168fXn/9dfztb3/DqFGjcOXKFbz//vto1aqV8M1SM2fOxKZNmxAVFYUXXngBbdu2xcWLF7Fz5058++23uHjxYoVtnD592u57oVatWlZPoQoKCsLDDz+MefPm4dKlS0hMTLSK9/HxwaxZszBkyBDExsaif//+lpJQYWFhGDt27B37IfoeKnu/v/HGG+jXrx98fHzQs2dPm2uJAWDChAn4/PPP8cQTT2DUqFGoV68ePvroIxw7dgyrV6/m05+IZKr6G/6JrJWVhPr555/t/j42NtaqJFSZ1atXKw8++KDi5+en+Pn5KeHh4crw4cOVAwcOWGL27t2rxMXFKbVq1VICAwOVF154Qfnll19sSsfcuHFDGTlypNKgQQNFp9NVWB7qTiWhbn/t999/r0RGRip6vV5p3ry5kpGRYbdkz5UrV5ShQ4cqAQEBSu3atZW+ffsqZ8+etSln8+effypDhgxRAgMDlVq1ainx8fHK/v37bUrZKIqiLF68WGnevLmlBFVZeajbS0IpiqLk5+db2tXr9UqHDh1syuuUld35xz/+YbNPbu+nrJJQinKzVNCsWbOUdu3aKQaDQalbt64SGRmpTJs2TSksLLTqw/Dhw4X6XFbCZ9WqVZZlZefa9u3blejoaMVoNCpNmzZVFixYYNNPZ/pU5l//+pfSsmVLxWAwKOHh4crSpUvtnhv79+9Xunbtqvj6+trs02+++UZp3769otfrldatWyuffvppuSWhyutHfn6+Mnz4cCU0NFTx8fFRQkJClEcffVRZtGiR3fhb3akklL1yVYsXL1YAKLVr17Yq+3SrlStXKvfee69iMBiUevXqKQMGDFB+//13qxhn3kOKoihvvfWW0rhxY8XLy8uqPJS999GRI0eUZ555RqlTp45iNBqVLl26KOvXr7eKsXc+KcqdS1URkTWdovDqa6KqlpqaavfyA9JWt27dcP78ebtT2UREVLk470BEREREmmNSSkRERESaY1JKRERERJpziaR04cKFCAsLg9FoRFRUFLZt23bH+FWrViE8PBxGoxEdOnSweTKOoiiYOnUqGjZsCF9fX8TFxdmtFUmkldTUVF5P6oJycnJ4PSkRkUY0T0pXrlyJ5ORkpKSkYOfOnejYsSPi4+PLfbTd1q1b0b9/fwwdOhS7du1CQkICEhISrP6QzJ49G//85z+RkZGB//znP/Dz80N8fDyuXbtWVZtFRERERA7Q/O77qKgodO7cGQsWLABw87F3oaGhGDlyJCZMmGATn5iYiOLiYqxfv96y7P7770dERAQyMjKgKAoaNWqEV199FePGjQNws55kcHAwli1bVmHtSSIiIiKqepoWzy8tLcWOHTswceJEyzIvLy/ExcVZPT3jVrm5uUhOTrZaVvbUHAA4duwY8vLyrB6ZFxAQgKioKOTm5tpNSktKSqyeHmI2m3Hx4kXUr1+/Sh5rSERERM5TFAWXLl1Co0aNNHmwwbVr11BaWiqlLb1eD6PRKKWt6kLTpPT8+fMwmUwIDg62Wh4cHGx5rvDt8vLy7Mbn5eVZfl+2rLyY26WlpWHatGmqtoGIiIhcy6lTp3DXXXdV6TqvXbuGZs2alZtrOCokJATHjh3zqMSUjxkFMHHiRKvR18LCQjRp0gRGABWNk4YLrqO1YJxoe7LX21wvGNhKpDHBtlpIjgsTjAsVjBP9PKsh2sGWgnGiO7CJYFwjwTjbx4/aV7/iEABAgGCc7aMd7RP9YJb5sXZDME70evViwbhCwbgLgnH2r9G39Ydg3EnBuKOCcYI3ot44Ihb3u+BqTwnGHReME+yeUJzorrN9iqr95gQH7w4Irtb+sJH6OFnrVXDz3ajm0cjOKi0tRV5eHk6dOgZ/f3+n2ioqKkJoaDOUlpYyKa0qgYGB8Pb2Rn5+vtXy/Px8hISE2H1NSEjIHePL/pufn4+GDRtaxURERNht02AwwGAw2CzXoeKk1LuC35cRzflET72agnGib0t/0asURDbYR7At211un69gnGhuI7xTBONqiJ4FsneM7LNFdAfWEowT3dGi7blyUir6DhedTjQJxokmw6LngOg+Fj1HRc95wfeQ6OEQPfVET3nRzyCZu0X0Y0Xws1t0l8g+U0TfGZI3V9NL7/z9/Z1OSj2Vpnff6/V6REZGIjs727LMbDYjOzsb0dHRdl8THR1tFQ8AGzdutMQ3a9YMISEhVjFFRUX4z3/+U26bRERERHLckPTjeTSfvk9OTkZSUhI6deqELl26ID09HcXFxRgyZAgAYNCgQWjcuDHS0tIAAKNHj0ZsbCzmzp2LHj16YMWKFdi+fTsWLVoE4Oa3ozFjxuDtt99Gy5Yt0axZM0yZMgWNGjVCQkKCVptJREREHkFGUsmkVBOJiYk4d+4cpk6diry8PERERCArK8tyo9LJkyet7qCLiYnB8uXLMXnyZEyaNAktW7bE2rVr0b59e0vMa6+9huLiYgwbNgwFBQV48MEHkZWV5VHXZRARERFVJ5rXKXVFRUVFCAgIgC8qvn6lrWCbrh7XQvRaKJE7rO4WbEv0vh/Z9wc1FYwTvSGqhsjdX4D4bWeiO1B0Q0Tv2AquOAQAECgYV0cwzh2uKRW9tvOyYFyBYNx5wbj8ikMAiN8hdEIw7rBgnOBtLjcE7+oRvYFJdDMk368lFCe66wTvJDpSUnEMAOwVXK2rxikAruLmDctVfV1nWe5QWHhCyo1OAQFNNdkOLWk+UkpERETkPkxwfvpd9GZH98Kk1Emip41o3HXBONFxmSuCcaWC36L1BQJBotVsRNoCgHOCcaK3l8q+bbSR4OiN8M2grj7qJjpSKjoCypJQtkSPregxE62beFow7rhgnOBQpCL4HhId8BWNE/1sEY0rEIwT+YwUbEv0s1v0b4HomSz6t0r230hyb0xKiYiIiKThjU5qMSklIiIikoZJqVqa1iklIiIiIgI4UkpEREQkEUdK1WJSSkRERCSNCc7fuuWZt35x+p6IiIiINMeRUieZBeNEy2fIjpNdqEZ/SSCoSLCxi4JxdQTjRMu2iFYWEo0T/WoXIlo6SvSoCdaCEa5Tc0Ewrr5gnGhJKBbPVx8nesxES0KJ1lI6LhamnBSLE+2e7BJOsks9iX6miXxGinzWQvyzW3apJ9lxon9LqwfWKVWLSSkRERGRNLymVC1O3xMRERGR5jhSSkRERCQNR0rVYlJKREREJA2TUrU4fU9EREREmuNIKREREZE0vPteLSalRERERNJw+l4tJqVVRLQGW6lgnOz6o1cE4+qKNFgo2FiBYJxo7T8/wbiagnGy65SKaiBY27GG6FlQIDmujuQ42YVjXblOqVbHTLT+qGDB0BtnxeK0qisqWr5V9LNFNK5AME7kM1Lww1v0s1t2PVPRv1XuVX+UKhuTUiIiIiJpOFKqFpNSIiIiImmYlKrFu++JiIiISHMcKSUiIiKShiOlajEpJSIiIpKGJaHU4vQ9EREREWmOI6VERERE0nD6Xi0mpU4SHWAXjbsuGCe7KqJoDbvSkopj9JcEGxONKxCME60/ahCM0wvGyZ5vEC3sV0+wVqTvZcEGCwTj6gjG1RKMEy0w6w51SkXfabKP2XmxsKuCVS9F63aK1gt19XqmBYJxEj/7RD5rAfn1R0XjRP9Wyf4bWT0wKVWL0/dEREREpDmOlBIRERFJw5FStZiUEhEREUnDpFQtTt8TERERkeY4UkpEREQkDeuUqsWklIiIiEgaE5xPKpmUkgqi1XtEy2fIjhMt7yFYCEao/Ii+SLCxAsG42oJxoiWhZFcg0qokVKlgXB3Bo+t3UixOL1heSLgklOiOdoeSUKJxgiWhSgWPrWjdoALBuELBONkloUTbE+1fgUZxAp+RoodM9LNbdqkn2XGiH3vk3piUEhEREUnDG53UYlJKREREJA2TUrV49z0RERERaY4jpURERETS8O57tZiUEhEREUnD6Xu1OH1PRERERJrjSCkRERGRNBwpVYtJaRURvTpEdv3REsE4mTXx6ooWzrskGFcgGCe7/qiPYJy3YJwo2XVKRU8Wf8E4g2jdU8E42ftZJ3ECSBE8GLLf4KJvSNE3uGjtYNH3pOw6pRclx2lVz1R0/wm8NWTWjgbETxXZ9Uw988pIJqVqcfqeiIiIiDTHkVIiIiIiaThSqhaTUiIiIiJpWBJKLU7fExEREZHmOFJKREREJM0NOH8HLKfviYiIiMgpTErV4vQ9EREREWmOI6VOEr0UWTROtPSkK9czLRYs7egnWtNPdj1T0TqlesE4UaIngVZ1SkULI9YUjBM9bqL7WfQrtLfoDhTgLsdMtJilu9Qp1ar+qGCcyGekq9cfFT2VZf+NrB44UqoWk1IiIiIiaXj3vVqcviciIiIizTEpJSIiIpLmhqQfxy1cuBBhYWEwGo2IiorCtm3b7hifnp6O1q1bw9fXF6GhoRg7diyuXRO92EM+Tt8TERERSXMDzo/5OZ6Urly5EsnJycjIyEBUVBTS09MRHx+PAwcOICgoyCZ++fLlmDBhApYsWYKYmBgcPHgQgwcPhk6nw7x585zsvzocKSUiIiKq5ubNm4cXXngBQ4YMQdu2bZGRkYGaNWtiyZIlduO3bt2KBx54AM8++yzCwsLw2GOPoX///hWOrlYmJqVERERE0sibvi8qKrL6KSmxX2+htLQUO3bsQFxcnGWZl5cX4uLikJuba/c1MTEx2LFjhyUJPXr0KDIzM9G9e3fnNt8JnL4nIiIiksYE5++ev/n60NBQq6UpKSlITU21iT5//jxMJhOCg4OtlgcHB2P//v121/Dss8/i/PnzePDBB6EoCm7cuIGXXnoJkyZNcrLv6jEprSKiZQxF42TXKZVZFlG0vp6faI1F0RqBovVHReNE5xG0mm9w9TqlovVHZR8PZ8sD3kp2nVLRYyF6bLWqUyoaJ1pXtEAwTrTuqWh7onGi2yt4PEQOh+ihlf0ZL/q3RfbfNLLv1KlT8Pf3t/zbYDBIazsnJwczZszAe++9h6ioKBw+fBijR4/GW2+9hSlTpkhbjyOYlBIRERFJI69Oqb+/v1VSWp7AwEB4e3sjPz/fanl+fj5CQkLsvmbKlCkYOHAgnn/+eQBAhw4dUFxcjGHDhuGNN96Al1fVj7jwmlIiIiIiaaq+JJRer0dkZCSys7Mty8xmM7KzsxEdHW33NVeuXLFJPL29b041KYri0Ppl4UgpERERUTWXnJyMpKQkdOrUCV26dEF6ejqKi4sxZMgQAMCgQYPQuHFjpKWlAQB69uyJefPm4d5777VM30+ZMgU9e/a0JKdVjUkpERERkTQ3AOgktOGYxMREnDt3DlOnTkVeXh4iIiKQlZVlufnp5MmTViOjkydPhk6nw+TJk3H69Gk0aNAAPXv2xPTp053su3o6RasxWhdWVFSEgIAA+KLi0yq0gt9rHddIg/ZE1xnkIxgYXHEIAKCBYJz9y2ts1ZO83jqS11tbcpyfYBxvdLLFG53sk32j0znJ682TvN78ikMA4KzA3UR/CK5SNO6URu1VdZwC4CqAwsJCoWsxZSrLHQoLn4C/v+gfuPLauo6AgH9rsh1a4jWlRERERKQ5Tt87SXQgRTROdIBENE6L0lGigzdXBTvn6+olobS59Eb+wdWqJJTogIJzAw+VS/RYiMZpNVIqu0xbgWCc7BFV0TjJI8Oin2kiu1mrUk+y/wbJ/htZPWgzfe8OmJQSERERSWOC80mpe6Xpojh9T0RERESa40gpERERkTQypt49c/reJUZKFy5ciLCwMBiNRkRFRWHbtm13jF+1ahXCw8NhNBrRoUMHZGZmWn53/fp1vP766+jQoQP8/PzQqFEjDBo0CH/8IXpPIREREZFaVV88311onpSuXLkSycnJSElJwc6dO9GxY0fEx8fj7NmzduO3bt2K/v37Y+jQodi1axcSEhKQkJCA3377DcDNJxTs3LkTU6ZMwc6dO7FmzRocOHAAvXr1qsrNIiIiIiIHaF6nNCoqCp07d8aCBQsA3HwsVmhoKEaOHIkJEybYxCcmJqK4uBjr16+3LLv//vsRERGBjIwMu+v4+eef0aVLF5w4cQJNmjSpsE+O1CnVog6oI+01lrxekTjRMqCi5Ud9AwQDReuFisaJ1gutL7m9OoJxovVHRe+W16pOKe++t8W77+27ILk90fqjgnFXC8XiRMqZipZQlV3P9LRgnFZ1TytqzzXqlHaBv79zV0cWFd1AQMA21imtSqWlpdixYwfi4uIsy7y8vBAXF4fc3Fy7r8nNzbWKB4D4+Phy44GbJ6dOp0OdOnXs/r6kpARFRUVWP0RERESOM8H5qXvPvPte0xudzp8/D5PJZHkEVpng4GDs37/f7mvy8vLsxufl2f9eee3aNbz++uvo379/ud820tLSMG3aNBVbIE70gS+icVrUHwXEBlJkD8poVqdUdGROdp1SrU6CEsE40QMsup9FR1Td4YlOoiOgWh0zrUZKReNkP3FKcpzMgWvRtrSqZyr7bxoR4ALXlFam69evo2/fvlAUBe+//365cRMnTkRhYaHl59Qp0YkEIiIiolvxRie1NB0pDQwMhLe3N/Lzra+wyc/PR0iI/SsTQ0JChOLLEtITJ07gu+++u+M1GQaDAQaDQeVWEBEREZW5gZtXtzrDM6fvNR0p1ev1iIyMRHZ2tmWZ2WxGdnY2oqOj7b4mOjraKh4ANm7caBVflpAeOnQI3377LerXF70ThYiIiIi0oHnx/OTkZCQlJaFTp07o0qUL0tPTUVxcjCFDhgAABg0ahMaNGyMtLQ0AMHr0aMTGxmLu3Lno0aMHVqxYge3bt2PRokUAbiakzzzzDHbu3In169fDZDJZrjetV68e9HrRC9iIiIiIHMWRUrU0T0oTExNx7tw5TJ06FXl5eYiIiEBWVpblZqaTJ0/Cy+uvAd2YmBgsX74ckydPxqRJk9CyZUusXbsW7du3BwCcPn0a69atAwBERERYrWvTpk3o1q1blWwXEREReSImpWppXqfUFVVGnVIt6oUC4jVDZbYnWn9UtG/1RS8y0apOaR2N4rSqUyp6+TXvvrfFu++rJk6jOqUXBM8DkRqkIrVMRdsCxOuFunp91OpRp7QF/P2d+yAqKjIhIOCIx9Up1XyklIiIiMh9mOD8SKlnFtNiUuok0dNGdMBFqwfDyKx1J9qW6OBNTcGdLFzPVKv6o6JETxYtDi4gvx6s6H52hyc6iR5bLQoMA+JvStl1QAsE4wSfmCS77ulVwc8g0d0n83NUNE7040L23yDZkw3VA5NStdy6TikRERERVQ8cKSUiIiKS5gacH/PzzJFSJqVERERE0jApVYvT90RERESkOY6UEhEREUnDkVK1mJQSERERSWOC80mlZ5aQZ1JaRUTLYmhVOkpm+RHR6jOy63QbBRvUiZYWEo2TfRGM7Boqsk8C2SWhRPefKz8hWLTejugxk10PSHbpqCLBONmloySvVxHcXtHuyfzsk31otSr15JnPJSK1mJQSERERSXMDFT8PsiIcKSUiIiIipzApVYt33xMRERGR5jhSSkRERCQNR0rVYlJKREREJItidj6n9MyclNP3RERERKQ9jpQSERERyWKG82VKPbN2PpNSZ8kuKSnanmhZRC3qmYrWzSsWjBMtdylaS9BftJigaJ1Sb8E42URPKtGTRXb9UdH9J1p/1JXndWQfCy0KDAPihTZF37yy65RKLhgqu/6o6G6R+Tkq+zNe9BTV6m9ftWCC8xvkVjtEnCt/zBMRERGRh+BIKREREZEsHClVjUkpERERkSy8plQ1Tt8TERERkeY4UkpEREQkC6fvVWNSSkRERCQLp+9V4/Q9EREREWmOI6VVRHQkXrSWnOiXKC3KHcounSi7xKKP4E7xFS1iKLtOqeyDK/skEK1TKrv+qOh+lvlVW3aRRdn1TEWPWYlgnGihTRevZ3pV8JwX3QwtPtNkr1Orvy0eOQtthvMbzpFSIiIiIiJtcKSUiIiISBbe6KQak1IiIiIiWXijk2qcviciIiIizXGklIiIiEgWTt+rxqSUiIiISBYmpaoxKXWS6HkjO060YoxonMyqQbJLOIlWIPKRHKcXLKPjLVrORpTsk0X04IruaNHyQqI7WnZJKJklumQfC1eu5eZInOzSUYLvIZPguSe7YpXszZX5OSp6qsj+m6HV3z5yb0xKiYiIiGThjU6qMSklIiIikoXT96rx7nsiIiIi0hxHSomIiIhkUeD89LsioyPVD5NSIiIiIlk4fa8ap++JiIiISHMcKSUiIiKShSOlqjEpdTGil6HILncoM052SUTR2n96wTjZ5TP9BTsoXD5Tq5qXsmtjyi4cK0qLOqWitHhDAi5fp9QkGFckuFrZdUVl1h8VjdPqVJH9N8gjsSSUapy+JyIiIiLNcaSUiIiISBZO36vGpJSIiIhIFialqnH6noiIiIg0x5FSIiIiIll4o5NqTEqJiIiIZDHD+el3D01KOX1PRERERJrjSGkVEf3SI/rlSnZtulKJ7ckuiSg7TrQspmi5S9E4P9F6plrVKRU9CWQXhBXdga78FVqrN7joMSsRjJNcp9QkuF53qCsqO06Lz25H4mR/TLkVTt+rxqSUiIiISBbefa+aK489EBEREZGH4EgpERERkSwcKVWNSSkRERGRLLymVDVO3xMRERGR5jhSSkRERCQLp+9V40gpERERkSwmST8qLFy4EGFhYTAajYiKisK2bdvuGF9QUIDhw4ejYcOGMBgMaNWqFTIzM9WtXAKOlDpJdnlCV69nKlJfzyixLUfiRGsdin4TE40TLbMpyk+wtqO36MkiWshQ9MCJ1ikV3YGy65lqQfYbUvaxFX0TCZ57JsHtkF1/VLQ92XGyP6tEdrNoW65ef5T1TKvOypUrkZycjIyMDERFRSE9PR3x8fE4cOAAgoKCbOJLS0vxP//zPwgKCsKXX36Jxo0b48SJE6hTp07Vd/7/Y1JKREREJIsC57NsxfGXzJs3Dy+88AKGDBkCAMjIyMDXX3+NJUuWYMKECTbxS5YswcWLF7F161b4+NwcHQgLC3Om107j9D0RERGRLBKn74uKiqx+Skrsj7OXlpZix44diIuLsyzz8vJCXFwccnNz7b5m3bp1iI6OxvDhwxEcHIz27dtjxowZMJm0u6CVSSkRERGRCwoNDUVAQIDlJy0tzW7c+fPnYTKZEBwcbLU8ODgYeXl5dl9z9OhRfPnllzCZTMjMzMSUKVMwd+5cvP3229K3QxSn74mIiIhkkVin9NSpU/D397csNhgMTjZ8yyrMZgQFBWHRokXw9vZGZGQkTp8+jX/84x9ISUmRth5HMCklIiIikkViSSh/f3+rpLQ8gYGB8Pb2Rn5+vtXy/Px8hISE2H1Nw4YN4ePjA2/vv+4ibdOmDfLy8lBaWgq9XvSuVnk4fU9ERERUjen1ekRGRiI7O9uyzGw2Izs7G9HR0XZf88ADD+Dw4cMwm/8a1j148CAaNmyoSUIKcKTU5Yh+uZIdJ7OsiOwqNaJxopWFZFcgkl2pSPSYGQUPmmicTvQkEN2BsktHuUNJKNmlngSPmSJY6kmL0keA/NJRsrdDdpzI4dWq1JPsOI+kUfH85ORkJCUloVOnTujSpQvS09NRXFxsuRt/0KBBaNy4seW61JdffhkLFizA6NGjMXLkSBw6dAgzZszAqFGjnOy8ekxKiYiIiGSReE2pIxITE3Hu3DlMnToVeXl5iIiIQFZWluXmp5MnT8LL669RgNDQUGzYsAFjx47FPffcg8aNG2P06NF4/fXXney8ekxKiYiIiNzAiBEjMGLECLu/y8nJsVkWHR2Nn376qZJ7JY5JKREREZEsGk3fuwOXuNHJ0We1rlq1CuHh4TAajejQocMdn9P60ksvQafTIT09XXKviYiIiG5jhvOF86vJc1dv3LiBb7/9Fh988AEuXboEAPjjjz9w+fJlVe1pnpSWPas1JSUFO3fuRMeOHREfH4+zZ8/ajd+6dSv69++PoUOHYteuXUhISEBCQgJ+++03m9ivvvoKP/30Exo1alTZm0FERETkMU6cOIEOHTqgd+/eGD58OM6dOwcAmDVrFsaNG6eqTc2T0luf1dq2bVtkZGSgZs2aWLJkid34d999F48//jjGjx+PNm3a4K233sJ9992HBQsWWMWdPn0aI0eOxGeffWZ5pisRERFRpTJL+nFxo0ePRqdOnfDnn3/C19fXsvxvf/ubVWkqR2ialKp5Vmtubq5VPADEx8dbxZvNZgwcOBDjx49Hu3btKuxHSUmJzfNliYiIiBwm47n31eCa0s2bN2Py5Mk2NU3DwsJw+vRpVW1qeqPTnZ7Vun//fruvycvLq/DZrrNmzUKNGjWEa22lpaVh2rRpDvbeMaJferSoK+pInEh9PdFafUbBONl1SkXLZ7p6nVLRc0o0zkewqKRonE70gLBOqQ1F8A0p+/2tVd1O2fVHReujyt4OmftZcqlazeqeVoMBP1LJbDbDZLI9E37//XfUrl1bVZuaT9/LtmPHDrz77rtYtmwZdDqd0GsmTpyIwsJCy8+pU6cquZdERETkljxk+v6xxx6zuolcp9Ph8uXLSElJQffu3VW1qelIqZpntYaEhNwxfvPmzTh79iyaNGli+b3JZMKrr76K9PR0HD9+3KZNg8EAg8Hg5NYQERGRx/OQklBz5szB448/jrZt2+LatWt49tlncejQIQQGBuLzzz9X1aamI6VqntUaHR1tcwHtxo0bLfEDBw7Er7/+it27d1t+GjVqhPHjx2PDhg2VtzFEREREHiI0NBS//PIL3njjDYwdOxb33nsvZs6ciV27diEoKEhVm5oXz3f0Wa2jR49GbGws5s6dix49emDFihXYvn07Fi1aBACoX78+6tevb7UOHx8fhISEoHXr1lW7cURERORZPGCk9Pr16wgPD8f69esxYMAADBgwQEq7mieljj6rNSYmBsuXL8fkyZMxadIktGzZEmvXrkX79u212gQiIiKim2RcE+ri15T6+Pjg2jXRWwDF6RRFUaS3Ws0VFRUhICAAvgAqulWqjmCbonEBktuTHSfSP9F77kTXKdqeaJy/YFxNwTg/ye2JViWQHSd6s7xoHO++t4N339slevd9seT2RIv/XZIcVyCxrUKJ66yMuKrunwLgKoDCwkL4+4t+2stRljsUzgb8fSuOv2NbV4GA17TZDlEzZszAwYMH8eGHH6JGDTljnJqPlBIRERG5jbLHjDrbhov7+eefkZ2djW+++QYdOnSAn5/1EM2aNWscbpNJqZNEzzvZNd1kr1dmTTzZozKy64+KjrZoRfaxFY2TvZ+9BE8E2fVgZd69qdX7UXS9smtZyh4BFa0XKnsEVLQ9V65nKrN2NOA+f6uqBQ+YvgeAOnXq4Omnn5baJpNSIiIiInLI0qVLpbfJpJSIiIhIFg+4+76yMCklIiIiksVDktJmzZrd8cmZR48edbhNJqVERERE5JAxY8ZY/fv69evYtWsXsrKyMH78eFVtMiklIiIiksVDbnQaPXq03eULFy7E9u3bVbWp6WNGiYiIiNyKSdJPNfXEE09g9erVql7LkVIXI7sch2hZEZntyS7oLbucjexa7bK/2ck+B0TjREtCiR4P2fvPg2rnC7fnLsXzXT1Oi880LT67HWmvGudMVMm+/PJL1KtXT9VrHU5Kmzdvjp9//tnm+fIFBQW47777VF3YSkREROQWPORGp3vvvdfqRidFUZCXl4dz587hvffeU9Wmw0np8ePHYTLZ7q2SkhKcPn1aVSeIiIiI3IIC568JrQYPgO/du7dVUurl5YUGDRqgW7duCA8PV9WmcFK6bt06y/9v2LABAQF/PQXdZDIhOzsbYWFhqjpBRERERNVHamqq9DaFk9KEhAQAgE6nQ1JSktXvfHx8EBYWhrlz50rtHBEREVG14iHT997e3jhz5gyCgoKsll+4cAFBQUF2Z9UrIpyUms03x6KbNWuGn3/+GYGBgQ6vjIiIiMiteUhJKEWxf41BSUkJ9Hq9qjYdvqb02LFjqlZERERERNXbP//5TwA3Z84//PBD1KpVy/I7k8mEH374ofKvKSUiIiKiCrj59P0777wD4OZIaUZGBry9/yrYp9frERYWhoyMDFVtMymtIrLPL9H2RONKBeO0qFMqWktQq7qYrv4ECtFzQLROqeikjOz97A51SmXXM5X5vnUkztXrirpDPVPRtkTPAdnnqCgXzq0qj5snpWUz5g8//DDWrFmDunXrSmubSSkREREROWTTpk3S22RSSkRERCSLh9zoBAC///471q1bh5MnT6K01Hrcft68eQ63x6SUiIiISBY3n74vk52djV69eqF58+bYv38/2rdvj+PHj0NRFNx3332q2pR6OZyXlxceeeQR7NixQ2azRERERORCJk6ciHHjxmHPnj0wGo1YvXo1Tp06hdjYWPTp00dVm1KT0iVLlqBr164YPny4zGaJiIiIqgcz/hotVftTDabv9+3bh0GDBgEAatSogatXr6JWrVp48803MWvWLFVtSp2+Hzx4MIDKefQUERERkcvzkGtK/fz8LNeRNmzYEEeOHEG7du0AAOfPn1fVpuqk9PDhwzhy5Ai6du0KX19fKIoCnU6ntjkiIiIiqibuv/9+bNmyBW3atEH37t3x6quvYs+ePVizZg3uv/9+VW06nJReuHABiYmJ+O6776DT6XDo0CE0b94cQ4cORd26dTF37lxVHaGbZNcxlB0nUhNPtG6eaO0/0fqZWtVilF2nVPRYiJ4ronVFZe9nrerGakH2+1a0PdnnvOh7t0Rye65cL9SRONH+yfwc1epvQTW4D0c7HnKj07x583D58mUAwLRp03D58mWsXLkSLVu2VHXnPaDi7+nYsWNRo0YNnDx5EjVr1rQsT0xMRFZWlqpOEBEREbkFs6QfF2YymfD777+jSZMmAG5O5WdkZODXX3/F6tWr0bRpU1XtOjxS+s0332DDhg246667rJa3bNkSJ06cUNUJIiIiIqoevL298dhjj2Hfvn2oU6eOtHYdHiktLi62GiEtc/HiRRgMBimdIiIiIqqWnL3zXsb0fxVo3749jh49KrVNh5PShx56CB9//LHl3zqdDmazGbNnz8bDDz8stXNERERE1YqHJKVvv/02xo0bh/Xr1+PMmTMoKiqy+lHD4en72bNn49FHH8X27dtRWlqK1157Df/9739x8eJF/Pjjj6o6QURERETVR/fu3QEAvXr1sqq+VFaNyWRyPLN2OClt3749Dh48iAULFqB27dq4fPkynnrqKQwfPhwNGzZ0uANEREREbsND6pRu2rRJepuq6pQGBATgjTfekN2XasnVy3HILi0j0j+tyqzILkEkSnapItmfRVqVIZJd6skdSkK5S+komaWPAPESU7JLR7lyKSrRc0B0nVqdo7LjqoWyJzo524aLi42Nld6mqqT02rVr+PXXX3H27FmYzdZ7rlevXlI6RkRERESua/Pmzfjggw9w9OhRrFq1Co0bN8Ynn3yCZs2a4cEHH3S4PYeT0qysLAwaNMjuI6TUXkNARERE5BZMcH5KrhqkUqtXr8bAgQMxYMAA7Ny5EyUlN+c8CgsLMWPGDGRmZjrcpsO7beTIkejTpw/OnDkDs9ls9cOElIiIiDyaBxTPB27efZ+RkYHFixfDx+evZwI+8MAD2Llzp6o2HU5K8/PzkZycjODgYFUrJCIiIqLq7cCBA+jatavN8oCAABQUFKhq0+Gk9JlnnkFOTo6qlRERERG5NQ+pUxoSEoLDhw/bLN+yZQuaN2+uqk2HryldsGAB+vTpg82bN6NDhw5WQ7YAMGrUKFUdISIiIqr2PKQk1AsvvIDRo0djyZIl0Ol0+OOPP5Cbm4tx48ZhypQpqtp0OCn9/PPP8c0338BoNCInJ8eqYKpOp2NSSkREROTmJkyYALPZjEcffRRXrlxB165dYTAYMG7cOIwcOVJVmzpFURRHXhASEoJRo0ZhwoQJ8PKSXfHRNRQVFSEgIAC+AHQVxPoJtllbclxNye2JbofIemW2VRntGSWvV7Q90TifikMqJU4vuT3ZdUpFyWxP9gya7BqQojUqReuFyq5nKrvuqWhcseT2rkher0h7MttypL1Lktcr2p5oXEXboQC4ipt3gPv7+wu2KkdZ7lAYD/iLflCW19Z1IGCDNtvhqNLSUhw+fBiXL19G27ZtUatWLdVtOTxSWlpaisTERLdNSImIiIhU85CSUGX0ej1q166N2rVrO5WQAip2W1JSElauXOnUSomIiIio+rpx4wamTJmCgIAAhIWFISwsDAEBAZg8eTKuXxedJ7Hm8EipyWTC7NmzsWHDBtxzzz02NzrNmzdPVUeIiIiIqj0Fzt+o5NCFldoYOXIk1qxZg9mzZyM6OhoAkJubi9TUVFy4cAHvv/++w206nJTu2bMH9957LwDgt99+s/rdrTc9EREREXkcEyq+IUWkDRe3fPlyrFixAk888YRl2T333IPQ0FD079+/apLSTZs2ObwSIiIiInIfBoMBYWFhNsubNWsGvV701llrvFuJiIiISBYPKZ4/YsQIvPXWW5Zn3gNASUkJpk+fjhEjRqhqU2ik9KmnnsKyZcvg7++Pp5566o6xa9asUdURIiIiomrPQ4rn79q1C9nZ2bjrrrvQsWNHAMAvv/yC0tJSPProo1b5omhuKJSUBgQEWK4XDQgIcLTfBPEvPVrVMZS5Xq22QTROtI6laA1D2WSfK7KPh2jNS3WTN+Vz5Wkd2X8/RPex7GPmLnVKtdoO0eMh8p6U2RagzTY4EkfVT506dfD0009bLQsNDXWqTaGkdOnSpXjzzTcxbtw4LF261KkVEhEREbktD7nRqTLyQeHBh2nTpuHy5cvSO0BERETkNsySfjyQ8N33Dj6NlIiIiIjc1IULFzB16lRs2rQJZ8+ehdlsnUlfvHjR4TYdKgnFOqREREREd+Ah0/cDBw7E4cOHMXToUAQHB0vJER1KSlu1alXhStVkxkRERERuwQznk8pqMH2/efNmbNmyxXLnvQwOJaXTpk3j3fdEREREHi48PBxXr16V2qZDSWm/fv0QFBQktQNEREREbsMM56fvq8FI6XvvvYcJEyZg6tSpaN++PXx8fKx+7+/v73CbwkkpryetGrJrxGlRC1R27T/REhE+FYcAEK9TKjtOlFb1R0XjRI+H7Lqx7lCn1NXrFbt6PVOt4rTYL65cQ9WR9jySjOtBq8E1pXXq1EFRUREeeeQRq+WKokCn08FkcnwjePc9ERERETlkwIAB8PHxwfLly6v+Rqfbb/UnIiIiott4yEjpb7/9hl27dqF169bS2nTlGTEiIiKi6kXD4vkLFy5EWFgYjEYjoqKisG3bNqHXrVixAjqdDgkJCcLr6tSpE06dOqWuo+Vw6EYnIiIiInI9K1euRHJyMjIyMhAVFYX09HTEx8fjwIEDd7xJ/fjx4xg3bhweeughh9Y3cuRIjB49GuPHj0eHDh1sbnS65557HN4GncKLRW0UFRUhICAAvqj4Bjo/wTZrCsa5Q3uy+2aQ3J5RozjRG7Fkx+kltyc6vSL7RjFXntbhjU5VE3dNo7hiwbgSie1dkdiWJ7WnALgKoLCwUNXd384oyx0K7wb8nbwDtsgEBBx2bDuioqLQuXNnLFiwAMDNyy5DQ0MxcuRITJgwwe5rTCYTunbtir///e/YvHkzCgoKsHbtWqH1eXnZfirrdLqqudGJiIiIiCogsSRUUVGR1WKDwQCDwXaoprS0FDt27MDEiRMty7y8vBAXF4fc3NxyV/Pmm28iKCgIQ4cOxebNmx3q4rFjxxyKF8Gk1M3JHpkRiRNtS3RURqsSRLJpdd267FE82SObHCmt/PZE3xuy2xMdidSqPa1KR8n8HJUdR64lNDTU6t8pKSlITU21iTt//jxMJhOCg4OtlgcHB2P//v12296yZQv+9a9/Yffu3ar61rRpU1WvuxMmpURERESyyChW9P/bOHXqlNX0vb1RUjUuXbqEgQMHYvHixQgMDFTdzieffIKMjAwcO3YMubm5aNq0KdLT09GsWTP07t3b4fZcefCBiIiIqHoxSfrBzaci3fpTXlIaGBgIb29v5OfnWy3Pz89HSEiITfyRI0dw/Phx9OzZEzVq1ECNGjXw8ccfY926dahRowaOHDlS4Wa+//77SE5ORvfu3VFQUGC5hrROnTpIT0+v8PX2MCklIiIiqsb0ej0iIyORnZ1tWWY2m5GdnY3o6Gib+PDwcOzZswe7d++2/PTq1QsPP/wwdu/ebXPZgD3z58/H4sWL8cYbb8Db+68LsTp16oQ9e/ao2g7Nk1JHa2qtWrUK4eHhMBqN6NChAzIzM21i9u3bh169eiEgIAB+fn7o3LkzTp48WVmbQERERHSTRnVKk5OTsXjxYnz00UfYt28fXn75ZRQXF2PIkCEAgEGDBlluhDIajWjfvr3VT506dVC7dm20b98een3FdVuOHTuGe++912a5wWBAcbFo3QVrmialZTW1UlJSsHPnTnTs2BHx8fE4e/as3fitW7eif//+GDp0KHbt2oWEhAQkJCTgt99+s8QcOXIEDz74IMLDw5GTk4Nff/0VU6ZMgdEoWriHiIiISCWJ0/eOSExMxJw5czB16lRERERg9+7dyMrKstz8dPLkSZw5c8a5bbtFs2bN7N4klZWVhTZt2qhqU9M6pY7W1EpMTERxcTHWr19vWXb//fcjIiICGRkZAIB+/frBx8cHn3zyiep+uVOdUtFUXIs6pbLrgMquUypat9PV2xO9u112e7z73hbvvq+e7YmO+cisj6pVHVDRbWCdUluWOqUhgL+TH1hFZiAgT5vtqMibb76JcePGYfny5UhNTcXcuXMxdOhQfPjhhzhy5AjS0tLw4Ycfol+/fg63rdnnfFlNrbi4uL86U0FNrdzcXKt4AIiPj7fEm81mfP3112jVqhXi4+MRFBSEqKioCgvBlpSUoKioyOqHiIiIyGEajZRWlWnTpuHy5ct4/vnnMWvWLEyePBlXrlzBs88+i/fffx/vvvuuqoQU0LAklJqaWnl5eXbj8/LyAABnz57F5cuXMXPmTLz99tuYNWsWsrKy8NRTT2HTpk2IjY21225aWhqmTZumaju0qiUn+3yVOTIju16o6BOJZNc9FaXVNzvRc0B0BNTV65RqVV9WhFZPapK9Xlev7ym7PdE42U/OEomTfa6Icpe/aZqSWBLKFd06wT5gwAAMGDAAV65cweXLl+/4OFMRblWn1Gy+eRR79+6NsWPHAgAiIiKwdetWZGRklJuUTpw4EcnJyZZ/FxUVCd15RkRERORpdDrrixtr1qyJmjVFL94rn2ZJqaM1tQAgJCTkjvGBgYGoUaMG2rZtaxXTpk0bbNmypdy+lPfYLiIiIiKHmHHz4lZnaHa3j5hWrVrZJKa3u3jxosPtapaU3lpTKyEhAcBfNbVGjBhh9zXR0dHIzs7GmDFjLMs2btxoqcGl1+vRuXNnHDhwwOp1Bw8erJTHYRERERFZMaPiu6Qr4uJJ6bRp0xAQECC9XU2n75OTk5GUlIROnTqhS5cuSE9Pt6mp1bhxY6SlpQEARo8ejdjYWMydOxc9evTAihUrsH37dixatMjS5vjx45GYmIiuXbvi4YcfRlZWFv7v//4POTk5WmwiERERkVvp16+f09eP2qNpUpqYmIhz585h6tSpyMvLQ0REhE1NLS+vv26XiImJwfLlyzF58mRMmjQJLVu2xNq1a9G+fXtLzN/+9jdkZGQgLS0No0aNQuvWrbF69Wo8+OCDVb59RERE5GFMcOuR0oqm7Z1qW8s6pa7KkTqlWtQBrYz2ZG6HVvVHRa8KFl2v6F3/ousVbU/0bnnZcbLvgufd9+rjePe9fSWS2xOtySm6Xpn1TGXXC5Vdk9VV66i6RJ1SX8DfybytSAECrrpmnVIvLy/k5eW530gpEREREVUfZZWOKgOT0mrKlUdSREcERdcpu/6o7BE3V37SECB+rmhVz1R2XVstyK6xKHtEVasnOmlVf1SrkVyZx032MdOq7qlH8oAbnSoLk1IiIiIiWdz8mtLK5OqDPERERETkAThSSkRERCQLR0pVY1JKREREJIsCj00qncXpeyIiIiLSHEdKiYiIiCQxwfmKHLIrelQXTEpdjFblPWSWmNKqhJOrF2HX4lgA2pV6kn18RUvkuDIWz7dPq1JPrry9rvwZ70icJ2JSqh6n74mIiIhIcxwpJSIiIpLEDOcfQuCpDzFgUkpEREQkCafv1eP0PRERERFpjiOlRERERJJw+l49JqVEREREknD6Xj1O3xMRERGR5jhSWkVk15ITpUXdU63q5onWCNSqTqlWRPefXnJ7oly9vqwIrWo7ir6HROuAirbnafVHtagZ6up1RbX6m1YdmOH8/vbE/QYwKSUiIiKShteUqsfpeyIiIiLSHEdKiYiIiCThjU7qMSklIiIikoRJqXqcviciIiIizXGklIiIiEgS3uikHpNSIiIiIkk4fa8ek9JqSqsadiJxstcpu/6oaO1EdyF6jY7oOSXanujxkF0b05VpUTe4Mtpz9fqjov2T/Rnkyp+jss8BosrApJSIiIhIEk7fq8eklIiIiEgSPtFJPd59T0RERESa40gpERERkSS80Uk9JqVEREREkvCaUvU4fU9EREREmuNIqZsT/bYlM06rUk+y2xPdDlf/Zie7NJPsElOix82V97PsUQ2tSkLJfu9qUUoJkH/uyW5PZDu0+OwmOTh9rx6TUiIiIiJJmJSq58qDD0RERETkIThSSkRERCQJb3RSj0kpERERkSScvleP0/dEREREpDmOlBIRERFJosD56XdFRkeqISalRERERJJw+l49JqVOkl1fT3acjwbr1ap2oui1KKWCce5C9vGQXfdUlCtfa6RVnVKt6pnKrj8qu+6p6Htcq89lmfWeZScvrv43jdwbk1IiIiIiSThSqh6TUiIiIiJJWBJKPVeeESMiIiIiD8GRUiIiIiJJOH2vHpNSIiIiIkmYlKrH6XsiIiIi0hxHSomIiIgk4Y1O6jEpraZcuZ6p7LqiorUJ3YVWNS9l1x/VavpJdDtEaLUNWtX6FV2vVnVKReuPyq57Knt7ZX6OimK90KpjhvP70VOTUk7fExEREZHmOFJKREREJAmn79VjUkpEREQkCe++V4/T90RERESkOY6UEhEREUnCkVL1mJQSERERScJrStXj9D0RERERaY4jpS5Gdl1R2fUORWpAiq5Tdt9Ev2G5Sz1OUaL7RYtzpTLIHGHQ6tjKriuqVXuy63tq9dmiRf+0OraiXP1zT0ucvlePSSkRERGRJExK1eP0PRERERFpjiOlRERERJIocP5yCUVGR6ohJqVEREREknD6Xj1O3xMRERGR5jhSSkRERCQJ65Sqx6S0mpJ9wsosPyK7hJPsUk/uMi0iuzST6Dkgu8SUbDKPr6tvg6uXPpJdOkqrElOy95/Mz1Gt/haQa1q4cCH+8Y9/IC8vDx07dsT8+fPRpUsXu7GLFy/Gxx9/jN9++w0AEBkZiRkzZpQbXxU4fU9EREQkiUnSj6NWrlyJ5ORkpKSkYOfOnejYsSPi4+Nx9uxZu/E5OTno378/Nm3ahNzcXISGhuKxxx7D6dOnVaxdDp2iKJ56k1e5ioqKEBAQAF8AugpiRYvYGwXjRNvTa7RekTgt1gmI7xPR9kRHGGVvh+h6ZY8ga7Ve2WR+0+ZIadWsV6uRUtH2SiW3JxJ3TYN1OrJe2ftE1vYqAK4CKCwshL+/v2CrcpTlDu8A8HWyrasAxsKx7YiKikLnzp2xYMECAIDZbEZoaChGjhyJCRMmVPh6k8mEunXrYsGCBRg0aJATvVePI6VERERELqioqMjqp6SkxG5caWkpduzYgbi4OMsyLy8vxMXFITc3V2hdV65cwfXr11GvXj0pfVdD86R04cKFCAsLg9FoRFRUFLZt23bH+FWrViE8PBxGoxEdOnRAZmam1e8vX76MESNG4K677oKvry/atm2LjIyMytwEIiIiIgB/3ejk7A8AhIaGIiAgwPKTlpZmd53nz5+HyWRCcHCw1fLg4GDk5eUJ9fv1119Ho0aNrBLbqqbpjU5l1z9kZGQgKioK6enpiI+Px4EDBxAUFGQTv3XrVvTv3x9paWl48sknsXz5ciQkJGDnzp1o3749ACA5ORnfffcdPv30U4SFheGbb77BK6+8gkaNGqFXr15VvYlERETkQWTWKT116pTV9L3BYHCyZftmzpyJFStWICcnB0aj6EV48mk6Ujpv3jy88MILGDJkiGVEs2bNmliyZInd+HfffRePP/44xo8fjzZt2uCtt97CfffdZ7l+AriZuCYlJaFbt24ICwvDsGHD0LFjxwpHYImIiIhcib+/v9VPeUlpYGAgvL29kZ+fb7U8Pz8fISEhd1zHnDlzMHPmTHzzzTe45557pPVdDc2SUjXXP+Tm5toMK8fHx1vFx8TEYN26dTh9+jQURcGmTZtw8OBBPPbYY+X2paSkxOa6DSIiIiJHmeH8nfeO3mSp1+sRGRmJ7Ozsv/phNiM7OxvR0dHlvm727Nl46623kJWVhU6dOjm4Vvk0m76/0/UP+/fvt/uavLy8Cq+XmD9/PoYNG4a77roLNWrUgJeXFxYvXoyuXbuW25e0tDRMmzbNia2pmKvfxSt6Z7jMdcqun6lV/VHZ3+xEt0O02oBsWtVPFN3PWpwHsrdV9jbIPmaid15rVR1Aq3qmrvxZpdXnoyfWPdWqeH5ycjKSkpLQqVMndOnSBenp6SguLsaQIUMAAIMGDULjxo0t16XOmjULU6dOxfLlyxEWFmbJpWrVqoVatWo5uQXquF3x/Pnz5+Onn37CunXr0LRpU/zwww8YPnz4HS/enThxIpKTky3/LioqQmhoaFV1mYiIiMgpiYmJOHfuHKZOnYq8vDxEREQgKyvLMph38uRJeHn99fX+/fffR2lpKZ555hmrdlJSUpCamlqVXbfQLClVc/1DSEjIHeOvXr2KSZMm4auvvkKPHj0AAPfccw92796NOXPmlJuUGgyGSrt4mIiIiDyHzBudHDVixAiMGDHC7u9ycnKs/n38+HGVa6k8ml1Tqub6h+joaKt4ANi4caMl/vr167h+/brVNwEA8Pb2htnsiZMIREREVJVkloTyNJpO3zt6/cPo0aMRGxuLuXPnokePHlixYgW2b9+ORYsWAbh5l1psbCzGjx8PX19fNG3aFN9//z0+/vhjzJs3T7PtJCIiIqI70zQpdfT6h5iYGCxfvhyTJ0/GpEmT0LJlS6xdu9ZSoxQAVqxYgYkTJ2LAgAG4ePEimjZtiunTp+Oll16q8u0jIiIiz6Ll9H11p1MURdG6E66m7Pm1vgB0FcTKfi676JWtsp/zLvN59Vo9q97Tnmkvuv+0ela97PY0f/zcHfDue+fak31Xvey772WvV2T/afVMe5nbAAD2H4qpvr2K+qfg5nPjHXlmvCxlucMbEP+bWp5rAKZDm+3Qkit/zhMRERGRh3C7klDVnavXMxX5FiM6QiZznY7EiXL1b2yiIwuy68HKrhcq+3xxZbK3Qas6oFqNgMreDq3iRPaf7HXK5qk34ojQqk6pO2BSSkRERCRJ2ROdnG3DE7n6YBAREREReQCOlBIRERFJwrvv1WNSSkRERCQJrylVj9P3RERERKQ5jpQSERERScLpe/WYlFZTssspySwbpFXfRMvKiJJd/F2U6LSNaNF+2eV7RPeL7FJUnkSrEk6itCrhpFUpKi1Kb4muU6uyYFQ+Tt+rx+l7IiIiItIcR0qJiIiIJOH0vXpMSomIiIgkYVKqHqfviYiIiEhzHCklIiIikkSB8zcqKTI6Ug0xKSUiIiKShNP36nH6noiIiIg0x5FSkk52vUvZ3xhFp1VKBeP0ajtSRWR/85RdF1GUVnVjRWh1jorSqg6o7PWKvie1qvGp1XuDXAtHStVjUkpEREQkCYvnq8fpeyIiIiLSHEdKiYiIiCTh9L16TEqJiIiIJOH0vXqcviciIiIizXGklIiIiEgSTt+rx6SUiIiISBIznE8qPXX6nkmpk2TXpRO9nkI0TnYtUJH1ivZN9j4RJbs9d/lGK3oua3XNjyd9SMuuiyl7vVrV7dSqjqrsOC3W6er7zpPe31Q+JqVEREREkvBGJ/WYlBIRERFJYoLzM0ruMgPnKN59T0RERESa40gpERERkSQcKVWPSSkRERGRJLymVD1O3xMRERGR5jhSSg4R+fYm+g1PtFyVVu25C9FpINH9J0r2fnblb9BanVNalCoC5JcXui65PdnHQ4v2PO1zyp1w+l49JqVEREREknD6Xj1XHnwgIiIiIg/BkVIiIiIiSfiYUfWYlBIRERFJYgKgk9CGJ+L0PRERERFpjiOlRERERJLwRif1mJQSERERScLpe/WYlLo50RNb9DoOkVqWstcpyl2uRZFdV1R2/VGt6p560oe0q9cVlV0vVKv2XD3OVddJVFmYlBIRERFJwpFS9ZiUEhEREUnCa0rVc5cZTyIiIiKqxjhSSkRERCQJp+/VY1JKREREJIkC56ffFRkdqYY4fU9EREREmuNIKREREZEkMqbeOX1P1Yro1IBorUiZ9f9Eh99l310ouy4mOUf0+HrSdI1Wd9Sy/mjVrFfm8dVinZXRnidiUqqeJ/09ICIiIiIXxZFSIiIiIknMcP7ue08dsWZSSkRERCQJp+/V4/Q9EREREWmOI6VEREREknCkVD0mpURERESS8JpS9Th9T0RERESa40ipi5Fdf1QLsr/hiX5z8tTpDllE958rn3uuTvY5Krs91h91jla1RbXgDttQWWTsG0/dv0xKiYiIiCRhUqoep++JiIiISHMcKSUiIiKSxARAcbINTx0pZVJKREREJAmTUvU4fU9EREREmuNIKREREZEkvNFJPSalbk52mR+ZJVREh+k99c3pKdyhFJVW5chkvzc8rdSTViWcRNbr6uXDqHycvleP0/dEREREpDmOlBIRERFJYobzI6XOvr66YlJKREREJIkZgM7JNjw1KdV8+n7hwoUICwuD0WhEVFQUtm3bVm7sf//7Xzz99NMICwuDTqdDenq6020SERERuQNH859Vq1YhPDwcRqMRHTp0QGZmZhX11D5Nk9KVK1ciOTkZKSkp2LlzJzp27Ij4+HicPXvWbvyVK1fQvHlzzJw5EyEhIVLaJCIiIpLFJOnHUY7mP1u3bkX//v0xdOhQ7Nq1CwkJCUhISMBvv/2mYu1y6BRF0WyUOCoqCp07d8aCBQsAAGazGaGhoRg5ciQmTJhwx9eGhYVhzJgxGDNmjLQ2yxQVFSEgIAC+qHgIXvSOYB/BONntiX7r0EtsT6t9Irqtou152npFafVNlnff2+Ld99VzvdcltiXat1LJ7YlsAyB3n4i0pwC4CqCwsBD+/v6CrcpRljvUhJzp+ytwbDsczX8SExNRXFyM9evXW5bdf//9iIiIQEZGhpNboI5m15SWlpZix44dmDhxomWZl5cX4uLikJubW6VtlpSUoKSkxPLvwsJCAGLXdIhm9KJxoh8Isv8YaVF+QnZyI/vblWh7otsh+1yRvV5RWiWlrlwiRau+MSl1Ls6Vt1f2tmoVJ/tzr6K4st9rON4m5TO3rI2ioiKr5QaDAQaDwSZeTf6Tm5uL5ORkq2Xx8fFYu3atU313hmZJ6fnz52EymRAcHGy1PDg4GPv376/SNtPS0jBt2jSb5ddU9YKIiIi0dOHCBQQEBFTpOvV6PUJCQpCXlyelvVq1aiE0NNRqWUpKClJTU21i1eQ/eXl5duNl9V8N3n0PYOLEiVbfFgoKCtC0aVOcPHmyyk9qslZUVITQ0FCcOnWqyqdiyBqPhWvh8XAdPBauo7CwEE2aNEG9evWqfN1GoxHHjh1DaanoxRJ3pigKdDrrCwHsjZK6E82S0sDAQHh7eyM/P99qeX5+frk3MVVWm+UNhwcEBPADxkX4+/vzWLgIHgvXwuPhOngsXIeXlzYXGhmNRhiNxipfr5r8JyQkRGoOJoNmd9/r9XpERkYiOzvbssxsNiM7OxvR0dEu0yYRERGRK1OT/0RHR1vFA8DGjRs1zZc0nb5PTk5GUlISOnXqhC5duiA9PR3FxcUYMmQIAGDQoEFo3Lgx0tLSANy8kHfv3r2W/z99+jR2796NWrVq4e677xZqk4iIiMjdOJpTjR49GrGxsZg7dy569OiBFStWYPv27Vi0aJFm26BpUpqYmIhz585h6tSpyMvLQ0REBLKysiwX3p48edJqCP6PP/7Avffea/n3nDlzMGfOHMTGxiInJ0eoTREGgwEpKSluf+1GdcBj4Tp4LFwLj4fr4LFwHZ58LBzNqWJiYrB8+XJMnjwZkyZNQsuWLbF27Vq0b99eq03Qtk4pERERERHgAo8ZJSIiIiJiUkpEREREmmNSSkRERESaY1JKRERERJrz2KR04cKFCAsLg9FoRFRUFLZt23bH+FWrViE8PBxGoxEdOnRAZmZmFfXU/TlyLBYvXoyHHnoIdevWRd26dREXF1fhsSNxjr4vyqxYsQI6nQ4JCQmV20EP4+jxKCgowPDhw9GwYUMYDAa0atWKn1WSOHos0tPT0bp1a/j6+iI0NBRjx47FtWt8eLWzfvjhB/Ts2RONGjWCTqcTek57Tk4O7rvvPhgMBtx9991YtmxZpfeTVFI80IoVKxS9Xq8sWbJE+e9//6u88MILSp06dZT8/Hy78T/++KPi7e2tzJ49W9m7d68yefJkxcfHR9mzZ08V99z9OHosnn32WWXhwoXKrl27lH379imDBw9WAgIClN9//72Ke+5+HD0WZY4dO6Y0btxYeeihh5TevXtXTWc9gKPHo6SkROnUqZPSvXt3ZcuWLcqxY8eUnJwcZffu3VXcc/fj6LH47LPPFIPBoHz22WfKsWPHlA0bNigNGzZUxo4dW8U9dz+ZmZnKG2+8oaxZs0YBoHz11Vd3jD969KhSs2ZNJTk5Wdm7d68yf/58xdvbW8nKyqqaDpNDPDIp7dKlizJ8+HDLv00mk9KoUSMlLS3Nbnzfvn2VHj16WC2LiopSXnzxxUrtpydw9Fjc7saNG0rt2rWVjz76qLK66DHUHIsbN24oMTExyocffqgkJSUxKZXI0ePx/vvvK82bN1dKS0urqosew9FjMXz4cOWRRx6xWpacnKw88MADldpPTyOSlL722mtKu3btrJYlJiYq8fHxldgzUsvjpu9LS0uxY8cOxMXFWZZ5eXkhLi4Oubm5dl+Tm5trFQ8A8fHx5caTGDXH4nZXrlzB9evXUa9evcrqpkdQeyzefPNNBAUFYejQoVXRTY+h5nisW7cO0dHRGD58OIKDg9G+fXvMmDEDJpOpqrrtltQci5iYGOzYscMyxX/06FFkZmaie/fuVdJn+gv/flcvmj7RSQvnz5+HyWSyecJTcHAw9u/fb/c1eXl5duPz8vIqrZ+eQM2xuN3rr7+ORo0a2XzokGPUHIstW7bgX//6F3bv3l0FPfQsao7H0aNH8d1332HAgAHIzMzE4cOH8corr+D69etISUmpim67JTXH4tlnn8X58+fx4IMPQlEU3LhxAy+99BImTZpUFV2mW5T397uoqAhXr16Fr6+vRj0jezxupJTcx8yZM7FixQp89dVXMBqNWnfHo1y6dAkDBw7E4sWLERgYqHV3CIDZbEZQUBAWLVqEyMhIJCYm4o033kBGRobWXfM4OTk5mDFjBt577z3s3LkTa9aswddff4233npL664RuTSPGykNDAyEt7c38vPzrZbn5+cjJCTE7mtCQkIciicxao5FmTlz5mDmzJn49ttvcc8991RmNz2Co8fiyJEjOH78OHr27GlZZjabAQA1atTAgQMH0KJFi8rttBtT895o2LAhfHx84O3tbVnWpk0b5OXlobS0FHq9vlL77K7UHIspU6Zg4MCBeP755wEAHTp0QHFxMYYNG4Y33njD6vnjVLnK+/vt7+/PUVIX5HHvDL1ej8jISGRnZ1uWmc1mZGdnIzo62u5roqOjreIBYOPGjeXGkxg1xwIAZs+ejbfeegtZWVno1KlTVXTV7Tl6LMLDw7Fnzx7s3r3b8tOrVy88/PDD2L17N0JDQ6uy+25HzXvjgQcewOHDhy1fDgDg4MGDaNiwIRNSJ6g5FleuXLFJPMu+LCiKUnmdJRv8+13NaH2nlRZWrFihGAwGZdmyZcrevXuVYcOGKXXq1FHy8vIURVGUgQMHKhMmTLDE//jjj0qNGjWUOXPmKPv27VNSUlJYEkoSR4/FzJkzFb1er3z55ZfKmTNnLD+XLl3SahPchqPH4na8+14uR4/HyZMnldq1aysjRoxQDhw4oKxfv14JCgpS3n77ba02wW04eixSUlKU2rVrK59//rly9OhR5ZtvvlFatGih9O3bV6tNcBuXLl1Sdu3apezatUsBoMybN0/ZtWuXcuLECUVRFGXChAnKwIEDLfFlJaHGjx+v7Nu3T1m4cCFLQrkwj0xKFUVR5s+frzRp0kTR6/VKly5dlJ9++snyu9jYWCUpKckq/osvvlBatWql6PV6pV27dsrXX39dxT12X44ci6ZNmyoAbH5SUlKqvuNuyNH3xa2YlMrn6PHYunWrEhUVpRgMBqV58+bK9OnTlRs3blRxr92TI8fi+vXrSmpqqtKiRQvFaDQqoaGhyiuvvKL8+eefVd9xN7Np0ya7fwPK9n9SUpISGxtr85qIiAhFr9crzZs3V5YuXVrl/SYxOkXhXAIRERERacvjriklIiIiItfDpJSIiIiINMeklIiIiIg0x6SUiIiIiDTHpJSIiIiINMeklIiIiIg0x6SUiIiIiDTHpJSIiIiINMeklIjcQk5ODnQ6HQoKCu4YFxYWhvT09CrpExERiWNSSkRVZvDgwdDpdNDpdNDr9bj77rvx5ptv4saNG063HRMTgzNnziAgIAAAsGzZMtSpU8cm7ueff8awYcOcXh8REclVQ+sOEJFnefzxx7F06VKUlJQgMzMTw4cPh4+PDyZOnOhUu3q9HiEhIRXGNWjQwKn1EBFR5eBIKRFVKYPBgJCQEDRt2hQvv/wy4uLisG7dOgDAn3/+iUGDBqFu3bqoWbMmnnjiCRw6dMjy2hMnTqBnz56oW7cu/Pz80K5dO2RmZgKwnr7PycnBkCFDUFhYaBmZTU1NBWA7fX/y5En07t0btWrVgr+/P/r27Yv8/HzL71NTUxEREYFPPvkEYWFhCAgIQL9+/XDp0iXhbc7JyYFer8fmzZsty2bPno2goCCrdREReTImpUSkKV9fX5SWlgK4Ob2/fft2rFu3Drm5uVAUBd27d8f169cBAMOHD0dJSQl++OEH7NmzB7NmzUKtWrVs2oyJiUF6ejr8/f1x5swZnDlzBuPGjbOJM5vN6N27Ny5evIjvv/8eGzduxNGjR5GYmGgVd+TIEaxduxbr16/H+vXr8f3332PmzJnC29itWzeMGTMGAwcORGFhIXbt2oUpU6bgww8/RHBwsCO7i4jIbXH6nog0oSgKsrOzsWHDBowcORKHDh3CunXr8OOPPyImJgYA8NlnnyE0NBRr165Fnz59cPLkSTz99NPo0KEDAKB58+Z229br9QgICIBOp7vjlH52djb27NmDY8eOITQ0FADw8ccfo127dvj555/RuXNnADeT12XLlqF27doAgIEDByI7OxvTp08X3t63334bGzduxLBhw/Dbb78hKSkJvXr1En49EZG740gpEVWp9evXo1atWjAajXjiiSeQmJiI1NRU7Nu3DzVq1EBUVJQltn79+mjdujX27dsHABg1ahTefvttPPDAA0hJScGvv/7qVF/27duH0NBQS0IKAG3btkWdOnUs6wRuTvmXJaQA0LBhQ5w9e9ahden1enz22WdYvXo1rl27hnfeecepvhMRuRsmpURUpR5++GHs3r0bhw4dwtWrV/HRRx/Bz89P6LXPP/88jh49ioEDB2LPnj3o1KkT5s+fX8k9Bnx8fKz+rdPpYDabHW5n69atAICLFy/i4sWLUvpGROQumJQSUZXy8/PD3XffjSZNmqBGjb+uIGrTpg1u3LiB//znP5ZlFy5cwIEDB9C2bVvLstDQULz00ktYs2YNXn31VSxevNjuevR6PUwm0x370qZNG5w6dQqnTp2yLNu7dy8KCgqs1inDkSNHMHbsWCxevBhRUVFISkpSldgSEbkrJqVE5BJatmyJ3r1744UXXsCWLVvwyy+/4LnnnkPjxo3Ru3dvAMCYMWOwYcMGHDt2DDt37sSmTZvQpk0bu+2FhYXh8uXLyM7Oxvnz53HlyhWbmLi4OHTo0AEDBgzAzp07sW3bNgwaNAixsbHo1KmTtG0zmUx47rnnEB8fjyFDhmDp0qX49ddfMXfuXGnrICKq7piUEpHLWLp0KSIjI/Hkk08iOjoaiqIgMzPTMn1uMpkwfPhwtGnTBo8//jhatWqF9957z25bMTExeOmll5CYmIgGDRpg9uzZNjE6nQ7/+7//i7p166Jr166Ii4tD8+bNsXLlSof6vWzZMuh0unJ/P336dJw4cQIffPABgJvXpC5atAiTJ0/GL7/84tC6iIjclU5RFEXrThARVWcpKSn4/vvvkZOTo3VXiIiqLZaEIiJy0r///W8sWLBA624QEVVrHCklIiIiIs3xmlIiIiIi0hyTUiIiIiLSHJNSIiIiItIck1IiIiIi0hyTUiIiIiLSHJNSIiIiItIck1IiIiIi0hyTUiIiIiLSHJNSIiIiItLc/wNEi0hL4OfQrQAAAABJRU5ErkJggg=="
      },
      "metadata" : { },
      "output_type" : "display_data"
    } ],
    "execution_count" : 4
  }, {
    "metadata" : {
      "ExecuteTime" : {
        "end_time" : "2025-03-17T18:08:23.577161Z",
        "start_time" : "2025-03-17T18:08:23.572046Z"
      }
    },
    "cell_type" : "code",
    "source" : "# fn.dprint(print_shape=True, print_memory_map=True)",
    "id" : "6a4ce34cb39e8cb5",
    "outputs" : [ ],
    "execution_count" : 5
  }, {
    "metadata" : {
      "ExecuteTime" : {
        "end_time" : "2025-03-17T18:08:26.824757Z",
        "start_time" : "2025-03-17T18:08:23.617272Z"
      }
    },
    "cell_type" : "code",
    "source" : [ "nt_test = np.array(50)\n", "alpha_test = np.array(1.0)\n", "fn.trust_input = True\n", "%timeit fn(nt_test, alpha_test)\n", "fn.trust_input = False" ],
    "id" : "cf9033d5598071da",
    "outputs" : [ {
      "name" : "stdout",
      "output_type" : "stream",
      "text" : [ "3.89 ms ± 104 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n" ]
    } ],
    "execution_count" : 6
  }, {
    "metadata" : {
      "ExecuteTime" : {
        "end_time" : "2025-03-17T18:08:38.878998Z",
        "start_time" : "2025-03-17T18:08:26.835721Z"
      }
    },
    "cell_type" : "code",
    "source" : [ "# with pytensor.config.change_flags(optimizer_verbose=True, optimizer_verbose_ignore=\"constant_folding,MergeOptimizer\"):\n", "fn_numba = build_heat_equation_fn(mode=\"NUMBA\")\n", "fn_numba(nt_test, alpha_test);" ],
    "id" : "83c05e38d2e02ff7",
    "outputs" : [ {
      "name" : "stderr",
      "output_type" : "stream",
      "text" : [ "/tmp/tmp5xa7er2j:21: NumbaWarning: \u001B[1m\u001B[1mCannot cache compiled function \"scan\" as it uses dynamic globals (such as ctypes pointers and large global arrays)\u001B[0m\u001B[0m\n", "  tensor_variable_9 = scan(nt, tensor_variable_8, tensor_variable_5)\n" ]
    } ],
    "execution_count" : 7
  }, {
    "metadata" : {
      "ExecuteTime" : {
        "end_time" : "2025-03-17T18:08:38.891948Z",
        "start_time" : "2025-03-17T18:08:38.888553Z"
      }
    },
    "cell_type" : "code",
    "source" : "# fn_numba.dprint(print_memory_map=True, print_shape=True)",
    "id" : "fdc0cb93efa8c937",
    "outputs" : [ ],
    "execution_count" : 8
  }, {
    "metadata" : {
      "ExecuteTime" : {
        "end_time" : "2025-03-17T18:08:40.840805Z",
        "start_time" : "2025-03-17T18:08:38.935833Z"
      }
    },
    "cell_type" : "code",
    "source" : [ "fn_numba.trust_input = True\n", "%timeit fn_numba(nt_test, alpha_test)\n", "fn_numba.trust_input = False" ],
    "id" : "51391ca854a2ecde",
    "outputs" : [ {
      "name" : "stdout",
      "output_type" : "stream",
      "text" : [ "2.33 ms ± 87.1 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n" ]
    } ],
    "execution_count" : 9
  }, {
    "metadata" : {
      "ExecuteTime" : {
        "end_time" : "2025-03-17T18:08:40.855892Z",
        "start_time" : "2025-03-17T18:08:40.853318Z"
      }
    },
    "cell_type" : "code",
    "source" : "",
    "id" : "5eaa889dcc5d584a",
    "outputs" : [ ],
    "execution_count" : null
  } ],
  "metadata" : {
    "kernelspec" : {
      "display_name" : "Python 3",
      "language" : "python",
      "name" : "python3"
    },
    "language_info" : {
      "codemirror_mode" : {
        "name" : "ipython",
        "version" : 2
      },
      "file_extension" : ".py",
      "mimetype" : "text/x-python",
      "name" : "python",
      "nbconvert_exporter" : "python",
      "pygments_lexer" : "ipython2",
      "version" : "2.7.6"
    }
  },
  "nbformat" : 4,
  "nbformat_minor" : 5
}