S1+S2 euclidean.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2023-07-27 22:37:09.800765: I tensorflow/core/platform/cpu_feature_guard.cc:193] This TensorFlow binary is optimized with oneAPI Deep Neural Network Library (oneDNN) to use the following CPU instructions in performance-critical operations:  AVX2 AVX512F AVX512_VNNI FMA\n",
      "To enable them in other operations, rebuild TensorFlow with the appropriate compiler flags.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CQT kernels created, time used = 0.0196 seconds\n",
      "CQT kernels created, time used = 0.0208 seconds\n",
      "CQT kernels created, time used = 0.0177 seconds\n"
     ]
    }
   ],
   "source": [
    "import auraloss\n",
    "from kymatio.torch import TimeFrequencyScattering\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import sys\n",
    "import torch\n",
    "import torchnmf\n",
    "import tqdm\n",
    "\n",
    "import cnn\n",
    "import synth\n",
    "\n",
    "sys.path.append('../dcase2023_task7_eval_fad')\n",
    "import audioset_model\n",
    "\n",
    "# duration = 4\n",
    "# event_duration = 2**(-4)\n",
    "# sr = 2**13\n",
    "# fmin = 2**8\n",
    "# fmax = 2**11\n",
    "# Q = 24\n",
    "# hop_length = 2**6\n",
    "# n_events = 2**6\n",
    "\n",
    "n_densities = 7\n",
    "n_slopes = 7\n",
    "n_folds = n_densities * n_slopes\n",
    "batch_size = 1\n",
    "\n",
    "dataset = cnn.ChirpTextureDataModule(\n",
    "    n_densities=n_densities,\n",
    "    n_slopes=n_slopes,\n",
    "    n_folds=n_folds,\n",
    "    batch_size=batch_size)\n",
    "dataset.setup()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2023-07-27 22:37:19.031971: I tensorflow/core/platform/cpu_feature_guard.cc:193] This TensorFlow binary is optimized with oneAPI Deep Neural Network Library (oneDNN) to use the following CPU instructions in performance-critical operations:  AVX2 AVX512F AVX512_VNNI FMA\n",
      "To enable them in other operations, rebuild TensorFlow with the appropriate compiler flags.\n",
      "/Users/user/miniconda3/envs/scrapl/lib/python3.10/site-packages/tensorflow/python/keras/engine/base_layer_v1.py:1694: UserWarning: `layer.apply` is deprecated and will be removed in a future version. Please use `layer.__call__` method instead.\n",
      "  warnings.warn('`layer.apply` is deprecated and '\n",
      "/Users/user/miniconda3/envs/scrapl/lib/python3.10/site-packages/tensorflow/python/keras/legacy_tf_layers/core.py:332: UserWarning: `tf.layers.flatten` is deprecated and will be removed in a future version. Please use `tf.keras.layers.Flatten` instead.\n",
      "  warnings.warn('`tf.layers.flatten` is deprecated and '\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO:tensorflow:Restoring parameters from ../dcase2023_task7_eval_fad/data/vggish_model.ckpt\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2023-07-27 22:37:19.417425: I tensorflow/compiler/mlir/mlir_graph_optimization_pass.cc:357] MLIR V1 optimization pass is not enabled\n"
     ]
    }
   ],
   "source": [
    "mss_distance = auraloss.freq.MultiResolutionSTFTLoss()\n",
    "\n",
    "def model_distance(x, y, model):\n",
    "    x_np = x[384:-384].detach().numpy()\n",
    "    y_np = y[384:-384].detach().numpy()\n",
    "    x_logmel = model.extract_features(x_np)\n",
    "    y_logmel = model.extract_features(y_np)\n",
    "    x_logmel = np.stack(x_logmel, axis=0).mean(axis=0)\n",
    "    y_logmel = np.stack(y_logmel, axis=0).mean(axis=0)\n",
    "    x_vgg = model.process_batch(x_logmel)\n",
    "    y_vgg = model.process_batch(y_logmel)\n",
    "    return np.linalg.norm(x_vgg - y_vgg)\n",
    "\n",
    "checkpoint = '../dcase2023_task7_eval_fad/data/vggish_model.ckpt'\n",
    "vggish = audioset_model.AudioSetModel(checkpoint)\n",
    "\n",
    "sc = TimeFrequencyScattering(\n",
    "    shape=(2**15),\n",
    "    J=6,\n",
    "    Q=(24, 2),\n",
    "    Q_fr=2,\n",
    "    J_fr=5,\n",
    "    T='global',\n",
    "    F='global',\n",
    "    format='time',\n",
    ")\n",
    "\n",
    "y = synth.generate_chirp_texture(\n",
    "    theta_density=torch.tensor(0.5),\n",
    "    theta_slope=torch.tensor(0.5),\n",
    "    duration=dataset.train_ds.duration,\n",
    "    event_duration=dataset.train_ds.event_duration,\n",
    "    sr=dataset.train_ds.sr,\n",
    "    fmin=dataset.train_ds.fmin,\n",
    "    fmax=dataset.train_ds.fmax,\n",
    "    n_events=dataset.train_ds.n_events,\n",
    "    Q=dataset.train_ds.Q,\n",
    "    hop_length=dataset.train_ds.hop_length,\n",
    "    seed=0,\n",
    ")\n",
    "\n",
    "Sy = sc(y)\n",
    "\n",
    "def scattering_audio_distance(x, y, Sy, sc):\n",
    "    # Skip order 0, which is not guaranteed to be nonnegative\n",
    "    Sx = sc(x)[1:]\n",
    "    return torchnmf.metrics.kl_div(Sx/Sx.sum(), Sy[1:]/Sy[1:].sum()), Sx\n",
    "\n",
    "distances = {\n",
    "    'mss': lambda x,y: mss_distance(\n",
    "        x.unsqueeze(0).unsqueeze(0), y.unsqueeze(0).unsqueeze(0)),\n",
    "    'vgg': lambda x, y: model_distance(x, y, vggish),\n",
    "    'jtfs': lambda x, y: scattering_audio_distance(x, y, Sy, sc),\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "49it [00:01, 30.89it/s]\n",
      "49it [00:01, 33.63it/s]\n",
      "49it [00:03, 12.35it/s]\n",
      "49it [00:04, 11.70it/s]\n",
      "49it [14:28, 17.73s/it]\n",
      "49it [14:15, 17.47s/it]\n"
     ]
    }
   ],
   "source": [
    "out = {}\n",
    "SX = {}\n",
    "\n",
    "for distance_key, distance in distances.items():\n",
    "\n",
    "    # 0: matched seed between x and y\n",
    "    # 1: difference seed between x and y\n",
    "    for seed in [0, 1]:\n",
    "\n",
    "        X = {}\n",
    "        sp_df = pd.DataFrame(columns=[\n",
    "            \"density_idx\", \"slope_idx\", \"distance\", \"density_grad\", \"slope_grad\"])\n",
    "        for _, row in tqdm.tqdm(dataset.df.iterrows()):\n",
    "            theta_density = torch.tensor(row[\"density\"], requires_grad=False)\n",
    "            theta_slope = torch.tensor(row[\"slope\"], requires_grad=False)\n",
    "            x = synth.generate_chirp_texture(\n",
    "                theta_density=theta_density,\n",
    "                theta_slope=theta_slope,\n",
    "                duration=dataset.train_ds.duration,\n",
    "                event_duration=dataset.train_ds.event_duration,\n",
    "                sr=dataset.train_ds.sr,\n",
    "                fmin=dataset.train_ds.fmin,\n",
    "                fmax=dataset.train_ds.fmax,\n",
    "                n_events=dataset.train_ds.n_events,\n",
    "                Q=dataset.train_ds.Q,\n",
    "                hop_length=dataset.train_ds.hop_length,\n",
    "                seed=seed,\n",
    "            )\n",
    "            if distance_key == 'jtfs':\n",
    "                xy_dist, Sx = distance(x, y)\n",
    "                SX[(distance_key, seed, row[\"density\"], row[\"slope\"])] = Sx\n",
    "            else:\n",
    "                xy_dist = distance(x, y)\n",
    "            # if hasattr(xy_dist, 'backward'):\n",
    "            #     xy_dist.backward()\n",
    "            # if theta_density.grad is None:\n",
    "            #     theta_density.grad = torch.zeros_like(theta_density)\n",
    "            # if theta_slope.grad is None:\n",
    "            #     theta_slope.grad = torch.zeros_like(theta_slope)\n",
    "            sp_row = {\n",
    "                \"density_idx\": int(row['density_idx']),\n",
    "                \"slope_idx\": int(row['slope_idx']),\n",
    "                \"distance\": xy_dist.item(),\n",
    "                #\"density_grad\": -theta_density.grad.item(),\n",
    "                #\"slope_grad\": -theta_slope.grad.item(),\n",
    "            }\n",
    "            sp_df.loc[len(sp_df)] = sp_row\n",
    "        out[(distance_key, seed)] = sp_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/user/miniconda3/envs/scrapl/lib/python3.10/site-packages/matplotlib/quiver.py:628: RuntimeWarning: Mean of empty slice.\n",
      "  amean = a[~self.Umask].mean()\n",
      "/Users/user/miniconda3/envs/scrapl/lib/python3.10/site-packages/numpy/core/_methods.py:189: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbcAAAKXCAYAAAAB0wRAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABkLElEQVR4nO3deVwTd/4/8FdACFeIoCigKAgtKp60aj2qtR6IVrHa2lpbwXst1lW/diu79aqr9PLoscW1h2CtWrVq3dqKaAFXW7daFbXrqngiqHgCAQmQfH5/+DA/UsAmIcOQ8fV8PObxMDOT97wT8+admcx8RiWEECAiIlIQJ7kTICIisjc2NyIiUhw2NyIiUhw2NyIiUhw2NyIiUhw2NyIiUhw2NyIiUhw2NyIiUhw2NyIiUhw2tzqgUqmwYMECudMgInposLlZITk5GSqVCiqVCvv27auyXAiBoKAgqFQqPPPMMzJkaL0lS5Zg27ZtD30OVPeUVE95eXlYsGABjh49+lDnUJ+wudnAzc0N69atqzI/MzMTly9fhlqtNpt/9+5dvPnmm3WVnlXqQ2OpDzmQfKytp/ooLy8PCxculL25yZ1DfcLmZoPBgwdj06ZNqKioMJu/bt06PPbYY/D39zeb7+bmhgYNGtR6u8XFxbWO4cjbt4bRaERpaancaZAFrK0nJSgpKZE7BYuVlpbCaDTKnYbV2NxsMHr0aNy8eRNpaWmmeWVlZdi8eTNeeumlKutX95tbbm4uJkyYgMDAQKjVaoSEhGDq1KkoKysD8P8P2WRmZuLVV19FkyZN0Lx5c9PzP/nkE0RERECtViMwMBDx8fG4c+eO2TbOnDmDkSNHwt/fH25ubmjevDlefPFFFBQUmPIqLi5GSkqK6fBQXFwcAGDBggVQqVT473//i5deegk+Pj7o1asXAOCpp57CU089VeV1xsXFITg42Gye0WjEBx98gPbt28PNzQ1+fn4YNGgQDh069Ic5VBevcm6/f4+nTZuGr776yvS+7Ny50/Rejx8/Hk2bNoVarUZERAS++OKLKnFJHtbWk9FoxIoVKxAREQE3Nzc0bdoUU6ZMwe3bt83WO3ToEKKiotC4cWO4u7sjJCQE48ePN1unuLgY//d//4egoCCo1WqEh4fj/fffx+9vlpKWloZevXqhYcOG8PLyQnh4OP76178CADIyMtClSxcAwLhx40yf4+TkZAD36qVdu3b49ddf0bt3b3h4eJieW9Pv8cHBwaY6uO/OnTuYOXMmgoODoVar0bx5c4wdOxY3btz4wxyqi3c/t8q1nJGRAZVKhQ0bNuDNN99Es2bN4OHhgcLCQgDAf/7zHwwaNAharRYeHh7o06cP9u/fXyVufVD73YmHUHBwMLp3747169cjOjoaAPDDDz+goKAAL774Ij788MMHPj8vLw9du3bFnTt3MHnyZLRu3Rq5ubnYvHkzSkpK4Orqalr31VdfhZ+fH+bNm2fac1qwYAEWLlyI/v37Y+rUqTh16hSSkpJw8OBB7N+/Hy4uLigrK0NUVBT0ej1ee+01+Pv7Izc3F9999x3u3LkDrVaLL7/8EhMnTkTXrl0xefJkAEBoaKhZrs8//zweeeQRLFmypErBW2LChAlITk5GdHQ0Jk6ciIqKCvz73//GgQMH8Pjjj1uUg6V+/PFHbNy4EdOmTUPjxo0RHByMa9eu4YknnjA1Pz8/P/zwww+YMGECCgsLMWPGDJu2RfZjbT1NmTIFycnJGDduHKZPn47z58/j448/xpEjR0yf//z8fAwcOBB+fn6YM2cOGjZsiAsXLmDLli2mOEIIDBs2DOnp6ZgwYQI6deqE1NRUvP7668jNzcXy5csBAL/99hueeeYZdOjQAW+99RbUajWys7NNf9TbtGmDt956C/PmzcPkyZPx5JNPAgB69Ohh2tbNmzcRHR2NF198ES+//DKaNm1q1Xuk0+nw5JNP4uTJkxg/fjwiIyNx48YNbN++HZcvX7YoB2ssWrQIrq6umD17NvR6PVxdXfHjjz8iOjoajz32GObPnw8nJyesXr0aTz/9NP7973+ja9euNm1LMoIstnr1agFAHDx4UHz88cdCo9GIkpISIYQQzz//vOjbt68QQoiWLVuKIUOGmJ4HQMyfP9/0eOzYscLJyUkcPHiwyjaMRqPZtnr16iUqKipMy/Pz84Wrq6sYOHCgMBgMpvkff/yxACC++OILIYQQR44cEQDEpk2bHviaPD09RWxsbJX58+fPFwDE6NGjqyzr06eP6NOnT5X5sbGxomXLlqbHP/74owAgpk+fXuPrfFAOv4/3+9wqAyCcnJzEb7/9ZjZ/woQJIiAgQNy4ccNs/osvvii0Wq3p/4/qni319O9//1sAEF999ZVZrJ07d5rN37p1qyl2TbZt2yYAiL///e9m85977jmhUqlEdna2EEKI5cuXCwDi+vXrNcY6ePCgACBWr15dZVmfPn0EALFy5coqy37/t+G+li1bmtXEvHnzBACxZcuWKuver6UH5fD7eJVzq1zL6enpAoBo1aqVWW0YjUbxyCOPiKioKLPaLSkpESEhIWLAgAFVYsuNhyVtNGrUKNy9exffffcdioqK8N1331V7COX3jEYjtm3bhqFDh+Lxxx+vsvz3h9smTZoEZ2dn0+Pdu3ejrKwMM2bMgJOTk9l63t7e2LFjBwBAq9UCAFJTU2t1fP9Pf/qTzc/95ptvoFKpMH/+/CrLfv867aFPnz5o27at6bEQAt988w2GDh0KIQRu3LhhmqKiolBQUIDDhw/bPQ+ynqX1tGnTJmi1WgwYMMDs//Oxxx6Dl5cX0tPTAQANGzYEAHz33XcoLy+vdpvff/89nJ2dMX36dLP5//d//wchBH744QezWN9++63Nvz2p1WqMGzfOpucC92qpY8eOePbZZ6ssk6KWYmNj4e7ubnp89OhRnDlzBi+99BJu3rxpet+Li4vRr18/7N27t979LsfmZiM/Pz/0798f69atw5YtW2AwGPDcc8/94fOuX7+OwsJCtGvXzqLthISEmD2+ePEiACA8PNxsvqurK1q1amVaHhISglmzZuGzzz5D48aNERUVhX/84x+m39ss9fvtW+Ps2bMIDAyEr6+vzTGs8ftcr1+/jjt37mDVqlXw8/Mzm+7/ocnPz6+T3OjBLK2nM2fOoKCgAE2aNKnyf6rT6Uz/n3369MHIkSOxcOFCNG7cGDExMVi9ejX0er0p1sWLFxEYGAiNRmO2jTZt2piWA8ALL7yAnj17YuLEiWjatClefPFFbNy40ao/5s2aNTP7ucFaZ8+etfhvhj38vpbOnDkD4F7T+/37/tlnn0Gv11v9t0Vq/M2tFl566SVMmjQJV69eRXR0tOkbnj1V/vZkraVLlyIuLg7ffvstdu3ahenTpyMxMREHDhwwOznF2u2rVKpqf38zGAw251qdmr6R1rSd3+d6/4/Pyy+/jNjY2Gqf06FDh1pkSPZkST0ZjUY0adIEX331VbUx/Pz8ANz77GzevBkHDhzAv/71L6SmpmL8+PFYunQpDhw4AC8vL4vzcnd3x969e5Geno4dO3Zg586d+Prrr/H0009j165dZkdWHhTDGnVZS9XlX1Mtvffee+jUqVO1sax5T+sCm1stPPvss5gyZQoOHDiAr7/+2qLn+Pn5wdvbGydOnLBpmy1btgQAnDp1Cq1atTLNLysrw/nz59G/f3+z9du3b4/27dvjzTffxE8//YSePXti5cqV+Pvf/w7AtkMaPj4+OHfuXJX597/p3hcaGorU1FTcunXrgXtvNeXg4+NT5QzQ6rZTEz8/P2g0GhgMhirvC9U/ltRTaGgodu/ejZ49e1rUMJ544gk88cQTWLx4MdatW4cxY8Zgw4YNmDhxIlq2bIndu3ejqKjIbO/tf//7H4D/X2sA4OTkhH79+qFfv35YtmwZlixZgr/97W9IT09H//79bT40WN1nvKysDFeuXKnyuv/ob8aDcnhQLVX+O1KT+yd5eXt7O0wt8bBkLXh5eSEpKQkLFizA0KFDLXqOk5MThg8fjn/961+m0+Erq26PqLL+/fvD1dUVH374odm6n3/+OQoKCjBkyBAAQGFhYZXrhtq3bw8nJyezQzOenp7VfugfJDQ0FP/73/9w/fp107ysrKwqpwSPHDkSQggsXLiwSozKudeUQ2hoKAoKCnDs2DHTvCtXrmDr1q0W5ens7IyRI0fim2++qfYPQ+X8SX6W1NOoUaNgMBiwaNGiKssqKipMn6Pbt29XqaX7exz3P/+DBw+GwWDAxx9/bLbe8uXLoVKpTGdu3rp1q8q2fh/L09MTAGyqpb1795rNW7VqVZU9t5EjRyIrK6vaz/791/mgHEJDQ3HgwAHTpUbAvd8jc3JyLMrzscceQ2hoKN5//33odLoqy+tjLXHPrZZqOtz1IEuWLMGuXbvQp08fTJ48GW3atMGVK1ewadMm7Nu374GHN/38/JCQkICFCxdi0KBBGDZsGE6dOoVPPvkEXbp0wcsvvwzg3mnx06ZNw/PPP49HH30UFRUV+PLLL01/8O977LHHsHv3bixbtgyBgYEICQlBt27dHpj/+PHjsWzZMkRFRWHChAnIz8/HypUrERERYboeBgD69u2LV155BR9++CHOnDmDQYMGwWg04t///jf69u2LadOmPTCHF198EW+88QaeffZZTJ8+HSUlJUhKSsKjjz5q8Ykgb7/9NtLT09GtWzdMmjQJbdu2xa1bt3D48GHs3r272j9cJJ8/qqc+ffpgypQpSExMxNGjRzFw4EC4uLjgzJkz2LRpEz744AM899xzSElJwSeffIJnn30WoaGhKCoqwqeffgpvb28MHjwYADB06FD07dsXf/vb33DhwgV07NgRu3btwrfffosZM2aY9lbeeust7N27F0OGDEHLli2Rn5+PTz75BM2bNzdd+xkaGoqGDRti5cqV0Gg08PT0RLdu3f7wN+uJEyfiT3/6E0aOHIkBAwYgKysLqampaNy4sdl6r7/+OjZv3oznn38e48ePx2OPPYZbt25h+/btWLlyJTp27PjAHCZOnIjNmzdj0KBBGDVqFM6ePYu1a9dafNmNk5MTPvvsM0RHRyMiIgLjxo1Ds2bNkJubi/T0dHh7e+Nf//qXRbHqjFynaTqiyqcuP8gfXQoghBAXL14UY8eOFX5+fkKtVotWrVqJ+Ph4odfrLdrWxx9/LFq3bi1cXFxE06ZNxdSpU8Xt27dNy8+dOyfGjx8vQkNDhZubm/D19RV9+/YVu3fvNovzv//9T/Tu3Vu4u7sLAKbThe+fbl/T6c9r164VrVq1Eq6urqJTp04iNTW12lP3KyoqxHvvvSdat24tXF1dhZ+fn4iOjha//vrrH+YghBC7du0S7dq1E66uriI8PFysXbu2xksB4uPjq8312rVrIj4+XgQFBQkXFxfh7+8v+vXrJ1atWlXt+lQ3bK0nIYRYtWqVeOyxx4S7u7vQaDSiffv24i9/+YvIy8sTQghx+PBhMXr0aNGiRQuhVqtFkyZNxDPPPCMOHTpkFqeoqEjMnDlTBAYGChcXF/HII4+I9957z+x09z179oiYmBgRGBgoXF1dRWBgoBg9erQ4ffq0Waxvv/1WtG3bVjRo0MDslPw+ffqIiIiIal+bwWAQb7zxhmjcuLHw8PAQUVFRIjs7u9pT92/evCmmTZsmmjVrJlxdXUXz5s1FbGys2WUuNeUghBBLly4VzZo1E2q1WvTs2VMcOnSoxksBarqE6MiRI2LEiBGiUaNGQq1Wi5YtW4pRo0aJPXv2VLu+nFRC2HBlLhERUT3G39yIiEhx2NyIiEhx2NyIiEhx2NyIiEhx2NyIiEhx2NyIiEhxHPoibqPRiLy8PGg0GklGxiayByEEioqKEBgYaHYnh/qG9UT1nTW15NDNLS8vD0FBQXKnQWSRnJwciweslgPriRyFJbXk0M3t97eqkMLvh8GRgre3t6Tx7485JxWp36NHHnlE0vjAvXt1ScVoNOL69et18nmtjbrIz5bh6qzVvn17SeMHBARIGj83N1fS+L/99puk8YF7w/9JwWg0Ijc316LPqkM3t7o4dFIXh5Gk3oYlt+SojQYNpP0Y1eY+WJaqi//n+n6ory7yq4v/y9rcJsoSHh4eksaXOn8l1JMln9X6+wMAERGRjdjciIhIcdjciIhIcdjciIhIcdjciIhIcWRtbkVFRZgxYwZatmwJd3d39OjRAwcPHpQzJSIiUgBZm9vEiRORlpaGL7/8EsePH8fAgQPRv39/ya/zICIiZZOtud29exfffPMN3n33XfTu3RthYWFYsGABwsLCkJSUJFdaRESkALJdxF1RUQGDwQA3Nzez+e7u7ti3b1+1z9Hr9dDr9abHhYWFkuZIpGSsJ1Iy2fbcNBoNunfvjkWLFiEvLw8GgwFr167Fzz//jCtXrlT7nMTERGi1WtPEcfCIbMd6IiWT9Te3L7/8EkIINGvWDGq1Gh9++CFGjx5d49AtCQkJKCgoME05OTl1nDGRcrCeSMlkHVsyNDQUmZmZKC4uRmFhIQICAvDCCy+gVatW1a6vVquhVqvrOEsiZWI9kZLVi+vcPD09ERAQgNu3byM1NRUxMTFyp0RERA5M1j231NRUCCEQHh6O7OxsvP7662jdujXGjRsnZ1pEROTgZN1zKygoQHx8PFq3bo2xY8eiV69eSE1NhYuLi5xpERGRg5N1z23UqFEYNWqUnCkQEZEC1Yvf3IiIiOyJzY2IiBSHzY2IiBSHzY2IiBSHzY2IiBRH1rMl7cXX17fGIbvsEVtq3t7eksb39PSUNH6TJk0kjV8XYx42bdpUstgGgwHXrl2TLL69Pfroo3B2dpYkdkhIiCRxK5P686jVaiWNX1BQIGn8xo0bSxofAMLDwyWJW1FRYfEwcdxzIyIixWFzIyIixWFzIyIixWFzIyIixWFzIyIixWFzIyIixWFzIyIixZG1uRkMBsydOxchISFwd3dHaGgoFi1aBCGEnGkREZGDk/Ui7nfeeQdJSUlISUlBREQEDh06hHHjxkGr1WL69OlypkZERA5M1ub2008/ISYmBkOGDAEABAcHY/369fjll1/kTIuIiBycrIcle/TogT179uD06dMAgKysLOzbtw/R0dFypkVERA5O1j23OXPmoLCwEK1bt4azszMMBgMWL16MMWPGVLu+Xq+HXq83PS4sLKyrVIkUh/VESibrntvGjRvx1VdfYd26dTh8+DBSUlLw/vvvIyUlpdr1ExMTodVqTVNdDKhLpFSsJ1IyWZvb66+/jjlz5uDFF19E+/bt8corr2DmzJlITEysdv2EhAQUFBSYJktHhyaiqlhPpGSyHpYsKSmpcqsaZ2dnGI3GatdXq9VQq9V1kRqR4rGeSMlkbW5Dhw7F4sWL0aJFC0RERODIkSNYtmwZxo8fL2daRETk4GRtbh999BHmzp2LV199Ffn5+QgMDMSUKVMwb948OdMiIiIHJ2tz02g0WLFiBVasWCFnGkREpDAcW5KIiBSHzY2IiBSHzY2IiBSHzY2IiBSHzY2IiBSHzY2IiBRH1ksB7MXX1xfOzs6SxG7YsKEkcSvTaDSSxvf09JQ0vp+fn6Txw8PDJY0PAP7+/pLFrqiowPHjxyWLb2+PP/44XF1dJYndrFkzSeJW1qhRI0nje3t7O3R8rVYraXwAaNGihSRxy8rKLF6Xe25ERKQ4bG5ERKQ4bG5ERKQ4bG5ERKQ4bG5ERKQ4bG5ERKQ4sje33NxcvPzyy2jUqBHc3d3Rvn17HDp0SO60iIjIgcl6ndvt27fRs2dP9O3bFz/88AP8/Pxw5swZ+Pj4yJkWERE5OFmb2zvvvIOgoCCsXr3aNC8kJETGjIiISAlkPSy5fft2PP7443j++efRpEkTdO7cGZ9++qmcKRERkQLI2tzOnTuHpKQkPPLII0hNTcXUqVMxffp0pKSkVLu+Xq9HYWGh2UREtmE9kZLJ2tyMRiMiIyOxZMkSdO7cGZMnT8akSZOwcuXKatdPTEyEVqs1TUFBQXWcMZFysJ5IyWxqbufOnbPLxgMCAtC2bVuzeW3atMGlS5eqXT8hIQEFBQWmKScnxy55ED2MWE+kZDadUBIWFoY+ffpgwoQJeO655+Dm5mbTxnv27IlTp06ZzTt9+jRatmxZ7fpqtRpqtdqmbRGROdYTKZlNe26HDx9Ghw4dMGvWLPj7+2PKlCn45ZdfrI4zc+ZMHDhwAEuWLEF2djbWrVuHVatWIT4+3pa0iIiIANjY3Dp16oQPPvgAeXl5+OKLL3DlyhX06tUL7dq1w7Jly3D9+nWL4nTp0gVbt27F+vXr0a5dOyxatAgrVqzAmDFjbEmLiIgIQC1PKGnQoAFGjBiBTZs24Z133kF2djZmz56NoKAgjB07FleuXPnDGM888wyOHz+O0tJSnDx5EpMmTapNSkRERLVrbocOHcKrr76KgIAALFu2DLNnz8bZs2eRlpaGvLw8xMTE2CtPIiIii9l0QsmyZcuwevVqnDp1CoMHD8aaNWswePBgODnd65UhISFITk5GcHCwPXMlIiKyiE3NLSkpCePHj0dcXBwCAgKqXadJkyb4/PPPa5UcERGRLWxqbmlpaWjRooVpT+0+IQRycnLQokULuLq6IjY21i5JEhERWcOm39xCQ0Nx48aNKvNv3brFgY+JiEh2NjU3IUS183U6nc0XdBMREdmLVYclZ82aBQBQqVSYN28ePDw8TMsMBgP+85//oFOnTnZN0BKenp5wdnaWJLZGo5EkbmWenp6Sxvfy8pI0vp+fn6Tx6+Ks22+++Uay2GVlZZLFlkLjxo0lG7mkcePGksStTOqalXpUF6nrtVGjRpLGB+6NGyyF0tJSi9e1qrkdOXIEwL09t+PHj8PV1dW0zNXVFR07dsTs2bOtCUlERGR3VjW39PR0AMC4cePwwQcfwNvbW5KkiIiIasOmsyUr3zmbiIiovrG4uY0YMQLJycnw9vbGiBEjHrjuli1bap0YERGRrSxublqtFiqVyvRvIiKi+sri5lb5UCQPSxIRUX1m03Vud+/eRUlJienxxYsXsWLFCuzatctuiREREdnKpuYWExODNWvWAADu3LmDrl27YunSpYiJiUFSUpLFcRYsWACVSmU2tW7d2paUiKAr02FR5iIELQ+C81vOCFoehEWZi6Ar08mdGpFDKTWW4rvC7/CXK3/B5NzJ+MuVv+C7wu9QarT8OjO52Xwn7ieffBIAsHnzZvj7++PixYtYs2YNPvzwQ6tiRURE4MqVK6Zp3759tqREDzldmQ59kvtgQeYCXC68DKMw4nLhZSzIXIA+yX3Y4IgsVGosxXs33sP2ou24bbwNAYHbxtvYXrQd7914z2EanE3NraSkxDQKwK5duzBixAg4OTnhiSeewMWLF62K1aBBA/j7+5umuhjBgJRn+c/LcfTqURiF+cgIRmHE0atHsfzn5TJlRuRYdut2I6c8BwLmwywKCOSU52C3brdMmVnHpuYWFhaGbdu2IScnB6mpqRg4cCAAID8/3+oLu8+cOYPAwEC0atUKY8aMwaVLl2pcV6/Xo7Cw0GwiAoBVh1dVaWz3GYURqw6vquOM6j/WE1Vnb/HeKo3tPgGBvcV76zgj29jU3ObNm4fZs2cjODgY3bp1Q/fu3QHc24vr3LmzxXG6deuG5ORk7Ny5E0lJSTh//jyefPJJFBUVVbt+YmIitFqtaQoKCrIlfVKgvKK8Wi1/GLGeqDp3jHdqtby+sKm5Pffcc7h06RIOHTqEnTt3mub369cPy5dbfvgnOjoazz//PDp06ICoqCh8//33uHPnDjZu3Fjt+gkJCSgoKDBNOTk5tqRPChSoCazV8ocR64mq09CpYa2W1xc2Db8FwPQbWWVdu3atVTINGzbEo48+iuzs7GqXq9VqyUfkJsc0OXIyFmQuqPbQpJPKCZMjJ8uQVf3GeqLq9Pbsje1F26s9NKmCCr09e8uQlfVs2nMrLi7G3Llz0aNHD4SFhaFVq1Zmk610Oh3Onj2LgIAAm2PQw2lm95no5N8JTirzj7STygmd/DthZveZMmVG5Fj6e/VHkEsQVFCZzVdBhSCXIPT36i9TZtaxac9t4sSJyMzMxCuvvIKAgADTsFzWmj17NoYOHYqWLVsiLy8P8+fPh7OzM0aPHm1TPHp4ebl6ITMuE8t/Xo5Vh1chrygPgZpATI6cjJndZ8LLVdp7ZBEphZuTG15v/Dp263Zjb/Fe3DHeQUOnhujt2Rv9vfrDzckxbkhtU3P74YcfsGPHDvTs2bNWG798+TJGjx6Nmzdvws/PD7169cKBAwckv/klKZOXqxfm9pmLuX3myp0KkUNzc3LDM97P4BnvZ+ROxWY2NTcfHx/4+vrWeuMbNmyodQwiIqLfs+k3t0WLFmHevHlm40sSERHVFzbtuS1duhRnz55F06ZNERwcDBcXF7Plhw8ftktyREREtrCpuQ0fPtzOaRAREdmPTc1t/vz59s6DiIjIbmz6zQ24d6ubzz77DAkJCbh16xaAe4cjc3Nz7ZYcERGRLWzaczt27Bj69+8PrVaLCxcuYNKkSfD19cWWLVtw6dIl073e6kpFRQWEqH6gz9oqLXWM2zs8iMFgkDR+fn6+pPHT0tIkjQ/cG0BAKuXl5ZLFlkK7du3g7u4uSezaDPJgKa1WK2l8V1dXSeO7uUl7HZlU/7eVeXlJc12pNScx2rTnNmvWLMTFxeHMmTNm/xGDBw/G3r2OMWI0EREpl03N7eDBg5gyZUqV+c2aNcPVq1drnRQREVFt2NTc1Gp1tfd+On36NEcXISIi2dnU3IYNG4a33nrL9FuCSqXCpUuX8MYbb2DkyJF2TZCIiMhaNjW3pUuXQqfTwc/PD3fv3kWfPn0QFhYGjUaDxYsX2ztHIiIiq9h0tqRWq0VaWhr279+PrKws6HQ6REZGon9/x7gVAhERKZvVzc1oNCI5ORlbtmzBhQsXoFKpEBISAn9/fwghbL79DRERkb1YdVhSCIFhw4Zh4sSJyM3NRfv27REREYGLFy8iLi4Ozz77rM2JvP3221CpVJgxY4bNMYiIiAAr99ySk5Oxd+9e7NmzB3379jVb9uOPP2L48OFYs2YNxo4da1USBw8exD//+U906NDBqucRERFVx6o9t/Xr1+Ovf/1rlcYGAE8//TTmzJmDr776yqoEdDodxowZg08//RQ+Pj5WPZeIiKg6VjW3Y8eOYdCgQTUuj46ORlZWllUJxMfHY8iQIRadjKLX61FYWGg2EZFtWE+kZFY1t1u3bqFp06Y1Lm/atClu375tcbwNGzbg8OHDSExMtGj9xMREaLVa0xQUFGTxtojIHOuJlMyq5mYwGNCgQc0/0zk7O6OiosKiWDk5Ofjzn/+Mr776yuKBQhMSElBQUGCacnJyLHoeEVXFeiIls+qEEiEE4uLioFarq12u1+stjvXrr78iPz8fkZGRpnkGgwF79+7Fxx9/DL1eD2dnZ7PnqNXqGrdNRNZhPZGSWdXcYmNj/3AdS8+U7NevH44fP242b9y4cWjdujXeeOONKo2NiIjIUlY1t9WrV9ttwxqNBu3atTOb5+npiUaNGlWZT0REZA2b78RNRERUX9k0tqRUMjIy5E6BiIgUgHtuRESkOGxuRESkOGxuRESkOGxuRESkOGxuRESkOPXqbElblZWVSXbRd3FxsSRxKysvL5c0fllZmaTxr127Jmn8X3/9VdL4wL27U0jF0iHp6gtfX194eHhIEtvV1VWSuJUZjUaHji/1e+Tt7S1pfAAWD6loLWvee+65ERGR4rC5ERGR4rC5ERGR4rC5ERGR4rC5ERGR4rC5ERGR4rC5ERGR4sja3BITE9GlSxdoNBo0adIEw4cPx6lTp+RMiYiIFEDW5paZmYn4+HgcOHAAaWlpKC8vx8CBA+vkwmkiIlIuWUco2blzp9nj5ORkNGnSBL/++it69+4tU1ZEROTo6tVvbgUFBQDuDf9DRERkq3oztqTRaMSMGTPQs2dPtGvXrtp19Ho99Hq96XFhYWFdpUekOKwnUrJ6s+cWHx+PEydOYMOGDTWuk5iYCK1Wa5qCgoLqMEMiZWE9kZLVi+Y2bdo0fPfdd0hPT0fz5s1rXC8hIQEFBQWmKScnpw6zJFIW1hMpmayHJYUQeO2117B161ZkZGQgJCTkgeur1Wqo1eo6yo5I2VhPpGSyNrf4+HisW7cO3377LTQaDa5evQoA0Gq1cHd3lzM1IiJyYLIelkxKSkJBQQGeeuopBAQEmKavv/5azrSIiMjByX5YkoiIyN7qxQklRERE9sTmRkREisPmRkREisPmRkREisPmRkREisPmRkREilNvBk6uDSGEZJcVlJWVSRK3MoPBIGl8qV+DSqWSNH5d/B9IOWhwRUWFZLGlEBwcDC8vL0li+/n5SRK3Mjc3N0nju7i4SBq/vLxc0vgNGzaUND4AaDQaSeIWFRVZvC733IiISHHY3IiISHHY3IiISHHY3IiISHHY3IiISHHY3IiISHHY3IiISHFkbW579+7F0KFDERgYCJVKhW3btsmZDhERKYSsza24uBgdO3bEP/7xDznTICIihZF1hJLo6GhER0fLmQIRESmQQw2/pdfrodfrTY+lHDKJSOlYT6RkDnVCSWJiIrRarWkKCgqSOyUih8V6IiVzqOaWkJCAgoIC05STkyN3SkQOi/VESuZQhyXVajXUarXcaRApAuuJlMyh9tyIiIgsIeuem06nQ3Z2tunx+fPncfToUfj6+qJFixYyZkZERI5M1uZ26NAh9O3b1/R41qxZAIDY2FgkJyfLlBURETk6WZvbU089JdkdtImI6OHF39yIiEhx2NyIiEhx2NyIiEhx2NyIiEhxHOoi7t+7fzKK0WiUbBtOTuz/f6SiokLS+OXl5ZLGB6R9DQaDAQDq/clT9/PT6XSSbaMuxq8sKyuTNL6Li4uk8aX+vN+9e1fS+ABQVFQkSdz7n01Lasmhm9v9N/DcuXMyZ0L0x4qKiqDVauVOo0b36+mpp56SNxGiP2BJLalEff86+QBGoxF5eXnQaDRQqVR/uH5hYSGCgoKQk5MDb29vSXKSehuML/82rI0vhEBRURECAwPr9ZGA+lZPjC//NupbfGtqyaH33JycnNC8eXOrn+ft7S3Zh62utsH48m/Dmvj1eY/tvvpaT4wv/zbqU3xLa6n+fo0kIiKyEZsbEREpzkPV3NRqNebPny/pbT6k3gbjy7+NungNjsDR32dHj18X23Dk+A59QgkREVF1Hqo9NyIiejiwuRERkeKwuRERkeKwuRERkeKwuT1kVCoVtm3bJncaRA6PtVS/sbkpzPXr1zF16lS0aNECarUa/v7+iIqKwv79++VOjcihsJYcm0MPv0VVjRw5EmVlZUhJSUGrVq1w7do17NmzBzdv3pQ7NSKHwlpycIIU4/bt2wKAyMjIqHEdAGLr1q2mx8eOHRN9+/YVbm5uwtfXV0yaNEkUFRWZlsfGxoqYmBixYMEC0bhxY6HRaMSUKVOEXq83rWMwGMSSJUtEcHCwcHNzEx06dBCbNm2S5DUS1QXWkuPjYUkF8fLygpeXF7Zt2wa9Xv+H6xcXFyMqKgo+Pj44ePAgNm3ahN27d2PatGlm6+3ZswcnT55ERkYG1q9fjy1btmDhwoWm5YmJiVizZg1WrlyJ3377DTNnzsTLL7+MzMxMu79GorrAWlIAubsr2dfmzZuFj4+PcHNzEz169BAJCQkiKyvLtByVvm2uWrVK+Pj4CJ1OZ1q+Y8cO4eTkJK5evSqEuPdt09fXVxQXF5vWSUpKEl5eXsJgMIjS0lLh4eEhfvrpJ7M8JkyYIEaPHi3hKyWSFmvJsXHPTWFGjhyJvLw8bN++HYMGDUJGRgYiIyORnJxcZd2TJ0+iY8eO8PT0NM3r2bMnjEYjTp06ZZrXsWNHeHh4mB53794dOp0OOTk5yM7ORklJCQYMGGD6tuvl5YU1a9bg7Nmzkr5WIimxlhwbTyhRIDc3NwwYMAADBgzA3LlzMXHiRMyfPx9xcXF239b9277v2LEDzZo1M1v2sA8sTI6PteS4uOf2EGjbti2Ki4urzG/Tpg2ysrLMlu3fvx9OTk4IDw83zcvKysLdu3dNjw8cOAAvLy8EBQWhbdu2UKvVuHTpEsLCwsymoKAgaV8YUR1jLTkONjcFuXnzJp5++mmsXbsWx44dw/nz57Fp0ya8++67iImJqbL+mDFj4ObmhtjYWJw4cQLp6el47bXX8Morr6Bp06am9crKyjBhwgT897//xffff4/58+dj2rRpcHJygkajwezZszFz5kykpKTg7NmzOHz4MD766COkpKTU5csnshvWkgLI/aMf2U9paamYM2eOiIyMFFqtVnh4eIjw8HDx5ptvipKSEiGE7acvz5s3TzRq1Eh4eXmJSZMmidLSUtM6RqNRrFixQoSHhwsXFxfh5+cnoqKiRGZmZp29diJ7Yi05Pt7PjR4oLi4Od+7c4TBDRLXEWqpbPCxJRESKw+ZGRESKw8OSRESkONxzIyIixWFzIyIixWFzIyIixWFzIyIixWFzIyIixWFzIyIixWFzIyIixXHoW94YjUbk5eVBo9FApVLJnQ5RtYQQKCoqQmBgIJyc6u/3SdYT1XfW1JJDN7e8vDzeCoIcRk5ODpo3by53GjViPZGjsKSWHLq5aTQaAEB4eDicnZ0l2cbjjz8uSdzKysvLJY0v9SA0wcHBksbXarWSxgeAEydOSBa7vLwcGzduNH1e66v7+XXv3h0NGkjzp+Hdd9+VJG5l+fn5ksbfv3+/pPEHDRokafyGDRtKGh8Avv/+e0nilpaW4u9//7tFteTQze3+oRNnZ2fJmpurq6skcSuT+hCQ1M1N6rsEu7m5SRofUMb/c23dz69BgwaSNTcvLy9J4lZW3c1E7Unqz7vU71Fd/B9IXbOW1FL9/QGAiIjIRmxuRESkOGxuRESkOGxuRESkOGxuRESkOLI2t6KiIsyYMQMtW7aEu7s7evTogYMHD8qZEhERKYCszW3ixIlIS0vDl19+iePHj2PgwIHo378/cnNz5UyLiIgcnGzN7e7du/jmm2/w7rvvonfv3ggLC8OCBQsQFhaGpKQkudIiIiIFkO0i7oqKChgMhioX+7m7u2Pfvn3VPkev10Ov15seFxYWSpojkZKxnkjJZNtz02g06N69OxYtWoS8vDwYDAasXbsWP//8M65cuVLtcxITE6HVak0Tx8Ejsh3riZRM1t/cvvzySwgh0KxZM6jVanz44YcYPXp0jaM9JyQkoKCgwDTl5OTUccZEysF6IiWTdWzJ0NBQZGZmori4GIWFhQgICMALL7yAVq1aVbu+Wq2WfFw3oocF64mUrF5c5+bp6YmAgADcvn0bqampiImJkTslIiJyYLLuuaWmpkIIgfDwcGRnZ+P1119H69atMW7cODnTIiIiByfrnltBQQHi4+PRunVrjB07Fr169UJqaipcXFzkTIuIiBycrHtuo0aNwqhRo+RMgYiIFKhe/OZGRERkT2xuRESkOGxuRESkOGxuRESkOGxuRESkOLKeLWkvFy5cgEqlkiR2RESEJHEr+/3g0faWn58vaXwPDw9J40udPwDs379fstgGg0Gy2FK4e/cuGjRw3D8NJSUlksbXaDSSxpf6vX/kkUckjQ/c+5sshbKyMovX5Z4bEREpDpsbEREpDpsbEREpDpsbEREpDpsbEREpDpsbEREpDpsbEREpjqzNzWAwYO7cuQgJCYG7uztCQ0OxaNEiCCHkTIuIiBycrFdqvvPOO0hKSkJKSgoiIiJw6NAhjBs3DlqtFtOnT5czNSIicmCyNreffvoJMTExGDJkCAAgODgY69evxy+//CJnWkRE5OBkPSzZo0cP7NmzB6dPnwYAZGVlYd++fYiOjq52fb1ej8LCQrOJiGzDeiIlk3XPbc6cOSgsLETr1q3h7OwMg8GAxYsXY8yYMdWun5iYiIULF9ZxlkTKxHoiJZN1z23jxo346quvsG7dOhw+fBgpKSl4//33kZKSUu36CQkJKCgoME05OTl1nDGRcrCeSMlk3XN7/fXXMWfOHLz44osAgPbt2+PixYtITExEbGxslfXVajXUanVdp0mkSKwnUjJZ99xKSkrg5GSegrOzM4xGo0wZERGREsi65zZ06FAsXrwYLVq0QEREBI4cOYJly5Zh/PjxcqZFREQOTtbm9tFHH2Hu3Ll49dVXkZ+fj8DAQEyZMgXz5s2TMy0iInJwsjY3jUaDFStWYMWKFXKmQURECsOxJYmISHHY3IiISHHY3IiISHHY3IiISHHY3IiISHHY3IiISHFkvRTAXlxdXaFSqSSJ3bJlS0niVtapUydJ4x84cEDS+BEREZLG7969u6TxAeDEiROSxa6oqMDZs2cli29v+fn5VUYOspe0tDRJ4lYWEBAgafwmTZpIGj8oKEjS+NevX5c0PgAUFRVJEresrMzidbnnRkREisPmRkREisPmRkREisPmRkREisPmRkREisPmRkREiiN7c8vNzcXLL7+MRo0awd3dHe3bt8ehQ4fkTouIiByYrNe53b59Gz179kTfvn3xww8/wM/PD2fOnIGPj4+caRERkYOTtbm98847CAoKwurVq03zQkJCZMyIiIiUQNbDktu3b8fjjz+O559/Hk2aNEHnzp3x6aefypkSEREpgKzN7dy5c0hKSsIjjzyC1NRUTJ06FdOnT0dKSkq16+v1ehQWFppNRGQb1hMpmU3N7dy5c3bZuNFoRGRkJJYsWYLOnTtj8uTJmDRpElauXFnt+omJidBqtaZJ6jHYiJSM9URKZlNzCwsLQ9++fbF27VqUlpbavPGAgAC0bdvWbF6bNm1w6dKlatdPSEhAQUGBacrJybF520QPO9YTKZlNze3w4cPo0KEDZs2aBX9/f0yZMgW//PKL1XF69uyJU6dOmc07ffp0jSPxq9VqeHt7m01EZBvWEymZTc2tU6dO+OCDD5CXl4cvvvgCV65cQa9evdCuXTssW7bM4lsqzJw5EwcOHMCSJUuQnZ2NdevWYdWqVYiPj7clLSIiIgC1PKGkQYMGGDFiBDZt2oR33nkH2dnZmD17NoKCgjB27FhcuXLlgc/v0qULtm7divXr16Ndu3ZYtGgRVqxYgTFjxtQmLSIiesjVqrkdOnQIr776KgICArBs2TLMnj0bZ8+eRVpaGvLy8hATE/OHMZ555hkcP34cpaWlOHnyJCZNmlSblIiIiGy7iHvZsmVYvXo1Tp06hcGDB2PNmjUYPHiw6e69ISEhSE5ORnBwsD1zJSIisohNzS0pKQnjx49HXFxcjbd0b9KkCT7//PNaJUdERGQLm5pbWloaWrRoYdpTu08IgZycHLRo0QKurq6IjY21S5JERETWsOk3t9DQUNy4caPK/Fu3bnFsSCIikp1NzU0IUe18nU4HNze3WiVERERUW1Ydlpw1axYAQKVSYd68efDw8DAtMxgM+M9//oNOnTrZNUEiIiJrWdXcjhw5AuDentvx48fh6upqWubq6oqOHTti9uzZ9s3QAkajESqVSpLYERERksStLC8vT9L4Ug+IW9NwafZSF5eHSHkPwfLycsliS0Gn01X5Pd1e9Hq9JHEra9iwoaTxjUajpPGleu/r0t69eyWJa817b1VzS09PBwCMGzcOH3zwAYfrISKiesmmsyUr31yUiIiovrG4uY0YMQLJycnw9vbGiBEjHrjuli1bap0YERGRrSxublqt1vS7llarlSwhIiKi2rK4uVU+FMnDkkREVJ/ZdFrO3bt3UVJSYnp88eJFrFixArt27bJbYkRERLayqbnFxMRgzZo1AIA7d+6ga9euWLp0KWJiYpCUlGTXBImIiKxl8524n3zySQDA5s2b4e/vj4sXL2LNmjX48MMPLY6zYMECqFQqs6l169a2pEQEXZkOizIXIWh5EJzfckbQ8iAsylwEXZlO7tSIHIquTIf3f3kf7b9oD7+P/dD+i/Z4/5f3HaqWbLoUoKSkBBqNBgCwa9cujBgxAk5OTnjiiSdw8eJFq2JFRERg9+7d/z+hBjalRA85XZkOfZL74OjVozCKexd6Xi68jAWZC7Dt1DZkxmXCy9VL5iyJ6j9dmQ7DtgzD8RvHTbWUV5yHd355B9+f+x7bR2x3iFqyac8tLCwM27ZtQ05ODlJTUzFw4EAAQH5+vtUXdjdo0AD+/v6mqXHjxrakRA+55T8vN2ts9xmFEUevHsXyn5fLlBmRY1l5dKVZY7vPKIw4fuM4Vh5dKVNm1rGpuc2bNw+zZ89GcHAwunXrhu7duwO4txfXuXNnq2KdOXMGgYGBaNWqFcaMGfPAoZz0ej0KCwvNJiIAWHV4VZVivM8ojFh1eFUdZ1T/sZ6oOiknUh5YSyknUuo4I9vY1Nyee+45XLp0CYcOHcLOnTtN8/v164flyy3/htytWzckJydj586dSEpKwvnz5/Hkk0+iqKio2vUTExOh1WpNU1BQkC3pkwLlFT14fM4/Wv4wYj1Rda6WXK3V8vrC5hE6/f390blzZ7NBPrt27WrVCSHR0dF4/vnn0aFDB0RFReH777/HnTt3sHHjxmrXT0hIQEFBgWnKycmxNX1SmEBNYK2WP4xYT1Qdfw//Wi2vL2w6e6O4uBhvv/029uzZg/z8/CojNZ87d86mZBo2bIhHH30U2dnZ1S5Xq9VQq9U2xSZlmxw5GQsyF1R7OMVJ5YTJkZNlyKp+Yz1RdWLbxeKdX96psZZi28XKkJX1bGpuEydORGZmJl555RUEBATY7XYzOp0OZ8+exSuvvGKXePTwmNl9Jrad2lblpBInlRM6+XfCzO4zZcyOyHH8qdOf8P2576ucVOKkckL7xu3xp05/kjE7y9nU3H744Qfs2LEDPXv2rNXGZ8+ejaFDh6Jly5bIy8vD/Pnz4ezsjNGjR9cqLj18vFy9kBmXieU/L8eqw6uQV5SHQE0gJkdOxszuMx3i1GWi+sDL1QvbR2zHyqMrkXIiBVdLrsLfwx+x7WLxp05/cphasqm5+fj4wNfXt9Ybv3z5MkaPHo2bN2/Cz88PvXr1woEDB+Dn51fr2PTw8XL1wtw+czG3z1y5UyFyaF6uXpjddTZmd637m0/bi03NbdGiRZg3bx5SUlLg4eFh88Y3bNhg83OJiIhqYlNzW7p0Kc6ePYumTZsiODgYLi4uZssPHz5sl+SIiIhsYVNzGz58uJ3TICIish+bmtv8+fPtnQcREZHd2HwR9507d/DZZ58hISEBt27dAnDvcGRubq7dkiMiIrKFTXtux44dQ//+/aHVanHhwgVMmjQJvr6+2LJlCy5dumS61xsREZEcbGpus2bNQlxcHN59913TrW8AYPDgwXjppZfslpylhBCSxd62bZtkse/r27evpPEvXLggafzIyEhJ49+4cUPS+ADscmlLTcrKyiSLLYWAgAA4OztLErtp06aSxK1M6s/LtWvXJI2fnp4uafzaXp9siVGjRkkSV6/XY8WKFRata9NhyYMHD2LKlClV5jdr1gxXrzrGoJpERKRcNjU3tVpd7e0xTp8+zQuwiYhIdjY1t2HDhuGtt95CeXk5AEClUuHSpUt44403MHLkSLsmSEREZC2bmtvSpUuh0+ng5+eHu3fvok+fPggLC4NGo8HixYvtnSMREZFVbDqhRKvVIi0tDfv370dWVhZ0Oh0iIyPRv39/e+dHRERkNaubm9FoRHJyMrZs2YILFy5ApVIhJCQE/v7+EELY7fY3REREtrLqsKQQAsOGDcPEiRORm5uL9u3bIyIiAhcvXkRcXByeffZZqfIkIiKymFXNLTk5GXv37sWePXtw5MgRrF+/Hhs2bEBWVhZ2796NH3/80eYLuN9++22oVCrMmDHDpucTERHdZ1VzW79+Pf76179We9Hx008/jTlz5uCrr76yOomDBw/in//8Jzp06GD1c4mIiH7PquZ27NgxDBo0qMbl0dHRyMrKsioBnU6HMWPG4NNPP4WPj49VzyUiIqqOVc3t1q1bDxw+p2nTprh9+7ZVCcTHx2PIkCEWnWmp1+tRWFhoNhGRbVhPpGRWNTeDwYAGDWo+wdLZ2RkVFRUWx9uwYQMOHz6MxMREi9ZPTEyEVqs1TUFBQRZvi4jMsZ5Iyay6FEAIgbi4OKjV6mqX6/V6i2Pl5OTgz3/+M9LS0uDm5mbRcxISEjBr1izT48LCQhYkkY1YT6RkVjW32NjYP1xn7NixFsX69ddfkZ+fbzaivMFgwN69e/Hxxx9Dr9dXGZlcrVbX2FiJyDqsJ1Iyq5rb6tWr7bbhfv364fjx42bzxo0bh9atW+ONN96Q7JYbRESkfDYNv2UPGo0G7dq1M5vn6emJRo0aVZlPRERkDZsGTiYiIqrPZNtzq05GRobcKRARkQJwz42IiBSHzY2IiBSHzY2IiBSHzY2IiBSHzY2IiBSnXp0taSuDwSDZHcA7d+4sSdzK7ty5I2n80tJSSePn5eVJGn/Tpk2Sxgfu3auQ7omOjpZs5BJPT09J4lbm5CTtd/YbN25IGt/FxUXS+Ldu3ZI0PgCcP39ekrjl5eUWr8s9NyIiUhw2NyIiUhw2NyIiUhw2NyIiUhw2NyIiUhw2NyIiUhw2NyIiUhxZm1tiYiK6dOkCjUaDJk2aYPjw4Th16pScKRERkQLI2twyMzMRHx+PAwcOIC0tDeXl5Rg4cCCKi4vlTIuIiBycrCOU7Ny50+xxcnIymjRpgl9//RW9e/eWKSsiInJ09eo3t4KCAgCAr6+vzJkQEZEjqzdjSxqNRsyYMQM9e/ZEu3btql1Hr9dDr9ebHhcWFtZVekSKw3oiJas3e27x8fE4ceIENmzYUOM6iYmJ0Gq1pikoKKgOMyRSFtYTKVm9aG7Tpk3Dd999h/T0dDRv3rzG9RISElBQUGCacnJy6jBLImVhPZGSyXpYUgiB1157DVu3bkVGRgZCQkIeuL5arZbsVhxEDxvWEymZrM0tPj4e69atw7fffguNRoOrV68CALRaLdzd3eVMjYiIHJishyWTkpJQUFCAp556CgEBAabp66+/ljMtIiJycLIfliQiIrK3enFCCRERkT2xuRERkeKwuRERkeKwuRERkeKwuRERkeKwuRERkeLUm4GTa2PixImSjbTg7+8vSdzK/Pz8JI3v7OwsafyAgABJ49fFDWyjo6Mli11eXo7du3dLFt/eLly4ABcXF0liP/bYY5LErUzqz3vLli0ljX/79m1J45eXl0saHwA2btwo+Tb+CPfciIhIcdjciIhIcdjciIhIcdjciIhIcdjciIhIcdjciIhIcWRtbnv37sXQoUMRGBgIlUqFbdu2yZkOEREphKzNrbi4GB07dsQ//vEPOdMgIiKFkfUi7ujoaEkvniUioocTf3MjIiLFcajht/R6PfR6velxYWGhjNkQOTbWEymZQ+25JSYmQqvVmqagoCC5UyJyWKwnUjKHam4JCQkoKCgwTTk5OXKnROSwWE+kZA51WFKtVks2+j/Rw4b1REoma3PT6XTIzs42PT5//jyOHj0KX19ftGjRQsbMiIjIkcna3A4dOoS+ffuaHs+aNQsAEBsbi+TkZJmyIiIiRydrc3vqqacghJAzBSIiUiCHOqGEiIjIEmxuRESkOGxuRESkOGxuRESkOGxuRESkOA51Effv3T/TsvL4ePZ29+5dyWLfV1JSImn80tJSSeNL/R5J+f97X3l5uWSxKyoqAKDenxl8Pz8p3wupP+sA4OzsLGl8qT/vjl6vdcGSWlKJ+l5xD3D58mWOh0cOIycnB82bN5c7jRqxnshRWFJLDt3cjEYj8vLyoNFooFKp/nD9wsJCBAUFIScnB97e3pLkJPU2GF/+bVgbXwiBoqIiBAYGwsmp/v4SUN/qifHl30Z9i29NLTn0YUknJyebvgl7e3tL9mGrq20wvvzbsCa+VquVLA97qa/1xPjyb6M+xbe0lurv10giIiIbsbkREZHiPFTNTa1WY/78+ZLe5kPqbTC+/Nuoi9fgCBz9fXb0+HWxDUeO79AnlBAREVXnodpzIyKihwObGxERKQ6bGxERKQ6b20NGpVJh27ZtcqdB5PBYS/Ubm5vCXL9+HVOnTkWLFi2gVqvh7++PqKgo7N+/X+7UiBwKa8mxOfQIJVTVyJEjUVZWhpSUFLRq1QrXrl3Dnj17cPPmTblTI3IorCUHJ0gxbt++LQCIjIyMGtcBILZu3Wp6fOzYMdG3b1/h5uYmfH19xaRJk0RRUZFpeWxsrIiJiRELFiwQjRs3FhqNRkyZMkXo9XrTOgaDQSxZskQEBwcLNzc30aFDB7Fp0yZJXiNRXWAtOT4ellQQLy8veHl5Ydu2bRbdJqa4uBhRUVHw8fHBwYMHsWnTJuzevRvTpk0zW2/Pnj04efIkMjIysH79emzZsgULFy40LU9MTMSaNWuwcuVK/Pbbb5g5cyZefvllZGZm2v01EtUF1pICyN1dyb42b94sfHx8hJubm+jRo4dISEgQWVlZpuWo9G1z1apVwsfHR+h0OtPyHTt2CCcnJ3H16lUhxL1vm76+vqK4uNi0TlJSkvDy8hIGg0GUlpYKDw8P8dNPP5nlMWHCBDF69GgJXymRtFhLjo17bgozcuRI5OXlYfv27Rg0aBAyMjIQGRmJ5OTkKuuePHkSHTt2hKenp2lez549YTQacerUKdO8jh07wsPDw/S4e/fu0Ol0yMnJQXZ2NkpKSjBgwADTt10vLy+sWbMGZ8+elfS1EkmJteTYeEKJArm5uWHAgAEYMGAA5s6di4kTJ2L+/PmIi4uz+7Z0Oh0AYMeOHWjWrJnZsod97EVyfKwlx8U9t4dA27ZtUVxcXGV+mzZtkJWVZbZs//79cHJyQnh4uGleVlaW2a3pDxw4AC8vLwQFBaFt27ZQq9W4dOkSwsLCzCbe1ZmUhrXkONjcFOTmzZt4+umnsXbtWhw7dgznz5/Hpk2b8O677yImJqbK+mPGjIGbmxtiY2Nx4sQJpKen47XXXsMrr7yCpk2bmtYrKyvDhAkT8N///hfff/895s+fj2nTpsHJyQkajQazZ8/GzJkzkZKSgrNnz+Lw4cP46KOPkJKSUpcvn8huWEsKIPePfmQ/paWlYs6cOSIyMlJotVrh4eEhwsPDxZtvvilKSkqEELafvjxv3jzRqFEj4eXlJSZNmiRKS0tN6xiNRrFixQoRHh4uXFxchJ+fn4iKihKZmZl19tqJ7Im15Ph4yxt6oLi4ONy5c4fDDBHVEmupbvGwJBERKQ6bGxERKQ4PSxIRkeJwz42IiBSHzY2IiBSHzY2IiBSHzY2IiBSHzY2IiBSHzY2IiBSHzY2IiBSHzY2IiBTHoe/nZjQakZeXB41GA5VKJXc6RNUSQqCoqAiBgYFwcqq/3ydZT1TfWVNLDt3c8vLyeJ8jchg5OTlo3ry53GnUiPVEjsKSWnLo5qbRaAAAfn5+kn0jnjZtmiRxKysqKpI0fuPGjSWNL/XeSNeuXSWNDwC9e/eWfBv3P6/11f38unfvjgYNpPnT8NJLL0kSt7IrV65IGl/qEQsDAgIkjS/V/21lkydPljS+JbXk0M3t/qETJycnyf7Aurm5SRK3srKyMknju7u7Sxpf6ubm5eUlafy6Ut8P9d3Pr0GDBpL9AZT6swhIX7NSNzep36O6aG5Ss6SW6u8PAERERDZicyMiIsVhcyMiIsVhcyMiIsVhcyMiIsWRtbkVFRVhxowZaNmyJdzd3dGjRw8cPHhQzpSIiEgBZG1uEydORFpaGr788kscP34cAwcORP/+/ZGbmytnWkRE5OBka253797FN998g3fffRe9e/dGWFgYFixYgLCwMCQlJcmVFhERKYBsV/NVVFTAYDBUueDS3d0d+/btq/Y5er0eer3e9LiwsFDSHImUjPVESibbnptGo0H37t2xaNEi5OXlwWAwYO3atfj5559rHD4nMTERWq3WNHEcPCLbsZ5IyWT9ze3LL7+EEALNmjWDWq3Ghx9+iNGjR9c4nFNCQgIKCgpMU05OTh1nTKQcrCdSMlkHGQsNDUVmZiaKi4tRWFiIgIAAvPDCC2jVqlW166vVaqjV6jrOkkiZWE+kZPXiOjdPT08EBATg9u3bSE1NRUxMjNwpERGRA5N1zy01NRVCCISHhyM7Oxuvv/46WrdujXHjxsmZFhEROThZ99wKCgoQHx+P1q1bY+zYsejVqxdSU1Ph4uIiZ1pEROTgZN1zGzVqFEaNGiVnCkREpED14jc3IiIie2JzIyIixWFzIyIixWFzIyIixWFzIyIixZH1bEl7eeaZZ+Dq6ipJ7BdffFGSuJUVFRVJGv/ChQuSxv/94Nf21rZtW0njA8Dw4cMli11eXo4dO3ZIFt/eAgICJKunuvi/9Pf3lzT+1atXJY0v9aVQUr8/ANCuXTtJ4hoMBpw8edKidbnnRkREisPmRkREisPmRkREisPmRkREisPmRkREisPmRkREisPmRkREiiNrczMYDJg7dy5CQkLg7u6O0NBQLFq0CEIIOdMiIiIHJ+tF3O+88w6SkpKQkpKCiIgIHDp0COPGjYNWq8X06dPlTI2IiByYrM3tp59+QkxMDIYMGQIACA4Oxvr16/HLL7/ImRYRETk4WQ9L9ujRA3v27MHp06cBAFlZWdi3bx+io6OrXV+v16OwsNBsIiLbsJ5IyWTdc5szZw4KCwvRunVrODs7w2AwYPHixRgzZky16ycmJmLhwoV1nCWRMrGeSMlk3XPbuHEjvvrqK6xbtw6HDx9GSkoK3n//faSkpFS7fkJCAgoKCkxTTk5OHWdMpBysJ1IyWffcXn/9dcyZM8c08n779u1x8eJFJCYmIjY2tsr6arUaarW6rtMkUiTWEymZrHtuJSUlcHIyT8HZ2RlGo1GmjIiISAlk3XMbOnQoFi9ejBYtWiAiIgJHjhzBsmXLMH78eDnTIiIiBydrc/voo48wd+5cvPrqq8jPz0dgYCCmTJmCefPmyZkWERE5OFmbm0ajwYoVK7BixQo50yAiIoXh2JJERKQ4bG5ERKQ4bG5ERKQ4bG5ERKQ4bG5ERKQ4bG5ERKQ4sl4KYC/+/v5wc3OTJHbDhg0liVtZRUWFpPGlem/uCw8PlzS+i4uLpPEBICgoSLLYZWVlksWWQkBAgGTDcjVq1EiSuJWVlpZKGj83N1fS+G3atJE0/uOPPy5pfADo16+fJHH1ej1Onjxp0brccyMiIsVhcyMiIsVhcyMiIsVhcyMiIsVhcyMiIsVhcyMiIsWRvbnl5ubi5ZdfRqNGjeDu7o727dvj0KFDcqdFREQOTNbr3G7fvo2ePXuib9+++OGHH+Dn54czZ87Ax8dHzrSIiMjBydrc3nnnHQQFBWH16tWmeSEhITJmRERESiDrYcnt27fj8ccfx/PPP48mTZqgc+fO+PTTT+VMiYiIFEDW5nbu3DkkJSXhkUceQWpqKqZOnYrp06cjJSWl2vX1ej0KCwvNJiKyDeuJlMym5nbu3Dm7bNxoNCIyMhJLlixB586dMXnyZEyaNAkrV66sdv3ExERotVrTJOV4gERKx3oiJbOpuYWFhaFv375Yu3ZtrQYpDQgIQNu2bc3mtWnTBpcuXap2/YSEBBQUFJimnJwcm7dN9LBjPZGS2dTcDh8+jA4dOmDWrFnw9/fHlClT8Msvv1gdp2fPnjh16pTZvNOnT6Nly5bVrq9Wq+Ht7W02EZFtWE+kZDY1t06dOuGDDz5AXl4evvjiC1y5cgW9evVCu3btsGzZMly/ft2iODNnzsSBAwewZMkSZGdnY926dVi1ahXi4+NtSYuIiAhALU8oadCgAUaMGIFNmzbhnXfeQXZ2NmbPno2goCCMHTsWV65ceeDzu3Tpgq1bt2L9+vVo164dFi1ahBUrVmDMmDG1SYuIiB5ytWpuhw4dwquvvoqAgAAsW7YMs2fPxtmzZ5GWloa8vDzExMT8YYxnnnkGx48fR2lpKU6ePIlJkybVJiUiIiLbLuJetmwZVq9ejVOnTmHw4MFYs2YNBg8eDCene70yJCQEycnJCA4OtmeuREREFrGpuSUlJWH8+PGIi4tDQEBAtes0adIEn3/+ea2SIyIisoVNzS0tLQ0tWrQw7andJ4RATk4OWrRoAVdXV8TGxtolSSIiImvY9JtbaGgobty4UWX+rVu3ODYkERHJzqbmJoSodr5Op4Obm1utEiIiIqotqw5Lzpo1CwCgUqkwb948eHh4mJYZDAb85z//QadOneyaIBERkbWsam5HjhwBcG/P7fjx43B1dTUtc3V1RceOHTF79mz7ZmgBb29vyfYY3d3dJYlbWVlZmaTx79y5I2l8g8Egafzf/7YrhbCwMMli12aIOjkIIWo8OlNbdfFeSL2NixcvShq/cePGksav/HdbKhUVFZLEteZvjVXNLT09HQAwbtw4fPDBBxyuh4iI6iWbzpasfHNRIiKi+sbi5jZixAgkJyfD29sbI0aMeOC6W7ZsqXViREREtrK4uWm1WqhUKtO/iYiI6iuLm1vlQ5E8LElERPWZTaeh3b17FyUlJabHFy9exIoVK7Br1y67JUZERGQrm5pbTEwM1qxZA+DeaeZdu3bF0qVLERMTg6SkJIvjLFiwACqVymxq3bq1LSkRERGZ2Hwn7ieffBIAsHnzZvj7++PixYtYs2YNPvzwQ6tiRURE4MqVK6Zp3759tqREBF2ZDosyFyFoeRCc33JG0PIgLMpcBF2ZTu7UiByKrkyHxP2JCPskDJ7veiLskzAk7k90qFqy6VKAkpISaDQaAMCuXbswYsQIODk54YknnrD6AscGDRrA39/fljSITHRlOvRJ7oOjV4/CKIwAgMuFl7EgcwG2ndqGzLhMeLl6yZwlUf2nK9Nh4LqByMrPMtVSblEu/r7/7/jXmX9h10u7HKKWbNpzCwsLw7Zt25CTk4PU1FQMHDgQAJCfn2/1hd1nzpxBYGAgWrVqhTFjxuDSpUu2pEQPueU/LzdrbPcZhRFHrx7F8p+Xy5QZkWP56OBHZo3tPqMwIis/Cx8d/EimzKxjU3ObN28eZs+ejeDgYHTr1g3du3cHcG8vrnPnzhbH6datG5KTk7Fz504kJSXh/PnzePLJJ1FUVFTt+nq9HoWFhWYTEQCsOryqSjHeZxRGrDq8qo4zqv9YT1Sdz7M+f2AtfZ7lGPfptOmw5HPPPYdevXrhypUr6Nixo2l+v3798Oyzz1ocJzo62vTvDh06oFu3bmjZsiU2btyICRMmVFk/MTERCxcutCVlUri8orxaLX8YsZ6oOld0V2q1vL6weURaf39/dO7c2WxQ265du9bqbMeGDRvi0UcfRXZ2drXLExISUFBQYJpycnJs3hYpS6AmsFbLH0asJ6pOgFdArZbXFzY1t+LiYsydOxc9evRAWFgYWrVqZTbZSqfT4ezZswgIqP7NU6vV8Pb2NpuIAGBy5GQ4qar/ODupnDA5cnIdZ1T/sZ6oOhM6TnhgLU3oWPWoWn1k02HJiRMnIjMzE6+88goCAgJMw3JZa/bs2Rg6dChatmyJvLw8zJ8/H87Ozhg9erRN8ejhNbP7TGw7ta3KSSVOKid08u+Emd1nypgdkeN4rctr+NeZf1U5qcRJ5YSOTTritS6vyZid5Wxqbj/88AN27NiBnj171mrjly9fxujRo3Hz5k34+fmhV69eOHDgAPz8/GoVlx4+Xq5eyIzLxPKfl2PV4VXIK8pDoCYQkyMnY2b3mQ5x6jJRfeDl6oVdL+3CRwc/wudZn+OK7goCvAIwoeMEvNblNYepJZuam4+PD3x9fWu98Q0bNtQ6BtF9Xq5emNtnLub2mSt3KkQOzcvVCwk9E5DQM0HuVGxm029uixYtwrx588zGlyQiIqovbNpzW7p0Kc6ePYumTZsiODgYLi4uZssPHz5sl+SIiIhsYVNzGz58uJ3TICIish+bmtv8+fPtnQcREZHd2HwR9507d/DZZ58hISEBt27dAnDvcGRubq7dkiMiIrKFTXtux44dQ//+/aHVanHhwgVMmjQJvr6+2LJlCy5dumS61xsREZEcbGpus2bNQlxcHN59913TrW8AYPDgwXjppZfslpylLly4ALVaXefbtZcGDWz6b7CYTiftPZgMBoOk8Y3G6gdxtaeahnyzh7KyMsliS+HOnTtwdXWVJPbly5cliVuX27h27Zqk8U+dOiVpfKnzB6T7m2bN3xqbDksePHgQU6ZMqTK/WbNmuHr1qi0hiYiI7Mam5qZWq6u9Pcbp06c5uggREcnOpuY2bNgwvPXWWygvLwcAqFQqXLp0CW+88QZGjhxp1wSJiIisZVNzW7p0KXQ6Hfz8/HD37l306dMHYWFh0Gg0WLx4sb1zJCIisopNv/pptVqkpaVh//79yMrKgk6nQ2RkJPr372/v/IiIiKxmdXMzGo1ITk7Gli1bcOHCBahUKoSEhMDf3x9CCJtvf0NERGQvVh2WFEJg2LBhmDhxInJzc9G+fXtERETg4sWLiIuLw7PPPitVnkRERBazqrklJydj79692LNnD44cOYL169djw4YNyMrKwu7du/Hjjz/afAH322+/DZVKhRkzZtj0fCIiovusam7r16/HX//6V/Tt27fKsqeffhpz5szBV199ZXUSBw8exD//+U906NDB6ucSERH9nlXN7dixYxg0aFCNy6Ojo5GVlWVVAjqdDmPGjMGnn34KHx8fq55LRERUHaua261bt9C0adMalzdt2hS3b9+2KoH4+HgMGTLEojMt9Xo9CgsLzSYisg3riZTMquZmMBgeOGaYs7MzKioqLI63YcMGHD58GImJiRatn5iYCK1Wa5qCgoIs3hYRmWM9kZJZdSmAEAJxcXE1DlKs1+stjpWTk4M///nPSEtLg5ubm0XPSUhIwKxZs0yPCwsLWZBENmI9kZJZ1dxiY2P/cJ2xY8daFOvXX39Ffn4+IiMjTfMMBgP27t2Ljz/+GHq9Hs7OzmbPUavVDj36P1F9wnoiJbOqua1evdpuG+7Xrx+OHz9uNm/cuHFo3bo13njjjSqNjYiIyFLS3kjsATQaDdq1a2c2z9PTE40aNaoyn4iIyBo2DZxMRERUn8m251adjIwMuVMgIiIF4J4bEREpDpsbEREpDpsbEREpDpsbEREpDpsbEREpTr06W9JWWq1WspEWjEajJHErs2Y8TlvcvXtX0vjWDLtmi7r4P2jUqJFksUtLSyWLLYW8vDy4uLhIEjs4OFiSuJVduXJF0vhRUVGSxo+JiZE0fmBgoKTxAaBNmzaSxLXmbxn33IiISHHY3IiISHHY3IiISHHY3IiISHHY3IiISHHY3IiISHHY3IiISHFkbW6JiYno0qULNBoNmjRpguHDh+PUqVNypkRERAoga3PLzMxEfHw8Dhw4gLS0NJSXl2PgwIEoLi6WMy0iInJwso5QsnPnTrPHycnJaNKkCX799Vf07t1bpqyIiMjR1avhtwoKCgAAvr6+1S7X6/VmQz0VFhbWSV5ESsR6IiWrNyeUGI1GzJgxAz179kS7du2qXScxMRFardY0BQUF1XGWRMrBeiIlqzfNLT4+HidOnMCGDRtqXCchIQEFBQWmKScnpw4zJFIW1hMpWb04LDlt2jR899132Lt3L5o3b17jemq1WrLR/4keNqwnUjJZm5sQAq+99hq2bt2KjIwMhISEyJkOEREphKzNLT4+HuvWrcO3334LjUaDq1evArh3fzZ3d3c5UyMiIgcm629uSUlJKCgowFNPPYWAgADT9PXXX8uZFhEROTjZD0sSERHZW705W5KIiMhe2NyIiEhx2NyIiEhx2NyIiEhx2NyIiEhx2NyIiEhx6sXwW7XVvn17eHh4SBK78qjpUrlz546k8e/evStp/Pt3c5BKXfwf1DRYtz2UlJRIFlsK48aNk6yeAgMDJYlbWbdu3SSNv3//fknjX7t2TdL4bm5uksYHgEaNGkkS15pa4p4bEREpDpsbEREpDpsbEREpDpsbEREpDpsbEREpDpsbEREpjqzNbe/evRg6dCgCAwOhUqmwbds2OdMhIiKFkLW5FRcXo2PHjvjHP/4hZxpERKQwsl7EHR0djejoaDlTICIiBeJvbkREpDgONfyWXq83G4qpsLBQxmyIHBvriZTMofbcEhMTodVqTVNQUJDcKRE5LNYTKZlDNbeEhAQUFBSYppycHLlTInJYrCdSMoc6LKlWq6FWq+VOg0gRWE+kZLI2N51Oh+zsbNPj8+fP4+jRo/D19UWLFi1kzIyIiByZrM3t0KFD6Nu3r+nxrFmzAACxsbFITk6WKSsiInJ0sja3p556CkIIOVMgIiIFcqgTSoiIiCzB5kZERIrD5kZERIrD5kZERIrD5kZERIrjUBdx/979My1LSkok20ZdjLen0+kkjV9aWippfKnzr4v/Ayk/Q3fv3gWAen9mMOvJMvf/P6XCeqqZNbWkEvW94h7g8uXLHA+PHEZOTg6aN28udxo1Yj2Ro7Cklhy6uRmNRuTl5UGj0UClUv3h+oWFhQgKCkJOTg68vb0lyUnqbTC+/NuwNr4QAkVFRQgMDISTU/39JaC+1RPjy7+N+hbfmlpy6MOSTk5ONn0T9vb2luzDVlfbYHz5t2FNfK1WK1ke9lJf64nx5d9GfYpvaS3V36+RRERENmJzIyIixXmomptarcb8+fMlvc2H1NtgfPm3URevwRE4+vvs6PHrYhuOHN+hTyghIiKqzkO150ZERA8HNjciIlIcNjciIlIcNjciIlKch6q56XQ6TJs2Dc2bN4e7uzvatm2LlStX2i1+XFwcVCqV2TRo0CC7xS8vL8cbb7yB9u3bw9PTE4GBgRg7dizy8vLsto1//OMfCA4OhpubG7p164ZffvnFLnETExPRpUsXaDQaNGnSBMOHD8epU6fsEhsAFixYUOW9b926td3iGwwGzJ07FyEhIXB3d0doaCgWLVpU78eLlIrUtQRIW0+spQdTRD2Jh8ikSZNEaGioSE9PF+fPnxf//Oc/hbOzs/j222/tEj82NlYMGjRIXLlyxTTdunXLLrGFEOLOnTuif//+4uuvvxb/+9//xM8//yy6du0qHnvsMbvE37Bhg3B1dRVffPGF+O2338SkSZNEw4YNxbVr12odOyoqSqxevVqcOHFCHD16VAwePFi0aNFC6HQ6O2QuxPz580VERITZe3/9+nW7xBZCiMWLF4tGjRqJ7777Tpw/f15s2rRJeHl5iQ8++MBu23AkUteSENLWE2vpwZRQTw9Vc4uIiBBvvfWW2bzIyEjxt7/9zS7xY2NjRUxMjF1iWeqXX34RAMTFixdrHatr164iPj7e9NhgMIjAwECRmJhY69i/l5+fLwCIzMxMu8SbP3++6Nixo11iVWfIkCFi/PjxZvNGjBghxowZI9k26zOpa0mIuq8n1tL/p4R6eqgOS/bo0QPbt29Hbm4uhBBIT0/H6dOnMXDgQLttIyMjA02aNEF4eDimTp2Kmzdv2i12dQoKCqBSqdCwYcNaxSkrK8Ovv/6K/v37m+Y5OTmhf//++Pnnn2uZZVUFBQUAAF9fX7vFPHPmDAIDA9GqVSuMGTMGly5dslvsHj16YM+ePTh9+jQAICsrC/v27UN0dLTdtuFI6qKWgLqtJ9aSOYevJ7u1SQdQWloqxo4dKwCIBg0aCFdXV5GSkmK3+OvXrxfffvutOHbsmNi6dato06aN6NKli6ioqLDbNiq7e/euiIyMFC+99FKtY+Xm5goA4qeffjKb//rrr4uuXbvWOn5lBoNBDBkyRPTs2dNuMb///nuxceNGkZWVJXbu3Cm6d+8uWrRoIQoLC+0S32AwiDfeeEOoVCrRoEEDoVKpxJIlS+wS2xFJXUtC1G09sZbMKaGeFNvc1q5dKzw9PU3T3r17xXvvvSceffRRsX37dpGVlSU++ugj4eXlJdLS0uwS//fOnj0rAIjdu3fb7TXcV1ZWJoYOHSo6d+4sCgoKbIpfWV0W5J/+9CfRsmVLkZOTY9e4ld2+fVt4e3uLzz77zC7x1q9fL5o3by7Wr18vjh07JtasWSN8fX1FcnKyXeLXZ1LXUk3b+L3a1BNrqXYcsZ4U29wKCwvFmTNnTFNJSYlwcXER3333ndl6EyZMEFFRUXaJX53GjRuLlStX2u01CHGvGIcPHy46dOggbty4YVPs39Pr9cLZ2Vls3brVbP7YsWPFsGHD7LINIYSIj48XzZs3F+fOnbNbzJo8/vjjYs6cOXaJ1bx5c/Hxxx+bzVu0aJEIDw+3S/z6TOpaqmkb1bG1nlhLtedo9eTQ93N7EI1GA41GY3pcWFiI8vLyKje4c3Z2htForHX86ly+fBk3b95EQECA1fFr2kZ5eTlGjRqFM2fOID09HY0aNbIp9u+5urrisccew549ezB8+HAA925euWfPHkybNq3W8YUQeO2117B161ZkZGQgJCSk1jEfRKfT4ezZs3jllVfsEq+kpMRunx1HI3UtVbeN6tSmnlhLteOQ9WS3NukA+vTpIyIiIkR6ero4d+6cWL16tXBzcxOffPJJrWMXFRWJ2bNni59//lmcP39e7N69W0RGRopHHnlElJaW2iH7e98yhw0bJpo3by6OHj1qdpquXq+vdfwNGzYItVotkpOTxX//+18xefJk0bBhQ3H16tVax546darQarUiIyPDLO+avqFb6//+7/9ERkaGOH/+vNi/f7/o37+/aNy4scjPz7dL/NjYWNGsWTPTqctbtmwRjRs3Fn/5y1/sEt/RSFlLQkhfT6ylB1NCPT1Uze3KlSsiLi5OBAYGCjc3NxEeHi6WLl0qjEZjrWOXlJSIgQMHCj8/P+Hi4iJatmwpJk2aZJcP833nz58XAKqd0tPT7bKNjz76SLRo0UK4urqKrl27igMHDtglbk15r1692i7xX3jhBREQECBcXV1Fs2bNxAsvvCCys7PtEluIe4e1/vznP4sWLVoINzc30apVK/G3v/3NLn8IHZGUtSSE9PXEWnowJdQTb3lDRESK81Bd50ZERA8HNjciIlIcNjciIlIcNjciIlIcNjciIlIcNjciIlIcNjciIlIcNreHjEqlwrZt2+ROg8jhsZbqNzY3hbl+/TqmTp2KFi1aQK1Ww9/fH1FRUdi/f7/cqRE5FNaSY1PswMkPq5EjR6KsrAwpKSlo1aoVrl27hj179kh+01QipWEtOTi7DeRFsrt9+7YAIDIyMmpcB4DZrTiOHTsm+vbtK9zc3ISvr6+YNGmSKCoqMi2PjY0VMTExYsGCBaJx48ZCo9GIKVOmmI0BZzAYxJIlS0RwcLBwc3MTHTp0EJs2bZLkNRLVBdaS4+NhSQXx8vKCl5cXtm3bBr1e/4frFxcXIyoqCj4+Pjh48CA2bdqE3bt3V7ktx549e3Dy5ElkZGRg/fr12LJlCxYuXGhanpiYiDVr1mDlypX47bffMHPmTLz88svIzMy0+2skqgusJQWQu7uSfW3evFn4+PgINzc30aNHD5GQkCCysrJMy1Hp2+aqVauEj4+P0Ol0puU7duwQTk5OptHXY2Njha+vryguLjatk5SUJLy8vITBYBClpaXCw8Ojyl2HJ0yYIEaPHi3hKyWSFmvJsXHPTWFGjhyJvLw8bN++HYMGDUJGRgYiIyORnJxcZd2TJ0+iY8eO8PT0NM3r2bMnjEYjTp06ZZrXsWNHeHh4mB53794dOp0OOTk5yM7ORklJCQYMGGD6tuvl5YU1a9bg7Nmzkr5WIimxlhwbTyhRIDc3NwwYMAADBgzA3LlzMXHiRMyfPx9xcXF235ZOpwMA7NixA82aNTNbplar7b49orrEWnJc3HN7CLRt2xbFxcVV5rdp0wZZWVlmy/bv3w8nJyeEh4eb5mVlZeHu3bumxwcOHICXlxeCgoLQtm1bqNVqXLp0CWFhYWZTUFCQtC+MqI6xlhwHm5uC3Lx5E08//TTWrl2LY8eO4fz589i0aRPeffddxMTEVFl/zJgxcHNzQ2xsLE6cOIH09HS89tpreOWVV9C0aVPTemVlZZgwYQL++9//4vvvv8f8+fMxbdo0ODk5QaPRYPbs2Zg5cyZSUlJw9uxZHD58GB999BFSUlLq8uUT2Q1rSQHk/tGP7Ke0tFTMmTNHREZGCq1WKzw8PER4eLh48803RUlJiRDC9tOX582bJxo1aiS8vLzEpEmTRGlpqWkdo9EoVqxYIcLDw4WLi4vw8/MTUVFRIjMzs85eO5E9sZYcn0oIIeRusFR/xcXF4c6dOxxmiKiWWEt1i4cliYhIcdjciIhIcXhYkoiIFId7bkREpDhsbkREpDhsbkREpDhsbkREpDhsbkREpDhsbkREpDhsbkREpDhsbkREpDhsbkREpDj/Dx8hCzg0n1ZYAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 500x750 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "side_length = 2.5\n",
    "fig, axes = plt.subplots(3, 2, figsize=(2*side_length, 3*side_length), sharex=True, sharey=True)\n",
    "\n",
    "for distance_id, distance_key in enumerate(distances.keys()):\n",
    "    for seed in [0, 1]:\n",
    "        ax = axes[distance_id, seed]\n",
    "        df = out[(distance_key, seed)]\n",
    "\n",
    "        distance_field = np.zeros((n_densities, n_slopes))\n",
    "        for density_idx in range(n_densities):\n",
    "            for slope_idx in range(n_slopes):\n",
    "                distance_field[density_idx, slope_idx] = df[\n",
    "                    (df['density_idx'] == density_idx) &\n",
    "                    (df['slope_idx'] == slope_idx)\n",
    "                ]['distance'].item()\n",
    "\n",
    "        ax.imshow(distance_field, cmap='gray_r')\n",
    "        ax.quiver(\n",
    "            df['slope_idx'], df['density_idx'],\n",
    "            df['slope_grad'], df['density_grad'],\n",
    "            color='blue'\n",
    "        )\n",
    "        # show green dot at the center\n",
    "        ax.scatter([5], [3], color='green', marker='o', s=30)\n",
    "        ax.set_xlabel('Slope')\n",
    "        ax.set_xticks(np.arange(n_slopes))\n",
    "        ax.set_xticklabels([\n",
    "            \"{:.0f}\".format(10*xt) for xt in\n",
    "            np.linspace(-1, 1, n_slopes+2)[1:-1]])\n",
    "        ax.set_yticks(np.arange(n_slopes))\n",
    "        ax.set_yticklabels([\n",
    "            \"{:.0f}\".format(10*xt) for xt in\n",
    "            np.linspace(0, 1, n_densities+2)[1:-1]])\n",
    "        if seed==0:\n",
    "            ax.set_ylabel('Density')\n",
    "        ax.invert_yaxis()\n",
    "\n",
    "        if distance_id==0:\n",
    "            if seed==0:\n",
    "                ax.set_title('Microstructure')\n",
    "            else:\n",
    "                ax.set_title('Mesostructure')\n",
    "\n",
    "    #plt.savefig('gradient_field.png', dpi=300, bbox_inches='tight')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 228,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "49it [00:00, 228.03it/s]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fe6c8e37b80>"
      ]
     },
     "execution_count": 228,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOcAAADsCAYAAACPM7HjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAANvElEQVR4nO3dTWxU9f7H8c9Mh05BppUHW56q1IQESaU8lBJSrg+hgQUs2BgXkBASNSFTbWVDGqPEEBhcaBBCeDKCiWBZEZEETFOhuSQ0xfZiaIiIceHw0BbUTFuUae2c/+J/O/f28tSZTs/50nm/krPoYcbfd6jvnmE6Z47PcRxHAMzxez0AgAcjTsAo4gSMIk7AKOIEjCJOwCjiBIwiTsCogNsLJhIJ3bx5U6FQSD6fz+3lAU85jqOenh7NmDFDfv+jj42ux3nz5k0VFxe7vSxgSjQa1axZsx55G9fjDIVCkqSZM2c+9idHpixdutSVdSSpsLDQtbUkady4ca6t1d/f79paktTV1eXaWi0tLa6sk0gkdP369WQHj+J6nINPZf1+v2txuvk/cDAYdG0tyd3H5tb3a9BYfmzD+ScdLwgBRhEnYBRxAkYRJ2AUcQJGESdgVFpx7t27V7Nnz1ZeXp6WLl3q2u+IgGyScpzHjx/X5s2btXXrVrW1tamsrEyrVq1y9RfGQDZIOc5PPvlEb775pjZu3Kh58+Zp//79mjBhgj7//PPRmA/IWinF2dfXp9bWVlVVVf3nP+D3q6qqShcuXMj4cEA2S+nte3fu3NHAwICKioqG7C8qKtKPP/74wPvE43HF4/Hk193d3WmMCWSfUX+1NhKJqKCgILlxRgowPCnFOXXqVOXk5Kizs3PI/s7OTk2bNu2B96mrq1MsFktu0Wg0/WmBLJJSnLm5uVq8eLEaGxuT+xKJhBobG7Vs2bIH3icYDCo/P3/IBuDxUj5lbPPmzdqwYYPKy8tVUVGhXbt26e7du9q4ceNozAdkrZTjfP3113X79m198MEH6ujo0IIFC3TmzJn7XiQCMDJpnWxdXV2t6urqTM8C4L/w3lrAKOIEjCJOwCjiBIwiTsAo4gSMIk7AKNc/VHqQz+dz7Vops2fPdmUdSY/9iP1MGz9+vGtrufkhz5LU0dHh2lptbW2urDMwMDDs23LkBIwiTsAo4gSMIk7AKOIEjCJOwCjiBIwiTsAo4gSMIk7AKOIEjCJOwCjiBIwiTsAo4gSMIk7AKOIEjCJOwCjiBIwiTsAo4gSMIk7AKOIEjCJOwCjiBIwiTsAozy7HkJ+fr5ycHFfWev75511ZR5IWLFjg2lqS1NfX59pabl+Owc1LTbi1FpdjAMYA4gSMIk7AKOIEjCJOwCjiBIwiTsAo4gSMIk7AqJTijEQiWrJkiUKhkAoLC7V27VpdvXp1tGYDslpKcTY1NSkcDqu5uVkNDQ3q7+/XypUrdffu3dGaD8haKb239syZM0O+PnLkiAoLC9Xa2qqXXnopo4MB2W5E/+aMxWKSpMmTJ2dkGAD/kfZZKYlEQrW1taqsrFRpaelDbxePxxWPx5Nfd3d3p7skkFXSPnKGw2G1t7ervr7+kbeLRCIqKChIbsXFxekuCWSVtOKsrq7WqVOndPbsWc2aNeuRt62rq1MsFktu0Wg0rUGBbJPS01rHcfT222/rxIkTOnfunEpKSh57n2AwqGAwmPaAQLZKKc5wOKxjx47p66+/VigUUkdHhySpoKDA1bPWgWyQ0tPaffv2KRaL6ZVXXtH06dOT2/Hjx0drPiBrpfy0FoA7eG8tYBRxAkYRJ2AUcQJGESdgFHECRhEnYJRn10qJx+Nj8lop5eXlrq0lSf39/a6t9ffff7u2liTdu3dvzK3FtVKAMYA4AaOIEzCKOAGjiBMwijgBo4gTMIo4AaOIEzCKOLNEb1+vtp/frpI9JcqL5KlkT4m2n9+u3r5er0fDQ3j29j24p7evVyu+XKEfOn9QwklIkm703NC2f27TyZ9OqnF9oybmTvR4SvwvjpxZ4NOWT4eEOSjhJPRD5w/6tOVTjybDoxBnFvjsX5/dF+aghJPQZ//6zOWJMBzEmQVu9d4a0Z/DG8SZBaZPnD6iP4c3iDMLvLHwDfl9D/5W+31+vbHwDZcnwnAQZxaoqahRWVHZfYH6fX6VFZWppqLGo8nwKMSZBSbmTlTj+ka9/4/3NTM0U36fXzNDM/X+P97n1yiG8XvOLDExd6LeW/6e3lv+ntejYJg4cgJGESdgFHECRhEnYBRxAkYRJ2AUcQJGefZ7zlgsJr/fnZ8Nt2/fdmUdSWpvb3dtLUlKJB58tsmTvpYk3bhxw7W1envdOek8lb9DjpyAUcQJGEWcgFHECRhFnIBRxAkYRZyAUcQJGEWcgFEjinPnzp3y+Xyqra3N0DgABqUd58WLF3XgwAHNnz8/k/MA+Le04uzt7dW6det06NAhTZo0KdMzAVCacYbDYa1evVpVVVWZngfAv6V8Vkp9fb3a2tp08eLFYd0+Ho8rHo8nv+7u7k51SSArpXTkjEajqqmp0dGjR5WXlzes+0QiERUUFCS34uLitAYFsk1Kcba2tqqrq0uLFi1SIBBQIBBQU1OTdu/erUAgoIGBgfvuU1dXp1gsltyi0WjGhgfGspSe1q5YsUKXL18esm/jxo2aO3eutmzZopycnPvuEwwGFQwGRzYlkIVSijMUCqm0tHTIvqeeekpTpky5bz+AkeEdQoBRI/4MoXPnzmVgDAD/iyMnYBRxAkYRJ2AUcQJGESdgFHECRhEnYJRn10qJx+Py+XxeLT9qfv/9d1fX++uvv1xby61r2wxy6/olkvuPbTjsTQRAEnECZhEnYBRxAkYRJ2AUcQJGESdgFHECRhEnYBRxAkYRJ2AUcQJGESdgFHECRhEnYBRxAkYRJ2AUcQJGESdgFHECRhEnYBRxAkYRJ2AUcQJGESdgFHECRnl2OYZYLObaWrdu3XJtrc7OTtfWkqR79+65ttaECRNcW0ty9+/yt99+c2Udx3GGfVuOnIBRxAkYRZyAUcQJGEWcgFHECRhFnIBRxAkYRZyAUSnHeePGDa1fv15TpkzR+PHj9eKLL+r7778fjdmArJbS2/f++OMPVVZW6tVXX9Xp06f1zDPP6Nq1a5o0adJozQdkrZTi/Oijj1RcXKzDhw8n95WUlGR8KAApPq09efKkysvL9dprr6mwsFALFy7UoUOHRms2IKulFOcvv/yiffv2ac6cOfr222+1adMmvfPOO/riiy8eep94PK7u7u4hG4DHS+lpbSKRUHl5uXbs2CFJWrhwodrb27V//35t2LDhgfeJRCL68MMPRz4pkGVSOnJOnz5d8+bNG7LvhRde0K+//vrQ+9TV1SkWiyW3aDSa3qRAlknpyFlZWamrV68O2ffTTz/pueeee+h9gsGggsFgetMBWSylI+e7776r5uZm7dixQz///LOOHTumgwcPKhwOj9Z8QNZKKc4lS5boxIkT+uqrr1RaWqpt27Zp165dWrdu3WjNB2StlD9DaM2aNVqzZs1ozALgv/DeWsAo4gSMIk7AKOIEjCJOwCjiBIxy/XIMqXwcfaa4eckCv9/dn3duPjafz+faWtL/nzThFrf+vxxcZzjr+RyXa7l+/bqKi4vdXBIwJxqNatasWY+8jetxJhIJ3bx5U6FQaNg/ibu7u1VcXKxoNKr8/PxRntBdPLYnU7qPzXEc9fT0aMaMGY99luX601q/3//YnxgPk5+fP+a+yYN4bE+mdB5bQUHBsG7HC0KAUcQJGPVExBkMBrV169YxeV4oj+3J5MZjc/0FIQDD80QcOYFsRJyAUcQJGEWcgFFPRJx79+7V7NmzlZeXp6VLl6qlpcXrkUYsEoloyZIlCoVCKiws1Nq1a+/7ZMOxYufOnfL5fKqtrfV6lIxw62Je5uM8fvy4Nm/erK1bt6qtrU1lZWVatWqVurq6vB5tRJqamhQOh9Xc3KyGhgb19/dr5cqVunv3rtejZdTFixd14MABzZ8/3+tRMmLwYl7jxo3T6dOndeXKFX388cejczEvx7iKigonHA4nvx4YGHBmzJjhRCIRD6fKvK6uLkeS09TU5PUoGdPT0+PMmTPHaWhocF5++WWnpqbG65FGbMuWLc7y5ctdWcv0kbOvr0+tra2qqqpK7vP7/aqqqtKFCxc8nCzzYrGYJGny5MkeT5I54XBYq1evHvL9e9K5eTEv03HeuXNHAwMDKioqGrK/qKhIHR0dHk2VeYlEQrW1taqsrFRpaanX42REfX292traFIlEvB4lo9K5mFe6XD8rBfcLh8Nqb2/X+fPnvR4lI6LRqGpqatTQ0KC8vDyvx8modC7mlS7TR86pU6cqJydHnZ2dQ/Z3dnZq2rRpHk2VWdXV1Tp16pTOnj2b9ql01rS2tqqrq0uLFi1SIBBQIBBQU1OTdu/erUAgoIGBAa9HTFs6F/NKl+k4c3NztXjxYjU2Nib3JRIJNTY2atmyZR5ONnKO46i6ulonTpzQd999N6auEL5ixQpdvnxZly5dSm7l5eVat26dLl26pJycHK9HTFs6F/NKmysvO41AfX29EwwGnSNHjjhXrlxx3nrrLefpp592Ojo6vB5tRDZt2uQUFBQ4586dc27dupXc/vzzT69HGxVj5dXalpYWJxAIONu3b3euXbvmHD161JkwYYLz5ZdfZnwt83E6juPs2bPHefbZZ53c3FynoqLCaW5u9nqkEZP0wO3w4cNejzYqxkqcjuM433zzjVNaWuoEg0Fn7ty5zsGDB0dlHU4ZA4wy/W9OIJsRJ2AUcQJGESdgFHECRhEnYBRxAkYRJ2AUcQJGESdgFHECRhEnYNT/AYX15P4aZJ9zAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 250x250 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "offset = 0\n",
    "exponent = 1\n",
    "Sy = SX[('jtfs', 0, 0.5, 0.)][offset:]**exponent\n",
    "D = np.zeros((n_densities, n_slopes))\n",
    "seed = 1\n",
    "\n",
    "for _, row in tqdm.tqdm(dataset.df.iterrows()):\n",
    "    Sm = 0.5 * (SX[('jtfs', seed, row[\"density\"], row[\"slope\"])][offset:] + Sy)\n",
    "    Sx = SX[('jtfs', seed, row[\"density\"], row[\"slope\"])][offset:]**exponent\n",
    "    xm_dist = torchnmf.metrics.kl_div(Sx/Sx.sum(), Sm/Sm.sum())\n",
    "    ym_dist = torchnmf.metrics.kl_div(Sy/Sy.sum(), Sm/Sm.sum())\n",
    "    xy_dist = torchnmf.metrics.kl_div(Sx/Sx.sum(), Sy/Sy.sum())\n",
    "    js_dist = 0.5 * (xm_dist + ym_dist)\n",
    "\n",
    "    neuc_dist = torchnmf.metrics.euclidean(Sx/Sx.sum(), Sy/Sy.sum())\n",
    "    euc_dist = torchnmf.metrics.euclidean(Sx, Sy)\n",
    "    is_dist = torchnmf.metrics.is_div(Sx, Sy)\n",
    "    logeuc_dist = torchnmf.metrics.euclidean(torch.log1p(Sx), torch.log1p(Sy))\n",
    "    l1_dist = torch.linalg.vector_norm(Sx-Sy, ord=1)\n",
    "    D[int(row['density_idx']), int(row['slope_idx'])] = euc_dist.item()\n",
    "\n",
    "\n",
    "side_length = 2.5\n",
    "fig, ax = plt.subplots(1, 1, figsize=(side_length, side_length), sharex=True, sharey=True)\n",
    "ax.imshow(np.log1p(D), cmap='gray_r')\n",
    "ax.scatter([3], [3], color='green', marker='o', s=30)\n"
   ]
  }
 ],
 "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
}