hierarchy_plot.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6c3d6c21-3247-4e49-b74e-c27eb602ae11",
   "metadata": {},
   "outputs": [],
   "source": [
    "import argparse\n",
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import random \n",
    "import networkx as nx\n",
    "\n",
    "from tensorflow import keras\n",
    "from tqdm import tqdm\n",
    "from pathlib import Path\n",
    "\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "from tensorflow.keras.utils import pad_sequences\n",
    "from dataloader import LSSTSourceDataSet, load, get_augmented_data, get_static_features, ts_length, get_ts_upto_days_since_trigger, ts_flag_value, static_flag_value, augment_ts_length_to_days_since_trigger\n",
    "from loss import WHXE_Loss\n",
    "from taxonomy import get_taxonomy_tree, source_node_label\n",
    "from vizualizations import make_gif, plot_reliability_diagram, plot_data_set_composition, plot_day_vs_class_score, plot_lc, make_z_plots\n",
    "from interpret_results import get_conditional_probabilites, save_all_cf_and_rocs, save_leaf_cf_and_rocs, save_all_phase_vs_accuracy_plot\n",
    "from train_RNN import default_batch_size\n",
    "\n",
    "from LSST_Source import LSST_Source"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 313,
   "id": "2ddb8331-5270-4e87-b837-13282ab29050",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading testing data from disc...\n"
     ]
    }
   ],
   "source": [
    "# This step takes a while because it has load from disc to memory...\n",
    "print(\"Loading testing data from disc...\")\n",
    "test_dir = Path(\"processed/test\")\n",
    "max_class_count = 12000\n",
    "\n",
    "X_ts = load(f\"{test_dir}/x_ts.pkl\")\n",
    "X_static = load(f\"{test_dir}/x_static.pkl\")\n",
    "Y = load(f\"{test_dir}/y.pkl\")\n",
    "astrophysical_classes = load(f\"{test_dir}/a_labels.pkl\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 314,
   "id": "b4e54d03-4fdd-4185-b82e-3ac0f041e1a2",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 463528/463528 [00:01<00:00, 255836.98it/s]\n"
     ]
    }
   ],
   "source": [
    "for i in tqdm(range(len(X_static))):        \n",
    "    X_static[i] = get_static_features(X_static[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 512,
   "id": "d667f901-415a-41c9-9fa1-3faea1471fcb",
   "metadata": {},
   "outputs": [],
   "source": [
    "rare_class = \"TDE\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 513,
   "id": "7ed94f50-194c-4192-a04d-0cd8445980c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_ts_balanced = []\n",
    "X_static_balanced = []\n",
    "Y_balanced = []\n",
    "astrophysical_classes_balanced = []\n",
    "\n",
    "for c in np.unique(astrophysical_classes):\n",
    "\n",
    "    if c == rare_class:\n",
    "        idx = list(np.where(np.array(astrophysical_classes) == c)[0])\n",
    "        \n",
    "        if len(idx) > max_class_count:\n",
    "            idx = random.sample(idx, max_class_count)\n",
    "    \n",
    "        X_ts_balanced += [X_ts[i] for i in idx]\n",
    "        X_static_balanced += [X_static[i] for i in idx]\n",
    "        Y_balanced += [Y[i] for i in idx]\n",
    "        astrophysical_classes_balanced += [astrophysical_classes[i] for i in idx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 514,
   "id": "591b91af-cbbb-4cf0-b8ed-fb35a2499cf5",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_dir =  \"models/lsst_alpha_0.5/\"\n",
    "tree = get_taxonomy_tree()\n",
    "model = keras.models.load_model(f\"{model_dir}/best_model.h5\", compile=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 515,
   "id": "8186aa43-b3cd-4694-ae92-7fcaad01b596",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0,  1,  5, 18])"
      ]
     },
     "execution_count": 515,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "true_idx = np.where(Y_balanced[0]==1)[0]\n",
    "true_idx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 516,
   "id": "5c21c9c8-a1de-41fa-ae4a-f243f25b9964",
   "metadata": {},
   "outputs": [],
   "source": [
    "N = len(X_static_balanced) # Number of samples\n",
    "D = 50 # Number of days\n",
    "days = np.arange(1, D, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 517,
   "id": "9b0d5dcc-a4ac-4f77-bf56-c2001de52949",
   "metadata": {},
   "outputs": [],
   "source": [
    "L1 = np.zeros((N, D))\n",
    "L2 = np.zeros((N, D))\n",
    "L3 = np.zeros((N, D))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 518,
   "id": "4e0f19ed-486d-4688-ab0b-e212eb60d279",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1/1 [==============================] - 12s 12s/step\n",
      "Day 1\n",
      "1/1 [==============================] - 10s 10s/step\n",
      "Day 2\n",
      "1/1 [==============================] - 10s 10s/step\n",
      "Day 3\n",
      "1/1 [==============================] - 10s 10s/step\n",
      "Day 4\n",
      "1/1 [==============================] - 9s 9s/step\n",
      "Day 5\n",
      "1/1 [==============================] - 9s 9s/step\n",
      "Day 6\n",
      "1/1 [==============================] - 9s 9s/step\n",
      "Day 7\n",
      "1/1 [==============================] - 9s 9s/step\n",
      "Day 8\n",
      "1/1 [==============================] - 9s 9s/step\n",
      "Day 9\n",
      "1/1 [==============================] - 9s 9s/step\n",
      "Day 10\n",
      "1/1 [==============================] - 9s 9s/step\n",
      "Day 11\n",
      "1/1 [==============================] - 9s 9s/step\n",
      "Day 12\n",
      "1/1 [==============================] - 10s 10s/step\n",
      "Day 13\n",
      "1/1 [==============================] - 9s 9s/step\n",
      "Day 14\n",
      "1/1 [==============================] - 9s 9s/step\n",
      "Day 15\n",
      "1/1 [==============================] - 9s 9s/step\n",
      "Day 16\n",
      "1/1 [==============================] - 8s 8s/step\n",
      "Day 17\n",
      "1/1 [==============================] - 8s 8s/step\n",
      "Day 18\n",
      "1/1 [==============================] - 8s 8s/step\n",
      "Day 19\n",
      "1/1 [==============================] - 7s 7s/step\n",
      "Day 20\n",
      "1/1 [==============================] - 7s 7s/step\n",
      "Day 21\n",
      "1/1 [==============================] - 7s 7s/step\n",
      "Day 22\n",
      "1/1 [==============================] - 7s 7s/step\n",
      "Day 23\n",
      "1/1 [==============================] - 7s 7s/step\n",
      "Day 24\n",
      "1/1 [==============================] - 7s 7s/step\n",
      "Day 25\n",
      "1/1 [==============================] - 7s 7s/step\n",
      "Day 26\n",
      "1/1 [==============================] - 7s 7s/step\n",
      "Day 27\n",
      "1/1 [==============================] - 7s 7s/step\n",
      "Day 28\n",
      "1/1 [==============================] - 7s 7s/step\n",
      "Day 29\n",
      "1/1 [==============================] - 9s 9s/step\n",
      "Day 30\n",
      "1/1 [==============================] - 8s 8s/step\n",
      "Day 31\n",
      "1/1 [==============================] - 7s 7s/step\n",
      "Day 32\n",
      "1/1 [==============================] - 7s 7s/step\n",
      "Day 33\n",
      "1/1 [==============================] - 7s 7s/step\n",
      "Day 34\n",
      "1/1 [==============================] - 7s 7s/step\n",
      "Day 35\n",
      "1/1 [==============================] - 7s 7s/step\n",
      "Day 36\n",
      "1/1 [==============================] - 7s 7s/step\n",
      "Day 37\n",
      "1/1 [==============================] - 7s 7s/step\n",
      "Day 38\n",
      "1/1 [==============================] - 9s 9s/step\n",
      "Day 39\n",
      "1/1 [==============================] - 9s 9s/step\n",
      "Day 40\n",
      "1/1 [==============================] - 8s 8s/step\n",
      "Day 41\n",
      "1/1 [==============================] - 7s 7s/step\n",
      "Day 42\n",
      "1/1 [==============================] - 7s 7s/step\n",
      "Day 43\n",
      "1/1 [==============================] - 7s 7s/step\n",
      "Day 44\n",
      "1/1 [==============================] - 7s 7s/step\n",
      "Day 45\n",
      "1/1 [==============================] - 9s 9s/step\n",
      "Day 46\n",
      "1/1 [==============================] - 8s 8s/step\n",
      "Day 47\n",
      "1/1 [==============================] - 7s 7s/step\n",
      "Day 48\n",
      "1/1 [==============================] - 7s 7s/step\n",
      "Day 49\n"
     ]
    }
   ],
   "source": [
    "for i in days:\n",
    "    x1, x2, y_true, _ = augment_ts_length_to_days_since_trigger(X_ts_balanced, X_static_balanced, Y_balanced, astrophysical_classes_balanced, i)\n",
    "    logits = model.predict([x1, x2], batch_size=12100)\n",
    "    _, pseudo_probabilities = get_conditional_probabilites(logits, tree)\n",
    "\n",
    "    print(\"Day\", i)\n",
    "    L1[:, i] = pseudo_probabilities[:, true_idx[1]]\n",
    "    L2[:, i] = pseudo_probabilities[:, true_idx[2]]\n",
    "    L3[:, i] = pseudo_probabilities[:, true_idx[3]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 519,
   "id": "74f446c6-95bc-4956-ab2a-5f917791a952",
   "metadata": {},
   "outputs": [],
   "source": [
    "threshold = 0.90\n",
    "flag_value = -9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 520,
   "id": "b45752c8-b4c6-4298-8c9e-1cf4e156985e",
   "metadata": {},
   "outputs": [],
   "source": [
    "L1_days = []\n",
    "for i in range(N):\n",
    "    \n",
    "    above_threshold = np.where(L1[i, :]>threshold)[0]\n",
    "    if len(above_threshold) > 0:\n",
    "        days_to_threshold = above_threshold[0]\n",
    "        L1_days.append(days_to_threshold)\n",
    "    # else:\n",
    "    #     L1_days.append(flag_value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 521,
   "id": "28a70784-4774-448c-90bb-58ee3c508b58",
   "metadata": {},
   "outputs": [],
   "source": [
    "L2_days = []\n",
    "for i in range(N):\n",
    "    \n",
    "    above_threshold = np.where(L2[i, :]>threshold)[0]\n",
    "    if len(above_threshold) > 0:\n",
    "        days_to_threshold = above_threshold[0]\n",
    "        L2_days.append(days_to_threshold)\n",
    "    # else:\n",
    "    #     L2_days.append(flag_value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 522,
   "id": "98bd4534-7d69-4897-8827-1597cba16aed",
   "metadata": {},
   "outputs": [],
   "source": [
    "L3_days = []\n",
    "for i in range(N):\n",
    "    \n",
    "    above_threshold = np.where(L3[i, :]>threshold)[0]\n",
    "    if len(above_threshold) > 0:\n",
    "        days_to_threshold = above_threshold[0]\n",
    "        L3_days.append(days_to_threshold)\n",
    "    # else:\n",
    "    #     L3_days.append(flag_value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 523,
   "id": "07ce070a-5c56-4bae-889c-b7cbe0c8db51",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(11703, 8966, 7705)"
      ]
     },
     "execution_count": 523,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(L1_days), len(L2_days), len(L3_days)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 524,
   "id": "7d0c17d9-d1e5-465d-a704-5f2e82414cac",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAACC3UlEQVR4nO3dd3xT5f4H8E/SkY403XuCzDKlDFlShmxkKiBbUYECKtcB/pQhIF4HolDEqzJUKFyGDNkCBVSUoch0IKNldNDSPZM8vz+w5zZ0N2lPmnzer1debc45ec73nJM03z7rKIQQAkRERERU5ynlDoCIiIiITIOJHREREZGFYGJHREREZCGY2BERERFZCCZ2RERERBaCiR0RERGRhWBiR0RERGQhmNgRERERWQhbuQOoi/R6PW7fvg0XFxcoFAq5wyEiIiILJoRAZmYmAgICoFSWXyfHxK4abt++jeDgYLnDICIiIisSHx+PoKCgcrdhYlcNLi4uAO6fYI1GI18gy9sCWQmA2g+YcVq+OAiDvhmE5NxkeDt6Y9fQXXKHY7X+7tcf2qQk2Pr44KG9e+QOh4jIJDIyMhAcHCzlH+VhYlcNRc2vGo1G3sTOwQYoVNz/KWccBFsnW9jABrZOtvK+J6yc2tYWWhsb2NryOhCR5alM9y8OniAiIiKyEEzsqiA6Ohrh4eFo166d3KEQERERlaAQQgi5g6hrMjIy4OrqivT0dDb3kAGdTofCwkK5wyAiojrIzs4ONjY2JZZXJe9gHzsiE8nKysLNmzfB/5WIiKg6FAoFgoKCoFarq10GEzsiE9DpdLh58yacnJzg7e3N+Q2JiKhKhBBITk7GzZs30bBhw1Jr7iqDiR2RCRQWFkIIAW9vbzg6OsodDhER1UHe3t64fv06CgsLmdhZpdh3gLwMwEEDRM6WOxqrFnM5Bk0VTZGal4pAx0C5w7FahUlJgE4H2NjAzsdH7nCIiKrEFK09HBVbl51ZB/wUff8nyWr/jf3ILsxGRn6G3KHUiokTJ2LRokUAgPXr12PQoEG1un+1Wo3bt2+XWK5LvQdtSgp0qfdqNR5zcvz4cbRq1UruMAAAYWFh+P77701e7oPvuePHj+Ohhx6CWq3GqVOn0KxZM5w4caJG91mX1dR1IfPAxI7IwoWFhcHZ2RnZ2dnSspycHLi4uCAsLMzo8seMGYNdu0x7t41mzZpBrVZDrVZDqVTC0dFRen78+HFkZWUhICDApPt80Nq1a9GrV69qvbYoVrVaDYVCAWdnZ+l5XFyciSM11LVrV/z2229Gl1OZ49fpdFiyZAkaNWoEZ2dn1KtXDzNmzMDdu3eN3n95HnzPzZs3D6+99hqysrLQrl07XLx4ER07djRqHwqFAjdv3ixzn3LSarUYMWIEgoODoVAocP36dblDIjPCxI7ICgQGBmL79u3S8x07dsDf31++gCpw8eJFZGVlISsrC40aNcLevXul5127dpU7vAoVxZqVlQWVSmVwPCEhIQDud5TW6/UyR2qc559/HmvXrsWaNWuQlpaGX375BYGBgTh58mStxhEXF4fw8PBa3afcunbtiv/+979QqVRyh0JmholdFXCCYqqrRo8ejfXr10vPv/76a4wZM8Zgm7i4OAwYMACenp5o2rQp9u3bJ637+++/0blzZ7i4uGD48OHIzc2V1j1YszNjxgwEBATAzc0NvXv3NqihUigU+OSTT1CvXj14eXlhyZIl1Tqe4rUpkZGRmDt3LiIiIuAV0QZT5s7FnaQk9OjRAxqNBsOHD0d+fr702i1btqBZs2bw8PDA448/jqSkpBLlX716FVOmTEFsbCzUajWaNWsG4P79ofv37w93d3eEh4djx44dVYp74sSJmD59Onr06AEnJyf8/fffWL16NRo1agQXFxe0bNkSsbGx0vaRkZGYN28e2rZtC41Gg5EjR0rH8ueff6JLly7QaDTw9fXFK6+8AgCIjY1FgwYNpDLKu65hYWH44IMP0LRpU7i5uWH69OnlHn9xv//+O1avXo2YmBh07twZdnZ2cHd3x+zZs9G/f/8S2//8889o164dNBoNQkNDsXz5cmndTz/9hIcffhgajQaBgYH48MMPy11e/D3XrFkzXL16Fb1794afn590XEVNjVqtFnPnzkVoaChcXV0RGRkp7XfYsGHw8fGBh4cHnnjiCaSmpgIAevfuDQBo3LixVEv84Pv8+PHjePjhh+Hm5oZu3brh8uXL0jpTvc/LYmtrixdeeKHStZJ79+5FgwYN4OHhgQULFhis+/bbb9GiRQu4uLigYcOG2Lx5MwDg5s2b0Gg0BjX9a9askc4NmTFBVZaeni4AiPT0dHkDeb+JEPM093+SrEZ+M1J8d/I7cenOJblDKSE0NFTExsaKoKAgkZSUJJKSkkRQUJD4/vvvRWhoqBBCCJ1OJ1q2bCk++ugjUVhYKH788Ufh5eUlEhIShBBCtGvXTsyZM0fk5+eLb775Rtja2oqFCxcKIYRYs2aN6Nmzp7S/mJgYkZaWJnJycsSkSZPE4MGDpXUAxIgRI0RmZqY4f/68UKlU4sqVK+XG37hxY3HkyBGDZQBEfHy8EEKIbt26ifDwcBEXFyeuHTsufDw9RZvmzcWFCxdERkaGaN68uVi9erUQQoiff/5ZBAYGinPnzomCggLxyiuviOHDh5e63wePSwghOnfuLF5++WWRl5cnjhw5ItRqtfjrr7/KjV+lUolr164JIYSYMGGC8PT0FKdPnxaFhYWioKBA7N69W8TFxQmtViv+85//CF9fX5GXlycdW5MmTcT169fFvXv3RHh4uHQsI0eOFG+//bbQ6/UiKytL/Pzzz0IIIY4cOSIeeughIUTF1zU0NFR06dJFJCcni/j4eOHt7S0OHz5c5vEXt3LlSun9U5bQ0FBx/PhxIYQQZ86cEWfOnBE6nU6cOnVKaDQa8csvvwghhOjQoYP4+uuvhRBCpKamVrj8wdiK7+fB54sWLRLt2rUTN27cEFqtVhw9elTa7quvvhJZWVkiLS1N9OnTR7zwwgvSuuLvsQf3effuXeHm5ia2bt0qCgoKxLvvvisaNGggCgsLpddW9n0+depU4erqWupj6tSp5Z5fIQzfX6VJSkoSarVa7Nq1S+Tn54tXXnlF2NjYSOcnNjZW/P7770Kn04ndu3cLZ2dncefOHSGEEN27dxcbNmyQynrsscfEmjVrKoyJqi83N1dcunRJ5ObmGiyvSt7BUbFENenHFcCJ6Iq3828FPLXRcNmGUcCdcvpKdYwCOk2vVBg2NjYYMWIENm3aBAAYPny4wVD6kydPIjc3FzNnzrxfdMeO6NatG/bu3Yvu3bvj3LlzOHbsGOzt7TFkyBB06NChzH2NGjVK+v21115D586dDdbPnj0barUazZs3R8uWLXH+/Hk89NBDlTqOsjzzzDMIDg5GXnYOurZtCx9vb6mWqX///lKfs9WrV2PatGlo0aIFAODNN9+Eh4cHtFotbG3L/3MYHx+P06dP4+DBg1CpVIiMjMTAgQOxefNmzJkzp9KxDh8+HBEREdLz4rVbzz77LObOnYu//voLzZs3BwBMnjwZoaGhAIABAwZIx2JnZ4cbN24gISEB/v7+aN++fYl9lXddJ06cCAB48cUX4eXlBeB+DeFvv/2G7t27V3gcKSkpUg1ZZbRp00b6vW3btujfvz9++OEHPPzww7Czs8OVK1eQmpoKDw8PuLu7S8dY2vKqWLt2LVauXCk1gT/66KPSurFjx0q/v/TSS/i///u/SpW5Z88etGzZEsOGDQMA/Otf/8KyZctw6tQpqRatsu/zlStXYuXKlVU+rsras2cPIiIiMHDgQADA/PnzsWzZMml9t27dpN/79++PFi1a4PTp0xg4cCDGjh2LmJgYjB49GomJifjxxx+xZcuWGouVTINNsUQ1KT8TyLxd8SOnlM7mOXfLf01+ZpVCGTNmDDZs2ID169eX2gx77do1uLm5SY99+/bhzp07uHPnDry9veHg4CBtHxwcXOZ+Fi9ejAYNGkCj0aB9+/ZISUkxWO/r6yv97uTkhKysrCodR2l8ik1t4qhSwcfD83/PHR2lfcTFxWHx4sXSMQYHB8PW1hYJCQkV7uP27dsl5ikMDQ0tdXRueYKCggyeb9++HW3atJFiSkpKMjhnZZ2vd999FwUFBWjdujUefvjhUjv2l3ddKyq/Ip6enpU6b0UuXryIxx57DN7e3nB1dcW2bduk4/z8889x8eJFNGjQAF26dJFGtJa1vCpu3ryJevXqlViu1Wrx4osvIjQ0FBqNBiNGjCjxXi3L7du3pUQRAJRKJYKDgw3eCzXxPq+OO3fuGHxenZyc4On5v8/H999/j86dO8PDwwNubm44ffq0dB5GjBiBo0eP4t69e9i8eTP69u3L22jWAayxM0efdgOySvb7KSEr4X8/P2ha8fZqH+D5o8bFRlWjcgFcKjF608mr9GXlvVblUqVQ2rZtK/UhateuHX766SdpXWBgIJo2bYpz586VeN2NGzdw9+5d5OXlScldfHx8qf2ujh49ipUrV+LIkSNo2LAh/vzzTzRp0qRKcdakwMBALFy4ELNmzapw2wfnkwoICEBycrLBeYiLi5Nq/yqreLn5+fkYPXo0tm3bht69e8PGxgb+/v6Vui2dv78/Vq9eDSEEdu7ciSeffBL37hlO81Leda1KnKXp3r07oqKicPbsWbRu3brC8qZPn46uXbti586dcHR0xOjRo6XjbNy4Mf773/9Cq9Vi1apVGD16NK5fv17m8qoIDg7G9evXDfodAvenL4mNjcWPP/6IwMBA7N+/H88//3ylygwICMCePXuk50IIxMfHV2uk9pQpU/D111+Xum7s2LFYtWpVlcsszt/f36BfZW5urkECO27cOMyZMwcTJ06Evb09OnbsKF0XjUaDPn36YOvWrYiJicHLL79sVCxUO1hjZ46ykipXyyP+GVEn9JXbvjLJIplWp+nAvy5X/HiwGRa4v6y811SyGba4bdu2Ydu2bSWWd+jQAXq9Hp988gkKCgpQUFCA48ePIy4uDqGhoWjevDkWLVqEwsJC7Ny5s8xRj5mZmbCzs4OXlxeys7Olue7MxaRJk7BixQqpOTM1NbXMARA+Pj64efMmtFotgPsJQps2bTBv3jwUFBTg2LFj2LVrF0aMGFHtePLz81FQUCDVOH700UdITk6u1Gu3bNmC27dvQ6FQwM3NDQqFokQyVt51rciDx/+gJk2a4Omnn8bo0aNx4sQJaLVapKen47333jNIeopkZmbCzc0NDg4OOH78OHbv3i2tW79+PVJSUmBrawsXFxepm0BZy6ti4sSJeOONNxAfHw+dTodjx45J8Tg4OMDd3R13797F+++/X+L4y0oi+/Xrh99++w07duyAVqvFhx9+CEdHR7Rt27bK8a1atcpgFHXxR3lJXX5+PvLy8kr8/qD+/fvjzJkz2LNnDwoKCrBgwQKD0diZmZnw9PSEnZ0dtm7dijNnzhi8fuzYsfjwww9x+fLlUgfFkPlhYmfOFArA0aPsh409oLS9/7O87Xjf0hrX3LM5VDYqONqZ9+3EwsPDS50WwtbWFrt378b+/fsRGBiIgIAALF68WPoC2LBhAw4dOgQPDw+sXbsWQ4cOLbX8vn37onPnzggNDUWLFi3QqVOnGj2eBymdnQA7Oyjs7Upd36lTJ7z//vsYP348NBoN2rRpgx9++KHUbXv06IGwsDB4e3ujZcuWAICNGzfit99+g4+PD55//nmsW7cODRs2rHa8Go0G7733Hvr06QM/Pz+kpKSUqFkqy8mTJxEREQG1Wo2pU6ciJiamxNQXFV3X8pR2/A/69NNPMX78eIwfPx6urq5o3bo14uPjS+3v9+9//xvR0dHQaDRYtmwZHn/8cWndnj170LhxY7i4uODjjz/Gl19+We7yqnjllVfQs2dPdOrUCZ6enpg/fz4AYPz48XB3d4evry+6du2Kvn37Grxu7ty5GD58ONzc3EpM5uvl5YXt27dj3rx58PT0xDfffIPt27fDzq70911NaNy4MRwdHZGfn48mTZqUeStDb29vxMTEYMaMGfD19YWjo6NBd4Dly5dj5syZcHd3x/79+w363AH3k9jExEQMHTqUU6vUEQpRmTp/MpCRkQFXV1ekp6fXTH+DD5rer2Fz9ACGmKBT7fZpQG7q/Wa9f12ueHuqsry8PFy7dg316tUz6ItGRFTXNW/eHB999BF69uwpdygWr6zvkqrkHayxIyIiolIdPHgQOTk5lRopTeaBgyeIiIiohJEjR+K7777DunXroFSyHqiuYGJHREREJRTNe0l1CxO7KoiOjkZ0dDR0Op3codx3aCGQlw44uAI935Q7Gqv2f9//Hx53fxz2mfZ4yMG4yXap+vKvXYPQaqGwtYWqlLnLiIgsHetWqyAqKgqXLl3CqVOn5A7lvsw7QMbN+z9JVreybkGr16JAVyB3KFZN5BdA5OdD5PM6EJF1YmJHREREZCGY2BERERFZCCZ2RERERBaCiR0RVdnEiROl24WtX78egwYNqtX9q9Vqgxuuk/l4++23MX161W93Z2rXr1+HrW3NjA+cMmUK3n33Xen58uXL4ePjAy8vL8TFxcHNza3G92mpiv9toephYkdk4cLCwuDs7Izs7GxpWU5ODlxcXBAWFmZ0+WPGjMGuXbuMLqe4Zs2aQa1WQ61WQ6lUwtHRUXp+/PhxZGVlVeuG61Wxdu1a9OrVq9qvDwsLK3Erqtpw/Phx6Vw5OTlBqVRKz5s1a1bj+3/99dexYsUKo8upzBd8SkoKnn/+efj7+0Oj0aB58+Z4//33a3zmglWrVuHVV18FABQUFODVV1/F999/j7t37yIkJARpaWlGlR8bG1vi1nLF9ym3S5cuoXfv3nB3dy/1b8j8+fPRrFkzKJVKrF27ttbjs3ZM7IisQGBgILZv3y4937FjB/z9/eULqAIXL16UboTeqFEj7N27V3retWtXucMza127dpXO1eHDhxESEiI9v3jxorSdVquVMUrj5ebmIjIyErdu3cL333+P9PR07Nq1C+fOnUNmZmatxZGUlITCwkI0atSo1vYpNzs7O4waNQofffRRqesbNGiApUuXokuXLrUcGQFM7IiswujRo7F+/Xrp+ddff40xY8YYbBMXF4cBAwbA09MTTZs2xb59+6R1f//9Nzp37gwXFxcMHz4cubm50roHa7ZmzJiBgIAAuLm5oXfv3oiLi5PWKRQKfPLJJ6hXrx68vLywZMmSah2PQqHAzZs3AQCRkZGYO3cuIiIi4BXRBlPmzsWdpCT06NEDGo0Gw4cPR35+vvTaLVu2oFmzZvDw8MDjjz+OpKSkEuVfvXoVU6ZMQWxsrEFNV3x8PPr37w93d3eEh4djx44dVY69vDIiIyMxb948tG3bFhqNBiNHjjSI/a233oKPjw/CwsKwfPnyKjc1hoWF4d1330XTpk2lGqHqXq9vv/0WjRs3lmp+N27cCOB+bc3kyZOl7Y4ePYqIiAi4ubkhMjISf//9N4D/NZV+9tln8Pf3h5+fH9atWwcAWLduHdavX4+FCxdCrVZjypQpJY7lyy+/xN27d7F582Y89NBDUCgUqFevHr788stSm0JXr16NRo0awcXFBS1btkRsbKzButDQULi4uKBx48bSurKWF9UmxsXFoXHjxtDpdFCr1XjiiSdKNAEnJyfjqaeekppqZ8+eDeD+Z+rRRx+Fm5sbAgIC8PrrrwMAdDod+vXrh6tXr0o1rcX3WWTlypWoX78+vL29MXbsWKSnpwP4X23fW2+9BQ8PD4SFhWH//v1lvSWqpWHDhnj66afLTGbHjh2LPn36wMnJqcKyyvvbUtY5AmBwPQDg2rVrcHNzQ15eXvUPzEIwsSOyAj169MD58+eRnJyM5ORknDt3ziAZ0+v1GDRoEPr06YPExESsXr0a48aNQ2JiIoD7iWG3bt2QkpKCcePGYdu2bWXuq3Pnzrh8+TLu3LmDoKAgzJw502D94cOHcf78ecTGxmLBggXSF70xtm7diu3bt+PCvv3Yd+wYRkRNw/Lly3Hr1i38+eef2LBhAwDg5MmTePHFF7Fx40YkJiaiSZMmmDZtWony6tevj1WrViEyMtKgpmv06NFo1qwZEhISsHLlSowdOxZXrlypUqwVlfHf//4XW7duRVxcHC5cuCDFvmvXLnz++ec4ceIEfvvtN3z77bfVOlfbtm1DbGwsLl++DKD612vy5MlYvXo1MjMz8dNPP6Fly5Yl9hUfH48RI0Zg2bJlSElJwfDhwzFq1ChpvU6nw4ULF3Djxg18+eWXiIqKQmZmJiZMmIAxY8bgzTffRFZWFlatWlWi7MOHD6NPnz5wdHSs1HH7+fnh0KFDSEtLw4wZMzBq1Cjk5+cjOzsbL774Ir777jtkZmbiwIEDCA0NLXN5cSEhIbh48SJsbGyQlZWFzZs3l9jvmDFj4OTkhL///hvx8fEYPHiwtG7hwoW4e/cujh49iq+//hrbt2+HjY0N9u7di/r160s1rQ86ePAgFi5ciG+//RbXr19Hbm4uXnjhBWn99evXoVKpkJSUhNdffx3PPfdcmeelZcuWcHNzK/VR9N6rSRX9bSntHAH3z2tMTIy0XUxMDIYMGQIHB4caj9nc8c4TdVnzYYA2H7BVyR2J1RvVeBQ0hRq4O7obLF93cR2+vPRlha8P9wjH8p7LDZbNODQDl1Ivlfma8eHjMaHZhErFZ2NjgxEjRki3CBo+fDhsbGyk9SdPnkRubq70pd6xY0d069YNe/fuRffu3XHu3DkcO3YM9vb2GDJkCDp06FDmvop/cb/22mvo3LmzwfrZs2dDrVajefPmaNmyJc6fP4+HHjLubh3PPPMMgoODoXV2RrcuXeDj4yPVsvXv3x+//fYbgPs1MNOmTUOLFi0AAG+++SY8PDyg1WorrP2Kj4/H6dOncfDgQahUKkRGRmLgwIHYvHkz5syZU6k4K1PG5MmTpQRiwIABUuxbt27F5MmTpXP16quv4tChQ1U8U8ALL7wAX19f6Xl1r5ednR0uX76MVq1awc/PD35+fiX2tX79egwdOlRqPp8xYwbeeustXL9+Xdpm7ty5sLe3R+/evaUEqHXr1hUeR0pKCtq2bVvp4+7fv7/0+7PPPou5c+fir7/+Qr169aBQKHDx4kWEhIRI5z47O7vU5VVx69YtxMbG4t69e3B2dgZw/7MFAA899JB0LRs2bIgxY8bg+++/x5AhQyosd+PGjXjuuecQHh4O4P6AldatW2PNmjUAAGdnZ7zyyitQKpUYO3Ysnn/+eaSlpZVak3nu3LkqH5ep3Lhxo9y/LeWdo7Fjx6J9+/ZYsWIF7OzsEBMTgw8++ECuQzErrLGryxr0ApoMuP+TZNW3Xl+o7dRwU7kZLM8uzEZSTlKFj9T81BJlpuanlvua7MLsEq8pz5gxY7BhwwasX7++1GbYoqaMose+fftw584d3LlzB97e3gb/CQcHB5e5n8WLF6NBgwbQaDRo3749UlJSDNYXTyqcnJxKrZGoKh8fHwCArYcHnNzc4BcSIq1zdHSU9hEXF4fFixdLxxgcHAxbW1skJCRUuI/bt2/D29vboIYoNDS0SqNzK1NGWecnISEBQUFB0rriv1fFg6+r7vXasmULvvnmGwQFBaFv375SDWBxcXFx+OqrrwzeV9nZ2bh16xaA+/9weHp6llp+RTw9PSt13Yps374dbdq0keJISkpCSkoKnJ2dERMTg48//hi+vr544okncPv27TKXV8XNmzfh4+MjJXXF3bp1C0OHDoWfnx9cXV2lWs3KuH37NkKKvcdDQ0ORl5eH1NT7f0e8vb2hVN7/ei9qDjXF58zUKvrbUt45ql+/Pho1aoQDBw7gwoULSE5ORs+ePWv9GMwREzuiGuRs5wwfJ58KHx4qjxKv9VB5lPsaZ7uSXxbladu2LVJTU3Hv3j20a9fOYF1gYCCaNm2KtLQ06ZGVlYU5c+bA398fd+/eNei7Eh8fX+o+jh49ipUrV2LPnj1IT0/HyZMnqxRjTQsMDMTChQsNjjM3N7fUJEmhUBg8DwgIQHJyssF5iIuLq9LoXGPK8PPzkxIiAFIfw6oqflzGXK8OHTpg9+7dSExMRKtWrTB16tQS2wQGBuLZZ581ON85OTklagUrirM0PXr0wIEDBwz6ZJUlPz8fo0ePxuLFi5GSkoK0tDT4+PhACAHgfm3e4cOHcfPmTahUKqkvV1nLKys4OBjJycnIyckpse6NN96Au7s7/vzzT6Snp+PFF1+U4qno2AMCAgz6QsbFxcHBwQEeHiX/jlSk+Aj0Bx/F++XWhIr+tpR3joD7ffliYmIQExODkSNHGrRCWDMmdkQ1aEKzCTj0xKEKHw82wwLA8p7Ly31NZZthi9u2bVup/eM6dOgAvV6PTz75BAUFBSgoKMDx48cRFxeH0NBQNG/eHIsWLUJhYSF27txZZgKQmZkJOzs7eHl5ITs72+zmo5o0aRJWrFghNW+mpqaWOQDCx8cHN2/elEaPBgcHo02bNpg3bx4KCgpw7Ngx7Nq1CyNGjChzfwUFBcjLy5MegYGBVS6jyLBhw/DFF1/g6tWryMjIwPvvv1+NM2CouteroKAAGzZsQEZGBuzs7KBWq0v9Un3qqaewefNmHD9+HHq9HpmZmdiyZUul9uHj42PQZPug8ePHw9PTE6NGjcLVq1cB3G/amzhxYonpRvLz81FQUCDV7H700UdITk4GACQmJuLbb79Fbm4uVCoVnJycYGNjU+byqggICEC3bt3w0ksvISsrC7m5ufjpp58A3D/3Li4uUKvVuHDhAr7++muDY09OTjaYoqi4kSNH4rPPPsPly5eRnZ2N//u//8OTTz5ZYUJYmuIj0B98PFizX0QIgby8PBQUFBj8XqSwsBB5eXnQ6/UGvz+oor8t5Z2jovOwZ88erF+/Hk899VSVj91SMbGry3LvATkp93+SrFJzU6HT61CoK5Q7lHKFh4dL/XKKs7W1xe7du7F//34EBgYiICAAixcvlv4Yb9iwAYcOHYKHhwfWrl2LoUOHllp+37590blzZ4SGhqJFixbo1KlTjR7Pg/SFhRB6PUQZ85h16tQJ77//PsaPHw+NRoM2bdrghx9+KHXbHj16ICwsDN7e3tLAgI0bN+K3336Dj48Pnn/+eaxbtw4NGzYsM56ePXvC0dFRemzYsKHKZRR5/PHHMXHiRHTo0AEtWrTAY489BpXKuP61xlyvdevWITQ0FO7u7jh48GCpc9fVq1cPGzduxCuvvAIPDw80adKk0iOJn376afz8889wc3MrdYCLo6MjYmNj4evri06dOsHFxQX9+/dHixYt4OLiYrCtRqPBe++9hz59+sDPzw8pKSnSqGC9Xo93330Xvr6+8PHxwa1bt7Bo0aIyl1fV+vXrkZaWhrCwMISEhGDnzp0A7vctPHLkCDQaDWbOnInhw4dLr2natCkGDx6M4ODgUvvF9e7dG3PmzEH//v0RGhoKOzs7LFu2rMqxVdeNGzfg6OiIbt26IS4uDo6Ojujdu7e0/tlnn4WjoyMOHjyI5557Do6Ojjh27FipZZX3t6W8cwTcb47v0qULbG1ty+33a20Uoni9JpUrOjoa0dHR0Ol0UtWwRqMx/Y4+aApk3gYcPYAhK8vebvs0IDe18tu5BAD/KtkPhow3avsoPBv4LAKCA9DUr6nc4VitvN//gNAWQmFrB4cmjeUOp0bt378fM2bMwJ9//il3KESymT59Otzd3bFw4UK5QzGJvLw8XLt2DfXq1TPoe5iRkQFXV9dK5R2ssauCqKgoXLp0CadOnZI7FCKyQtu3b0dhYSESEhLw1ltvGUydQWRtEhIS8N///heTJk2SOxSzwsSOiKiOWLZsGTw8PNCiRQs0bNgQc+fOlTskIlmsWLECDRo0wNSpU1G/fn25wzErnMeOiKiOKD7TPpE1mz59OqZPny53GGaJNXZEREREFoKJHREREZGFYGJHREREZCGY2BERERFZCCZ2RERERBaCiR0RVdnEiROlWfjXr1+PQYMG1er+1Wp1lW/Ibi7i4uJKvZtAkeLntjyxsbHS3ROq4/jx42jVqlW1X28q/fr1w6ZNm+QOwywoFIpq3wO4PG+//bbBCNItW7YgMDAQarUaSUlJNfJ5enCflqqyn9faxMSuLuvxBtD/vfs/SVaLOi+Ct5M3AtWBcodSQlhYGJydnQ3uO5mTkwMXFxeEhYUZXf6YMWOwa9cuo8sprviNyZVKJRwdHaXnx48fR1ZWFgICAkq8zr5eGFQNGsC+XpjRMaxduxa9evWq1muL30hdoVDA2dlZeg6gxL1MTe369etQKBRQq9VwdnZGUFAQRo0ahfPnz0vbdO3aVbpnbm0p7Zzu3bsXI0eONOl+ih9/8Ufx4ze1+fPnY/LkyeVuk5ubi1dffRWhoaFQq9Vo1KgR3njjjTLvCWsqr7/+usEt31599VWsWbMGWVlZ8PHxKfPzVFnXr1+Hra3h7GkP7lNOu3fvRqdOneDq6oqAgADMmjULhYX/u/3j/Pnz0axZMyiVSqxdu1a+QE2EiV1dpgkAXIPv/yRZBbkEwU5pB3tbe7lDKVVgYCC2b98uPd+xYwf8/f3lC6gCxW9M3qhRI+zdu1d63rVr1zJfp1SpoHRwgNLIe6gaq/iN1FUqlcHxhISE1EoMNjY20j5//PFHBAUFoUOHDvj111+rXaZWqzVhhDWr+PEXPVq0aCFbPEIIDB48GMeOHcOuXbuQkZGBY8eOITc3F3///XetxhIXF1fqPaMtVUZGBubPn4+EhAT89ttvOHXqFN577z1pfYMGDbB06VJ06dJFxihNh4kdkRUYPXo01q9fLz3/+uuvMWbMGINt4uLiMGDAAHh6eqJp06bYt2+ftO7vv/9G586d4eLiguHDhyM3N1da92AtzIwZMxAQEAA3Nzf07t0bcXFx0jqFQoFPPvkE9erVg5eXF5YsWVKt4yneZBUZGYm5c+ciIiICarUazzzzDO7cuYMePXpAo9Fg+PDhyM/Pl167ZcsWNGvWDB4eHnj88ceRlJRUovyrV69iypQpiI2NhVqtRrNmzQAA8fHx6N+/P9zd3REeHl7pG9oX92DtRnnn9u+//8ajjz4KNzc3BAQE4PXXX6/y/hQKBUJCQvD++++jb9++eOuttwAYNuXq9XrMnDkTXl5ecHNzQ7t27XD37l3p9StWrEC9evXQvXv3UpuAi1+PsLAwvPvuu2jUqBE8PT3x8ssvQ6/Xl3lOIyMj8fXXX0txzJs3D8HBwfD398fMmTOla7d27Vr06NEDU6dOhUajQXh4OH755Zcqn49FixaVuAVV9+7dpRjOnz+PRx99FO7u7oiIiMDp06cNjrO0929sbCzefvttrFu3Dmq1Gv369Sux3wMHDuD777/Htm3b0LJlSyiVSvj5+eGDDz5Ay5YtS2z/7bffokWLFnBxcUHDhg2xefNmg3WNGzeWat03btxY7vLitYlqtRo6nQ6NGzdG+/btpeMqun7Z2dmYNm0aAgIC4O7ujnHjxgEA7t27h759+8LLywve3t547rnnpGvTu3dv6HQ6qWY0Li6uRA3mtm3b0LRpU7i7u2PgwIG4desWgP99Hj777DP4+/vDz88P69atq9zFrKTRo0ejd+/ecHR0hLe3N8aNG4cTJ05I68eOHYs+ffrAycmpwrKq+3lt3LixwQTj165dg5ubG/Ly8kxzkMUwsSOyAj169MD58+eRnJyM5ORknDt3ziAZ0+v1GDRoEPr06YPExESsXr0a48aNQ2JiIoD7fxi7deuGlJQUjBs3Dtu2bStzX507d8bly5dx584dBAUFYebMmQbrDx8+jPPnzyM2NhYLFiwwSW3F1q1bsX37dly5cgW7d+/G448/juXLl+PWrVv4888/sWHDBgDAyZMn8eKLL2Ljxo1ITExEkyZNMG3atBLl1a9fH6tWrUJkZCSysrJw8eJF6Tw0a9YMCQkJWLlyJcaOHYsrV64YFXtF53bhwoW4e/cujh49iq+//tqg5rWqHn/8cfzwww8llh84cAA//vgjrl69ipSUFHz66acGNyA/ePAgfvvtNxw4cKBS+4mJicGxY8dw/vx57N27F2vWrCnznBb3xRdfYMuWLThx4gQuXLiAM2fOGCT/x48fx6OPPop79+5h2LBheOmll6p8DkaOHIkdO3ZITXEJCQk4efIkBg8ejKysLPTt2xcvvPAC7t69izfffBPDhg0z+PIt7f0bGRmJ119/HRMmTEBWVhb27t1bYr+HDx9G+/btK93k6eLigi1btiA9PR0fffQRJk2ahISEBADA5MmTsXr1amRmZuKnn36SEsOylheXlZUFAPjjjz9w8uTJEutffPFFxMXF4dy5c0hKSsLzzz8P4P7fiKioKNy6dQvnzp3D6dOn8cknnwC4//4pXkP6YK3077//jkmTJuGzzz5DQkIC6tevj7Fjx0rrdTodLly4gBs3buDLL79EVFQUMjMzSz0vAwcOhJubW6mPd955p1Ln9tixY9I/FlVV3c/rmDFjEBMTI20XExODIUOGGHzOTIW3FKvLrn8P6AoAG3sgzDKqkOuq2PhYeBd6I6Mgw+CDmrJmLVIr0WfDITwcwZ+sNFgWP3Ua8i5dKvM1HhMnwnPSxErFZ2NjgxEjRkid1IcPHw4bGxtp/cmTJ5GbmyslYR07dkS3bt2wd+9edO/eHefOncOxY8dgb2+PIUOGoEOHDmXua9SoUdLvr732Gjp37mywfvbs2VCr1WjevDlatmyJ8+fP46GHHqrUcZTlmWeeQXBwMLRpaXi0Y0f4+PhIf7j79+8v9SVbvXo1pk2bJjXJvfnmm/Dw8IBWqy3RR+hB8fHxOH36NA4ePAiVSoXIyEgMHDgQmzdvxpw5c6oV940bN8o9tw899JB0bho2bIgxY8bg+++/x5AhQ6q1Pz8/P9y7d6/Ecjs7O2RmZuL3339Hu3bt0KZNG4P1s2fPhkajqfR+XnjhBfj5+QEAXnrpJWzatAnPPPNMha/buHEjXn75ZQQFBQEA5s6di5kzZ2L+/PkAgCZNmmD06NEAgKeeeqrcPlw6na7EIJWUlBQ0bNgQYWFhOHDgAAYMGIAtW7agT58+cHFxwcaNG9GsWTMMHz4cADBkyBAsWrQIP/30EyIjI6VzUZ33b0pKinROKqNbt27S7/3790eLFi1w+vRpDBw4EHZ2drh8+TJatWoFPz8/qdyylleWXq/HV199hfPnz8PLywsApOZJT09PaZCUv78/nn/+eRw8eBAvvvhiheVu3rwZQ4YMkcp6++234e7ujjt37kjbzJ07F/b29ujduzecnJzw999/o3Xr1iXK+vbbb6t0TA/aunUrDh06VK3+pcZ8XseOHYv27dtjxYoVsLOzQ0xMDD744AOjjqUsrLGry85uAE7+5/5PktXai2uRlp+GlJwUg+X6rCxoExMrfOhSU0uUqUtNLfc1+n/+866sMWPGYMOGDVi/fn2pzbBFTQNFj3379uHOnTu4c+cOvL29DRLW4ODgMvezePFiNGjQABqNBu3bt0dKiuE58fX1lX53cnKSahCM4ePjAwDQJiRCpdfDy8FRWufo6CjtIy4uDosXL5aOMTg4GLa2tlJNSHlu374Nb29vODr+r+zQ0FCjRhNWdG5v3bqFoUOHws/PD66urli2bFmJ81nV/bm7u5dY3rNnT0yZMgXPPfcc/P398fLLLxt0Li9KtCqr+DEEBwcbfIGX5/bt2wa1PQ+e36q8d2xsbJCWlmbwKPpnZtSoUdI/OZs2bZIGb8TFxeHo0aMGn4PLly9XO4biPD09K/U+K/L999+jc+fO8PDwgJubG06fPi1d+y1btuCbb75BUFAQ+vbti8uXL5e7vLKSk5ORn5+PevXqlViXmZmJ8ePHIygoCBqNBrNmzar0e/HB66pWq+Hp6SmdVxsbG3h6ekrrTfV34UFHjhzB1KlTsWvXLulvRlUY83mtX78+GjVqhAMHDuDChQtITk5Gz549jT+oUjCxI6pBSrUatr6+FT5sPDxKvNbGw6Pc1yj/GWFZWW3btkVqairu3buHdu3aGawLDAxE06ZNDb4Es7KyMGfOHPj7++Pu3bsGzVHx8fGl7uPo0aNYuXIl9uzZg/T09FKbeuQUGBiIhQsXGhxnbm5uqYmLQqEweB4QEIDk5GSD8xAXF2fUaMKKzu0bb7wBd3d3/Pnnn0hPT8eLL74IIUS19/ftt9+WqEEt8tJLL+Hs2bM4deoU9u/fb9Ans/i5cHZ2NuhXVNRcX1zxY4iPj5cG6jx4Th8UEBBg0CfT2PNblieffBI7d+7E1atXcfbsWQwcOBDA/fdHnz59DN4f2dnZeOqppyoss6Jj69GjB06ePFnp5G7cuHGYMGECEhISkJaWhrZt20rXvkOHDti9ezcSExPRqlUrTJ06tdzlleXt7Q2VSoXr16+XWLd06VIkJyfj7NmzyMjIwNKlS6V4qnpds7OzkZKSUq1r269fvxKjnYseb7/9dpmv+/nnn/Hkk0/iv//9L9q2bVvl/QLGf17Hjh2LmJgYxMTEYOTIkQatJqbExI6oBnlOmoiGR2MrfDzYDAsAwZ+sLPc1lW2GLW7btm2l9o/r0KED9Ho9PvnkExQUFKCgoADHjx9HXFwcQkND0bx5cyxatAiFhYXYuXNnmQlbZmYm7Ozs4OXlhezsbLOb32nSpElYsWKF1AyTmppa5gAIHx8f3Lx5UxoJGhwcjDZt2mDevHkoKCiQRjeOGDGi2vFUdG4zMzPh4uICtVqNCxcuSB38q0IIgZs3b2L27NnYs2cP3nij5PRIp0+fxqlTp6DVauHi4gI7O7syv3QaNWqEe/fu4ejRo8jPz8fChQtLbLN8+XIkJibizp07WLZsGZ588kkAJc/pg0aOHIkPPvgAt27dQmpqKhYuXGjQtG8qYWFhaNq0KZ599ln0798fzs7OAO733/r111+xfft2aLVa5ObmYt++fUhPT6+wTB8fH9y4caPMxLt3797o3Lkzhg8fjvPnz0Ov1yMpKQmvvfYazp07V2L7zMxMeHp6ws7ODlu3bsWZM2cAAAUFBdiwYQMyMjJgZ2cHtVoNGxubMpdXhVKpxPjx46XauMLCQqlPZmZmJpycnODq6oobN25g5cr//c3y8vKCXq8vcw6+ESNGYPv27fjxxx9RUFCAN954A506darWyPziI+QffJQ1uOj8+fMYNGgQvvjiC6lJvbjCwkLk5eVBr9cb/P4gYz+vI0eOxJ49e7B+/fpK/bNQXUzsiKxIeHh4qdMc2NraYvfu3di/fz8CAwMREBCAxYsXS3/cNmzYgEOHDsHDwwNr167F0KFDSy2/b9++6Ny5M0JDQ9GiRQt06tSpRo+nqjp16oT3338f48ePh0ajQZs2bUodTADcr2EJCwuDt7e31Al948aN+O233+Dj44Pnn38e69atQ8OGDY2KqbxzO3fuXBw5cgQajQYzZ86U+n5VRvFRiu3bt8e1a9dw4sQJRERElNg2PT0dTz/9NNzc3NC4cWN07ty5zC8eV1dXfPTRR3jyySdRr169ErW/wP0asa5du6J58+Z47LHHpFGopZ3T4p555hkMHToU7du3R3h4OFq1alXt/ovFj7/oUXyk98iRI3H48GEp6Sw6tt27d2P58uXw8fFBWFgY/vOf/1RqfyNGjEBWVpY06vNBCoUCO3bsQKdOnTBgwABoNBp07twZ9vb2pU40vXz5csycORPu7u7Yv3+/QZ+7devWITQ0FO7u7jh48KDU17Cs5VWxdOlSBAQEoFmzZvD19ZWO/4UXXpCa8ocPH27wPnV2dsbs2bPRunVruLm5GdTOAUDTpk3x+eefY9KkSfD19cUff/xRrX9Sqmvp0qVISUnBU089Jb0Xio9cfvbZZ+Ho6IiDBw/iueeeg6OjI44dO1ZqWcZ8Xj09PdGlSxfY2tqW20/ZWAphTL2+lcrIyICrqyvS09Or1KG40j5oCmTeBhw9gCEla3Ik26cBuamV384lAPhX1fpcUOWM2j4KzwY+i4DgADT1ayp3OFYr7/c/ILSFUNjawaFJY7nDsUphYWH4+uuvLWZOMCJTmj59Otzd3Uut6QaAvLw8XLt2DfXq1TPoy1eVvIOjYqsgOjoa0dHR0Ol0codCREREdUhCQgL++9//4qeffqrR/bAptgqioqJw6dIlnDp1Su5QiIiIqI5YsWIFGjRogKlTp6J+/fo1ui/W2BERkcmUNqKSyNpNnz4d06dPr5V9scaOiIiIyEKwxq4uc3Qz/EmycVO5wUZhAxtlzcxLRJWj+OfuEYoK7iJBRGSOTDGelX/96rI+ZU/GSLVrea/luHLlCnxVvhVvTDVG1cC4W5MREcmpoKAAAIyavJiJHZEJ2NrawsnJCcnJybCzs4NSyV4ORERUeXq9HsnJyXBycqrw3tXlYWJHZAIKhQL+/v64du0abty4IXc4RERUBymVSoSEhFR4m7byMLEjMhF7e3s0bNhQqkonIiKqCnt7e6NbfJjY1WUnPwMKsgB7NdD+WbmjsWoLTixAen46XFWumNdxntzhWK07c+dBl54OG1dX+L+1QO5wiIhqHRO7uuz2r/+7pRjJ6tjNY0jKSYKPk4/coVi1rKNHoU1MhK0vB7EQkXViD28iIiIiC8HEjoiIiMhCMLEjIiIishBM7IiIiIgsBBM7IiIiIgvBxI6IiIjIQjCxIyIiIrIQTOyIiIiILAQnKK7LQjsBBdmAvbPckVi9/vX6I6MgAxp7jdyhWDXNgAHQZaTDRuMqdyhERLJgYleXPTxW7gjoH/9q+y+5QyAAvq++IncIRESyYmJnhka62uCuWwCgUAK/LDa+QG8nQDjAS9hgk/GlERERkZliYmeG7iqBJJt/Lk1BuvEF2igBKAGdML4sIiIiMltM7MyYUgho7IzvP5dRmA29QmGCiIiIiMicMbEzY156gTcferLM9e2PR8M+LxMFDi442TWqzO0WXlqNJBsmdjVp0DeDkJybDG9Hb+waukvucKzW3/36Q5uUBFsfHzy0d4/c4RAR1TpOd1KH2WgLYKsrgI22QO5QrF6ONgfZhdnI0ebIHYpV0+fkQJ+dDX0OrwMRWScmdkREREQWgokdERERkYWw2sRu6NChcHd3x4gRI+QOhYiIiMgkrDaxe+GFF/Dll1/KHQYRERGRyVhtYhcZGQkXFxe5wyAiIiIyGbNM7I4dO4ZBgwYhICAACoUC27dvL7FNdHQ0wsLC4ODggA4dOuDkyZO1HygRERGRGTHLxC47OxutWrVCdHR0qes3bdqEWbNmYd68efjll1/QqlUr9OnTB0lJSdI2rVu3RvPmzUs8bt++XVuHQURERFSrzHKC4n79+qFfv35lrl+6dCmeffZZTJo0CQCwatUq7N69G6tXr8bs2bMBAGfPnjVZPPn5+cjPz5eeZ2RkmKxsIiIiIlMxy8SuPAUFBThz5gzmzJkjLVMqlejVqxdOnDhRI/tcsmQJFixYUCNlG+PPZgOg1GuhV9a5y2hx5j4yF3m6PDjYOMgdilXzmz8PIj8fCpVK7lCIiGRR5zKCu3fvQqfTwdfX12C5r68vfv/990qX06tXL/z222/Izs5GUFAQNm/ejI4dO5a67Zw5czBr1izpeUZGBoKDg6t3ACaU4tNI7hDoH92Cu8kdAgFw6d5d7hCIiGRV5xI7U/nuu+8qva1KpYKKNQBERERk5sxy8ER5vLy8YGNjg8TERIPliYmJ8PPzkykqIiIiIvnVucTO3t4eEREROHTokLRMr9fj0KFDZTalWip1+m1o7sVDnc6RvnK7mHIRZ5PO4mLKRblDsWq5Fy4i59dfkXuB14GIrJNZNsVmZWXhypUr0vNr167h7Nmz8PDwQEhICGbNmoUJEyagbdu2aN++PZYtW4bs7GxplKy1aPHLJqjyM5GvcsGJ7i/JHY5Vm3l4JpJykuDj5INDTxyq+AVUI25GRUGbmAhbX180PBordzhERLXOLBO706dPo3uxTtBFAxcmTJiAtWvXYuTIkUhOTsbcuXORkJCA1q1bY9++fSUGVJhadHQ0oqOjodPpanQ/RERERNVhloldZGQkhBDlbjN9+nRMnz69liK6LyoqClFRUcjIyICrq2ut7puIiIioInWujx0RERERlY6JHREREZGFYGJHREREZCGY2BERERFZCCZ2VRAdHY3w8HC0a9dO7lCIiIiISmBiVwVRUVG4dOkSTp06JXcoRERERCUwsSMiIiKyEGY5jx1Vzsmu0+QOgf6xc8hOCCGgUCjkDsWq1d+9G4AAwOtARNaJiV0dprNVyR0C/cPZzlnuEAiAjZrXgYisG5tiiYiIiCwEE7sq4KhYIiIiMmdsiq0Cc7tXbNC1E7DV5kNrq8LNeh3lDseqrbu4DtmF2XC2c8aEZhPkDsdqpaxZC31WFpRqNTwnTZQ7HCKiWsfErg4Lvv4TVPmZyFe5MLGT2ZeXvkRSThJ8nHyY2Mkode1aaBMTYevry8SOiKwSm2KJiIiILAQTOyIiIiILwcSOiIiIyEIwsSMiIiKyEEzsiIiIiCwEE7sq4Dx2REREZM6Y2FVBVFQULl26hFOnTskdChEREVEJTOyIiIiILAQnKK7DMjX+yCvQoNCeNz6XW7hHOPyc/eCh8pA7FKvmEB4OnZ8fbDx4HYjIOjGxq8MuRIySOwT6x/Key+UOgQAEf7JS7hCIiGTFplgiIiIiC8HEjoiIiMhCMLEjIiIishDsY1eHNT+zEXYF2Si0d2Z/O5nNODQDqfmp8FB5sL+djOKnToMuNRU2Hh7sb0dEVomJXRVER0cjOjoaOp1O7lAAAC4Zd6DKz0S+ykXuUKzepdRLSMpJgo+Tj9yhWLW8S5egTUyEra+v3KEQEcmCTbFVwAmKiYiIyJwxsSMiIiKyEGyKNUOvrtPBOQdQCkBjt6HM7ZTdcgAHQJmeA6+Xyt5uYaEeeoUe2U4Anq6BgImIiMgsMLEzQy7ZgHtW0bOcsjfUC+mnzb2yt3P/56dSmCI6IiIiMldM7MyYXgHAWVX2BgqF9FOvLme77HwmdURERFaAiZ0ZS3cG8MyjZa7XJ28G9DnQO6uQXs52+OJgsRpAIiIislQcPEFERERkIZjYEREREVkINsXWYbedwmErCqFV2MkditUbHz4e2YXZcLZzljsUq+YxcSL0WVlQqtVyh0JEJAsmdnXYHedmcodA/5jQbILcIRAAz0kT5Q6BiEhWbIqtgujoaISHh6Ndu3Zyh0JERERUAhO7KuAtxYiIiMicsSm2DlPqC6GAgIACeiX72ckpuzAbQggoFAr2s5ORLisbgACggI2a14GIrA8Tuzrs4ZTtUOlzkK90whnvJ+QOx6o9vv1xJOUkwcfJB4eeOCR3OFbr6oAB0CYmwtbXFw2PxsodDhFRrWNTLBEREZGFYGJHREREZCGY2BERERFZCCZ2RERERBaCiR0RERGRhWBiR0RERGQhmNgRERERWQgmdkREREQWgokdERERkYXgnSfqsN/dekAhdBAKm0ptr1cAPTf3NNn+vRy9sGngJpOVV5d93ONjFOoKYWfDW7vJKSg6GqKwAAo7e7lDISKSBRO7KoiOjkZ0dDR0Op3coQAAsu08q/yapJykGoiEmnk2kzsEAuDYnNeBiKwbE7sqiIqKQlRUFDIyMuDq6ip3OJWm/Oem6EoBuKncjC4vPT8dAsLocoiIiMi0mNhZAQ+dgA4KeOkFFnZeaHR5b/7wJtLy04wPjIiIiEyKiV0d5p4fD6XQQa+wwT1VcK3v/27uXZP12avr/fWOxh9Fni4PDjYO6BbcTe5wrFbmkSMQ+flQqFRw6d5d7nCIiGodE7s6rH7GT1Dpc5CvdMIZ79pP7PRCzz57/3jrp7eQlJMEHycfHAo+JHc4Vith/gJoExNh6+vLxI6IrBITO6oyF3sXk5XF/npERESmw8SOquzVdq+arCz21yOqWdeGj4D27l2Tl2vr5YV6W7eYvFwiMg4TOzILpuyvB9T9Pntk3UyZjGkTE01STmnl/tUt0mTlMVEkMg0mdmQW2F+P6H+0d+/WSEJm4+ZmdBm6tDTp95pKGomo+pjYkaxM2V8PYJ89sjAKBWxMNGemUqOB7+zZRpeT+M470GdkmCCi+3Tp6YAQ0CYnswaQyASY2JGsTNlfD2CfPbIsNq6u8H/7bbnDMGCK5LC4O6+/fr8WUK9nDSCRCTCxoyrzmLYEynum+4/dlP6dnwE99Mh2SQaekDsaIqqIUqMxaXlFNYBE1oqJnRXR5QJeo+cYXY7N3TTjg6kh7v/8VCrM436+ZB1MPfJUm5xssrLMXY3VABJZKaMSu23btmHw4MGwsbExVTxUBTqFLbQKO+gUlbyMQmHypEyvcTZpeUbLzIbShP+sV/YL227YPTg6KWCXdrdS/YTY/6dmKJ2coHR2htLJqdztTJ6IsQmRiMyEUYndiBEjEBAQgKeffhqTJ09GSEiIqeKiSjjrNbRS29k46qGAgFAoUGhrusEKerUTMqeMMFl5JvHeKrhnmi6zq+zoxA8/KfYa8EteLg/t3VOp7Wpq1ClgmpGnRUzdTElEls+oxC4qKgrr16/HokWLsGTJEvTt2xfPP/88BgwYAIVCYaoYyUhBfbPu33rMVoMzjaLkDqduMtHoRPb/MTMmHHUKmG7kKRFRdRmV2C1fvhzvvfceNm7ciE8//RS7d+/Gnj17EBgYiMmTJ2Py5MkICAgwVaxEsjHV6ET2/zEv5jjqlIjIGEYPnnBwcMDEiRMxceJEXLhwAZ9++inWr1+P+fPnY9GiRRgwYACef/559O3b1xTxElWKa6beJHNiWVMndrJMTye8hxS96Uexeyo1WO33isnLJSLjmHRUbPPmzaVavE2bNuGNN97Azp07sXPnToSEhCAqKgpTp06Fs7OZdbivo0IzT8NWXwCt0h43XNrKHY5ZUYra7dC+Im07MvU5cFE6YbrbkFrbLxlKfPc96DLSYaNxhe+rTDoAIEWfgWRdutxhEFEtMfl0J9nZ2Vi/fj3+85//4NatWwCA1q1b48qVK3j11Vfx0UcfYdeuXWjdurWpd211vPKu3e87p3RiYvePdGcl9NBBCSW8bIzveJ6qy4QeAhmOmZhy+80yt0vRZUAPASUUOJhzpsztFuky4A4gNS/V6NiopIzdu6FNTIStry8TuwcooIAryh8tXBnpyOHdXYjMmMkSu19//RWffvopYmJikJmZCScnJ0yePBnTpk1D69atkZWVhZUrV2LevHmYOXMmjh07ZqpdE0nenqRGGrKhhAKeJpiGJ1mnAPDPQKBK1HroIcqtHdFDf78owXn2qGymbD5N0d0vxxVO+LfjU0aX91ruBqQh2+hyiKhmGJXY5eTkICYmBp9++inOnDkDIQSaNm2KKVOmYMKECdAUG6qvVqvx6quvIj4+Hl988YXRgcshOjoa0dHR0On4pWzuKkqwqsMNZXchKKrFqLhWhE1i5iA1LxUaAHd1GeXWxFaVqfqdsfmUiKrLqMQuICAAmZmZsLGxwfDhwzFt2jRERkaW+5rAwEDk5eUZs1vZREVFISoqChkZGXA14RQJZDoahSNqopVIo3DE/zmUPW9gUS1GxbUiq2CqAE09yW4Rc5w8ubLHWjTYpaIbyqszCgHcr0E15wTKVM2nwD+fDSKyeEYldi4uLvjXv/6FZ599Fn5+fpV6zbRp0zB69GhjdktUpvKSL3NiilG71nS3gypPKFzBDeWVxX4vrya2sopqbFN0GRhsghpAUzef1gRTHWsRjrIlMg2jErsbN25AqVRWvGExGo3GoImWyBqZetSuKe52UDR5ckW1XVVl0hrACiYUliaArmC7u7oM6KFHprONSfud1UQXAHNl7sdqyvexOdZiE5XFqMSuV69emDhxIsaPH1/mNl9//TVWr16Nw4cPG7MrIotg6lG7gOnudiBNnlxBbZecKppQuOgYKtpuyu03kaxLhxuc8W8TxFWTXQDMjamPtcZG2Zrx+5ioJhmV2MXGxlbYp+7GjRs4evSoMbshshhFo3a9bVyxI2Ch3OEYMPV9Sa3p9ml1pQuAKZj6WE09ytaU7+M6UYtN9ACTz2P3oNzcXNja1vhuiMhIpr7HKW+fRnIw5fu4LtRiEz3I6IxLoVCUulwIgbi4OOzZswfBwcHG7oZKcU8VBFt9PrRKldyhWL0WNsHIFvlwVvBayMmheXPos7Oh5N1tyARYi011UZUTO6VSaZDMzZ8/H/Pnzy9zeyEEXn/99WoFR+W7qukodwj0j7H2XeUOgQC4P2WeI0ipbmItNtVFVU7sHn30USmxO3bsGEJCQhAWFlZiOxsbG3h6eqJnz56YPHmy0YESEZHl4vQpRKZR5cQuNjZW+l2pVGLSpEmYO3euKWMisnim/BLjFxhZAnOfPoWorjCqj921a9fgZoL5s4isDb/EiO6rM9OnENURRiV2oaGhpoqDqqFlyrew0+eiUOmIc54D5Q7Hqi3O+wYZIrfCW4+Z8kuMX2AlJb7zDvQZGSab249qnrlPn0JU11QpsXvrrbegUCgQFRUFDw8PvPXWW5V6nUKhwJtvmq7vBN1np8+FSp8jdxgEIEPk3v8yqSDPMuWXGL/AStJnZLBzOpk9zotHNalKid38+fOhUCgwcuRIeHh4lDsatjgmdkTWyxRfYgVJiVDi/q3AppTTN3GRLgPuldiu6F6sRLLgvHhUg6qU2B05cgQAEBISYvCciKhMJvgSK7ojtR76cvsm6qGv1HZEcuC8eFQbqpTYdevWrdznRFT7zHWaCFN+id3VZUAPPdKcFXBD2ZMPK5EBQECJ8rcrYo73YiXLxXnxqDbwXl9EdZy5jrA15ZfYlNtvIlmXDjc449+OZU9C7Kr4EkA2XBVO5W5HRGSpjErsrl+/jkuXLqFbt25w/ucWPlqtFgsXLsT27dvh7OyMV155BUOHWs8NsolqS01NE2GuNYBERFQxoxK7BQsWYOfOnUgs1n9m0aJFWLhwofT8ySefxPHjx/HII48YsysiekBNTRNhrjWARERUMaMSuxMnTqBnz56wtb1fjF6vx8qVK9GkSRMcOHAACQkJ6NWrFz788ENs2rTJJAETUc3gRLFERHWfUYldYmKiwSTFZ8+exd27dzFv3jwEBQUhKCgIgwcPxvHjx40OlIhqFieKJSKq+4xK7AoLC6FQKKTnP/zwAxQKBXr06CEtCwoKwp07d4zZDZXhhjoCSqGFXsExMHIbbtceBdDCnuORZJXT+xEoCrUQdrwORGSdjPrrFxQUhHPnzknP9+zZAy8vLzRt2lRalpSUBI2J5+6h++461pc7BPpHe9sGcodgtkwxGKOyEwoXtmxk1H7IcphyEBAHAFFdYlRiN3DgQHz44Yd4+eWX4eDggIMHD2LSpEkG2/z555+8pyyRFeNgDJID33dkrYxK7F599VVs374dS5cuBQAEBgZiwYIF0vqkpCScOHECM2fONC5KMgl7bSYi/nzPZOUV2qpxrv5Uk5VHlsXUgzGkMonKYcr3HQcAUV1kVGLn4+OD8+fP49ChQwDu34nCxcVFWn/37l2899576NOnj3FRUqkctOlQQA8BJfJsXSvcXgEBlZb3yKwJCfo06KGHEkr4Kd3kDscsmHowRmUo794D9AJQKqD3cq/1/ZP8TPm+qysDgExxP+Yitl5eqLd1i0nKInkY3cPY0dERAwcOLHVdeHg4wsPDjd0FlaHZvQNQ6XOQr3TCGe8nytyuUGnaWg57fS4U/C/WwIf5e5CG7ArvjEA1y2XtLigzsqHXOCP95fFyh0MWwuwn7TbB/ZjJcljl0LH4+HiMGzcOSUlJsLW1xZtvvoknnig7MarrznmWnnhXV0TyZqj0OSYtk4jIXJlrfz1T3o9Zl54OCP7DbgmMTuxSU1OxevVqnDx5Evfu3YNOpyuxjUKhkJprzYGtrS2WLVuG1q1bIyEhAREREejfv790WzQiIiJzn7TblPdjvvP669ClpZmsvGvDR0B7967JyivCpuKKGZXY/f7774iMjERycjJEOZl+8bnuzIG/vz/8/f0BAH5+fvDy8kJqaioTOyIiknDS7urT3r3L5mGZKI158csvv4ykpCS89tpruHr1KgoLC6HX60s8SqvFK8+xY8cwaNAgBAQEQKFQYPv27SW2iY6ORlhYGBwcHNChQwecPHmyWsdw5swZ6HQ6BAcHV+v1REREVAaFAjZubkY/YGYVRObMqBq748ePY8CAAXj77bdNFQ8AIDs7G61atcLTTz+NYcOGlVi/adMmzJo1C6tWrUKHDh2wbNky9OnTB3/88Qd8fHwAAK1bt4ZWqy3x2gMHDiAgIADA/Wbk8ePH47PPPjNp/ERERATYuLrC3wQ5gqmbii2ZUYmdEKJGRr3269cP/fr1K3P90qVL8eyzz0qTIa9atQq7d+/G6tWrMfufPgdnz54tdx/5+fkYMmQIZs+ejU6dOlW4bX5+vvQ8I4NThhARUfWY/ShbqtOMSuwiIiLwxx9/mCqWSikoKMCZM2cwZ84caZlSqUSvXr1w4sSJSpUhhMDEiRPRo0cPjBs3rsLtlyxZYjDxMhERUXWZ6yhbwHRz4mmTk40PhqrFqMRu7ty56NOnD2JjYxEZGWmikMp39+5d6HQ6+Pr6Giz39fXF77//XqkyfvjhB2zatAktW7aU+u999dVXaNGiRanbz5kzB7NmzZKeZ2RksE8eERFVibmPsgXAOfEsgFGJXXx8PAYPHozevXtj9OjRiIiIgJubW6nbjh9vPpOFdunSBXq9vtLbq1QqqFSqGoyIiIgsnTmPsjXlnHi1US6VzajEbuLEiVAoFBBC4KuvvsJXX31VYmoTIQQUCoXJEjsvLy/Y2Ngg8YH/KBITE+Hn52eSfdQV5zwGQAEBAY4WktvrDkOkW4qRfDKeHy7dUoyIKs+Uc+KRvIxK7NasWWOqOCrN3t4eEREROHToEIYMGQIA0Ov1OHToEKZPn17r8cip0MZJ7hDoH64KXgtzIFw4FyURWTejErsJEyaYKg4DWVlZuHLlivT82rVrOHv2LDw8PBASEoJZs2ZhwoQJaNu2Ldq3b49ly5YhOztbGiVLREREZI3M8l6xp0+fRvfu3aXnRQMXJkyYgLVr12LkyJFITk7G3LlzkZCQgNatW2Pfvn0lBlSYWnR0NKKjo6s84TIREREZz1SjdotY4i3KTJLYJScnY+vWrbh8+TKys7Px+eefS8uvXbuGFi1awNHRsdLlRUZGlnuLMgCYPn16rTe9RkVFISoqChkZGXB1da3VfZfGN+dPKEUh9Ao7JDo1kjscq3ZMexn5ohAqhR0etW0qdzhWy/70JSjyCyFUdihoa/o5NolIZhy1WyGjE7svvvgCM2fORF5enjRQoiixS0xMRMeOHfGf//wHzzzzjNHBkqGg7N+g0ucgX+nExE5muwt/RRqy4QZnJnYycow9DWVGNvQaZyZ2RBbE1KNrdenpQAUVSHWVUYndwYMH8dxzz6Fly5ZYsGAB9u/fj1WrVknrmzdvjmbNmmH79u1M7IiIiKhaTD1qt+gWZaZs2jWXZl2jErt///vf8Pf3x9GjR6HRaPDrr7+W2KZly5aVviMEERERUa2xwKZdoxK706dPY9SoUdCUU0UaFBSEhIQEY3ZjNjh4goiIqO4zZdOuuTXrGpXYFRQUwNm5/Hmj0tLSYGNjY8xuzIa5DZ4gIiKiqjNl025Rs665MGqa/LCwMJw5c6bcbX7++Wc0btzYmN0QERERUSUYldgNHjwYx48fx+bNm0tdv2bNGpw7dw7Dhw83ZjdEREREVAlGNcW++uqr2LhxI0aPHo0tW7YgPT0dALBixQocP34c27ZtQ8OGDTFjxgyTBEtERESGUnQZGHz7TZOU5anUYLXfKyYpi+RhVGLn7u6Oo0ePYvz48Qa1djNnzgQAdO3aFRs2bKiwHx4RERFVjx4Cybp0ucMgM2H0BMUhISGIjY3FuXPncOLECaSkpMDV1RWPPPIIIiIiTBEjlSHPRgOdwg6Fysrf1YNqhq9SA0dhB42C10JOOk9XCJU99GpeB7J8GoUjYKLBmOnIgTBVYSQrk90rtmXLlmjZsqWpijNL5jbdyUWPPrLu316biYg/3zNJWYW2apyrP9UkZclhlmqg3CEQgKxJg+UOgajW/J/DUJOV9VruBqQh22TlkXxMktjduHEDycnJUCgU8Pb2RkhIiCmKNTuc7sSQAgIqbYbcYRAREdE/qp3Y3b17F2+//TZiYmKQlJRksM7X1xdjxozBnDlz4OHhYXSQZF5M2fRrr8+FgtX/REREJlGtxO6vv/7CY489hvj4eAghYGtrC09PTwghkJqaioSEBCxduhRbt27Fd999h/r165s6bpLROU/TNTtGJG+GSp9jsvKIiEh+Tye8hxS96Vt0OGq3YlVO7PR6PcaMGYO4uDhERkbijTfeQJcuXWBvbw8AyM/Px/Hjx7F48WIcPXoUY8eOxY8//mjywAlomH4Mtvp8aJUq/OX6qNzhWLUvCg4jS+RBrXDAM/Y95A7Hajlt+Q7K7DzonR2QM6KX3OEQWa0UfQZH6sqkyondgQMHcPr0aTz55JOIiYmBQqEwWK9SqdCrVy/07NkTI0eOxNatW3Hw4EE89thjJgua7tMUJEKlz0G+0knuUKzen7oEpCEbbuDUPnKyu34byoxs6DW8DkTmQAEFXGH8dxRH7VZelRO7rVu3QqVSYfny5SWSuuIUCgVWrFiBnTt3YsuWLUzsiIiIrIwrnPBvx6eMLoejdiuvyrcU++WXX9C5c2d4e3tXuK2Pjw+6dOmCX375pVrBEREREVHlVbnGLj4+Hl26dKn09s2aNUNMTExVd2OWzG0eOyIiIlMy1e3JUnScCksuVU7sMjIy4ObmVunt3dzckJmZWdXdmCXOY0dERJaMtyer+6qc2BUUFMDGxqbS2yuVShQUFFR1N0RERFRLTHl7shLlUq2q1jx25Q2aICIiorrFlLcnI3lVK7GbP38+5s+fb+JQiIiIiMgY1UrshKhafS1r+IiIiIhqXrXuPEHmIdGxIWxFIbQKO7lDsXpdbBsjVxTAUWEvdyhWLT+iKRR5BRAOvA5EZJ2qVWNH5uGmurXcIdA/BtlFyB0CAcjr3k7uEIiIZMXEjsyCvTYTEX++Z7LyCm3VOFd/qsnKIyIiqguY2FUBJyiuOQoIqLSc0JKIiMgYTOyqgBMUm16h0rRzHNnrc6HgjaKJiMhKMbGrwyKSN0Olz0G+0glnvJ+QO5xqOec50KTlFZ2T2lZ0g2o3OJvkhtdUPa7vfwllRjb0Gmekvzxe7nCIyMRMdcuzIp5KDVb7vWKy8swBEzsiIiKqE3jLs4oxsSMiIiKzZupbnqUjB8JCu+0wsSMiIiKzZupbnhV1n7FETOyIiIjIKpmiz94iXQbcAaTmpZomKCMxsSMiIiKrZIo+e3rcvyOXTpjHVGhM7IiIiMiqmLbPnnkN5mBiR0RERFbFtH32VsGkIzuMpJQ7ACIiIiIyDSZ2VRAdHY3w8HC0a8cbjRMREZH5YVNsFZjbLcX+cu0KpdBBr7CROxSr97R9JLTQwRa8FnLKHt4T0OoAW14HIrJOTOzqsAx7P7lDoH80tgmQOwQCoK0XKHcIRESyYlMsERERkYVgYkdERERkIdgUW4dpChKkPnZslpXXH7rbUh87NsvKx/baLamPHZtlicgaMbGrwxqmH4dKn4N8pRPOeD8hdzhWbXVBLNKQDTc449+OT8kdjtVy3noIyoxs6DXOSH95vNzhEBHVOjbFEhEREVkIJnZEREREFoKJHREREZGFYB87onK0vPoJ7LRZFW5n768BbJWw12Yi4s/3Kty+0FaNc/WnmiJEIiIiCRM7onLYabOg0mZUYks17leA6yu5PRERkekxsSOqBAEFCpSO5WyhkH7mK53K3MpenwsFhEljIyIiKsLEjixSZZtEK1MOABQoHcudUqZA+SOAfBQonXDGu1eZ20Ukb4ZKn2N0XERERKVhYlcF0dHRiI6Ohk6nkzsUqoACgk2iRERkdZjYVUFUVBSioqKQkZEBV1dXucOhUhSW21xqfuUSERGZEhO7Oox3myjpnOdAWfb777xOsuyXDPFuE0Rk7TiPHREREZGFYGJHREREZCHYFEskA1ON2gU42TEREf0PE7s6LCjrLGxFIbQKO9xUt5Y7HKu2y/YachVaOApbDNLWq3B7jtqtGQ5HTkGRVwDhYI+87u3kDoeIqNYxsavDfHP/gkqfg3ylExM7mX1vewdpyny46VXlJnamHF3LyY5LUp25DGVGNvQaZyZ2RGSVmNgR1SJTjtrlZMdERPQgDp4gIiIishBM7IiIiIgsBBM7IiIiIgvBxI6IiIjIQjCxIyIiIrIQTOyIiIiILAQTOyIiIiILwXns6rAMe1/Y6vOhVarkDsXqNdK7IUsUQC3s5Q7FqhWGBUCZnQe9s4PcoRARyYKJXR32l+ujcodA/3imIFzuEAhAzohecodARCQrNsUSERERWQgmdkREREQWgokdERERkYVgH7sqiI6ORnR0NHQ6ndyhAACape6HnT4XhUpHXPToI3c4Vm2p6ldkKAqgEfaYlf+w3OFYLfWaHVBm5UKvdkTWpMFyh0NEVOtYY1cFUVFRuHTpEk6dOiV3KAAAB10GnHTpcNBlyB2K1UtU5OKOMgeJily5Q7FqNinpsEm+B5uUdLlDISKSBRM7IiIiIgvBxI6IiIjIQjCxIyIiIrIQTOyIiIiILARHxRLVcfbaTET8+Z7cYZSp0FaNc/Wnyh0GEZFVYGJHVMcpIKDScmQ0ERExsSOqswqVjnKHUC57fS4UEHKHQURkVZjYEdVR5zwHyh1CuSKSN0Olz5E7DCIiq8LErg676dwKSlEIvcJO7lCs3oDCMOQrtFAJfqTklBvZFor8QghV3f5MtLz6Cey0WSYtk30diawDv4XqsESnRnKHQP94VBcgdwgEoKBtuNwhmISdNov9JomoWpjYERGZKQEFCozsS8m+jkTWhYkdEZGZKlA64oz3E0aVwb6ORNaFiV0dZqfLgQICAgoU2jjJHY5VS0c+9AoBpVDAFSq5w7FaisxsQC8ApQLCxVnucIiIah0TuzqsZepuqPQ5yFc6Gf1fPRnnbYczSFPmw02vwr/zOskdjsWp7GCCG9ucoMtRwsZJj9BhFddScUABEVkaJnZEZPYqO5hAIRz/+clJm4nIOjGxI6I6o6LBBAIK6We+suzuCRxQQESWiokdEdUZFQ0mcLU5BiXyUWjjiDPefcrcrmhAganvs8umXSKSGxM7IqpRpkie7LWZJorGEO+zS0SWhokdEdUoc0yeTH2fXTbtEpG5YGJHRDXC1MmTKcs09X12Td20W1M1lERk+ZjYEVGNMHXyVBeYY+0kEVkXJnZEREaqidrJmiyXiCwXEzsiIiNZY+0kEZknJnZ12EX33lBADwGl3KFYvZfyW/9zJRRyh2LVModFSLcUIyKyRkzs6rA8W1e5Q6B/+Aneq9cc6N15f1gism6s6iEiIiKyEFaZ2KWlpaFt27Zo3bo1mjdvjs8++0zukIiIiIiMZpVNsS4uLjh27BicnJyQnZ2N5s2bY9iwYfD09JQ7tCrxyr0KpdBCr7DFXcf6codj1U7aJKIAOtjDBu11vnKHY7Xsfr8DhVYHYWuDwib+codDRFTrrDKxs7GxgZPT/T5R+fn5EEJAiLo3a3xo1hmo9DnIVzoxsZPZVru/kabMh5texcRORk4//AVlVj70ahXSmdgRkRUyy6bYY8eOYdCgQQgICIBCocD27dtLbBMdHY2wsDA4ODigQ4cOOHnyZJX2kZaWhlatWiEoKAivvPIKvLy8TBQ9ERERkTzMssYuOzsbrVq1wtNPP41hw4aVWL9p0ybMmjULq1atQocOHbBs2TL06dMHf/zxB3x8fAAArVu3hlarLfHaAwcOICAgAG5ubvjtt9+QmJiIYcOGYcSIEfD1ZU0LEVF5Wl79BHbaLJOXW2irxrn6U01eLpG1McvErl+/fujXr1+Z65cuXYpnn30WkyZNAgCsWrUKu3fvxurVqzF79mwAwNmzZyu1L19fX7Rq1QrHjx/HiBEjSt0mPz8f+fn50vOMDN4yiIisk502i7dNIzJjZpnYlaegoABnzpzBnDlzpGVKpRK9evXCiRMnKlVGYmIinJyc4OLigvT0dBw7dgxTp5b9n+KSJUuwYMECo2MnIpKLvTYTEX++Z5JyAEBAgQIT3PLMXp8LBepeH2cic1XnEru7d+9Cp9OVaDb19fXF77//Xqkybty4geeee04aNDFjxgy0aNGizO3nzJmDWbNmSc8zMjIQHBxcvQMgIpKBAsKkNW0FSkec8X7C6HIikjdDpc8xQUREBNTBxM4U2rdvX+mmWgBQqVRQqVQ1FxARUQ0pNEGtWm2Ua6oaxSLss0fWqs4ldl5eXrCxsUFiYqLB8sTERPj5+ckUFRGReTrnOVDuECrF1DWKRNbKLKc7KY+9vT0iIiJw6NAhaZler8ehQ4fQsWNHGSMjIqKqKlQ6Il/pZLKHgELuQyKSlVnW2GVlZeHKlSvS82vXruHs2bPw8PBASEgIZs2ahQkTJqBt27Zo3749li1bhuzsbGmUbE2Jjo5GdHQ0dDpdje6nsoqaQmqqqYUqTyPsAf0/P0k2eid7g59k/kxdo8g+e2TtzDKxO336NLp37y49Lxq4MGHCBKxduxYjR45EcnIy5s6di4SEBLRu3Rr79u2r8XnooqKiEBUVhYyMDLi6utboviqjrjSxWIP/y28rdwgEIHP0I3KHQEQkK7NM7CIjIyu8xdf06dMxffr0WoqIiIiIyPzVuT52RERERFQ6JnZEREREFsIsm2KpcupnnICtPh9apQpXNRwRLKev7f5AtqIQzsIOYwsbyx2O1XI6dAmKvEIIBzvk9AyXOxwiolrHxK4KzG1UrHv+Taj0OchXOskditU7b5OCNGU+3PQqoFDuaKyX3fW7UGblQ6/mhOJEZJ3YFFsFUVFRuHTpEk6dOiV3KEREREQlMLEjIiIishBM7IiIiIgsBPvYERERlaLl1U9gp80yebmFtmqcqz/V5OUSAUzsiIiISmWnzYJKmyF3GERVwsSOiIioHAIKFJjgntz2+lwoUP5dlYiMxcSuCsxtuhMiIqp5BUpHnPF+wuhyIpI3Q6XPMUFERGXj4Ikq4HQnREREZM5YY1eH3XWoB1t9AbRKe7lDsXrtdD7I0WvhJPiRklNBIz8o8gshVHZyh0JEJAt+C9VhN1zayh0C/WNEYQO5QyAAuV0byR0CEZGs2BRLREREZCGY2BERERFZCCZ2RERERBaCfezqsNZ3v4G9PhcFSkec9RoqdzhWba7Dz0hX5MNVqPBWXge5w7Fami9/gDI7H3pnFTLGd5Y7HCKiWscauyqIjo5GeHg42rVrJ3coAAAboYWtKISN0ModitXLhw55Ch3ywTkO5aQo1EFRoIOikNeBiKwTE7sq4Dx2REREZM6Y2BERERFZCPaxIyIii2OvzUTEn+8ZXQZRXcPEjoiILI4CAipthtxhENU6JnZERGQxCpWOdaJMoprCxI6IiCzGOc+BcodAJCsmdkRERGRSLa9+AjttlsnLLbRV41z9qSYv15IwsSMiIiKTstNmsY+jTJjYVUF0dDSio6Oh05nH5KdXNY9AKXTQK2zkDsXqjSlohEKFHnaCMwjJKad7U0CnB2x4HYjMgYACBSboo2ivz4UCwgQRWT4mdlUQFRWFqKgoZGRkwNXVVe5wcE8VLHcI9I+Wei+5QyAAhfW95Q6BiIopUDrijPcTRpcTkbwZKn2OCSKyfPy3loiIiMhCMLEjIiIishBsiq3DnAtToBA6CIUNsu085Q7Hqt1QZEKr0MNWKBEqXOQOx2rZJGYAej2gVELnq5E7HCKiWsfErg5rknYYKn0O8pVOJunDQNW3UnUeacp8uOlV+HdeJ7nDsVrqb89CmZUPvVqF9GcelTscIqJax6ZYIiIiIgvBxI6IiIjIQjCxIyIiIrIQTOyIiIiILAQTOyIiIiILwcSuCqKjoxEeHo527drJHQoRERFRCZzupArM7ZZiRERkvVpe/QR22iyTllloq8a5+lNNWqY5MuW5OyPU938RepOUZywmdkRERHWQnTYLKm2G3GHUSaY9d0WJnYmKMxITOyIiojpMQIECpaNRZdjrc6Ewl8ykFpni3JkbJnZ12K+eQ6CAgIBC7lCs3oK89hAAr4TM0sd1Ai8EWZsCpaPRdx+KSN4MlT7HRBHVHaY4d8BBk8RiKkzs6jC90k7uEOgfDvwomQd7Xgcism4cFUtERERkIfjvLREREQEA7LWZiPjzPZOUUxPMPT5zwMSuDvPPvghbUQitwg53nJvJHY5VO2gbhzyFDg7CBo9pQ+QOx2qpfrkBRYEWwt4W+W1C5Q6HqM5RQJj1SFtzj88cMLGrwwJyLkGlz0G+0omJncy+s72JNGU+3PQqJnYycvj1BpRZ+dCrVUzsiKqgsIZGhpqqXHOPz5wwsSMiIrJy5zwHyh1Cucw9PnPCwRNEREREFoKJHREREZGFYFMsERFRLeLITqpJTOyqIDo6GtHR0dDpdHKHQkREdRRHdlJNYmJXBVFRUYiKikJGRgZcXV3lDoeIiOoQjuyk2sDEjoiIqBZwZCfVBg6eICIiIrIQrLGrw7JtPVAgnFGoUMkditUL0avhLlRwEXZyh2LVtN4uUKodoHfkdSAi68TErg773b2n3CHQP6IKWsodAgHIfvxhuUMgIpIVm2KJiIiILAQTOyIiIiILwcSOiIiIyEKwj10d1uTeIdiJfBQqVOxvJ7No+3PIVBTCRdixv52MnHf+CmVuIfSOduxvR0RWiYldHeasTYVKn4N8pZPcoVi9OGUW0pT5cNNzhLKcbJMzoczKh17N60BE1olNsUREREQWgokdERERkYVgYkdERERkIZjYEREREVkIJnZEREREFoKJHREREZGFYGJHREREZCE4j101CCEAABkZGTVSfrZOBzsdkK0DkFtQ5nYZeXqohEC+Qo/scrajmqcVOuiUOmj1OmTn8VrIxUang1Kng16n42eCiGpHse/smsoLisotyj/KoxCV2YoM3Lx5E8HBwXKHQURERFYkPj4eQUFB5W7DxK4a9Ho9bt++DRcXFygUCpOXn5GRgeDgYMTHx0Oj0Zi8fKo8XgvzwOtgPngtzAevhXmojesghEBmZiYCAgKgVJbfi45NsdWgVCorzJhNQaPR8MNqJngtzAOvg/ngtTAfvBbmoaavg6ura6W24+AJIiIiIgvBxI6IiIjIQjCxM0MqlQrz5s2DSqWSOxSrx2thHngdzAevhfngtTAP5nYdOHiCiIiIyEKwxo6IiIjIQjCxIyIiIrIQTOyIiIiILAQTOzMUHR2NsLAwODg4oEOHDjh58qTcIVm8Y8eOYdCgQQgICIBCocD27dsN1gshMHfuXPj7+8PR0RG9evXCX3/9JU+wFmzJkiVo164dXFxc4OPjgyFDhuCPP/4w2CYvLw9RUVHw9PSEWq3G8OHDkZiYKFPElumTTz5By5YtpXm5OnbsiL1790rreQ3k884770ChUODFF1+UlvF61I758+dDoVAYPJo0aSKtN5frwMTOzGzatAmzZs3CvHnz8Msvv6BVq1bo06cPkpKS5A7NomVnZ6NVq1aIjo4udf27776Ljz/+GKtWrcLPP/8MZ2dn9OnTB3l5ebUcqWU7evQooqKi8NNPP+HgwYMoLCxE7969kZ2dLW3z0ksvYdeuXdi8eTOOHj2K27dvY9iwYTJGbXmCgoLwzjvv4MyZMzh9+jR69OiBwYMH4+LFiwB4DeRy6tQpfPrpp2jZsqXBcl6P2tOsWTPcuXNHenz//ffSOrO5DoLMSvv27UVUVJT0XKfTiYCAALFkyRIZo7IuAMQ333wjPdfr9cLPz0+899570rK0tDShUqlETEyMDBFaj6SkJAFAHD16VAhx/7zb2dmJzZs3S9tcvnxZABAnTpyQK0yr4O7uLj7//HNeA5lkZmaKhg0bioMHD4pu3bqJF154QQjBz0RtmjdvnmjVqlWp68zpOrDGzowUFBTgzJkz6NWrl7RMqVSiV69eOHHihIyRWbdr164hISHB4Lq4urqiQ4cOvC41LD09HQDg4eEBADhz5gwKCwsNrkWTJk0QEhLCa1FDdDodNm7ciOzsbHTs2JHXQCZRUVEYMGCAwXkH+JmobX/99RcCAgJQv359jBkzBnFxcQDM6zrwXrFm5O7du9DpdPD19TVY7uvri99//12mqCghIQEASr0uRevI9PR6PV588UV07twZzZs3B3D/Wtjb28PNzc1gW14L0zt//jw6duyIvLw8qNVqfPPNNwgPD8fZs2d5DWrZxo0b8csvv+DUqVMl1vEzUXs6dOiAtWvXonHjxrhz5w4WLFiArl274sKFC2Z1HZjYEZFZioqKwoULFwz6sFDtady4Mc6ePYv09HRs2bIFEyZMwNGjR+UOy+rEx8fjhRdewMGDB+Hg4CB3OFatX79+0u8tW7ZEhw4dEBoaiv/+979wdHSUMTJDbIo1I15eXrCxsSkxiiYxMRF+fn4yRUVF557XpfZMnz4d3377LY4cOYKgoCBpuZ+fHwoKCpCWlmawPa+F6dnb26NBgwaIiIjAkiVL0KpVK3z00Ue8BrXszJkzSEpKQps2bWBrawtbW1scPXoUH3/8MWxtbeHr68vrIRM3Nzc0atQIV65cMavPBRM7M2Jvb4+IiAgcOnRIWqbX63Ho0CF07NhRxsisW7169eDn52dwXTIyMvDzzz/zupiYEALTp0/HN998g8OHD6NevXoG6yMiImBnZ2dwLf744w/ExcXxWtQwvV6P/Px8XoNa1rNnT5w/fx5nz56VHm3btsWYMWOk33k95JGVlYW///4b/v7+ZvW5YFOsmZk1axYmTJiAtm3bon379li2bBmys7MxadIkuUOzaFlZWbhy5Yr0/Nq1azh79iw8PDwQEhKCF198EYsWLULDhg1Rr149vPnmmwgICMCQIUPkC9oCRUVFYcOGDdixYwdcXFykvimurq5wdHSEq6srnnnmGcyaNQseHh7QaDSYMWMGOnbsiEceeUTm6C3HnDlz0K9fP4SEhCAzMxMbNmxAbGws9u/fz2tQy1xcXKQ+pkWcnZ3h6ekpLef1qB0vv/wyBg0ahNDQUNy+fRvz5s2DjY0NRo8ebV6fi1odg0uVsnz5chESEiLs7e1F+/btxU8//SR3SBbvyJEjAkCJx4QJE4QQ96c8efPNN4Wvr69QqVSiZ8+e4o8//pA3aAtU2jUAINasWSNtk5ubK6ZNmybc3d2Fk5OTGDp0qLhz5458QVugp59+WoSGhgp7e3vh7e0tevbsKQ4cOCCt5zWQV/HpToTg9agtI0eOFP7+/sLe3l4EBgaKkSNHiitXrkjrzeU6KIQQonZTSSIiIiKqCexjR0RERGQhmNgRERERWQgmdkREREQWgokdERERkYVgYkdERERkIZjYEREREVkIJnZEREREFoKJHREREZGFYGJHZKFOnz6Nxx57DF5eXlAoFGjdurXcIZmFv/76C0OHDoWfnx8UCgXc3NwAAPPnz4dCoUBsbGyly4qMjIRCoaiZQM1AWFgYwsLC5A6jhNjYWCgUCsyfP7/EurLe9xMnToRCocD169drNdbymOv5pbqN94olq/Tgl7G9vT00Gg2Cg4PRpk0bDB8+HL1794aNjY1MERonIyMDAwYMQF5eHsaNGwcvLy/4+fnJHZbsdDodhgwZgitXrmDcuHEICgqCg4OD3GFV2fXr11GvXj1MmDABa9eulTscs2Fu7/vIyEgcPXoUvMET1SYmdmTV5s2bB+D+F35aWhouXryIr776Cl988QXatm2L9evXo1GjRjJHWXUnT55EUlISFi9ejNdff13ucMzGtWvXcOnSJTz77LP4z3/+Y7Bu+vTpGDVqFEJCQmSKjiqrffv2uHz5Mry8vAyWl/e+X7JkCWbPno3AwMDaDLVchw4dkjsEskBM7MiqldaUk5iYiBkzZmDz5s3o1asXTp8+DR8fn9oPzgi3b98GAAQEBMgciXkp77x4eXmVSBTIPDk5OaFJkyYllpd3ff39/eHv71/jsVXFQw89JHcIZIkEkRUCIMp7++t0OhEZGSkAiBdeeMFg3enTp8XMmTNFy5Ythbu7u1CpVKJBgwZi1qxZIjU11WDbVatWCQBi/vz5pe7nzp07wtbWVjRv3lxalpGRId566y3RrFkz4eLiItRqtahfv7548sknxenTp8s9rmvXrknH9uBjzZo10na3b98W06ZNE6GhocLOzk54eXmJoUOHllr+mjVrpNfv3btXdOvWTWg0mnLPX3HZ2dninXfeEREREUKtVgtnZ2fRpEkTMWPGDJGQkGCwbXXjOnz4sOjWrZtQq9XCxcVF9O/fX1y6dMlg+7LOy7x584QQQsybN08AEEeOHCmxr5iYGNGmTRvh4OAgvL29xdixY8WtW7dEt27dyjwP+/btE/369ROenp7C3t5e1K9fX7z88svi3r17JbYNDQ0VoaGhIisrS7z88ssiODhY2Nvbi4ceeki88847Qq/XS9sWxVnRNS5PfHy8mDFjhmjQoIFwcHAQ7u7uol27duKtt94qNa7i0tLSxLvvviu6d+8uAgMDpes0aNAg8eOPP5a6v2PHjomBAweKwMBAYW9vL3x9fUWHDh1KfC4SEhLEv/71L9GoUSPh5OQkXF1dRaNGjcSECRPE33//LW135MgRg2tXmff9hAkTBABx7dq1EvH9/PPP4sknnxQBAQHC3t5e+Pn5iccee0xs2rTJYLs1a9aIYcOGiXr16gkHBwfh4uIiOnXqJL766iuD7cqLp1u3buWeXyGEyMvLE0uWLBHNmzcXjo6OwsXFRXTp0qVEPMX3NWHCBHHt2jUxcuRI4enpKVQqlYiIiBC7du0q9ZqQ5WKNHVEplEol3njjDcTGxiImJgYffvih1C/vs88+wzfffINu3bqhV69e0Ov1OHPmDJYuXYq9e/fi559/houLCwBgzJgxePXVV/HFF1/gjTfeKNFnb/Xq1dBqtXj++ecBAEII9O3bFz/++CM6duyIyZMnw9bWFjdv3sSRI0fQtWtXRERElBm3m5sb5s2bh7Nnz2LHjh0YPHiw1Hm86Oe1a9fQpUsX3L59Gz169MDo0aMRHx+PzZs3Y/fu3di6dSsGDhxYouwtW7Zg37596NevH6ZMmYIbN25UeB7v3buH7t2747fffkPjxo3x9NNPw97eHn///TfWrFmDYcOGwdfX16i4vv32W+zYsUOK69KlS9izZw9OnTqFS5cuSbVw8+bNw/Xr17Fu3Tp069YNkZGRACD9LMuHH36IWbNmwc3NDePHj4ebmxv279+PTp06wdXVtdTXLFiwAPPnz4eHhwcGDhwIHx8fnDt3Du+//z727NmDEydOQKPRGLymsLAQffr0we3bt9GvXz/Y2tpi+/btmD17NvLy8qRuA5GRkUhLS8NHH32EVq1aYciQIVIZlRkgc/r0afTp0wepqal49NFHMWzYMOTk5ODSpUuYP38+3nzzzXJff/nyZfzf//0fHn30UQwYMADu7u6Ii4vDzp07sXfvXuzatQt9+/aVtt+3bx8GDBgAjUaDxx9/HIGBgUhNTcXly5excuVK6bhycnLQuXNn/P3333jssccwaNAgCCFw48YN7NixAyNGjED9+vVLjaky7/uyfPbZZ5g6dSpsbGzw+OOPo2HDhkhKSsLp06excuVKPPnkk9K2U6dORbNmzfDoo4/C398fKSkp2LNnD8aNG4c//vgDCxcuNIhn7dq1uHHjhnSMACocLFFQUIA+ffrg6NGjaNKkCaKiopCTk4MtW7Zg5MiROHv2LN5+++0Sr7tx4wbat2+P+vXrY9y4cUhNTcWmTZswePBgfPfdd+jevXu5+yULIndmSSQHVFBjJ8T9/5ptbW0FAHH16lVp+fXr14VWqy2x/eeffy4AiHfeecdgeVRUlABQ4j9nvV4v6tWrJ5ycnERaWpoQQohz584JAGLIkCElytfpdCVqBMtSvDbrQb179xYAxKJFiwyW//DDD8LGxkZ4eHiIzMzMEmUpFAqxd+/eSu2/yOjRowUAMWXKFKHT6QzWZWZmSsdtTFw2Njbiu+++M3jN7NmzBQDx73//22D5gzU9xZVWY3ft2jVhZ2cn3N3dDWp6dDqdGDZsWKnvo8OHDwsAomPHjiVq54pifvHFFw2Wh4aGCgCiX79+IicnR1qemJgoXF1dhaurqygoKDCIC//U0lRFfn6+CAsLEwDE+vXrS6yPj48vEVdpNXbJycmlvtbf3180adLEYHnReTp79myJ1xQvZ+fOnaWem6K4MzIypOdlXcfy3vel1dhdvHhR2NraCnd3d3HhwoVSj6m4K1eulBpbjx49hK2trbh586bBuvJqdIUo/fy+/fbb0nuhsLBQWp6YmCi9T3744QdpefHawQdrQPft2yeVRdaD050QlUGlUsHT0xMAkJycLC0PDQ0tdbTs008/DY1Gg/379xssnzp1KgDg008/NVh+4MABXLt2DSNHjixR8+Po6FiifKVSCXd39+odzD9u3ryJAwcOICQkBK+++qrBuk6dOmH06NFITU3Ftm3bSrx28ODBBjUxFUlKSsKmTZvg7++P999/H0ql4Z8btVotHbcxcY0aNQo9e/Y0WPbcc88BuN+Z3hjr169HYWEhZsyYYVDTolQq8d5775U4JgD4+OOPAdyvCSqaSqXIxIkT0bp1a6xfv77U/X388ccG197HxweDBw9Geno6/vjjD6OOBQB27dqF69ev4/HHH8dTTz1VYn1QUFCFZbi6upbaFzEoKAgjRozA77//jri4uBLrS3tPl1ZOadvZ29tLteCm9Mknn0Cr1eLNN99Es2bNSqx/8HyU1ifO3t4eUVFR0Gq1JhkMsXr1aigUCixduhS2tv9rVPPx8ZFqUz///PMSrwsNDcUbb7xhsKxPnz4ICQkx+nNAdQsTO6JyiH+mKSg+PUphYSFWrFiBLl26wMPDAzY2NlAoFFAqlcjIyMCtW7cMyihqutm7dy/i4+Ol5UWjMqdMmSItCw8PR+vWrRETE4POnTvj3XffxY8//oiCggKTHM+vv/4KAOjatSvs7OxKrO/Ro4fBdsW1b9++Svs6deoU9Ho9Hn30UTg7O9dYXG3bti2xLDg4GMD9pmBj/PLLLwCAbt26lVhXv359aT/FnThxAnZ2dti8eTPmz59f4lFQUIDk5GSkpKQYvM7V1RUNGjSosWMBgJ9++gkA0K9fP6PK+eGHH/Dkk08iODgYKpUKCoUCCoUCy5cvBwCDz8CYMWMAAB06dMCUKVOwadMm3Lx5s0SZ3bp1Q2BgIN555x307dsXH3/8Mc6cOQOdTmdUrOWp6vmIi4tDVFQUmjRpAicnJ+m4hw8fDgAlPvtVlZmZiStXriAgIKDUwSHlfQ5at25d6j+cwcHBJnnvUN3BPnZEZcjLy0NqaioAwNvbW1o+cuRIfPPNN6hfvz4GDx4MPz8/qFQqAMCyZcuQn59foqxp06bh2LFj+Pzzz7FgwQIkJCRg586daN26tUHCZGNjg8OHD+Ott97Cli1b8NprrwEAXFxcMGHCBCxZsgRqtbrax5Seng4AZY4OLFqelpZWYl1V5wMrKqMy00sYE9eDtWIApJoOY5OCoriK+gE+yM/Pr0Rfw5SUFGi1WixYsKDcsrOysqQaYaD04wBMdyxA1a5JWb755huMGDECDg4OeOyxx/DQQw/B2dkZSqUSsbGxOHr0qMFnYNiwYfj222/xwQcfYPXq1VLNdUREBJYsWYLHHnsMAKDRaPDTTz9h3rx52Llzp1Tz7eXlhWnTpuGNN94oNek3RlXOx9WrV9G+fXvcu3cPXbt2Re/eveHq6gobGxup72Zpn/2qMPXnALj//tHr9UbFRXULEzuiMnz//ffQarXw9fWVmuFOnz6Nb775Br169cLevXsNmkr0ej3efffdUssqGiTwxRdfYO7cuSUGTRTn7u6ODz/8EB9++CGuXLmCo0eP4tNPP8WKFSuQlpaGr776qtrHVNT0mZCQUOr6O3fuGGxXXFXvsFD0RVOZWgxj4qpJRftLTEwstamutHhdXV2h1+ulfwrMSVWuSVnefPNN2Nvb4/Tp02jatKnBuueffx5Hjx4t8ZoBAwZgwIAByM7Oxs8//4xvv/0Wn3zyCQYOHIhff/0V4eHhAO43fX7xxRcQQuDSpUs4fPgwoqOj8dZbb0Gv10uDE0yl+PkorYasuKVLlyIlJQVr1qzBxIkTDdbFxMRg3bp1Rsdjrp8DqlvYFEtUCr1ej8WLFwOAQV+kK1euAAAef/xxg6QOuN+fKzc3t9Ty7OzsMHnyZNy6dQu7du3C559/DrVaLTVTlaVBgwZ45plncPToUajVauzYscOYw8LDDz8M4H9J64OOHDkCAGjTpo1R+wHuN90qlUocO3YM2dnZZhNXVRTtr7Rk5erVqwZN60UeeeQR3Lt3DxcvXqyxuIqa3Kpai/fII48AAPbu3VvtfV+5cgXh4eElkjq9Xo/vv/++3Nc6OzujR48eWLp0KV5//XUUFBSUGotCoUCzZs0wY8YMHDx4EACwffv2asdclqqcj6LPflGza3GlvT+Aql8nFxcXPPTQQ7h16xb++uuvEuvl+hxQ3cLEjugBSUlJGDVqFGJjYxESEmIwg31Rzd2D9xNNSkpCVFRUueU+99xzsLGxwfTp03Ht2jU89dRTJTqEX7t2DVevXi3x2nv37iE/P7/UjuVVERQUhMceewzXr1/HsmXLDNb9/PPP2LBhA9zd3TF06FCj9gPcb74eNWoU7ty5g5dffrlEc1BWVpbU9FSbcVXFmDFjYGdnh+XLlxvcY1Sv1+OVV14ptYnrpZdeAgA8++yz0oS5xWVnZ0t9u6rL3d0dCoWi1EEK5Rk0aBDCwsKwc+dOxMTElFhfWt+3B4WFheGvv/4yODYhBObPn49Lly6V2P7YsWOlJuuJiYkA7k82DAAXL16UlpW3nSlNnToVtra2WLhwYamxFz8fZX329+/fX+pgBgBSU3tVrtPTTz8NIQReeeUVg4Tw7t27Uo3l008/XenyyPqwKZasWtGdJ/R6vXRLse+//x4FBQVo37491q9fbzByr127dujcuTO2bduGTp06oUuXLkhMTMTevXvRuHHjcu/0EBISggEDBmDnzp0AUGoz7G+//YZhw4ahXbt2aNq0KQICApCcnIwdO3agsLBQ6nNnjFWrVqFz58545ZVXcODAAbRt21aaL06pVGLNmjUmG4G4YsUKXLhwAatWrUJsbCz69OkDe3t7XLt2Dfv378fOnTuleeRqM67KCgsLwzvvvIN//etfePjhh6URzPv370daWhpatmyJc+fOGbymZ8+eeOeddzBnzhw0bNgQ/fv3R7169ZCVlYUbN27g6NGj6NKlC/bt21ftuNRqNTp06IDjx49jzJgxaNSokTQPW8uWLct8nb29PTZv3ozevXvjqaeewqeffopHHnkEeXl5uHz5Mg4dOlRqElbcSy+9hClTpuDhhx/G8OHDYWdnhx9++AGXLl3CoEGDsGvXLoPtZ86ciVu3bqFz584ICwuDvb09zpw5g8OHDyM0NBSjRo0CABw8eBCvvPIKOnbsiEaNGsHHxwc3b97Ejh07oFQq8corr1T7fJUlPDwcK1eulI5n8ODBaNiwIVJSUnDq1CloNBqplmzatGlYs2YNnnjiCYwYMQIBAQG4cOEC9u3bhyeffBKbNm0qUX7Pnj2xefNmDBs2DP3794ejoyNCQ0Mxbty4MmN6+eWXsXfvXuzYsQOtWrVC//79kZOTg82bNyMpKQmvvvoqunTpYvJzQRZE3tlWiOSBB2aDt7e3F56enqJNmzZi8uTJYu/evSXmXSuSkpIipk6dKkJDQ4VKpRL169cXc+bMEdnZ2WXOJF9k+/btAoBo27Ztqevj4+PFnDlzRKdOnYSvr6+wt7cXgYGBom/fvmLPnj2VPr7y5vMSQoibN2+KKVOmiJCQEGFnZyc8PT3F4MGDxcmTJ6tcVkWysrLEokWLRIsWLYSjo6NQq9WiadOm4oUXXhCJiYk1FhcemOVfiKrPY1dkw4YN4uGHHxYqlUp4eXmJMWPGVHjniePHj4snnnhC+Pv7S3dnaNWqlXjppZfEqVOnDLYt731TVlx//fWXGDhwoPDw8BAKhaJK1+jGjRti6tSpIiwsTNjZ2QkPDw/Rvn17sXjx4krFtWbNGtGqVSvh5OQkPD09xZAhQ8S5c+dKjXXTpk1i1KhRokGDBsLZ2Vm4uLiIZs2aiddff10kJSVJ2126dEm89NJLIiIiQnh5eQl7e3sRGhoqhg8fbjBvmxCmm8euyI8//iiGDRsmvL29hZ2dnfD39xd9+vQRmzdvNtjuhx9+EN27dxdubm5CrVaLzp07i2+++abMeLRarZgzZ46oV6+eNCdmZe48kZubKxYvXiyaNWsmHBwcpH1t2LChxLYVzWlY0Vx6ZHkUQvwznwMR1bj58+djwYIF+Pzzz/HMM8/IHQ4REVkYJnZEtSQzMxMNGzZEYWEh4uPja6TPEBERWTf2sSOqYbt378Yvv/yCXbt2ITExEe+//z6TOiIiqhFM7Ihq2ObNm7Fu3Tr4+vpizpw50qhJIiIiU2NTLBEREZGF4Dx2RERERBaCiR0RERGRhWBiR0RERGQhmNgRERERWQgmdkREREQWgokdERERkYVgYkdERERkIZjYEREREVkIJnZEREREFuL/AVelvXCbjQUoAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "bins = np.arange(0, D+1, 2)\n",
    "\n",
    "L1_median = int(np.median(L1_days))\n",
    "L2_median = int(np.median(L2_days))\n",
    "L3_median = int(np.median(L3_days))\n",
    "\n",
    "plt.hist(L1_days, color=\"C1\", alpha=0.3, bins=bins, density=True)\n",
    "plt.hist(L1_days, color=\"C1\", bins=bins, density=True, histtype='step', linewidth=2)\n",
    "\n",
    "plt.hist(L2_days, color=\"C2\", alpha=0.3, bins=bins, density=True)\n",
    "plt.hist(L2_days, color=\"C2\", bins=bins, density=True, histtype='step', linewidth=2)\n",
    "\n",
    "plt.hist(L3_days, color=\"C3\", alpha=0.3, bins=bins, density=True)\n",
    "plt.hist(L3_days, color=\"C3\", bins=bins, density=True, histtype='step', linewidth=2)\n",
    "\n",
    "\n",
    "plt.axvline(L1_median, color=\"C1\", linestyle='--', linewidth=2, label=f\"Median Time to Transient Classification = {L1_median} day\")\n",
    "plt.axvline(L2_median, color=\"C2\", linestyle='--', linewidth=2, label=f\"Median Time to Long Transient Classification = {L2_median} day\")\n",
    "plt.axvline(L3_median, color=\"C3\", linestyle='--', linewidth=2, label=f\"Median Time to Tidal Disruption Event Classification = {L3_median} day\")\n",
    "\n",
    "plt.legend(fontsize='small')\n",
    "\n",
    "plt.xlabel(\"Days for confident classification\", fontsize='x-large')\n",
    "plt.ylabel(\"Density\",  fontsize='x-large')\n",
    "\n",
    "plt.yscale('log')\n",
    "plt.tight_layout()\n",
    "#plt.savefig(f'{model_dir}/TDE_levels.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 525,
   "id": "70b88671-ccd6-45fb-95c8-3b00176be1e0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 525,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAicAAAGdCAYAAADJ6dNTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAYxElEQVR4nO3dfWyVd9348U+BtYyNljG0pQJjRtmskxKh1KqLbuskaMgejNkfi3aYYNTObKuasDtxuETD4haGuJMbHzL5S0FMxqKLRtIN0IWNrgy3WfeAwa0JtEgWWui2srTX7w+zan8w7hV6er6n5/VKTkLPdXqdD9826TvXOee6yrIsywIAIBFTCj0AAMB/EycAQFLECQCQFHECACRFnAAASREnAEBSxAkAkBRxAgAkZVqhBxir4eHhOHz4cMycOTPKysoKPQ4A8B5kWRYnTpyI2tramDLl7MdGii5ODh8+HPPnzy/0GADAOeju7o558+ad9TFFFyczZ86MiH//5yorKws8DQDwXvT398f8+fNH/o6fTdHFyTsv5VRWVooTACgy7+UtGd4QCwAkRZwAAEkRJwBAUsQJAJAUcQIAJEWcAABJEScAQFLECQCQFHECACRFnAAASREnAEBSxAkAkBRxAgAkpWBXJX7jjTfiIx/5SHzpS1+KBx54oFBjnObBnS8XeoQxu+v6RYUeAQDGTcGOnPzwhz+MT3ziE4V6egAgUQWJk1deeSVefPHFWLlyZSGeHgBI2JjjZM+ePbFq1aqora2NsrKy2LFjx2mPyeVysXDhwpg+fXo0NjbGvn37Rm3/zne+E+vXrz/noQGAyWvMcTIwMBD19fWRy+XOuH3btm3R1tYW69ati/3790d9fX2sWLEijh49GhERjz76aCxatCgWLXpv75MYHByM/v7+UTcAYPIa8xtiV65cedaXYzZs2BBr1qyJ1atXR0TE5s2b47HHHouHH3441q5dG0899VRs3bo1tm/fHidPnoy33347Kisr45577jnj/tavXx/33nvvWMcEAIrUuL7n5NSpU9HZ2RnNzc3/eYIpU6K5uTn27t0bEf+Oje7u7vjnP/8ZDzzwQKxZs+ZdwyQi4u67746+vr6RW3d393iODAAkZlw/Snzs2LEYGhqK6urqUfdXV1fHiy++eE77rKioiIqKivEYDwAoAgU7z0lExG233VbIpwcAEjSuL+vMmTMnpk6dGr29vaPu7+3tjZqamvF8KgBgkhrXOCkvL4+lS5dGe3v7yH3Dw8PR3t4eTU1N4/lUAMAkNeaXdU6ePBkHDx4c+frQoUNx4MCBmD17dixYsCDa2tqipaUlli1bFsuXL4+NGzfGwMDAyKd3zlUul4tcLhdDQ0PntR8AIG1lWZZlY/mGXbt2xTXXXHPa/S0tLbFly5aIiHjooYfi/vvvj56enliyZEls2rQpGhsbx2Xg/v7+qKqqir6+vqisrByXff4319YBgPE3lr/fY46TQhMnpxMnAKRuLH+/C3bhPwCAMxEnAEBSxAkAkJSiiZNcLhd1dXXR0NBQ6FEAgDwqmjhpbW2Nrq6u6OjoKPQoAEAeFU2cAAClQZwAAEkRJwBAUsQJAJAUcQIAJKVo4sRHiQGgNBRNnPgoMQCUhqKJEwCgNIgTACAp4gQASIo4AQCSIk4AgKSIEwAgKeIEAEhK0cSJk7ABQGkomjhxEjYAKA1FEycAQGkQJwBAUsQJAJAUcQIAJEWcAABJEScAQFLECQCQlKKJEydhA4DSUDRx4iRsAFAaiiZOAIDSIE4AgKSIEwAgKeIEAEiKOAEAkiJOAICkiBMAICniBABIijgBAJIiTgCApBRNnLi2DgCUhqKJE9fWAYDSUDRxAgCUBnECACRFnAAASREnAEBSxAkAkBRxAgAkRZwAAEkRJwBAUsQJAJAUcQIAJEWcAABJEScAQFLECQCQFHECACSlaOIkl8tFXV1dNDQ0FHoUACCPiiZOWltbo6urKzo6Ogo9CgCQR0UTJwBAaRAnAEBSxAkAkBRxAgAkRZwAAEkRJwBAUsQJAJAUcQIAJEWcAABJEScAQFLECQCQFHECACRFnAAASREnAEBSxAkAkBRxAgAkRZwAAEkRJwBAUsQJAJCUoomTXC4XdXV10dDQUOhRAIA8Kpo4aW1tja6urujo6Cj0KABAHhVNnAAApUGcAABJEScAQFLECQCQFHECACRFnAAASREnAEBSxAkAkBRxAgAkRZwAAEkRJwBAUsQJAJAUcQIAJEWcAABJEScAQFLECQCQFHECACRFnAAASREnAEBSxAkAkBRxAgAkRZwAAEkRJwBAUsQJAJAUcQIAJEWcAABJEScAQFLECQCQFHECACRlwuPk+PHjsWzZsliyZElcddVV8fOf/3yiRwAAEjZtop9w5syZsWfPnpgxY0YMDAzEVVddFTfffHNceumlEz0KAJCgCT9yMnXq1JgxY0ZERAwODkaWZZFl2USPAQAkasxxsmfPnli1alXU1tZGWVlZ7Nix47TH5HK5WLhwYUyfPj0aGxtj3759o7YfP3486uvrY968efHd73435syZc87/AQBgchlznAwMDER9fX3kcrkzbt+2bVu0tbXFunXrYv/+/VFfXx8rVqyIo0ePjjxm1qxZ8de//jUOHToUv/rVr6K3t/fc/wcAwKQy5jhZuXJl/OAHP4ibbrrpjNs3bNgQa9asidWrV0ddXV1s3rw5ZsyYEQ8//PBpj62uro76+vr485///K7PNzg4GP39/aNuAMDkNa7vOTl16lR0dnZGc3Pzf55gypRobm6OvXv3RkREb29vnDhxIiIi+vr6Ys+ePXHFFVe86z7Xr18fVVVVI7f58+eP58gAQGLGNU6OHTsWQ0NDUV1dPer+6urq6OnpiYiIV199Na6++uqor6+Pq6++Or71rW/Fxz72sXfd59133x19fX0jt+7u7vEcGQBIzIR/lHj58uVx4MCB9/z4ioqKqKioyN9AAEBSxvXIyZw5c2Lq1KmnvcG1t7c3ampqxvOpAIBJalzjpLy8PJYuXRrt7e0j9w0PD0d7e3s0NTWN51MBAJPUmF/WOXnyZBw8eHDk60OHDsWBAwdi9uzZsWDBgmhra4uWlpZYtmxZLF++PDZu3BgDAwOxevXqcR0cAJicxhwnzzzzTFxzzTUjX7e1tUVEREtLS2zZsiVuueWW+Ne//hX33HNP9PT0xJIlS+KPf/zjaW+SHatcLhe5XC6GhobOaz8AQNrKsiI7d3x/f39UVVVFX19fVFZWjvv+H9z58rjvM9/uun5RoUcAgLMay9/vCb+2DgDA2YgTACAp4gQASIo4AQCSUjRxksvloq6uLhoaGgo9CgCQR0UTJ62trdHV1RUdHR2FHgUAyKOiiRMAoDSIEwAgKeIEAEiKOAEAkiJOAICkFE2c+CgxAJSGookTHyUGgNJQNHECAJQGcQIAJEWcAABJEScAQFLECQCQFHECACSlaOLEeU4AoDQUTZw4zwkAlIaiiRMAoDSIEwAgKeIEAEiKOAEAkiJOAICkiBMAICniBABIijgBAJJSNHHiDLEAUBqKJk6cIRYASkPRxAkAUBrECQCQFHECACRFnAAASREnAEBSxAkAkBRxAgAkRZwAAEkRJwBAUsQJAJCUookT19YBgNJQNHHi2joAUBqKJk4AgNIgTgCApIgTACAp4gQASIo4AQCSIk4AgKSIEwAgKeIEAEiKOAEAkiJOAICkiBMAICniBABIijgBAJJSNHGSy+Wirq4uGhoaCj0KAJBHRRMnra2t0dXVFR0dHYUeBQDIo6KJEwCgNIgTACAp4gQASIo4AQCSIk4AgKSIEwAgKeIEAEiKOAEAkiJOAICkiBMAICniBABIijgBAJIiTgCApEwr9ACcvwd3vlzoEcbsrusXFXoEABLlyAkAkBRxAgAkRZwAAEkRJwBAUsQJAJAUcQIAJKVo4iSXy0VdXV00NDQUehQAII+KJk5aW1ujq6srOjo6Cj0KAJBHRRMnAEBpECcAQFLECQCQFHECACRFnAAASREnAEBSxAkAkBRxAgAkRZwAAEkRJwBAUsQJAJAUcQIAJEWcAABJEScAQFLECQCQFHECACRFnAAASREnAEBSphV6AErTgztfLvQIY3bX9YsKPQJASXDkBABIijgBAJIiTgCApIgTACAp4gQASIo4AQCSIk4AgKSIEwAgKeIEAEiKOAEAkjLhcdLd3R2f/exno66uLhYvXhzbt2+f6BEAgIRN+LV1pk2bFhs3bowlS5ZET09PLF26ND7/+c/HRRddNNGjAAAJmvA4mTt3bsydOzciImpqamLOnDnx+uuvixMAICLO4WWdPXv2xKpVq6K2tjbKyspix44dpz0ml8vFwoULY/r06dHY2Bj79u074746OztjaGgo5s+fP+bBAYDJacxxMjAwEPX19ZHL5c64fdu2bdHW1hbr1q2L/fv3R319faxYsSKOHj066nGvv/56fOUrX4mf/exn5zY5ADApjfllnZUrV8bKlSvfdfuGDRtizZo1sXr16oiI2Lx5czz22GPx8MMPx9q1ayMiYnBwMG688cZYu3ZtfPKTnzzr8w0ODsbg4ODI1/39/WMdGQAoIuP6aZ1Tp05FZ2dnNDc3/+cJpkyJ5ubm2Lt3b0REZFkWt912W1x77bXx5S9/+f/c5/r166Oqqmrk5iUgAJjcxjVOjh07FkNDQ1FdXT3q/urq6ujp6YmIiCeffDK2bdsWO3bsiCVLlsSSJUvi+eeff9d93n333dHX1zdy6+7uHs+RAYDETPindT796U/H8PDwe358RUVFVFRU5HEiACAl43rkZM6cOTF16tTo7e0ddX9vb2/U1NSM51MBAJPUuMZJeXl5LF26NNrb20fuGx4ejvb29mhqahrPpwIAJqkxv6xz8uTJOHjw4MjXhw4digMHDsTs2bNjwYIF0dbWFi0tLbFs2bJYvnx5bNy4MQYGBkY+vXOucrlc5HK5GBoaOq/9wLl6cOfLhR5hzO66flGhRwAYs7Isy7KxfMOuXbvimmuuOe3+lpaW2LJlS0REPPTQQ3H//fdHT09PLFmyJDZt2hSNjY3jMnB/f39UVVVFX19fVFZWjss+/1sx/gGCdyNOgFSM5e/3mOOk0MQJvHfiBEjFWP5+T/hViQEAzkacAABJEScAQFLECQCQlKKJk1wuF3V1ddHQ0FDoUQCAPCqaOGltbY2urq7o6Ogo9CgAQB4VTZwAAKVBnAAASREnAEBSxnxtHYB8KsazNDsTL4wvR04AgKQUTZz4KDEAlIaiiRMfJQaA0lA0cQIAlAZxAgAkRZwAAEnxUWKYxIrxY7kAjpwAAEkRJwBAUoomTpznBABKQ9HEifOcAEBpKJo4AQBKgzgBAJIiTgCApIgTACAp4gQASIo4AQCSIk4AgKS4tg7AeSrGaxjddf2iQo8A76pojpw4QywAlIaiiRNniAWA0lA0cQIAlAZxAgAkRZwAAEkRJwBAUsQJAJAUcQIAJEWcAABJEScAQFLECQCQFHECACSlaOLEtXUAoDQUTZy4tg4AlIaiiRMAoDSIEwAgKeIEAEiKOAEAkiJOAICkiBMAICniBABIijgBAJIiTgCApEwr9AAATLwHd75c6BHG7K7rFxV6BCaIIycAQFLECQCQFHECACSlaOIkl8tFXV1dNDQ0FHoUACCPiiZOWltbo6urKzo6Ogo9CgCQR0UTJwBAaRAnAEBSxAkAkBRxAgAkRZwAAEkRJwBAUlxbB4Ci4HpApcOREwAgKeIEAEiKOAEAkiJOAICkiBMAICniBABIijgBAJIiTgCApIgTACAp4gQASIo4AQCSIk4AgKQUzYX/crlc5HK5GBoaKvQoAPCeFOPFCiMKf8HCojly0traGl1dXdHR0VHoUQCAPCqaOAEASoM4AQCSIk4AgKSIEwAgKeIEAEiKOAEAkiJOAICkiBMAICniBABIijgBAJIiTgCApIgTACAp4gQASIo4AQCSMq3QA4xVlmUREdHf35+X/b81cDIv+wWAYpGPv7Hv7POdv+NnU3RxcuLEiYiImD9/foEnAYDJ6X/yuO8TJ05EVVXVWR9Tlr2XhEnI8PBwHD58OGbOnBllZWXv6Xv6+/tj/vz50d3dHZWVlXmekHdY98Kw7oVh3QvDuhfGuax7lmVx4sSJqK2tjSlTzv6ukqI7cjJlypSYN2/eOX1vZWWlX94CsO6FYd0Lw7oXhnUvjLGu+/91xOQd3hALACRFnAAASSmJOKmoqIh169ZFRUVFoUcpKda9MKx7YVj3wrDuhZHvdS+6N8QCAJNbSRw5AQCKhzgBAJIiTgCApIgTACApJREnuVwuFi5cGNOnT4/GxsbYt29foUeaVPbs2ROrVq2K2traKCsrix07dozanmVZ3HPPPTF37ty48MILo7m5OV555ZXCDDtJrF+/PhoaGmLmzJnx/ve/P2688cZ46aWXRj3mrbfeitbW1rj00kvj4osvji9+8YvR29tboIknh//93/+NxYsXj5x4qqmpKf7whz+MbLfmE+O+++6LsrKyuPPOO0fus/bj7/vf/36UlZWNul155ZUj2/O55pM+TrZt2xZtbW2xbt262L9/f9TX18eKFSvi6NGjhR5t0hgYGIj6+vrI5XJn3P6jH/0oNm3aFJs3b46nn346LrroolixYkW89dZbEzzp5LF79+5obW2Np556Knbu3Blvv/12fO5zn4uBgYGRx9x1113xu9/9LrZv3x67d++Ow4cPx80331zAqYvfvHnz4r777ovOzs545pln4tprr40bbrgh/va3v0WENZ8IHR0d8dOf/jQWL1486n5rnx8f/ehH48iRIyO3v/zlLyPb8rrm2SS3fPnyrLW1deTroaGhrLa2Nlu/fn0Bp5q8IiJ75JFHRr4eHh7Oampqsvvvv3/kvuPHj2cVFRXZr3/96wJMODkdPXo0i4hs9+7dWZb9e40vuOCCbPv27SOP+fvf/55FRLZ3795CjTkpXXLJJdkvfvELaz4BTpw4kX34wx/Odu7cmX3mM5/J7rjjjizL/L7ny7p167L6+vozbsv3mk/qIyenTp2Kzs7OaG5uHrlvypQp0dzcHHv37i3gZKXj0KFD0dPTM+pnUFVVFY2NjX4G46ivry8iImbPnh0REZ2dnfH222+PWvcrr7wyFixYYN3HydDQUGzdujUGBgaiqanJmk+A1tbW+MIXvjBqjSP8vufTK6+8ErW1tfHBD34wbr311njttdciIv9rXnQX/huLY8eOxdDQUFRXV4+6v7q6Ol588cUCTVVaenp6IiLO+DN4ZxvnZ3h4OO6888741Kc+FVdddVVE/Hvdy8vLY9asWaMea93P3/PPPx9NTU3x1ltvxcUXXxyPPPJI1NXVxYEDB6x5Hm3dujX2798fHR0dp23z+54fjY2NsWXLlrjiiiviyJEjce+998bVV18dL7zwQt7XfFLHCZSC1tbWeOGFF0a9Fkz+XHHFFXHgwIHo6+uL3/72t9HS0hK7d+8u9FiTWnd3d9xxxx2xc+fOmD59eqHHKRkrV64c+ffixYujsbExLrvssvjNb34TF154YV6fe1K/rDNnzpyYOnXqae8e7u3tjZqamgJNVVreWWc/g/y4/fbb4/e//3088cQTMW/evJH7a2pq4tSpU3H8+PFRj7fu56+8vDw+9KEPxdKlS2P9+vVRX18fP/7xj615HnV2dsbRo0fj4x//eEybNi2mTZsWu3fvjk2bNsW0adOiurra2k+AWbNmxaJFi+LgwYN5/32f1HFSXl4eS5cujfb29pH7hoeHo729PZqamgo4Wem4/PLLo6amZtTPoL+/P55++mk/g/OQZVncfvvt8cgjj8Tjjz8el19++ajtS5cujQsuuGDUur/00kvx2muvWfdxNjw8HIODg9Y8j6677rp4/vnn48CBAyO3ZcuWxa233jryb2uffydPnox//OMfMXfu3Pz/vp/3W2oTt3Xr1qyioiLbsmVL1tXVlX3ta1/LZs2alfX09BR6tEnjxIkT2bPPPps9++yzWURkGzZsyJ599tns1VdfzbIsy+67775s1qxZ2aOPPpo999xz2Q033JBdfvnl2ZtvvlngyYvXN77xjayqqirbtWtXduTIkZHbG2+8MfKYr3/969mCBQuyxx9/PHvmmWeypqamrKmpqYBTF7+1a9dmu3fvzg4dOpQ999xz2dq1a7OysrLsT3/6U5Zl1nwi/fendbLM2ufDt7/97WzXrl3ZoUOHsieffDJrbm7O5syZkx09ejTLsvyu+aSPkyzLsp/85CfZggULsvLy8mz58uXZU089VeiRJpUnnngii4jTbi0tLVmW/fvjxN/73vey6urqrKKiIrvuuuuyl156qbBDF7kzrXdEZL/85S9HHvPmm29m3/zmN7NLLrkkmzFjRnbTTTdlR44cKdzQk8BXv/rV7LLLLsvKy8uz973vfdl11103EiZZZs0n0v8fJ9Z+/N1yyy3Z3Llzs/Ly8uwDH/hAdsstt2QHDx4c2Z7PNS/Lsiw7/+MvAADjY1K/5wQAKD7iBABIijgBAJIiTgCApIgTACAp4gQASIo4AQCSIk4AgKSIEwAgKeIEAEiKOAEAkiJOAICk/D9+AKiUBPpgwwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(L1_days, alpha=0.5, bins=10)\n",
    "plt.yscale('log')\n",
    "np.median(L1_days)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 526,
   "id": "31a81983-93c8-4c4b-b72c-0c60123d93cf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "11.0"
      ]
     },
     "execution_count": 526,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiYAAAGdCAYAAAAmK7htAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAATcUlEQVR4nO3db2yVd/n48avAWoYrdBtaqMCGUbZUpItQSGOIUZhkMcvcfMCDPSAz0TjPzAbGZPslDpeYQDTZdObEGY3ySJkz2YwmGgnbihq2lX9us8qGIRsJ/9wDCqsDFvr5PTBrvnw3/I6tp/d1Dq9XchJ67sM5F1eb8M7duz1tpZQSAAAJTKl6AACAtwgTACANYQIApCFMAIA0hAkAkIYwAQDSECYAQBrCBABIY1rVA1yssbGxOHz4cHR2dkZbW1vV4wAA70IpJU6dOhU9PT0xZcqFz4s0XZgcPnw45s+fX/UYAMB7cOjQoZg3b94FjzddmHR2dkbEf/5hM2fOrHgaAODdOHnyZMyfP3/8//ELaboweevbNzNnzhQmANBk/q/LMFz8CgCkIUwAgDSECQCQhjABANIQJgBAGsIEAEhDmAAAaQgTACANYQIApCFMAIA0hAkAkIYwAQDSECYAQBpN9+7CjfTQtpeqHuE9WX/joqpHAIAJ4YwJAJCGMAEA0hAmAEAawgQASEOYAABpCBMAIA1hAgCkIUwAgDSECQCQhjABANIQJgBAGsIEAEhDmAAAaTRNmNTr9ejt7Y3+/v6qRwEAGqRpwqRWq8Xw8HAMDQ1VPQoA0CBNEyYAQOsTJgBAGsIEAEhDmAAAaQgTACANYQIApCFMAIA0hAkAkIYwAQDSECYAQBrCBABIQ5gAAGkIEwAgDWECAKQhTACANIQJAJCGMAEA0hAmAEAawgQASEOYAABpCBMAIA1hAgCkIUwAgDSECQCQhjABANIQJgBAGsIEAEhDmAAAaQgTACANYQIApCFMAIA0hAkAkIYwAQDSECYAQBrCBABIQ5gAAGkIEwAgDWECAKQhTACANIQJAJCGMAEA0hAmAEAawgQASEOYAABpCBMAIA1hAgCkIUwAgDSECQCQRtOESb1ej97e3ujv7696FACgQZomTGq1WgwPD8fQ0FDVowAADdI0YQIAtD5hAgCkIUwAgDSECQCQhjABANIQJgBAGsIEAEhDmAAAaQgTACANYQIApCFMAIA0hAkAkIYwAQDSECYAQBrCBABIQ5gAAGkIEwAgDWECAKQhTACANIQJAJCGMAEA0hAmAEAawgQASEOYAABpCBMAIA1hAgCkIUwAgDSECQCQhjABANIQJgBAGsIEAEhDmAAAaQgTACANYQIApCFMAIA0hAkAkIYwAQDSECYAQBrCBABIQ5gAAGkIEwAgDWECAKQhTACANIQJAJCGMAEA0hAmAEAawgQASEOYAABpCBMAIA1hAgCkIUwAgDSECQCQhjABANIQJgBAGsIEAEhDmAAAaQgTACANYQIApCFMAIA0hAkAkIYwAQDSECYAQBrCBABIQ5gAAGkIEwAgDWECAKQhTACANIQJAJCGMAEA0hAmAEAawgQASEOYAABpCBMAIA1hAgCkIUwAgDSECQCQhjABANIQJgBAGsIEAEhj0sPkxIkTsWzZsrjhhhti8eLF8ZOf/GSyRwAAkpo22S/Y2dkZO3bsiBkzZsTo6GgsXrw4brvttrj66qsnexQAIJlJP2MyderUmDFjRkREnDlzJkopUUqZ7DEAgIQuOkx27NgRN998c/T09ERbW1s88cQTb3tMvV6Pa6+9NqZPnx4rVqyI55577rzjJ06ciL6+vpg3b15885vfjNmzZ7/nfwAA0DouOkxGR0ejr68v6vX6Ox5/9NFHY8OGDbFx48bYs2dP9PX1xZo1a+L48ePjj+nq6oq//vWvcfDgwfjFL34Rx44de+//AgCgZVx0mNx0003xne98J2699dZ3PP7ggw/Gl7/85bjjjjuit7c3HnnkkZgxY0b87Gc/e9tju7u7o6+vL/70pz9d8PXOnDkTJ0+ePO8GALSmCb349ezZs7F79+647777xu+bMmVKrF69Onbu3BkREceOHYsZM2ZEZ2dnjIyMxI4dO+LOO++84HNu2rQpHnjggYkcs+U8tO2lqke4aOtvXFT1CAAkNKEXv7722mtx7ty56O7uPu/+7u7uOHr0aEREvPLKK7Fy5cro6+uLlStXxte//vX4xCc+ccHnvO+++2JkZGT8dujQoYkcGQBIZNJ/XHj58uWxb9++d/34jo6O6OjoaNxAAEAaE3rGZPbs2TF16tS3Xcx67NixmDNnzkS+FADQgiY0TNrb22Pp0qWxffv28fvGxsZi+/btMTAwMJEvBQC0oIv+Vs7rr78eBw4cGP/44MGDsW/fvrjqqqtiwYIFsWHDhli3bl0sW7Ysli9fHt///vdjdHQ07rjjjgkdHABoPRcdJrt27YrPfOYz4x9v2LAhIiLWrVsXW7ZsibVr18a//vWvuP/+++Po0aNxww03xB/+8Ie3XRALAPC/tZUm+33wJ0+ejFmzZsXIyEjMnDlzQp+7GX/stln5cWGAS8u7/f970t8rBwDgQoQJAJCGMAEA0hAmAEAak/6bXwH+m2a8CN3F3DBxmuaMSb1ej97e3ujv7696FACgQZomTGq1WgwPD8fQ0FDVowAADdI0YQIAtD5hAgCkIUwAgDSECQCQhjABANIQJgBAGsIEAEhDmAAAaQgTACANYQIApCFMAIA0hAkAkEbThIl3FwaA1tc0YeLdhQGg9TVNmAAArU+YAABpCBMAIA1hAgCkIUwAgDSECQCQhjABANKYVvUAAM3uoW0vVT3CRVt/46KqR4B35IwJAJCGMAEA0hAmAEAarjGhEr4nD8A7ccYEAEhDmAAAaQgTACCNpgmTer0evb290d/fX/UoAECDNE2Y1Gq1GB4ejqGhoapHAQAapGnCBABofcIEAEhDmAAAaQgTACANYQIApCFMAIA0vFcOtLBmfE8i4NLmjAkAkIYwAQDSECYAQBrCBABIQ5gAAGkIEwAgDWECAKQhTACANIQJAJCGMAEA0miaMKnX69Hb2xv9/f1VjwIANEjThEmtVovh4eEYGhqqehQAoEGaJkwAgNbn3YXhXfJOvQCN54wJAJCGMAEA0hAmAEAawgQASEOYAABpCBMAIA1hAgCkIUwAgDSECQCQhjABANIQJgBAGsIEAEhDmAAAaQgTACANYQIApCFMAIA0hAkAkIYwAQDSECYAQBpNEyb1ej16e3ujv7+/6lEAgAZpmjCp1WoxPDwcQ0NDVY8CADTItKoHAGDyPbTtpapHuGjrb1xU9QhMgqY5YwIAtD5hAgCkIUwAgDSECQCQhjABANIQJgBAGsIEAEhDmAAAaQgTACANYQIApCFMAIA0hAkAkIYwAQDSECYAQBrCBABIQ5gAAGkIEwAgDWECAKQhTACANIQJAJCGMAEA0hAmAEAawgQASEOYAABpNE2Y1Ov16O3tjf7+/qpHAQAapGnCpFarxfDwcAwNDVU9CgDQIE0TJgBA6xMmAEAawgQASEOYAABpCBMAIA1hAgCkIUwAgDSECQCQhjABANIQJgBAGsIEAEhDmAAAaQgTACANYQIApCFMAIA0hAkAkIYwAQDSECYAQBrCBABIQ5gAAGkIEwAgDWECAKQhTACANIQJAJCGMAEA0hAmAEAawgQASEOYAABpCBMAIA1hAgCkIUwAgDSECQCQhjABANIQJgBAGk0TJvV6PXp7e6O/v7/qUQCABmmaMKnVajE8PBxDQ0NVjwIANEjThAkA0PqECQCQhjABANKYVvUAAPBuPLTtpapHuGjrb1xU9QhNxxkTACANYQIApCFMAIA0hAkAkIYwAQDSECYAQBrCBABIQ5gAAGkIEwAgDWECAKQhTACANIQJAJCGMAEA0hAmAEAawgQASEOYAABpCBMAIA1hAgCkIUwAgDSECQCQhjABANIQJgBAGsIEAEhDmAAAaQgTACANYQIApCFMAIA0hAkAkIYwAQDSECYAQBrTqh4AAFrVQ9teqnqEi7b+xkWVvr4zJgBAGsIEAEhDmAAAaQgTACANYQIApCFMAIA0hAkAkIYwAQDSECYAQBrCBABIQ5gAAGkIEwAgDWECAKQhTACANJomTOr1evT29kZ/f3/VowAADdJWSilVD3ExRkZGoqurKw4dOhQzZ86c0OeuP3lgQp8PAJpN7bMfbcjznjx5MubPnx8nTpyIWbNmXfBx0xry6g106tSpiIiYP39+xZMAQOv5fw1+/lOnTv3XMGm6MyZjY2Nx+PDh6OzsjLa2tnf1d96qtEacZeHC7L0a9l4Ne6+GvVfjvey9lBKnTp2Knp6emDLlwleSNN0ZkylTpsS8efPe09+dOXOmL9wK2Hs17L0a9l4Ne6/Gxe79v50peUvTXPwKALQ+YQIApHFJhElHR0ds3LgxOjo6qh7lkmLv1bD3ath7Ney9Go3ce9Nd/AoAtK5L4owJANAchAkAkIYwAQDSECYAQBqXRJjU6/W49tprY/r06bFixYp47rnnqh6ppezYsSNuvvnm6Onpiba2tnjiiSfOO15Kifvvvz/mzp0bl19+eaxevTpefvnlaoZtEZs2bYr+/v7o7OyMD33oQ/GFL3wh9u/ff95jTp8+HbVaLa6++uq44oor4otf/GIcO3asoolbw49+9KNYsmTJ+C+VGhgYiN///vfjx+18cmzevDna2trinnvuGb/P7ifet7/97Whrazvvdv31148fb9TOWz5MHn300diwYUNs3Lgx9uzZE319fbFmzZo4fvx41aO1jNHR0ejr64t6vf6Ox7/73e/Gww8/HI888kg8++yz8YEPfCDWrFkTp0+fnuRJW8fg4GDUarV45plnYtu2bfHmm2/G5z73uRgdHR1/zPr16+O3v/1tPPbYYzE4OBiHDx+O2267rcKpm9+8efNi8+bNsXv37ti1a1d89rOfjVtuuSX+9re/RYSdT4ahoaH48Y9/HEuWLDnvfrtvjI9//ONx5MiR8duf//zn8WMN23lpccuXLy+1Wm3843PnzpWenp6yadOmCqdqXRFRHn/88fGPx8bGypw5c8r3vve98ftOnDhROjo6yi9/+csKJmxNx48fLxFRBgcHSyn/2fFll11WHnvssfHH/P3vfy8RUXbu3FnVmC3pyiuvLD/96U/tfBKcOnWqfOxjHyvbtm0rn/70p8vdd99dSvH13igbN24sfX1973iskTtv6TMmZ8+ejd27d8fq1avH75syZUqsXr06du7cWeFkl46DBw/G0aNHz/sczJo1K1asWOFzMIFGRkYiIuKqq66KiIjdu3fHm2++ed7er7/++liwYIG9T5Bz587F1q1bY3R0NAYGBux8EtRqtfj85z9/3o4jfL030ssvvxw9PT3xkY98JG6//fZ49dVXI6KxO2+6N/G7GK+99lqcO3cuuru7z7u/u7s7/vGPf1Q01aXl6NGjERHv+Dl46xjvz9jYWNxzzz3xqU99KhYvXhwR/9l7e3t7dHV1nfdYe3//XnjhhRgYGIjTp0/HFVdcEY8//nj09vbGvn377LyBtm7dGnv27ImhoaG3HfP13hgrVqyILVu2xHXXXRdHjhyJBx54IFauXBkvvvhiQ3fe0mECl4JarRYvvvjied/7pXGuu+662LdvX4yMjMSvf/3rWLduXQwODlY9Vks7dOhQ3H333bFt27aYPn161eNcMm666abxPy9ZsiRWrFgR11xzTfzqV7+Kyy+/vGGv29Lfypk9e3ZMnTr1bVcJHzt2LObMmVPRVJeWt/bsc9AYd911V/zud7+Lp556KubNmzd+/5w5c+Ls2bNx4sSJ8x5v7+9fe3t7fPSjH42lS5fGpk2boq+vL37wgx/YeQPt3r07jh8/Hp/85Cdj2rRpMW3atBgcHIyHH344pk2bFt3d3XY/Cbq6umLRokVx4MCBhn69t3SYtLe3x9KlS2P79u3j942NjcX27dtjYGCgwskuHQsXLow5c+ac9zk4efJkPPvssz4H70MpJe666654/PHH48knn4yFCxeed3zp0qVx2WWXnbf3/fv3x6uvvmrvE2xsbCzOnDlj5w20atWqeOGFF2Lfvn3jt2XLlsXtt98+/me7b7zXX389/vnPf8bcuXMb+/X+vi6dbQJbt24tHR0dZcuWLWV4eLh85StfKV1dXeXo0aNVj9YyTp06Vfbu3Vv27t1bIqI8+OCDZe/eveWVV14ppZSyefPm0tXVVX7zm9+U559/vtxyyy1l4cKF5Y033qh48uZ15513llmzZpWnn366HDlyZPz273//e/wxX/3qV8uCBQvKk08+WXbt2lUGBgbKwMBAhVM3v3vvvbcMDg6WgwcPlueff77ce++9pa2trfzxj38spdj5ZPqfP5VTit03wje+8Y3y9NNPl4MHD5a//OUvZfXq1WX27Nnl+PHjpZTG7bzlw6SUUn74wx+WBQsWlPb29rJ8+fLyzDPPVD1SS3nqqadKRLzttm7dulLKf35k+Fvf+lbp7u4uHR0dZdWqVWX//v3VDt3k3mnfEVF+/vOfjz/mjTfeKF/72tfKlVdeWWbMmFFuvfXWcuTIkeqGbgFf+tKXyjXXXFPa29vLBz/4wbJq1arxKCnFzifT/w4Tu594a9euLXPnzi3t7e3lwx/+cFm7dm05cODA+PFG7bytlFLe3zkXAICJ0dLXmAAAzUWYAABpCBMAIA1hAgCkIUwAgDSECQCQhjABANIQJgBAGsIEAEhDmAAAaQgTACANYQIApPH/ARztXYpH0ZYbAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(L2_days, alpha=0.5, bins=10)\n",
    "plt.yscale('log')\n",
    "np.median(L2_days)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 527,
   "id": "1248f284-dcf4-4bde-8325-9545dc5cbbe9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "21.0"
      ]
     },
     "execution_count": 527,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj4AAAGdCAYAAAASUnlxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAbg0lEQVR4nO3de2zV9f348VehFkS5qGhL5aJm4uyUdnJp0BG1lBE0RHTZTHSzskQzV53SOYO/ZLIlRsyWKXM2Ot2FxWSjY1F2yW5dnTANSil04rqBODLJsEXi5FIVTM/n94exsV/qpYV6aN+PR3ISz+d8zue8+m5Nn/mcz6EFWZZlAQCQgGH5HgAA4OMifACAZAgfACAZwgcASIbwAQCSIXwAgGQIHwAgGcIHAEhGYb4HOJbkcrnYtWtXjB49OgoKCvI9DgDwEWRZFvv374/S0tIYNuyDz+kIn/fYtWtXTJo0Kd9jAAD9sHPnzpg4ceIH7iN83mP06NER8c7CjRkzJs/TAAAfxb59+2LSpEndv8c/iPB5j3ff3hozZozwAYBB5qNcpuLiZgAgGcIHAEiG8AEAkiF8AIBkCB8AIBnCBwBIhvABAJIhfACAZAgfACAZwgcASIbwAQCSIXwAgGQIHwAgGf46+8fo/sZt+R6hz5bMm5rvEQDgqHHGBwBIhvABAJIhfACAZAgfACAZwgcASIbwAQCSIXwAgGQIHwAgGcIHAEiG8AEAkiF8AIBkCB8AIBnCBwBIhvABAJIhfACAZAgfACAZwgcASIbwAQCSIXwAgGQIHwAgGcIHAEiG8AEAkiF8AIBkCB8AIBnCBwBIhvABAJIhfACAZAgfACAZwgcASIbwAQCSIXwAgGQIHwAgGcIHAEiG8AEAkiF8AIBkCB8AIBnCBwBIhvABAJIhfACAZAgfACAZwgcASIbwAQCSIXwAgGQIHwAgGcIHAEiG8AEAkiF8AIBkCB8AIBnCBwBIhvABAJIhfACAZAgfACAZwgcASIbwAQCSIXwAgGQIHwAgGcIHAEiG8AEAkiF8AIBkCB8AIBnCBwBIhvABAJIhfACAZAgfACAZwgcASIbwAQCSIXwAgGQIHwAgGcIHAEiG8AEAkiF8AIBkCB8AIBnCBwBIhvABAJIhfACAZAgfACAZwgcASIbwAQCSIXwAgGQMqfB5/fXXY8aMGVFRURHnnXdePProo/keCQA4hhTme4CjafTo0bFu3boYNWpUdHZ2xnnnnRdXXXVVnHLKKfkeDQA4BgypMz7Dhw+PUaNGRUTEwYMHI8uyyLIsz1MBAMeKYyp81q1bFwsXLozS0tIoKCiINWvWHLZPfX19nHHGGTFy5MiorKyMDRs29Hj89ddfj/Ly8pg4cWJ84xvfiPHjx39M0wMAx7pjKnw6OzujvLw86uvre328oaEh6urqYtmyZbFp06YoLy+P+fPnx+7du7v3GTduXPz973+PHTt2xM9//vPo6Oj4uMYHAI5xx1T4LFiwIO6+++648sore338vvvuixtuuCEWL14cZWVl8fDDD8eoUaPiJz/5yWH7FhcXR3l5efztb39739c7ePBg7Nu3r8cNABi6jqnw+SCHDh2KlpaWqK6u7t42bNiwqK6ujvXr10dEREdHR+zfvz8iIvbu3Rvr1q2Lc845532PuXz58hg7dmz3bdKkSQP7RQAAeTVowmfPnj3R1dUVxcXFPbYXFxdHe3t7RET85z//iTlz5kR5eXnMmTMnbrnlljj//PPf95h33nln7N27t/u2c+fOAf0aAID8GlIfZ581a1a0trZ+5P1HjBgRI0aMGLiBAIBjyqA54zN+/PgYPnz4YRcrd3R0RElJSZ6mAgAGk0FzxqeoqCimT58eTU1NsWjRooiIyOVy0dTUFDfffHN+hwMGhfsbt+V7hD5bMm9qvkeAIeWYCp8DBw7E9u3bu+/v2LEjWltb4+STT47JkydHXV1d1NTUxIwZM2LWrFmxYsWK6OzsjMWLF+dxagBgsDimwmfjxo1x6aWXdt+vq6uLiIiamppYuXJlXH311fHqq6/GXXfdFe3t7VFRURF//OMfD7vgGQCgN8dU+FxyySUf+icmbr75Zm9tAQD9MmgubgYAOFLCBwBIhvABAJIhfACAZAgfACAZwgcASIbwAQCSIXwAgGQIHwAgGcIHAEiG8AEAknFM/a0uAHq6v3FbvkfolyXzpuZ7BOiVMz4AQDKEDwCQDOEDACRD+AAAyRA+EVFfXx9lZWUxc+bMfI8CAAwg4RMRtbW10dbWFs3NzfkeBQAYQMIHAEiG8AEAkiF8AIBkCB8AIBnCBwBIhvABAJIhfACAZAgfACAZwgcASEZhvgfg2HZ/47Z8j9BnS+ZNzfcIAByjnPEBAJIhfACAZAgfACAZwgcASIbwAQCSIXwAgGQIHwAgGcIHAEiG8AEAkiF8AIBkCB8AIBnCBwBIhvABAJIhfACAZAgfACAZwgcASIbwiYj6+vooKyuLmTNn5nsUAGAACZ+IqK2tjba2tmhubs73KADAABI+AEAyhA8AkAzhAwAkQ/gAAMkQPgBAMoQPAJAM4QMAJKMw3wPA0XZ/47Z8j9BnS+ZNzfcIAElwxgcASIbwAQCSIXwAgGQIHwAgGcIHAEiG8AEAkiF8AIBkCB8AIBnCBwBIhvABAJIhfACAZAgfACAZwgcASIbwAQCSIXwAgGQIHwAgGcIHAEiG8AEAkiF8AIBkCB8AIBnCBwBIhvABAJIhfCKivr4+ysrKYubMmfkeBQAYQMInImpra6OtrS2am5vzPQoAMIAK8z0AMDjd37gt3yMA9JkzPgBAMoQPAJAM4QMAJEP4AADJED4AQDKEDwCQDOEDACRD+AAAyRA+AEAyhA8AkAzhAwAkQ/gAAMkQPgBAMoQPAJAM4QMAJEP4AADJED4AQDIK8z0AEHF/47Z8jwBH1WD8mV4yb2q+R+Bj4IwPAJAM4QMAJEP4AADJED4AQDKEDwCQDOEDACRD+AAAyRA+AEAyhA8AkIwhGz47d+6MSy65JMrKymLatGmxevXqfI8EAOTZkP2TFYWFhbFixYqoqKiI9vb2mD59elx22WVxwgkn5Hs0ACBPhmz4TJgwISZMmBARESUlJTF+/Ph47bXXhA8AJKxfb3X997//jS9+8YtxyimnxPHHHx/nn39+bNy48agNtW7duli4cGGUlpZGQUFBrFmzptf96uvr44wzzoiRI0dGZWVlbNiwodf9WlpaoqurKyZNmnTUZgQABp8+h8///ve/uOiii+K4446LP/zhD9HW1hbf+9734qSTTup1/2eeeSbefvvtw7a3tbVFR0dHr8/p7OyM8vLyqK+vf985Ghoaoq6uLpYtWxabNm2K8vLymD9/fuzevbvHfq+99lpcd9118cgjj/ThqwQAhqKCLMuyvjxh6dKl8cwzz8Tf/va3D903l8vFBRdcEGeffXasWrUqhg8fHhERW7dujYsvvjjq6urijjvu+OABCwriiSeeiEWLFvXYXllZGTNnzowHH3yw+7UmTZoUt9xySyxdujQiIg4ePBjz5s2LG264Ib70pS996Lz79u2LsWPHxt69e2PMmDEfun9f3d+47agfE4CjY8m8qfkegX7qy+/vPp/x+c1vfhMzZsyIz3/+83HaaafFpz/96Xj00Ud7P/iwYfH73/8+Nm/eHNddd13kcrl46aWXoqqqKhYtWvSh0fN+Dh06FC0tLVFdXd3jtaqrq2P9+vUREZFlWVx//fVRVVX1odFTX18fZWVlMXPmzH7NAwAMDn0On3//+9/x0EMPxdlnnx1/+tOf4qabboqvfe1r8bOf/azX/UtLS+PJJ5+Mp59+Oq655pqoqqqK6urqeOihh/o99J49e6KrqyuKi4t7bC8uLo729vaIeOcttoaGhlizZk1UVFRERUVFbNmypdfj1dbWRltbWzQ3N/d7JgDg2NfnT3XlcrmYMWNG3HPPPRER8elPfzpeeOGFePjhh6OmpqbX50yePDkee+yxuPjii+Oss86KH//4x1FQUHBkk3+Iz3zmM5HL5Qb0NQCAwaXPZ3wmTJgQZWVlPbade+658fLLL7/vczo6OuLGG2+MhQsXxhtvvBFLlizp+6TvMX78+Bg+fPhhF0d3dHRESUnJER0bABi6+hw+F110UWzdurXHtm3btsWUKVN63X/Pnj0xd+7cOPfcc+Pxxx+PpqamaGhoiNtvv71/E0dEUVFRTJ8+PZqamrq35XK5aGpqitmzZ/f7uADA0Nbnt7qWLFkSF154Ydxzzz3xhS98ITZs2BCPPPJIrx8Xz+VysWDBgpgyZUo0NDREYWFhlJWVRWNjY1RVVcXpp5/e69mfAwcOxPbt27vv79ixI1pbW+Pkk0+OyZMnR0REXV1d1NTUxIwZM2LWrFmxYsWK6OzsjMWLF/f1SwIAEtHnj7NHRPzud7+LO++8M1588cU488wzo66uLm644YZe921sbIw5c+bEyJEje2zfvHlznHrqqTFx4sTDnvPUU0/FpZdeetj2mpqaWLlyZff9Bx98ML773e9Ge3t7VFRUxAMPPBCVlZV9/XK6+Tg7QLp8nH3w6svv736Fz1AlfADSJXwGrwH9d3wAAAYr4QMAJEP4AADJED4AQDKEDwCQDOEDACRD+AAAyRA+AEAyhA8AkAzhAwAkQ/gAAMkQPgBAMoQPAJAM4QMAJEP4AADJED4AQDKEDwCQDOEDACRD+AAAyRA+AEAyhA8AkAzhAwAkQ/gAAMkQPgBAMoQPAJAM4QMAJEP4AADJED4AQDKEDwCQDOEDACRD+AAAyRA+AEAyhA8AkAzhAwAkQ/gAAMkQPhFRX18fZWVlMXPmzHyPAgAMIOETEbW1tdHW1hbNzc35HgUAGEDCBwBIhvABAJIhfACAZAgfACAZwgcASIbwAQCSIXwAgGQIHwAgGcIHAEiG8AEAkiF8AIBkCB8AIBnCBwBIhvABAJIhfACAZAgfACAZwgcASIbwAQCSIXwAgGQIHwAgGcIHAEiG8AEAkiF8AIBkCB8AIBnCBwBIhvABAJIhfACAZAgfACAZwgcASIbwAQCSIXwAgGQIHwAgGcIHAEiG8AEAkiF8AIBkCB8AIBnCBwBIhvABAJIhfACAZAgfACAZwgcASEZhvgcAgGPB/Y3b8j1Cny2ZNzXfIww6zvgAAMkQPgBAMoZs+OzcuTMuueSSKCsri2nTpsXq1avzPRIAkGdD9hqfwsLCWLFiRVRUVER7e3tMnz49LrvssjjhhBPyPRoAkCdDNnwmTJgQEyZMiIiIkpKSGD9+fLz22mvCBwASdkRvdd17771RUFAQt91221Ea5x3r1q2LhQsXRmlpaRQUFMSaNWt63a++vj7OOOOMGDlyZFRWVsaGDRt63a+lpSW6urpi0qRJR3VOAGBw6Xf4NDc3xw9/+MOYNm3aB+73zDPPxNtvv33Y9ra2tujo6Oj1OZ2dnVFeXh719fXve9yGhoaoq6uLZcuWxaZNm6K8vDzmz58fu3fv7rHfa6+9Ftddd1088sgjH+GrAgCGsn6Fz4EDB+Laa6+NRx99NE466aT33S+Xy0VtbW1cc8010dXV1b1969atUVVVFT/72c96fd6CBQvi7rvvjiuvvPJ9j33ffffFDTfcEIsXL46ysrJ4+OGHY9SoUfGTn/yke5+DBw/GokWLYunSpXHhhRf24ysFAIaSfoVPbW1tXH755VFdXf3BBx82LH7/+9/H5s2b47rrrotcLhcvvfRSVFVVxaJFi+KOO+7o19CHDh2KlpaWHq8/bNiwqK6ujvXr10dERJZlcf3110dVVVV86Utf+sDj1dfXR1lZWcycObNf8wAAg0Ofw2fVqlWxadOmWL58+Ufav7S0NJ588sl4+umn45prromqqqqorq6Ohx56qM/DvmvPnj3R1dUVxcXFPbYXFxdHe3t7RLzzFltDQ0OsWbMmKioqoqKiIrZs2dLr8Wpra6OtrS2am5v7PRMAcOzr06e6du7cGbfeems0NjbGyJEjP/LzJk+eHI899lhcfPHFcdZZZ8WPf/zjKCgo6POwffGZz3wmcrncgL4GADC49OmMT0tLS+zevTsuuOCCKCwsjMLCwli7dm088MADUVhY2OM6nvfq6OiIG2+8MRYuXBhvvPFGLFmy5IiGHj9+fAwfPvywi6M7OjqipKTkiI4NAAxdfQqfuXPnxpYtW6K1tbX7NmPGjLj22mujtbU1hg8ffthz9uzZE3Pnzo1zzz03Hn/88WhqaoqGhoa4/fbb+z10UVFRTJ8+PZqamrq35XK5aGpqitmzZ/f7uADA0Nant7pGjx4d5513Xo9tJ5xwQpxyyimHbY94J0YWLFgQU6ZMiYaGhigsLIyysrJobGyMqqqqOP3003s9+3PgwIHYvn179/0dO3ZEa2trnHzyyTF58uSIiKirq4uampqYMWNGzJo1K1asWBGdnZ2xePHivnxJAEBCBvRfbh42bFjcc889MWfOnCgqKureXl5eHn/5y1/i1FNP7fV5GzdujEsvvbT7fl1dXURE1NTUxMqVKyMi4uqrr45XX3017rrrrmhvb4+Kior44x//eNgFzwAA7yrIsizL9xDHin379sXYsWNj7969MWbMmKN+/Psbtx31YwKQriXzpuZ7hGNCX35/D9m/zg4A8H8JHwAgGcIHAEiG8AEAkiF8AIBkCB8AIBnCBwBIhvABAJIhfACAZAgfACAZwgcASIbwAQCSIXwAgGQIHwAgGcIHAEiG8AEAkiF8AIBkCB8AIBnCBwBIhvABAJIhfACAZAgfACAZwgcASIbwAQCSIXwAgGQIHwAgGcIHAEhGYb4HAAD65/7Gbfkeoc+WzJua19d3xgcASIbwAQCSIXwAgGQIHwAgGcIHAEiG8AEAkiF8AIBkCB8AIBnCBwBIhvABAJIhfACAZAgfACAZwgcASIbwAQCSIXwAgGQU5nuAY0mWZRERsW/fvgE5/ludBwbkuAAwWAzE79h3j/nu7/EPInzeY//+/RERMWnSpDxPAgBD0/8bwGPv378/xo4d+4H7FGQfJY8SkcvlYteuXTF69OgoKCj4yM/bt29fTJo0KXbu3BljxowZwAl5lzXPD+ueH9Y9P6x7fvRn3bMsi/3790dpaWkMG/bBV/E44/Mew4YNi4kTJ/b7+WPGjPE/x8fMmueHdc8P654f1j0/+rruH3am510ubgYAkiF8AIBkCJ+jYMSIEbFs2bIYMWJEvkdJhjXPD+ueH9Y9P6x7fgz0uru4GQBIhjM+AEAyhA8AkAzhAwAkQ/gAAMkQPkeovr4+zjjjjBg5cmRUVlbGhg0b8j3SkLJu3bpYuHBhlJaWRkFBQaxZs6bH41mWxV133RUTJkyI448/Pqqrq+PFF1/Mz7BDyPLly2PmzJkxevToOO2002LRokWxdevWHvu89dZbUVtbG6ecckqceOKJ8bnPfS46OjryNPHQ8NBDD8W0adO6/+G22bNnxx/+8Ifux635wLv33nujoKAgbrvttu5t1n1gfOtb34qCgoIet09+8pPdjw/UugufI9DQ0BB1dXWxbNmy2LRpU5SXl8f8+fNj9+7d+R5tyOjs7Izy8vKor6/v9fHvfOc78cADD8TDDz8czz33XJxwwgkxf/78eOuttz7mSYeWtWvXRm1tbTz77LPR2NgYb7/9dnz2s5+Nzs7O7n2WLFkSv/3tb2P16tWxdu3a2LVrV1x11VV5nHrwmzhxYtx7773R0tISGzdujKqqqrjiiiviH//4R0RY84HW3NwcP/zhD2PatGk9tlv3gfOpT30qXnnlle7b008/3f3YgK17Rr/NmjUrq62t7b7f1dWVlZaWZsuXL8/jVENXRGRPPPFE9/1cLpeVlJRk3/3ud7u3vf7669mIESOyX/ziF3mYcOjavXt3FhHZ2rVrsyx7Z52PO+64bPXq1d37/POf/8wiIlu/fn2+xhySTjrppOxHP/qRNR9g+/fvz84+++yssbExu/jii7Nbb701yzI/6wNp2bJlWXl5ea+PDeS6O+PTT4cOHYqWlpaorq7u3jZs2LCorq6O9evX53GydOzYsSPa29t7fA/Gjh0blZWVvgdH2d69eyMi4uSTT46IiJaWlnj77bd7rP0nP/nJmDx5srU/Srq6umLVqlXR2dkZs2fPtuYDrLa2Ni6//PIe6xvhZ32gvfjii1FaWhpnnXVWXHvttfHyyy9HxMCuuz9S2k979uyJrq6uKC4u7rG9uLg4/vWvf+VpqrS0t7dHRPT6PXj3MY5cLpeL2267LS666KI477zzIuKdtS8qKopx48b12NfaH7ktW7bE7Nmz46233ooTTzwxnnjiiSgrK4vW1lZrPkBWrVoVmzZtiubm5sMe87M+cCorK2PlypVxzjnnxCuvvBLf/va3Y86cOfHCCy8M6LoLH+AD1dbWxgsvvNDjvXcGzjnnnBOtra2xd+/e+NWvfhU1NTWxdu3afI81ZO3cuTNuvfXWaGxsjJEjR+Z7nKQsWLCg+7+nTZsWlZWVMWXKlPjlL38Zxx9//IC9rre6+mn8+PExfPjww64w7+joiJKSkjxNlZZ319n3YODcfPPN8bvf/S7++te/xsSJE7u3l5SUxKFDh+L111/vsb+1P3JFRUXxiU98IqZPnx7Lly+P8vLy+P73v2/NB0hLS0vs3r07LrjggigsLIzCwsJYu3ZtPPDAA1FYWBjFxcXW/WMybty4mDp1amzfvn1Af96FTz8VFRXF9OnTo6mpqXtbLpeLpqammD17dh4nS8eZZ54ZJSUlPb4H+/bti+eee8734AhlWRY333xzPPHEE/Hkk0/GmWee2ePx6dOnx3HHHddj7bdu3Rovv/yytT/KcrlcHDx40JoPkLlz58aWLVuitbW1+zZjxoy49tpru//bun88Dhw4EC+99FJMmDBhYH/ej+jS6MStWrUqGzFiRLZy5cqsra0tu/HGG7Nx48Zl7e3t+R5tyNi/f3+2efPmbPPmzVlEZPfdd1+2efPm7D//+U+WZVl27733ZuPGjct+/etfZ88//3x2xRVXZGeeeWb25ptv5nnywe2mm27Kxo4dmz311FPZK6+80n174403uvf5yle+kk2ePDl78skns40bN2azZ8/OZs+encepB7+lS5dma9euzXbs2JE9//zz2dKlS7OCgoLsz3/+c5Zl1vzj8t5PdWWZdR8oX//617Onnnoq27FjR/bMM89k1dXV2fjx47Pdu3dnWTZw6y58jtAPfvCDbPLkyVlRUVE2a9as7Nlnn833SEPKX//61ywiDrvV1NRkWfbOR9q/+c1vZsXFxdmIESOyuXPnZlu3bs3v0ENAb2seEdlPf/rT7n3efPPN7Ktf/Wp20kknZaNGjcquvPLK7JVXXsnf0EPAl7/85WzKlClZUVFRduqpp2Zz587tjp4ss+Yfl/8bPtZ9YFx99dXZhAkTsqKiouz000/Prr766mz79u3djw/UuhdkWZYd2TkjAIDBwTU+AEAyhA8AkAzhAwAkQ/gAAMkQPgBAMoQPAJAM4QMAJEP4AADJED4AQDKEDwCQDOEDACRD+AAAyfj/Gh4j2irtAk8AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(L3_days, alpha=0.5, bins=10)\n",
    "plt.yscale('log')\n",
    "np.median(L3_days)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 528,
   "id": "ea4e95ec-16fb-44da-bae3-f2e0ccd06b40",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "11703"
      ]
     },
     "execution_count": 528,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(L1_days)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 529,
   "id": "262c6840-cf61-4ad1-9934-42e9544184c4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8966"
      ]
     },
     "execution_count": 529,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(L2_days)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 530,
   "id": "d6af9999-0269-4643-9909-ed6dc3e376c3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7705"
      ]
     },
     "execution_count": 530,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(L3_days)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}