{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Runtime Benchmark (MNIST)\n",
    "\n",
    "This notebook compares compute costs of UMAP (=$k$-NN), the exact $k$-MST (=KDTree-based boruvka) and approximate $k$-MST (=NNDescent-based boruvka) algorithms. The dataset samples and generated graphs are stored for re-analysis and visualization. On MNIST, the approximate $k$-MST is roughly two orders of magnitude faster than the exact $k$-MST algorithm!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "from tqdm import tqdm\n",
    "from scipy.sparse import save_npz\n",
    "\n",
    "from sklearn.datasets import fetch_openml\n",
    "from sklearn.utils.random import sample_without_replacement\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\", category=UserWarning)\n",
    "\n",
    "from umap import UMAP\n",
    "from multi_mst import KMST, KMSTDescent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# Trigger numba compilation\n",
    "_ = KMSTDescent().fit(np.random.rand(100, 2))\n",
    "_ = KMST().fit(np.random.rand(100, 2))\n",
    "_ = UMAP(force_approximation_algorithm=True, transform_mode=\"graph\").fit_transform(\n",
    "    np.random.rand(100, 2)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "## Timed algorithms\n",
    "\n",
    "Implement parameter sweep, output logging, and timing code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def time_task(task, *args, **kwargs):\n",
    "    \"\"\"Outputs compute time in seconds.\"\"\"\n",
    "    start_time = time.perf_counter()\n",
    "    result = task(*args, **kwargs)\n",
    "    end_time = time.perf_counter()\n",
    "    return end_time - start_time, result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def run_dmst(data, n_neighbors):\n",
    "    mst = KMSTDescent(num_neighbors=n_neighbors)\n",
    "    compute_time, umap = time_task(lambda: mst.fit(data).umap(transform_mode=\"graph\"))\n",
    "    return compute_time, umap.graph_\n",
    "\n",
    "\n",
    "def run_kmst(data, n_neighbors):\n",
    "    mst = KMST(num_neighbors=n_neighbors)\n",
    "    compute_time, umap = time_task(lambda: mst.fit(data).umap(transform_mode=\"graph\"))\n",
    "    return compute_time, umap.graph_\n",
    "\n",
    "\n",
    "def run_umap(data, n_neighbors):\n",
    "    umap = UMAP(n_neighbors=n_neighbors, transform_mode=\"graph\")\n",
    "    compute_time, umap = time_task(lambda: umap.fit(data))\n",
    "    return compute_time, umap.graph_\n",
    "\n",
    "\n",
    "mains = {\"dmst\": run_dmst, \"kmst\": run_kmst, \"umap\": run_umap}\n",
    "\n",
    "\n",
    "def run(data, algorithm, n_neighbors):\n",
    "    return mains[algorithm](data, n_neighbors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def compute_and_evaluate_setting(\n",
    "    data,\n",
    "    algorithm=\"dmst\",\n",
    "    repeat=0,\n",
    "    frac=1.0,\n",
    "    n_neighbors=5,\n",
    "):\n",
    "    compute_time, graph = run(data, algorithm, n_neighbors)\n",
    "    save_npz(\n",
    "        f\"./data/generated/mnist/graph_{algorithm}_{n_neighbors}_{frac}_{repeat}.npz\",\n",
    "        graph.tocoo(),\n",
    "    )\n",
    "    return (\n",
    "        algorithm,\n",
    "        frac,\n",
    "        n_neighbors,\n",
    "        repeat,\n",
    "        graph.nnz,\n",
    "        compute_time,\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def init_file(path):\n",
    "    handle = open(path, \"w\", buffering=1)\n",
    "    handle.write(\n",
    "        \"algorithm,sample_fraction,n_neighbors,repeat,num_edges,compute_time\\n\"\n",
    "    )\n",
    "    return handle\n",
    "\n",
    "\n",
    "def write_line(handle, *args):\n",
    "    handle.write(\",\".join([str(v) for v in args]) + \"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "repeats = 5\n",
    "algorithms = ['dmst', 'umap', 'kmst']\n",
    "fraction = np.exp(np.linspace(np.log(0.1), np.log(1), 5)).round(2)\n",
    "n_neighbors = [2, 3, 6]\n",
    "\n",
    "total = len(algorithms) * len(fraction) * len(n_neighbors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 7000., 12600., 22400., 39200., 70000.])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fraction * 70000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(70000, 784)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df, target = fetch_openml(\"mnist_784\", version=1, return_X_y=True)\n",
    "df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 5/5 [5:04:38<00:00, 3655.75s/it]  \n"
     ]
    }
   ],
   "source": [
    "output = init_file(\"./data/generated/mnist/metrics.csv\")\n",
    "for repeat in tqdm(range(repeats)):\n",
    "    pbar = tqdm(desc=\"Compute\", total=total, leave=False)\n",
    "    for frac in fraction:\n",
    "        sample_idx = sample_without_replacement(df.shape[0], int(df.shape[0] * frac))\n",
    "        np.save(f\"./data/generated/mnist/sampled_indices_{frac}_{repeat}.npy\", sample_idx)\n",
    "        X = df.iloc[sample_idx, :]\n",
    "        for algorithm in algorithms:\n",
    "            for k in n_neighbors:\n",
    "                result = compute_and_evaluate_setting(X, algorithm, repeat, frac, k)\n",
    "                write_line(output, *result)\n",
    "                pbar.update()\n",
    "    pbar.close()\n",
    "output.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results\n",
    "\n",
    "NNDescent-based $k$-MST is more expensive that NNDescent based $k$-NN. Scaling appears a bit steeper but still usable. Definately a lot quicker than KDTree-based MSTs!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "import matplotlib.lines as ml\n",
    "from lib.plotting import *\n",
    "\n",
    "configure_matplotlib()\n",
    "\n",
    "import warnings\n",
    "warnings.simplefilter(action=\"ignore\", category=FutureWarning)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "values = pd.read_csv(\"./data/generated/mnist/metrics.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEKCAYAAAAMzhLIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAABcSAAAXEgFnn9JSAABHPklEQVR4nO3deXhb13kn/i8AAtwkEgC1cieuFsqyLQkktThxm0ag3CRjJ7EBMrblbaYiknaeiWeaEGJnJnaS35QGk0z1e/rMkwLu1LHkJRLoJJ06rWNAatJMQkkkIdmxbEkWwE27ROCSFHcCd/64vJcEsZAgsZLv53n02Lw8wD2gKLw45z3nvBKO4zgQQgghYUiT3QFCCCGpjQIFIYSQiChQEEIIiYgCBSGEkIgoUBBCCImIAgUhhJCIKFAQQgiJiAIFIYSQiChQEEIIiYgCBSGEkIgoUBBCCImIAgUhhJCIKFAQQgiJiAIFIYSQiDKS3YFEcTqd6OjoAAC4XC4AgNlsTmaXCCEkLUhWSj2K2tpaGI1G6PV6AIBEIoHdbodOp0tyzwghJLWtmBGFxWKBRqMBALAsCwBQq9VJ7BEhhKSHFTOiELAsC5PJBIZh0NjYmOzuEEJIykuLEQXLsnA4HDCZTLDb7eLIYDar1YrOzk4wDAOXywWDwRA0rcSyLKxWKzweDxiGiXu/OY7DCovDZB4SiQQSiSTZ3SAkKikfKJxOJxwOBzQaDdxutzhtNJvVaoXdbofNZhOvMQwDm80GrVYrXlMqleIogmEY9Pf3xyWhzXEcbt68icHBQfj9/pg/P0lfEokECoUChYWFyMrKSnZ3CFmQlF8eq9Vq0djYGDHpbDabUV9fH3BNr9ejubkZAD+SUKlUAUFGq9XC6XTGpc9DQ0NgWZaCBAnCcRzGx8fR29uLsbGxZHeHpJHhsUlMTPmScu+UH1HMh2VZuN1uKJXKgOsMw6ClpQUAP5LQaDQBbZxOZ9yWx96+fRsAnyxfu3ZtXO5B0tPU1BSuXbuGsbExXL9+PeQ0KiFzfevEebQ6rwEA/vbJHXh0R3FC75/2gULYGzF3BdPcr81mM0wmEwoKCuByuWAymcSlsvPZvn172O9duHAh4GuO4zA5OQkAKCgogFSa8oM2kkAKhQJFRUVwuVyYmJgAx3GUsyDzeu/CTShkEmRmSGH9TRcFimgJ00lzRxSzv69UKqHT6RKyZ2J28pqCBAklI4P/ZycsdqBAQSL58CqLEnUOLt4YglTK4TuP3pfwPqR9oAgXIDweT8zuMXvUwLJsQK5DCESEEBIr41M+/PLDGzja1oPzfax4PUMqxWu/60ZNRUFC+5P2H3nnbqITzDfSWKwjR46goqJC/HPkyJGYPn8imUymoCR/OK2trWAYBlVVVSG/73a7IZFIYDAY4Ha7xccYDAYYjUbU1tbCarUCmNklbzQaoVKpUFVVBZPJBIPBgNbW1qDnFpY1t7a2QqVShWwzl9VqhUqlEvsSb62trQm7F1m+rrGjaHnvIh5qPoX/cuIDnO9jkSGVYHWmDBvzs7BmlQK3BkYS3q+0H1FoNBpx6ezspbAul2vBOYhovPDCCwFTWCUlJTG/R6IYjUa0trYuKJjq9XrY7XZYrVY4nc6AnzUwc26W0WiERqOBw+EIWrJsNBoB8IFCWKbc0dEBo9GIhoYGcb/MXIcOHRKfx2KxLOi1NTQ0wGQyRWzT0tICrVYbkylJvV4Pg8GAV199lUaYJCocx+F3V/rxels3Tn5yC/7p2et1qzPxuLYIT+0pxdePduD2vQn035tAXvbqhPcxbUYUkaaSTCZT0BuIw+FAU1NTzPvx2muv4eGHHxb/vPbaazG/R6I4HI6o3iSFXI+w7FgQavRmsViCRh9CMAm3I16pVAYF95aWlqClz7Fy/PjxmD6f0WicNzgRIhgcm8Rrv+vC/v/5Gxz832dg/5gPElWlSvxQ/yB+/a3P4fAXtqFUnYtn9pbBM8wvfvj/63cmvK8pP6Jwu92wWCzingeTyQStVov6+nrxU21DQwMA/h+qsDPbYrEEfeqNhRdffBHPP/+8+HU0nx5HJ3yY8MV3b4VCJkW2Qragtna7XXwTtlqtcLlcAT/XUMxmM6qqquB2u8VpP6vVCpPJhBMnTojtNBoNzGYz1Go1tFpt0PLkhbJYLOjs7Az5PYfDAYvFArPZLI4qW1tb4XK50NHREXFKzel0wu12w2azwel0QqlUwmKxwGg0wmKxoLq6WpxGY1lW/J0S7hvquk6nQ21t7YJHPWRlunRzCEfbuvHzc9cwMsHvi8hRyPDFBzbimX2leLBIGbTA4cm95Xhyb3kSestL+UAhvOHMRwgWqeq7/3QBr/++WxxWxotUAjz3UDleejT8kl6Bw+GA2WxGS0sLGhsbYTQag6bw5tJqtdBqtTCbzeIbosvlCtoPIPydNTc3i2/ETU1NUZ+vFWqPDDCTExCmpFiWhcFgEIMKy7KoqKgI+7wajQZqtRq1tbXQ6XRioLDZbLDZbGhtbRW/BvjpstbWVmg0mpDXhZGQELBofwSZbdLnx/sXbuH1tm6c7ZqZHSkvyEF9TQnqqktQsCoziT2MLOUDRao5cuQIvvvd74pfv/TSS3j55ZfnfdzRtp64BwkA8HP8veYLFMInYuETNBA4/280GsURQl1dXcD3zGYzamtrYTabceLECTH3MNfsAG+1WsVj3hf6Jup0OkO2bW5uhkajCch/nDhxAtXV1eLXc4NLqNczdyOmWq0Wcyz9/f1gWRYtLS1QKpVQKpVwu91ob28PeV2g0WjC9pusPLcGx/DWmV68fbYXt4fGAQAyCfBHW9bi6T2l+NzWdciQpX4GgAJFlBabzH52X1lCRhQyqQTP7iubt53D4YBer4fZbIbBYIDZbA54XRaLJewUik6ng1arRXNzM1iWDTmam/upWkguR3O0u1KpDJmb0ul0sFqt4kgImClGFU6k1xOKMLISnl94jQaDIeR1gcfjoWT2CsdxHM52eXD0dA9+9dFNTE3/o1fnyvGVncV4Zm8pKtauSnIvo0OBIkqvvfbaokYULz26HY2PVKZMjsJut6O2thYajQY6nQ5utxtOpxMejydsgnv2nL8wqrDb7SHbWiwWNDU1iW+abrdbnOJZKI1GEzLPUFNTIy65FQ52rK2tFUc5SqUSDodjQct+wy2SEJ5r9muIdF1A004r1/D4FH5x/hqOtfXg4s0h8fqDRfn42u4SfHlnEXIz0/MtNz17nURLSWZnK2TIxsISzfHGsizq6uoAAPX19bBYLGBZNmwOobW1FQ6HQ/wUr9PpoNfrxaDS0tIClmVhsVigVqtRU1ODQ4cOQa1Wi0e6z54qAvjpKCHPUF1dHTI3otPpxOW4DocDHR0dUCqV0Gq16OzsxP79+8U+CYl2jUYDk8kk5hPCvakLezmMRiO0Wq343BqNRlyuu3//flRXV4s5lnDXlUolWJaFWq2mQLHCuO7cw7G2HrzTeRVD41MAgCy5FI9s34Bn9pWhqlSV9rvvV1zhonjz+/24dOkSAGDr1q10jMcSCfsx0qG+ubChcL6FFfQ7kv58fg4nP7mFY6d78NtP74rXi5XZqKspxtdqSrEub/kcI08jiijNPcJDSGiS+NDpdLDb7Sk/pcOybNAGQ7L89N8bx/GOPrx5uhfX2FEAgATAQ0wBnt5bitptGyDPWH6BnwJFlBa76oksntlsFotXpaqOjg4KEsvY+T4WR9u68e6HNzAxxecZ87Mz8OWdRXhmbyk2r89Lcg/ji6aeojTfiIKmFch86HckPYxN+vBPH1zHsdM9+PDqgHi9csNqPLm7FI9ri7A6S57EHiYOjSiiRFNNhCxvfZ4RvHGmByfa++Ad4WvLyGUS6Latx7P7yrBXU5D2yeloUaCIEuUoCFl+/H4O//bpHRxr68GpS7chzLOsz8uEoaoET+8pxUZldnI7mUQUKKJEOQpClo+BkUnYOvvwxukedPfPHN+9u1yFg3vL8Kf3b4RiGSano0WBIkpL2UdBCEkNF64P4FhbD35x/hrGJvnkdK5Chkd3FOLZfWW4rzA/yT1MLRQookRTTYSkp4kpP/7loxs41taDjh6veJ1Zm4uv1ZTCUF0MZY4iiT1MXTSmWsHSrcLdXImuYrdQVO0utdwYGMWP3r+Eh145hW/+9Dw6eryQSSXQbVuHn7xQDft//iMc+iMNBYlIOBIVr9fLdXV1iX+8Xm/A930+H/fxxx9zH3/8Mefz+ZLTyQVyuVycRqNZcPuGhgYOANfZ2Rn2e3a7neM4jrPb7VxDQ0NQG47jOLPZLF7TarWcxWLhOI7/2dpstqDn1uv1YfukVCo5l8u14NcQL2azWXztHMf3ee7vhiCdfkfSld/v53535Q739WMdnKbpl1yZ6V2uzPQup/3e+9xf//JjrufuvWR3Ma3QiCJKy6lm9kqvcBdLc6vlUbW72BgenxI3uC3EvfEpHGvrxoG/+Tc89eoZ/MtHN+Hzc9hVosQP9A/id4c/j6YvbkNpQW4ce738UI4iSktKZk+MAL6JmPcpgEwBKHIW1DQdK9zNV8UuVPU5lmXR3NyMgoICtLe3w2g0iifmWiwWsb1wEm6o52hpaRELFlksFng8HvF55lbLEw5NpGp3S/N7112Y37sEdY4cf1O/M+LU0JXbQ/zBfM5ruDd9MF+2XIYvPrABz+4rw44SVaK6vSxRoIjSopPZ/3IYOGsBuPgeMw6JFNhtBL7wyrxN063C3XxV7JxOZ8jqc263GwUFBWhsbATLsnC73XC73TAajWJwaGlpifgcjY2NMJlM4s/sxIkTqK2tBcdxQdXyBFTtbmnO9bLw+zncvTeBrrvD2FUaGCimfH44PrmFo209+L2rX7xequarxj25uwTq3NStGpdOKFAkyllr/IMEwN/jrHXeQJGOFe7mq2J3/PjxkNXntFqtWDvDYDCIRZRqa2vFxwoBLNxzCPfT6/VQKpVoaGgQA6swWpo7aqJqd0vz2I5CdN0ZxkZlFh4omlmuemdoHMfbe/HmmV7cGBgDwJcAfnjzWhzcW4r9levoWJQYo0CRKLsbEjSikPH3mkc6Vribr4pduKp0AOD1enHixAlYLBbY7faw/Yj0HAACHjffyJKq3S3Nq7914x3nVcikEjz24EYMjk/haFsP/vkPNzDp47dOK3Pk+OrOIjy7ryztqsalEwoUifKFV4D930mZHEU6Vribr4pduOpzwtRaQ0MDdDodjEYj6uvrcejQoagr2EV6859bLY+mnZbmF85r4MCPFr726hn0D8/827m/MA9P7SnFV3YVIUdBb2PxRj/hRFLkAFhYojne0rHCnU6ni1jFLlz1OYCvdV1TU4P+/n7YbDZxSqmqqgparRZqtRoWiyXsc5w4cUJ8fbOn1ISvZ1fLE3IhVO1u8brvDiMnU4Z7E1OY8HHoH56AIkOKR+5bj+ceKkN1eUGyu7ii0DHjUaJjxhMrnSrczRap2h39joTm83P4zeXbONrWg99cviMezJchlaCqNB9/+3QV1q1ePlXj0gn9hkZpOe2jSAfCiCWddjoLS23nK4lKeN7hCVh+48Lnfviv+Pc/6cCvL/FBIkchw/q8TBSrsvFAcT4FiSSiqaco0aGAiZcOFe5mo2p3C/OHqwM42taN//PBdYxPb6pblZmBL+8sxHP7yvDNt524MTSBu/cm4BmZSnJvVzYKFFGiQwGTI5od5MmWTn1NtPEpH/75DzdwtK0H53pZ8frmdavw1J4SGKpKsSqLf1u6v1iFTzqvQiYBXnioIswzkkSgQEEIibtr7CjePN2D4+194uqlDKkE+7etx3P7SrGPWRNUNe4Hhh34r1/ahhyFDIoMWTK6TaatmEAhHNcA8Bu5amtro94lTAhZOI7j8Lsr/Tja1g3HJ7fgn05Or12VCUN1MZ7ZV4aN+ZGrxtGJrqlhxQSK2cc1AADDMOIOW0JI7AyOTeJnnVdx9HQP3HeGxevVZSo8vbcU/+6BjZDTCCGtxDRQnDp1SjwaQjjSoLq6Gnl5ebG8zaIIB8gJ+QWtVksrUwiJoUs3h3C0rRs/P3cNIxM+APzBfI/u2Ijn9pVjexFVjUtXSw4Ug4ODaGxsxKuvvgphS4ZSqRT3GkgkEuh0Ohw+fBh/8id/stTbLZrX6w342ul0hj2jiBCyMJM+P96/cAtH27pxpmtmZ3p5QQ6e2lOK+ppS5GfLk9hDEgtLChTvvPMODh06hLq6Opw4cQIajQa7du0KaON0OnHy5EkcOnQIDMPAZrNFPcJgWRYOhwMmkwl2uz3kMkmr1YrOzk4wDAOXywWDwRB29YlQRY1yFIQszu3BMbx9tg9vne3BrcFxAIBMIsEfb12DZ/aW4XNb1wUlp0kaW2zFI6vVytXV1XEsyy74MXa7nTtw4AA3MDCw4Md0dnZyZrOZs9lsYaurWSyWoCpoGo0mZFuXy8XpdLqw1ceWKp2qlzU2NnJKpXJBPwubzcZpNBpOp9OF/Z7w97SQditZOv2OzOb3+7kz7n7uL97s5JhZVeN2fe997v979wJVjVvGFh0oHA7Hoh7Hsix37ty5qB/n9XrDBgqNRhNUQrOxsTEoeLhcroBroZ5rqdLpTWAxpVCVSiXX2NgY8nvRtlup0ul3hOM4bnh8knvzdA/3yN/8RgwOZaZ3uUf/9rfc8fYebmxiKtldJHG26Kmn/fv3R/z++fPn0dHRgbq6uoCppvz8fOzcuXOxtw0iFKKZuwmOYRixGA3AL481mUwwm81gWRYejwcOhyNikZ7lbjGlUE+ePImqqqqgIj2LaUdSm/vOPRw73YPWzqsYGuN3RmdmSPGF+zfg+YfKsbOUqsatFDFb9XT48GF0dXXBbDbD6/Wiuroa+fn5MJvNsNvtKC8vj9WtAnR0dABAUH2BuV/X1taKZTQFCz1mYfv27WG/d+HChYV2FaNTo5j0Ty64/WLIpXJkZ0Remy5YTClUrVYLi8WC2trakJXtom1HUovPz+HUxds42taN3356V7xerMrGkzUleHJPGdS5tLdhpYlZoHA6nXj//fcBAK+88gqUSiX6+/vBsiyamprw4x//OFa3CiCsrgp3rIawJHa+ojfxZj5rxlsX34I/zoWLpBIpnqp8CqbdpnnbLqYUKsCfiCrUs4j0c11oO5J8/ffGcbyjD2+e7sU1dhQAIAHwmU1r8Oy+Mui2rYdUSsnplSpmgWL2J0aHwyHWOoj3uUgLLSKzFLNHDXOPGZ+9NyOSty++HfcgAQB+zo+3L749b6BYSilUgB+NMQwDg8EQcWS20HYkOc73sTja1o13P7yBiemD+fKyMvCVXUV4bl85mHVUNY7EMFAIb8xdXV1wu90B9YjjeUS0EKBmv3nP/jrWgerIkSP47ne/K3790ksv4eWXX573cU9WPpmQEYVMIsOTlU/O224ppVAFdrsdDMOItReW2o4kxtikD+9+eAPH2rrxwdUB8XrlhtV4ek8pntAWIydzxRzaQBYgZr8NdXV1UKvVGBgYAMMwePzxxzEwMACHw4H+/v5Y3SaIRqOBRqMJmjJxuVzQ6/Uxv98LL7wQ8IZaUlKyoMeZdpvwn7T/KWVyFEsthQrwP3ubzQaDwRDws19oO5JYfZ4RvHGmByfa++Ad4X8P5TIJDty3Hs/tK0dNhZr2PpDQYrmEyuVyca2treLeCqfTybW2ti56Ke3c50aEfRRz1+2H20exVC+99BIHQPzz0ksvBXw/XZY+zt5L0tnZyTU0NETc3xBpD4SwHyOaditZIn9HfD4/9+tLt7n/8JOzXPnhmaWtu/+HnfvRry5yNwdG4np/sjykfClU4dRXp9MpLufUarVBq3Oi2Zm9FD09Pejr6xO/LikpQVlZmfj1Si1z6Xa7F7SyaaHtlrNE/I4MjEzC1tmHN073oLt/RLy+u0KNZ/aW4Qv3b0CGbGX8bpKlW3Sg+PrXv45vfOMb2LFjR1SP++EPfwitVovPf/7zi7lt0r388ssRcxQrNVCQhVvo78iNgVGMTPjArF14Qvnj64M4drobvzh3HaOT/MF8uQoZHttZhOf2laFyY/IP6CTpZ9GBgmVZ6HQ61NbW4q/+6q+wevXqiO1/9rOfobm5GUajEX/2Z3+2qM6mgrmrnuZWvKNAQeazkN+RizcG8cw/nMWUz4//+qVt0FeFz4VNTPnx3oWbOPr7bnT0zBx+qVmbi6d3l6KuugSr6WA+sgSLTmYrlUp0dHTAZDIhPz8fDMNAq9VCo9GAYRiwLIv+/n44nU50dHRAo9HAarUGHRpICAn23kc3cWeIP2zvRHtfyEBxc2AMb53pwVtn+3D33vTBfFIJPl+5Ds/uLcNnNwdXjSNkMZa86slsNqOpqQlWqxUnTpyA3W4PqkfxyiuvLJsAsdjlsYREw3VnSPx/YQMcwFeNO+324Njpbvzqwi34psvGrVmlwBPaYjy7rwxFqpyE95csbymfzE41NPVElmohvyOHWz/ETzv4RRMlqiz8y4t/jJ87r+LY6R5cvnVPbLerRImn95bisR2FVFeaxA3tqonS3MAw1+yhvt/vp0BBgkxN8QfsSSSSsFNDXXdnRhR3701g71+fxL1x/nHZchm++MAGPLevHA+WKOPeX0IoUERpvhGFRCKBXC7H5OQk+vv7sXbt2sR3kqSsqakpXLt2DQCgUChCBoopnx/9wxOQSyWY9HMYnfQD8KNUnYMnd5fgyd2lUObQwXwkcWjqKUrzLY8F+PKwwpsBIaHIZDKUlpYiKytLvHZnaBzH23vx1pleXB8YE68X5MphfuJBfL6SDuYjyUGBIkrzjSgAPuF48+ZNDA4Owu+P/0GAJH1IJBIoFAoUFhYiKysLHMfB2cviWFs3fvmHG5j08f8c87IysDozA8ocOb73lftRVaae55kJiZ+YB4pTp06BZVk8/vjjAPhP19HWyF4uOL6CYLK7QVKIkJcYnfDh/3xwDUfbenDh+qD4/e2FeXh6byke31mELAXNDJPUELNAce7cOezfvx8sy4JhGHz66acAgHfeeQcSiUQMHOluISMKQsLp6R/GG6d7cKLjKgZG+YP5FDIpHtm+Hs89VI7qcho5kNQTsyU5JpMJTU1NcLlceOKJJ8TrTzzxBI4fPx6r2yTdkSNHUFFRIf45cuRIsrtEUpzfz+FfL97GC6+dxed++Gu8+tsuDIxOYmN+Fv6zbjN+3/R5/O1TWgoSJGXFtHDRt7/9bQAIWsnhdDpjdZuke/HFF/H888+LX9NogoTDjkzA1sHvfej1zBzMt09TgGf2leKR7Rsho+Q0SQNxmQT1emfOmzl37lw8bpE0NNVE5vPRtQEcbevGP56/jvHpqnGrMzPw2M5CPLevHFs2RD4XjZBUE7NAUVVVhaamJjQ1NYkjipMnT6Kurk4si0rIcjU+5cM//+EGjrb14FwvK17fvG4Vnt5TCkNVCXKzKDlN0lNMVz3V1tbi1KlTAdcqKirgdDqXzconSmaT2a6xo3jrTA9+erYP/cMTAIAMqQT7t63Dc3vLsW9TAR3MR9JezJfHOp1OnDx5Ev39/aipqQlIbC8HC9lwR5Y3juPwuyv9ONrWDccntzB9Lh/WrspEXXUxDu4tw0bl/OVoCUkXCdlwd/78eezcuTPet0kIGlGsXENjk3ink09Ou+4Mi9erylQ4uLcU/+6BjZDTwXxkGUpIoPjGN76BH//4x/G+DSFxcfnWEI62deNnzmsYmeCrxuUoZPjSAxvx3L4y3F+sTG4HCYmzmAWK7u5uGAyGsEthfT5fLG5DSEJM+vywf3wLR9u6cdrtEa+XF+Tgqd2lqKspoYP5yIoRs2UYer0ebrcbhw4dQlVVFdRq2jxE0s/twTG8fbYPb53twa3B6apxEgn+aMsaPLuvDJ/buo6S02TFidmIQq1Ww+l0ory8PBZPR0jCcByH9m4vjrZ1472PbmJqOjutzlXg8V1FeHZfGUoLcpPcS0KSJ2YjCp1OF/Z7yymZTZaPkYkp/OLcdRxt68bFmzOFgh4oyscz+0rx2INFyFJQcpqQmI0oBgYGcPjwYbS0tGD16sCdp8spmU2rntKf+849vHG6F7bOPgyN8VXjMjOk+NP7N+D5feXYVaZKcg8JSS0xCxRdXV2ora1FV1dXyO8vl2Q27aNITz4/h1MXb+NoWzd+++ld8XqRMhtfqynBU3tKUbAqM4k9JCR1xSxQHDhwAB0dHWhoaADDMGIyu7+/H1arFR0dHbG4TdLRiCK9eIYncLy9D2+c7sE1dhQAIAHwmU1rcHBfKWor10Mmo7rmhEQSsxxFR0dH2GT2clolQoEhPZzvY3G0rRvvfngDE9MH8+VlZeAr08npTevoYD5CFiohyeza2tpY3YaQsMYmfXj3wxs41taND64OiNcrN6zGU7tLoa8qRk4mHcxHSLQSksyur69PieJFbrcbZrMZVVVVaGhoSHZ3SIz0eUbw5pleHG/vhXeErxonl0mg27Yez+0rwx4NHcxHyFLE7ONVRUUFBgYGYLVaY/WUMeV0OtHR0QGPxxOQYyDpye/n8Nsrd3GsrRsnL96G8HFnfV4mDFUlOLi3FBvy6WA+QmIhZoFi//79kEgkqK+vD7jOcRxeeeWVWN1m0bRaLbRaLTo7O5PdFRKBd3gCZ7r6sbNEhQ35WUHfHxidRGvnVbxxugddd2cO5qspV+Hg3jJ88f4NdDAfITEWs0Dx9a9/HWq1Grt27Qr63uyKd4REUm9pg+vuMNasUuDMX83kvT6+Pohjp7vxi3PXMTrJL7XOVcjw6I5CPLevDNsK85PVZUKWvZiOKMKpqalZ0nOzLAuHwwGTyQS73Q6NRhPUxmq1orOzEwzDwOVywWAwREywk9T06e174ADcGhzH2MQU3v/kNo61daO9e+bDhmZNLp7cXYL6mlLkZcuT11lCVoglBYrBwUGxct3g4GDYds3NzYtOZjudTjgcDmg0Grjd7pD5BavVCrvdDpvNJl5jGAY2mw1arXZR9yXJIaysUMgkePgHv8adoemD+aQS/MnWtXhmXxke3rQGUintfSAkURYdKA4cOIDOzk709/cDAMrLyzEwMDDPo6In5BYiJaDNZjPMZnPANb1ej+bm5oDgQVIXx3E47fZAIZNgwsdhwsfhztA4CnIVeKKqGM/sLUOJOifZ3SRkRVp0oDCZTAFv3nV1dVAqlUF7JjiOw+HDhxfdwfmwLAu32x20CY5hGLS0tMTkHtu3bw/7vQsXLsTkHivVvfEp/PzcNRxr68blW/fE61kZUnz/y9vx2M4iZMopOU1IMi06UMzNSTz55JOoqKhAaWlpUNumpqbF3mZewtEgc+tfzP1amMJyOBziY/R6fch8B4m/K7eHcKytB+84r+HeOH8wX5ZcCrlMirwsOTJkEuirS2j/AyEpIGbJ7MuXL6O9vR3f+ta3Aq6fO3cODocDTzzxRKxuFUAY1YQ7VoNlWSiVSnEKq7GxMep7zB41zD3rSXh+Mr8pnx+OT27j2Olu/O5Kv3i9RJ2Nr9WU4ms1JTjiuIx//sMNMMpVFCQIEXS+Bbz7HwGpDGjsAjJXJfT2S05mA0BeXh4KCgpC5ig8Hg+sVmvcjhkP9ybt8XhCXl+qI0eO0OmxUbp7bxzH2/vw5ukeXB8YAwBIJcBnN63BwX1l2L91nXgw3/e+fD++qdsCNZUZJSsd2wf0neH/dP4E4HyAVAH8z/uApt6EdmVJgaK9vR21tbUBn/xMJlNQu3iuPBKmjuYmu+cbaSzWiy++iOeff178mkYToXEcB2cvi2Nt3fjlH25g0sevZ8rPluOr0wfzadYGfyqSSCRYQ8d9k5XGNwnc/APQdxboO83/d/BacDuJDJAn/sSBJQWK/fv3w2KxQK1W4+zZs+jq6gramQ1EPjBwqTQajbh0dnZAcrlc0Ov1cbsvCW10wod/+uA6Xm/rxoXrM0um79uYh6f2lODxXXQwHyEY9QJXO4De0/yI4VonMDkS2EYiA9ZWAoU7AfUWwPkPgCwL+IvTCe9uzA4FBICTJ09G3Hi3FG63GwzDoLOzM2iEYrVaYbPZYLfbxWvx2kdBhYtC6+kfxhune3Ci4yoGRvmD+RQyKWq3r8ez+8qwu1xNOQeyMnEc4HHPBIW+M8Cdi8HtMlcDG3YAxVVA8V6guBrIUQGy5G8qjWmgiAe32w2LxSKuWtLpdNBqtaivrw8IAonamd3T04O+vj7x65KSEpSVlcX8PunA7+fwm8t3cLStG7++fEc8mG9jfhYMVSV4ak8JHcxHVp7JMeDGB/wUUu90YBi5G9wuv5QfLRTXACV7gfXb+WmlFNxMmvKBItXQiAJgRyZg67iKN870oKd/Zri8V6PGwT1leGT7BsgzUu+XnZC4uHd7ZqTQewa4cR7wTQS2kSmAddv5wFCyByjdC+QVARnpsWiDAkWUVnIp1I+uDeBoWzf+8fx1jE9XjVuVmYHHdhTimX2l2LaRDuYjy5zfz08bCYGh7ww/rTRXtpoPCkVVfGAo0gKZ+YAsPfNz6dlrkjDjUz78yx9u4vW2bpzrZcXrm9atwpM1JdBXFyM/Oz0+FREStYlhPuksBoZ2YHzuNgAJUMAAG3cBJdVAyT5gzRZAngUsk7xcQgLF+fPnsXPnzkTcKu5Wyj6K6+wo3jzTg5+e7UP/MD+Mzph1MN9n6WA+shwNXJ1OOk8vU735Eb9/YTZ5NrD+AX6UULwbKN0D5K5Lm2mkxYj51FN3d3fA1yzLwmQy4Ve/+lUsb5M0y3nqieM4/N7Vj6Nt3bB/fAv+6d+Mtasy8URVEQ7uKUMxHcxHlgvfFHDro8D8wuDV4HarNgCFu2amkQp3AIpV/C7pFSJmI4pz585Bp9MFbXzjOI6WRaa4obFJ/Mx5DUfbuuG6M1M1TluqxFO7S/GlBzciW0GzlCTNjbLT00in+VFD2L0LW6dXI+3mk84qDT9aWMHvYzH713/o0CGoVCpYrdZlfdDecpp6unxrCEfbuvFz5zUMT/DD62y5DF98YAOefagcDxblU5An6UnYuzB7p/PtTzBT8WSaYjWw8UF+tFC8GyjZDeSoU2LvQiqJ2dSTVCqF2+1GeXl5LJ4uZaX71NOkzw/7x7dwtK0bp90z52GVqXPwtd2lqK8uhpqO0CCpwv1vQN5GYM3myO2mxoHr5wNXIw3fCW6XX8JPIxVXT+9duD9l9y6kkpiNKLRabdDR3oLZlfDSXboFBsHtoTG8faYPb53twa3B6apxEgke3rwGB/eW4k9mHcxHSEp47YtAz+/4/3/aBmw+MPO9e3eAq2dndjtfPxe8d0EqB9bdN5N0LtsL5BUv66RzvMRsRNHV1QWr1Yrm5uag733jG9+I2+mxiZZOIwqO49DR48XRth6899HMwXyqHP5gvoN7S6FZuzrJvSQkjO+tAfz8cTAo3g3sfGpmKink3gUVsHEnUFTNTyEVV6X13oVUErNAUV1dja6uLgDBRYPcbjd8Pl+oh6WddNiZPTIxhX88fx2v/74bF28OidfvL8zDk7tL8eVdhViVSXOwJEVNDPOJ5nf/C8D2AP4pgPMHt1MzgdNIayuX1d6FVBKzQFFXVweWZYMO4eM4Dq+++mrc6kMkWiqPKLruCgfz9WFojK8al5khxYHt6/HM3jLU0MF8JBUNXJtJOPee5o/bnrt3QZbJJ50LtfwS1ZI9wKrlvXchlcRsTFZbWwudToeKioqg761ZsyZWt0m6VAoMAODzc/jXi7dx9HQP/u3yTPKuSJmNuupiPLm7FOvyspLYQ0Jm8U0Bty9MH5Y3HRwG+oLb5a4DJoeBjCw+2bxmG/D0CUo6J0lMl8eGwzBMrG5DpnmGJ3Ciow9vnO7BVe8oAEACYB9TgKf3lKJ223oo5CtnQxBJUeLehTMzdRcm7gW2kUiBNVuBol3Texf2AWoN8OPPAlPDwNQYkFNAQSKJYjb19LOf/Szs95qbm9He3h6L26x4H/SxONrWg3/68Dompg/my8vKwGM7C/H0nlJUbsij6SWSHBwHeLundzlPr0YKuXdhFbBxx8xO5+LdQG5B8E7n6+eA17/MTzH9xVkKFEkUs0CxadMmdHV1Ye7TCW9ayyWZnYwcxdikD7/88AaOtnXjg6szB5JtXb8aX6spwePaIuRTjWmSaFPjwI0PZ3Y6950Fhm8Ht5u7d2HDg5R0TjMxm3rS6/WoqamBTqdDfj5/3LTT6cThw4dhtVpjdZukS+TO7D7PCN4804sTHX3wTB/MJ5dJ8PnKdTi4txT7NGuQQXsfSKIM3w08F+n6OcA3HthGJgfWbpve6bwHKNvDBwra6ZzWYrqPoqCgIGhjndPpRFNTEx0KGMbw+BTO9bLYXpgHVa4Cfj+H/3vlLo629eDUxZmD+datzsQTVcV4ancJStS5S3sRhMzH7wfuXg5cjeRxBbfLUvGH5BXV8KeoFlUDWfk0TbTMJKRwkUwmWzZTT7Gm+9Gv4b47jDWrFDD+8Sa8cboHXXdnDuarLlPhyd2l+ML9G5CTSRuHSJxMjPCJZvEIjLPAGBvcTqWZlXTey48e5HTky3IX13eewcFBWK3WlFpOmmqu3BlGVoYUd4Ym8P13PwYA5Cpk+OIDG3FwbxkeKMqHVEpzuSTGBm/Mquk8vXfBPxXYJiOLPwtJSDqX7gZWbaSdzitQzP7G1Wo1BgbmVn7iN9wZjcZY3WbZkQAYm169VLEmF/U1JTBoi1Gwmj6lkRgR9i4IU0h9Z4GB3uB2uetmJZ338BXbMlfRNBKJXaDQ6XRQq9Wora0NuK7RaLBr165Y3WbZ4cCfvZQpl+GX//Eh5GTR6iWyRGMDwNX2mcAQdu/CFv7AvKIaoGwfP61E00gkhJgFCqPRiOrqanHF03IVj+Wxfg7w+TjAT1NMZNrYIPA39wOTY8Cz/wiU7wvdbvbeBWE10u2PEXLvwoYH+YPyiqePwAi1d4GkrKtDV5GTkQN1duhTuuMpZoFi//79Yb9XX1+P48ePx+pWSRWP5bEDo/wJmZmZ9I+WTPuHPwXGp6dyjz0G/Pfp41mmJoAbHwTWXbh3K/jxecWBVdrWPwAosmnvQpqyfGDB/zr/vyCRSPCTR36CXesTO0uzpFVPp06dgtPpBMMw+OpXvxqyTV1dHd55551ls+op1iOK4+29ePW3bjz6YCG+qduy9A6S5eFvHpjJI8iygL1fnz4Cwxm8d0GaMavuwh5+maqylPYuLCOfeeszGJwcBADsXLsTx754LKH3X/SI4tVXXw1IUh84cADvvfee+PX58+dhMBjgcrnQ0NCwtF6mkFjvxK6vKUV9TWnMno+kMb8f6P+UDwhrNgND1/mVSL4x4HdHZtplKfnRQtF00rmoiq/FQEnnZWu1fLUYKCpVlQm//6JHFJs2bYJWq0VTUxM4jsPhw4fxyCOP4C//8i/xwx/+ECaTCfn5+Xj11VfxxBNPxLrfhKS/iRF+d7OwTPXqWWDUG9xOlglUfgkoruEPzFt3HyWdV5gdr++AVCKFTCpDQWYBfmVI7AbmJQWKK1euBFw7cOAAJBIJ7HY7tFotTp48mVLJbaPRCIZh0N/fj5qaGuj1+mR3iawkgzcCcws3Pgizd2E7MDUJjNzmv161EfgP74V+TrLscByHWyO3cMlzCZe9l3HJewmnek9h0j+J7IxsFOcW42dfCX8IazwseuppboEiANi1axd+8IMfwGw249vf/rZ4PRVqZhuNRhiNRrHftbW10Gq10Gg0Se3XhG8Cn3g+AZPPYJViVVL7QmLI7wNuXZjZ5dx3GmBD7V1Yy5fvLKnh8wuFWn7vgkQC2P49cO8GUJfY+WiSOKNTo3CxLjEoCH8GJwZDtpdKpPj+Z76f4F4uYUQRaiXTO++8E3LfRFNTU8ha2okkkUgCTra1Wq3o7OyExWJJYq+AV86+gnO3z6F4VTF+9LkfJbUvZAnGBmf2LvSd5mswhNq7ULAFKNo5nXTey5fzpGmkZY/jONwcvolL3ulRwnRg6B3qhT9EmVeZRIaS1SWoyK9ARX4FNis3Y5NyE8ryypAtz054/xc9oujs7MSf//mfBx0rHsqJEyeSGiicTmfQNbVajY6OjiT0JtDtEf5Y5rujd+Hn/JBKKCGZ8jiOr+XcO+tcpNsXgus6y3P58p1FVUDJbj7xnLOGjsBY5kanRnHFe0WcNhJGCUMTQyHb52fmQ5OvQUVeBTRKDbaotoBRMsjNyIVcJodcKk96jZlF/8a63W783d/93YLaLvVFsiwLh8MBk8kEu90ecrpIGCEwDAOXywWDwQCdTgcA8Hg8IVcqzV7mmix/sfMvcLL3JPZs3ENBIlVNTfBnIfWdntnUdu9mcLu8Iv4IjKKq6b0LDwKKHFqNtExxHIcbwzfEEcIl7yV86v0UPYM94OZueETgKEGTr8Fm5WZsUW/BhpwNUMgUyJBmQJaiGyAXHSh0Oh3ef//9BbU9cODAYm8Dp9MJh8MBjUYDt9sd8s3darXCbrfDZrOJ1xiGgc1mg1arhVqtDvm4VDiskFEyYJRUKjaljHhmppD6zvJHYEyNBbaRZvAnpxbu4kcLpXsBZRmQQUewLEcjkyO4wl4RRweXPHxQGJoMPUpQZirFgMAoGWxVbcUm5SZkZ2QjQ5YBuTS99rgsOlCYTKYFt13KoYBarRZarTbip3+z2Qyz2RxwTa/Xo7m5WQwWc7ndblRXVy+6X2SZ4Djg7vTeBSEw3L0c3C4rH9iwgz8wr3g3UFLN12KgaaRlheM4XB++HpRc7h3sDTlKyJBkoCSvRJw62qLagsqCSqzPWY8MaUZKTBvFwqJ/yyMd2TFXPPdRsCwLt9sdNDpgGAYtLS3i13q9Hg6HQ5yOam9vR1NTU9z6RVLU5Ci/d0E4RbXvDDDqCW6nquBXIxVXT9dduI8v30nTSMvGyOQIPmU/DRghXPZexr3JeyHbqzJV4ihhk2oTtqq2gslnkCPPSelpo1hI+49DQkJarQ48KGvu1zabDUajEU6nE/39/eLy2IXYvn172O9duHAhyh6ThBq6OZNX6Dsdeu+CLJPfu1C4i084l+wB8gppGmmZ4DgO1+5dE5PLn3o/xSXPJfQN9c07Stik3ITNqs2oVFViXc46yGVyZEjT/m0zamn/ioUpqXD5BpZlxe8leyksiTO/D7j9yayCPGf41UlzCXsXhLORincBijyaRkpBI5MjkEvlkC/w3CphlDB36mh4cjhke2GUsEm5CZtUm1CpqgSTzyBbnr3ge64Eaf8vI1yA8HhCTCcs0uxRw9xDAWcHIhIDA9eAc28C9z0GrJvnTJvxocC6C1c7gLlLECVSoGDzzKa20r2AiqFppDRwsvckvvN/v4O8rDz89Es/RX7mzCkPUY8SpBkoXV3KrzZSbcYW1RZsVW3F+tz1yyaPEE9pHyiEpbJzk93zjTQWKx7HjJNZjj4GDFwFOl4FvvXpzHWO43c2C6uRes+E2buQw9ddEFYjlewGctby00j0ZpBWvtf2PQxODmJwchCvnHkFu9bvWtAoQZ2l5kcJ+fy00Vb1VnGUsBKnjWIh7X9qGo1GXDo7O+fgcrnicpbTCy+8ICbEAaCkpCTm91jR7t0GpsaBES9wtTPwbKShG8Ht8wr5kp2Fu4CSvfwGN0UOHbGdhvycnx8lePhAMD41DplEBh/nw7td7+LdrncD2sulcpSsLgGjZLBZuRmbVZuxTb1NzCWQ2EmbQBFpKslkMsFisQQEBofDEbCvIlZee+01GlHEg7B3QSoDZArANwH8/ecD20gzgLWVs85G2guoSoGMTKrUlmaGJ4fFVUaz8wkjUyMh2+cp8rBFtUVMLm9VbcUmFb8vgTaqxt+SChclgtvthsViETfe6XQ6aLVa1NfXB4wgIu3MjqV4lEJdcTgO6HcF7nS+eym4nbB3oXAXHxiKqvm6CxmZNI2UJvycH9eGrgUcZ3HJcwlX710N2V4ulaM0rxRMPoP2m+2Qy+RQSBVYl70Or3/x9QT3nghSPlCkGgoUizA5Flh3IezehXL+cD2pDJDnAcZTfM6BlqmmBWGUMPt47E+9n4YdJRRkFYjJ5c3KzahUV0Kj1CBHngMAeKT1EQxPDWPSN4kvlX8J3/nMdxL5csgsFCii9PLLL9PU03yGbgXmFq6fB/yTgW1kCn7vwsZd/Ka2kj1A3kbA8jBf6W1qDHjxAiCjKaVU4+f8uDp0dWaUMJ1TiDRKKMsrA5PPYLN6MyqVlags4PclRFptdNlzGd/8129iQ+4G/P2Bv1/WG9pSHQWKKNGIYg5x78J0UOg9HXrvQk4Bn1so3MXXdN64E1Cs5pepzn6z+P4awDcdVP5bP5CRNmm0ZenexL2AfQnCKGF0ajRk+zXZa8SNasIS1E2qTcjKyEpwz0ks0b9CEp3xIf6QPGEK6Wo7MD63yIqEr/m8YQe/qa1kD1CwCZBnzz+N5Js18qDRRML4OT/6hvoCDr277L2Ma/euhWwvjBLEgKDeiq2qrfOOEkh6okARpRW1j4LjgIG+mQ1tfaf5qm0h9y48MH02Ug3/Z9UaQJYV/W7nh/8SaP8HYMsXKGEdJ0MTQ3wuQaiV4LmMT9nIowQmn8EmFR8UKtX87uXMDCq4tFLQ1FOUlvXUk28SuPnhrMBwFhi6HtxudSGwcQdQuJM/AmPjA3yRnoxM2u2cQnx+H67euxo0bRRplFCeVw5GyYgBoVJdibU5axPcc5JqaESxko14po/AmF6ieq0TmPupUiKb3rsgTCPtBVRlQEYWrUZKIUMTQ0G1EhY6SqhUVYorjhQy+jslwShQRCltp57EvQuz6i7cuRjcLjOPDwrCEduF1UCucnHTSGRJOI6Dn/MHrPaJdpSgkCpQnl/Ob1SbXoJaWVCJNdlrEvUyyDJAU09RSpupp8kx4Mb5wLoLI3eD2ynLplcj7eSTzuu2ARnZNI2UZHdH7sJoN2JocgiPlD2CocmheUcJa7PXglEy2KTchK3qrdim2gaNUkPHWZAlo0CxXNy7PbM8NdLehXXb+aBQqOV3O+cV8qMFOjQvaXx+H3qHegOmjc7dPoeBiYGQ7TNlmWIuQaiVsE29DQU5BQnuOVkpaC4hSjEfUfgmgTuXALWGP8xuIfx+ftpo9k5nb1dwu2w1v29h406gqIr//6w8frRAnzKTYnBiUNygJgSGK+wVjPnGQraXSWTIkefgK5u+gq2qrdiq3gpNngYKyg+RBKJAEaWY5yhOfg+4cpKf8tH/79Btxu/xiWZxt/PZ0HsXCjZNr0aaPmJbvQmQZ/KJZ9rVmlDCKGH2zuXL3su4MRziBFzMjBKEQ+/ar7ejd7gXUokUCokCjTWNCX4FhMygqacoxXxE8T+KgckhQJIB/Pc7fF6A7Qs8AuPmRwDnC3ycPBtY/wAfFAqnj8HIXTu9GokOzUukgfGBgDoJlz2XI44S1mWvE3MJW9RbUKmqREV+RcC+hIO/PIiL3ouY8E0gQ5oB5zPORL0cQoLQiCJKMU9eTw7xx2fLFIDteeBaBzAYYgXLqg18bkGYRtrwIJCZy9d7pmmIhPD5fegZ6hGDgbDq6ObwzZDtM2WZ4u7lzUq+gM59BfdBnaWed/eyQqbAuG+cfx4pbWwjyUWBIhX4p/g/n/wj//XsvQvCaEFVzgcEWqaaELNHCcJS1CvsFfHNe6512eugUWqwWblZ3KxWnl++6DOOHmUeRfutdgDAZ4s+u+jXQUgs0DtOKpAp+FFF1fNAYRVQXMXXXZBl0jLVOJvyT6F3sDeoXsKtkVsh28/NJWxVbUWluhLqbHVMC+h8dfNXAQCeUQ8Obj8Ys+clZDEoRxGlmOco/uELwL1bQF4xcLCVppHiaGB8IKCa2iXvJbhYV/hRQs46fvfyrFxCeV45suR0EipZWShQRCku9Sj63YCylKaUYmTKP4WewZ6AaaNL3ku4PXI7ZPssWRbK88v5egkqPpdQqaqEKktFNRAIAQWKqKXNzuw0NumbXPBuYnaMnTkFdTowuFgXJvwTIdtvyNmAivwKbFZtxiblJlSqp1ccyTLpeGxCwqBAQVLKt3/zbZzsPYkH1zyIn3zhJ+L1uaOES95LuOS5hDujd0I+z+xRgjh1pK6EKlNFR1oQEiUKFCSl7Dq2C3KJHFKpFH++48/F6moRRwm5G6DJ14DJZ7BFzVdVK8srQ6Ysk6aOCIkBChQkaab8U+ge6A5YcdR2vQ2+uZsLp2VnZKM8r5yfOppehrpVvRWqLBXkUjlNHRESJ5Q9JQnhHfMGJZfdrDvsKCFDkoHdG3fPJJhVW1GaV4pMWSZNHRGSYDSiiBIlsyOb9E+Ko4TZgSFcLiE7IxsVeRWoyK/AJuUmHL90HHKZHFKJFD9/9OeQZ1BQICTZaEQRpbQtXBQHnjFP4CjBcwnuATcm5x5vPm1j7kZU5FcE1F8uXV2KrIwscero7UtvY2hiCBO+iZhuYCOELB6NKKK0EkcUwihh9s7ly97LuDsaohASZkYJQmW1raqt2KzaDGWmEhmyDMil4UcJD731EIYmh6CQKtBxsIPyDoSkAAoUJED/aH/AtNElLz9KmPJPhWxfmFuIivyKmQSzmh8lCLmEaEcFH/d/jPe638PDRQ+jZkNNLF4SIWSJKFCsUJP+SXQNdIkV1T7xfIJPvZ+if6w/ZPvsjGw+IORVgFEy2KLagi2qLcjPzIdcKqcEMyHLGOUokuza0DW8duE1fHXzV7G9YHtc7tE/2o9LXj4gXPRcxCXvJXQNdIUcJUggEXMJwg7mrcqtKMkrEQNChpR+bQhZSWhEkWQ1b9RgzDcGmUSG88+eX9JzTfom4R5wi1NHFz0Xcdl7GZ4xT8j2ORk5qMivgCZfA02+BltUW8AoGeRn5kMhU9DeBEIIABpRJJ1QBc3H+TA6OYpsefaCHnd39O5MHmE6l9A92B12lFC4qlCcOhKOtChZVQKFTDFvgpkQsrJRoEgyBRSYAL/pLFSQmD1K+MTziZhT8I57Qz5frjxX3JcgJJg3qzcjT5HHTx1J5XSsBSEkKhQokmwCE8iSZSFXnoub927CPeDGRc9Fcdqoe7A75JEWs0cJs6eOilYViSuOaOqIEBILlKNIsgdefwBZsixM+Cfg5/wh2+TKc1GeVw5NvkbcwbxFtQWrFashl8pp6ogQElc0okgyKaRinkIKKTau2hgwStis2oyiVUVicpmmjgghiUaBIsn2l+7HZe9lrMteh+9+5rtYpVgFuVROq44IISmDpp6SbHRqFJc8l8AoGaxWrE52dwghJAidurYELMvi5ZdfDjj7KVrZGdkoV5TjR3/9o0U9z2L7EO3jtm/fju3b47MhcCWIxe9KsqRC3xPVh3jcJ1bPuZTnWXIfOLJoXV1dHACuq6srac+z2MdG+7j77ruPu++++6LuH+HF6nclGVKh74nqQzzuk87vEwIaURBCCImIAgUhhJCIKFAsgVKpxEsvvbTkehRLeZ7FPjZWfScLk84/71Toe6L6EI/7pPP7hIBWPZEFERLZFy5cSHJPCCGJRiMKQgghEdGIghBCSEQ0oiCEEBIRBQpCCCERUaAghBASEQUKQgghEVGgIDHT2toKk8mU7G4QQmKMAgWJCZZlKUgQskxRoCAx0dzcDL1en+xuEELigAIFWbLW1lYYjUYUFBQkuyuEkDigQLGCsSyL1tZWMAwDt9sdso3VaoXRaERLSwuMRiMcDkfQc3g8Hmg0mkR0mRCSBFQKdYVyOp1wOBzQaDRwu90hC5pYrVbY7XbYbDbxGsMwsNls0Gq1APgpJ7PZnKhuE0KSgI7wWOFYloVKpUJnZ6f45i9gGAZmszkg92AymeB2u2Gz2dDa2gq73Q6GYQAAdrsdHo8H9fX1aGxsTOjrIITED40oSEgsy8LtdgcdS8wwDFpaWgAAer0+IIi4XC5otVoKEoQsM5SjICF1dHQAANRqdcD1uV8LWlpa0NHRAYfDIQYSQsjyQCMKEpKQswhX6IRl2YDvNTY20kiCkGWKRhQkpHABwuPxJLYjhJCko0BBQhKWu85dDTXfSIMQsvxQoCAhaTQacensbC6Xi3ZgE7LCUKBY4SJNJZlMJlgsloBrDocDTU1N8e4WISSF0D6KFcrtdsNisYgb73Q6HbRaLerr6wP2U1itVnR2doJhGLhcLhgMBuh0uiT2nBCSaBQoCCGERERTT4QQQiKiQEEIISQiChSEEEIiokBBCCEkIgoUhBBCIqJAQQghJCIKFIQQQiKiQEEIISQiChSEEEIiokBBCCEkIgoUhBBCIqJAQRKitbUVtbW1kEgkMJlMIds4nU4YjUZIJBIYDAY4nc4E9zI8q9Uq9t9gMERsy7IsJBIJVCoVjEYj3G73ol9/a2srDAaD+HwmkymgRkhLSwsMBgNqa2vBMAwMBgMcDgcAwGAwgGEYVFVVoaqqChKJBBKJBFVVVaitrRWvx5PD4UBtbW1c70ESgCMkQVwuFweAA8DZbLaw7TQaTQJ7tXCdnZ1i/71eb9h2jY2NHABOp9MFXF/K6wfANTQ0BFzT6XSc2WwOundjYyPHcVxQe61WyymVyoBrc9vEms1m45RKZcSfVyR2u33RjyWxQyMKklB6vR4ajQYGgyGoKJIgVavnKZVKsW/Nzc1h27W2toZ9DbF6/a2trejo6AiqU240GlFQUAAAQZ/k1Wo11Gp1wLV4jyj0ej28Xu+i/k6FURxJPgoUJKHUajXsdjuA4DeydFBXVweNRoOWlpaQ329paYHRaAx6QxbE6vW3t7eHvIdGoxGDx0IqETY0NCy6D/EmvL5U/eCwklCgIAmn0Whgs9ngdrvnne93OBwwGAxQqVTiNbfbDZPJBIZh0NraGtSOZVnx/2tra+F2uwOuVVVVhf00vxBCjiFUsLBYLEGf8ueK5vWHU1BQIP4c5tY1XwohJ1JVVSX+zCQSCRiGCfl6rVYrTCYTjEYjqqqqAvIvQs5F+DsR2gvP73a7xdxJVVVVUE6qvb1d7JPVagXLsmL/jEYjamtrw+Z7SIwle+6LrBwulytgTryhoYEDwFksloB2Wq1W/H+v1yu2m33NZrMFzPV3dnZyer1enMt3uVycy+XilEolp9Fogq7Nvsdi+q9UKoPm+202m5gf0Gg0IXMU0b5+AULkKLRarZjzEO43N2cxm06nmzf/43K5OL1ezymVSq6hoYHr7OzkXC4XZzabOQDi6+M4Ph8yu09er5dTKpXi6/Z6vWK+RsgzuFwuTqfTcUqlkmtsbOS8Xi/ncrk4jUYT8PO02+3i36fZbObMZjN37ty5oJ+5Xq+P+HpIbFCgIAkz942S42be7Do7OwOuzSYEhdm8Xm9QUlhoNzv5OfeNiuO4oMCzmP5bLJag+2u1WvE+CwkUwmPme/0cFzpQcBwfIM1ms/jmLgSNUBYSKDhu5mc21+yfm5CYn91vjpv5uQjX7XZ7yL+TUEEWAOdyuYKuRWoz+/9J/NDUE0mqkydPQqlUYv/+/VFNoUSat579PYZhgq7FIoHb0NAApVIpJrUdDgeqq6ujnk9f7OsXaLVaNDY2wmazwev1orGxEW63W1wiuxgFBQUhX4eQU2FZVpwm0mg0AW2Eeurz3X9ufmXu84Si0+mg0WjAMAwYhoHRaIzptBsJjwIFSSqlUgmbzSbOh8daqIRvuERztJqamuB0OuFwOGA2mxc1Xx7r19/U1AQgdq9xLmHll8fjAQDxv4J4JqCVSiVcLhdsNht0Oh1OnDgRMrdBYo8CBUk6nU6HxsZGOByOiG+2sz89psKbQ2NjI5RKJYxGI5RK5YI+FYey0Nc/m9VqDfkzcLvdUCqV0Gq1i+pLJMePHxdXSVVXVwMIHjkIgUMYWcSS0+kEy7LQ6/WwWCzistuOjo6Y34sEokBBEmruJ1CB2WyGTqdDS0tL0Iok4Q3YarUC4N8wFvqGGu5+c5lMpgXtuJ6roaEBbrdb/CQ/330X8/rD9WXuXg63241Dhw7h5MmTYR+30J+HsMJIIPy/2WwGwE956XS6oFVXZrMZZrN50UEz1M/Y4XDA6XQGTamxLAu1Wo26urpF3YtEIdlJErIy2Gw2cWewsNplLmHVzNxEJ8fxiVQhWavT6cRkqlar5Ww2G2exWMTEsF6v51wuV8A1YQWP0Je51/R6fVDyeTbhueb23+v1BjxOWDWE6dVIwmqrxb5+m83G6XQ6DkDAYzs7O8XVQ1qtltPr9VxDQ0PI5xWS6EKf9Ho9Z7fbw75WYYWT8Pw6nS7sTvKGhgZOq9WKO8JnP6/wmoV7Col34e9xdn+FXe86nU78O/F6vQF/f3a7ndNqtZxOp+MaGhrEny2JPwnHcVwS4hMhYTmdzrhMnaSLZL/+lpYWNDc3w+v1Jq0PJLXQ1BNJOSs5SAD0+knqoUBBCAnQ39+f7C6QFEOBghASgGVZ8Q8hAEA5CkKISDhTSSCcdktWNgoUhBBCIqKpJ0IIIRFRoCCEEBIRBQpCCCERUaAghBASEQUKQgghEVGgIIQQEhEFCkIIIRFRoCCEEBIRBQpCCCERUaAghBASEQUKQgghEVGgIIQQEhEFCkIIIRFRoCCEEBLR/wMs5O8Co4/SpQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 410.433x253.648 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib.ticker import FixedLocator\n",
    "\n",
    "# Fit robust linear regression in log-log space\n",
    "fig = sized_fig(1/2)\n",
    "ax = plt.gca()\n",
    "for i, alg in enumerate([\"kmst\", \"dmst\", \"umap\"]):\n",
    "    alg_values = values[values[\"algorithm\"] == alg]\n",
    "    sns.regplot(\n",
    "        x=np.log10(alg_values[\"sample_fraction\"] * df.shape[0]),\n",
    "        y=np.log10(alg_values[\"compute_time\"]),\n",
    "        ci=95,\n",
    "        order=1,\n",
    "        robust=True,\n",
    "        color=f\"C{i}\",\n",
    "        units=alg_values[\"repeat\"],\n",
    "        scatter_kws={\"edgecolor\": \"none\", \"linewidths\": 0, \"s\": 2},\n",
    "        line_kws={\"linewidth\": 1},\n",
    "        ax=ax,\n",
    "    )\n",
    "ax.set_xlabel(\"Num. MNIST points\")\n",
    "ax.set_ylabel(\"Run time (s)\")\n",
    "\n",
    "# Draw log y-ticks\n",
    "y_ticks = np.array([0.0, 1.0, 2.0, 3.0])\n",
    "plt.ylim(-1, plt.ylim()[1])\n",
    "ax.set_yticks(y_ticks)\n",
    "ax.get_yaxis().set_major_formatter(lambda x, pos: f\"$10^{{{int(x)}}}$\")\n",
    "ax.get_yaxis().set_minor_locator(\n",
    "    FixedLocator(locs=np.concat((\n",
    "        np.log10(np.arange(0.2, 1, 0.1) * 10.0 ** y_ticks[0]),\n",
    "        np.log10(np.arange(2, 10) * 10.0 ** y_ticks[None].T).ravel())\n",
    "    ))\n",
    ")\n",
    "\n",
    "# Draw log x-ticks\n",
    "x_ticks = np.array([4.0])\n",
    "ax.set_xticks(x_ticks)\n",
    "ax.get_xaxis().set_major_formatter(lambda x, pos: f\"$10^{{{int(x)}}}$\")\n",
    "ax.get_xaxis().set_minor_locator(\n",
    "    FixedLocator(locs=np.log10(np.array(\n",
    "        [0.6, 0.7, 0.8, 0.9, 1, 2, 3, 4, 5, 6, 7, 8, 9, 20, 30, 40, 50, 60, 70, 80]\n",
    "    ) * 10.0 ** x_ticks[None].T).ravel())\n",
    ")\n",
    "plt.xlim(np.log10([6000, 80000]))\n",
    "\n",
    "# Legend\n",
    "adjust_legend_subtitles(\n",
    "    plt.legend(\n",
    "        loc=\"upper left\",\n",
    "        handles=[\n",
    "            ml.Line2D([], [], color=f\"C{j}\", label=f\"{v}\")\n",
    "            for j, v in enumerate(['$k$-MST (kd-tree)', '$k$-MST (descent)', '$k$-NN'])\n",
    "        ]\n",
    "    )\n",
    ")\n",
    "plt.subplots_adjust(left=0.17, right=0.9, top=0.95, bottom=0.24)\n",
    "plt.savefig(\"./images/mnist_scaling.pdf\", pad_inches=0)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "work",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}