jacobians.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 313,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import torch\n",
    "\n",
    "def load_jacobian(density_idx, slope_idx, seed, data_dir):\n",
    "    filepath = f\"jacobian_density-{density_idx}_slope-{slope_idx}_seed-{seed}.pt\"\n",
    "    jac_path = os.path.join(data_dir, \"jacobians_7x7\", filepath)\n",
    "    return torch.load(jac_path).double()\n",
    "\n",
    "def load_jtfs(density_idx, slope_idx, seed, data_dir):\n",
    "    filepath = f\"jtfs_density-{density_idx}_slope-{slope_idx}_seed-{seed}.pt\"\n",
    "    jtfs_path = os.path.join(data_dir, \"jtfs_7x7\", filepath)\n",
    "    return torch.load(jtfs_path).double()\n",
    "\n",
    "data_dir = \"icassp2024_scrapl_data\"\n",
    "\n",
    "y_coords = (4, 3)\n",
    "Sy = load_jtfs(*y_coords, seed=0, data_dir=data_dir)[:, 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 358,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "\n",
    "SX, JX, SY = [], [], []\n",
    "for seed in range(1, 10):\n",
    "    y_coords = (4, 3)\n",
    "    x_coords = (5, 3)\n",
    "    Sx = load_jtfs(*x_coords, seed=seed, data_dir=data_dir)[:, 0]\n",
    "    Sy = load_jtfs(*y_coords, seed=seed, data_dir=data_dir)[:, 0]\n",
    "    Jx = load_jacobian(*x_coords, seed=seed, data_dir=data_dir)\n",
    "    SX.append(Sx)\n",
    "    SY.append(Sy)\n",
    "    JX.append(Jx)\n",
    "SX = torch.stack(SX, dim=0)\n",
    "SY = torch.stack(SY, dim=0)\n",
    "JX = torch.stack(JX, dim=0)\n",
    "\n",
    "G_scrapl = []\n",
    "for _ in range(100):\n",
    "    Sx_scrapl, Sy_scrapl, Jx_scrapl = [], [], []\n",
    "    for path in range(SX.shape[-1]):\n",
    "        random_seed = torch.randint(SX.shape[0], (1,))\n",
    "        Jx_scrapl.append(JX[random_seed, path, :])\n",
    "        Sx_scrapl.append(SX[random_seed, path])\n",
    "        Sy_scrapl.append(SY[random_seed, path])\n",
    "    Jx_scrapl = torch.stack(Jx_scrapl, dim=0).squeeze()\n",
    "    Sx_scrapl = torch.stack(Sx_scrapl, dim=0).squeeze()\n",
    "    Sy_scrapl = torch.stack(Sy_scrapl, dim=0).squeeze()\n",
    "    g_scrapl = Jx_scrapl.T @ (Sy_scrapl - Sx_scrapl)\n",
    "    G_scrapl.append(g_scrapl)\n",
    "G_scrapl = torch.stack(G_scrapl, dim=0).squeeze()\n",
    "\n",
    "G_jtfs = []\n",
    "for seed in range(9):\n",
    "    g_jtfs = JX[seed, :, :].T @ (Sy - SX[seed, :])\n",
    "    G_jtfs.append(g_jtfs)\n",
    "G_jtfs = torch.stack(G_jtfs, dim=0).squeeze()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 362,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Histogram density')"
      ]
     },
     "execution_count": 362,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkgAAAGwCAYAAABSN5pGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABUgUlEQVR4nO3deVwVVf8H8M+97MgOwQWlUDF3RUGvqLkkCmopRqZmKcqDlZIlWYYproWZC5Uaj5lgTxJkmZkahiRqiZK45UZqGKYsLgkBynLv/P7g58S9gHLhwrB83q/XvJo5c+bMd6YBv5w5MyMTBEEAEREREYnkUgdARERE1NgwQSIiIiLSwgSJiIiISAsTJCIiIiItTJCIiIiItDBBIiIiItLCBImIiIhIi6HUATRVarUa169fh6WlJWQymdThEBERUQ0IgoB//vkHLi4ukMur7ydiglRL169fh6urq9RhEBERUS1cvXoVbdq0qXY9E6RasrS0BFB+gq2srCSOhoiIiGoiPz8frq6u4r/j1WGCVEv3b6tZWVkxQSIiImpiHjY8hoO0iYiIiLQwQSIiIiLSwgSJiIiISAvHIBEREdWBWq1GSUmJ1GHQ/zMyMoKBgUGd22GCREREVEslJSXIyMiAWq2WOhSqwMbGBgqFok7vKWSCREREVAuCICArKwsGBgZwdXV94EsHqWEIgoCioiLk5uYCAJydnWvdFhMkIiKiWigrK0NRURFcXFxgbm4udTj0/8zMzAAAubm5cHR0rPXtNqa7REREtaBSqQAAxsbGEkdC2u4nrKWlpbVugwkSERFRHfB7nI2PPv6fMEEiIiIi0sIEiYiIiEgLB2kTERHp0drE3xt0f3OGP96g+2sp2INERETUggQGBsLf3x8ymeyB0+LFi3HlypUq173wwgsAygeqr1ixAp06dYKZmRns7OygVCqxadMmiY+y7tiDRERE1AJlZWWJ8/Hx8QgPD0d6erpYZmFhgZs3bwIA9u3bh65du4rr7j9Kv2TJEvz3v//FunXr4OXlhfz8fBw7dgx///13Ax1F/WGCRERE1AIpFApx3traGjKZTKMMgJgg2dvbV1oHADt37sTMmTMxfvx4saxnz571FHHD4i02IiIiqhWFQoGffvoJN27ckDoUvWOCRERERA/Uv39/WFhYiNOJEycAAGvWrMGNGzegUCjQo0cPvPzyy/jhhx8kjlY/eIuNiKiCotIiKGOVAICjzx+FuRE/IUEUHx+Pzp07i8uurq4AgC5duuDMmTNIS0vDL7/8goMHD+Lpp59GYGBgkx+ozR4kIiIieiBXV1e4u7uLk4mJibhOLpejT58+eP3117F9+3bExMTgs88+Q0ZGhoQR1x0TJCIiItKbLl26AAAKCwsljqRueIuNiIiIauXZZ5/FgAED0L9/fygUCmRkZCAsLAyPP/44OnXqJHV4dcIEiYiISI8a+5ut1Wo1DA3188+/r68vvvzyS0RERCAvLw8KhQJPPvkkFi9erLd9SKVpR09EREQ6yc3Nhbu7u0ZZYGAgAgMDK9V1c3ODIAjVthUcHIzg4GB9h9gocAwSERFRC/D3339j165dSE5Oho+Pj9ThNHrsQSIiqsBAboAnWj8hzhM1F9OnT8evv/6KN954A2PHjpU6nEaPCRIRUQUmBibY4LNB6jCI9O7bb7+VOoQmhbfYiIiIiLQwQSIiIiLSwgSJiKiCotIi9N3aF3239kVRaZHU4RCRRDgGiYhIy92yu1KHQEQSaxQ9SOvXr4ebmxtMTU2hVCqRmppabd1PP/0UTzzxBGxtbWFrawsfH59K9QVBQHh4OJydnWFmZgYfHx9cvHhRo87t27cxefJkWFlZwcbGBkFBQSgoKKiX4yMiIqKmRfIepPj4eISGhiIqKgpKpRKRkZHw9fVFeno6HB0dK9VPTk7GpEmT0L9/f5iamuL999/HiBEjcPbsWbRu3RoAsHLlSnz00UfYsmUL2rZti4ULF8LX1xfnzp2DqakpAGDy5MnIyspCYmIiSktLMW3aNMyYMQOxsbENevxEzdL+iLq3MTSs7m0QEdWSTHjQKzIbgFKpRJ8+fbBu3ToA5a9Ad3V1xauvvoq33377odurVCrY2tpi3bp1mDJlCgRBgIuLC9544w3MnTsXAJCXlwcnJyfExMRg4sSJOH/+PLp06YJff/0VXl5eAICEhASMGjUKf/31F1xcXB663/z8fFhbWyMvLw9WVlZ1OANEzVATTpCKSougjFUCAI4+fxTmRuaSxEGN371795CRkYG2bduKf3wD0M/1r4ta/KzcuHED4eHh2L17N3JycmBra4uePXsiPDwcAwYMAACcOHEC7733Hg4ePIi8vDy4urpiyJAhePPNN/H444/jypUraNu2rdimra0tunfvjuXLl+OJJ56otM+XXnoJmzZtQlxcHMaPH6+xbvHixViyZAkAwMDAAG3atMG4ceOwbNkyWFhYiPs6ceIEPDw8Hnp81f6/Qc3//Zb0FltJSQnS0tI03ugpl8vh4+ODlJSUGrVRVFSE0tJS2NnZAQAyMjKQnZ2t0aa1tTWUSqXYZkpKCmxsbMTkCAB8fHwgl8tx9OjRKvdTXFyM/Px8jYmIiKgpCggIwIkTJ7Blyxb8/vvv2LlzJ4YMGYJbt24BAHbt2oV+/fqhuLgYW7duxfnz5/HFF1/A2toaCxcu1Ghr3759yMrKwsGDB+Hi4oKnnnoKOTk5GnWKiooQFxeHt956C5s3b64ypq5duyIrKwtXrlzB+++/j40bN+KNN96onxNQA5LeYrt58yZUKhWcnJw0yp2cnHDhwoUatTFv3jy4uLiICVF2drbYhnab99dlZ2dXun1naGgIOzs7sY62iIgIMbslIiJqqu7cuYNDhw4hOTkZgwcPBgA89thj6Nu3L4DyZGbatGkYNWqUxssl27ZtC6VSiTt37mi0Z29vD4VCAYVCgfnz5yMuLg5Hjx7FmDFjxDrbtm1Dly5d8Pbbb8PFxQVXr16Fq6urRjuGhoZQKBQAgAkTJiApKQk7d+7Ef//73/o4DQ/VKAZp19aKFSsQFxeHb7/9tlIXmr6FhYUhLy9PnK5evVqv+yMiachlcng5ecHLyQtyWZP+FUlUJQsLC1hYWGDHjh0oLi6utH7v3r24efMm3nrrrSq3t7GxqbL87t27+PzzzwEAxsbGGus+++wzvPDCC7C2tsbIkSMRExPz0DjNzMxQUlLy0Hr1RdKffgcHBxgYGFTqisvJyRGzyOqsWrUKK1aswI8//ogePXqI5fe3e1CbCoUCubm5GuvLyspw+/btavdrYmICKysrjYmImh9TQ1NE+0Uj2i8apob1+4cXkRQMDQ0RExODLVu2wMbGBgMGDMD8+fNx+vRpABCf+u7UqVON2uvfvz8sLCzQqlUrrFq1Cp6enhg2bJi4/uLFizhy5AgmTJgAAHjhhRcQHR2NBw2BTktLQ2xsLJ588snaHmadSZogGRsbw9PTE0lJSWKZWq1GUlISvL29q91u5cqVWLZsGRISEjTGEQHlXYAKhUKjzfz8fBw9elRs09vbG3fu3EFaWppY56effoJarYZSqdTX4RERETVKAQEBuH79Onbu3Ak/Pz8kJyejd+/eiImJeWDiUpX4+HicOHEC33zzDdzd3RETEwMjIyNx/ebNm+Hr6wsHBwcAwKhRo5CXl4effvpJo53ffvsNFhYWMDMzQ9++feHt7S0+wCUFyR/zDw0NxdSpU+Hl5YW+ffsiMjIShYWFmDZtGgBgypQpaN26NSIiyp8KeP/99xEeHo7Y2Fi4ubmJY4budxnKZDK8/vrrWL58OTp06CA+5u/i4gJ/f38AQOfOneHn54fg4GBERUWhtLQUISEhmDhxYo2eYCMiImrqTE1NMXz4cAwfPhwLFy7Ef/7zHyxatAiRkZEAgAsXLjyws+I+V1dXdOjQAR06dEBZWRnGjRuHM2fOwMTEBCqVClu2bEF2djYMDf9NOVQqFTZv3qzR09SxY0fs3LkThoaGcHFxqXSbrqFJfoN9woQJWLVqFcLDw+Hh4YGTJ08iISFBHGSdmZmJrKwssf4nn3yCkpISPPvss3B2dhanVatWiXXeeustvPrqq5gxYwb69OmDgoICJCQkaIxT2rp1Kzp16oRhw4Zh1KhRGDhwIDZu3NhwB05EjVJRaREGxQ3CoLhB/NQItShdunRBYWEhRowYAQcHB6xcubLKetqDtCt69tlnYWhoiA0bNgAA9uzZg3/++QcnTpzAyZMnxenLL7/E9u3bNdoyNjaGu7s73NzcJE+OgEbQgwQAISEhCAkJqXJdcnKyxvKVK1ce2p5MJsPSpUuxdOnSauvY2dnxpZBEVKW/i/+WOgSienPr1i2MHz8e06dPR48ePWBpaYljx45h5cqVGDt2LFq1aoVNmzZh/PjxGDNmDGbPng13d3fcvHkTX331FTIzMxEXF1dl2zKZDLNnz8bixYvx0ksv4bPPPsPo0aPRs2dPjXpdunTBnDlzsHXrVsyaNavGsaenp1cq69q1q8YtPX2RvAeJiIiIGo6FhQWUSiXWrl2LQYMGoVu3bli4cCGCg4PFMT9jx47F4cOHYWRkhOeffx6dOnXCpEmTkJeXh+XLlz+w/alTp6K0tBQff/wxdu/ejYCAgEp15HI5xo0bh88++0yn2CdOnIhevXppTNoPZemL5G/Sbqr4Jm2iB+CbtKkFeNDbmklaTf5N2kRERESNERMkIiIiIi1MkIiIiIi0NIqn2IiIGgu5TI6u9l3FeSJqmZggERFVYGpoirinqn6EmYhaDv55RERERKSFCRIRERGRFiZIREQV3C27C9+vfeH7tS/ult2VOhwikgjHIBERVSAIAq4XXhfniahlYg8SERERkRYmSERERC1IYGAgZDKZONnb28PPzw+nT5/WS/uLFy+Gh4eHXtqSEhMkIiKiFsbPzw9ZWVnIyspCUlISDA0N8dRTT0kdVqPCBImIiKiFMTExgUKhgEKhgIeHB95++21cvXoVN27cAABcvXoVzz33HGxsbGBnZ4exY8fiypUr4vbJycno27cvWrVqBRsbGwwYMAB//vknYmJisGTJEpw6dUrsoYqJiZHmIOuIg7SJiIj0qKi0qNp1BnIDmBiY1KiuXCaHqaHpQ+uaG5nXIsp/FRQU4IsvvoC7uzvs7e1RWloKX19feHt749ChQzA0NMTy5cvF23ByuRz+/v4IDg7Gl19+iZKSEqSmpkImk2HChAk4c+YMEhISsG/fPgCAtbV1neKTChMkIqIKZDIZ2lu3F+eJdKWMVVa77onWT2CDzwZxechXQ6p9nYSXkxei/aLFZb9v/PB38d+V6v029TedY9y1axcsLCwAAIWFhXB2dsauXbsgl8sRGxsLtVqNTZs2iT8D0dHRsLGxQXJyMry8vJCXl4ennnoK7duX/6x07txZbNvCwgKGhoZQKBQ6x9WYMEEiIqrAzNAMO/x3SB0GUb0aOnQoPvnkEwDA33//jQ0bNmDkyJFITU3FqVOncOnSJVhaWmpsc+/ePVy+fBkjRoxAYGAgfH19MXz4cPj4+OC5556Ds7OzFIdSb5ggERER6dHR549Wu85AbqCxnPxccrV1tT+WnBCQUKe4KmrVqhXc3d3F5U2bNsHa2hqffvopCgoK4Onpia1bt1ba7pFHHgFQ3qM0e/ZsJCQkID4+HgsWLEBiYiL69euntxilxgSJiIhIj3QZE1RfdXUlk8kgl8tx9+5d9O7dG/Hx8XB0dISVlVW12/Tq1Qu9evVCWFgYvL29ERsbi379+sHY2BgqlareYm0ofIqNiKiCu2V34b/DH/47/PmpEWq2iouLkZ2djezsbJw/fx6vvvoqCgoK8PTTT2Py5MlwcHDA2LFjcejQIWRkZCA5ORmzZ8/GX3/9hYyMDISFhSElJQV//vknfvzxR1y8eFEch+Tm5oaMjAycPHkSN2/eRHFxscRHWzvsQSIiqkAQBFzOuyzOEzVHCQkJ4pghS0tLdOrUCdu2bcOQIUMAAAcPHsS8efPwzDPP4J9//kHr1q0xbNgwWFlZ4e7du7hw4QK2bNmCW7duwdnZGbNmzcJLL70EAAgICMD27dsxdOhQ3LlzB9HR0QgMDJToSGuPCRIREVELEhMT89B3EykUCmzZsqXKdVZWVvj222+r3dbExARff/11XUJsFHiLjYiIiEgLEyQiIiIiLUyQiIiIiLQwQSIiIiLSwkHaREQVyGQyuLRyEeeJHoZPOzY++vh/wgSJiKgCM0Mz7H12r9RhUBNgYFD+VuySkhKYmZlJHA1VVFRU/mFfIyOjWrfBBImIiKgWDA0NYW5ujhs3bsDIyAhyOUetSE0QBBQVFSE3Nxc2NjZiElsbTJCIiIhqQSaTwdnZGRkZGfjzzz+lDocqsLGxgUKhqFMbkidI69evxwcffIDs7Gz07NkTH3/8Mfr27Vtl3bNnzyI8PBxpaWn4888/sXbtWrz++usaddzc3Kq8UGfOnIn169cDAIYMGYIDBw5orH/ppZcQFRWln4MioibrXtk9BCYEAgBi/GJgamgqbUDUqBkbG6NDhw4oKSmROhT6f0ZGRnXqObpP0gQpPj4eoaGhiIqKglKpRGRkJHx9fZGeng5HR8dK9YuKitCuXTuMHz8ec+bMqbLNX3/9VeMjeWfOnMHw4cMxfvx4jXrBwcFYunSpuGxuXn8fASSipkMtqHH21llxnuhh5HI5TE2ZSDc3kiZIa9asQXBwMKZNmwYAiIqKwu7du7F582a8/fbbler36dMHffr0AYAq1wPAI488orG8YsUKtG/fHoMHD9YoNzc316n7rbi4WOODe/n5+TXeloiIiJoWyUaUlZSUIC0tDT4+Pv8GI5fDx8cHKSkpetvHF198genTp1d6XHfr1q1wcHBAt27dEBYWJo54r05ERASsra3FydXVVS8xEhERUeMjWQ/SzZs3oVKp4OTkpFHu5OSECxcu6GUfO3bswJ07dyp9Rfj555/HY489BhcXF5w+fRrz5s1Deno6tm/fXm1bYWFhCA0NFZfz8/OZJBERETVTkg/Srk+fffYZRo4cCRcXF43yGTNmiPPdu3eHs7Mzhg0bhsuXL6N9+/ZVtmViYgITE5N6jZeIiIgaB8lusTk4OMDAwAA5OTka5Tk5OXV+NA8A/vzzT+zbtw//+c9/HlpXqVQCAC5dulTn/RIREVHTJ1kPkrGxMTw9PZGUlAR/f38AgFqtRlJSEkJCQurcfnR0NBwdHTF69OiH1j158iQAwNnZuc77JZLc/oi6tzE0rO5tNGG2JrZSh0BEEpP0FltoaCimTp0KLy8v9O3bF5GRkSgsLBSfapsyZQpat26NiIjyX/glJSU4d+6cOH/t2jWcPHkSFhYWcHd3F9tVq9WIjo7G1KlTYWioeYiXL19GbGwsRo0aBXt7e5w+fRpz5szBoEGD0KNHjwY6ciJqrMyNzHFw4kGpwyAiiUmaIE2YMAE3btxAeHg4srOz4eHhgYSEBHHgdmZmpsar269fv45evXqJy6tWrcKqVaswePBgJCcni+X79u1DZmYmpk+fXmmfxsbG2Ldvn5iMubq6IiAgAAsWLKi/AyUiIqImRSbwM8S1kp+fD2tra+Tl5cHKykrqcIj+1RhusTWGGIiIqlDTf7/5ZT0iogruld3DtIRpmJYwDffK7kkdDhFJpFk/5k9EpCu1oMaxnGPiPBG1TOxBIiIiItLCBImIiIhICxMkIiIiIi1MkIiIiIi0MEEiIiIi0sKn2IiItJgZmkkdAhFJjAkSEVEF5kbmSJ2cKnUYRCQx3mIjIiIi0sIEiYiIiEgLEyQiogqKVcWYuW8mZu6biWJVsdThEJFEOAaJiKgClVqFQ9cOifMwkDggIpIEe5CIiIiItDBBIiIiItLCBImIiIhICxMkIiIiIi1MkIiIiIi0MEEiIiIi0sLH/ImIKjA3MsdvU3+TOgwikhh7kIiIiIi0MEEiIiIi0sIEiYiogmJVMUKTQxGaHMpPjRC1YEyQiIgqUKlVSPwzEYl/JpZ/aoSIWiQmSERERERamCARERERaWGCRERERKSFCRIRERGRFr4okkif9kfUvY2hYXVvg4iI6oQ9SERERERaJE+Q1q9fDzc3N5iamkKpVCI1NbXaumfPnkVAQADc3Nwgk8kQGRlZqc7ixYshk8k0pk6dOmnUuXfvHmbNmgV7e3tYWFggICAAOTk5+j40ImqCzAzNcPT5ozj6/FGYGZpJHQ4RSUTSBCk+Ph6hoaFYtGgRjh8/jp49e8LX1xe5ublV1i8qKkK7du2wYsUKKBSKatvt2rUrsrKyxOnnn3/WWD9nzhx8//332LZtGw4cOIDr16/jmWee0euxEVHTJJPJYG5kDnMjc8hkMqnDISKJSJogrVmzBsHBwZg2bRq6dOmCqKgomJubY/PmzVXW79OnDz744ANMnDgRJiYm1bZraGgIhUIhTg4ODuK6vLw8fPbZZ1izZg2efPJJeHp6Ijo6GocPH8aRI0f0foxERETU9EiWIJWUlCAtLQ0+Pj7/BiOXw8fHBykpKXVq++LFi3BxcUG7du0wefJkZGZmiuvS0tJQWlqqsd9OnTrh0UcffeB+i4uLkZ+frzERUfNToirBOz+/g3d+fgclqhKpwyEiiUiWIN28eRMqlQpOTk4a5U5OTsjOzq51u0qlEjExMUhISMAnn3yCjIwMPPHEE/jnn38AANnZ2TA2NoaNjY1O+42IiIC1tbU4ubq61jpGImq8ytRl2Hl5J3Ze3okydZnU4RCRRCQfpK1vI0eOxPjx49GjRw/4+vpiz549uHPnDr766qs6tRsWFoa8vDxxunr1qp4iJiIiosZGsvcgOTg4wMDAoNLTYzk5OQ8cgK0rGxsbPP7447h06RIAQKFQoKSkBHfu3NHoRXrYfk1MTB447omIiIiaD517kAYPHozPP/8cd+/erdOOjY2N4enpiaSkJLFMrVYjKSkJ3t7edWq7ooKCAly+fBnOzs4AAE9PTxgZGWnsNz09HZmZmXrdLxERETVdOidIvXr1wty5c6FQKBAcHFynJ79CQ0Px6aefYsuWLTh//jxeeeUVFBYWYtq0aQCAKVOmICzs37cKl5SU4OTJkzh58iRKSkpw7do1nDx5UuwdAoC5c+fiwIEDuHLlCg4fPoxx48bBwMAAkyZNAgBYW1sjKCgIoaGh2L9/P9LS0jBt2jR4e3ujX79+tT4WIiIiaj50vsUWGRmJVatWYefOndiyZQsGDRoEd3d3TJ8+HS+++GKlQdcPMmHCBNy4cQPh4eHIzs6Gh4cHEhISxDYyMzMhl/+bw12/fh29evUSl1etWoVVq1Zh8ODBSE5OBgD89ddfmDRpEm7duoVHHnkEAwcOxJEjR/DII4+I261duxZyuRwBAQEoLi6Gr68vNmzYoOupICIiomZKJgiCUJcGcnNzsXHjRrz77rtQqVQYNWoUZs+ejSeffFJfMTZK+fn5sLa2Rl5eHqysrKQOhxqLxvAtNsZQJ0WlRVDGKgEAR58/CnMjc0niIKL6UdN/v+s0SDs1NRXR0dGIi4uDo6MjAgMDce3aNTz11FOYOXMmVq1aVZfmiYganJmhGQ5MOCDOE1HLpHOClJubi//973+Ijo7GxYsX8fTTT+PLL7+Er6+v+Fr+wMBA+Pn5MUEioiZHJpPBztRO6jCISGI6J0ht2rRB+/btMX36dAQGBmqM7bmvR48e6NOnj14CJCIiImpoOidISUlJeOKJJx5Yx8rKCvv37691UEREUilRlWDlrysBAG/1eQvGBsbSBlTX8VwSjeUiaup0fsx/0aJFuHPnTqXy/Pz8Zj8wm4iavzJ1GeLT4xGfHs9PjRC1YDonSAcOHEBJSeUPON67dw+HDh3SS1BEREREUqrxLbbTp08DAARBwLlz5zQ+7KpSqZCQkIDWrVvrP0IiIiKiBlbjBMnDwwMymQwymazKW2lmZmb4+OOP9RocERERkRRqnCBlZGRAEAS0a9cOqampGk+vGRsbw9HREQYGBvUSJBEREVFDqnGC9NhjjwEo/6AsERERUXNWowRp586dGDlyJIyMjLBz584H1h0zZoxeAiMiIiKSSo0SJH9/f2RnZ8PR0RH+/v7V1pPJZFCpVPqKjYiowZkamiIhIEGcJ6KWqUYJUsXbarzFRkTNmVwmR2sLPpFL1NLp/B6kqlT14kgiIiKipkrnBOn9999HfHy8uDx+/HjY2dmhdevWOHXqlF6DIyJqaKWqUqw+thqrj61GqapU6nCISCI6J0hRUVFwdXUFACQmJmLfvn1ISEjAyJEj8eabb+o9QCKihlSqLkXM2RjEnI1BqZoJElFLpfPHarOzs8UEadeuXXjuuecwYsQIuLm5QalU6j1AIiIiooamcw+Sra0trl69CgBISEiAj48PgPJPkPAJNiIiImoOdO5BeuaZZ/D888+jQ4cOuHXrFkaOHAkAOHHiBNzd3fUeIBEREVFD0zlBWrt2Ldzc3HD16lWsXLkSFhYWAICsrCzMnDlT7wESERERNTSdEyQjIyPMnTu3UvmcOXP0EhARERGR1HROkADg4sWL2L9/P3Jzcyu9ODI8PFwvgRERERFJRecE6dNPP8Urr7wCBwcHKBQKyGQycZ1MJmOCRERNmqmhKb4d8604T0Qtk84J0vLly/Huu+9i3rx59REPEZGk5DI53G35wAlRS6fzY/5///03xo8fXx+xEBERETUKOidI48ePx48//lgfsRARSa5UVYoNJzdgw8kN/NQIUQum8y02d3d3LFy4EEeOHEH37t1hZGSksX727Nl6C46IqKGVqkvxyalPAACBXQNhZGD0kC2IqDnSOUHauHEjLCwscODAARw4cEBjnUwmY4JERERETZ7OCVJGRkZ9xEFERETUaNTqPUgAUFJSgoyMDLRv3x6GhrVuhmpobeLvDbavOcMfb7B9NbT6Po/9Mm+J897t7Ot1X9T01eR6rHhN1caRsvJ9NOefa6L6oPMg7aKiIgQFBcHc3Bxdu3ZFZmYmAODVV1/FihUr9B4gERERUUPTOUEKCwvDqVOnkJycDFPTf1+i5uPjg/j4eJ0DWL9+Pdzc3GBqagqlUonU1NRq6549exYBAQFwc3ODTCZDZGRkpToRERHo06cPLC0t4ejoCH9/f6Snp2vUGTJkCGQymcb08ssv6xw7ERERNU86J0g7duzAunXrMHDgQI23aHft2hWXL1/Wqa34+HiEhoZi0aJFOH78OHr27AlfX1/k5uZWWb+oqAjt2rXDihUroFAoqqxz4MABzJo1C0eOHEFiYiJKS0sxYsQIFBYWatQLDg5GVlaWOK1cuVKn2ImIiKj50nnw0I0bN+Do6FipvLCwUCNhqok1a9YgODgY06ZNAwBERUVh9+7d2Lx5M95+++1K9fv06YM+ffoAQJXrASAhIUFjOSYmBo6OjkhLS8OgQYPEcnNz82qTrKoUFxejuLhYXM7Pz6/xtkTUdJgYmODL0V+K80TUMuncg+Tl5YXdu3eLy/eTok2bNsHb27vG7ZSUlCAtLQ0+Pj7/BiOXw8fHBykpKbqGVa28vDwAgJ2dnUb51q1b4eDggG7duiEsLAxFRUUPbCciIgLW1tbi5OrqqrcYiajxMJAboJtDN3Rz6AYDuYHU4RCRRHTuQXrvvfcwcuRInDt3DmVlZfjwww9x7tw5HD58uNJ7kR7k5s2bUKlUcHJy0ih3cnLChQsXdA2rSmq1Gq+//joGDBiAbt26ieXPP/88HnvsMbi4uOD06dOYN28e0tPTsX379mrbCgsLQ2hoqLicn5/PJImIiKiZ0jlBGjhwIE6ePIkVK1age/fu+PHHH9G7d2+kpKSge/fu9RFjrc2aNQtnzpzBzz//rFE+Y8YMcb579+5wdnbGsGHDcPnyZbRv377KtkxMTGBiwu52ouauVFWKL85/AQB4ofMLfJM2UQtVqxcYtW/fHp9++mmdduzg4AADAwPk5ORolOfk5Og0Nqg6ISEh2LVrFw4ePIg2bdo8sK5SqQQAXLp0qdoEiYhahlJ1KdakrQEATOg4gQkSUQtVozFI+fn5NZ5qytjYGJ6enkhKShLL1Go1kpKSdBrLpE0QBISEhODbb7/FTz/9hLZt2z50m5MnTwIAnJ2da71fIiIiaj5q1INkY2NT4yfUVCpVjXceGhqKqVOnwsvLC3379kVkZCQKCwvFp9qmTJmC1q1bIyIiAkD5wO5z586J89euXcPJkydhYWEBd3d3AOW31WJjY/Hdd9/B0tIS2dnZAABra2uYmZnh8uXLiI2NxahRo2Bvb4/Tp09jzpw5GDRoEHr06FHj2ImIiKj5qlGCtH//fnH+ypUrePvttxEYGCj29KSkpGDLli1iIlNTEyZMwI0bNxAeHo7s7Gx4eHggISFBHLidmZkJufzfTq7r16+jV69e4vKqVauwatUqDB48GMnJyQCATz4p/wr3kCFDNPYVHR2NwMBAGBsbY9++fWIy5urqioCAACxYsECn2ImIiKj5qlGCNHjwYHF+6dKlWLNmDSZNmiSWjRkzBt27d8fGjRsxdepUnQIICQlBSEhIlevuJz33ubm5QRCEB7b3sPWurq46PW1HRERELY/O70FKSUmBl5dXpXIvL68HfiaEiIiIqKnQOUFydXWt8gm2TZs28b1ARERE1Czo/Jj/2rVrERAQgB9++EF8PD41NRUXL17EN998o/cAiYgakomBCTb7bhbniahl0rkHadSoUbh48SLGjBmD27dv4/bt23j66afx+++/Y9SoUfURIxFRgzGQG6CPog/6KPrwUyNELVitXhTZpk0bvPvuu/qOhYiIiKhRqFWCRETUXJWqS/H1718DAJ59/FkYyfkmbaKWiAkSEVEFpapSvHf0PQDA2PZjmSARtVA6j0EiIiIiau6YIBERERFp4S02IqLqHFwNyOvwa3JomP5iIaIGpfNP/q1btxAeHo79+/cjNzcXarVaY/3t27f1FhwRERGRFHROkF588UVcunQJQUFBcHJygkwmq4+4iIiIiCSjc4J06NAh/Pzzz+jZs2d9xENEREQkOZ0TpE6dOuHu3bv1EQsRkeSMDYyxftj68vnfD0kcDRFJReen2DZs2IB33nkHBw4cwK1bt5Cfn68xERE1ZYZyQwxqMwiD2gyCoYwP+hK1VDr3INnY2CA/Px9PPvmkRrkgCJDJZFCpVHoLjoiIiEgKOidIkydPhpGREWJjYzlIm4ianVJ1KXb/sRsAMFpQw4i9SEQtks4J0pkzZ3DixAl07NixPuIhIpJUqaoUC39ZCAAY0eY5JkhELZTOP/leXl64evVqfcRCRERE1Cjo3IP06quv4rXXXsObb76J7t27w8hI80OOPXr00FtwRERERFLQOUGaMGECAGD69OlimUwm4yBtIiIiajZ0TpAyMjLqIw4iIiKiRkPnBOmxxx6rjziIiIiIGo1af6b63LlzyMzMRElJiUb5mDFj6hwUERERkZR0TpD++OMPjBs3Dr/99ps49giA+D4kjkEioqbM2MAYqwavKp//45jE0RCRVHR+zP+1115D27ZtkZubC3Nzc5w9exYHDx6El5cXkpOT6yFEIqKGYyg3hK+bL3zdfPmpEaIWTOcepJSUFPz0009wcHCAXC6HXC7HwIEDERERgdmzZ+PEiRP1EScRERFRg9H5zyOVSgVLS0sAgIODA65fvw6gfPB2enq6fqMjImpgZeoy7L2yF3uv7EWZoJY6HCKSiM49SN26dcOpU6fQtm1bKJVKrFy5EsbGxti4cSPatWtXHzESETWYElUJ5h6YCwA42uY53mYjaqF0TpAWLFiAwsJCAMDSpUvx1FNP4YknnoC9vT3i4+P1HiARERFRQ9P5TyNfX18888wzAAB3d3dcuHABN2/eRG5uLp588kmdA1i/fj3c3NxgamoKpVKJ1NTUauuePXsWAQEBcHNzg0wmQ2RkZK3avHfvHmbNmgV7e3tYWFggICAAOTk5OsdOREREzZNOCVJpaSkMDQ1x5swZjXI7OzvxMX9dxMfHIzQ0FIsWLcLx48fRs2dP+Pr6Ijc3t8r6RUVFaNeuHVasWAGFQlHrNufMmYPvv/8e27Ztw4EDB3D9+nUx6SMiIiLSKUEyMjLCo48+qrd3Ha1ZswbBwcGYNm0aunTpgqioKJibm2Pz5s1V1u/Tpw8++OADTJw4ESYmJrVqMy8vD5999hnWrFmDJ598Ep6enoiOjsbhw4dx5MgRvRwXERERNW0632J75513MH/+fNy+fbtOOy4pKUFaWhp8fHz+DUYuh4+PD1JSUuqtzbS0NJSWlmrU6dSpEx599NEH7re4uBj5+fkaExERETVPOg/SXrduHS5dugQXFxc89thjaNWqlcb648eP16idmzdvQqVSwcnJSaPcyckJFy5c0DWsGreZnZ0NY2Nj2NjYVKqTnZ1dbdsRERFYsmRJreIiIiKipkXnBMnf378ewmj8wsLCEBoaKi7n5+fD1dVVwoiIqD4YGRhh2YBl5fOZZyWOhoikonOCtGjRIr3s2MHBAQYGBpWeHsvJyal2ALY+2lQoFCgpKcGdO3c0epEetl8TE5Nqxz0RUfNhJDeCv7t/+cLV85LGQkTSkewNaMbGxvD09ERSUpJYplarkZSUBG9v73pr09PTE0ZGRhp10tPTkZmZWev9EhERUfOicw+Sra1tlY/0y2QymJqawt3dHYGBgZg2bdpD2woNDcXUqVPh5eWFvn37IjIyEoWFheK2U6ZMQevWrREREQGgfBD2uXPnxPlr167h5MmTsLCwgLu7e43atLa2RlBQEEJDQ2FnZwcrKyu8+uqr8Pb2Rr9+/XQ9HUTUzJSpy3D4+mEAQH9BzTdpE7VQOidI4eHhePfddzFy5Ej07dsXAJCamoqEhATMmjULGRkZeOWVV1BWVobg4OAHtjVhwgTcuHED4eHhyM7OhoeHBxISEsRB1pmZmZDL//3ldP36dfTq1UtcXrVqFVatWoXBgwcjOTm5Rm0CwNq1ayGXyxEQEIDi4mL4+vpiw4YNup4KImqGSlQlmJU0CwA/NULUkumcIP38889Yvnw5Xn75ZY3y//73v/jxxx/xzTffoEePHvjoo48emiABQEhICEJCQqpcdz/puc/NzQ2CINSpTQAwNTXF+vXrsX79+oe2RURERC2Pzn8a7d27V+MdQvcNGzYMe/fuBQCMGjUKf/zxR92jIyIiIpKAzgmSnZ0dvv/++0rl33//Pezs7AAAhYWFsLS0rHt0RERERBLQ+RbbwoUL8corr2D//v3iGKRff/0Ve/bsQVRUFAAgMTERgwcP1m+kRERERA1E5wQpODgYXbp0wbp167B9+3YAQMeOHXHgwAH0798fAPDGG2/oN0oiIiKiBqRzggQAAwYMwIABA/QdCxEREVGjUKsESaVSYceOHTh/vvwts127dsWYMWNgYGCg1+CIiBqakYER5ivnl89f58MmRC2VzgnSpUuXMGrUKFy7dg0dO3YEUP4hV1dXV+zevRvt27fXe5BERA3FSG6ESZ0mlS9kRUgbDBFJRuen2GbPno327dvj6tWrOH78OI4fP47MzEy0bdsWs2fPro8YiYiIiBqUzj1IBw4cwJEjR8RH+gHA3t4eK1as4LgkImryVGoVjuceBwD0FtQw4Ju0iVoknRMkExMT/PPPP5XKCwoKYGxsrJegiIikUqwqxvS90wGUf2rEnAkSUYuk80/+U089hRkzZuDo0aMQBAGCIODIkSN4+eWXMWbMmPqIkYiIiKhB6ZwgffTRR2jfvj28vb1hamoKU1NTDBgwAO7u7oiMjKyHEImIiIgals632GxsbPDdd9/h0qVL4mP+nTt3hru7u96DIyIiIpKCzj1IS5cuRVFREdzd3fH000/j6aefhru7O+7evYulS5fWR4xEREREDUrnBGnJkiUoKCioVF5UVIQlS5boJSgiIiIiKemcIAmCAJlMVqn81KlTGo/+ExERETVVNR6DZGtrC5lMBplMhscff1wjSVKpVCgoKMDLL79cL0ESETUUI7kRQj1Dy+dzr0kcDRFJpcYJUmRkJARBwPTp07FkyRJYW1uL64yNjeHm5gZvb+96CZKIqKEYGRhhWrdp5Qv7+akRopaqxgnS1KlTAQBt27bFgAEDYGhYq+/cEhERETV6Oo9BsrS0FB/vB4DvvvsO/v7+mD9/PkpKSvQaHBFRQ1OpVThz8wzO3DwDlaCWOhwikojOCdJLL72E33//HQDwxx9/YMKECTA3N8e2bdvw1ltv6T1AIqKGVKwqxqTdkzBp9yQUM0EiarF0TpB+//13eHh4AAC2bduGwYMHIzY2FjExMfjmm2/0HR8RERFRg6vVY/5qdflfVfv27cOoUaMAAK6urrh586Z+oyMiIiKSgM4JkpeXF5YvX47//e9/OHDgAEaPHg0AyMjIgJOTk94DJCIiImpoOidIkZGROH78OEJCQvDOO++I32D7+uuv0b9/f70HSERERNTQdH5Wv0ePHvjtt98qlX/wwQcwMDDQS1BEREREUtLby4xMTU311RQRERGRpGqUINnZ2eH333+Hg4OD+MmR6ty+fVtvwRERNTQjuRFe6flK+Tx/nxG1WDVKkNauXQtLS0sA5WOQiIiaKyMDI8z0mFm+wE+NELVYNUqQ7n9mRHueiIiIqDmq8Rik/Pz8GtWzsrKqdTBERFJTC2r8cecPAEA7QYD8AUMKiKj5qvFj/jY2NrC1ta12ur++NtavXw83NzeYmppCqVQiNTX1gfW3bduGTp06wdTUFN27d8eePXs01stksiqnDz74QKzj5uZWaf2KFStqFT8RNR/3yu5h3M5xGLdzHO4JKqnDISKJ1LgHaf/+/eK8IAgYNWoUNm3ahNatW9cpgPj4eISGhiIqKgpKpRKRkZHw9fVFeno6HB0dK9U/fPgwJk2ahIiICDz11FOIjY2Fv78/jh8/jm7dugEAsrKyNLb54YcfEBQUhICAAI3ypUuXIjg4WFy+P86KqMXj2BsiauFqnCANHjxYY9nAwAD9+vVDu3bt6hTAmjVrEBwcjGnTpgEAoqKisHv3bmzevBlvv/12pfoffvgh/Pz88OabbwIAli1bhsTERKxbtw5RUVEAAIVCobHNd999h6FDh1aK1dLSslJdIiIiIp3fpK1PJSUlSEtLg4+Pj1gml8vh4+ODlJSUKrdJSUnRqA8Avr6+1dbPycnB7t27ERQUVGndihUrYG9vj169euGDDz5AWVlZtbEWFxcjPz9fYyIiIqLmSW8viqyNmzdvQqVSVfqGm5OTEy5cuFDlNtnZ2VXWz87OrrL+li1bYGlpiWeeeUajfPbs2ejduzfs7Oxw+PBhhIWFISsrC2vWrKmynYiICCxZsqSmh0ZERERNWJ0SpAe9MLKx2Lx5MyZPnlzpTd+hoaHifI8ePWBsbIyXXnoJERERMDExqdROWFiYxjb5+flwdXWtv8CJiIhIMjVOkLR7YO7du4eXX34ZrVq10ijfvn17jXfu4OAAAwMD5OTkaJTn5ORUOzZIoVDUuP6hQ4eQnp6O+Pj4h8aiVCpRVlaGK1euoGPHjpXWm5iYVJk4ERERUfNT4zFI1tbWGtMLL7wAFxeXSuW6MDY2hqenJ5KSksQytVqNpKQkeHt7V7mNt7e3Rn0ASExMrLL+Z599Bk9PT/Ts2fOhsZw8eRJyubzKJ+eIqOUwkhshsGsgArsGwqgJ9JITUf2ocQ9SdHR0vQQQGhqKqVOnwsvLC3379kVkZCQKCwvFp9qmTJmC1q1bIyKi/LHj1157DYMHD8bq1asxevRoxMXF4dixY9i4caNGu/n5+di2bRtWr15daZ8pKSk4evQohg4dCktLS6SkpGDOnDl44YUXav0uJyJqHowMjPCG1xvlC3zdAVGLJekgbQCYMGECbty4gfDwcGRnZ8PDwwMJCQniQOzMzEzI5f92dPXv3x+xsbFYsGAB5s+fjw4dOmDHjh3iO5Dui4uLgyAImDRpUqV9mpiYIC4uDosXL0ZxcTHatm2LOXPmaIwxIiIiopZL8gQJAEJCQhASElLluuTk5Epl48ePx/jx4x/Y5owZMzBjxowq1/Xu3RtHjhzROU4iav7UghpZheUvm3Xmp0aIWqxGkSARETUW98ruwe8bPwDA0TbPwVzGX5NELZGkL4okIiIiaoyYIBERERFpYYJEREREpIUJEhEREZEWJkhEREREWpggEREREWnh86tERBUYyg0xoeOE8vkiiYMhIskwQSIiqsDYwBgL+i0oX+CnRohaLN5iIyIiItLCHiQiogoEQcDfxX8DAGwFATJ+aoSoRWKCRJWsTfy9wfY1Z/jjDbavhpbyx61abXekrG7nv19m7farC+929vW+D6ncLbuLwfGDAfBTI0QtGW+xEREREWlhgkRERESkhQkSERERkRYmSERERERamCARERERaWGCRERERKSFz68SEVVgKDfEmPZjyudL+A4kopaKCRIRUQXGBsZ4d+C75Qv81AhRi8VbbERERERa2INERFSBIAi4W3YXAGDGT40QtVhMkIiIKrhbdhfKWCUAfmqEqCXjLTYiIiIiLUyQiIiIiLQwQSIiIiLSwgSJiIiISAsTJCIiIiItTJCIiIiItPD5VSKiCgzkBhj+2PDyeTXfgUTUUjFBIiKqwMTABGuGrClf4KdGiFqsRnGLbf369XBzc4OpqSmUSiVSU1MfWH/btm3o1KkTTE1N0b17d+zZs0djfWBgIGQymcbk5+enUef27duYPHkyrKysYGNjg6CgIBQUFOj92IiIiKjpkTxBio+PR2hoKBYtWoTjx4+jZ8+e8PX1RW5ubpX1Dx8+jEmTJiEoKAgnTpyAv78//P39cebMGY16fn5+yMrKEqcvv/xSY/3kyZNx9uxZJCYmYteuXTh48CBmzJhRb8dJRERETYfkCdKaNWsQHByMadOmoUuXLoiKioK5uTk2b95cZf0PP/wQfn5+ePPNN9G5c2csW7YMvXv3xrp16zTqmZiYQKFQiJOtra247vz580hISMCmTZugVCoxcOBAfPzxx4iLi8P169er3G9xcTHy8/M1JiJqfopKi9B9S3d039IdReoyqcMhIolImiCVlJQgLS0NPj4+YplcLoePjw9SUlKq3CYlJUWjPgD4+vpWqp+cnAxHR0d07NgRr7zyCm7duqXRho2NDby8vMQyHx8fyOVyHD16tMr9RkREwNraWpxcXV11Pl4iIiJqGiRNkG7evAmVSgUnJyeNcicnJ2RnZ1e5TXZ29kPr+/n54fPPP0dSUhLef/99HDhwACNHjoRKpRLbcHR01GjD0NAQdnZ21e43LCwMeXl54nT16lWdj5eIiIiahmb5FNvEiRPF+e7du6NHjx5o3749kpOTMWzYsFq1aWJiAhMTE32FSERERI2YpD1IDg4OMDAwQE5OjkZ5Tk4OFApFldsoFAqd6gNAu3bt4ODggEuXLoltaA8CLysrw+3btx/YDhEREbUMkvYgGRsbw9PTE0lJSfD39wcAqNVqJCUlISQkpMptvL29kZSUhNdff10sS0xMhLe3d7X7+euvv3Dr1i04OzuLbdy5cwdpaWnw9PQEAPz0009Qq9VQKpX6OTgiqpu6voNoaJh+4iCiFknyp9hCQ0Px6aefYsuWLTh//jxeeeUVFBYWYtq0aQCAKVOmICzs3190r732GhISErB69WpcuHABixcvxrFjx8SEqqCgAG+++SaOHDmCK1euICkpCWPHjoW7uzt8fX0BAJ07d4afnx+Cg4ORmpqKX375BSEhIZg4cSJcXFwa/iQQERFRoyL5GKQJEybgxo0bCA8PR3Z2Njw8PJCQkCAOxM7MzIRc/m8e179/f8TGxmLBggWYP38+OnTogB07dqBbt24AAAMDA5w+fRpbtmzBnTt34OLighEjRmDZsmUaY4i2bt2KkJAQDBs2DHK5HAEBAfjoo48a9uCJqNExkBvgidZPlM/L+KkRopZK8gQJAEJCQqq9pZacnFypbPz48Rg/fnyV9c3MzLB3796H7tPOzg6xsbE6xUlEzZ+JgQk2+GwoX+CnRohaLMlvsRERERE1NkyQiIiIiLQ0iltsRESNRVFpEYZ8NQQAkKx4GuZy/pokaon4k09EpOVu2V2pQyAiifEWGxEREZEWJkhEREREWpggEREREWlhgkRERESkhQkSERERkRY+xUZEVIFcJoeXk1f5vMSxEJF0mCAREVVgamiKaL/o8gV+aoSoxeIfSERERERa2INERFRfHtID1S/zVgMFQkS6YoJERFRBUWkR/L7xAwAkPDKcnxohaqH4k09EpOXv4r+lDoGIJMYxSERERERamCARERERaWGCRERERKSFCRIRERGRFiZIRERERFr4FBsRUQVymRxd7buWz0scCxFJhwkSEVEFpoamiHsqrnyBnxoharH4BxIRERGRFiZIRERERFp4i42IqIK7ZXfhv8MfALDDdiDM+KkRohaJP/lERBUIgoDrhdfL520lDoaIJMNbbERERERamCARERERaWGCRERERKSFCRIRERGRlkaRIK1fvx5ubm4wNTWFUqlEamrqA+tv27YNnTp1gqmpKbp37449e/aI60pLSzFv3jx0794drVq1gouLC6ZMmYLr169rtOHm5gaZTKYxrVixol6Oj4iIiJoWyROk+Ph4hIaGYtGiRTh+/Dh69uwJX19f5ObmVln/8OHDmDRpEoKCgnDixAn4+/vD398fZ86cAQAUFRXh+PHjWLhwIY4fP47t27cjPT0dY8aMqdTW0qVLkZWVJU6vvvpqvR4rETV+MpkM7a3bo711e8ikDoaIJCP5Y/5r1qxBcHAwpk2bBgCIiorC7t27sXnzZrz99tuV6n/44Yfw8/PDm2++CQBYtmwZEhMTsW7dOkRFRcHa2hqJiYka26xbtw59+/ZFZmYmHn30UbHc0tISCoWiHo+OiJoaM0Mz7PDfUb7AT40QtViS9iCVlJQgLS0NPj4+YplcLoePjw9SUlKq3CYlJUWjPgD4+vpWWx8A8vLyIJPJYGNjo1G+YsUK2Nvbo1evXvjggw9QVlZWbRvFxcXIz8/XmIiIiKh5krQH6ebNm1CpVHByctIod3JywoULF6rcJjs7u8r62dnZVda/d+8e5s2bh0mTJsHKykosnz17Nnr37g07OzscPnwYYWFhyMrKwpo1a6psJyIiAkuWLNHl8IiIiKiJkvwWW30qLS3Fc889B0EQ8Mknn2isCw0NFed79OgBY2NjvPTSS4iIiICJiUmltsLCwjS2yc/Ph6ura/0FT0SSuFt2F5N2TQIAfGnpyU+NELVQkv7kOzg4wMDAADk5ORrlOTk51Y4NUigUNap/Pzn6888/8dNPP2n0HlVFqVSirKwMV65cQceOHSutNzExqTJxIqLmRRAEXM67XD5v6SlxNEQkFUnHIBkbG8PT0xNJSUlimVqtRlJSEry9vavcxtvbW6M+ACQmJmrUv58cXbx4Efv27YO9vf1DYzl58iTkcjkcHR1reTRERETUXEjedxwaGoqpU6fCy8sLffv2RWRkJAoLC8Wn2qZMmYLWrVsjIqL8aZLXXnsNgwcPxurVqzF69GjExcXh2LFj2LhxI4Dy5OjZZ5/F8ePHsWvXLqhUKnF8kp2dHYyNjZGSkoKjR49i6NChsLS0REpKCubMmYMXXngBtrb8OiUREVFLJ3mCNGHCBNy4cQPh4eHIzs6Gh4cHEhISxIHYmZmZkMv/7ejq378/YmNjsWDBAsyfPx8dOnTAjh070K1bNwDAtWvXsHPnTgCAh4eHxr7279+PIUOGwMTEBHFxcVi8eDGKi4vRtm1bzJkzR2OMEREREbVckidIABASEoKQkJAq1yUnJ1cqGz9+PMaPH19lfTc3NwiC8MD99e7dG0eOHNE5TiIiImoZJH+TNhEREVFj0yh6kIiIGguZTAaXVi7l8xLHQkTSYYJERFSBmaEZ9j67t3yBnxoharF4i42IiIhICxMkIiIiIi28xUZEVMG9snsITAgEAMSYdYUpPzVC1CLxJ5+IqAK1oMbZW2fL59t0lTgaIpIKb7ERERERaWGCRERERKSFCRIRERGRFiZIRERERFqYIBERERFp4VNsRERabE1spQ6BiCTGBImIqAJzI3McnHiwfIGfGiFqsXiLjYiIiEgLEyQiIiIiLbzFRkRUwb2ye3hl3ysAgE+M2vFTI0QtFH/yiYgqUAtqHMs5Vj7fpp3E0RCRVHiLjYiIiEgLe5CIqHmq7RNo6jL9xkFETRJ7kIiIiIi0MEEiIiIi0sIEiYiIiEgLxyAREWkxkxlIHQIRSYwJEhFRBeZyQ6S6TpA6DCKSGG+xEREREWlhgkRERESkhQkSEVEFxYIKM3OTMTM3GcWCSupwiEgiHINERFSBShBw6N51cR4yiQMiIkmwB4mIiIhICxMkIiIiIi2NIkFav3493NzcYGpqCqVSidTU1AfW37ZtGzp16gRTU1N0794de/bs0VgvCALCw8Ph7OwMMzMz+Pj44OLFixp1bt++jcmTJ8PKygo2NjYICgpCQUGB3o+NiIiImh7JE6T4+HiEhoZi0aJFOH78OHr27AlfX1/k5uZWWf/w4cOYNGkSgoKCcOLECfj7+8Pf3x9nzpwR66xcuRIfffQRoqKicPToUbRq1Qq+vr64d++eWGfy5Mk4e/YsEhMTsWvXLhw8eBAzZsyo9+MlIiKixk/yBGnNmjUIDg7GtGnT0KVLF0RFRcHc3BybN2+usv6HH34IPz8/vPnmm+jcuTOWLVuG3r17Y926dQDKe48iIyOxYMECjB07Fj169MDnn3+O69evY8eOHQCA8+fPIyEhAZs2bYJSqcTAgQPx8ccfIy4uDtevX2+oQyciIqJGStKn2EpKSpCWloawsDCxTC6Xw8fHBykpKVVuk5KSgtDQUI0yX19fMfnJyMhAdnY2fHx8xPXW1tZQKpVISUnBxIkTkZKSAhsbG3h5eYl1fHx8IJfLcfToUYwbN67SfouLi1FcXCwu5+XlAQDy8/N1P/BauFfYPG//NdT5u6++z2Ph3eKHV3qIusaojxgeJr/w3sMrNVFF6jKo7pY/3p9feA9l8vr7NdkQ/6/uX08N/bNG1Fjd/1kQBOGB9SRNkG7evAmVSgUnJyeNcicnJ1y4cKHKbbKzs6usn52dLa6/X/agOo6OjhrrDQ0NYWdnJ9bRFhERgSVLllQqd3V1re7wqAbmSx1Ao7RO6gDo/znjvNQh6EH59cSfNSJN//zzD6ytratdz/cg1VBYWJhGz5Varcbt27dhb28Pmaz+XpSSn58PV1dXXL16FVZWVvW2n+aI5672eO5qj+eu9njuao/nruYEQcA///wDFxeXB9aTNEFycHCAgYEBcnJyNMpzcnKgUCiq3EahUDyw/v3/5uTkwNnZWaOOh4eHWEd7EHhZWRlu375d7X5NTExgYmKiUWZjY/PgA9QjKysrXvS1xHNXezx3tcdzV3s8d7XHc1czD+o5uk/SQdrGxsbw9PREUlKSWKZWq5GUlARvb+8qt/H29taoDwCJiYli/bZt20KhUGjUyc/Px9GjR8U63t7euHPnDtLS0sQ6P/30E9RqNZRKpd6Oj4iIiJomyW+xhYaGYurUqfDy8kLfvn0RGRmJwsJCTJs2DQAwZcoUtG7dGhEREQCA1157DYMHD8bq1asxevRoxMXF4dixY9i4cSMAQCaT4fXXX8fy5cvRoUMHtG3bFgsXLoSLiwv8/f0BAJ07d4afnx+Cg4MRFRWF0tJShISEYOLEiQ/tciMiIqLmT/IEacKECbhx4wbCw8ORnZ0NDw8PJCQkiIOsMzMzIZf/29HVv39/xMbGYsGCBZg/fz46dOiAHTt2oFu3bmKdt956C4WFhZgxYwbu3LmDgQMHIiEhAaampmKdrVu3IiQkBMOGDYNcLkdAQAA++uijhjvwGjIxMcGiRYsq3d6jh+O5qz2eu9rjuas9nrva47nTP5nwsOfciIiIiFoYyV8USURERNTYMEEiIiIi0sIEiYiIiEgLEyQiIiIiLUyQGrExY8bg0UcfhampKZydnfHiiy9qfEz3ypUrkMlklaYjR45IGHXj8LBzBwCnT5/GE088AVNTU7i6umLlypUSRdt4XLlyBUFBQWjbti3MzMzQvn17LFq0CCUlJRp1eN1VVpNzB/C6q8q7776L/v37w9zcvNoX8FZ1zcXFxTVsoI1QTc5dZmYmRo8eDXNzczg6OuLNN99EWVlZwwbaBEn+mD9Vb+jQoZg/fz6cnZ1x7do1zJ07F88++ywOHz6sUW/fvn3o2rWruGxvb9/QoTY6Dzt3+fn5GDFiBHx8fBAVFYXffvsN06dPh42NDWbMmCFx9NK5cOEC1Go1/vvf/8Ld3R1nzpxBcHAwCgsLsWrVKo26vO401eTc8bqrWklJCcaPHw9vb2989tln1daLjo6Gn5+fuNyQXzNorB527lQqFUaPHg2FQoHDhw8jKysLU6ZMgZGREd577z0JIm5CBGoyvvvuO0EmkwklJSWCIAhCRkaGAEA4ceKEtIE1AdrnbsOGDYKtra1QXFws1pk3b57QsWNHqUJstFauXCm0bdtWXOZ1V3Pa547X3YNFR0cL1tbWVa4DIHz77bcNGk9TUt2527NnjyCXy4Xs7Gyx7JNPPhGsrKw0rkOqjLfYmojbt29j69at6N+/P4yMjDTWjRkzBo6Ojhg4cCB27twpUYSNV1XnLiUlBYMGDYKxsbFYz9fXF+np6fj777+lCrVRysvLg52dXaVyXncPp33ueN3VzaxZs+Dg4IC+ffti8+bNEPgav4dKSUlB9+7dxZcvA+XXXH5+Ps6ePSthZI0fE6RGbt68eWjVqhXs7e2RmZmJ7777TlxnYWGB1atXY9u2bdi9ezcGDhwIf39//mP1/x507rKzszV+YQAQl7Ozsxs0zsbs0qVL+Pjjj/HSSy+JZbzuaqaqc8frrvaWLl2Kr776ComJiQgICMDMmTPx8ccfSx1Wo8drrg6k7sJqaebNmycAeOB0/vx5sf6NGzeE9PR04ccffxQGDBggjBo1SlCr1dW2/+KLLwoDBw5siENpcPo8d8OHDxdmzJih0f7Zs2cFAMK5c+ca9Lgagq7nThAE4a+//hLat28vBAUFPbR9Xnc1O3ct6bqrzXl70C02bQsXLhTatGlTD5FLT5/nLjg4WBgxYoRGWWFhoQBA2LNnT30eRpPHQdoN7I033kBgYOAD67Rr106cd3BwgIODAx5//HF07twZrq6uOHLkCLy9vavcVqlUIjExUZ8hNxr6PHcKhQI5OTka295fVigUeo9darqeu+vXr2Po0KHo37+/+CHoB+F1V7Nz15KuO13Pm66USiWWLVuG4uLiZvf9MX2eO4VCgdTUVI2y5nrN6RsTpAb2yCOP4JFHHqnVtmq1GgBQXFxcbZ2TJ0/C2dm5Vu03dvo8d97e3njnnXdQWloqjktKTExEx44dYWtrq5+AGxFdzt21a9cwdOhQeHp6Ijo6WuNj0dXhdVfuYeeuJV13dfl5rYmTJ0/C1ta22SVHgH7Pnbe3N959913k5ubC0dERQPk1Z2VlhS5duuhlH80VE6RG6ujRo/j1118xcOBA2Nra4vLly1i4cCHat28v9h5t2bIFxsbG6NWrFwBg+/bt2Lx5MzZt2iRl6JKrybl7/vnnsWTJEgQFBWHevHk4c+YMPvzwQ6xdu1bi6KV17do1DBkyBI899hhWrVqFGzduiOvu/7XJ665qNTl3vO6qlpmZidu3byMzMxMqlQonT54EALi7u8PCwgLff/89cnJy0K9fP5iamiIxMRHvvfce5s6dK23gjcDDzt2IESPQpUsXvPjii1i5ciWys7OxYMECzJo1q1kml3ol9T0+qtrp06eFoUOHCnZ2doKJiYng5uYmvPzyy8Jff/0l1omJiRE6d+4smJubC1ZWVkLfvn2Fbdu2SRh141CTcycIgnDq1Clh4MCBgomJidC6dWthxYoVEkXceERHR1c75uE+XndVq8m5EwRed1WZOnVqledt//79giAIwg8//CB4eHgIFhYWQqtWrYSePXsKUVFRgkqlkjbwRuBh504QBOHKlSvCyJEjBTMzM8HBwUF44403hNLSUumCbiJkgsDnJImIiIgq4mP+RERERFqYIBERERFpYYJEREREpIUJEhEREZEWJkhEREREWpggEREREWlhgkRERESkhQkSERERkRYmSETU4GQyGXbs2CF1GA0uMDAQ/v7+4vKQIUPw+uuvSxYPEVWPCRIR6dWNGzfwyiuv4NFHH4WJiQkUCgV8fX3xyy+/SB1ao7N9+3YsW7ZMr21qJ2FEVDv8WC0R6VVAQABKSkqwZcsWtGvXDjk5OUhKSsKtW7ekDk0vSktLYWRkpJe27Ozs9NIOEekfe5CISG/u3LmDQ4cO4f3338fQoUPx2GOPoW/fvggLC8OYMWOq3e63337Dk08+CTMzM9jb22PGjBkoKCgQ19/vFVmyZAkeeeQRWFlZ4eWXX0ZJSYlYR61WIyIiAm3btoWZmRl69uyJr7/++oHxZmVlYfTo0TAzM0Pbtm0RGxsLNzc3REZGinVkMhk++eQTjBkzBq1atcK7774LlUqFoKAgcV8dO3bEhx9+qNG2SqVCaGgobGxsYG9vj7feegvan77UvsVWXFyMuXPnonXr1mjVqhWUSiWSk5PF9TExMbCxscHevXvRuXNnWFhYwM/PD1lZWQCAxYsXY8uWLfjuu+8gk8kgk8k0tieimmOCRER6Y2FhAQsLC+zYsQPFxcU12qawsBC+vr6wtbXFr7/+im3btmHfvn0ICQnRqJeUlITz588jOTkZX375JbZv344lS5aI6yMiIvD5558jKioKZ8+exZw5c/DCCy/gwIED1e57ypQpuH79OpKTk/HNN99g48aNyM3NrVRv8eLFGDduHH777TdMnz4darUabdq0wbZt23Du3DmEh4dj/vz5+Oqrr8RtVq9ejZiYGGzevBk///wzbt++jW+//faB5yIkJAQpKSmIi4vD6dOnMX78ePj5+eHixYtinaKiIqxatQr/+9//cPDgQWRmZmLu3LkAgLlz5+K5554Tk6asrCz079//wf8DiKhqAhGRHn399deCra2tYGpqKvTv318ICwsTTp06pVEHgPDtt98KgiAIGzduFGxtbYWCggJx/e7duwW5XC5kZ2cLgiAIU6dOFezs7ITCwkKxzieffCJYWFgIKpVKuHfvnmBubi4cPnxYYz9BQUHCpEmTqozz/PnzAgDh119/FcsuXrwoABDWrl2rEevrr7/+0OOeNWuWEBAQIC47OzsLK1euFJdLS0uFNm3aCGPHjhXLBg8eLLz22muCIAjCn3/+KRgYGAjXrl3TaHfYsGFCWFiYIAiCEB0dLQAQLl26JK5fv3694OTkJC5PnTpVYx9EVDscg0REehUQEIDRo0fj0KFDOHLkCH744QesXLkSmzZtQmBgYKX658+fR8+ePdGqVSuxbMCAAVCr1UhPT4eTkxMAoGfPnjA3NxfreHt7o6CgAFevXkVBQQGKioowfPhwjbZLSkrQq1evKuNMT0+HoaEhevfuLZa5u7vD1ta2Ul0vL69KZevXr8fmzZuRmZmJu3fvoqSkBB4eHgCAvLw8ZGVlQalUivUNDQ3h5eVV6Tbbfb/99htUKhUef/xxjfLi4mLY29uLy+bm5mjfvr247OzsXGWvFxHVDRMkItI7U1NTDB8+HMOHD8fChQvxn//8B4sWLaoyQdKH++OVdu/ejdatW2usMzExqXP7FZM3AIiLi8PcuXOxevVqeHt7w9LSEh988AGOHj1a630UFBTAwMAAaWlpMDAw0FhnYWEhzmsPEJfJZNUmXURUe0yQiKjedenSpdr3HnXu3BkxMTEoLCwUE5FffvkFcrkcHTt2FOudOnUKd+/ehZmZGQDgyJEjsLCwgKurK+zs7GBiYoLMzEwMHjy4RjF17NgRZWVlOHHiBDw9PQEAly5dwt9///3QbX/55Rf0798fM2fOFMsuX74szltbW8PZ2RlHjx7FoEGDAABlZWVIS0vT6LGqqFevXlCpVMjNzcUTTzxRo2OoirGxMVQqVa23J6JyHKRNRHpz69YtPPnkk/jiiy9w+vRpZGRkYNu2bVi5ciXGjh1b5TaTJ0+Gqakppk6dijNnzmD//v149dVX8eKLL4q314Dy22VBQUE4d+4c9uzZg0WLFiEkJARyuRyWlpaYO3cu5syZgy1btuDy5cs4fvw4Pv74Y2zZsqXK/Xbq1Ak+Pj6YMWMGUlNTceLECcyYMQNmZmaQyWQPPM4OHTrg2LFj2Lt3L37//XcsXLgQv/76q0ad1157DStWrMCOHTtw4cIFzJw5E3fu3Km2zccffxyTJ0/GlClTsH37dmRkZCA1NRURERHYvXv3A+OpyM3NDadPn0Z6ejpu3ryJ0tLSGm9LRP9iDxIR6Y2FhQWUSiXWrl2Ly5cvo7S0FK6urggODsb8+fOr3Mbc3Bx79+7Fa6+9hj59+sDc3BwBAQFYs2aNRr1hw4ahQ4cOGDRoEIqLizFp0iQsXrxYXL9s2TI88sgjiIiIwB9//AEbGxv07t272v0CwOeff46goCAMGjQICoUCEREROHv2LExNTR94nC+99BJOnDiBCRMmQCaTYdKkSZg5cyZ++OEHsc4bb7yBrKwsTJ06FXK5HNOnT8e4ceOQl5dXbbvR0dFYvnw53njjDVy7dg0ODg7o168fnnrqqQfGU1FwcDCSk5Ph5eWFgoIC7N+/H0OGDKnx9kRUTibw5jURNXKBgYG4c+dOvX+e5K+//oKrqyv27duHYcOG1eu+iKhxYw8SEbVYP/30EwoKCtC9e3dkZWXhrbfegpubmzhuiIhaLiZIRNRilZaWYv78+fjjjz9gaWmJ/v37Y+vWrXr7lAgRNV28xUZERESkhU+xEREREWlhgkRERESkhQkSERERkRYmSERERERamCARERERaWGCRERERKSFCRIRERGRFiZIRERERFr+D6a5zjiV0vvNAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(G_jtfs[:, 0].detach(), density=True, alpha=0.5, label='JTFS')\n",
    "plt.hist(G_scrapl[:, 0].detach(), density=True,  alpha=0.5, label='SCRAPL')\n",
    "\n",
    "y_coords = (4, 3)\n",
    "Sx_best = load_jtfs(*x_coords, seed=0, data_dir=data_dir)[:, 0]\n",
    "Jx_best = load_jacobian(*x_coords, seed=0, data_dir=data_dir)[:, 0]\n",
    "G_best = (Jx_best.T @ (Sy - Sx_best)).detach()\n",
    "plt.plot([G_best.item()]*2, [0, plt.gca().get_ylim()[1]], '--', label='Best')\n",
    "\n",
    "half_width = torch.max(torch.abs(G_jtfs - G_best)).detach().numpy()\n",
    "plt.xlim(G_best - 0.75*half_width, G_best + 0.75*half_width)\n",
    "plt.legend()\n",
    "plt.xlabel(\"Slope gradient\")\n",
    "plt.ylabel(\"Histogram density\")\n",
    "\n",
    "# Caution: does not work on every path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "scrapl",
   "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.11"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}