{ "cells": [ { "cell_type": "markdown", "id": "828d581c-2c0b-4950-a1dd-0cf96f55f648", "metadata": {}, "source": [ "# K-D tree fundamentals\n", "\n", "TODO: what is a K-D tree\n", "\n", "The [nanoflann](https://github.com/jlblancoc/nanoflann) library is used within teqpflsh due to its computational efficiency." ] }, { "cell_type": "raw", "id": "a4329548-602c-4174-b3f6-bb7cb4373e5a", "metadata": { "editable": true, "raw_mimetype": "text/restructuredtext", "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "The :py:class:`L2TreeHolder ` class makes a copy of the data for the tree to ensure that the lifetime of the data copied into the holder is longer than the tree itself. The underlying :py:class:`L2Tree ` object obtained via the ``.tree`` attribute then makes a reference to the data held in the holder class. " ] }, { "cell_type": "code", "execution_count": 1, "id": "53459aa8", "metadata": { "execution": { "iopub.execute_input": "2025-01-06T11:32:10.863718Z", "iopub.status.busy": "2025-01-06T11:32:10.863422Z", "iopub.status.idle": "2025-01-06T11:32:11.221441Z", "shell.execute_reply": "2025-01-06T11:32:11.221201Z" } }, "outputs": [], "source": [ "import numpy as np\n", "import teqpflsh\n", "import matplotlib.pyplot as plt\n", "\n", "def boxpoly(*, top, bottom, left, right, ptr):\n", " X = np.array([left, right, right, left, left])\n", " Y = np.array([bottom, bottom, top, top, bottom])\n", " return ptr.makeclosedpolygon(X, Y)" ] }, { "cell_type": "code", "execution_count": 2, "id": "19ebf7aa-32b1-48a2-af0e-468ec2f0d83c", "metadata": { "execution": { "iopub.execute_input": "2025-01-06T11:32:11.222931Z", "iopub.status.busy": "2025-01-06T11:32:11.222806Z", "iopub.status.idle": "2025-01-06T11:32:11.230280Z", "shell.execute_reply": "2025-01-06T11:32:11.230050Z" } }, "outputs": [], "source": [ "ptr = teqpflsh.GeometryFactoryHolder()\n", "\n", "# Polygon for the shifted circle\n", "t = np.linspace(0, 2*np.pi, 10000)\n", "X = 0.5 + 0.3*np.cos(t)\n", "Y = 0.3*np.sin(t)\n", "poly1 = ptr.makeclosedpolygon(X, Y)\n", "poly2 = boxpoly(left=0, right=1, bottom=0, top=1, ptr=ptr)\n", "\n", "# Polygon for the square [0,1]x[0,1] minus small circle\n", "poly = poly2.difference(poly1)\n", "X, Y = poly.getXY()\n", "\n", "def do_one(*, NKD, Nsample, plot=False, close=True):\n", " def get_random(NKD):\n", " \"\"\" Random points for the tree \"\"\"\n", " XX, YY = [], []\n", " while len(XX) < NKD:\n", " x_, y_ = np.random.random(2)\n", " pt = ptr.createPoint(float(x_), float(y_))\n", " if poly.containsPoint(pt):\n", " XX.append(x_)\n", " YY.append(y_)\n", " return XX, YY\n", " XX, YY = get_random(NKD)\n", " \n", " if plot:\n", " plt.plot(X, Y, 'k')\n", " plt.plot(XX, YY, '.', ms=5)\n", "\n", " holder = teqpflsh.L2TreeHolder(np.array(XX), np.array(YY), 10)\n", " tree = holder.tree\n", "\n", " xsample, ysample = get_random(Nsample)\n", " d2 = [tree.get_nearest_indexd2(x_, y_)[1] for x_, y_ in zip(xsample, ysample)]\n", " \n", " if plot:\n", " plt.axis('off')\n", " plt.axis('equal');\n", " if close:\n", " plt.close()\n", " \n", " return np.mean(np.array(d2)**0.5), tree.get_used_bytes()" ] }, { "cell_type": "markdown", "id": "ec39b27b-1ab0-43f0-9865-736e54b21e09", "metadata": {}, "source": [ "Here is a small number of \"lighthouse\" points randomly distributed in the domain. A random point is first pulled from [0,1]x[0,1] and checked whether it is within the domain or not. This so-called point-in-polygon problem is quite slow (relatively)." ] }, { "cell_type": "code", "execution_count": 3, "id": "f37f0986-2884-419d-bc03-944bc1eae624", "metadata": { "execution": { "iopub.execute_input": "2025-01-06T11:32:11.231535Z", "iopub.status.busy": "2025-01-06T11:32:11.231429Z", "iopub.status.idle": "2025-01-06T11:32:11.286869Z", "shell.execute_reply": "2025-01-06T11:32:11.286648Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAGFCAYAAABg2vAPAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIzhJREFUeJzt3QmQVcXVwPHDzLDMIIQouwRE2REZN1wgCClBYkUjRsElhkiinzERFbU00RA/qkxF1KhINJpIQhlNXKKlIkHUBEVBcCFgUEDBIIgoKAg67DNfneZ748DAMPPm3u7b3f9f1dR9EyNcZ+7rd/r06dMNKioqKgQAAESrwPUNAAAAtwgGAACIHMEAAACRIxgAACByBAMAAESOYAAAgMgRDAAAEDmCAQAAIkcwAABA5AgGAACIHMEAAACRIxgAACByBAMAAESOYAAAgMgRDAAAEDmCAQAAIkcwAABA5AgGAACIHMEAAACRIxgAACByBAMAAESOYAAAgMgRDAAAEDmCAQAAIlfk+gZiMGbMGHnuuedc3wYAeGno0KFy5513ur6NoDWoqKiocH0TISsrK5OmTZu6vg14rqCkhbQ89Qpp1L6bbFu9VNZNu0PKyza4vi3Ami+//FJKSkpc30awCAYsPMAHHHCAeT19+nQpLi52fUvw0C2vlcmidTulvEKkoIFI75aFcs2xDIwI2+bNm2XYsGHm9RdffMHEKkUsE1g0YMAAHmbkZczMGVJesdO81oBg1ZcFMnDgQNe3hYxbu2mrXP3oAlmwaoP07dBCbj27r7Rq1lh8mkzBDgoIAQ/oQF6oKQERcz2iQwvXtwQPaCDw8nvrZEPZdnPV74G9IRgAPKAzugFdWsrXSxqaq36Pr2a/oybPk9LxM8xVv8cumhHYqakkEXNduCr5OhN+/mGgZsBizQBrXkDy9ANIZ736YadZEw2Wpozu5/q2MsHGzybNv4Px0x4yAwC8ZmP26ysbGSV+/mGggBCA9/UUVWem1FN8RYsF086S8PMPA5kBwEOs036Fegq3+PmHgZqBlLHmhTSwTo4YMH7awzIB4CGX67S+710HUB3LBICHXPYdYO86EB6CAcBDLtdpqR4HwsMyAeAhG1Xi+0L1OBAeMgOoF6ra40P1OBAedhOkLPRqWKraAaQl9PEzS8gMoF5YPwYA/xEMoF44TQ8A/EcwgHph/RgA/EfNQMpY80LSaPqDWDB+2kNmAPAMTX8AJI1gAPAMRZsAkkYwAHiGok0ASSMYADxD0SaApFFAmDIKYADUFsWhu2P8tIfMAABkBMWhcIWDihANZl3IOopD4QqZAUSDWReyjuJQuEIwgGgw60LWURwKV1gmQFSzrqonLDLrQtbostX+Tv1kuQtpIDMQGB0o9Fjh0vEzzFW/xy7MuhAClruQBjIDgQ4UOvvNDRT7m2nEojazLiDrWO5CGsgMBIaBAggbRYZIA8FAYBgowsdSUNxY7kIa6EAYWAetXHGRZgQ0EKC4KDwaAFQthNQPBJY/ECI6ENpDzUBgqEYOH0tB1fFMA/XDMkGEqEb2G0tB1fFMA/VDMBChfGeWrFVnA2vG1ZEtAeqHZYII1aX5TtX0q85FN27eLjsrhG2LDrFFsjoaSgH1Q2YgQnWZWVZNv64v2xUI+D77IsMRHrIlQP2wmyBlvlfD6gemBgJ78rmKnWp8wA++j58+ITOAWher6UVnXr7PvlhfBoDdUTOAGukHfmh9C1hfBuqHrZzhYZkgZaS5sofGTIAfS22Mn/aQGUB0qMZHvpgR78JSW3ioGQCAWqK50S40vgoPwQAA1BIz4l3YyhkelgkA7BNp8d1RfLoLS23hITMAYJ9Ii++OGTFCRWYAwD6RFt8dM2KEiswAgH2iUAyIA8EAgH0iLQ7EgaZDKaNpBgDkh/HTHjIDAABEjmAAAIDIEQwAABA5ggEAACJHMAAAQOQIBgAAiBwdCBEM+ugDQH7IDCAY9NEHgPyQGQhI7DNj+ugDQH7IDAQk9pkxffQBID8EAwGJfWZMH30AyA/LBIHNjDUjoIFAjDNjjpcFwrLui23S+qwbpVH7bvI/Dy2U2885KqqlT5vIDASEmTGAkFz/1GJp0vlIKSxuLnOWfxbd0qdNZAYCwswYQEjeWr1RGhQUmtc7KyS6pU+byAwAADKpT/vmUlG+07wubCDRLX3aRDAAAMikm07vIVveny87yz6XEw49kKXPFDWoqKjYVX6OVHz55ZdywAEHmNdffPGFNG3a1PUtIUWx93oAksT4aQ+ZASBBsfd6AOAnggEgQbH3egDgJ4IBIEF0QQTgI4IBIEH0egDgIwoIU0YBDADkh/HTHpoOAagROySA8LFMACDvHRIaKIyaPE9Kx88wV/0egH8IBgDkvUOCrZRAGAgGAOS9Q4KtlOEh2xMnggEAee+QYCtleMj2xIkCQgB5n4apgYF+WGhGQAMBtlL6L6lsD4WnfiEYAJA3js0Oj35wa0ZAA4H6ZHtyGQb9c3IZBp6V7GKZAACQeOMs6kn8QmYAAJB4tiepDAPsIDMAAEgcrbn9QjvilNFOEwDyw/hpD5kBAAAiR80AAFjEljtkEZkBALCIpj7IIoIBALCILXfIIoIBALDIlxbOnFEQF4IBwAEG2nj5suWO5Yy4UEAIOECr1nj50sKZ5Yy4kBkAHGCgTQ9Zl7iWM5AMggHAAQba9JDejms5A8lgmQBwgKN/00PWJa7lDCSDYABwgIE2PRyQA9QdywRIFOu1cI30NlB3HFSUstgO2tAAoOqsTAdjZsAA8hHb+OkSmQGkul770tK1ZAgAIOMIBpBalbzSsICKbmDvWFZDVhAMeC5rg0luvfarcICKbsSlLu9JtkEiKwgGPJe1wSRXJT+wWyv20SNKdXlPsg0SWUEw4LmsDiZUdCNWdXlP0nwKWUGfgYD2VOeW6jU9qf+7fgDrTN0F9tEjVnXpc0DzKWQFWws93xqj65G5wURt3LxddlbsmmWwrQ+wr+p78gjHQbnv2FpoD8FAQA+zZgR0nTJHU/Tzxw1N7e8DgDQRDNhDzUBAWH8EAOSDYCAgFO0BAPLBMkHKSHMBQH4YP+0hMwAAQOQIBgAAiBzBAAAAkaPpEAB41L9AOxy6biqG8JAZAAAPZO0cEoSFYAAAPJDVc0gQBoIBAJk/Ghs0FUO6CAYAj6X1oU1KOntoKoY0UUAIeCz3oa1p49yHdhKHU8WYks56gR4ngSJNZAYAj6X1oR1jSppsCGJGMAB4LK0P7RhT0jFmQ4AclgkAj+mHtM5g9YNLA4GkPrRjTElrYJVbcoklGwLkcFBRyjhoA/CrZqBqYJWlmoEYMX7aQzCQMh5mAMgP46c91AwAABA5ggEAACJHASEAIJN9HkqYrlrDjxoA4Bx9HtwiGAAAOEefB7cIBgAAzsXY9TJLCAYAAM7F2PUyS+gzkDL2yQIIic0DnRg/7SEzAACoNQr9wsTWQsChrB+bC+yJQr8wkRkAHGKWBQ0IR02eJ6XjZ5irfp9lFPqFiWAAcIhZFnwLCCn0CxPLBJ4jzew3js2FbwFhjMdbx4DMgOd8m1Vgd8yyEFPa3bclkZiwtTBlaW+N0TeVBgI5+qEyf9zQRP8OAOln9zQjcETg2T0NAKpmwjQArinLwNZCe1gm8BxpZv+wtINY0+6+LYnEhGUCz6WdZiatlzyWdhCrmJZEfMMyQcp8T3PVNa2H/WNpB7Gq65KI7+OnT1gmQI1I6yUv1KUdlj+wPzEtifiGZQLUiLRe8kLdQcDyB+AvMgOokX5Q7ZnWQ/2EOjsiiwT4i2AAUX5wwd/lD5YjgOSxTADAq+UPliOA5JEZAOBVFonlCD+R0ck2MgNAwujNkC6KWv1ERifbCAaAhDHopSvU3RihI6OTbSwTAAlj0EsXRa1+CrW/RijIDAAJI40NVEdGJ9toR5wy2mnGJ6ZT6IA0MX7aQzCQMh5mAMgP46c91AwAQIDYyoe6oGYgMGxrA6DY1YK6IBgIDAMAAMWuFtQFwUBgGAAAKHa1oC4IBgLDAABAsZUPdcFugsCqYdnWBiAU7Cawh2AgZTzMAJAfxk97WCYAACByBAMAAESOYAAAgMgRDAAAEDmCAQAAIkcwAABA5AgGAACIHMEAAACRIxgAACByRa5vAPnjvHIgLrznkRYyAx7juGIgLrznkRYyAx7juGIgrozAS0vXSu4wGd7zSBKZAY9xXDEQV0ag6qlyvOeRJIIBj3FeeXizv1GT50np+Bnmqt8De2YBlU4BeM8jSRxhnDKO4ERtaQCgsz8d9HXWp4P9lNH9XN8WMsC3ZyOpQkfGT3vIDAAZQQ0IQskCUujoHwoIgYzQGVTV2R/rwcjRWXWWMwF7IrD1D5kBICN8m/0B+0Jxs3+oGUgZa14AYpOrGdCMgAYC1AxkH8FAyniY80OnNQCMn/awTIBMogAJAOyhgBASewESWQgAsSMzAC8KkHq0bZ5aQx6yEABiRzAALyrrVVof2GyDAhA7lgngxb5qzQik9YHN/n4AsSMzAIl93zL7+wHEjq2FKWNrTLb2LQPwB+OnPQQDKeNhBoD8MH7aQ80AMoHtfQDgDjUDyAS29wGAOwQDyAS29wGAOwQDyAROOQMAdwgGkAls7wMAd9hNkDKqYQEgP4yf9pAZAAAgcmwtBABUw3bfuJAZAABUw3bfuBAMAACqYbtvXAgGAADVsN03LgQDAIBq2O4bF7YWpoytMQCQH8ZPe8gMAAAQOYIBAAAiRzAAAEDkaDoEABlDwx/YRmYAADKGhj+wjcwAAC/ENFum4Q9sIzMAwAsxzZZp+APbCAaABGeuoybPk9LxM8xVv0dyYpot0/AHtrFMACQ8c9UPqtzMdcrofq5vK6jZcu7nG/psWZc/eHZgE5kBICExzVx9ny2TxckWfh/u0Y44ZbTTjIcOYlVnrvqBFdrsLpQivhh+VyH8Phg/7SEzACQkhnXeUIr4yOJkC78P96gZABISwzpvKIN2TPUHPuD34R6ZgQCx/oa0hLLlLYYsjk/4fbhHzUDKXKx5JbEeGsracCj/HVn7eWpGQAMBfp5IEzUD9hAMBPgwa0ZA13RzNNqeP25olAVWofx3ADEiGLCHZYIAJZHKDWVtOJT/DgBIE8FAgJJYfwtlbTiU/w4ASBPLBCnzNc0VytpwKP8dQIx8HT99RDCQMh5mAMgP46c9LBMAABA5mg4hOmw3BIDdkRlAdEJpqQsASSEzgOiw3XD/yJ4AcSEzgOhaKbPdcP/IngBxIRhAdB9C9EHfP7InQFxYJkCm2PgQiuF0wfriFDkgLmQGkCmk8LOxxEL2BIgLTYdSRtOMuqFjYP44lAmhYfy0h2UCZAop/Pyxzu83dnDAJZYJgECwxOI3dnDAJYIBIBCs8/tdx1HXzI6NbbiIBzUDKWPNCy7o2/qzzz6T1atXm69169bJxo0b5fPPPzdX/dq+fbuUl5dXfum/06RJEykpKTFfxcXF0qxZM2nVqpW0bt268qtly5ZSWFjo+j8xuDqOutZ8JF0jksVlCsZPe6gZADylH94rVqyQxYsXy9KlS2XJkiXmumzZMhMAbN2azkyxYcOGcsghh8hhhx0mhx56qPnq1auX9O3bV9q1aycNGuxaqkDdZvv64btn8WxSf3Zdlin0z8otU1C/Ew+CAcAT+sE/d+5ceeONN8zXm2++KevXr6/x39FZ/MEHH2yuX/va18xX8+bNzYy/cePGUlBQUPmltmzZIps3b5aysjLzpRmEtWvXyieffGK+NMOgGYV3333XfO3poIMOMkGBfp1wwgnSv39/ad++vcSqLv0a6lo8m3QvCApQ48YyQcpIcyFf//3vf2XmzJmVXxoM7G2W3rVrV+nevXvlV5cuXaRDhw5mlq4f+EnasWOHyTosX7688uu9996Tt956y2Qmdu7cWe3f0SyCBgUnnXSSnHLKKdKxY0eJRZpbZZP+s7O4NZXx0x6CgZTxMKO2dMb98ssvy9SpU+Xpp5+uNvPWdfrS0lI55phj5OijjzZfhx9+uDRq1EiyQLMKixYtkoULF8rrr78us2fPNq+1HqEqXVIYNmyY+dIAISv3H7ss9vhg/LSHYCBlPMyoybZt2+TZZ5+Vv/3tbzJt2jTZsOGr1GxRUZEce+yxMmjQIPN14oknVj5LvtBlhldffdUEOc8//7xZ5qgaHLRo0UK++93vytlnny1DhgwhMMBuGD/tIRhw8DBnsWoX9uiH4axZs+Shhx6SRx99dLd1f13bP/XUU+W0006ToUOHmvX9kOgOhxdeeEGmT59ugp81a9ZU/jOtZzjzzDPlwgsvlAEDBlCICIIBiwgGHDzMWVybQ/q0AO9Pf/qT3HfffWatPadt27ZyzjnnyFlnnSXHH398NNv2tL5AlxIeeeQR+fvf/y4fffRR5T/r1q2bjB49WkaNGmV+PogTwYA9BAMOHmZtEqJdxnK0Scz8cUMd3iXSom8vzQLcc8895gNP6wKUzvj1w/+8884zSwCxBAA1ZUt0KWHKlCny8MMPm/eN0p+L/pyuvPJKOe6446zfF1k8twgG7CEYSBmZgbhUfnis3CBtGm6RT6fdIW+88q/Kf96vXz+55JJLZOTIkaaxD6rbtGmTyRbcf//9MmfOnMr/XbcqalAwfPhwU09hA+9VtwgG7KEdsQO0jQ3X2Ifny0tLP5ENm7fL4g0iK9udZLr6XXTRRaY3gBbQ6Zo4gcC+aQ+EH/3oR2YJ4d///rdZKtDCQg0MRowYYbZP6nJLLsuSJvbeIxZkBlJGZBsH3X//wAMPyP8uLBFp/FXFf2PZLrPGnmja+CJ/Wmh49913m+UWbXykOnfuLDfccINccMEFpt9CGsgMuMX4aQ+ZAaCea926xq37/bXgbfOqxVJRvqvxTmEDkeO7tScQSIAWEY4fP940Yrr11lvNz/T99983GYQePXqYXRlpzGvI4iEWZAZSRmQbLt0/P2bMGHnttdcqW/Fefu0NsuTrx8mi1Zsy07gl1PfV73//e5kwYYLZpaG0D8Ntt91mdmQgDIyf9hAMpIyHOTzajve6664zywK5Ne6rr75arrjiiuD6AmSdvqc0U3DLLbeYsxSUbtPU/03PZIDfGD/tYZkAqENdgM48dQ98LhDQYkA9KXDcuHEEAg7oB8WNN95ofgf6u9BGRdrNsWfPnnLXXXft9awEANURDAC1sGDBApN+1gyAzlb09bx582Ty5Mk0xckAzQLo70J3bGg/At2eqEs4+nvS0x0B1IxgAN7v69eKb23kpFf9PunDd7RiXQ8H0g8abZn7xz/+UV555RVzbgCy5cgjjzRbEnXXgf6u9MAk/T3p71DPgQCwd9QMpIw1r3SlufXrP//5j5x77rnmqrRv/qRJk8zRwPBjO6LWcehuj1ygoMs7vXv3dn1rqCXGT3vIDMBraTSF0fhYP/Q1G6CBgG5je+yxx0w7YQIBf+jyjdYP6LbDAw88UObPn2+Ofb799turHasMxI5gIND0diy0X7xmBJRedTtffaxdu9acGHjZZZfJ1q1bzQmCCxculO9973sJ3TFs07MNNKj79re/bX6nY8eONccm6wmKAHYhGHBE+9drelsPLNKrfg+3TWG0IPCoo46SZ555Rho3biwTJ06UqVOnSps2bRK9Z9inGR39vWotgf5u9feqv+tcjwggdtQMOFrz4uTCbNGiwJ/+9KemyEy3DuqyQJ8+fVzfFlKgywVnn322LFu2zLQxvuOOO+QnP/mJ2ZaIbKFmwB4yA4Gkt5EfTRtffPHF5iAhDQQ0fayzRQKBcGkhoe4M0dMP9bAjDQIvvfRSKwcfAVlFMOAIPc/dW79+vQwbNkz+8Ic/mFnhTTfdJI8//jjNgyKg2w61IPTmm282v3ttbaz1IfpMADFimSBlpLmySQ+80cH/nXfeMe2EdfuZFpghPk899ZScd9555r2qxyNrPUGXLl1c3xYYP60iM4Do6DKAdqnTQEA7182aNYtAIGKnn366aSL1jW98Q5YsWSL9+/c3dQVATAgGEJWZM2fK4MGDzUl3ffv2lblz55or4qbPgO4mKS0tNc/GoEGD5KWXXnJ9W4A1BAOIxj/+8Q+TAdDU48knn2wyApxsh6pNijRYHDhwoGzcuFFOOeUUefrpp53cC31IYBvBAKKghYG6U0DPGtCmQjrIa60AsGdh4fTp080zos+K7jjQZ8c2+pDANoIBBE97BowYMcJsHRs5cqSpIm/SpInr20JGFRcXmwDg/PPPN0cg6zPz5JNPet9mG6gJwQCCpl3n9LAhHdRHjRolDz74oGk0A9SkqKhIpkyZYnYZ7NixwzQp0l0GttCHBLYRDEQktnXIf/7zn+ZMAR3MNSC4//77pbCw0PVtwRP6rGhAoJkBzSrpszRjxgwrfzd9SGAbfQYi2ieb5nG/WTNnzhwZMmSI+flrrYCeXEdGAPnIBZO63KTv3xdffNGcfoi4xs/QkRmISCzrkO+++64pANOBZOjQoaahEIEA6rNkoMtLugNFnyltVqXnGgAhIRiISAzrkOvWrTOD9aeffirHHHOMKQTTU+qA+mjUqJEpPM31IdBth3oFQkEwEJHQ1yF1K9gZZ5wh7733nnTq1MlsHyStiKTomRXaq6Jz584mM6DPmh50BYSAmoGUseZlhz7GF1xwgUnn6l7x2bNnS69evVzfFgK0dOlS0856w4YNMnr0aHP8Nccfp4Px0x4yAwjCxIkTTSCg67u6NEAggLR069bN1KEUFBTI5MmT5a677nJ9S0C9EQzAe9pD/qqrrjKvb731VvnWt77l+pYQOC1MveWWW8zrsWPHmm2sgM8IBuC11atXm+6C2lRIG8SMGTPG9S0hEldeeaVZmtJnT7cerlmzxvUtAXkjGIC3coPwxx9/LH369JH77ruPtVtYo8/avffea5493Vnw/e9/3zyTgI8IBuCt3/zmN2aJQAuMtE6A4iK4OMdA6wdKSkrkhRdeMM8k4COCAXjZ/ljPnv/Vr35lXv/ud7+TLl26JPZnA3XRs2dP8wyqcePGySuvvOL6loA6IxiAd8ewbtq0ydQHaEr2nHPOMeu2gEt6CJYuE5SXl8sPf/hDKSsrc31LQJ0QDMC79sfXXnutafrSsWNHueeee6gTgHP6DE6aNEk6dOhgml794he/cH1LQJ0QDMCr9sezZs0yAYD685//LC1ahNdSGX7SZlfagEjdeeedpp4F8AXBALxpf6zthn/84x+b1xdddJEMHjw4oTsFkqFnFuSe0QsvvJDlAniDdsQpo51mcq6//nr59a9/Le3atZO3336brAAyaePGjdK7d29ZtWqV/PKXv5Tx48e7viVvMX7aQ2YAXnjnnXdkwoQJ5rVWbhMIIMsHGukygbr55pvNkdpZ2okD7A3BADJPk1fa7W3Hjh1y2mmnyfDhw13fElAjfUZ1yWDbtm1y2WWXmWc4CztxgH0hGPBMjDOGZ555Rp599llp2LCh/Pa3v3V9O0CtdhfoAUaNGjUyz+4TTzyRiZ04wL4QDHgmthmDzqw0K6D0SnMh+KJr165yzTXXVG6H3b59u9OdOEBNCAY8E9uM4e677zb7ttu2bSs33HCD69sB6kSDgNatW5tnOLft0MVOHGB/2E3gWTWsLg1oRkADAZ0x6EAxZXQ/CZH+vA499FBZu3atOYRItxMCvtGC15/97GfSpk0bExTkxoNQ6FKlZih1oqIZDQ1cWjVrnMifzW4Ce8gMeCamGcPEiRNNIHDYYYeZFq+AjzSI1WdYT9e87bbbJDSxLV2GisxAyohs87Nhwwbp3Lmzuf7lL3+R888/3/UtAXl75JFHZOTIkWbb4YoVK4LaGqvFzBoI5OhEZf64oYn82Yyf9pAZQCbdfvvtJhDo1auXOYwI8NlZZ51lGhFpQ6LcCYehoNgxDAQDyBydAei2LHXjjTdKYWGh61sC6qWgoEB+/vOfVwa6OuMNRUxLlyFjmSBlpLnyqxW4/PLLzTbCxYsXEwwgCNo0q3v37rJ8+XLTLyO3ZdZWMZ6PGD/tITOATNG92Lkiq6uuuopAAMEoKiqS6667zrzWZ3xvfQcoxoMrBAPIlEcffVQ++OADadWqlYwaNcr17QCJ+sEPfmD6Dnz44Yfy5JNPSux9RJAdBAPIlEmTJpmr7ssuLi52fTtAoho3biwXX3yxeZ2ri6mKYjy4Qs1Ayljzqr2FCxdK3759TTp15cqVpusgEBo92viQQw6RnTt3yoIFC+SII46oVjOgGQENBKgZYPy0hcwAMuPee+811zPOOINAAMHq0KGDnHnmmeb1ntsM9YNfO4rqPn29xhwIwC6CAWSCRv0PPPCAeX3JJZe4vh0gVZdeeqm5/vWvf5WysjLXtwMQDCAbHnvsMdm0aZPZTjh48GDXtwOkauDAgWapQJ/5vRUSArYRDCATHnroIXPVHQTaoAUImT7jurNATZkyxfXtABQQpo0CmP1bs2aNHHzwwVJeXm5OddNDXYDQLVu2zGTCNDDQgtn27du7vqXMYfy0hykYMnGIiwYCxx13HIEAoqHPev/+/c2z//DDD7u+HUSOYADOaRGVOvfcc13fCmDViBEjzPWJJ55wfSuIHMsEKSPNVbPVq1ebJYLc63bt2rm+JcAa7bbZqVMnadCggVku0+6EtTmnIJYzDBg/7SEzAKemTZtmrv369SMQQHQ6duwoRx99tOic7Kmnnqr1OQWcYYCkEQzAqalTp5rrd77zHde3AjgxfPhwc3388cdrfU4BZxggaQQDcGbLli3y3HPPmdcEA4jV6aefbq4zZ86UrVu31uqcAs4wQNIIBuDMiy++aLqvac1AaWmp69sBnDj88MOlTZs2snnzZpkzZ07l/651AAO6tJSvlzQ0V/2+Nv8MyEdRXv8WkACdCakhQ4aYAiogRvrsn3zyyfLggw/K888/L4MGDdrtnIK9qemfAfkgMwDnwcBJJ53k+lYApzQYUBoMAC6QGQiEb1uNdMvQ66+/bl7nZkJA7MHAa6+9Jhs3bpTmzZtb+7t9GzuQDjIDgfBtq9Hs2bNlx44dZmuVHtgCxH6ssfYb0G6EuSDZFt/GDqSDYCAQvm010mBAffOb33R9K0AmaDtuNXfuXKt/r29jB9JBMBAI37YavfHGG+Z67LHHur4VIBOOP/54c3311Vet/r2+jR1IB8FAIHzbavTmm2+aq3ZfA7B7ZsBml3jfxg6kg7MJUkZv7eo+/vhjadu2rdlSpcVSuZ8PEDPtM6DvBa0b+PDDDznSmPHTKjIDcJYV6NGjB4EA8P+Ki4ulS5cu5vWiRYtc3w4iQzAA695++21z7dOnj+tbATKld+/eu71HAFsIBmDd0qVLzbV79+6ubwXIZDBAZgC2EQzAuiVLlpgrwQCwu169epkrwQBsIxiAs8xAt27dXN8KkCm598T777/v+lYQGYIBWKUVwR999JF53bVrV9e3A2SKduRU+h6pepwxkDaCAVilW6aU9l5v0YLmJkBVLVu2NLsK1KpVq1zfDiJCMACrVq9eba7soQaq094buezAihUrXN8OIkIwAKsIBoCa5YKBDz74wPWtICIEA3ASDLRr1871rQCZ1Lp1a3P99NNPXd8KIkIwAOutiJW2IwZQ3YEHHmiuBAOwiWAAVulZBIriQWDvDjroIHP97LPPXN8KIkIwACfBgO4mALDvzADBAGwiGICTYKBZs2aubwXIpFzWbP369a5vBREhGIBVmzZtMlcyA8DeNWnSxFy3bdvm+lYQEYIBWFVWVmauJSUlrm8FyKRGjRqZK8EAbCIYgFU7duww14YNG7q+FSCTCAbgQpGTvzXCmbDq2bOnNG3aVGK2fPlycy0q4tEDagoGFixYsNd+HNqlsC7f+/zvlJeXV/v/Ix2MyCmrqKiofL1y5Uqn95IVBQUF0rlzZ9e3AWSSHu2tAYFmBtasWeP6djJBO5bmaimQjgYVVT+tkDh9Q0+YMMG87t+/vxQWFkrsOnXqZL4A7N26desqu3Xm7DlU1/d7X/5M/b53796V/ReQDoIBAAAiRwEhAACRIxgAACByBAMAAESOYAAAgMgRDAAAEDmCAQAAIkcwAABA5AgGAACIHMEAAACRIxgAACByBAMAAESOYAAAgMgRDAAAEDmCAQAAIkcwAABA5AgGAACIHMEAAACRIxgAACByBAMAAESOYAAAgMgRDAAAEDmCAQAAIkcwAABA5AgGAACQuP0fea62JrZ300YAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "do_one(NKD=100, Nsample=100, plot=True, close=False);" ] }, { "attachments": {}, "cell_type": "markdown", "id": "799a54e4-dc2d-4551-a040-9ea475124465", "metadata": {}, "source": [ "As you increase the number of points $N$ inside the domain, the distance to the nearest point goes down like $N^{-1/2}$ and in general the scaling should be like $N^{-1/D}$ where $D$ is the number of spatial dimensions (I think). \n", "\n", "The required memory is linear with the number of points in the K-D tree(!)" ] }, { "cell_type": "code", "execution_count": 4, "id": "eafaec7f-87e7-4a6c-ac95-f908a4c42527", "metadata": { "execution": { "iopub.execute_input": "2025-01-06T11:32:11.288160Z", "iopub.status.busy": "2025-01-06T11:32:11.288045Z", "iopub.status.idle": "2025-01-06T11:32:15.361678Z", "shell.execute_reply": "2025-01-06T11:32:15.361430Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkUAAAG5CAYAAACAxkA+AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAWQJJREFUeJzt3QlYlGXXB/A/IJsoKLjgjrmGG26gibmnZrZoZqtbWZ+pLb5W+lqZb4uvrVb6alm2aWpaauWaW+aWpuIuliupqIiCooLA813nnmaYAQYeYJj1/7uuueBZ5pkHmpjjuc99bi9N0zQQEREReThvR98AERERkTNgUERERETEoIiIiIjIgEEREREREYMiIiIiIgMGRUREREQMioiIiIgMyvzzlXTIzs7GmTNnUL58eXh5eTn6doiIiEgHacl45coVVK9eHd7e1vNBDIqKQAKiWrVqOfo2iIiIqBgSEhJQs2ZNq8cZFBWBZIiMv9Tg4GBH3w4RERHpkJqaqpIaxs9xaxgUFYFxyEwCIgZFRERErqWw0hcGRU4sK1vD9uPJOH/lBqqUD0B03VD4eLOWiYiIqDQwKNJh+vTp6pGVlWW311y5/ywm/XQQZ1NumPZVCwnAxL6R6NW0mt3ug4iIyFN4aVKSTbrHJENCQpCSkmKz4bP8skG/HEzEiDm7kPs/jDFHNOPRVgyMiIiIbPz5zUyRA+WXDQoP9seNzOw8AZHQ/gmM5Dk9IsM5lEZERGRDbN7owIBIskHmAZFITE3H5Ws3rT5PAiN5jmSXiIiIyHYYFDmADJlJtqck45Yy3EZERES2w6DIASTLkztDhCKWdkn9EREREdkOa4ocIL8sz+N/LMVdh37DssYdsKJRLE6HVMn3uVJFFB5iKMgmIiIi22GmyAHyy/LceXgTWp6Nx8vrZ2PzzGEYvXlennOMZdUyLZ9F1kRERLbFoMgBJMsjPYfMw5oR947Hq92fwu+1miIbXjgWcSvCg3OCp1su/o2WWZc4HZ+IiKiUsE+Rg/oUGWefCfP/ABIoVb6ajNcf74TuLWqZehh1mDASlZYvBdq0AQYMAB58EKhdu4Q/ERERkftLZZ8i5ybZHsn65OlTJF2rH+2Onv9kg9rXCzMUYftkA97ewB9/GB4hIcBTT9n0nrisCBEReTJmipywo7XVQOT8eeCHH4DvvwfmzgWqmBVjT58uN2jIItWvX+T74LIiRETk6Z/fDIocHBTZhKZBq1sXXidPqs20W5shYPpH8OnSWVcAZhzK47IiRETkjjh85kFW7fkbu1r1QwevDbjt5B4EHdqHQT/E4+HKjUzBzC/bj+LVtSfyZIJe6XMrXl92iMuKEBGRx2OmyMUzRbmzPBWvpeD247vwY2RnwMtLZXlEWO/uKJd+Dcsax2J5o1gcC6upAh69//HnDW9nqG8qAtYoERGRM2CmyEOXC7lUNgRLm3RR30v48dqPB1D++lX8fPYI/LMyceuFExj72xw8MvANbI6IKrVlRVijREREroZ9itxtuRAz2j8LzP550xcxI7/GC72fwYa6rZEcGIw/akZanNszfgvqJSXYZFkRq4vdptxQ++U4ERGRs2GmyIUVJXtzOTAYC5vfoR7+mRlIL+NnOibb7y7/AOUzriO+Um01vPZNqz5ILhtS5GVFClrsNneNkuDwGhEROQsGRS6suIvCmgdEouK1VOyo2QSxJ+LQKOkUGiTNw7dRvQpcVsRavZCe7JUcn7buL8zfcYrDa0RE5DQYFLnBciEyLJVfZkbCmKrB/uq7c6nWz/GqVRM3Fv+IPou2o9muX9EgKQEXyoXmBClNwoG77gJat1Z9kFZqoZj086F8A5r0zGxd9/7BmiN59hmH19gCgIiIHIGzz9xk9ll+y4UI4+yzws6RIMTqbLHdu4FWhuuIv0Jr4sGHJyMpqGKeaz3XvQE+WPNnsX8e43Ddppe6ciiNiIjs+vnNQms3WS5EAglzsm0MdvScIyQIkWn390TVUF9NQYl0yP7mG2h39UWGjy/8s24iqWwFi2vVvJyomkjO235KLWRb3HDGOLwmwRkREZE9cfjMDUhQI4XLBRUt6znHqvLlgUcfxbb2vfFk/bWoLQGQV87zwtIuY8OnT+JExWqqSDt08MN4OcUbXl5eeTJTWim1ACAiIiopBkVuwpjlKek5hQUqV/yDcKBqPYv9Tc8dRaa3D+oln8borQuQGfcTKm/ch4m/HM+z2O2DbWvpGl4rbhE5ERFRcXlcUHTfffdhw4YN6NatGxYtWuTo23Ep1gKVX29pjVaj56Lb0R24M34TWjWLwB1tbkG3VnVVZupCcipafTUd1e5/CGhZH/N3JBRYHF6UFgBERES24nE1Rc8++yy+/vprR9+GS892y2/ALc2/LH6K7IRJQ15H2II5Fpmpuy8cQs3p78GnbRv4NGqIuUeXwEvLznOdgloAFEaKxLcevYilcafVV9kmIiIqCo8Lijp37ozyUiNDRSaBigQsosCAxifX2yosDOjfHwgMBI4exS0Hd+J/j7XJKfyWCZCalqfwuygz8GKnrMNDs7bh2flx6qtss3M2ERG5bFC0ceNG9O3bF9WrV1dFukuWLMlzzvTp0xEREYGAgADExMRg+/btDrlXT6V3JpuFmBhAhirPnwcWLAAmTFDnybR7WWh2dttAHP52JDZfXIFe1xIMQZLODBCXFCEiIresKUpLS0OLFi0wbNgw9OvXL8/xBQsWYMyYMZg5c6YKiKZOnYqePXsiPj4eVapUUedERUUhMzMzz3NXr16tgq2iSE9PVw/zPgdUgpls5coBDzyQt/D7m9+Av08B771neMTGAr/9VuiiskVZUoQ9j4iIyKWCot69e6uHNe+//z6GDx+OoUOHqm0JjpYtW4bZs2dj3Lhxal9cXJzN7mfy5MmYNGmSza7nTko6k83Ciy8CTZsaskk//wy0aWPKAElwE3zjKupd/Bu7qzey6HodEuina0kRCd5sdq9EROS2nGr4rCAZGRnYuXMnunfvbtrn7e2ttrdu3Voqrzl+/HjV/dL4SEiwvoo8lUDZssD99wPz56shtqx/T7DIAPWK34LFc8Zi84xhmLB2FqLOxKvjian6ehmx5xEREblVUJSUlISsrCxUrVrVYr9sJyYm6r6OBFEDBgzA8uXLUbNmzQIDKn9/f9UO3PxBpaxsWWxP9bLIAFW4cQVX/QJR48oFPPHHUgz//Xt1PPlqztBmQdjziCgfkpVt1Aho0AD47LP8z4mIAJo3l7oEoEuXnP333QdUrGj4x4zR5csqy6vOlczvrFml/zMQufPwmT2sWbPG0bdARczsfBrTH1+1ugudju9SfZCka7YIDfJTNUY1D+zEnYc3YVnjWOyscSs0L0Osz55HRFZI3eWYMcD69UBIiGGxZwl0ZKZoblu2GOoBzT37LDBsGPDVVzn7ZFbvxo2GzG9amiEwktrQ/K5J5KRcJlNUqVIl+Pj44Ny5cxb7ZTs8PNxh90W2l19mJ93XH6sbtsdzfV9QX0V4SKAquu63fx2G7vwJi+a+hK3/G4JhO5aWqOcRuYiLFwGZYHHiBNzOgw8aJh2UFpm126QJUKOGIeCRWs7Vq/U/v3NnQxBkzsfHEBAJmaDyT6sNIlfiMkGRn58fWrdujbVr15r2ZWdnq+327Q0fkqVF2gBERkaibdu2pfo6VHiTSCH75bicJ7PQ6j89BCtadkeqX1mEX01GmexMyxYB168DWVl2/imo1L35JnDPPYYhnuKaPt3w/IAAQ+uIwlp8vPaaYd0/80fjxjnHr1wBnnsOqFPH0JfrttuAHTssryGvl/sa8hg5Muecl182/HwpKZbP7dTJcO68eZb7P/4YKMrs2jNnDAGRkXx/+nTe8+S15DXlb9/cuYVfV4bQWrQAatYEXnhB/jWr/56InIBTDZ9dvXoVf/31l2n7+PHjajZZaGgoateurabjDx48GG3atEF0dLSaki/T+I2z0UrLyJEj1UOm5IdIqpns0iRSZpnlXkQ2vwxQ2ycfRNYTA7Hj8FlgzS/o2bo1XmrfJCdDJP/inj4d2f364VCHnvirURSqVAjSvyAuOZ9r14DPPwdWrSr+NaRnlgwhzZxpCIimTgV69gTi4w0ZKGskw2I+DF/G7M/oE08A+/cD33xjCFLmzJFCRuDgwZwgRIIk8yBdzu/RAxgwIGefDD3Vq2d4vjFYkqzL7t1AtWrA998DDz2Uc/7OnUCrVjnbUteTT2sSlQ0qSvC0aZPhvs+eNfwczZoZaoysqVAB2LNHUviGoTOpOcpVB0rk1DQnsn79evn8y/MYPHiw6ZyPP/5Yq127tubn56dFR0dr27Zts9v9paSkqPuRr1T6Vuw7o7V7a41W56WfTQ/Zlv1F0rGjMZGvHjOi+xX/Wq6sRg1Nmz7dct/mzZoWGKhpJ07Y5jWysjTtrbc0LSJC0wICNK15c01buDDn+Pnzmla1qqa9+ablPfj6atqaNYbtTp00beRIwyM4WNPCwjTt5Zc1LTs75zlyzcqVS3av0dGG1zC/9+rVNW3yZOvPmThR01q0yP/YtWua5uOjaT//bLm/VStNmzDB+jWffVbT6tWz/PnEpEmaFhubsx0fb3gPy3/DsmU1LS0t51izZpr26quabvI7v/dey3uYO7fg54wdq2lffJGzvX69pvXvb/38ESMs/9sTOZDez2+nCoqcHYMi+8vMyta2/JWkLdn9t/oq20W1ctcJbcj9E7XvmnbXLvsHaQMe/q8KiiL+eWyZ85OmrVunaZmZNntNp9Svn6YNGZKzLR/Cbdtq2r//nfdcCVqCggp+nDyZ93lvvKFpjRtr2sqVmnb0qOFD1N9f0zZsyDln2TJDELRjh6alpmraLbdo2vPP5xyXoKhcOcMH9eHDmjZnjiEI+PTTnHOeeUbTevUq/u8iPd0QwCxebLl/0CBNu/vugoMiuZdq1TStbl1Ne/jhnN+D/CwStBiDO6MOHQw/k7X7kKDPPEg0WrFC0/z8NO3GDcP2t98aAk15Tni4pi1aZNh//bqmlSmjaUuW6P/5b97UtPr1Ne3vvzXtyhVNa9hQ05KSLM+5etXwMwk5R4K77dutB0WJiTnnX76saU2aaNrevfrvicgJPr+daviMyNZNIqXr9Wsr/8LZem2xvl5b/DtrJDK9fSy6Xme//gYQ/7saMjnVuRdernIbNgbVzLeLtktr185ytpAM8UjvrfHj8577f/9n0X08X7mHYaS49q23DENLxjq/W24xDMF88omhNkXceScwfDjwyCOGKdxBQdIp1fJatWoBH3xgqGmRaeP79hm25Xni5MmiDQPllpRkGMLKPbQj24cPW3+eDLN9+aXhnmRISZq7duxoGAKTwmP5uV9/Hbj1VsO1pPZH2n7Ur5//9WQpI6nDGTIk7zH5+TIyAGk5IjVKu3YZhq78/AwzxaTZqawpKMNVMlRmPnxWGBnyk2FlmWafnW1ooGqcJSZDb9IEV4bA5HWE/K7kd2+sq5ShNHldmWUm9UMLFxoKrZ98MicnO3q0YbiNyIUwKNJZaC0P6ZNErkW6WZv3PLrp42txXLKlR8tWQnSFivA9fx61v/saYX0qAU1zgiLzLtrGwEiCrSIvc+IMQZF0fr961RBs/PvfwBtv5J1uLUJDDY+ikHpAqfWR+hhz8sHesqXlvnffNdTNyIep1MP4++e9V7lHIwk25ENc/h+UD18pnpfiaHPys02ZUvA9HjpkWRhdVOYd9yVAkSBJApbvvgMef9wQaMpUdanDkfuUQEVqf+RnzI/URck18wvwpFBbyO9USFBkDHykXkceEojK/sqVDYFkUdx9t+GRm3FVAAloJfDJj7XWJjZcUYDIERgU6cBCa9dVaDdrLy+8escIBN/3MdZ/8h1idq3HmgYxFqeM2fgNKl27jOUXuqLHjOfxy5GkAtdkc1rSi8bb2/AhKh9q8kFqbZKCZHzkURApHq5dO2dbgi2xbJnlzCaRO+g5etQwA0qyFDKlvqgZBZnVdOmS5b5//Sv/jIs5+aA3Pl+CllwtPtR2UVp8SGFxw4aGgFBIcfSvvxoyKLJWohRFDxyY87rmJNsl/x1++CH/aycnG77Kfych/92MxdUyJd7X11BonrvImoiKjUERuTW93awvZmhYWqUplvZqarHfOzsLD+5dhcppl/HQnlW4vmgyFtzxrBqOM5dfNsnpMkrSQ0aCD5m5JN2Gly83BEn5Kc7wWWSkIfg5dSpnqCw/kjl69FFDsCDDUDJjS4bHzGd8/f675XO2bTN0XpZARkjmSWZmmZPgwRhAFEaGoCRIlBYf995r2CcBmmyPGgXdJBCUAO+xxyz3y5CgPCRwk8Dl7bfzPveLLww/c58++V9bhuRkaEoCuGPHDMNsxuBHhr8kyyP/LeV3V8CakUSkH4Mi8oieRxK05NdGztj1OrRcrkzGPzQvLzx311j0ObwJPY9sQVjqZRyvaBkM1Eg5j3PlQpHlU0ZlkHpEhqvARxa1dbqMkgxLSU8b6e8j2QZrijN8JjU1Y8cCzz9vCDBiYw19djZvBmSJnMGDDedNmGDY/9FHhqE7Cc5kyEmWnTCSwEqmyz/1lCFDIvds3sxQps5LLZQEHbLcRHHI9eWepK4pOtowJV8yPObZs2nTgMWLDcGSkJ+vb1/DkJlkuiZONARqxgyOBEBSTyPBnmSPpFePDNflzsjJ70eCInl98yn95n77DbjjDsP3kg2SQE6GHI2knkiCMRlek98pEZUYgyJya3p7HoUE+uX7fFkyZHNElHq8cscItDhzBCdCLYeGZi5+EzVTzmN1g3ZY3jgWO440weVMTb1m7kDMWkbJbqSxngy7vPNO6VxfiowlWyOF05LdkOElyW5I/ZLYsMEQfMjyEsa1BKUOR+5rxgxgxAjDvkGDDHVDEqxI0CHLSkgRr5FkvOS6UssjgVNxSKbqwgXg1VcNxcxSYLxypWXxtRRkSybI6O+/DQGQdNOWn1MCP8liGTNUEuxJsCbnSVApgYs0YZTfuTkZNpPAT4LB/Ny4YSjClvsREhhKQCSBkZHUbkmNlWTeOHxGZBNeMgXNNpfyjELrI0eOICUlhYvDupjCsjYyzBU7ZV2BGaWQQF9cvn7TYn/wjatYO+v/UPnaZdO+TTPm4YVLlS1eK7/s1KaXutp/KE1mG8kHaGkuIVFSksGSAEWCp4JI7ZJkYmSYydowoKuSAFEyVEVZeoOIrDLWBBf2+c1MkQ4stHZ9EvjIsJa1+h49GaWhHSLwwZo/La6bGlAOMSO/QkzCAbVYbduEAzjSqBXOrso5b9DOn3CqQrjKNsnsN7m2BExyLyVpN6CbDNVIRkRmOv35J7B0KdyC1OLIzyPLUxR15pWzk8ySDBkSkV0xKCKPUVjPIwmcZFgrd0ZJsjoSMElQNX9HQp5sUra3D7bWaY5tdZojPNgfL1YIMh0LSr+GCetnwz/rJlL8g/BLg3Z4r+OjOBtcWQVndinElpXLu3Y11LZIYa47ZTllnTF3JMXnRGR3HD4rhfQbubaCAhUZhpNsEqxkkySokvqkh2ZtU9thaZcxest89D6yBVWvJuOmtw/ajJqDlMDyeL57Q8zfcQpnL1839eSxVojtVLPYiIjc9PObQVERMCii4tYnydT+1qcPodGFk5jbqg9Cyvoi5dpNQMvGqs9HYV94PVWkvSmiFTLK+FoUYjvlLDYiIhfCoKgUMCgivZmbgjJKsl2hrC8uX7uJVn8fwg9zXzAdT/Uri+7DZ8KnRnVViP3LwcR8Z7GZZ6YYGBER2ebz282mbBDZtz7pnqga6mvuoSxjfZLUI5mT7ee7N1ABkdhdoxH6P/I2Zre+G2fLhSGxfCWcLxeqskLbjl5UGaJ2J/fA/2a6xXWMQZIclwCNiIhKjoXWOnDtM7LljLef956x6IO0s2akerze7QlUvpqzdMXWY0nwPnUK8+ZPwFW/QKyr1xbLGsVibf1oZPqUsf8sNiIiN8egSAdOySdbzniztvSIBEjny5uf64WaKedwunxl1LhyAXcf2ojbj+9ShdpFWt+NiIh0YVBE5KRLj0gwNa12M8SO+BxRZ46oPkjpZfxUlsgo4OYN3DbxOeCxgcCddxrWNyMiomJhUETkpEuPtLslzBQ87a7RWD2Q69x+Z/ei8o+LAHlIQCRrmslCqVY6PHNqPxGRdQyKiBygsEaRxhllhQVPvR/qDtS8DixcCJw4YehcbR4QyeRSWTA0KIhT+4mICsEp+UXAKflka3oyN7qCGU1D1o4/cDDhEo7Va5JzrT92qLXEEmO74o3AJqpY+5pfoOk6nNpPRJ4glX2KbI9BETmKnr5I+QVOX51ZjYbT3jbtO1ilLu4c+rHzLFBLRGQHXBDWhjgln5x53TZjo8jc/7qRWqSeZTvizRntcemrb1Wh9tp60RbnSAuA9qf2qP2c2k9Eno6ZoiJgpoicjXFJEfMMkTnJ+4QE+uLydVlSRINfVqZaRsRoyB8/4rW1nyLdxxfJHbui2ugngX79WJBNRG6FmSIiDyCBi7WASMi/eFRAJLy8LAIikeYXgKOhNVAv+TSqbVgF1KqClQ3bsyCbiDwSgyIiF6a3cWOFQF+kXL+ZZ4htYfM7sKhZD3S4fhZfB5/EH/VaWgzFdTm6Aw/HrcCKRrF48fxF4PGODIyIyG0xKCJyYda6Y+c2tENdTF1zJP+p/V5eePTJu6BFhuPZKeugISfQuufgBvT4a7t6pPuUwbdxDyBr7Rw1lKZniI3DcETkShgUEXlAd+xRXeujUXi5AvsibT16Mc9Q3PR2D+BkheqqSLvBxQQc8w1WQU7K9Qx1rfSz55Dp7YPUgHJ5htjYF4mIXA0LrYuAhdbkjIyzz2ClwaN5D6KCMjdL407j2flxVl+nwYWTOF8uFP27N8MXm0+o13rx1y/xxPYl2BQRheWNY/FLg3aY8sTt6vz8ZsSxLxIROQILrYk8hN7u2IVN7S9sKO7PynXU1yVxZ0zBToOkU/DLzkTXY3+ox4yY+/FaFbm+V76ZK+2fwEjutUdkOIfSiMipMCgicgMS+EiQUZL6HT1DcRWDfJGclmHaN7z/q6iXlIA+8Ztw5+FNWNY4Fomp6abjnY/uQJWrl7C6YTtcDjT860yuLcEb+yIRkbPh8FkRcPiMPH0obliHCHy++YTu630779+47dRe3PT2wdbazTGj3QBsrdNcHfvwwSjcE1VD97VYtE1ExcXhMxtiR2vyFIUNxYUE+ukPijQNv9VtiQo3riDy/HHcfmI35kb1thyuk3+TeRUe2LBom4jsgZmiImCmiDyFtayMsYN2QUNsVYP91XfnUnPOqZt8Gr3jN+PzNvcgw9c/Z621J4cDJ08iu//92NmqM874lcuTBbK2jAmLtolILy4IWwoYFBHpm+0mCp0R16gSUKUKcPmy2pfp5Y3Rd7+IFY1jTVkgqZMqbBkTLmZLRLb6/PZGCVStWhUdO3bE008/jf/973/47bffcOnSpZJckoicnHGITYIRc7JtzNroOQe+vtg4dxmmdBqMfVXroYyWjd3VG6vzJBMlQdXCGd8j40yi1XsxL9omIiqpEmWKpMYmPj4e+/fvV499+/bh4MGDuH79Opo0aYIVK1bAnTBTRJSjpB2tcy9mW/VKEs6Vr2R6rpemYd3nI1A7+Qx+r9UUyxt1wM+3djTNYjNX1KJtIvIsqfYotPbx8UFkZKR6PPDAA9i6dasKhBYvXoyLFy+W5NJE5OQK6nmk55zci9maB0SifHoaUn0D4aNlqxls8jhY9RbsqhFstccSZ6gRUUmUKChKSkrCqlWrsGzZMuzevRutW7dGr169sG7dOlSuXLlEN0ZEnr2YrSwdcs/gDxB54yI67tmANn8fwO7qjSzOGb9+Nq5WrY7ogJacoUZEjg2KpKaoRYsWGDt2LL755huVOSIisuVitj3vao+pAWH4NKa/RdF26LUUPL5jiapF0n6ejtAakfC781mgYnXTOcbaJM5QIyI9SlRo/c4776Bly5b48MMPUb16dbRp0wZDhgzBu+++i5UrV5bk0kTk5owdtK0Nbsn+av8sZptf0XZYhSD89a9XoLVrp+qPmiX+hfNBoRbnlE2/pr5KBkmG1oiICmLTKfnHjx83FV0fOHAAc+bMgTthoTWR8y1mu/XoRfzrvR8Ree441jSIMV2jTFYmfp8+CMdDa2B5o1j0mTQarWOb2fXnIyLnwD5FpYBBEZHtlbQWaGncaTw7Py7P/panD2PxnLGmbc3LC16nTwPVOIxG5GlS7TH7TGqKGjZsiGbNmqFp06amrxUrVizJZYnIg5R0MVtrtUm7azRGuxFfoveRzaqb9q1hgShfrZpFxqnJyu9R98F74BNRx8Y/FRG5IvYpKuLaZ0eOHGGmiMiJ6Fl6RHW9fq4Dfjl6yZSVqncxAWs/G6HOudy8FSoMehgYORII0FcATkSuwyHDZ+Z9igIDA7F9+3a4Ew6fEbn+0iPG483PHsGE9bPRNuEAvKHhRuWqCEg8A3h7s98RkZuxS1BkrU9Rz5493bJPEYMiItesTSpoDbXKV5PR+8gWhASUwXNLPsIvBxPVddLOJWH2oklY3aAd/mjTFU8O7panxonBE5FrsEtQJH2JjH2KBg4c6PZ9ihgUETm3gmaoPTRrW6HPf757Q0xdc0Rlk+7ftwbvLp9qOrajZiQuLluNXs0Ny4mwWSSR67DLgrDsU0REzsS4rIisgyZfjVmbwrpnG32x+bhpeG39LW3w8h1PY0vt5sjy8sblgPKYtOywCrwkIBr59Q6UOXnC4vnGZpFynIhcj+5MUfv27VUAFBUVpR7NmzdHQK6CRPYpIiJnpDdTZE1Y2mW1FtuJ0BqY+3gMxi7ag4i92zFv/r+xr2o9LG8ci2WNYnGqYrWcwu6XuqqgjENsRG44Jb9Pnz7Yu3cv3nvvPRw9ehReXl5o0KCBKUiSgEmG0vr27aseRETO1j27oBlqIYG+uHz9Zr7PvxhUQT3E1mNJasis14XjKoPU7NxR9ehzeBPuGvKhur4cl0Ao5XoGh9iIXIju4bOXX34Z3333nZqCv2XLFlVILT2JJDj68ssvVXF1tWrV1DAaEZEzkcyMBCIid47GuD20Q4TOqxme8UWbe9B21DcY13MUNka0xI+33m5x1o51f2D/0y8h4PhRi/0cYiNyXsWqKRoxYoTq2yNB0ty5c1Vvop9//lkFRUOHDrX9XRIRlZBkZvJbQ022Zf+org10rcUmtUpGyWVDMD+qFwYNfF0tWGvhu+8w9rdvsH7WU1gxexQe2/Wz2m3MVHE9NiLnU6yO1ocOHVJDZubuvPNO/O9//8MHH3xgq3sjIrJr92zJJkkWx8tKvyM53u6WsEKH4ioG+WJPSA38WrcVbju5B7deOIFGF06azlFDbJevq/uQIIt1R0TOoVhT8jt16oQuXbrgtddey1NoLUNqaWlpcEcstCZyf3qm2hfWLHJYhwh8vtkwMy3k+hXc8edWHKxaDweq1jOde/fBDXjj0E84d0dfvOJ3K7b5V7X6ekTkxH2Kdu7cia5du6Jfv354/vnnVSCUkZGh+hX99NNPOHky519E7oRBEZFn0JO5KSh4Cgn0K3S224zFb6mmkUZzonrj5Z4j83TiZmBE5OQLwkrn6t9//x2jRo1Sw2i+vr7Izs5GmTJl8Pnnn5fkvomInKbfUXGH4iSoKmyI7f0HxuKPQ1vRYfd6xJ6IQ1z1Rqbj8pxalxPx9WfL0eP9YZzaT2QnJV777NSpU4iLi4O3t7cKlqTY2l0xU0REehU2xPZc9wb4YM2f6vvgG1eR4VMGN3xzisAnrvkEQ3f+hGv1GuJMjz6YENgcv/tVLnCIjYETkQMyReZq166tHkRElHe2W+4htvB/gpn0zGzTvtSAcnme75+ZgXSfMih79AjqHz2Cdh0exu+xD+eZ2m8cYuOyI0ROkCnyJMwUEZE912OTLtr3JOxEpz2/4r+dhuBopVqmY/32r8UtyWfwe5uueOjxuzBy3u48Q3WsTSKy44KwnoZBERHZMliKnbKu0Kn9yWn5d9le8vUYRJ09or4/GVYDb90+GKsa3pbvdcyXHTG+NofZyJOk2mNBWE8hjSojIyPRtm1bR98KEXlQl+37omrk/2RNw+w292BVg3ZI9/FFnYuncb2Mv8UpQenX1Hnmy44IGWaTYEyyVM/Oj1NfZZsdtomYKSoSZoqIyNZKOrVfgp+uR3dgRaMOyPTJKROdvmQympw7pvYvaxyL4c/2h7+vj6pD4jAbeZpUDp/ZHoMiIioN1oazijvEViYrE9unD0Lo9VTTviPv/A+Dbza0CL4KGmbjEBu5E7vNPiMiotLpi2QcYito6ZE37mmK15cdsgicJGMU+3+fqwxSn8ObcPvJOFy6vSvO/vCX6fld/9qOpKAK2BveAPDyshhmS7mewZls5JGYKSoCZoqIyBEKm25fWE+kmQ80wQ1vX1VDZDhJw6aZw1Az9QISQqpieaMO+LrVXTgdUkUtUfLF5hMcYiO3wkwREZGHLGRbWE+knk2rqRYARuUzrmF39cZqeK1Wyjk8tf0HLG7SRR1bEncm36E62SevJq8h98IhNnJHDIqIiNxg6ZHCAif53rj0yBX/IIy+5yUE3LyBzsd2IjrhAOIrRyBU1SdlqPM/XjoFieXDsLxRLHbLEiQcYiMPwOGzIuDwGRG5ssKG2WTo7PPNJxCemoRtM4aYjp8uXxmDHviPqXkkh9jI1bBPERERWTAOs8mwmjnZlv3dI8PV9qWywXjyvglYEtkJV/0CVWfthAqGY8YhtjrJp+Gl5SxVIoxBkmSQZGhNyFcZulsad1p9Ne4nckYcPiMi8iAFDbNJwGIYYgNWN2yvHrIGW/2kU8go42tqAZB2+Qp++/JZpPoHYWWj21QfpJ01boXm5c0hNnJpzBQREXlofdI9UTXUV2PdUX5dttPL+OFAeH2LLtsNLiYgy8sb1a5exNCdP+G7ueNQKS3F4jV+OZiohupy90UyLmTLDtrkjBgUERFRkYbY9ofXR5vRc/F4/1fwfZMuWNMgBhfKVcw5WdNQ7b230DZhP7yzswodYiNyFiy0LgIWWhORpyhJl+3bUk5g7sxRavt8UEWsaHQbJnceihu+loHWvOHtCpxRR2QrLLQmIiK7DLEZGbdjm9bCoqbdkOIfhCppl9S0/xu5FqyVbJIEXILF2OQsWGhNRERFUlizSLWQbYIPfLNuosOJOATeTFd9jowqpV3Cz18+C/x9D7Z3uxPPnw3G6Ss567exGJschcNnRcDhMyKikg+xPRy3Am+tmm7aPlGhGro8+YmavSbY74hsjcNnRETklENsC5v1wPb/zcWPbXrjUkB57KnW0BQQiTL/ZJjeWLLPYiiNw2xU2jh8RkREdh9iywr0wzMnQ1Cm81MITk+zeG7siTh8sWgSLv4YjKT996Dq8MFYWTWSPY+o1DEoIiIiuzeKlGyPyPQpg+SyIRbPC7uWguTAYIRdTwUWfINLRw5jxB0T8wzFGXsecZiNbIVBERER2X0hWwmQrFnUrDsWN+mCdqf24T0cwWeZVSwCogYXTuLxP5ZieaMO2FqnhcogSfBlrGeytiiukZ5zyDMxKCIiIruTQMSwpEj+xdjZ3j441qIdjg54Cp999rvFsb6HNuLBvavVQ2qS5rS8E9sfiNK1rIh00uYwHFnDQmsiIrI7Pf2O5HjS1fQ8z11fry3mRPXGhbIVUPHGFQRkplssKyLrtUk7gNzLisiDS49QQTglvwg4JZ+IyLYKy9zILLOHZm3L97myhEhMwgGcLR+G1Fp1kZyWofY/GLcS4zd8gdUN2mN54w7YHBGF0IrlVLiVmGoZEJkHYlIEvumlrhxK8+DPbw6fERGRUxZjFzbMJkNs2+o0R8UgX1NAJGJP7kFIehoG7F+jHrI+27/u+leB9yHXlsBM7oNLj3gujxo+S0hIQOfOnREZGYnmzZtj4cKFjr4lIiKPZ63fkd5htvuialjsf6bvWAx8aDK+atVHrb22tn6MxfFbzx9Dt79+h19mThdtIy494tk8avjs7NmzOHfuHKKiopCYmIjWrVvjyJEjCAoK0vV8Dp8RETnfMJtaVqSAITYJnrK8fUz7/rviI1WknepXFmsaxGBei57YUaupaZFaPQXb5Fo4fJaPatWqqYcIDw9HpUqVkJycrDsoIiIi5xtmkyxOQUNsql4oWBak9cK51Bs4Vy4MieVCEX41Gf0OrMfBKnXxR62mqqboUloGRn67iz2RPJRTDZ9t3LgRffv2RfXq1eHl5YUlS5bkOWf69OmIiIhAQEAAYmJisH379mK91s6dO5GVlYVatWrZ4M6JiMhZlxURr93dBK/dbThnasdH0P7pL9H/kbcxu/XdWNkoVu1/pc+teH3ZQQzbvhgf/vgOeh7ZAn9ZzPafmiMhGSQOpbkvp8oUpaWloUWLFhg2bBj69euX5/iCBQswZswYzJw5UwVEU6dORc+ePREfH48qVaqoc2RoLDMzM89zV69erYItIdmhQYMGYdasWQXeT3p6unqYp9+IiMj1lhUxZnfMz9lZM1I9JMs0459hONk/YN8aNE46iXsO/YqrfoGY2P3/8H2zbnmKsdkE0v04bU2RZIoWL16Me++917RPAqG2bdti2rRpajs7O1tlekaPHo1x48bpuq4EOT169MDw4cPx2GOPFXjua6+9hkmTJuXZz5oiIiLnVJKO1lJU/ey83Wh5Jh53xm9C7/jNqJl6Afc/MgV/1Gxiev6nnSoDlcIwcc0J1h25WU2RywRFGRkZKFu2LBYtWmQRKA0ePBiXL1/G0qVLC72m/KgPP/wwGjVqpAKewuSXKZIgjEEREZH7ydMTSdPQ4uwR7K3WAJpXTrXJhnVvo3LcdtVEclmjDlhfrw1u+AaYhupy1x0xo+R4bldonZSUpGqAqlatarFftg8fPqzrGps3b1ZDcDId31iv9M0336BZs2b5nu/v768eRETk/vL0RPLywp7qjUzHJYypHuQDr+PHEHTzBu46/Jt6PHvXv7C0SRf1HDnHfC02LiviWlwmKLKF2NhYNeRGRESUm7FgW2aZSXBjPoxizOs80P4WdLo6A80T/8Sdhzeh29EdFn2Q5Dm3b1yKo0EncLJdZ4z44TBnsrkQlwmKZPq8j4+P6jNkTrZlej0REVFpF2ynZ2arDNLeag3V479dhlk83yc7Cy9s/BqVVqagtq8/ptdtjfc7Poq/KtU2nZNfRolDbM7BZYIiPz8/1Wxx7dq1ppoiyfrI9qhRo0r1taUNgDxk+I6IiDy3J5LUHRUk4GY6vmveA48mbEfw6VO488gWTOk8xOIcLy1b1SgZZ7KxWaTzcKpC66tXr+Kvv/5S37ds2RLvv/8+unTpgtDQUNSuXVvVA0lh9SeffILo6Gg1Jf+7775TNUW5a41KAztaExF5NsnoxE5Zl2+jSPOFZV/s2QifffwDYk7tx+y291ic8/03Y5EUVAHLGsei+iMD8MnuC3mulV/RNrNJHjb7bMOGDSoIyk0CoS+//FJ9L9Px33nnHbVMh/Qk+uijj9RUfXtgUERERFI8LTVBsFJ3JIGMtaVHal1OxG+fPGHaTi/ji7sHfYD4yhFWA6xNL3XFLwcTmU3ytKDI2TEoIiIiUdisMqsZJU3DrReOo8/hzbjryGYEXb+KmJFfqeVIjNqf3It94fVx1b+s2n6+e0NMXXNEVzaJ8segqBQwKCIiIqPChrMKyygNu60Ofly9GxfKhZqOlUu/hp0fP6KesbFuKyxvFIstUZ1xLjP/VbnMs0kcSiv557dTrX3mrKTIOjIyUnXTJiIiKmgtttwz2SRoMSfbsr97k2oWAZGokXoeCSFV4Z+ViR5/bceUFR/h+tVrVu/BfOkRY6AmxeDSnVu+cp22omGmqAiYKSIiIltllAoaYmuYdFINsVXJTMP4Lk+aDnlnZ+G9ZR9gU0RL/NIgBqkB5dT+Dx+Mgn8Zb9YdWcHhs1LAoIiIiGypsCG257o3wAdr/jTtj07Yj+++Naz1meFdBpsiovD0veMwondz1h0VgMNnRERETq6wIbZRXRuobI8xuJGhtQ86PIz4SrXhl52J6qkXUCEsBPO2nzIFRFKXZGTcJxkk41Aah9isY6aoCJgpIiIiexdtW8sm1U9KQNj1FNw29D5TNqlmyjms+/QpbKnTAssad8AvDdrhcqDh82re8HYe2ygylcNnpdPR+siRIwyKiIjIaVoAyNIjz86PU/se2LMab6/8yHTOhaAKiB75teqgPaxDBL7YfEL3EFuWGzWLZFBUCpgpIiIiR7EWpMgQmHmjyFsu/o3e8ZvRJ34T9oY3wLjez6j9oUF+uHEpBfcc/BWrGrZHctkQq1P7VxbSh8nVMCgqBQyKiIjIlZYe8c26iUwfX1QM8kVy2k30OfQbpv84BZle3thWuxmWRnbCwuZ3mM6f988QmwzXuVPRNgutiYiIPIBkdiSDI3IPbklAJO6LqqG+3vQpg/1V66GMlo3Yk3vQ99BvFucnplxXGaL8siX5FW27mzKOvgEiIiKyzSy23ENeMiQmAZOsxfb55hNY3bC9etS5dAZ3xm/Gn2G1La7js2M73ps5Ccsbx2Jlw/ZICqpotVmkNKx0p7ojweGzIuDwGRERuVWjyFw1RV8f/A4N5nxquJaXNxY2626qSTLnas0iOXxGRETkYawtPVLQEJtxe2LfSFx54v/wVuehiKvWAD5aNs7lWoYk+MZVVL56CSeSrqm6I/OASEjQJfulUNvIlfoiMVOkA6fkExGROyhsVlmWWUapRso5pPv44UK5nCG0YTuW4uX1nyGuTjMsqd8eKxp1yLN+m/lMtl8OJjpFNomzz0oBh8+IiMjVFVYHtLKApUdeXzUdj8atMO1b1LQbxvZ5Pt/Xeb57Q6dZekTv5zcLrYmIiDxwiK04RduVvvkcqxMS8Pu7s3Dn4U34uXGsxXNlbbYm545ieaMO+GKzr9VZbBIYyfV7RIabap6coWCbmaIiYKaIiIg8RZbOZpHmPvrxbdx9aKP6/o8at2JW9H1Y1fA2q69hr6VHWGhNRERENi/ajq4barFIrbkttZtjRw1DQXeb04dQOe2yxXEvLdtiW2qO9BZs2wMzRUXATBEREREKrDsSE6KC8fdnc/DzrR0teh29tOFLxCTsw7JGsVjRuANuVKuJ5LSMfF8j99IjJcFMEREREZWKXv/UHUnQYk62Zf/QB2KxqtsDuGje/FHT0OvIZrQ6E49X1n+OLTOG4a6N31t9DfNGkfbCQmsiIiIqVmAkhdLWCqSlJkiySV7GbJKXFwY+9F/0PrJFddOO/vsAdtY0DLUZBdy8gRu+loGWXNtemCnSQXoURUZGom3bto6+FSIiIqevO7KWTTpfPgyruw7ApeWrsXPLfhyocgvMlcu4jtwk2LIX1hQVAWuKiIiISmfpEZ/sLPhm3TRlihxRU8ThMyIiIrJ7XyTj0iPGIbYsbx/1yL30iD37FXH4jIiIiJyyYNveC8syU0REREROW7BtTwyKiIiIyKmXHrEXDp8RERERMVNUNMaJelLFTkRERK7B+Lld2IR7BkVFcOXKFfW1Vq1ajr4VIiIiKsbnuEzNt4Z9ioogOzsbZ86cQfny5eHllX8BmDR43LFjh9Vr5HdcIlgJtBISElyu/1FhP68zvk5xr1XU5+k9X895xXlfufJ7y5PeV458b/F95Rqv5W5/s9o64H0loY4ERNWrV4e3t/XKIWaKikB+kTVr1izwHB8fnwL/YxZ0XPa70h8YPT+vM75Oca9V1OfpPV/PeSV5X7nie8uT3leOfG/xfeUar+Vuf7N8HPS+KihDZMRCaxsbOXJkiY67Gnv9PLZ8neJeq6jP03u+nvP4vnL+1ynJtRz13uL7yjVey93+Zo104vcVh8+cAJcPodLC9xaVBr6vyF3fV8wUOQF/f39MnDhRfSWyJb63qDTwfUXu+r5ipoiIiIiImSIiIiIiAwZFRERERAyKiIiIiAwYFBERERExKCIiIiIyYFBERERExKCIiIiIyIBBERERERGDIiIiIiIDBkVEREREDIqIiIiIDBgUERERETEoIiIiIjJgUERERETEoIiIiIjIgEEREREREYMiIiIiIgMGRUREREQMioiIiIgMGBQRERERMSgiIiIiMmBQRERERMSgiIiIiMiAQRERERERgyIiIiIiAwZFRERERADKOPoGXEl2djbOnDmD8uXLw8vLy9G3Q0RERDpomoYrV66gevXq8Pa2ng9iUFQEEhDVqlXL0bdBRERExZCQkICaNWtaPc6gqAgkQ2T8pQYHBzv6doiIiEiH1NRUldQwfo5bw6CoCIxDZhIQMSgiIiKyjaxsDduPJ+P8lRuoUj4A0XVD4eNt+zKVwkpfGBQRERGRw6zcfxaTfjqIsyk3TPuqhQRgYt9I9Gpaza73wtlnRERE5LCAaMScXRYBkUhMuaH2y3F7YlBEREREDhkykwyRls8x4z45LufZC4MiIiIisjupIcqdITInoZAcl/PshUERERER2Z0UVdvyPFtgUERERER2J7PMbHmeLTAoIiIiIruTafcyy8zaJHnZL8flPHthUERERER2J32IZNq9yB0YGbfleGn0K7KGQRERERE5hPQhmvFoK4SHWA6Rybbst3efIjZvJCIiIoeRwKdHZLhdOloXhkEREREROZQEQO3rhTn2Jjh8RkREROQGQdH06dMRERGBgIAAxMTEYPv27VbP/eGHH9CmTRtUqFABQUFBiIqKwjfffGPX+yUiIiLn5bJB0YIFCzBmzBhMnDgRu3btQosWLdCzZ0+cP38+3/NDQ0MxYcIEbN26FXv37sXQoUPVY9WqVXa/dyIiInI+Xpqm2W9RERuSzFDbtm0xbdo0tZ2dnY1atWph9OjRGDdunK5rtGrVCn369MHrr7+u6/zU1FSEhIQgJSUFwcHBJbp/IiIisg+9n98umSnKyMjAzp070b17d9M+b29vtS2ZoMJIHLh27VrEx8fj9ttvt3peenq6+kWaP4iIiMg9uWRQlJSUhKysLFStWtViv2wnJiZafZ5EiOXKlYOfn5/KEH388cfo0aOH1fMnT56sIkvjQzJRRERE5J5cMigqrvLlyyMuLg47duzAm2++qWqSNmzYYPX88ePHq0DK+EhISLDr/RIREZH9uGSfokqVKsHHxwfnzp2z2C/b4eHhVp8nQ2z169dX38vss0OHDqlsUOfOnfM939/fXz2IiIjI/blkpkiGv1q3bq3qgoyk0Fq227dvr/s68hypGyIiIiJyyUyRkKGvwYMHq95D0dHRmDp1KtLS0tQ0ezFo0CDUqFFDZYKEfJVz69WrpwKh5cuXqz5FM2bMcPBPQkRERM7AZYOigQMH4sKFC3j11VdVcbUMh61cudJUfH3q1Ck1XGYkAdPTTz+Nv//+G4GBgWjcuDHmzJmjrkNERESlIytbc4p1zdy6T5EjsE8RERGRfiv3n8Wknw7ibMoN075qIQGY2DdSLQRrL27dp4iIiIicPyAaMWeXRUAkElNuqP1y3NkwKCIiIiKbD5lJhii/oSjjPjku5zkTBkVERERkU1JDlDtDZE5CITku5zkTBkVERERkU1JUbcvz7IVBEREREdmUzDKz5Xn2wqCIiIiIbEqm3cssM2sT72W/HJfznAmDIiIiIrIp6UMk0+5F7sDIuC3Hna1fEYMiIiIisjnpQzTj0VYID7EcIpNt2W/PPkVu39GaiIiInFuvptXQIzLcZTpaMygiIiKiUiMBUPt6YXAFHD4jIiIiYlBEREREZMCgiIiIiIhBEREREZEBC62JiIioWGRBV1eZWaYHgyIiIiIqspX7z6qV7s0XfpUu1dKU0Rl7EOnB4TMiIiIqckA0Ys4ui4BIJKbcUPvluCsq9UxRxYoV4eWlL5WWnJxc2rdDREREJRwym/TTQWj5HJN98okvx6Vpo6sNpZV6UDR16lTT9xcvXsQbb7yBnj17on379mrf1q1bsWrVKrzyyiulfStERERUQtuPJ+fJEOUOjOS4nOcqTRvtFhQNHjzY9H3//v3xn//8B6NGjTLte+aZZzBt2jSsWbMGzz//fGnfDhEREZXA+Ss3bHqex9YUSUaoV69eefbLPgmKiIiIyLlVKR9g0/M8NigKCwvD0qVL8+yXfXKMiIiInFt03VA1y8xatZDsl+Nynqux65T8SZMm4YknnsCGDRsQExOj9v3+++9YuXIlZs2aZc9bISIiomLw8fZS0+5llpkEQOYF18ZASY67WpG13TNFQ4YMwebNmxEcHIwffvhBPeT7TZs2qWNERETk/Ho1rYYZj7ZCeIjlEJlsy35X7VPkpWlafrPqKB+pqakICQlBSkqKCuaIiIg8WZaLdLTW+/lt947WR48exRdffIFjx46p6fpVqlTBihUrULt2bTRp0sTet0NERETF5OPt5XLT7p1m+OzXX39Fs2bNVB3R999/j6tXr6r9e/bswcSJE+15K0RERESOC4rGjRunmjf+8ssv8PPzM+3v2rUrtm3bZs9bISIiInJcULRv3z7cd999efbLEFpSUpI9b4WIiIjIcUFRhQoVcPZs3kXidu/ejRo1atjzVoiIiKiQIuqtRy9iadxp9VW23Z1dC60ffPBBvPTSS1i4cKFaJDY7O1tN0R87diwGDRpkz1shIiIiK1buP6sWdTVf40waMkr/IVedbu90maK33noLjRs3Rq1atVSRdWRkJG6//XbcdtttePnll+15K0RERGQlIBoxZ1eeRV8TU26o/XLcXdmtT5G8TEJCAipXrqzqh6S+SAKjli1bokGDBnAF7FNERETuLCtbQ+yUdXkCIiOvfxo0bnqpq1P2I3KZPkUSFNWvXx8HDhxQQZBki4iIiMh5bD+ebDUgEpJFkeNynjv1J7L78Jm3t7cKhi5evGivlyQiIqIiOH/lhk3PczV2rSn673//ixdeeAH79++358sSERGRDlXKB9j0PFdj16BIZpht374dLVq0QGBgIEJDQy0eRTV9+nREREQgICAAMTEx6trWzJo1Cx07dkTFihXVo3v37gWeT0RE5Gmi64aqWWbWqoVkvxyX89yRXafky1pntrJgwQKMGTMGM2fOVAGRXLtnz56Ij49XzSBz27BhAx566CE1002CqClTpuCOO+5QNU7skURERARVPC3T7mWWmQRA5jOxjIGSHHelImunnH1maxIItW3bFtOmTVPb0vNIirdHjx6tlhMpTFZWlsoYyfP19kji7DMiIvIEK92sT5HTzD6TGzHegHxfEL2BRkZGBnbu3Inx48dbFHLLkNjWrVt1XePatWu4efNmgcN26enp6mFU2P0TERG5g15Nq6FHZLiaZSZF1VJDJENm7pohsltQJNkYWdpDhrRkmQ/pZJ2bJKtkv2Rv9JA+R3Ju1apVLfbL9uHDh3VdQzprV69eXQVS1kyePBmTJk3SdT0iIiJ34uPt5ZbT7h0aFK1bt86UjVm/fj2cgcyCmz9/vqozkvoiayQTJXVL5pki9lciIiJyT6UeFHXq1Cnf70uiUqVK8PHxwblz5yz2y3Z4eHiBz3333XdVULRmzRo0b968wHP9/f3Vg4iIyJ26VnvasJhTzj7buHFjgcdlHTQ9/Pz80Lp1a6xduxb33nuvqdBatkeNGmX1eW+//TbefPNNrFq1Cm3atCni3RMREbk2dyugdumgqHPnznn2mdcY6a0pEjKsNXjwYBXcREdHqyn5aWlpGDp0qDouM8pkqr3UBQmZgv/qq6/i22+/Vb2NEhMT1f5y5cqpBxERkScs9Jp7yrlxodcZj7by+MDIrs0bL126ZPE4f/48Vq5cqabWr169ukjXGjhwoBoKk0AnKioKcXFx6lrG4utTp06pAm+jGTNmqFlr999/P6pVq2Z6yDWIiIjcfchMMkT59eAx7pv000F1nidzij5Fv/76q8r8yDR7Z8Y+RURE5Iq2Hr2Ih2ZtK/S8ecPbueWMM72f33bNFFkj2R3pRE1ERES25+kLvTplTdHevXsttiVJJUNcMhtMhsCIiIjI9jx9oVenDIok8JHC6twjdu3atcPs2bPteStEREQet9CrFFXnVzMjU57C3XihV6cMio4fP26xLUtzVK5cucAGikRERFQynr7Qq1MGRXXq1LHnyxEREdE/ZLq9TLvP3adIMkTsU+SAoOijjz7Sfe4zzzxTqvdCRETkaTx1oVennJJft25dXLhwQa1QL4vDisuXL6Ns2bJqGM10U15eOHbsGJwNp+QTERG5Hqecki9LbEix9aFDh5CcnKwe8n2rVq3wxhtvqJojeThjQEREROTspPmi9CRaGndaffX0ZoxOnSmqV68eFi1ahJYtW1rsl6aN0mk6dyG2s2GmiIiInBXXNXOxTJH0JMrMzMyzX9Y8y73iPRERERVtXTPzgMh8XTM5ToWza1DUrVs3PPXUU9i1a5dFlmjEiBHo3r27PW+FiIjILXBdMxcNiqRBY3h4uFrZ3t/fXz1khXtZ5uOzzz6z560QERG5BZlJljtDZE5CITku55ETTcmXGWbLly/HkSNHcPjwYbWvcePGaNiwoT1vg4iIyG1wXTMXDYqMIiIi1FIfUnhdpoxDboGIiMgtcF0zFx0+k/5Ejz/+uOpL1KRJE5w6dUrtHz16tFoUloiIiIq3rpm19ouyX457+rpmThcUjR8/Hnv27MGGDRss1juTIusFCxbY81aIiIjcal0zkTsw4rpmThwULVmyBNOmTUNsbKzqWm0kWaOjR4/a81aIiIjcbl0zWcfMnGzLfk/vU6SXXQt6ZImPKlWq5NmflpZmESQRERGRJZlSX9CaZVzXzMWCIpmKv2zZMlVDJIyBkEzHb9++vT1vhYiIyO26VUsA1L5emIPu0vXZNSh666230Lt3bxw8eFB1tv7www/V91u2bMGvv/5qz1shIiJyqW7VuVsvGrtVc3jMRWuKpJZICq0lIGrWrBlWr16thtO2bt2K1q1b2/NWiIiInB67VbtppujmzZtqiY9XXnkFs2bNstfLEhEReUS3ag6buVCmyNfXF99//729Xo6IiMjlsVu1Gw+f3XvvvWpaPhERERWO3arduNC6QYMG+M9//oPNmzerGqKgoCCL488884w9b4eIiMglulVLUXV+VUNe//QiYrdq2/DSZBEyO6lbt671G/HywrFjx+DMUlNTERISgpSUFAQHBzv6doiIyINmnwnzD2xj9yHOPrPd57ddM0XHjx+358sRERG5fGNGY7fq3H2KwvPpU0QlwyXqiYiInLwxI7tVu+Hwmavj8BkREdmjMSOHxhzz+W3X2WdERERkwMaMzodBERERkZM3ZiT7YFBERETkAGzM6HwcVmidlpaGBQsW4Pr167jjjjtUDyMiIiJPmVnGxoweGhSdOnUKjz32GHbt2oV27drh888/R48ePfDnn3+q44GBgVixYgVuv/12e9wOERGRw2eWyWwyNmb0wOGzsWPHIiMjAzNnzkTZsmXRs2dPlRk6e/Yszp07h969e+O1116zx60QERHZdWZZ7rohCYJk/y8HE1VwJHJPrDduy3FOu3ezKfnh4eH48ccfER0djeTkZFSqVEkt9dG+fXt1fM+ePejWrRuSkpLgzDgln4iI9A6ZxU5ZZ7WQ2pgF2vRSVxUcFdaniNyoo/X58+dRp04d9X1oaKjKFlWtWtUiaLp06ZI9boWIiMipZpaxMaMHzj6Ttc3y+74kpk+fjoiICAQEBCAmJgbbt2+3eu6BAwfQv39/db68/tSpU21yD0RERCWdWSYBUPt6Ybgnqob6yoDIzWefvfrqqypDJKS+6M0331SpLHHt2rUiX09mro0ZM0bVKUlAJEGO1CrFx8ejSpUqec6X17jlllswYMAAPP/88zb4iYiIiPLHmWWuyS41RZ07d9aVHVq/fr3ua0og1LZtW0ybNk1tZ2dno1atWhg9ejTGjRtX4HMlW/Tcc8+pR0HS09PVw3xMUl6DNUVERKSnpqiwmWVSU8SskIfVFG3YsMGm15NM086dOzF+/HjTPm9vb3Tv3h1bt2612etMnjwZkyZNstn1iIjIM1a2l69SKC2zzGSPeWDEmWXOy2HNG0tCZqllZWVZFGsL2T58+LDNXkeCLhmiy50pIiIiz6Z3ZXtZ0DX3eZIh4swyDw2KJKh4/fXXERQUZBFg5Of999+HM/H391cPIiKiwla2N/YfMl/ZnjPLXEupB0W7d+/GzZs3Td9bU5QZadLnyMfHRzV+NCfbMr2fiIjIESvbyyeZHJdAyHwoTWaUkfMr9aDIvHi6KIXUBfHz80Pr1q2xdu1a3HvvvaZCa9keNWqUTV6DiIioJP2HGAi5HpesKRIyFDd48GC0adNGdcqWKfmyyOzQoUPV8UGDBqFGjRqqWNpYnH3w4EHT96dPn0ZcXBzKlSuH+vXrO/RnISIi1yii5sr27s0uQdGwYcN0nTd79mzd1xw4cCAuXLig+h8lJiYiKioKK1euNBVfyyK0MiPN6MyZM2jZsqVp+91331WPTp062Xx2HBERuWcRNfsPuTe79CmS4ESW+ZCgpKCXW7x4MZwZ1z4jIvLMImpj1ev0h1vi9WWH2H/IxThVn6IRI0Zg3rx5OH78uBreevTRR9UaaERERK5URC0B0St9IjHyW/Yfckd2WftM1ig7e/YsXnzxRfz000+q188DDzyAVatWFZg5IiIicrYi6opBfmravWSEzMm2+XR8cj12K7SWfj8PPfSQepw8eRJffvklnn76aWRmZqrFWqXgmYiIyBWKqGXhVvYfcj8OmX0mNUbSl0iyRNKZmoiIyNWKqNl/yP3YZfhMyMKqUlfUo0cPNGzYEPv27VOLucosMWaJiIjIXkXUuYfIjJ2oL6WlqwDJWq5H9stxyQiRe7JLpkiGyebPn69qiWR6vgRH0pWaiIjIHlhETU41Jb927dpqSn5By3n88MMPcGackk9E5Jr1QluPXsRDs7YV+vx5w9sh5XpGoYu9kmtxqin50l26KGubERER2bJeKD0zW9c1WETt2ewSFMlMMyIiIketXP9c9wa6rsMias9mt0JrIiIiR9QLiXnbTyE8mEXUVDAGRURE5DLBj9QGLY07rb7Ktt6mi4mp6Xgourbazh0YsYiaHNqniIiIyN71QhGVyqqO07mvI52oWURNgkERERF5TL2Q1AmxiJqsYVBEREROOY1eb38hY73QudSCV6431guxiJqsYVBEREROOSwmw1l664We794QU9ccYdNFKhEWWhMRkVMuuyHH9S7SaqwX4sr1VBLMFBERkd2HxvQMi8nxdwe00PU6rBciW2BQREREdh8aCwn0K3RYTB3XDM+R7BHrhai0cfiMiIjsPjS25mCiruskpaWrIEqwvxCVNgZFRERk02aKejpML447res1ZAhM6oFYL0T2wOEzIiIq0jR5WwyNJafdRGiQHy6lZegaFpPAh/VCVNoYFBERUZGmyRfWTHFYhwhdr3NvVHV8sfmE7mn0rBei0sbhMyIiD2Nt2EtPLdDyvWdsNjQmmR8Oi5EzYaaIiMiDFJQFkiClsGnyLy/dr4a+bDU0JtkfDouRs2BQRETkQXVAha0hpqcWyNZDYxwWI2fBoIiIyM0Ln2UYSk+zRAlibEWyP3KPXJGeXAmDIiIiFw94CssASX2Onhlhl6/rywJxaIzcFYMiIiIXDnimP9wSry87VOhyGS/2aqzrfisE+iLl+s0CA55X+kRi5Le7ODRGbodBERGRCwc8egqf5fWTr6br+pmGdqhb6GrzqpmidysOjZHbYVBERFTEYMeZAh69hc8y5KVnDbFRXeujUXi5QgMeNlMkd+SlaVp+/39QPlJTUxESEoKUlBQEBwc7+naIqBjBTEmzO8Zz8gt4jFcxBjzWanjkvIpBvroDGluYN7wdUq5nqPuGlSyQeW8gPb9LInf7/GZQVAQMiohcP5gpTnbHPGiQ7EjslHVOFfDoKXze9FJX9fvS83sicjcMikoBgyLyFLYIUvSeY69gRpQ0uyPBxbsDWuCRz36HMwU8xsJnPRkgwSwQeZpUBkW2x6CIivqhUth5trqOswUpRTnHHsFM1WB/9V1iasmzO6O61Me09X/B2QIeZoCIrGNQ5CJBUUEfUkX515wtPhCd8fXsHTToOU/vh4+e7IYtruNsQYrec/QMQ9kymLGVUV3qYdr6o04Z8DADRJQ/BkUuEBQV9IdO6P1Xn62GH5zt9ewdNOi978I+7PWc9+TtdfHpxuMlvo4zBil6znHEMJStzH08BmMX7Sl0FhcDHiLn4RFB0fTp0/HOO+8gMTERLVq0wMcff4zo6Gir5y9cuBCvvPIKTpw4gQYNGmDKlCm48847HRIUFfRhZ+0/SH71Abb40NT7AW3P1yvsGrYOGvRcS2/Nya8vdEGnd9YX2D1YPsvMFiYv1nVcPUix5TCUvQuWfzmYqGsWFwMeIufg9kHRggULMGjQIMycORMxMTGYOnWqCnri4+NRpUqVPOdv2bIFt99+OyZPnoy77roL3377rQqKdu3ahaZNm9o1KJI/gAX9i7wg5n+YRUn/Za/3A9qer6f3Q9xWQYPea+kdpnmlz60qeCopW13HWYMUvcNQtghmjO+Dc6m2ye4IBjxErkPv57fLNm98//33MXz4cAwdOlRtS3C0bNkyzJ49G+PGjctz/ocffohevXrhhRdeUNuvv/46fvnlF0ybNk09Nz/p6enqYf5LtQX5A1mcgMi8O61cQxS2llFiauFdbK0FKI56vcKuYbynb7aeKPH9FOVaeutWTiZf03Weva6jn33/fdT+lkr4ftfpAoehihrMWOvC/NrdTdRXye7YqlOz3uaFXOaCyHW4ZFCUkZGBnTt3Yvz48aZ93t7e6N69O7Zu3Zrvc2T/mDFjLPb17NkTS5Yssfo6klWaNGkSbE3+gDrDNZz59ewdNNjyWnVCyzrVdewdpOgNZNrVC1OBRkGBiq2DGcn02LJTMwMeIvfikkFRUlISsrKyULVqVYv9sn348OF8nyN1R/mdL/utkaDLPJCSTFGtWrVKfP/yR9YZruHMr2fvoEHvtfTUnDzWPgKfbTpuNbgQ8vkqA9cluY6zBil6AxkJKFQw86j9ghlmd4jI7YIie/H391cPW5M/wgWtQVQQ4weYXEMUtpZRYR+Iej+g7fl6ej/EbRU06L2W3tXB/cp4FxpcDO9oKDYv6XWcOUjRc44jghkGPETkVkFRpUqV4OPjg3Pnzlnsl+3w8PB8nyP7i3J+aZI/ygV92Bm3C/sAM35f0g9EPR/Q9nw9vR/itgoa9F6rKMM0eoKLlrUr2uQ6zhqk2HoYisEMEZU2l519JjPOZPq9TMMX2dnZqF27NkaNGpVvofXAgQNx7do1/PTTT6Z9t912G5o3b2610Do39imy7+s5Y58iT+hoTUTkbjxiSv7gwYPxySefqOBIpuR/9913qqZIaoVkun6NGjVUsbRxSn6nTp3w3//+F3369MH8+fPx1ltvOWRKvjl2tHa9jtZERORa3D4oEjKd3ti8MSoqCh999JHKIInOnTsjIiICX375pel86WP08ssvm5o3vv322w5r3khERET24RFBkb3JL7NChQpISEhgUEREROQijLPHL1++rIIjtyq0dpQrV66or7aYlk9ERET2/xwvKChipqgIpJj7zJkzKF++PLy88q8xadu2LXbs2GH1GvkdN0awrpiBKuzndcbXKe61ivo8vefrOa847ytXfm950vvKke8tvq9c47Xc7W9WWwe8ryTUkYCoevXqqtmzNcwUFYH8ImvWrFngOdIqoKD/mAUdl/2u9AdGz8/rjK9T3GsV9Xl6z9dzXkneV6743vKk95Uj31t8X7nGa7nb3ywfB72vCsoQGVkPl6hYRo4cWaLjrsZeP48tX6e41yrq8/Ser+c8vq+c/3VKci1Hvbf4vnKN13K3v1kjnfh9xeEzJ8BZbVRa+N6i0sD3Fbnr+4qZIicgS4lMnDixVJYUIc/G9xaVBr6vyF3fV8wUERERETFTRERERGTAoIiIiIiIQRERERGRAYMiIiIiIgZFRERERAYMilzAfffdh4oVK+L+++939K2Qm5A2+p07d0ZkZCSaN2+OhQsXOvqWyA3IYptt2rRBVFQUmjZtilmzZjn6lsjNXLt2DXXq1MHYsWNL5fqcku8CNmzYoNZs+eqrr7Bo0SJH3w65gbNnz+LcuXPqwysxMRGtW7fGkSNHEBQU5OhbIxeWlZWF9PR0lC1bFmlpaSow+uOPPxAWFuboWyM3MWHCBPz1119qjbR3333X5tdnpsgFyL/oZRFaIlupVq2aCohEeHg4KlWqhOTkZEffFrk4WdNKAiIhwZH8m5v/7iZb+fPPP3H48GH07t0bpYVBUSnbuHEj+vbtq1bm9fLywpIlS/KcM336dERERCAgIAAxMTHYvn27Q+6VPPN9tXPnTvUvfPmXF3k2W7yvZAitRYsWavHsF154QQXcRBtt8N6SIbPJkyeX6n0yKCplkkKWPxDyHzs/CxYswJgxY1Rr8127dqlze/bsifPnz9v9Xsnz3leSHRo0aBA+/fRTO905ufv7qkKFCtizZw+OHz+Ob7/9Vg3TEqWV8L21dOlSNGzYUD1KldQUkX3Ir3vx4sUW+6Kjo7WRI0eatrOysrTq1atrkydPtjhv/fr1Wv/+/e12r+T+76sbN25oHTt21L7++mu73i+5/98roxEjRmgLFy4s9Xsl939vjRs3TqtZs6ZWp04dLSwsTAsODtYmTZpk83tjpsiBMjIy1NBF9+7dTfu8vb3V9tatWx16b+Te7yv5uzRkyBB07doVjz32mAPvltzpfSVZIZkUImSlcxkyadSokcPumdznvTV58mQ1a/bEiROqwHr48OF49dVXbX4vDIocKCkpSdVyVK1a1WK/bMuMICN5YwwYMADLly9X4/QMmKik76vNmzerdLWM60vBtTz27dvnoDsmd3lfnTx5Eh07dlRDH/J19OjRaNasmYPumNzts9AeytjlVahE1qxZ4+hbIDcTGxuL7OxsR98GuZno6GjExcU5+jbIzQ0ZMqTUrs1MkQPJrAyZwpq7EFG2ZZo0UXHwfUWlge8r8oT3FoMiB/Lz81NN89auXWvaJ/96l+327ds79N7IdfF9RaWB7yvyhPcWh89K2dWrV1X3TSOZpirp5dDQUNSuXVtNQRw8eLBqjS+p56lTp6qpi0OHDnXofZNz4/uKSgPfVwRPf2/ZfD4b5ZlKL7/m3I/Bgwebzvn444+12rVra35+fmpa4rZt2xx6z+T8+L6i0sD3FXn6e4trnxERERGxpoiIiIjIgEEREREREYMiIiIiIgMGRUREREQMioiIiIgMGBQRERERMSgiIiIiMmBQRERERMSgiIiIiMiAQRERERERgyIicgedOnWCl5cX5s2bZ7H/448/RvXq1R12X0TkWhgUEZFLk+Ubd+/ejWrVquH777+3OLZz5060atXKYfdGRK6FQRERubQ///wTV65cwcsvv4wVK1bg2rVrpmO7du1C69atHXp/ROQ6GBQRkUuTbFBAQACeeOIJBAcHq8BI3LhxA4cOHWKmiIh0Y1BERC5NskHNmzeHn58f7rvvPixatEjt37NnDzIzM4sVFJ05cwaPPPJIgedcvnwZn376abHvm4icD4MiInL5oMgY+PTr1w/Lli1Denq62l+5cmXUqlWryNeU4uy5c+cWeA6DIiL3w6CIiFyaed1Q586d4evri1WrVlkUWZ84cUJlkx544AHceuutGDx4sMoiiSlTpqBp06Zo1qyZKRCS89u0aWP6vkWLFuo58tyBAweq4u4JEybg4MGDiIqKwn/+8x9cvXoVvXr1UteRh9wDEbmWMo6+ASKi4jp27JjK2BiDnzJlyuDuu+9Ws9D27duH3r17m87dv38/Zs+erYKdhx9+GHPmzEGTJk3w3Xff4Y8//lAF2m3btkWXLl3yvI7UJsl0fwmK5PimTZvw5ptvIj4+Xj1XyGuGhYVh5cqVKmiS4m8ici3MFBGRy5JskNQSSabHqH///vjxxx9x4MABi3qi+vXrm7I/Dz74oApsNm/erM6XQu3Q0FB069YNO3bsyPM6jRo1QmRkpOqF1LJlS5U9yk2yQxs3bsSLL76Ibdu2qaJvInItDIqIyKWHziQgksDIqEePHsjKykJGRoZFUCQBjfn35tuF8ff3N33v4+Ojrp9bw4YNERcXp7JPY8aMwbRp04r5UxGRozAoIiKXNXnyZJUtyh3ApKamqiGsunXrWvQzkiBKLFiwALGxserxww8/qMLsS5cuYd26dYiOjtb12uXLl7cYIpMZa0FBQar26LnnnlMBEhG5FtYUEZFHkIySFFXv3btX1Q5JXZEUZQ8YMEAVakvmaNKkSaozdn7DY7lJ/ZBkomTYTK4RExODsWPHqkxSYGAgPv/8c7v8XERkO16a/HOKiMiNSZBz//33m4qiiYjyw+EzIiIiImaKiIiIiAyYKSIiIiJiUERERERkwKCIiIiIiEERERERkQGDIiIiIiIGRUREREQGDIqIiIiIGBQRERERGTAoIiIiImJQRERERGTAoIiIiIhAwP8D91WjxdEjLqsAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "Ntrees = np.geomspace(10, 10**4, dtype=int)\n", "\n", "fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)\n", "\n", "d2s = [do_one(NKD=NKD_, Nsample=10**3)[0] for NKD_ in Ntrees]\n", "ax1.plot(Ntrees, d2s, 'o')\n", "\n", "pf = np.polyfit(np.log(Ntrees), np.log(d2s), 1)\n", "xx = np.geomspace(np.min(Ntrees), np.max(Ntrees), 1000)\n", "ax1.plot(xx, np.exp(np.polyval(pf, np.log(xx))), dashes=[2,2], color='r')\n", "ax1.text(100, 0.07, rf'$y=\\exp({pf[1]:0.3f})N^{{{pf[0]:0.3f}}}$', color='r')\n", " \n", "ax1.set_xscale('log')\n", "ax1.set_yscale('log')\n", "ax1.set(ylabel=r'$d_{NN}$')\n", "\n", "MiBs = [do_one(NKD=NKD_, Nsample=10**3)[1]/1024**2 for NKD_ in Ntrees]\n", "ax2.plot(Ntrees, MiBs, 'o')\n", "ax2.set(xlabel=r'$N_{\\rm points}$', ylabel='MiB required');\n", "\n", "# del ptr" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.8" } }, "nbformat": 4, "nbformat_minor": 5 }