{ "cells": [ { "attachments": {}, "cell_type": "markdown", "id": "8414ba1c-c56d-46e3-a099-d00044905eda", "metadata": {}, "source": [ "# Parallelized Cubic Spline Interpolation" ] }, { "attachments": {}, "cell_type": "markdown", "id": "1cf56b82-e6a1-405c-ba54-9a31decb1f64", "metadata": {}, "source": [ "A part of the Fast EMRI waveforms package is parallelized cubic spline interpolation. This generally means fitting and evaluating many splines in parallel with the same input x array. This is available for GPUs and CPUs (not parallelized for CPU). The user can perform this operation entirely in Python while leveraging [CuPy](https://cupy.dev/) for GPUs. However, the evaluation will not be as efficient as when it is implemented properly in a customized kernel. The spline class ([CubicSplineInterpolant](https://bhptoolkit.org/FastEMRIWaveforms/html/user/sum.html#few.summation.interpolatedmodesum.CubicSplineInterpolant)) can provide an 1D flattened array of all spline coefficients for use in a custom CUDA kernel. " ] }, { "cell_type": "code", "execution_count": null, "id": "3c526d1a", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "import few\n", "\n", "# tune few configuration\n", "cfg_set = few.get_config_setter(reset=True)\n", "\n", "# Uncomment if you want to force CPU or GPU usage\n", "# Leave commented to let FEW automatically select the best available hardware\n", "# - To force CPU usage:\n", "# cfg_set.enable_backends(\"cpu\")\n", "# - To force GPU usage with CUDA 12.x\n", "# cfg_set.enable_backends(\"cuda12x\", \"cpu\")\n", "# - To force GPU usage with CUDA 11.x\n", "# cfg_set.enable_backends(\"cuda11x\", \"cpu\")\n", "\n", "cfg_set.set_log_level(\"info\");" ] }, { "cell_type": "code", "execution_count": null, "id": "19ba8ad4", "metadata": {}, "outputs": [], "source": [ "from few.trajectory.inspiral import EMRIInspiral\n", "from few.trajectory.ode import SchwarzEccFlux\n", "from few.amplitude.romannet import RomanAmplitude\n", "\n", "traj = EMRIInspiral(func=SchwarzEccFlux)\n", "amp = RomanAmplitude()\n", "\n", "# parameters\n", "M = 1e5\n", "mu = 1e1\n", "p0 = 10.0\n", "e0 = 0.7\n", "theta = np.pi / 3.0\n", "phi = np.pi / 2.0\n", "\n", "# get trajectory\n", "t, p, e, x, Phi_phi, Phi_theta, Phi_r = traj(M, mu, 0.0, p0, e0, 1.0)\n", "# Schwarzschild\n", "a = 0.0\n", "\n", "teuk_modes = amp(a, p, e, x)" ] }, { "cell_type": "code", "execution_count": 18, "id": "18b9882f-3e6d-4ca2-8730-c034ae96e649", "metadata": {}, "outputs": [], "source": [ "from few.summation.interpolatedmodesum import CubicSplineInterpolant\n", "\n", "# let's take the amplitudes from the last step and spline those.\n", "# We have to arange everything in the shape (ninterps, length)\n", "# We will split real and imaginary components\n", "\n", "interp_in = np.zeros((teuk_modes.shape[1] * 2, teuk_modes.shape[0]))\n", "\n", "interp_in[: teuk_modes.shape[1], :] = teuk_modes.T.real\n", "interp_in[teuk_modes.shape[1] :, :] = teuk_modes.T.imag\n", "\n", "spline = CubicSplineInterpolant(t, interp_in)" ] }, { "cell_type": "code", "execution_count": 19, "id": "04ebee73-84eb-4e73-945a-8d59eef9795e", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGdCAYAAAAfTAk2AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAWHZJREFUeJzt3Xl4VOX5//H3mUlmsm9kh4QEwr5DIAQsqESDaC2Wti60CqVqVVTEDay7VajFuleqrVL7VVH7U+tKi+AGhH2HEFkCYclCCMlkIdvM+f0RGI3skjCT4fO6rnOROec5z7nPTOLcPudZDNM0TURERER8iMXTAYiIiIi0NCU4IiIi4nOU4IiIiIjPUYIjIiIiPkcJjoiIiPgcJTgiIiLic5TgiIiIiM9RgiMiIiI+x8/TAXiCy+Vi3759hIaGYhiGp8MRERGRU2CaJpWVlSQmJmKxnLiN5pxMcPbt20dSUpKnwxAREZEfYffu3XTo0OGEZc7JBCc0NBRoeoPCwsI8HI2IiIicCofDQVJSkvt7/ETOyQTnyGOpsLAwJTgiIiJtzKl0L1EnYxEREfE5SnBERETE5yjBEREREZ9zTvbBORWmadLY2IjT6fR0KOIFrFYrfn5+mlZARKSNUIJzDPX19RQWFlJTU+PpUMSLBAUFkZCQgM1m83QoIiJyEkpwfsDlcpGfn4/VaiUxMRGbzab/az/HmaZJfX09+/fvJz8/ny5dupx0gikREfEsJTg/UF9fj8vlIikpiaCgIE+HI14iMDAQf39/du3aRX19PQEBAZ4OSURETkD/G3oc+j90+SH9ToiItB36L7aIiIj4HCU4PuT8889nypQpng7jjBiGwQcffOCRaz/88MP079/fI9cWEZGWpQTHh7z33ns89thjp1R2586dGIbB2rVrWzeoVqakREREjkWdjH1IVFSUR67b0NCAv7+/R64tIiJyLGrB8SHff0SVkpLCE088wW9/+1tCQ0NJTk7m5ZdfdpdNTU0FYMCAARiGwfnnn+8+9ve//50ePXoQEBBA9+7d+etf/+o+dqTl5+2332bkyJEEBATwxhtvMGfOHCIiIvjggw/o0qULAQEBZGdns3v37mYxvvTSS3Tu3BmbzUa3bt3417/+dcJ7uvfee+natStBQUF06tSJBx54gIaGBgDmzJnDI488wrp16zAMA8MwmDNnDgDl5eX87ne/IyYmhrCwMC688ELWrVvXrO6ZM2cSFxdHaGgokyZNora29rTebxEROdpX3+7n4Q838fH6fR6NQwnOKTBNk5r6xrO+maZ5RnE/9dRTpKens2bNGm6++WZuuukm8vLyAFi+fDkAn3/+OYWFhbz33nsAvPHGGzz44IM8/vjj5Obm8sQTT/DAAw/wz3/+s1nd06ZN4/bbbyc3N5fs7GwAampqePzxx3n99ddZvHgx5eXlXHXVVe5z3n//fW6//XbuvPNONm7cyI033sjEiRP54osvjnsPoaGhzJkzh82bN/Pss8/yyiuv8PTTTwNw5ZVXcuedd9KrVy8KCwspLCzkyiuvBOCXv/wlJSUlfPbZZ6xatYqBAwcyatQoysrKAHjnnXd4+OGHeeKJJ1i5ciUJCQnNEjkREflxDqz/HzHLZ3Jw9X88GoceUZ2CQw1Oej7437N+3c2PZhNk+/Ef0ZgxY7j55puBppaQp59+mi+++IJu3boRExMDQLt27YiPj3ef89BDD/HUU0/x85//HGhq6dm8eTN/+9vfuO6669zlpkyZ4i5zRENDAy+88AIZGRkA/POf/6RHjx4sX76cIUOGMGvWLCZMmOCOaerUqSxdupRZs2ZxwQUXHPMe7r//fvfPKSkp3HXXXcydO5d77rmHwMBAQkJC8PPza3YPixYtYvny5ZSUlGC32wGYNWsWH3zwAf/+97+54YYbeOaZZ5g0aRKTJk0C4I9//COff/65WnFERM5QbOVGLrLOZ0O1Z+cLUwuOD+vbt6/7Z8MwiI+Pp6Sk5Ljlq6ur2b59O5MmTSIkJMS9/fGPf2T79u3Nyqanpx91vp+fH4MHD3a/7t69OxEREeTm5gKQm5vL8OHDm50zfPhw9/Fjefvttxk+fDjx8fGEhIRw//33U1BQcML7XrduHVVVVbRr167ZfeTn57vvIzc3152IHZGZmXnCekVE5OSWt59An7p/8N8Ot3s0DrXgnIJAfyubH832yHXPxA87/hqGgcvlOm75qqoqAF555ZWjvvyt1uaxBAcHn1FspyInJ4fx48fzyCOPkJ2dTXh4OHPnzuWpp5464XlVVVUkJCTw5ZdfHnUsIiKidYIVEREAXId7V3h6mSMlOKfAMIwzelTkjY4sGPn91dLj4uJITExkx44djB8//rTrbGxsZOXKlQwZMgSAvLw8ysvL6dGjBwA9evRg8eLFzR51LV68mJ49ex6zviVLltCxY0f+8Ic/uPft2rXrqPv44YrvAwcOpKioCD8/P1JSUo5Zd48ePVi2bBnXXnute9/SpUtP/WZFROSYnIf7j1qU4IgnxMbGEhgYyLx58+jQoQMBAQGEh4fzyCOPcNtttxEeHs7o0aOpq6tj5cqVHDx4kKlTp56wTn9/f2699Vaee+45/Pz8mDx5MkOHDnUnPHfffTe/+tWvGDBgAFlZWXz00Ue89957fP7558esr0uXLhQUFDB37lwGDx7MJ598wvvvv9+sTEpKCvn5+axdu5YOHToQGhpKVlYWmZmZjB07lieffJKuXbuyb98+PvnkE6644grS09O5/fbbmTBhAunp6QwfPpw33niDTZs20alTp5Z5g0VEzlF9ij/gdf+PqSy9FDj2/8CeDeqDc47y8/Pjueee429/+xuJiYn87Gc/A+B3v/sdf//733nttdfo06cPI0eOZM6cOe5h5ScSFBTEvffeyzXXXMPw4cMJCQnh7bffdh8fO3Yszz77LLNmzaJXr1787W9/47XXXms2RP37Lr/8cu644w4mT55M//79WbJkCQ888ECzMuPGjWP06NFccMEFxMTE8NZbb2EYBp9++ikjRoxg4sSJdO3alauuuopdu3YRFxcHNI3AeuCBB7jnnnsYNGgQu3bt4qabbvqR76aIiBwRdaiAEdYNtKvbffLCrcgwz3QschvkcDgIDw+noqKCsLCwZsdqa2vJz88nNTVVK0afhjlz5jBlyhTKy8s9HUqr0e+GiMjJLX3pRoYWzyUn4Voyb3y+Res+0ff3D6kFR0RERFqOeXgwi4f74CjBERERkRZjuBOcMxsJfKaU4EiLmDBhgk8/nhIRkVPkTnA8m2IowREREZGWowRHREREfI2JQaNpAYsSHBEREfER/2k/lbS6/2NVx995NA4lOCIiItJiXO5BVBpFJSIiIj7CpaUaRERExNcML32HUf7LoPRqoLPH4lALzjnu4Ycfpn///qd1zvnnn8+UKVM8HoeIiHifDofyGG1dQUTtHo/GoRacc9xdd93FrbfeelrnvPfee/j7+7dSRCIi0pYZNHXCMT08TFwJzjnKNE2cTichISGEhISc1rlRUVGtFJWIiLR5h+fBMTRMXFpKXV0dt912G7GxsQQEBHDeeeexYsUKAL788ksMw+Czzz5j0KBB2O12Fi1adNSjocbGRm677TYiIiJo164d9957L9dddx1jx451l/nhI6qUlBSeeOIJfvvb3xIaGkpycjIvv/xys9juvfdeunbtSlBQEJ06deKBBx6goaGhNd8OERHxAC3V0BbVVx9/a6g9jbKHTl72R7jnnnv4f//v//HPf/6T1atXk5aWRnZ2NmVlZe4y06ZNY+bMmeTm5tK3b9+j6vjTn/7EG2+8wWuvvcbixYtxOBx88MEHJ732U089RXp6OmvWrOHmm2/mpptuIi8vz308NDSUOXPmsHnzZp599lleeeUVnn766R91nyIi4s2aRlEZekTVhjyRePxjXS6G8e9+9/rPadBQc+yyHc+DiZ989/qZPlBzoHmZhytOK7Tq6mpeeukl5syZwyWXXALAK6+8wvz58/nHP/7B4MGDAXj00Ue56KKLjlvP888/z/Tp07niiisAeOGFF/j0009Pev0xY8Zw8803A02tNU8//TRffPEF3bp1A+D+++93l01JSeGuu+5i7ty53HPPPad1nyIi4t0ML1mqQQmOj9i+fTsNDQ0MHz7cvc/f358hQ4aQm5vrTnDS09OPW0dFRQXFxcUMGTLEvc9qtTJo0CBcR2ZuOo7vtwYZhkF8fDwlJSXufW+//TbPPfcc27dvp6qqisbGRsLCwk77PkVExLu5ExwP98FRgnM67tt3/GM/fNZ497YTlP3Bhz5lw4+P6TQFBwe3Sr0/HFVlGIY7KcrJyWH8+PE88sgjZGdnEx4ezty5c3nqqadaJRYREfGcv8Y+xFdbCnm0Qz8yPBiH+uCcDlvw8Tf/gNMoG3jysqepc+fO2Gw2Fi9e7N7X0NDAihUr6Nmz5ynVER4eTlxcnLtjMoDT6WT16tWnHc/3LVmyhI4dO/KHP/yB9PR0unTpwq5du86oThER8U6NWKnDhuHn2elEzkqC8+KLL5KSkkJAQAAZGRksX778uGU3bdrEuHHjSElJwTAMnnnmmTOu81wQHBzMTTfdxN133828efPYvHkz119/PTU1NUyaNOmU67n11luZMWMG//nPf8jLy+P222/n4MGDZ7SmSJcuXSgoKGDu3Lls376d5557jvfff/9H1yciIt7LW5ZqaPUE5+2332bq1Kk89NBDrF69mn79+pGdnd2sf8b31dTU0KlTJ2bOnEl8fHyL1HmumDlzJuPGjeM3v/kNAwcOZNu2bfz3v/8lMjLylOu49957ufrqq7n22mvJzMwkJCSE7OxsAgICTn7ycVx++eXccccdTJ48mf79+7NkyRIeeOCBH12fiIh4rzHlc3na/0Wiy86s9f9MGaZ5ONVqJRkZGQwePJgXXngBAJfLRVJSErfeeivTpk074bkpKSlMmTLlqGUBzqROAIfDQXh4OBUVFUd1dK2trSU/P5/U1NQz+lL3FS6Xix49evCrX/2Kxx57zNPheJR+N0RETm7jjPPpXbeGlYOeJP2nN7Zo3Sf6/v6hVm3Bqa+vZ9WqVWRlZX13QYuFrKwscnJyzlqddXV1OByOZpsc265du3jllVf49ttv2bBhAzfddBP5+flcc801ng5NRETagCNLNXh6mHirXr20tBSn00lcXFyz/XFxcRQVFZ21OmfMmEF4eLh7S0pK+lHXPhdYLBbmzJnD4MGDGT58OBs2bODzzz+nR48eng5NRETaAM2DcxZNnz6dqVOnul87HA4lOceRlJTUbCSWiIjI6Tncydjq2aUaWjXBiY6Oxmq1Ulxc3Gx/cXHxcTsQt0addrsdu93+o64nIiIip85yZLFNXx5FZbPZGDRoEAsWLHDvc7lcLFiwgMzMTK+pU0RERFrK4bFLHl5ss9UfUU2dOpXrrruO9PR0hgwZwjPPPEN1dTUTJ04E4Nprr6V9+/bMmDEDaOpEvHnzZvfPe/fuZe3atYSEhJCWlnZKdbaEVh5cJm2QfidERE7OcLfg+HgfnCuvvJL9+/fz4IMPUlRURP/+/Zk3b567k3BBQQGW761XsW/fPgYMGOB+PWvWLGbNmsXIkSP58ssvT6nOM3FkyYGamhoCAwNPUlrOJTU1TYun/nBZChER+c5D4Y+Tu7eMpxOGeTSOVp8HxxudbBx9YWEh5eXlxMbGEhQU5PHniOJZpmlSU1NDSUkJERERJCQkeDokERGvdfkLi1i/p4LXJgzmgu6xLVr36cyDc06MojpdRzorn+szI0tzERERP7pzvIjIucLpamo38XTbgBKcYzAMg4SEBGJjY2loaPB0OOIF/P39sXp4yKOISFtwVc2bWP2KCKmMAlq2Bed0KME5AavVqi81ERGR03Be/WJS/Xax8dAEj8bh2S7OIiIi4lPcSzVYfHipBhERETm3GBzpg+PZJyBKcERERKTFuGcyVguOiIiI+Ap3C45FLTgiIiLiI470wVELjoiIiPiMIy04Fg+34GiYuIiIiLSYG20zKa2o5uV2PTwahxIcERERaTGlRFBIAIZfgEfj0CMqERERaTEuU0s1iIiIiI+5ofEtTL9qbDVdgXCPxaEER0RERFrMFa7/EeXnIL/+Lo/GoUdUIiIi0mK+G0WlYeIiIiLiIyyaB0dERER8jWYyFhEREZ9jMfWISkRERHyM+xGVoQRHREREfMR3j6g8m2JomLiIiIi0mMsb/4TT5eSN0ASPxqEER0RERFpMvhlPo2li8bN5NA49ohIREZEWo6UaRERExOdMtb6NCwNLw1DAcwtuKsERERGRFmG6XEz2+w8AB50zPRqLHlGJiIhIi3A6ne6fNQ+OiIiI+ASXy/XdC81kLCIiIr7A5fquBcdqVQuOiIiI+ADTpUdUIiIi4mNczRIcz45jUoIjIiIiLeL7fXC0VIOIiIj4BJdfEGPqnsDA5AMPz2SsBEdERERahAsLm80UACwaRSUiIiK+4MgyDQAWLdUgIiIivsBVV8XN1g9wGRYM41KPxnJWWnBefPFFUlJSCAgIICMjg+XLl5+w/Lvvvkv37t0JCAigT58+fPrpp82OT5gwAcMwmm2jR49uzVsQERGRkzDqKrnH/x3usr7j6VBaP8F5++23mTp1Kg899BCrV6+mX79+ZGdnU1JScszyS5Ys4eqrr2bSpEmsWbOGsWPHMnbsWDZu3Nis3OjRoyksLHRvb731VmvfioiIiJzAkVFULi/oAdPqEfzlL3/h+uuvZ+LEifTs2ZPZs2cTFBTEq6++eszyzz77LKNHj+buu++mR48ePPbYYwwcOJAXXnihWTm73U58fLx7i4yMbO1bERERkRM4MtGfiYc74NDKCU59fT2rVq0iKyvruwtaLGRlZZGTk3PMc3JycpqVB8jOzj6q/JdffklsbCzdunXjpptu4sCBA8eNo66uDofD0WwTERGRluU83ILj9PUWnNLSUpxOJ3Fxcc32x8XFUVRUdMxzioqKTlp+9OjRvP766yxYsIA//elPfPXVV1xyySXNVjH9vhkzZhAeHu7ekpKSzvDORERE5IfMwwmON7TgtMlRVFdddZX75z59+tC3b186d+7Ml19+yahRo44qP336dKZOnep+7XA4lOSIiIi0MNPVCIDLCxKcVm3BiY6Oxmq1Ulxc3Gx/cXEx8fHxxzwnPj7+tMoDdOrUiejoaLZt23bM43a7nbCwsGabiIiItCx3HxzDxx9R2Ww2Bg0axIIFC9z7XC4XCxYsIDMz85jnZGZmNisPMH/+/OOWB9izZw8HDhwgISGhZQIXERGR01Yf3J5xdQ8x2Zju6VBa/xHV1KlTue6660hPT2fIkCE888wzVFdXM3HiRACuvfZa2rdvz4wZMwC4/fbbGTlyJE899RSXXnopc+fOZeXKlbz88ssAVFVV8cgjjzBu3Dji4+PZvn0799xzD2lpaWRnZ7f27YiIiMhxuPwCWGV2I8ri2XWo4CwkOFdeeSX79+/nwQcfpKioiP79+zNv3jx3R+KCggIs31txdNiwYbz55pvcf//93HfffXTp0oUPPviA3r17A2C1Wlm/fj3//Oc/KS8vJzExkYsvvpjHHnsMu93e2rcjIiIix+E6vFKDp5dpADBM83sLR5wjHA4H4eHhVFRUqD+OiIhIC8nbns/brz6FERDCAw/MbPH6T+f7u02OohIRERHvY63ay4P+/6LEGQW0fIJzOjzfzVlERER8w7m0VIOIiIicG1zm4QTH14eJi4iIyLnjnFmLSkRERM4d3y3V4Pn0wvMRiIiIiE8wzaYWHJ9fqkFERETOHabzcAuOF/TB0TBxERERaRFVEd35df10oiPCecbDsSjBERERkRbRaAtjkasP3fxDPR2KHlGJiIhIyziyVIPh+S44asERERGRluFftZerrAsJbogDRng0FrXgiIiISIsILP+Wmf5/5+rauZ4ORQmOiIiItIwjE/25DKuHI1GCIyIiIi3ku4n+PN8JRwmOiIiItAjT1EzGIiIi4muOJDheMIxKCY6IiIi0DPdim+qDIyIiIj7C9KIWHM2DIyIiIi2iNGoAN9bfQWK79vTycCxKcERERKRFHAqI47+uwZwXGO3pUPSISkRERFqGy2xaq8ELnlApwREREZGWEVhZwM8si+hRt97ToSjBERERkZbR7sBKnrX9lUsr3/F0KEpwREREpIW4NNGfiIiI+Jjvhol7Pr3wfAQiIiLiG8ymif5QgiMiIiK+wtQjKhEREfE5hx9RecM4cSU4IiIi0iLcLThe8IhKMxmLiIhIi9gdNZS36m+mc2R3Bnk4Fs+nWCIiIuITDgal8B/XeewM7ufpUJTgiIiISMtwNa3UgMWiPjgiIiLiI0KqdnKxZQUdar/1dChKcERERKRldNq/gJdtTzO87H1Ph6IER0RERFqIe5i459OLsxLBiy++SEpKCgEBAWRkZLB8+fITln/33Xfp3r07AQEB9OnTh08//bTZcdM0efDBB0lISCAwMJCsrCy2bt3amrcgIiIiJ3FOLdXw9ttvM3XqVB566CFWr15Nv379yM7OpqSk5JjllyxZwtVXX82kSZNYs2YNY8eOZezYsWzcuNFd5sknn+S5555j9uzZLFu2jODgYLKzs6mtrW3t2xEREZHjMLyoBccwTdNszQtkZGQwePBgXnjhBQBcLhdJSUnceuutTJs27ajyV155JdXV1Xz88cfufUOHDqV///7Mnj0b0zRJTEzkzjvv5K677gKgoqKCuLg45syZw1VXXXXSmBwOB+Hh4VRUVBAWFtZCdyoiInJuW/LKHQzb+yrLY8Yx5JZXW7z+0/n+btUUq76+nlWrVpGVlfXdBS0WsrKyyMnJOeY5OTk5zcoDZGdnu8vn5+dTVFTUrEx4eDgZGRnHrbOurg6Hw9FsExERkZblbjPxghacVo2gtLQUp9NJXFxcs/1xcXEUFRUd85yioqITlj/y7+nUOWPGDMLDw91bUlLSj7ofEREROQEvekTl+QjOgunTp1NRUeHedu/e7emQREREfM6WiBFMb5hEXruskxduZa2a4ERHR2O1WikuLm62v7i4mPj4+GOeEx8ff8LyR/49nTrtdjthYWHNNhEREWlZe4J68JZzFIVhfT0dSusmODabjUGDBrFgwQL3PpfLxYIFC8jMzDzmOZmZmc3KA8yfP99dPjU1lfj4+GZlHA4Hy5YtO26dIiIi0vqOdMGxesFSDa2+mvjUqVO57rrrSE9PZ8iQITzzzDNUV1czceJEAK699lrat2/PjBkzALj99tsZOXIkTz31FJdeeilz585l5cqVvPzyywAYhsGUKVP44x//SJcuXUhNTeWBBx4gMTGRsWPHtvbtiIiIyHFE1uzkPMsGImsDgG4ejaXVE5wrr7yS/fv38+CDD1JUVET//v2ZN2+eu5NwQUEBFst3DUnDhg3jzTff5P777+e+++6jS5cufPDBB/Tu3dtd5p577qG6upobbriB8vJyzjvvPObNm0dAQEBr346IiIgcx+DS97jd9m+Wl0wERnk0llafB8cbaR4cERGRlpfz/EQyD7zH8qRJDJn0lxav32vmwREREZFzx5GlGgyL59MLz0cgIiIiPsFwr0Vl9XAkSnBERESkhRims+lfTfQnIiIiPuNIt149ohIRERGf4V6qwfOPqFp9mLiIiIicG1aEXsj8A9H0jxrs6VDUgiMiIiItY3NgOv9wjqE8so+nQ1GCIyIiIi3DdbgPjmF4fqkGJTgiIiLSImJqCxhofEtwwwFPh6IER0RERFrG2IOv8p79YZKLFpy8cCtTgiMiIiIt48goKi9YTVwJjoiIiLQIgyNLNXh+mLgSHBEREWkRR5ZqMLxgHhwlOCIiItIiDPdMxkpwRERExEcYOA//oD44IiIi4iOOPKKyWDy/UILnIxARERGfsDDgYuZXpzE0orunQ1GCIyIiIi3jK/tI1jr70S+ym6dD0SMqERERaRnm4U7GXjANjhIcERERaRlxjXvoaezEv7HK06HoEZWIiIi0jCmVT9PTvoWN+yOBzh6NRS04IiIi0iIsh2cytmiiPxEREfEVR5ZqwOL59MLzEYiIiIhPsJhai0pERER8jEHTKColOCIiIuIzLGbTUg0WPaISERERX+FNLTgaJi4iIiIt4kPrRfjXlnJBSAdPh6IER0RERFrGW5bLKGqsZURER0+HokdUIiIi0jJch5dqMLxgqQa14IiIiEiLSHAVEWjU4eeq83QoSnBERESkZcx2PkSCvZRdB7tBh1iPxqJHVCIiItIi3Es1eMEoKiU4IiIi0iKOJDiG5sERERERX3FkqQaL1fM9YJTgiIiISIuwHJnoz/B8etGqEZSVlTF+/HjCwsKIiIhg0qRJVFVVnfCc2tpabrnlFtq1a0dISAjjxo2juLi4WRnDMI7a5s6d25q3IiIiIifh7oPj6y0448ePZ9OmTcyfP5+PP/6Yr7/+mhtuuOGE59xxxx189NFHvPvuu3z11Vfs27ePn//850eVe+211ygsLHRvY8eObaW7EBERkVNh9aI+OK2WYuXm5jJv3jxWrFhBeno6AM8//zxjxoxh1qxZJCYmHnVORUUF//jHP3jzzTe58MILgaZEpkePHixdupShQ4e6y0ZERBAfH99a4YuIiMhpesOZhY16LguM8HQordeCk5OTQ0REhDu5AcjKysJisbBs2bJjnrNq1SoaGhrIyspy7+vevTvJycnk5OQ0K3vLLbcQHR3NkCFDePXVVzEPz554LHV1dTgcjmabiIiItBzTNJnReDWPNF4HQdGeDqf1WnCKioqIjW0+yY+fnx9RUVEUFRUd9xybzUZERESz/XFxcc3OefTRR7nwwgsJCgrif//7HzfffDNVVVXcdtttx6x3xowZPPLII2d2QyIiInJc329nsFo8v1bDabfgTJs27ZidfL+/bdmypTVidXvggQcYPnw4AwYM4N577+Wee+7hz3/+83HLT58+nYqKCve2e/fuVo1PRETkXOM0TeI5QCwHsZhOT4dz+i04d955JxMmTDhhmU6dOhEfH09JSUmz/Y2NjZSVlR2370x8fDz19fWUl5c3a8UpLi4+YX+bjIwMHnvsMerq6rDb7Ucdt9vtx9wvIiIiLcNlmiwNuBWAyroLIKS9R+M57QQnJiaGmJiYk5bLzMykvLycVatWMWjQIAAWLlyIy+UiIyPjmOcMGjQIf39/FixYwLhx4wDIy8ujoKCAzMzM415r7dq1REZGKokRERHxEJfT5f7ZG5ZqaLU+OD169GD06NFcf/31zJ49m4aGBiZPnsxVV13lHkG1d+9eRo0axeuvv86QIUMIDw9n0qRJTJ06laioKMLCwrj11lvJzMx0j6D66KOPKC4uZujQoQQEBDB//nyeeOIJ7rrrrta6FRERETkJl7PR/bPF4vl5cFo1gjfeeIPJkyczatQoLBYL48aN47nnnnMfb2hoIC8vj5qaGve+p59+2l22rq6O7Oxs/vrXv7qP+/v78+KLL3LHHXdgmiZpaWn85S9/4frrr2/NWxEREZETcLq+63djsXq+BccwTzS+2kc5HA7Cw8OpqKggLCzM0+GIiIi0eRUVFYQ/nQxAw7278Q9s+e/X0/n+9vxUgyIiItLmfb8Fx+oFfXCU4IiIiMgZczm/S3AML0hwPN8LSERERNo8l2HlzcYLsBomV1r9PR2OEhwRERE5c06/QO5rvB4/i8GVXtCCo0dUIiIicsZch4csWbxgmQZQC46IiIi0AJfTSRhV2AzvSC28IwoRERFp04yqQtYH3EC96QeM83Q4ekQlIiIiZ87V2DSKyoV3PKJSgiMiIiJnzHV4HhyX4R2phXdEISIiIm2aaTYlOE4vSS28IwoRERFp08zDLTiml6QW3hGFiIiItGlHVhN3eUlq4R1RiIiISJvmcrqa/vWS1ELDxEVEROSMNfqH8L5zONhCuMLTwaAER0RERFpAbXAidzTcQoeQQK9IcLyjHUlERETaNNNsWqvB6iVLNSjBERERkTPmdLqw0YAfLk+HAugRlYiIiLSA4OIVfBtwHQWH2gObPR2OWnBERETkzJkaJi4iIiK+xnQ1PZoysXo4kiZKcEREROSMma7DLThai0pERER8xZG1qPSISkRERHyHU6uJi4iIiI9xuRfbVB8cERER8RG1tnb815nOFntvT4cCaB4cERERaQGlkf24tWEqQ8OjGOvpYFALjoiIiLQAp6tpqQY/i3ekFt4RhYiIiLRpjS7vWotKj6hERETkjCUVvM82+4Ns2J8BzPN0OGrBERERkTNnOhvxM1xYDNPToQBKcERERKQFmIeHiWNomLiIiIj4iCOLbZpKcERERMRnuJTgiIiIiK85stimxTvGLynBERERkTPnbADA9PUEp6ysjPHjxxMWFkZERASTJk2iqqrqhOe8/PLLnH/++YSFhWEYBuXl5S1Sr4iIiLSucnsCXzv7UBLQydOhAK2Y4IwfP55NmzYxf/58Pv74Y77++mtuuOGGE55TU1PD6NGjue+++1q0XhEREWldue0u4tqG6SyLv8bToQCtNNFfbm4u8+bNY8WKFaSnpwPw/PPPM2bMGGbNmkViYuIxz5syZQoAX375ZYvWKyIiIq3ru5mMvaP3S6tEkZOTQ0REhDsJAcjKysJisbBs2bKzXm9dXR0Oh6PZJiIiIi3HvRaV1TuWamiVBKeoqIjY2Nhm+/z8/IiKiqKoqOis1ztjxgzCw8PdW1JS0o+OQURERI42ctcLrLdP4rx9r3k6FOA0E5xp06ZhGMYJty1btrRWrD/a9OnTqaiocG+7d+/2dEgiIiI+xd9ZQ5hxCH+cng4FOM0+OHfeeScTJkw4YZlOnToRHx9PSUlJs/2NjY2UlZURHx9/2kEe8WPrtdvt2O32H31dEREROTHD9K5h4qcVRUxMDDExMSctl5mZSXl5OatWrWLQoEEALFy4EJfLRUZGxo+LtBXrFRERkTNjHJ7oD6t3JDit0genR48ejB49muuvv57ly5ezePFiJk+ezFVXXeUe6bR37166d+/O8uXL3ecVFRWxdu1atm3bBsCGDRtYu3YtZWVlp1yviIiInH3GkcU2Lf6eDeSwVhvL9cYbb9C9e3dGjRrFmDFjOO+883j55ZfdxxsaGsjLy6Ompsa9b/bs2QwYMIDrr78egBEjRjBgwAA+/PDDU65XREREzj6LebgFx0seURmmaZqeDuJsczgchIeHU1FRQVhYmKfDERERafPWzvop/au+Zkn3+xh21b2tco3T+f72jtl4REREpE0r8k9ipasrdYGxJy98FijBERERkTP2ftRv+UX9w+yLv9DToQBKcERERKQFuGcytvjwTMYiIiJybjkn1qISERGRc8vtxfezzH4zifsXeToUQAmOiIiItIBQ50HijHKsXrJUgxIcEREROWPWw/PgWP18fKI/EREROXdYzKaWG6u/zcORNFGCIyIiImfsSAuOxU8JjoiIiPiII0s1+CnB8T019Y0crKrFdLk8HYqIiMhZZeVwHxx/7+iD4x0rYvmI+ZuL+eidf/C8//OUWqKptEVTGxCHMyQea1gi9nbtCUzNJLZDJ0LseutFRMR37CCJMlcwVnuop0MBlOC0qANV9SQYBwg06kky90HdPqgDKoC9TWVuWXgbn7iGEhbgx0+DN3Nj/escCojBGRyPEZaALaojIfGdiGjfBXtUMli9IxMWERE5kVss91N2qJ7/RXfxdCiAEpwW9dvzUqkdNIM9+26gvHgnNft301C+F7OyCP+aYoLr9lNqS4BacNQ2YjTsItl/BzTsgEqgqHl9f7BMYUPURSSEB9Dffw+D63KwtkshMLYzke3TiI7viNVq9ci9ioiIfF99Y1P3DJvVO3q/KMFpYQGBQXTo3IsOnXsd8/jbQFVdI4Xlh9hf2JEv9mZQX7YHV2Uh/tVFhNYWEtNYRHv2s/lQJOv3VLB+TwVx1oWk+/8T8r+rq870Z68lhjL/eD6Pv56G+AEkhAeQFOwkMdRKbGwCUSF2DMM71gURERHfVe9sSnD8/ZTgnLNC7H50iQulS1xv6N/7qOOmaXKwuo7Hyg9R6KhnX/khbAWl5BTuJ+TQXqIaiohz7cduNJBs7iO5fh+Pbi1m9bc7APiN9X885j+HKjOAb4mlzD+B6qD2mBHJ+LdLJaDzebRPTCQxIhCrlyyKJiIibZdpmvzXcjvYwF67AEjydEhKcLyRYRhEhQQQFRKAO/0ZlgJc5y7jbGxg/758DhZu5VBxPmNDfkK/aj8Ky2vpua8eDkGIUUs3CqCxABw0bQXw86UPs9rsir/V4JrQtVxmyaEuNAlrZApB8Z1p1z6NuOSu+NsDz/q9i4hI29PodJFqKQbA4ecdXSeU4LRRVj9/YpK7EpPcFYB+zY4Ogoa/UH9gJwf3bqOiaDv1+/OxlO8iqGYvluAUbOUW6p0uEqs2Mtjva6gBioEtTTW4TINCSzTPxz6GJaE3Ke2C6RrooGMoxKd0w67kR0REDquvr+PIkBh/W4BHYzlCCY6v8g/EFt+DuPgexP3g0L8Bp8uksOIQZXkBrNrZC9fBndgqdxNWu49YZxHBRh0J5n4+3+WkZFcBANP83mKE30c0mhb2WGI5YO9ATUgKRHUmMKEr4d1G0iGuHf5e0sFMRETOjvraGoIP/2wLCPJoLEcowTlHWS0GHSKD6DB0FAwd1eyY6XKxv3gP+3flcpe1B7vKath5oIa4XRZqau0EGXV0MIvoUFsEtSuhFPgWMv/7PCWWaDpEBnJNwBL6WndhtEsjMKEb0ck9iE/qrFFfIiI+yFlXAzS1/lv97R6OpokSHDmKYbEQk5BMTEIyPZsdeR3T5aK0uID9OzdTvS8PZ+k27I58QmoLqfCPxtlgsutADR38vyTTuhxKgNyms2tNfwqsCRwMSOLzrg+TEB9HartgUsMhoV0UFrX8iIi0SfW1TQlOHf4EesnIXSU4cloMi4XohBSiE1KAMc2ObTJNih115JdW49x4Jcv3dcJeuZPI2t0kOIsIMBpIdRWQWF3IL5YW42I/AM/7P8eFljUUWRM5GJhEXVgqlug0QhK7EpfSi5i49hgWJT8iIt7KWX8IgDrDhrf00FSCIy3GMAziwwOIDw+Azr8Dfuc+5mpsoHD3Vg4U5FJeWshEe2d2llaTf6CajhUlBBt1dHblQ3U+VH8NhcAGaDCt9HH9iw7tQkmNDibbWEpCkEloYndiUroTHZOo5EdExMMaXCbfutpzyBJMhKeDOUwJjpwVFj9/ElJ7kpDa9NDrvO8da6zPYV/BtxzYnUtt0bdwYDuBVbuIrttNlelPVQNsKapkS1ElN9r+Tn/LDljXdK7DDKLYL5GKoGTqwtMoHHA7qdFBpLQLJirYpkkORUTOgqrQTlxe/2cSwgPI8XQwhynBEY/zswWQmNaXxLS+Rx1raKjni4oG8kuryC+toWpdBpscYUTX7yHOLCXMqCHMuQ0qt5FfsYnx2853n/tGwJ+I86umIjCZhvBU/GPSCG3fjfiUnoS1iwclPyIiLeJQvROAQH/vGUiiBEe8mr+/jdRoG6nRhwcgnjfbfazuUBXFu7ZQtnsLtcXb2F/jYpjRjp2l1eyrqKWHuZ2oxiqo3Na01tceYE3TuXl0ZFrsS6S2CyYlOpghjauIiY4lLrUnIZFxSn5ERE5DTcPhBMemBEfkjNkDQ0junk5y93T3vp8e/re2vpHCbf9h5+4t1BZvxXpwB8HVu4iu30s8pRQ7w1hTUM6agnIAVtjvJsZwAOAgmBK/RCqDO9IYkYoloR+BfX9GSnQQQTb9yYiI/FDIroXMtz3KrkN9gJ94OhxACY74qACbH6k900ntmX7UseqqSuKKinixJoSdB6opKCmnaFsqzoY9xHOAMKoJa9wKFVuhAhbt6MW4L6MAiAuzM9uYCQEROCNT8I/pQkSHbsSl9CIgLFotPyJybqoppYtlLzVmoqcjcVOCI+ec4JBQuqWF0q3Z3q8BcFQ6KMrfQvneXOqKt+FXvoO8hgQiDvlTXtNApaOCAQEroA6oAHYCK5pqcBDM8uDzWZg2nU7RwaS0C6abaztxyV2xhUWfzVsUETmrnPWHJ/qzesckf6AER6SZsNAwwvoOgb5D3PsygUlAeU09O4sOsDT3WepLmpKfkJoC4hr2EkcZYVRTUlHDm8ualrYIoI4tARMBcBDCflt7qoKTcUV0wh7XhbBO6cR17qelLUSkzTMPz4NjKsERaXsigmz075QAnSY022+aJgfKyynamUuUw2RybTT5pdXUFG+nqCKKeKOMMKoIq8+D+jw4COTDO4tGcp/r93SIDKRLO39urXwWV1RnAmLTiEzqQXTHHliDozxyryIip8PVUNv0r593LLQJSnBEzphhGLSLjKRd5DB6AaPdRwZimr+gpKyMop1bcOzNo2H/NvzL8wmrKSDX1ZnGRpOdB2rwK9tDX/v/mpKf7XBkIgkHoZTa27Ml/nLKevya1OhgUtoFkhBsxeIlK/aKiNDQ1IKDEhyRc4NhGMS2a0dsu+EwaHizY71dJjdU1pJfWk3RnngWbr0Z/4qm5Ce+cR9xxkHCqCSsbgvvbevDC3kbAehoFLHAdhe7re3ZH9iJmohuWOJ7Epnan5S0XgQH2DxxqyJyDjMPt+Dg5y0LNSjBEfEYi8UgITyQhPBA6BwNIwe5jzldJnv2H6B4Vy5Ve/MIakjkwpoYdpZW0/XgXvwMFx1du+lYvRuqv4K9wCqoMe08ZbuO3A5X0i0+hB4xNnpEmiQnpeDv5z3zU4iIb6kikD1mNI0B3vNYXQmOiBeyWgw6xEXTIe4nwE8Y+b1jjY0/Yc/eayjbsZa6fRvxP7CFiMqtJDTsIsioI7/azue5xXyeW8wIyzpet/2JMjOU3f4plId0wRXTg5DkviR1G0hcTIyWsxCRM/ZhxK/5/e4sHknrxVBPB3NYqyU4ZWVl3HrrrXz00UdYLBbGjRvHs88+S0hIyHHPefnll3nzzTdZvXo1lZWVHDx4kIiIiGZlUlJS2LVrV7N9M2bMYNq0aa1xGyJex8/PSoeOaXTomAb84rsDLicH9+RxbXUg6WWQV1xFUv5SnA6DKKOSqMYNUL4ByoGtwAKYZtzB3g6X0CsxnIHt6ukZ3khi575Y/PT/PiJy6mrqz6GZjMePH09hYSHz58+noaGBiRMncsMNN/Dmm28e95yamhpGjx7N6NGjmT59+nHLPfroo1x//fXu16GhoS0au0ibZLESmdyTIcB3g9z7YNZPp3jHOg7kr6V+70YCDuYRXbOdaLOMjbXt2Li1lG+2lnKt9b9c7P9Pqs0A8u1dKY/sg7VDOtHdh5GS2kWPuETkuI4kOEG+nuDk5uYyb948VqxYQXp600yyzz//PGPGjGHWrFkkJh57psMpU6YA8OWXX56w/tDQUOLj41syZBGfZdiCiOueSVz3zGb7ayv283g5bCo6xMZ9FXTaZqG6yk6wUUvv+vVQvB6K34BVUGJG8HjU40SmDqB/UgT9OoSTEh2sx1siAsDvy57kPtsuXOWPAt4xm3GrJDg5OTlERES4kxuArKwsLBYLy5Yt44orrjij+mfOnMljjz1GcnIy11xzDXfccQd+J2hSr6uro66uzv3a4XCc0fVFfEFAeAz9wqFfxyN7nqGx4c/kb11L2bc5sHcVUeUbSGrIJ5oKFhTaqSrcCcB9fm9wkd8a9gb1pD6uP6GdM0jtPZToiDBP3Y6IeFByQz5plnw2Gw2eDsWtVRKcoqIiYmNjm1/Iz4+oqCiKiorOqO7bbruNgQMHEhUVxZIlS5g+fTqFhYX85S9/Oe45M2bM4JFHHjmj64qcC/z8/UntOZjUnoPd+8z6agq3refx+mTW7a5g7e6DDCzeRir7SK3ZB/mfQz7UzfdjsyWVkvDe5Pe/l36d4umdGI7NTzM1i/g6u9k0TNw/0Hu6jJxWgjNt2jT+9Kc/nbBMbm7uGQV0MlOnTnX/3LdvX2w2GzfeeCMzZszAbj/2FNHTp09vdp7D4SApKalV4xTxFYYtmMSemfwM+Fn/9gDUOz4hf9NiHNtysBWvIbF6M+FU0tPcSszBEibM+wWwA7ufhcciPiYp3J/gThmk9DufsOgEj96PiLS8ALNpoj9bUBtNcO68804mTJhwwjKdOnUiPj6ekpKSZvsbGxspKytr8b4zGRkZNDY2snPnTrp163bMMna7/bjJj4icPltYDKmZYyFzbNMO06SyeBv7Ni5ib0kpWQ1xrNp1kIM19Zxf+TGxVeWw9zX4BgosHSiKGIglZRiJ/UaRkNxFfXlE2rgg8xAYYAvynsfUp5XgxMTEEBMTc9JymZmZlJeXs2rVKgYNapq8bOHChbhcLjIyMn5cpMexdu1aLBbLUY/EROQsMgxC47vQLb4L3YALaVqja3uxg12LbiV/zypiHRtJdRWQ7NpDctkeKPuQtSs7M9b+JINTo8hIjWJ4u2pSO/fAogVIRdoMp9NJsNHUzzWgrSY4p6pHjx6MHj2a66+/ntmzZ9PQ0MDkyZO56qqr3COo9u7dy6hRo3j99dcZMqRpUGtRURFFRUVs27YNgA0bNhAaGkpycjJRUVHk5OSwbNkyLrjgAkJDQ8nJyeGOO+7g17/+NZGRka1xKyLyIxmGQVp8OPziLve+0pIidq1dQP32RbQ7sJKldd0pqazjk/WFfLE+n/X231FuhJAf1I/69kNp1+t8OvXKwM/f34N3IiIncqimkiMz3AWG+HiCA/DGG28wefJkRo0a5Z7o77nnnnMfb2hoIC8vj5qaGve+2bNnN+sMPGLECABee+01JkyYgN1uZ+7cuTz88MPU1dWRmprKHXfc0ax/jYh4r+jYeKIvHg+MByCprpH+eytYkV/G/rwcGoutROEgquYb2PoNbP0zle8HsiOwN3s7X0Xs4J/Tp0M4ds3JI+I1DtVU4zCjCKKO8IBgT4fjZpimaXo6iLPN4XAQHh5ORUUFYWHek22KnOsa6g6Rv34RB3O/IrBwGamHNhBKU+fFBxom8C/nxdj9LFycWMt429eEdBtJ5wEXEBgS7uHIRc5dO/ZXceFTXxEa4MeGh7Nb9Vqn8/2t+dhFxGv42wPpOvgiGHwRAM7GRnZsXk7Jpi+gti/t9tg4UF1P8N5FDPV/Ffa8SuPnFrb6p3EwOp3AtPNIGTiK0ChNBCpytlTWNgIQFuBdj5KV4IiI17L6+dGp7zA69R3GUOBR02T7/mp2rjzEii276OBYS4Kxny6N30LRt1D0JiyC+8IeJ6DLhQxJjWJIcihRYd7TbC7iaxy1TZP7hQZ4V0rhXdGIiJyAYRikxYaQNuYXMOYXmKbJvoKt7F67EHPnYuLK15Ds2sPHJbE4SvJ5dXE+d/vN5Zf+SyiMGAAdh9G+7wVEp/QDi0ZqibQE+84veN/2J/bU9wNGeDocNyU4ItJmGYZBYseuJHbsCvwegOLiYv5Y5GR5/gGW55cx6OBWYs1SYg/Oh4PzYe0jVBDKntA+NHYYStQFt9IhJlJz8Yj8SIZjLwMs2zDMaE+H0owSHBHxKXFxcVweB5f3a5qSouzgZ6xY8yXVWxcRUbqCrvVbCDcqCa9cwsHN6xm4ZiDx4UEMSY1iXMBKOiclkth7BIbde2ZkFfFmzkMVADT4e9ffjBIcEfFpUZGRRF14BVzYtMivo7qGlesW49jyNfsOVGBtsFJYUct/1u7lXvsTJK4to/EjC3vsXaiKG0xYtxG0738x1mDNtSVyTLVNCY7LpgRHRMRjwoKDSB92EQxrGqk1rt7JmoKDrNy2j+3rBuCqWUcHo5SUujwoyIOC/8M132Bp8AWsGzKLIalR9G4fjr9mWxYBwKhzAOCye9d0DUpwROScFmizMiwtmmFp0TD6PeoanazLy6V4wxf47VlKx8o1dDb2sqnCzozPtgAQ7u/k9aBnqU7IIKL3RXTpdx7+mm1ZzlGWwwmOJUAJjoiI17L7WenXqzf06g3cSqPTRe62rQQWlHFRoT8rdpbRo3YT/epWws6VsPNFKj4KZmtQf2qTRhDX/2I6d+uv9bTknGFrKAfA4mWPcZXgiIicgJ/VQo9u3ejRDa4BXC6THfmdWLbShb3gGzpXrybcqCb90GL4djF8O4OZxm8pSPs1wzpHM6xTFKkxIRqlJT7rkNNKpRmIX4hGUYmItFkWi0Fa5zTSOk8DpuFqbGDHxiUc2PA/gvcuIu3QRr6uS2PzhiI+3VDEFZZvuNX+MYURg/DrdB6pAy8mtn1HT9+GSIuZbptGvqOad9KGejqUZpTgiIicAYufP536j6RT/5EANNRW81jRIRZvP8jibaWM2LuJTuZuOh3cDas+gFWw20ikOHIQts7nkTz8V0RERHn2JkTOQMWhppmMI4JtHo6kOSU4IiItyD8gmEEpwQxKiea2UV04VPEam1fOo3rr17Tbv4KUxnyS2EdS2T4o+4iBi8KJj+/AkNQoRkUU0aNjAtHJPUCPtKQNcLlMymvqAYgI9K6O9kpwRERaUWB4ND1H/RpG/RqAirJStq+aT83Wr6kv20NZbRhlhQ42Fzq42P+PRFs3U2pEsjt0AI1JmcT1uZCkrgMwLFYP34nI0aocpbznfz8HzVDCAi72dDjNKMERETmLwqOiGXjR1XDR1QAsd9SyfGcZK3YcwLYxgLoGP6I5SLRjIWxaCJse5yChbAoZzuYhTzCoYxS924dh91PCI55XfaCI/pYdVJmBBNjUgiMiIofFhgVwWd9ELuubCGO/wFFVSe6ar6j+9mvCipeRVreZSKOS+ooinvi0aR4em5+F/wt+FiI6Yk8dSlK/C4hKSPHsjcg5qbqiFACHEUqIh2P5ISU4IiJeJCwklP4/uQx+chkA9XW15G1YTGWhg4sOxrF610H8qosYUpcDxTlQPBeWQrERzb7QvrjaDyaqz0V07J6OxaJ+PNK6ah37Aaiyhnk4kqMpwRER8WI2ewDd0kfRDfgZYJomuwqLWb5iJq6CpUQfXEeqcydxlBLnWAiOhczZsIan/H7HwORIhnQIZJT/BpL6nU9wVKKnb0d8TENlUwvOISU4IiJyJgzDICUxnpSf3QTcBEBF+UHy139D1bbFBBWvZpWrL5W1jXz17X5qtm7hFvuj8BUUWuIpCe8LSRnE9/oJsWkDMaze1W9C2hZndRkAdTbvWqYBlOCIiLR54RGR9B9xOYy4HIC+Thc3FlWyatdBqjfvYfveZFJdu0lwFZFwsAgO/g/WQw0BvBV/D/T+OYM6RtIzPhSbvzovy6lzHU5wGm3etUwDKMEREfE5flYLvduH07t9OAy7EbiRopJidq37mkM7lhBeupa0+lxCjUN8tMuPtTs3A/Ar/0VMsf2HA1H9sXYcTmL6ZUTEa9ZlOb5D9Y1UmoEQ5H2TVSrBERE5B8THxhF/0S+BXwJwqLaedZtWkl0eSbvdVawqOEi/+jwSnXtJ3L8X9n8CK+9juyWFXVHDMNKy6Dr4Itq3876+FuI5c0Ov47pd2TzSoyfDPB3MDyjBERE5BwUG2Og3aBj9Dr92uUx27unB1xu+pmFnDgkHltHduZXOrp10Lt0JpW+S8eULBMckMaJLDCNSQxnSNZEQu75GzmX7K+sAiAkL8HAkR9NvpoiIYLEYdEpOolPyeGA8AOX7C9mz6hPYNh9XRSH766Io3l/Njv3VXLBiJvuMMraEZODsnEXqwFH06RiLVUPTzyn7qw4nOKF2D0dyNMM0TdPTQZxtDoeD8PBwKioqCAtTc6uIyKmoONRAzvZSluQVMm3DJQRR6z5WbdpZbvRhf9xPiOg3hiH9+xMR5F2LL0rLW/VQBuWuQLr+7jWSOqa2+vVO5/tbCY4SHBGR01dTxv5183Bs+Izo4m8Idx50H/rG2ZvrGu8jvWMUF/aIZVSXCNISojC0gKhPqa6sIPipZACqpuYTEtb6HY1P5/tbj6hEROT0BUURk3kNMZnXgMtF4771FK3+CHPrfDY2DMFVDst3lpG/czu/WTiVHGtvDiaOJHrgZfTv219rafmA8uJdBANVZiDBoRomLiIivsZiwa9Dfzp06A88wE3AZWU1fJFXQs3Ktwg+UMcw1yrYswr2/IWd/0lgZ+QwgnqOptewSwgOCfXwDciPUbW/AIADlnaEeGHrnBIcERFpcUlRQVybmQJDp3Fo7xj2LPsQ647PSa7eQIpRSEr5/4Ml/487v5lMRdoVZPeKI6tHHJHB6rfTVtSV7QGgwj/aw5EcmxIcERFpPYZBYId+dOnQD3gA16EKtq/4FMeGz0goXcIXDb0pyy3m89xiJvj9j58Hb6Aq7TI6/+Qq4uK1dpY3qz+4F4Bqe5yHIzk2JTgiInLWWALD6TziahhxNabLxZslVczbWMR/NxXz0wOL6Vu3FTatonHjY6yx96eq80/pOlLJjjcyK5oSHFdIvIcjOTYlOCIi4hGGxUL3+DC6x4cxJasr+7bPYcXit2hX8CmdGncwoH415K6mYfMfWREwhNyRL3FJn0SvnHPlXFRd17RMgyW8vadDOSYNE9cwcRERr3Ng12Z2ffMmkTs/IbVxB/9zDuKGhjuxGDC0Uztuit1E7+GXEhntna0H54KfPr+IDXsreOU3A7moV8JZuaaGiYuISJvWrmNP2nX8I/BHSnZuojJ3L/13BLN2dzl7dmzmJ3un0rD6btYHDqCu6+V0vWA84ZHe2dnVV+0tPwRA+8hgD0dybEpwRETEq8Wm9GJcSi/GAbvLali5aB471qXSyZlP39qVsH4ldeseY1XoMIx+V9FrxM+x271vbSRfUlPfSFl1PQDtIwM9HM2xWVqz8rKyMsaPH09YWBgRERFMmjSJqqqqE5a/9dZb6datG4GBgSQnJ3PbbbdRUVHRrFxBQQGXXnopQUFBxMbGcvfdd9PY2NiatyIiIl4gKSqIKy7/OZ0eWMvu8d+Qk3ITOy3J2I0GBlV9xcDFN3H/jJn84f0NrNxZxjnYC+Os2J+/gXm2e3ne/hLhgf6eDueYWrUFZ/z48RQWFjJ//nwaGhqYOHEiN9xwA2+++eYxy+/bt499+/Yxa9Ysevbsya5du/j973/Pvn37+Pe//w2A0+nk0ksvJT4+niVLllBYWMi1116Lv78/TzzxRGvejoiIeJGkLn1J6tIXzBns2JhD6aJ/Elf8NR/W9qNuWQFvLCvg92GLOS/eSdLIiXTs3N3TIfsMx948+lh242f1zuQGWrGTcW5uLj179mTFihWkp6cDMG/ePMaMGcOePXtITDy1IX/vvvsuv/71r6mursbPz4/PPvuMyy67jH379hEX1zT2fvbs2dx7773s378fm+3kk0Spk7GIiG9yOl0szS/j/TV7+WzDPj7kDjpbCgHY6Neb8i4/p0fWtbRrF+PhSNu25XMfZ8iWJ1kV9BMG3fPxWbvu6Xx/t9ojqpycHCIiItzJDUBWVhYWi4Vly5adcj1HbsLPz89db58+fdzJDUB2djYOh4NNmza13A2IiEibY7VaGJ4Wzaxf9mPlfaNwpN/K5oABuEyD3o0bOS/3UUKe68HKP1/O6oX/psHp8nTIbZKrbBcA9aEdPBzJ8bXaI6qioiJiY2ObX8zPj6ioKIqKik6pjtLSUh577DFuuOGGZvV+P7kB3K+PV29dXR11dXXu1w6H45SuLyIibVdggI0Bl98Cl99C2b4d7Fg4h5j8D+jo3EV69Ve894WTG5aE8/OBHfjloA50idOaWKfKXrUbAEtUimcDOYHTbsGZNm0ahmGccNuyZcsZB+ZwOLj00kvp2bMnDz/88BnVNWPGDMLDw91bUlLSGccnIiJtR1RiJ9J//Sgd719HwS/msSzuKj6zXUxpVT0vf72DW555g7w/DmHpv5/F4Sj3dLheL7x2HwCBsZ09HMnxnXYLzp133smECRNOWKZTp07Ex8dTUlLSbH9jYyNlZWXEx594YqbKykpGjx5NaGgo77//Pv7+33Viio+PZ/ny5c3KFxcXu48dy/Tp05k6dar7tcPhUJIjInIuMgySe2eS3DuTgU4XX+bt552Vuxm+9XW6NebBxgep2jCDxZEXETL8d/QZNAKLxftWyvYk0+Ui1lkMBkR16OLpcI7rtBOcmJgYYmJO3jkrMzOT8vJyVq1axaBBgwBYuHAhLpeLjIyM457ncDjIzs7Gbrfz4YcfEhDQfC6DzMxMHn/8cUpKStyPwObPn09YWBg9e/Y8Zp12ux27XVN7i4jId/ytFi7qGcdFPeMoLZrF8vldab/jHdpTxPDyD+GTD9nyWWeK0q6iz5jf0y5Cg1IA9pcfZL8ZSypFxCZ19XQ4x9WqSzVccsklFBcXM3v2bPcw8fT0dPcw8b179zJq1Chef/11hgwZgsPh4OKLL6ampob333+f4ODvZkeMiYnBarXidDrp378/iYmJPPnkkxQVFfGb3/yG3/3ud6c8TFyjqERE5FhMl5Nty+dxaOk/6HHwS/wNJ6VmGCMb/0pWnw6Mz+jI4JRIDOPcbdVZubOMX8zOoUO4nUXTs87qtb1mqYY33niDyZMnM2rUKCwWC+PGjeO5555zH29oaCAvL4+amhoAVq9e7R5hlZaW1qyu/Px8UlJSsFqtfPzxx9x0001kZmYSHBzMddddx6OPPtqatyIiIucAw2Kly9BLYeilHDpYzKb/vcyKAgfVByz8Z+0+Ply7hzdDnsNIy6Ln6N8RFtHO0yGfdbsONH1nJ0eHeDiSE9Nim2rBERGRk9iwp4L/W7qL0nWf8g/rDABqTDvroscQN/puOnXp5eEIz56/zP+W5xZs5eohScz4ed+zeu3T+f5WgqMER0RETpGjrITc/75M3Na5pLiahko3mhaWB5+PfeRUBg45z+cfX616aixBFdvZ3u9uLht37Vm9tldM9CciIuJrwqJiybj6fjrev57ci/6PTYHp+BkuhtUsZNBnlzH5qVeZu7yA2ganp0NtNXHV39LDUkBiO++eN0gtOGrBERGRM1C4ZRkH/vsnzLKd/LTuUcCgXbCN2/vBz0aNIDzYd0bxNtTXYjyegJ/hYt9vV5OYfHbnwfGaTsYiIiK+LqF7Bgnd38NRXc0fVhUzZ8lOHOUHGLvqVkpWtWNxl0kMvfxGosKCT16Zl9u3YzMdDRfVZgDx7VM9Hc4J6RGViIhICwgLDub6EZ346u7zmX2RP1bDIM3Yw5htj1D3VB/mv/owJWUHPB3mGSnNXw/APv8kLFbvTiG8OzoREZE2xs9qYfiosQTek0ten7soMyJJMA5wUcHT+D/bj29euQvHgZKTV+SFbDs+B2B/mPePGlOCIyIi0gosQRF0G/cAkfflsiX9MQqtCUQalQzb83eue/4j/vbV9rbVGbmxntTSLwCoSbvcw8GcnPrgiIiItCLDP5Dul92GecnNbPj8Xyxbs4Y15fGs+WwLc5bs5LH0Wi64cDRWq9XToZ5QbY2DD52Z9OdbUgaO8nQ4J6UER0RE5CwwrH70yZ5Iz4smEL56D0/P/5YgxzYuWHQvO5Z2ov7CR+k1/DJPh3lci/Y4ua9+Iu3DA1gUF+7pcE5Kj6hERETOIqvF4JfpSSy863zuHQSHjAC6OLfTa/541j05mr3b1nk6xGN6Z2XTxIYX945vE5MZKsERERHxgAB/Kxf94kYabl7F0uif02ha6FeTQ+y/LmD57BupLN/v6RDdSr5dQfmWrzBwMT4j2dPhnBIlOCIiIh4UGdueoZNfY/dVC1kTkIG/4WRI0VzKnzmPucvycbo8Px9v0aczecf2KE/FfEparHfPYHyEEhwREREvkNpjAP3v/S9rR/6DnUYH5jSMYtr7m7n0uW9Ysq3UY3EV7sqj+8Gm0VNp5/3SY3GcLiU4IiIiXsIwDPpf8AsSp6+hQ/YUwgL82FJUySuvzmbdk9ns2Xr2++fs+eARbIaTTbZ+9M248Kxf/8dSgiMiIuJlbDYbE0d05au7L+C6ocnc6zeXfjVLifu/C1j61+upOFB8VuLYu30DA8o+A8B60YNn5ZotRQmOiIiIl4oMtvHI2D4EXPN/rA4cir/hZGjJO/D8AHLe/CN1dbWtd3HTpOS96fgZLtYGZNB9cFbrXasVKMERERHxcind+zPw3v+y7sJ/km/pSDjVZH77Z4pmDmTpwv9gmi3fEXnTpy8xoPobGk0LIZc+0uL1tzYlOCIiIm1EvxFjSb5vFSt6P0AZYXQ09zL7841c8dclfLGlpMUSnZ2l1Ty6wmCHK55vOlxPWp/MFqn3bDLM1kj7vJzD4SA8PJyKigrCwsI8HY6IiMhpq64oY/kHz3PTtsHUNjR9ld8T9TUju8XSY/RNWOxBP6re3WU1/Pofy9h1oIbhSQG8dv0IbDb/lgz9Rzud728lOEpwRESkDSuprOXv3+Tz3tItfG7cQoRRzQEi+DZ1PF2yfkt0+7RTq8jlYud/HuPtDZW8VHMBSVGB/L/fDyM2LKB1b+A0KME5CSU4IiLia8oqHKz+4Hl65L9Ge76bBXmHfxcOJl9MXObVtO/c+6hlFsxD5exa8SmHlv6DHjUrqTet3BD6En+64WfEeVFyA0pwTkoJjoiI+KpDh2rZ8N9/ELzpTbrXb8JqNH3NP9v4c/5pv4ZeiWEM8Cvg4rJ/4V9fQedDG/DDCUCt6c9nyXeR/Zu7CLJ533rcp/P97X3Ri4iIyI8WGBjAkLG3wNhb2LungO2L3yVkx2d83pBBWXU932wtJcKylt62r9zn5JsJ7IoaTvJFv+eKnoM9GH3LUQuOWnBEROQcUNfoJLewkm+LKqkv2UrigSXYbTbsXS+kR6/+BNu9v81DLTgiIiLSjN3PSv+kCPonRQBJQNtZduHH0Dw4IiIi4nOU4IiIiIjPUYIjIiIiPkcJjoiIiPgcJTgiIiLic5TgiIiIiM9RgiMiIiI+RwmOiIiI+BwlOCIiIuJzlOCIiIiIz1GCIyIiIj5HCY6IiIj4HCU4IiIi4nPOydXETdMEmpZdFxERkbbhyPf2ke/xEzknE5zKykoAkpKSPByJiIiInK7KykrCw8NPWMYwTyUN8jEul4t9+/YRGhqKYRgtWrfD4SApKYndu3cTFhbWonXLqdFn4Hn6DDxPn4Hn6TNoeaZpUllZSWJiIhbLiXvZnJMtOBaLhQ4dOrTqNcLCwvQL7WH6DDxPn4Hn6TPwPH0GLetkLTdHqJOxiIiI+BwlOCIiIuJzlOC0MLvdzkMPPYTdbvd0KOcsfQaep8/A8/QZeJ4+A886JzsZi4iIiG9TC46IiIj4HCU4IiIi4nOU4IiIiIjPUYIjIiIiPkcJTgt68cUXSUlJISAggIyMDJYvX+7pkLzS119/zU9/+lMSExMxDIMPPvig2XHTNHnwwQdJSEggMDCQrKwstm7d2qxMWVkZ48ePJywsjIiICCZNmkRVVVWzMuvXr+cnP/kJAQEBJCUl8eSTTx4Vy7vvvkv37t0JCAigT58+fPrpp6cdS1s0Y8YMBg8eTGhoKLGxsYwdO5a8vLxmZWpra7nlllto164dISEhjBs3juLi4mZlCgoKuPTSSwkKCiI2Npa7776bxsbGZmW+/PJLBg4ciN1uJy0tjTlz5hwVz8n+dk4llrbmpZdeom/fvu5J4DIzM/nss8/cx/X+n10zZ87EMAymTJni3qfPoI0zpUXMnTvXtNls5quvvmpu2rTJvP76682IiAizuLjY06F5nU8//dT8wx/+YL733nsmYL7//vvNjs+cOdMMDw83P/jgA3PdunXm5ZdfbqamppqHDh1ylxk9erTZr18/c+nSpeY333xjpqWlmVdffbX7eEVFhRkXF2eOHz/e3Lhxo/nWW2+ZgYGB5t/+9jd3mcWLF5tWq9V88sknzc2bN5v333+/6e/vb27YsOG0YmmLsrOzzddee83cuHGjuXbtWnPMmDFmcnKyWVVV5S7z+9//3kxKSjIXLFhgrly50hw6dKg5bNgw9/HGxkazd+/eZlZWlrlmzRrz008/NaOjo83p06e7y+zYscMMCgoyp06dam7evNl8/vnnTavVas6bN89d5lT+dk4WS1v04Ycfmp988on57bffmnl5eeZ9991n+vv7mxs3bjRNU+//2bR8+XIzJSXF7Nu3r3n77be79+szaNuU4LSQIUOGmLfccov7tdPpNBMTE80ZM2Z4MCrv98MEx+VymfHx8eaf//xn977y8nLTbrebb731lmmaprl582YTMFesWOEu89lnn5mGYZh79+41TdM0//rXv5qRkZFmXV2du8y9995rduvWzf36V7/6lXnppZc2iycjI8O88cYbTzkWX1FSUmIC5ldffWWaZtN9+vv7m++++667TG5urgmYOTk5pmk2JaoWi8UsKipyl3nppZfMsLAw9/t+zz33mL169Wp2rSuvvNLMzs52vz7Z386pxOIrIiMjzb///e96/8+iyspKs0uXLub8+fPNkSNHuhMcfQZtnx5RtYD6+npWrVpFVlaWe5/FYiErK4ucnBwPRtb25OfnU1RU1Oy9DA8PJyMjw/1e5uTkEBERQXp6urtMVlYWFouFZcuWucuMGDECm83mLpOdnU1eXh4HDx50l/n+dY6UOXKdU4nFV1RUVAAQFRUFwKpVq2hoaGh27927dyc5ObnZ59CnTx/i4uLcZbKzs3E4HGzatMld5kTv8an87ZxKLG2d0+lk7ty5VFdXk5mZqff/LLrlllu49NJLj3qf9Bm0fefkYpstrbS0FKfT2eyXHCAuLo4tW7Z4KKq2qaioCOCY7+WRY0VFRcTGxjY77ufnR1RUVLMyqampR9Vx5FhkZCRFRUUnvc7JYvEFLpeLKVOmMHz4cHr37g003bvNZiMiIqJZ2R++P8d6b44cO1EZh8PBoUOHOHjw4En/dk4llrZqw4YNZGZmUltbS0hICO+//z49e/Zk7dq1ev/Pgrlz57J69WpWrFhx1DH9DbR9SnBEznG33HILGzduZNGiRZ4O5ZzTrVs31q5dS0VFBf/+97+57rrr+Oqrrzwd1jlh9+7d3H777cyfP5+AgABPhyOtQI+oWkB0dDRWq/WoHu3FxcXEx8d7KKq26cj7daL3Mj4+npKSkmbHGxsbKSsra1bmWHV8/xrHK/P94yeLpa2bPHkyH3/8MV988QUdOnRw74+Pj6e+vp7y8vJm5X/4/vzY9zgsLIzAwMBT+ts5lVjaKpvNRlpaGoMGDWLGjBn069ePZ599Vu//WbBq1SpKSkoYOHAgfn5++Pn58dVXX/Hcc8/h5+dHXFycPoM2TglOC7DZbAwaNIgFCxa497lcLhYsWEBmZqYHI2t7UlNTiY+Pb/ZeOhwOli1b5n4vMzMzKS8vZ9WqVe4yCxcuxOVykZGR4S7z9ddf09DQ4C4zf/58unXrRmRkpLvM969zpMyR65xKLG2VaZpMnjyZ999/n4ULFx71OG/QoEH4+/s3u/e8vDwKCgqafQ4bNmxolmzOnz+fsLAwevbs6S5zovf4VP52TiUWX+Fyuairq9P7fxaMGjWKDRs2sHbtWveWnp7O+PHj3T/rM2jjPN3L2VfMnTvXtNvt5pw5c8zNmzebN9xwgxkREdGsd700qaysNNesWWOuWbPGBMy//OUv5po1a8xdu3aZptk0NDsiIsL8z3/+Y65fv9782c9+dsxh4gMGDDCXLVtmLlq0yOzSpUuzYeLl5eVmXFyc+Zvf/MbcuHGjOXfuXDMoKOioYeJ+fn7mrFmzzNzcXPOhhx465jDxk8XSFt10001meHi4+eWXX5qFhYXuraamxl3m97//vZmcnGwuXLjQXLlypZmZmWlmZma6jx8ZInvxxReba9euNefNm2fGxMQcc4js3Xffbebm5povvvjiMYfInuxv52SxtEXTpk0zv/rqKzM/P99cv369OW3aNNMwDPN///ufaZp6/z3h+6OoTFOfQVunBKcFPf/882ZycrJps9nMIUOGmEuXLvV0SF7piy++MIGjtuuuu840zabh2Q888IAZFxdn2u12c9SoUWZeXl6zOg4cOGBeffXVZkhIiBkWFmZOnDjRrKysbFZm3bp15nnnnWfa7Xazffv25syZM4+K5Z133jG7du1q2mw2s1evXuYnn3zS7PipxNIWHev9B8zXXnvNXebQoUPmzTffbEZGRppBQUHmFVdcYRYWFjarZ+fOneYll1xiBgYGmtHR0eadd95pNjQ0NCvzxRdfmP379zdtNpvZqVOnZtc44mR/O6cSS1vz29/+1uzYsaNps9nMmJgYc9SoUe7kxjT1/nvCDxMcfQZtm2GapumZtiMRERGR1qE+OCIiIuJzlOCIiIiIz1GCIyIiIj5HCY6IiIj4HCU4IiIi4nOU4IiIiIjPUYIjIiIiPkcJjoiIiPgcJTgiIiLic5TgiIiIiM9RgiMiIiI+RwmOiIiI+Jz/D4gzlc5YcScZAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# get new values\n", "t_new = np.linspace(t[0], t[-1], 1000)\n", "\n", "# notice the new shape is (ninterps, length) just like the inital input to the spline\n", "new_teuk_modes = spline(t_new)\n", "\n", "# (220) mode (real part)\n", "ind = amp.special_index_map[(2, 2, 0)]\n", "\n", "plt.plot(t_new, new_teuk_modes[ind], label=\"interpolated\")\n", "plt.plot(t, teuk_modes[:, ind].real, \"--\", label=\"original\")\n", "plt.legend()\n", "plt.show()" ] }, { "attachments": {}, "cell_type": "markdown", "id": "88352b10-edac-456c-a390-111716f38fee", "metadata": {}, "source": [ "To get the array of interp coefficients for CUDA implementations, do the following. The underlying shape of the array is (4, length, ninterps). It is flattened though for entry into GPU kernel. " ] }, { "cell_type": "code", "execution_count": 20, "id": "cc04fd9e-60ee-4d0f-b90b-e45d08b2cb26", "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "array([ 6.23856751e-06, 3.56609646e-06, -1.26237401e-06, ...,\n", " 3.36335764e-11, 3.08981746e-11, 8.40344230e-12])" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "spline.interp_array" ] } ], "metadata": { "kernelspec": { "display_name": "few-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.12.3" } }, "nbformat": 4, "nbformat_minor": 5 }