src/hawkes_experiments.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d20cbf57",
   "metadata": {},
   "source": [
    "Dependencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2e7f19ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "import sys\n",
    "from hawkes.hawkes import hawkes, hawkes_calculate, sampleHawkes, plotHawkes, iterative_sampling, extract_samples, sample_counterfactual_superposition, check_monotonicity_hawkes\n",
    "from sampling_utils import thinning_T, return_samples\n",
    "from counterfactual_tpp import sample_counterfactual, check_monotonicity, distance, covariance\n",
    "from multiprocessing import cpu_count, Pool\n",
    "from tqdm import tqdm\n",
    "import utils"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "81d38778",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sample_counterfactual_superposition_plot(mu0, alpha, new_mu0, new_alpha, all_events, lambda_max, w, T):\n",
    "    \"\"\"Generates samples from the counterfactual intensity, and return the counterfactuals.\n",
    "    This is done in 3 steps:\n",
    "        1. First we calculate the counterfactual basedon the history in each exponential (created by superposiiton.).\n",
    "        2. Then we determine the events that were rejected in original intensity and accepted in intevened intensity (rej_acc).\n",
    "        3. Then we create a new exponential and sample for each rej_acc event.\n",
    "    \"\"\"\n",
    "    def constant1(x): return mu0\n",
    "    def constant2(x): return new_mu0\n",
    "    all_counterfactuals = {}\n",
    "    all_counterfactuals['mu'] = []\n",
    "    all_counterfactuals['exp'] = []\n",
    "    all_counterfactuals['new'] = []\n",
    "    count = 0\n",
    "    counterfactuals = {}\n",
    "    for t_i, events in all_events.items():\n",
    "        sample = list(events.keys())\n",
    "        if count == 0:\n",
    "            lambdas = [constant1(s) for s in sample]\n",
    "            counterfactuals_events, counterfactual_indicators = sample_counterfactual(\n",
    "                sample, lambdas, lambda_max, list(events.values()), constant2)\n",
    "            all_counterfactuals['mu'].extend(counterfactuals_events)\n",
    "        else:\n",
    "            def f(t): return alpha * np.exp(-w * (t - t_i))\n",
    "            def g(t): return new_alpha * np.exp(-w * (t - t_i))\n",
    "            lambdas = [f(s) for s in sample]\n",
    "            counterfactuals_events, counterfactual_indicators = sample_counterfactual(\n",
    "                sample, lambdas, lambda_max, list(events.values()), g)\n",
    "            all_counterfactuals['exp'].extend(counterfactuals_events)\n",
    "        count += 1\n",
    "        counterfactuals.update(\n",
    "            {sample[s]: counterfactual_indicators[s] for s in range(len(sample))})\n",
    "\n",
    "    rej_acc_events = {}\n",
    "    for events in list(all_events.values()):\n",
    "        for t_i in list(events.keys()):\n",
    "            if events[t_i] == False and counterfactuals[t_i] == True:\n",
    "                rej_acc_events[t_i] = True\n",
    "\n",
    "    new_events = {}\n",
    "    iterative_sampling(new_events, rej_acc_events, new_mu0,\n",
    "                       new_alpha, w, lambda_max, T)\n",
    "    # These are the additional counterfactuals sampled from the new intensity.\n",
    "    sampled_counterfactuals = list(new_events.keys())\n",
    "    sampled_counterfactuals.sort()\n",
    "    all_counterfactuals['new'].extend(sampled_counterfactuals)\n",
    "\n",
    "    # Combine all counterfactuals\n",
    "    real_counterfactuals = [k for k, v in counterfactuals.items() if v == True]\n",
    "    real_counterfactuals.extend(sampled_counterfactuals)\n",
    "    real_counterfactuals.sort()\n",
    "    real_counterfactuals = list(dict.fromkeys(real_counterfactuals))\n",
    "    return real_counterfactuals, all_counterfactuals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "aebb5826",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotHawkes_plot(tevs, l_0, alpha_0, w, T, resolution, label, legend):\n",
    "\n",
    "    tvec = np.arange(0, T, step=T / resolution)\n",
    "\n",
    "    n = -1\n",
    "    l_t = np.zeros(len(tvec))\n",
    "    for t in tvec:\n",
    "        n += 1\n",
    "        m = l_0(t)\n",
    "        l_t[n] = m + alpha_0 * np.sum(np.exp(-w * (t - tevs[tevs < t])))\n",
    "\n",
    "    plt.plot(tvec, l_t, label=label)\n",
    "\n",
    "    return tvec, l_t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "889dfb28",
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_monotonicity_hawkes(mu0, alpha, new_mu0, new_alpha, all_events, sampled_events, real_counterfactuals, w):\n",
    "    count = 0\n",
    "    monotonic = 1\n",
    "    def constant1(x): return mu0\n",
    "    def constant2(x): return new_mu0\n",
    "    for t_i, events in all_events.items():\n",
    "        sample = list(events.keys())\n",
    "        if count == 0:\n",
    "            for s in sample: \n",
    "                if constant2(s) >= constant1(s) and s in sampled_events:\n",
    "                    if s not in real_counterfactuals:\n",
    "                        return('NOT  MONOTONIC, 1')\n",
    "                        monotonic = 0\n",
    "                if constant2(s) < constant1(s) and s not in sampled_events:\n",
    "                    if s in real_counterfactuals:\n",
    "                        return('NOT  MONOTONIC, 2')\n",
    "                        monotonic = 0\n",
    "        else:\n",
    "            for s in sample:\n",
    "                def f(t): return alpha * np.exp(-w * (t - t_i))\n",
    "                def g(t): return new_alpha * np.exp(-w * (t - t_i))\n",
    "                if g(s) >= f(s) and s in sampled_events:\n",
    "                    if s not in real_counterfactuals:\n",
    "                        return('NOT  MONOTONIC, 3')\n",
    "                        monotonic = 0\n",
    "                if g(s) < f(s) and s not in sampled_events:\n",
    "                    if s in real_counterfactuals:\n",
    "                        return('NOT  MONOTONIC, 4')\n",
    "                        monotonic = 0\n",
    "        count += 1\n",
    "        \n",
    "    if monotonic == 1:\n",
    "            return('MONOTONIC')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "766ee3aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(5)\n",
    "w = 1\n",
    "mu0 = 1\n",
    "new_mu0 = 1\n",
    "alpha = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "e361e103",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.4412274868850414\n"
     ]
    }
   ],
   "source": [
    "# run this cell only one time if you want to reproduce the results for \"Intervention increases self-excitation\"\n",
    "# run this cell four times if you want to reproduce the results for \"Intervention decreases self-excitation\"\n",
    "new_alpha = alpha + np.random.normal(loc=0.0, scale=1, size=1)[0]\n",
    "print(new_alpha)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "2b295296",
   "metadata": {},
   "outputs": [],
   "source": [
    "lambda_max = 1.5\n",
    "T = 5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df57029c",
   "metadata": {},
   "source": [
    "Run (uncomment) the following four cells in case you want to generate new data and save them. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b009d657",
   "metadata": {},
   "outputs": [],
   "source": [
    "def counterfactual1(_):\n",
    "    initial_sample, indicators = thinning_T(0, lambda t: mu0, lambda_max, T)\n",
    "    events = {initial_sample[i]: indicators[i] for i in range(len(initial_sample))}\n",
    "    all_events = {}\n",
    "    all_events[mu0] = events\n",
    "    iterative_sampling(all_events, events, mu0, alpha, w, lambda_max, T)\n",
    "    sampled_events = list(all_events.keys())[1:]\n",
    "    sampled_events.sort()\n",
    "    sampled_events = np.array(sampled_events)\n",
    "    counters = []\n",
    "    for counter in range(50):\n",
    "        real_counterfactuals = sample_counterfactual_superposition(mu0, alpha, new_mu0, new_alpha, all_events, lambda_max, w, T)\n",
    "        counters.append(real_counterfactuals)\n",
    "    return sampled_events, counters, all_events"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e2b2c4aa",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1000/1000 [06:24<00:00,  2.60it/s]\n"
     ]
    }
   ],
   "source": [
    "# Run in case you want to generate new data\n",
    "# with Pool(48) as pool:\n",
    "#     result = list(tqdm(pool.imap(counterfactual1, list(range(1000))), total = 1000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a0e7de06",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save\n",
    "# all_events_save = [result[i][2] for i in range(1000)]\n",
    "# import json\n",
    "# # json.dump(all_events, output_file)\n",
    "# with open('data_hawkes/allevents.json', 'w') as fout:\n",
    "#     json.dump(all_events_save , fout)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db982302",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# save = [result[i][0] for i in range(1000)]\n",
    "# np.save('data_hawkes/hawkes_1000_100.npy', save, allow_pickle=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f32a7d4",
   "metadata": {},
   "source": [
    "The following cells will read the same data we used (and generated) in our paper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e0a0fa27",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read\n",
    "import json\n",
    "with open(r'data_hawkes/allevents.json', \"r\") as read_file:\n",
    "    data = json.load(read_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "079dbfac",
   "metadata": {},
   "outputs": [],
   "source": [
    "new_data = []\n",
    "for dic in data:\n",
    "    new_dic = {}\n",
    "    for key, value in dic.items():\n",
    "        sub_dic = {}\n",
    "        for k, v in value.items():\n",
    "            sub_dic[float(k)] = v\n",
    "        new_dic[float(key)] = sub_dic\n",
    "    new_data.append(new_dic)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "431e8dff",
   "metadata": {},
   "outputs": [],
   "source": [
    "def counterfactual2(_):\n",
    "    all_events = new_data[_]\n",
    "    counters = []\n",
    "    all_counters = []\n",
    "    for counter in range(100):\n",
    "        real_counterfactuals, all_counterfactuals = sample_counterfactual_superposition_plot(mu0, alpha, new_mu0, new_alpha, all_events, lambda_max, w, T)\n",
    "        counters.append(real_counterfactuals)\n",
    "        all_counters.append(all_counterfactuals)\n",
    "        sampled_events = list(all_events.keys())[1:]\n",
    "        sampled_events.sort()\n",
    "        sampled_events = np.array(sampled_events)\n",
    "        if check_monotonicity_hawkes(mu0, alpha, new_mu0, new_alpha, all_events, sampled_events, real_counterfactuals, w) != 'MONOTONIC':\n",
    "            print('Not Monotonic')\n",
    "    return all_events, counters, all_counters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "74b42139",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|███████████████████████████████████████████████████████████████████████████████| 1000/1000 [02:49<00:00,  5.91it/s]\n"
     ]
    }
   ],
   "source": [
    "with Pool(48) as pool:\n",
    "    result = list(tqdm(pool.imap(counterfactual2, list(range(1000))), total = 1000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "88eb1b18",
   "metadata": {},
   "outputs": [],
   "source": [
    "b = np.load('data_hawkes/hawkes_1000_100.npy', allow_pickle=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "5e6259bf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1000\n",
      "59\n"
     ]
    }
   ],
   "source": [
    "def count_interval(start, end, samples):\n",
    "    return len(samples[(start <= samples) & (samples < end)])\n",
    "\n",
    "number_of_events = [count_interval(0, T, b[i]) for i in range(len(b))]\n",
    "print(len(number_of_events))\n",
    "number_of_groups = 3\n",
    "quantile_indices = pd.qcut(number_of_events, number_of_groups, labels = range(number_of_groups)).to_numpy()\n",
    "print(max(number_of_events))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "38cca3a9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "364.0\n",
      "348.0\n",
      "288.0\n"
     ]
    }
   ],
   "source": [
    "groups = [[] for i in range(number_of_groups)]\n",
    "samples_groups = [[] for i in range(number_of_groups)]\n",
    "for i in range(1000):\n",
    "    k = quantile_indices[i]\n",
    "    groups[k].extend(result[i][1])\n",
    "    samples_groups[k].append(b[i])\n",
    "for i in range(number_of_groups):\n",
    "    print(len(groups[i]) / 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e7b8d4e",
   "metadata": {},
   "source": [
    "In the following cells we start plotting the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "fc3f18a9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['5.04', '17.17', '36.17']\n",
      "0, 9\n",
      "******************\n",
      "10, 22\n",
      "******************\n",
      "25, 59\n"
     ]
    }
   ],
   "source": [
    "plot_groups = [[] for i in range(number_of_groups)]\n",
    "for i in range(1000):\n",
    "    k = quantile_indices[i]\n",
    "    plot_groups[k].append((b[i], result[i][1], result[i][2], result[i][0]))\n",
    "#***********************************************************************************    \n",
    "means = []\n",
    "for i in range(number_of_groups):\n",
    "    sum = 0\n",
    "    for j in range(len(samples_groups[i])):\n",
    "        sum += len(samples_groups[i][j])\n",
    "    means.append(sum / len(samples_groups[i]))\n",
    "means = [ '%.2f' % elem for elem in means ]\n",
    "print(means)\n",
    "#***********************************************************************************    \n",
    "print(min(len(samples_groups[0][j]) for j in range(len(samples_groups[0]))), end = ', ')\n",
    "print(max(len(samples_groups[0][j]) for j in range(len(samples_groups[0]))))\n",
    "print('******************')\n",
    "print(min(len(samples_groups[1][j]) for j in range(len(samples_groups[1]))), end = ', ')\n",
    "print(max(len(samples_groups[1][j]) for j in range(len(samples_groups[1]))))\n",
    "print('******************')\n",
    "print(min(len(samples_groups[2][j]) for j in range(len(samples_groups[2]))), end = ', ')\n",
    "print(max(len(samples_groups[2][j]) for j in range(len(samples_groups[2]))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "169d4c25",
   "metadata": {},
   "outputs": [],
   "source": [
    "d = []\n",
    "for i in range(number_of_groups):\n",
    "    for j in range(len(groups[i])):\n",
    "        d.append([means[i], len(groups[i][j])])\n",
    "        \n",
    "df = pd.DataFrame(d, columns=['Group', 'Number of Counterfactuals'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "74dd282d",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.to_csv('data_hawkes/higher_hawkes.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "318f49dd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.axis.YTick at 0x7f3040af8ca0>,\n",
       " <matplotlib.axis.YTick at 0x7f3040db0fa0>,\n",
       " <matplotlib.axis.YTick at 0x7f3040db04f0>,\n",
       " <matplotlib.axis.YTick at 0x7f30402e6b80>,\n",
       " <matplotlib.axis.YTick at 0x7f30402f9310>,\n",
       " <matplotlib.axis.YTick at 0x7f30402f9a60>,\n",
       " <matplotlib.axis.YTick at 0x7f30402fe1f0>]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAADkCAYAAAClrDJbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfqElEQVR4nO3dfZxVZb338c8A5kOgMGTGg6AS9zfL1JTM7NQJfdEtWppDWU75FHfouaVXZkmalZV5qjFLS285nMygHB8ySiowCVFPp0c1Tx7E351DjMKNog6giE8j+/5jrQ17hlnDzDB7r9l7vu/Xa16z1rXXrP3bzsv5ca3run5XXaFQwMzMrL8NyTsAMzOrTU4wZmZWFsPyDqA3Zs6cWbj++uvzDsPMzLLVFQ+qqgezYcOGvEMwM7MeqqoEY2Zm1cMJxszMysIJxszMysIJxsysDNra2rjwwgtpa2vLO5TcOMGYmZVBc3MzK1asoLm5Oe9QcuMEY2bWz9ra2li6dCmFQoGlS5cO2l6ME4yZWT9rbm5m69atAGzdunXQ9mKcYMzM+tny5ctpb28HoL29neXLl+ccUT6cYMzM+tnUqVMZNiwplDJs2DCmTp2ac0T5cIIxM+tnjY2NDBmS/HkdMmQIjY2NOUeUDycYM7N+Vl9fz7Rp06irq2PatGnU19fnHVIuqqrYpZlZtWhsbKS1tXXQ9l4gI8FIagK+DrwA3AEcCnwmIn5SwdjMzKyKZfVg3hcRcySdAqwGGoB7gW4TjKSRwA+AQ4AC8AkggFuAA9J7nRoRGyTVAVcDJwBbgLMi4oFd+zhmZgND6ULL2bNn5x1OLrLGYHZLv58I/DQiNvXwflcDd0TEm4DDgJXARcCyiJgMLEvPAaYDk9OvWcB1vQ/fzGzg8ULLRFaC+aWkR4AjgWWS9gVe7O5GkvYB3gNcDxARL0fERuBkYH562Xzgg+nxycCCiChExB+BkZLG7MJnMTMbELzQMpGVYC4FjgGmRMQrJI+wTtrJvQ4EngJukPRXST+Q9Fpgv4hYl17zBLBfejwOeLzk59ekbWZmVc0LLRNZCeYPEdEWEa8CRMTzwJKd3GsYcARwXUS8DXie7Y/DSO9TIBmb6TFJsyTdJ+k+72hpZtXACy0THQb5Jb2BpBexp6S3sX1v5b2BvXZyrzXAmoj4U3p+G0mCeVLSmIhYlz4CW5++vhbYv+Tnx6dtHUTEPGAeQENDQ6+Sk5lZHhobG1m6dCnghZal/ifwbZI/9t8Brky/LgC+0N2NIuIJ4HFJSpuOAx4GFgFnpm1nArenx4uAMyTVSToa2FTyKM3MrGp5oWWiQw8mIuYD8yXNiIif9eF+nwJulPQaYBVwNkkSu1XSTKAVODW9djHJFOVHScZ4zu7bRzAzG3i80BLqCoUdnzpJ2h2YQbJ2ZVsSioivVSyyLjQ0NBQWLlyYZwhmZta94tBK5kLL24FNwP3AS5WIyMzMaktWghkfEcdXNBIzM6spWdOUfy/prRWNxMzMakpWD+afgLMk/YPkEVkdUIiIQysWmZmZVbWsBDO9olGYmVnN6fIRWUS0kiyCPDY93pJ1rZmZWVe6TBqSLgU+D1ycNu3GTkr1m5mZlcrqlZxCUtzyeYCI+H/AiEoFZWZm1S8rwbxcWpgyrYpsZmbWY1kJ5lZJ/0ayR8sngd8C/165sMzMrNp1OYssIr4taRrwLCDgyxGxtKKRmZlZVesywUi6ALjFScXMzPoqax3MCOBOSW3ALcBPI+LJyoVlZmbVLmsdzFcj4i3AecAY4B5Jv61oZGZmVaytrY0LL7yQtra2vEPJzc4WT64HngCeAV5f/nDMzGpDc3MzK1asoLm5Oe9QcpO10PJ/S7obWAaMBj7pOmRmZj3T1tbGnXfeSaFQ4M477xy0vZisMZj9gfMj4sEKxmJmVhOam5t55ZVXAHjllVdobm5m9uzZOUdVeVljMBcDwyWdDSBpX0kHVjQyM7Mqddddd3V7PlhkTVO+FJhCsgbmBrbXIntXdzeTtBp4DngVaI+IKZLqSWaiHQCsBk6NiA2S6oCrgRNIimmeFREP7PpHMjPL1/Dhw3nhhRe2nY8YMTgrbZWjFtnUiDg8Iqak5xcByyJiMsmYzkVp+3Rgcvo1C7iu9+GbmQ08Tz31VIfz9evX5xRJvipRi+xkYH56PB/4YEn7gogoRMQfScrSjNmF9zEzswGkv2uRFUgWaN4vaVbatl9ErEuPnwD2S4/HAY+X/OyatM3MzGpAf9ci+6eIWCvp9cBSSY90um9BUqE3AaaJahbA2LFje/OjZma52H333XnppZc6nA9GWdOUSRNKr2qRRcTa9Pt6ST8HjgKelDQmItalj8CKDyPXkkyHLhqftnW+5zxgHkBDQ0OvkpOZWR5Kk0tX54NFv22DLOm1kkYUj4H3Af8NLALOTC87E7g9PV4EnCGpTtLRwKaSR2lmZlblMnswfbAf8HNJxfs2R8Qdkv5CMqYzE2gFTk2vX0wyRflRkmnKZ/djLGZmlrN+SzARsQo4rIv2Z4DjumgvkBTTNDOzGtQhwUh6iHRqcid1QMH1yMzMrKc692Den0sUZmZWczokmIhozSsQMzOrLVm1yI4Gvg8cDLwGGAo8HxF7VzA2M7Oq5HUwiaxpytcApwF/B/YE/hdwbaWCMjOrZsVS/Vnng0XmOpiIeBQYGhGvRsQNwPGVC8vMrHpt3bq12/PBImua8hZJrwEelNQErKMfF2WamVnty0oap5OMu8wmKdm/PzCjUkGZmVn1yyp2WZxN9gLw1cqFY2ZmtSJrFtk/6GLBZUQcVPaIzMysJmSNwUwpOd4D+DBQX/5wzMysVmQ9InumU9NVku4Hvlz+kMzMrBZkPSI7ouR0CEmPpj8rL5uZWY3LmkV2ZcnXN4Aj2F5m38z6QUtLCzNmzGDVqlV5h2JWFlm9kplp+f1tJB1YgXjMBo2mpia2bNlCU1MTc+fOzTscs36X1YO5rYdtZtYHLS0tPPbYYwC0tra6F2M1qfN+MG8C3gLsI6mh5KW9SWaTmVk/aGpq2uHcvRirNZ0fkYlkT5iRwAdK2p8DPlmhmMxqXrH3UtTa6p0yrPZ03g/mduB2Se+MiD/kFJNZzZswYUKHJDNx4sQcozErj6xB/nMlrYyIjQCSRgFXRsQndnZDSUOB+4C1EfH+dHLAzcBo4H7g9Ih4WdLuwALgSOAZ4CMRsXpXP5BZNZgzZw6zZ8/ucG5Wa7IG+Q8tJheAiNgAvK2H9/w0sLLk/FvAdyPijcAGYGbaPhPYkLZ/N73ObFCYNGkSEyZMAJLey0EHuQqT1Z6sBDMk7bUAIKmeHiy0lDQeOBH4QXpeBxzL9hlo84EPpscnp+ekrx+XXm82KMyZM4e99trLvRerWVlJ40rgD5J+CtQBHwIu78H9rgLmACPS89HAxohoT8/XAOPS43HA4wAR0S5pU3r906U3lDQLmAUwduzYHoRgVh0mTZrEz372s7zDMCubLnswEbGAZP+XJ4EngIaI+HF3N5L0fmB9RNzfnwFGxLyImBIRU0aNGrXzHzAzswGhuy2TVwC3AouAzZIm7ORe7wJOkrSaZFD/WOBqYKSkYk9pPLA2PV5LspEZ6ev7kAz2m5lZDegywUg6SdLfgX8A9wCrgSXd3SgiLo6I8RFxAPBR4K6I+BiwnOQRG8CZwO3p8aL0nPT1uyJihz1ozMysOmX1YC4Djgb+b0QcCBwH/LGP7/F54AJJj5KMsVyftl8PjE7bLwAu6uP9zcxsAMoa5H8lIp6RNETSkIhYLumqnt40Iu4G7k6PVwFHdXHNiyQbmZmZWQ3KSjAbJQ0H7gVulLQeeL5yYZmZWbXr8IgsXV0PyRqVF4DPAHcALXSsTWZmZtatzj2YP5BsLjY3Ik5P2+ZjZmbWS50TzGskNQLHdCrXD0BELKxMWGZmVu06J5hzgY+xY7l+gALgBGNmZj3SuVz/7yT9HlgTET0pDWNmZtalHdbBRMRWti+MNDMz65OshZbLJM1wdWMzM+urrHUw55Csrn9V0gskFZULEbF3xSIzM7Oq1mWCiYgRXbWbmZn1VJcJJn009jHgwIi4TNL+wJiI+HNFozMzs6qVNQbzf4B3Ao3p+Wbg2opEZGZmNSErwbwjIs4DXgSIiA3AayoWlZmZVb2sBPOKpKEkiyuRtC+wtWJRmQ0CbW1tXHjhhbS1teUdillZZCWY7wE/B14v6XLgd8A3KhaV2SDQ3NzMihUraG5uzjsUs7LoMsFExI3AHJKksg74YETcWsnAzGpZW1sbS5cupVAosHTpUvdirCZlbZn844h4JCKujYhrImKlpB9XOjizWtXc3MzWrclT561bt7oXYzUp6xHZW0pP0vGYI8sfjtngsHz5ctrb2wFob29n+fLlOUdk1v86rIORdDHwBWBPSc+SrOAHeBmY192NJO1BsgPm7ul9b4uISyUdCNwMjAbuB06PiJfTzc0WkCSuZ4CPRMTq/vpgZgPZ1KlT+c1vfkN7ezvDhg1j6tSpeYdk1u869GAi4hvpKv4rImLviBiRfo2OiIt3cq+XgGMj4jDgcOB4SUcD3wK+GxFvBDYAM9PrZwIb0vbvpteZDQqNjY0MGZL87zdkyBAaGxt38hNm1SerVMzFksYBE0uviYh7s24UEQWSBZkAu6VfBeBYti/YnA98BbiOZFvmr6TttwHXSKpL72NW0+rr65k2bRqLFy9m2rRp1NfX5x2SWb/LKhXzTeCjwMPAq2lzgeQRWKZ0rOZ+4I0kK/9bgI0R0Z5esgYYlx6PAx4HiIh2SZtIHqM93emes4BZAGPHju3FRzMb2BobG2ltbXXvxWpWVjXlUwBFxEu9uVlEvAocLmkkyTqaN+1aeBAR80jHfxoaGty7sZpRX1/PFVdckXcYZmWTNYtsFckjrj6JiI3AcpJ6ZiMlFRPZeGBterwW2B8gfX0fksF+s0GhpaWFGTNmsGrVqrxDMSuLrASzBXhQ0r9J+l7xq7sbSdo37bkgaU9gGrCSJNEUd8g8E7g9PV6UnpO+fpfHX2wwaWpqYsuWLTQ1NeUdillZZD0iW5R+9cYYYH46DjMEuDUifiXpYeBmSV8H/gpcn15/PfBjSY8CbSRjPmaDQktLC4899hgAra2trFq1ioMOOijnqMz6V12hUD2dhoaGhsLChQvzDsNsl51zzjnbEgzAxIkTmTt3bo4RWX+aPn36Dm1LlizJIZJcFNdPZs4i+wdpJeVSEeF/Ypn1g9LkAkkvxqzWZD0im1JyvAfwYcAT9c36yYQJE3bowZjVmh4/IpN0f0TkWo/Mj8isVrS0tDB79uxt59dee63HYHJw7rnnVn3vcQA+Xt3pI7IjSk6HkPRosno7ZtZLkyZNYty4caxdu5bx48c7ueSkXH+YB/kYzDZZSePKkuN2YDVwatmjMRtEhg8f3uG7Wa3JqkXm0q5mZdTW1kZEAPDII4/Q1tbmemQ1ZMmSJR16MYOx9wLZj8j2AS4F3pM23QN8LSI2VSows1rWuUTMlVdeyeWXX55TNGblkbWS/4fAcySPxU4FngVuqFRQZrXuwQcf7HD+wAMP5BOIlU2x1zJYey+QPQYzKSJmlJx/VdKDFYjHzMxqRFYP5gVJ/1Q8kfQu4IXKhGRmZrUgqwfzLyR1xfZJzzcAZ1UkIjMzqwlZs8geBA6TtHd6/mwlgzIzs+qXNYvsX4GmdF8XJI0CPhsRX6xgbGZmVsWyxmCmF5MLQERsAE6oSERmg0BdXV2352a1ICvBDJW0e/Ek3UBs94xrrUzOP/98pk+fzgUXXJB3KNbPOtcArKZtM8x6KmuQ/0ZgmaTi2pezgfmVCcmKiiu9V65cmXMkZma912UPJiK+BXwdODj9uiwivK9rBZ1//vkdzt2LMbNqk1khOSLuAO6oYCxWoth7KXIvxsyqTb+V4Je0P7AA2I9kN8x5EXG1pHrgFuAA0qrMEbFBUh1wNcnkgS3AWRHhehlmZjUia5C/L9pJpjK/GTgaOE/Sm4GLgGURMRlYlp4DTAcmp1+zgOv6MRYzM8tZhwQjaVn6/Vu9vVFErCv2QCLiOWAlMA44me0TBOYDH0yPTwYWREQhIv4IjJQ0pi8fohZJ6nB+8MEH5xSJmVnfdH5ENkbSMcBJkm6mZOtLgJ4+wpJ0APA24E/AfhGxLn3pCZJHaJAkn8dLfmxN2raupA1Js0h6OIwdO7Ynb18Trrrqqg77SXznO9/JMRozs97rnGC+DHwJGA90/otWAI7d2Q0lDQd+BpwfEc+W/ks8IgqSejXhPyLmAfMAGhoaBtViAUlEhHsvZlaVOiSYiLgNuE3SlyList7eTNJuJMnlxohYmDY/KWlMRKxLH4GtT9vXAvuX/Pj4tM1SV111Vd4hmJn1WVaxy8skncT2HS3vjohfdXejdFbY9cDKiCjt/SwCzgS+mX6/vaR9dvoo7h3AppJHaWZmVuWyil1+AziKZEU/wKclHRMRX+jmXu8CTgceKtmc7AskieVWSTOBVpIdMgEWk0xRfpRkmvLZu/A5zMxsgMlaB3MicHhEbAWQNB/4K0nC6FJE/I5OkwJKHNfF9QXgvF5Fa2ZmVaO7dTAjS473ybrIzMysK1k9mG8Af5W0nKRX8h62L5A0MzPbqaxilzeRrMZfSDIr7J0RcUslAzMzs+rWXbHLdSQzvczMzHqtP2uRmZmZbeMEY2ZmZbFDgpE0VNIjeQRjHbW0tDBjxgxWrVqVdyhmZr22Q4KJiFeBkDQhh3isRFNTE1u2bKGpyZuJmln1yRrkHwWskPRn4PliY0ScVJGojJaWFh577DEAWltbWbVqFQcddFDOUZmZ9VxWgvlSRaOwHXTutTQ1NTF37tycojEz672sYpf3SJoITI6I30raCxha2dAGt2Lvpai1tTWnSMzy0/jxRjY8syHvMHZJ6b5O1WjU6FE0/6S5Tz+bVezykySbfNUDk0g2AptLFzXFrDwmTJjQIclMnDgxx2jM8rHhmQ2ManSlqjxtaO57gs+apnweSXXkZwEi4u/A6/v8LtZrc+bM6fbczGygy0owL0XEy8UTScNIdrS0Cpk0aRJveMMbABgzZowH+M2s6mQlmHskfQHYU9I04KfALysXlgFs3rwZgOeeey7nSMzMei8rwVwEPAU8BJxDsjnYFysVlCXTlIsJZvPmzV5saWZVJ6ua8lZgPnAZ8FVgfrpBmFXIJZdc0u25mdlA12WCkXQi0AJ8D7gGeFRSdc+1qzKbNm3qcL5x48Z8AjEz66OshZZXAlMj4lEASZOAXwNLuruZpB8C7wfWR8QhaVs9cAtwALAaODUiNkiqA64GTgC2AGdFxAO7+oHMzGxgyBqDea6YXFKrgJ6MNP8IOL5T20XAsoiYDCxj+86Y04HJ6dcs4LoexmxmZlWgQw9GUkN6eJ+kxcCtJNOTPwz8ZWc3i4h7JR3Qqflk4L3p8XzgbuDzafuCdGznj5JGShqTbnQ26I0cObLDY7H6+vr8gjEz64POj8g+UHL8JPDP6fFTwJ59fI/9SpLGE8B+6fE44PGS69akbU4w7Djm0tbWlk8gg9y5555bsTI95SgpMnHiRNews9x0SDARcXY53ywiCpJ6NRtN0iySR2iMHTu2LHGZZSnXH+euksmSJd0OcZpVnaxaZAcCnyIZmN92TR/L9T9ZfPQlaQywPm1fC+xfct34tK2DiJgHzANoaGjwVGkzsyqRNcj/C5IZX98nmVFW/OqLRcCZ6fGZwO0l7WdIqpN0NLDJ4y82WHTurbj3YrUoa5ryixHxvd7eTNJNJAP6r5O0BrgU+CZwq6SZQCtwanr5YpIpyo+STFMu6+M5M6tOG5o37fwiG5CyEszVki4F7gReKjbubJ1KRJyW8dIOZf7T2WPn9TBOs5qzZMkSpk+f7t7LTrhcf752JcFnJZi3AqcDxwJb07ZCem5mZrZTWQnmw8BBpSX7zczMeiNrkP+/gZEVjMPMzGpMVg9mJPCIpL/QcQymL9OUzcxsEMpKMJdWNAozsy6MGj1ql/aEt103avSoPv9slwkmIu7p8x3NzPpJ80+a8w5hlwz2WYJZK/mfI5k1BvAaYDfg+YjYu1KBmZlVs2I5oMGcZLJ6MCOKx+m+LScDR1cqKDMzq35ZYzDbpAsif5EuvLxoZ9ebVdJpHzudjW1P5x3GLilHFeVKGVn/Om668cd5hzHgdP6dDtZeTNYjsoaS0yHAFODFikRk1gsb256mMOWsvMMYtDbe96O8Q7ABLKsHU7ovTDtJ4cuTyx6NmVkFVft+PzCw9/zJGoNx4Ukzq3ne76e8Om+Z/OVuri1ExGVljsfMzGpE5x7M811c81pgJjAacILppNq72AO5e91TdR4HMBuQOm+ZvG1TMUkjgE+T7NNyM33fcKymuYudv6od5C9NjFX6GZzcrTs7jMFIqgcuAD4GzAeOiAjXajAzs17pPAZzBdAAzAPeGhGbc4nKtm1GVXpuOxpZ/7ramCpbpZ9hZP3r8g7BBrDOPZjPklRP/iJwiaRiex3JIL9LxdiAUu5FfpUcYyuHWhhjq0Z77703zz777LbzkSNH5hdMjuoKhcLOryojSccDVwNDgR9ExDezrm1oaCgsXLiwV/c/47TTeGrjxl2K0fpu35EjWXDTTXmHMeB4jK32DeInEHXFg52WiiknSUOBa4FpwBrgL5IWRcTD/fUeT23cyGVPPdNft7Ne+lLeAZjlpNiLGay9F8g5wQBHAY9GxCoASTeTVAzotwQzdMgQvrTv6P66nfXS0CFZm6aa1bZbbrkl7xByl3eCGQc8XnK+BnhHf77Br3796/68XcUN1iJ5Zlb98k4wOyVpFjALYOzYsTlHs6NKDAKXs9quB4Hz4VmCNhjknWDWAvuXnI9P27aJiHkk06ZpaGjId0ZCF/zH2cysa3knmL8AkyUdSJJYPgo05huSWWW412K1LtcR2IhoB2YDvwFWArdGxIo8YzIzs/6Rdw+GiFgMLM47DjMz61+eQ2pmZmXhBGNmZmXhBGNmZmWRey2y3pD0FFC9lQf75nXA03kHYWXj329tG4y/36cj4niosgQzGEm6LyKm5B2HlYd/v7VtsP9+/YjMzMzKwgnGzMzKwglm4JuXdwBWVv791rZB/fv1GIyZmZWFezBmZlYWTjBmZlYWTjA5krQ57xis9yQVJP2k5HyYpKck/aqX97lb0pT0eLGkkf0cqpVJ5/93JZ0l6Zr0+FxJZ+zk57ddX8tyL3ZpVoWeBw6RtGdEvABMo9M+Rr0VESf0S2SWu4jwJlEpJ5gBRtLhwFxgL6AF+ASwG7AkIo6UdBjwIDAxIh6T1AK8NSK25BTyYLUYOBG4DTgNuAl4N4Ck1wLfBw4h+d19JSJul7QncANwGPAIsGfxZpJWA1OA4cCvIuKQtP1zwPCI+Iqku4G/pu/zWuAM4GLgrcAtEfHF8n5k6wlJXwE2R8S3Jb0duB7YCiwFphd/t8BYSXcAk4CfR8ScXAIuIz8iG3gWAJ+PiEOBh4BLI2I9sIekvUn+uNwHvFvSRGC9k0subgY+KmkP4FDgTyWvXQLcFRFHAVOBK9Kk8y/Alog4GLgUOLIP7/tyujJ8LnA7cB5JIjtL0ug+fxrrrT0lPVj8Ar6Wcd0NwDkRcTjwaqfXDgc+QvIPhI9I2p8a4wQzgEjaBxgZEfekTfOB96THvwfelZ7/a/r93cB/VDpOg4j4G3AASe+l835G7wMuSv/w3A3sAUwg+Z39pOTn/9aHt16Ufn8IWBER6yLiJWAVHbcft/J6ISIOL34BX+58QTqmNiIi/pA2NXe6ZFlEbIqIF4GHgYnlDDgPfkRWPe4lSSgTSf7l+nmgAPw6z6AGuUXAt4H3AqW9hzpgRkRE6cWSenLPdjr+w2+PTq+/lH7fWnJcPPf/z9Wl9Pf3KjX4+3MPZgCJiE3ABknvTptOB4q9mf8APg78PSK2Am3ACcDvKh6oFf0Q+GpEPNSp/TfApyTVAUh6W9p+L9CYth1C8mitsyeB10saLWl34P1lidzKLiI2As9Jekfa9NEcw8lFzWXMKrOXpDUl598BzgTmStqL5LHH2QARsTr9g3Vveu3vgPERsaGSAdt2EbEG+F4XL10GXAX8TdIQ4B8kieI64AZJK4GVwP1d3PMVSV8D/kwyM+2R8kRvFTIT+HdJW0n+sbgp53gqyqVizMzKRNLwiNicHl8EjImIT+ccVsW4B2NmVj4nSrqY5G9tK3BWvuFUlnswZmZWFh7kNzOzsnCCMTOzsnCCMTOzsvAgv1kvSBoPXAu8GRhKsor/s+lq+v56j/eSlIT5fXp+LkmJmQWSfkRSq+y2/no/s3JxD8ash9J1SAuBX0TEZGAyScHKpn5+q/cCxxRPImJuRCzo5/cwKzv3YMx67ljgxYi4ASAiXpX0GaBV0t+BN0XEbIB0b5hvR8Tdkq4D3k6SjG6LiEvTa1aT1Jv7AEnV5Q8DLwLnAq9K+jjwKeA40uq8pcFIOpJkce5w4GngrIhYV87/AGa94R6MWc+9hU6r7yPiWWA13f9j7ZK0AvKhwD9LKi0R83REHEGyyv9zEbGapFLyd9NCil0WM5W0G8mWAB+KiCNJytZc3qdPZVYm7sGYld+pkmaR/P82hmT8plhJeWH6/X6goRf3FEmZ/qVpEc2hgHsvNqA4wZj13MPAh0ob0j163gA8A/yPkpf2SF8/EPgc8PaI2JAO0pdWSC5ODuhtNd06knL97+zNBzCrJD8iM+u5ZSQFSs8AkDQUuBK4hqSg5eGShqQbRx2V/szeJFssb5K0HzC9B+/zHDBiJ9cEsK+kd6ax7CbpLb39QGbl5ARj1kMRUQBOAT6UDuo/A2yNiMuB/yRJMg+TVFh+IP2Z/yLZ5vgRkg2n/rMHb/VL4JR0t8R3d3VBRLxM0pv6lqT/ItlG+5iurjXLi2uRmfWRpGOAm4BTIuKBvOMxG2icYMzMrCz8iMzMzMrCCcbMzMrCCcbMzMri/wM5MmePeON9oAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 465.994x252 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df = pd.read_csv('data_hawkes/higher_hawkes.csv')\n",
    "width_pt = 397\n",
    "fig_height, fig_aspect = utils.get_fig_dim(width_pt, fraction=0.65)\n",
    "fig, ax = plt.subplots(figsize=(4*fig_aspect,3.5))\n",
    "plt.style.use(['science','no-latex'])\n",
    "plt.rcParams[\"axes.spines.right\"] = False\n",
    "plt.rcParams[\"axes.spines.top\"] = False\n",
    "plt.rcParams[\"xtick.top\"] = False\n",
    "plt.rcParams[\"ytick.right\"] = False\n",
    "sns.boxplot(x=\"Group\", y=\"Number of Counterfactuals\", data=df, palette=\"Set1\", whis=[5, 95]).set(\n",
    "    xlabel= \"Quantile\", \n",
    "    ylabel= \"Number of counterfactual events\"\n",
    ")\n",
    "sns.despine()\n",
    "ax.set_xticklabels(['Low', 'Medium', 'High'])\n",
    "ax.set_yticks(range(0,700,100))\n",
    "# fig.savefig('hawkes_high.pdf', bbox_inches = 'tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "66751c43",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAooAAABTCAYAAAAY2iZCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAIIklEQVR4nO3dbYhcZxXA8f8m4sbErBpIidCXVK0HjbEKbZbQQguNIOaTVarQWkO/GMGI9kOrWPuizYsxiZJEm6LRUkXQNkILEb9ECFjyZgWVgkeQtppoTSBKIiVR2vHD3Mhk8qS5zd65s5v5/2Bh73Ofneewh9k5e585d8Y6nQ6SJElSv1nDDkCSJEnTk4WiJEmSiiwUJUmSVPSG0mBELAIeBq7NzOsL52cB64CTwGJgZ2buH2CckiRJatn5rijeCDwFjJ3n/G3ARGauBe4FHo+I2QOIT5IkSUNSLBQz80m6VwvPZyWwr5p7HDgFLGk8OkmSJA1Nceu5hss4u5A8UY2dZdOmTZ3x8XEAli1bxuTk5EUuJ0mSpIacb8f4HBdbKB4F5vccT1RjZxkfH2fNmjUXuYQkSZKGqXbXc0TMi4iF1eFuYHk1vgCYAzzXfHiSJEkalvN1Pd8EfAp4e0TcB2wGVgFLgdXAz4APRsQDwJXAnZn5SisRS5IkqRXFQjEz9wJ7+4a/03P+VbrdzpIkSbpEecNtSZIkFQ20UDx8+DDbtm3jwIEDg1xGkiRJA3CxXc+1XH755XY9S5IkzVBuPUuSJKnIQlGSJElFFoqSJEkqslCUJElSkV3PkiRJKrLrWZIkSUVuPUuSJKnIQlGSJElFFoqSJEkqslCUJElSkV3PkiRJKrLrWZIkSUVuPUuSJKnIQlGSJElFFoqSJEkqslCUJElSkV3PkiRJKrLrWZIkSUVuPUuSJKmoeEUxIlYAtwJHgU5mPtR3fhWwGjhVDe3MzB8NME5JkiS17JxCMSLmAjuAJZl5OiJ2RcQtmbmnb+onM/OFNoKUJElS+0pXFJcDL2bm6er4GWAl0F8ofi4iXgLmAtsz8/jgwpQkSVLbSoXiZcDJnuMT1VivvcDuzDwWER8BngBu6X+gM13PAMuWLWNycrKRoCVJkjR4pULxKDC/53iiGvu/zHy+5/BXwNMRMTszX+mdZ9ezJEnSzFXqet4HXBUR49XxDcDuiFgQERMAEbE+Is4UmdcAL/QXiZIkSZrZzrmimJkvR8Rnga0RcQz4fWbuiYiNwHFgA/AS8EhEPA8sBe5oM2hJkiQN3lin0xnYg2/btq3j1rMkSdK0MlZ3ojfcliRJUpGf9SxJkqQiP+tZkiRJRW49S5IkqchCUZIkSUUWipIkSSqyUJQkSVKRXc+SJEkqsutZkiRJRW49S5IkqchCUZIkSUUWipIkSSqyUJQkSVKRXc+SJEkqsutZkiRJRW49S5IkqchCUZIkSUUWipIkSSqyUJQkSVKRXc+SJEkqsutZkiRJRW49S5IkqahYKEbEioj4bkQ8GBEPFM7PiYjtEfHliPhBRLy79DiHDx9uOl5NU769YHSY69FivkeHuR4dEXFz3bnnFIoRMRfYAXwxMx8E3h8Rt/RN+wLwl8xcD3wL2Fl68CNHjtSNQzPcwYMHhx2CWmKuR4v5Hh3meqTcXHdi6YricuDFzDxdHT8DrOybsxLYB5CZfwCujYiJ1x/n4AzzP6NhrT2q/w2OYq6HvfawjOrv21y79qVsVH/fMyXXpWaWy4CTPccnqrE6c070Tjp06ND+iDhTcL5QfbVlccvrTYe1h7UuwOLt27cPbW1GL9fDXNtcj9jaQ8r3Ykb09z3Mtc31yKw9p+7EUqF4FJjfczxRjb3eOWTm8rqBSJIkaXopbT3vA66KiPHq+AZgd0Qs6Nle3k13i5qIWAr8LjNPnPtQkiRJmqnGOp3OOYMR8SHg48Ax4L+Z+VBEbASOZ+aGiHgTsAn4O/AuYF1m/qnFuCVJkjRgxULx9YqIFcCtdLefO5n5UN/5OXQLyyPANcAGC8uZqUau7wUW0f0n4jrg/sz8Y+uBasoulOueebcDPwbmZ+a/WwxRDanxvB4Dznx6wmLgrZl5V6tBqhE1cn013dfrQ8AHgJ9k5tNtx6mpi4hFwMPAtZl5feH8LGAd3Z6TxcDOzNzfP2/KN9xu8nY6mt5q5vrNwN2ZuRHYBXyz3SjVhJq5JiLeA7y35fDUoJq5vgP4V2Zuzcy7gW+3G6WaUDPX9wC/zswNwDeAze1GqQbdCDwFjJ3n/G3ARGauBe4FHo+I2f2Tmvhklkvidjqq5YK5zsyvZuaZy9SzAK8wzUwXzHX1onMPULzSqBmjzt/w24EFEfH5iFiHz+uZqk6u/wEsrL5fCDzbUmxqWGY+ydl3qOnXW5sdB04BS/onNVEoTuV2OppZaucxIt4IfBq4r4W41Lw6uV4LfC0z/9NaVBqEOrm+iu6Vh63AY8AvS1ceNO3VyfUWYDIitgD3Az9sKTa1r9ZrehOFYmO309G0VyuPVZH4CPCVzPxzS7GpWa+Z64i4Angb8ImI+FI1fHdEXNdeiGpInef1CeAAQPX+8gngilaiU5Pq5Pox4PvVWww+Cvw0Iha0E55aVus1vYlC0dvpjI4L5rrajnwU2JKZz0bEx4YUq6bmNXOdmX/NzFWZuaF6LxN0c/6b4YSrKajzN3wP8A6Aamw28FLrkWqq6uT6CrrNiAD/BF6lmVpB00BEzIuIM28t6K3NFtC9Cfdz/T/TVNezt9MZETVy/XPgfcDfqh+ZV+q20vR3oVxXcxYCnwG+Xn09mpl+yPsMU+N5/RZgI/Ai8E5gV2b+YngR62LVyPWNdBtQfwtcDTybmTuGFrAuWkTcBNwJfJjuLt9m4C5gaWaurrqe1wMvA1cC3yt1PTdSKEqSJOnS4+VkSZIkFVkoSpIkqchCUZIkSUUWipIkSSr6HyBHFaPtbsVfAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 794.439x79.4439 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# some plotting params\n",
    "width_pt = 397\n",
    "fig_height, fig_aspect = utils.get_fig_dim(width_pt, fraction=0.65)\n",
    "fig, ax = plt.subplots(figsize=(fig_height*5,fig_height/2))\n",
    "plt.style.use(['science','no-latex'])\n",
    "plt.rcParams[\"axes.spines.right\"] = False\n",
    "plt.rcParams[\"axes.spines.top\"] = False\n",
    "plt.rcParams[\"xtick.top\"] = False\n",
    "plt.rcParams[\"ytick.right\"] = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "0dc3a03b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The following numbers indicate the samples we plotted in the paper.\n",
    "# i = 0, j = 0, k = 0\n",
    "# i = 1, j = 100, k = 20\n",
    "# i = 2, j = 0, k = 0\n",
    "# *************************\n",
    "# i = 0, j = 0, k = 0\n",
    "# i = 1, j = 0, k = 0\n",
    "# i = 2, j = 0, k = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "bb6d0a27",
   "metadata": {},
   "outputs": [],
   "source": [
    "# you can select j and k randomly\n",
    "i = 0\n",
    "k = 0\n",
    "j = 0 # from 0 to the length of the group\n",
    "all_events = plot_groups[i][j][3].copy()\n",
    "mu_events = [k for k, v in all_events[mu0].items() if v == True]\n",
    "s_plot = plot_groups[i][j][0]\n",
    "s = np.array(list(set(s_plot) - set(mu_events)))\n",
    "c = plot_groups[i][j][1][k]\n",
    "all_c = plot_groups[i][j][2][k].copy()\n",
    "all_c['new'] = list(set(all_c['new']) - set(all_c['exp']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "113b8ed3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'time')"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuwAAACaCAYAAADyxB0JAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABHfklEQVR4nO3deXxTVd748U+aNF3StKVQWpay0wsIZRFwV3BGFjdEcQUGH0flcWEUFZkRYRxHZ1BHHami4zg+/kRmlHF9HMFHZXFBERBZSy87dGFrC22aLmma+/sjC2npkqZJk6bf9+vVV5KTe889Oefc9Htvzj1Xp2kaQgghhBBCiPAUFeoCCCGEEEIIIRonAbsQQgghhBBhTAJ2IYQQQgghwpihsTcURRkD/At4SlXVt7zSjwG5Xot+p6rq40EroRBCCCGEEB1YgwG7oihTgRuB0gbe/lxV1duDWSghhBBCCCGEU2NDYjapqnobYGnLwgghhBBCCCHqavAMu6qq+U2sM1hRlFWAGcgBHlNVtSgYhRNCCCGEEKKja3QMexN2AQ/iPPv+LLBKUZSxqqrWmdD9L3/5ixYTEwPA2LFjOe+881pZVCGEEEIIIToMnftJiwN2VVXvcD9XFOUJoAwYC/zovVxMTAxz5szxv4hCCCGEEEKI1k3rqKqqFSgBegemOEIIIYQQQghvLQrYFUW5XFGU0V6vjUAnoDDQBRNCCCGEEEK0/Ax7L+AeRVHcY2rmAPupNxxGCCGEEEIIERgNBuyKopyrKMo6YATwW0VRPnS9tQbnAPhvFUVZD1wBXKOqak0blFUIIYQQQogOR6dpWvNL+WH+/Plaz549ZYYYIYQQQgghWs7/WWJ81bNnT5klRgghhBBCiFZq1SwxQgghhBBCiOCSgF0IIYQQQogwJgG7EEIIIYQQYSxoAXt+fj7Z2dn8+KPM+CiEEEIIIYS/5KJTIYQQQgghwpgMiWnA5s2bmThxItdddx1FRUV+5WGxWFi8eHGD7+Xn5/P666/jcDhaU0wRpqTt27+m2hB8b8fm8vHVXXfd5fevla1Zty3IPiGEEM2TgL0BBoOBN954g7vuuou1a9e2eH1N03jwwQeZMGECAHa7naysLE6ePAk4f30AeO211wJX6A6sfv2GUv22h7rlk7avK5zazq25NgTf9uGG8vHXiy++yNixY5tdbubMmXz44Yd10nxdN1RknxBCiOZJwN6AESNGkJGRwdChQ9m9e3eL1//888/R6XSMGjUKgAMHDmA2m0lNTfUsM2vWLJYtW8bx48cDVu6OqqH6DZX6bQ9nl0/a/oxwajs3X9oQmm/HhvLxV0JCAjqdrvkFA7xuW5F9QgghmiYBexMKCwvJzc1t8XqrVq1izJgxAGzfvp05c+bgcDiYOXMmTzzxBAAxMTEMGTKEL774IpBFDiuvbYbv8+qmfZ/nTA+Uxuq3NWpra1m+fDkzZ87klltuYfTo0ZSUlPi0rnfbN1a+9tD2q4pXcdWOqxi9ZTRX7biKVcWrAr6NcGs7N1/aEJrfh+vnA3DkyBHuuOMOZsyYwW233caWLVsAePjhhxk0aBC33HILNpuNa665hiuuuIL169fzxhtvcNFFF5GdnQ2Aw+Hg97//PbfeeivTp09nwYIFVFRU8Pzzz7N7925ef/11Zs6cybp16+qs+95773H55Zczd+5cFi1axNSpU7nrrruorq72lG/Hjh1cf/313HLLLZ5tTJo0idWrVwe1ztvDPiGEEKEUtItO3bPEjB07lvPOOy9Ymwma2tpa/vSnP1FaWtridXNycpg0aRIAWVlZTJgwAbvdzvz58+ss17t3b3JycgJS3nCUlQb3rYJXJsOFGc5g3f06YNtoon79oWkaDzzwAGlpabz55pscPXqUG264gZSUFJ/W9277psoXzm2/qngVTx15iiqtCoBjNcd46shTAEzuHLjGC7e2c/O1DaHpdqyfj91uZ/bs2fz6179m2rRp5ObmMmvWLFavXs3zzz9PQkICRUVFGI1GhgwZwsMPP0zXrl256KKL2Ldvnyefb7/9loKCAv71r38BcN9991FSUsLDDz/M1q1bmTp1Ktdffz0A48aN86x78803c+LECVasWMF//vMfzGYz1157LV9++SVXX301NpuN+++/n3nz5nH11Veze/dubrjhBp566il+8YtfBLXOm6tLIYTo6IJ2ht09S0x7DNYBVqxYgclkYsKECeTn57do3eLiYkwmk+e1qqooinLWciaTye+LWtuDCzOcwfl9q+D5H+oG74HUWP0CFBUVsWjRIqZNm8bChQv57rvvqKioQFVVli5detbyH374IcePH2fBggVER0eza9cuBg4c6HNZ6rd9Y+UL57Z/ufBlT7DuVqVV8XLhywHfVri03Z49e3jnnXcA39sQmm7H+vls27aNvLw8pkyZAsCgQYNIS0tj3bp1AMybN49du3Yxd+5cxo4dS9euXRvMNzExkT179rB+/XocDgcvvPAC3bt39+lzAgwfPpykpCSioqIYOHCg5/tt69atFBcXM3my86Bs8ODB9O/fv8m8AlXnEN77hBBChFrQzrC3xNjffcbugpafyfbV4B5JbPzzVT4vX1ZWRnZ2Nn/729/Izc0lNzfXc2GULzRNq/N6z549PPTQQ2ctp9PpIn5mhAszYMYwWLIRfjM28ME6NF6/AEuXLiUrK4spU6awY8cOXn75ZXbu3MnAgQN5/PHHz1r+3Xff5eabbyYqynksm5ubS2Zmps9lqd/2jZUvnNv+eE3D44gbS2+NcGm7DRs2eAJVX9sQmm7H+vm4x2ffcccdnjSbzYbFYgGcY80ffPBBFi5cyMKFCxst68iRI/njH//I3//+dx577DFuvvlmZs+e3cwnPCMhIcHzPCYmhpqaGgBOnjxJYmIier3e835ycnKTeQWqziG89wkhhAi1sAjYWxJMt4VXXnmF8ePHM2zYMDRN45tvvuGXv/ylz+unpKRgtVoBZ/BfVFREv379zlrOarXSpUuXgJU7HH2fB+/scAbr7+yAC3oGNmhvqn4BFi1a5Hl+7rnncvvttzeal91uZ+fOnTz33HOetI0bN3Lttdf6XB7vtm+qfOHc9mnRaRyrOdZgeiCFU9slJSV5Lij1tQ2h6Xasn096ejrR0dEsW7bMk1ZRUeEJdjVN4+uvv+YXv/gFTz31FC+88EKD+VosFsaOHctll13GkSNHuPPOO0lLS+OGG27w6bM2JjU1lbKyMux2OwaD81/D6dOnG10+kHUO4b1PCCFEqMlFp/UcOHCATz75xHM2LTMz0zMGdN++fUyYMIHly5ezcOFCsrOz+fe//8306dPrXGSlKAqFhYUA5OXlkZ6ejtFo5He/+12d5QoKChodDhAJvMesP3zBmeEx9S9EbY2m6tcfRqPRc93C2rVr2bp1q+eMoS/t7932TZUvnNv+/u73E6uLrZMWq4vl/u73B3Q74dR27mEq4HsbQtPtWD+f4cOH061bN8+FlXa7nfvuu49Dhw4B8K9//YvJkyfzxz/+kU2bNvHVV181mO+XX37Je++9B0CvXr1IS0vznJk2mUxUVlZy6NAhnnnmmRbV34gRI+jcuTMrV64EYPfu3Z6yNaaxOvdlX/GucwjvfUIIIUItaAG7+6LTcL5hR30bN27kzjvvRKfTef5ZPvroo/zwww+sXLmSAQMGkJyczE033cRVV12FxWLhxhtvJCMjo86ZtIkTJ3pmf8jIyMBoNDJ9+nR69OjhuRirpqaG7du3B2SO5nC1/XjdMevuMe3bAziyorH69YfBYODxxx/nkUceYfbs2RQUFKBpmifo86X9vdu+sfKFe9tP7jyZx3s9Tnp0Ojp0pEen83ivxwN6wSmEX9u5+dKG0Pw+XD8fvV7Pa6+9xnvvvceMGTP41a9+xVVXXcWgQYNYvHgxL774IocOHSI/P5+oqCgee+wxXn/9dd544w2+/fZbPvroI/79738zYsQIvv/+e371q18xbdo0evXq5Ql+b7jhBt5++20eeeQRLrvssjrrvvrqq3z00Ud8++23/POf/+S9997zvPfpp59iNBrJzs7mrbfe4tZbb+Wjjz4iKyur0Skhm6rzltS3L3UphBAdnqZpDf5lZmaOyczM3JeZmXl7vfQRmZmZP2RmZq7PzMz8NDMzs3ND6y9ZskSLNEVFRdpvfvMbTdM07Y033tDWrVunaZqmzZgxo85ydrtdu/XWW7Xc3NxG83rnnXe0Z555JniFFS1WUVGhVVVVeV5/+eWX2tSpUz2vfWl/afvQCETbufnShprWfDv6mk84OXXqVJ3XV155paeu6muqzltS35om+4QQQjTCE1c3eIZdUZSpwFygtF66EfgEeExV1YuALUCHuT3dtm3bGD58uOf5sGHDKC8vp6ioqM5MMnq9npdeeonly5c3mE9BQQGFhYXMnTu3TcotfPPBBx/w2WefAVBeXs7SpUuZNWuW531f2l/aPjQC0XZuzbUh+NaOvuQTbubNm+cZsrJz505Onjzpqbf6mqrzltS37BNCCNE8ndbAjAiKovRUVTVfUZR1wFuqqr7lSp8CvKiqaj/X6wzgEJCuqmqde4tnZ2drc+bMCW7phQigNWvW8NJLL5GYmIjdbuf666/nxhtvDHWxhA+k7QLjrbfe4n//93+Jj4/HZrMxd+5cLrjgggaXlToXQoig84xJbDBgd2sgYH8KGKWq6pVey1iAaaqq/p/3uhKwCyGEEEII4TdPwN7Si07TqDdMBjgNNHyHDyGEEEIIIUSr+DMPe0On5M+aRsA9SwzA2LFj2+0dT4UQQgghhAillgbsJ4CR9dKSXel19OzZExkSI4QQQgghROu0dEjMJsBzZwvXRafxwE+BLJQQQgghhBDCqaUB+yrAoCjKZa7XdwAf1p8hRgghhBBCCH8U2Urp9c10Ht/7P2y3HKCpCVI6isbmYT/XNUPMCOC3iqJ8CKCqajVwHbBYUZTvgNHAf7dJSYUQQgghRNg5WHGUvdaCgOVXVFNKXtVJbJqda3/+PUO+v5Pf73ubXeWHAraNJrdvqz+/Sug1Oa1ja8yfP1/r2bOnXHAqhBBCCBHBrt/6Bz46sZ5hCX2ZlnYJN6RdzBBTb3S6s+Yk8UlO+WGmbfsjORe9gaZpbCzNZcXxb1hx7GuSDCZuSr+Mm9MvQzFlBPiTQElNGZ3XTiPnwjcYnNAr4Pm3kKcC/Zklxidy0akQQgghROSLjTLy9tBH6RuXzgfHv2XylgWY9LHc0PVipqVdynBzvxYF7w7NQZQrVtXpdJyXPJjzkgfzXOZdbCjdzYpj33D55kfpEp3EtLSLua7rRQxN6OP3AYK3akcN6caUcAjW62jpGPYOYfPmzUycOJHrrruOoqKiVuVlsVhYvHhxo+/n5+fz+uuv43A4WrUdIYQQQohQcOAgWqfn4k5DeXHQPRy+5B3+39B52DQ7N2x7kgHf3c6je/7Oj6d349Caj3ccaETpzg5Ro3RRXJh8Dn8ddA95ly4ne9C9lNRYuObnRQz47nYeVv/Gd6d2UqvV+v9ZNI2oAAT+gSYBewMMBgNvvPEGd911F2vXrvU7H03TePDBB5kwYYInzW63k5WVxcmTzut0e/bsCcBrr73WukKLsCMHa+1XoNquuXzaC+mrQoimOIPcMyGlTqdjbNIgns28i30Xv8X7wxdi1Bm4fddf6PnNbdy960U+PfEDFbVVTeTXdNAcpYvi0pQsXhx0DwcveZsPhi/CrI/j/tyX6bbuFu7c9QL/ObmBqlpbyz4LZ87uhxMJ2BswYsQIMjIyGDp0KLt37/Y7n88//xydTseoUaM8aQcOHMBsNpOamupJmzVrFsuWLeP48eOtKndHVf8gKBzIwZpvIrntGsqnvZK+KoRoivcQlvp0Oh0jEwfw1MD/YvdF/+Dr0c8zyJTBC4c/JH3dLVz78yLeyF/FseqSM/nhIKoFIapOp2NEYn+eGPArtl7wGj+et4ShCX34y6H3Sf/6ZqZtfZL/V/AFJ6pP+fBZGj67H2rhV6IwUlhYSG5urt/rr1q1ijFjxnheb9++nTlz5uBwOJg5cyZPPPEEADExMQwZMoQvvviitUXukBo6CAo1OVjzTSS3XUP5tGcdva8K0VGdqD7FeRvm8OT+d9hcuqfBIS2NDWFpyEBTDx7qM421Y57j0KVvc0v6OL4q2cLg9XcydsMcnjqwnJ/L9rXqHHff+G482Pt61o35C3sv/h+uSj2PT09uIHP9HYzdMIcn9r3NplK14c/SxMFHKAUtYM/Pzyc7O5sff/wxWJsIqtraWv70pz+Rn5/vdx45OTlkZJy5gjkrK4sJEyZw3XXXsWzZMk/ADtC7d29ycnJaU+SwtHzDcvrM70PUXVH0md+H5RuWBzT/xg6C/LVo0SKGDh3Kzz//DMCDDz6IoigcPXq0Rfm0+4O1g8vh4z7wzyjn48HAthsEvu3Aud8uX76cmTNncssttzB69GhKSkqaX9FLoNqufj7vvfcel19+OXPnzmXRokVMnTqVu+66i+rqas8yR44c4Y477mDGjBncdtttbNmyBYDf/va3KIrC9OnTqa2t5dChQ1x88cVomsbx48e59tprmTRpUqOfNRD1ErZ9VQgRVMdsp9hYplJqtzJz5zN0+/oWfrXjWd49upaSmjIAajWHX+O+U6ITua3b5bybtYAT41aweOAdFNnKuD/3FSodLRvK0phUYzL/1WMi749YxIlxK3gm89dYa6uYtfM5un19C7N2PMuKY19zuqYccB586HX6gGw7kMJilpih6+9il/VwsIrCOabe7Lzo7y1aZ8WKFZhMJs477zzy8/M9Pwk3Z8+ePWzcuJEZM2ZQXFyMyWSq876qqlx55ZVnrWcymcjLy2tRGcPd8g3LuXvZ3VTYKgA4XHKYu5fdDcD086cHZBvugyC73c78+fNbnd+CBQv44IMPGDhwIAAPPfQQ69evp1u3bi3KJycnh0mTJvlUzrA7WDu4HDbeDbXOdqPisPM1QN/AtBsEvu00TeOBBx4gLS2NN998k6NHj3LDDTeQkpLSonwC1Xb187n55ps5ceIEK1as4D//+Q9ms5lrr72WL7/8kquvvhq73c7s2bP59a9/zbRp08jNzWXWrFmsXr2axYsX88MPP/DII4+g1+tZt24dJSUl7Nq1i6FDh3LNNddwwQUXNPhZA1UvzX1eIQLh9fzPSDYk8IuUkXQ2Joa6OAJnMD7c3I/nldk8r8zmYMVRVhVtYvnRNdyd8xJDE/pQWF3coiEsDYmOMnB555Fc3nkkLyr/jYPAXzNjjIpmfMoIxqeM4Dnlbs9nebvwK+7c9SIjzP0YlTgwLC86DVrA3hItDaaDraysjOzsbP72t7+Rm5tLbm6uzwH7hg0bmDx5MkCDd+bas2cPDz300FnpOp0u4i7oWvDRAk+w7lZhq2DBRwsCFrBD4wdBAEVFRSxZsoScnBwGDx7MxIkTGTVqFHl5eaxevZp77723zvL79+8nLS2NhIQET96ZmZk+lSNiDta2LTgTrLvVVjjTAxiwQ2Db7sMPP+T48eMsWbKEqKgodu3a5Tnwak4w2q6hfACGDx9OUlISAAMHDvT8irdt2zby8vKYMmUKAIMGDSItLY1169Zx9dVXc9lll7Fu3TpGjhzJjh07uOKKK/j6668ZOnQoOTk53HnnnQ2WI1D10tznFSIQZue8xCjzAO7c9SID47tzRedz+WXnkVycPJRYvTHUxeuQ6o8n7xvfjXt7Xcu9va6lqtbGN6e2s6ZkG8PN/QK2TZ1Oh57gn+X2/iyVtdWsLdnG/xVvZlKX0UHfdkuFRcAebl555RXGjx/PsGHD0DSNb775hl/+8pc+rZuUlOQZ55qSkoLVavW8V1ZWRlFREf36nd2prVYrXbp0CcwHCBNHSo60KN1fjR0EASxdupSsrCymTJnCjh07ePnll9m5cycDBw7k8ccfP2v5vXv31glmVFVFURSfyhExB2sVjbRPY+mtEMi2e/fdd7n55puJinL+Y8nNzfX5YCsYbdfYTencB4PgHGZSU1MD4Bkbfscdd3jet9lsWCwWAMaNG8df//pXZs+eTVxcHGPGjGH58uXMmDEDs9nc6PzDgaqX5j6vEIHQ1ZjMylFP0yk6gR9Lc/mq+GcW7XubHeUHOT9pkCeAH2HuH5YXBkYih6ahb6SuY/VGJnQZzYQwDHBbKk4fw5WpY7kydWyoi9IgCdjrOXDgAJ988gmfffYZAJmZmbz55psA7Nu3j3vvvZdZs2aRm5tL165dSU9P5+OPPyY7O5uUlBTP2TEARVEoLCz0vM7LyyM9PR2j0cjvfvc75s2b5/lZuqCggLFjw7OT+KtXSi8Ol5w91KlXSuBuRtDUQRA4x6S7nXvuudx+++1N5rd371569+7teb158+Y6wxqaEjEHa/G9nMNgGkoPoEC2nd1uZ+fOnTz33HOetI0bN3Lttdf6VJZgtF39fJqTnp5OdHQ0y5Yt86RVVFR4Au0LLriAuXPn8uGHH3LhhRdywQUXsGDBAj755BMuvfTSBvMMZL1AGPZVEXHc0/kZo6K5pNMwLuk0jD8M+BWlNVbWndrGV8U/c9v2P1NcY2F8ynDGdcrispSsVt1Vs6MqtpWxo/wgY5MU4vWxjS4XrhdhdjRy0amXjRs3cuedd6LT6fjqq68AePTRR/nhhx9YuXIlAwYMIDk5mZtuuomrrroKi8XCjTfeSEZGRoP/mCdOnOi5aAwgIyMDo9HI9OnT6dGjhydYr6mpYfv27REx/Zu3p6c+Tbwxvk5avDGep6c+HbBt1D8IaumFdPXt3bvXc2Z0y5YtbN682XOGfd++fUyYMIHly5ezcOFCsrOz+fe//8306dMpKSlp0cGadzkLCgp8PovfJoY/Dfq67YY+3pkeQIFuO6PRSGlpKQBr165l69atnjPJoWi7+vk0Z/jw4XTr1s1zUafdbue+++7j0KFDAMTFxTF27FiWLl3KJZdcQqdOnRg6dCivv/46F1xwQYvrpbk6AerUS3OfV4hAaGw6v6RoE1O6Xkj24PvIvfhNtlzwCld1Gcvmsr1c+/Pv6bruJqZtfZLsIx+z3XLAp5vzdHTvHlvH+M3z6LruJi7a+CCP7X2Tz4s2YbHXHRLZkhlgRPCExUWn4WLs2LGsWbOmTtqSJUs8z4uLi+nWrRvR0dHs2rWLiy66CHD+E/OeDcbtmmuuYcWKFZ5hFYmJiaxcufKs5VasWMF1113n8zj59sI9Tn3BRws4UnKEXim9eHrq0wEdv+59ENTYRXctsXfvXvbu3cv+/fuZPHkymqaxePFi3nzzzToHbD/99BNr1qxhzpw5/PTTT1it1jrbnjhxYp3ZNBorp/tgraEhHiHjHqe+bYFzGEx8L2ewHuDx64FsO4PBwOOPP84jjzxCv379uOSSS9A0zROwh6Lt6ufz6aef8tFHH1FdXc0///lP9Ho93377LTExMfTp04drrrmG1157jT/84Q+8/fbbOBwOrr/+egYNGuTJY9y4cVRXV2M2mz2vN23a1OBY+ebqxWw2+1wnvnxeIQLB17tMZsR2ZVaPCczq4TzRlVd1gq9LtvP1qR1kH/mE4poyLu00jMs6ZXFZpyyyzH3DcuaPUKrVapnTawqLB/6aH07v5utT21l88D02l+1hiKm3q/6GodPpwvIizA5H07Sg/C1ZskSLNKtXr9b+8Y9/aJqmaXPmzNGKi4s1i8WiTZo0ScvLy2twnRMnTmgLFy5sNM/8/Hzt2Wef1Ww2W1DKLHxXXl6unXPOOVp1dXWD7xcVFWm/+c1vNE3TtDfeeENbt26dpmmaNmPGjLOWtdvt2q233qrl5uY2uc133nlHe+aZZ1pZclFRUaFVVVV5Xn/55Zfa1KlTPa9D0Xa+5hNMTdVLS+pE06SviraRuHqKdtpW3up8CiqLtH8WrtFm7/qrNui7O7Sk1ddpEzf/TvvDvmXal0U/aaU1rd9Ge/fioQ+0B3YvPSu90l6tfV28Tfvj/ne0X256VIv/8mrtF5seDUEJheYVV8tvHC1w+eWXey4IW7JkCSkpKSQkJLBq1apGz46npqby5JNPNppnjx49mDdvHtHR0UEps/Ddvn376NGjB0ZjwzMRbNu2jeHDh3ueDxs2jPLycoqKis6ar1+v1/PSSy+xfHnj85cXFBRQWFjI3LlzA/chOqgPPvjAc91JeXk5S5cuZdasWZ73Q9F2vuQTbE3VS0vqRPpq+3Sw4ijHfbizYziwO2oB38+wN6d7bGdu7Tae14Y8wO6L/sHei/+HezOuoaK2mif3v0P3r29lxA//zb05S3in8CsOVBxt9ELxSOXQHA1eTBqrN3JpShaP95vOl6Of4fTlH/HhiEUN5CDaki5YHTQ7O1trb0NiRMf2/vvvs2bNGpYuXRrqoogWWrNmDS+99BKJiYnY7Xauv/56brzxxlAXK+SkXjq2Hl/fSmF1Mb1j0xibpHBe0iDOSxrEqMQBTV5k2Nas9koS1kwhM74n+yoKsfzi46CXz+aoYatlP9+fzuH70zmsP70Lu1bLhclDuDBpCGOTBnFu4kASDHFBLUeg1DjsROl0LRr285dD/+Z49SmeU+4OYslEK3mOXiVgF0IIIYJE07SQzV4y7Pu7WT7st8RGGfmxNJcfS3PZWKqyq/wQmaaenJc0yBPIDzJlhGyMd7GtjP7fzeLbMS9wsPIY13Zt/CLqYNE0jbyqk3x/ehffl+awqXQP2y0H6BuXzpgkhTGJmYxJUsgy9yUmKvzmg7/u59/zyckfuKxTFqMTMxmTlMnoxEz6xXVrtP89e3AFxTVlPJPZ8D0cRFjwNF7QLjp1zxIzduxYzjvvvGBtRgghhAhLR6uL6f71rZxj6s25iQMZlTiQcxMHMsLcv03O3Lrnz8409STT1JOZ3Z33E6mqtbHVsp+NpbmsLv6ZPx14lxO204x2lW9U4kBGJQ5gYHyPNpkdxIGDaJ2BYea+DDP3Dfr2GqLT6egV15VecV25pdt4wHkWfmf5ITaVqmwq28Pf8leyt6KAcxJ6ewL4sUlKSA923BxovHnOw/SI6czmsr28e2wdj6ivY62tZnTSQMYkKoxOHMiYJIUeMV2c91TAIReTtiMtDtgVRXkL6FMv+WpVVcu9E9rjLDFCCCFEoFjslfSOTePtYY+ypWwfP5Xt5Z9H17Cz/BC94royyjzAEyCPNA8gKbrh2X781dgUibF6I+cnD+b85MGetGJbGZvKVLaU7eOD49+yYN//cNJWynBzP0YlDmCU2RnEDzb1IjoqsOf6AjVuPdCMUdGug5eBzHalVdRW8XPZfjaW5fJl8Rb+dPBfHKs+RZa5LyPM/RlpHsAIc3+GJvRp0zuzOjSNVGPSWTcxOlZdwuayPWwq3cMbBZ8zO+clonRRjEnM5JS9nMs6ZbVZGUXr+LXXqao6LsDlEEIIISKKAwexXkHfnTjvGlvjsLPbeoSfyvaypWwf7x//ju3lB+gWk8Io80BGmPuRZe7HsIS+ZMSm+j2kpqk7VNbX2ZjIpC5jmNRljCftVI2Fn8v2scWyjy+Lf+KZQ+9xuPIEQxP6MCpxACPN/RnuCk5b84tBYwcW4SheH8tFnc7hok7neNJO1VjYZjnAVst+vju9k+wjn7CnIp8B8d0ZYe7vCuT7MyKxPynRiUEpV2N1mB6TwtWp53N16vnAmaE/7oOzyV7tLcKb3OlUCCGECALnmeOzg6joKANZrqD8v3pMBJxzYuda89hSto/tlgMsOfIx2y0HqXRUk5XQ17V8X7IS+vkcILd2yEOnaDOXdx7J5Z1HetLK7ZVssxxgi2UvG0pzeT1/JbuteXSLSWFYQh/PgcYwcx8GxPXAENX8UJFwPcPuq07RZsalDGdcynBPWrXDxq7yw2y17GerZT8fn/iebZYDJEebGGkeQJa5L0MT+jA0oQ+Z8T1b/auFr3XoPfTnhrRLWrXNSFZdU0tMdHjN2+9XD1EU5TXgHKAceFZV1bUBLZUQQgjRzrXklu56nZ5zEvpwTkKfOuknbafZYTnI9vKDbDjtDJBzyo/QLSbFE8APc603IL57ncDPoWkBP3OdYIg76wyz3VHLvsoCdlgOsaP8IMuPrmbH3kMcqy5hkCmDYea+ZCU4x6cPS+hLmrFTnV8NnGeH22/A3pCYKKPnlxU3h+bgQOVRfi7bz87yQ6w49g2Lyt/mSNUJBsR39wTwQxP6cE5Cb/rGpfs8Nj7S6lDTNLYeOkVGl3i6mNt2RqPTVhsZ97zP5sVXoXRPatNtN8WfgH03sEZV1U2KoowB1iqKcrGqqlu9F3JfdArIhadCCCHavcraauL0MT4vH4hbuqcak886y12r1bKvopDtloNsLz/AsqOrybEeJr+qiAHx3TknoTdDTL0ptVvb5My1IUrPIFMvBpl6cSOXetIt9gp2lR9mR/lBdpQf5NOTP7K9/AA6dAwx9WJIQm+GmHrRKdoc8os220KULooB8T0YEN+jTj1V1laTa81jZ/khdpYf4vX8lewsP0SRrZTBCb08QfwQU28GmzLoHZd2Vr9q7Nec9urAiXIu/f3nJMdHo9dHoXRPZFD3JOdjjySU7kl07xQXlBmYKm120pJiwypYBz8CdlVVn/F6vklRlM+A2cA93svJRadCCCEixXendnLJpodIN6YwyJSBYurp/IvPYJApg95xXc8KOoM11EOv06OYMlBMGQ0GfrvKD5FjPcLEzueSGh26oMNsiD/r4lZN0zhhO02O9TC7yg+TU36YHOsRRpj7haycoRanj2Fk4gBGJg6ok15mt5JTfoSd5c5fLv6vaDO51jxKaixkmnoyyJTBYFMvBpkyOFVjiagz7LaaWpTuiWz681UcL61CLSxFLSwjt6CUz7YUkFtYSmW1HaV7EpndE1G6JzGoh/OxT6oJfZT/By+1Dg19VPjVZSDGsB/BOTxGBInFYuGVV17ht7/9bYPv5+fns3LlSu68806iWtFJhRBCNKy8tpIrOo/iH+c8hGrNJ9eah2rNZ1XRJlRrPidsp+kf1w3FlOEJ6Gu1tr2YsrHAL5zodDrSYjqRFtOJ8SkjQl2csJZoMJ11wAPOXy5Uaz67rUfYbT3Cu8fWYdccdIvpHKKSwhfbCvn5YDH9083062qmX5qZZJP/s+S4g2adTkd6chzpyXFcNiS9zjKnrDbUwlJyC5zB/He5x1ELyzheWkmvLgn0T0ugf5qZgd0S6Z9mpn+amR4p8UQ1E4zXOjQM+vCLpfyZ1vFRVVWf9UpKAwoDV6TQs9vtjBo1itWrV5OamhrSsmiaxoMPPsh9993nSatfvp49ewLw2muvce+994aqqEIIEbEcmoZBpycjtisZsV35ZedRdd632ivZW1GIWuEM5N1nQweZMkJUYhGpzIZ4RidlMjopM9RF8Xj7m/3k5JcyNCOZ/ccsHDhhwWjQe4Lmfq6AuV9aAv3SzKQkND20rNbR/K9TnUxGzh+YyvkD68ZplTY7B0+Us++Yhf3HLfx8sIT3Nxxm3zELpRU2+qQmMCDdTP90Z5kGuB7Tk51DbCLpDPtDiqK8parqCUVR+gJTgGsDXK6zPPvsesaM6c748WduqrB27UE2bSrk0UcvCui2Dhw4gNlsDnmwDvD555+j0+kYNerMP4eGyjdr1izGjRvHDTfcQFpaWiiKKoQQEau5qQdNhjhGJDqn7hOio9HrdDw2dSjTzu8DOE82niyrYt8xCwdOlHPguIWVP+ez/7iF/ccsGPRRdQL43qkJ9ElNoE+qie6d4lsVNMcZDQzpmcyQnslnvVdeVcOB4+XsP25h37Eyfthzkne+PcC+YxaqbLX06ZpAF3NMWM5a5E/A/hfgI0VR7IAJuF9V1a8DW6yzjRnTnZtuep8VK6Yxfnxf1q496HkdSNu3b2fevHk4HA5mzpxJ//79eeKJJ/zOb9GiRXz44YcsW7aMkSNH8uCDD7Jq1SrWrVtHt27dml1/1apVjBlzZp7UxsoXExPDkCFD+OKLL5g5c6bf5RVCCHG29j71oBDBVFtv/9DpdHRNiqNrUhwXKl3rLKtpGkWWavYft3DguIUDx8v5etcx3jpZzuGTVk5Zq+mcEEP3lPiAlzMhNpqs3p3I6t3prPdKK2wcOlHOwRPlxMWE30XQ/lx0+hecQXuT3LPEBGqGmPHj+7JixTRuuul97rlnNK++utkTvAdSVlYWEyZMwG63M3/+/Fbnt2DBAj744AMGDnRO7fTQQw+xfv16n4J1gJycHCZNmuRT+Xr37k1OTk6ryyyEaFq1w4bdUYupDW4vL8JDrVYbURf1CRFILTkjrtPpSE2MJTUx9qzhLOAc0nL4pBWjoW3HkSfFGxneJ4XhfVLadLu+CtqNk4IxS8z48X25557R/PGP37Bw4aUBD9bdVFXlyiuvbPC9oqIilixZQk5ODoMHD2bixImMGjWKvLw8Vq9efdYY8v3795OWlkZCQoIn78xM38edFRcXYzLVvV11Y+UzmUzk5eX5nLcQwj//nbOEtwq/oJPBTI/YzvSM6ULP2FR6xnahZ0wXergee8amkmQwBWXqMREYx6pLsNgr6RHbmXh94/M9B2KKRiHaoyNFVs556BP6pJrokWKiR0pcvcd4rFX2Zi/m9FWc0cCgHuE1pWI4aFd3Ol279iCvvrqZhQsv5dVXNzN+fJ+gBO179uzhoYceavC9pUuXkpWVxZQpU9ixYwcvv/wyO3fuZODAgTz++ONnLb93717P2XVwBtuKovhcFk3TfC6fTqfD4XD4nLcQwj82Rw3Lhj7KpC5jyK8qoqC6iPyqk+RXF/Hd6V3kV52koLqYvKqTaGie4L1HTGdPYN8jpgvpMZ3oFpNCmrETxqjoUH+sDmnGjmdYXfIzsVFGYqOMdI/pTI/YzvSI6eJ8HtOZHrFdUK356CVgFx3QaauN/mlm3n/4MgpLKskvsVJQUsnuglK+2nGUwpJKjpVW0jWxbW9w1NG0m4Dde8z6+PF9GT++T53XgVJWVkZRURH9+jU8J+yiRYs8z88991xuv/32JvPbu3cvvXv39rzevHlznSEuzUlJScFqtfpUPqvVSpcuXXzOWwjhHwcaep2eLsYkuhiTGEHjFxqW2a3kVxV5gvj8qpNstezns5MbOWYr4Wh1CSdsp0kymEg3OgP49JgU56OxE91iOnsC+3RjipyxDzC9LorPR/2JCZ3P5ZTdQkFVMQXVRRRWF1NQVcT28oOsKtpEYXUx16SeH+riCtHmah0aCbEGBqQnMiA9MdTF6bDaTcC+aVNhneDcPaZ906bCgAbseXl5pKenYzQa+d3vfse8efNISfF/PNPevXvp1asXAFu2bGHz5s088MADAOzbt497772XWbNmkZubS9euXUlPT+fjjz8mOzublJQUFEWhsPDMrJlNla+goICxY8e24tMLIXzRkgsQEw0mhiSYGJLQu9FlHJqDIlspx2ynOFpdwrHqEo7ZTnGk6iQby1RXmvO9Gs1OurGTJ6hPNSbR1ZhMarTr0ZjsSesSnYQhKvwungonDs1BlM4533NKdCIp0YkMMwdnuKUQwZJbUEpesZUu5lg6m2PoYo4hPiYwIV6twxGW0xx2NO0mYG9o6kbnmfbAfrFmZGRgNBqZPn06F1xwQauCdXAG7Hv37mX//v1MnjwZTdNYvHgxb775JgMGDCA5OZmbbrqJn376iTVr1jBnzhx++uknrFYrKSkpTJw4kS+++KLZ8tXU1LB9+/YGh+UIIQKrqLySNSeOYTh0mJSEGFISjM5HcwzxRn2Lz4BH6aLoGtOJrjGdyGrmjo9WeyXHbKdcAXwxJ2tKOWkrZU9FAetP53DSdpoTttOcrCmluKaMJIOJ1OgkUo3JroA+yfPoHeR3MSbSOTqxww3NcaC16c2NRNNeWrmbIksVXcyxnv2qsznGs591MhlbdRfLSLXgX1v4cV8RfVITKLJUU2SpQq/T0SUxls4JMXRJdAbxnb0Ceu/gvktiLMnx0Q1+d9U6tICNTxf+C1rAHuhZYtpKYmIiK1euDEheVquVEydOsGXLFoxG5x2/brzxRs/7xcXFdOvWjejoaHbt2sVFFzkPSgoKCsjIcN5s45prrmHFihWese+NlW/FihVcd911npsoCSGCZ9/xMvbsLaQi/ggl5TZKyqtdfzY0NDonuAMMr2C+ieedTEaf/yGaDHH0N8TRP757s8vWarWcqinnpK3UGcR7BfOqNZ9vbTs5aSvlpCvtVE05xigDnaMTSYk2n/XYUFrnaDOdos1ER7Wb8z911LrOsAfLybIqYqP1JMQaZCiTDx5/92f++4pMjp2uJCf/NMWWM/tWSXk1ZZU1JMZFe/ahM8H8mcC+s1eg38lkJNlkJDa65QfS7YlOp+P12Rdw5UhnDKBpGtZquzN4L6uiyFJNsaXaE8wfOG7xvFdS7ky3VttJSXAG8N7fU5W2WjnDHgba1Swx7c2+ffvo0aOHJ1ivb9u2bQwfPtzzfOrUqZSXl1NUVER+fj49e/ZEr9fz0ksvkZ2dzZNPPtlgPgUFBRQWFjJ37tygfRYhxBkaGjdf0Je/XHrJWe9VVNvrBBhnAo5q8osr2H74VJ33nUGIncQ4A0nxzuAiKd5IUny053my13P3Msnx0Z7ncY2c1fceZz+YXs1/Lk2jvLaS4poySmosnkf384LqIrZbDrheWyhxvX/KXk68PoYUg5nORmcg38lgJjnaRLIhgSSDiWSDieTohDOvXe8lG0zE62NDFkw5NAf6IJ1hLyipYNCDHxNv1GOrdZzVlsnxRpJMrsf4aE9weaYfuJaNMxLdxlPchUpncwyPThlKaiMXMNY6HJyy2lz7Vf19zMbBE+We/arYUs0pq43TVhtAnf2pU516PlPf3umefS4uOuzPMDu0utMq6nQ6EmKjSYiNpk9qgk952Oy1nqD+lNf3U3G5jatGycnAUGufp0Taib1799K/f+MXo11++eWe50uWLPE8X7VqVZ3lUlNTGw3WAXr06MG8efNaUVIhOq4n97/Di4c/pH98N5IMJpIMJhIN8Z7n3n/udFtUFdGN/CwfH2MgPsZAz86mBt9viL3WQVllDaetNkoraiitcAYZpyucr09bbaiFZV7pZ5Ypraih1qF5gjt38JcUbyQxPhpzXDSJcc5/3O7n5ljDWenmOAP6qCjMhnjMhnj6xKX7XH6H5sBirxvon7aXc7qmnNN2K6ft5eRVnaTU9fx0jevRXk5pjRWbZncF9QkkR5vqPPcO+N1tYDbEOx/1rkdDHAn6OL+mXQzmdI3Wajv908xsfe4abPZaSitqOGW1ebWdjdNWZ9opqzPY9G7X067nZZU1xBn1ngA/Kd7Zdolx7rZzt2u91/We+zNcq605mrklvT4qii7mWLqYWzYjSaXN7qnXs9vAxrHTleQWlLr2Oa82sNoor7JjjjO49i0jSSZnXbvrNTG+7usG6z82mpjoqKDVf61DQ9/KvI0GPd06xdOtU+BvWCRaTwL2IJo2bRrTpgX2TqxCiMDaVX6YOb2mcE3q+ZTarXX+yuwVFFQXkWM9TGmNlbLaCkprrFTpraTqA3dzDYM+yvMTtD+qbLXOIKOibsBvqayhrLKG8io7eUVWyiprsFTVYKl0/VXZPc/Lq+zEGfWeIMMd1Jvj6gYjCbEGEuOiMcUaMMU4/+JjDCTERhMfY6JnbDKZJgOmGL3PY41tjhpnMF9j9QrqzwT7pXYreysKOW0vx2KvoMxegaW20vVYgcVeibW2CpM+FrMhjkS9M6g36+POCvDNhrg6z4tryoI2JMbh0DDonXkbDXpSE/WNnjluLh9LVY0ngCytcAbx7rZzt+ux05Wu1/az3rNU1mCzO+oEks529WpnV5op1kCCq12929nZxgYSXGmN/bLTGq25JX1T4owG4owGup19g0sfyuRw7Vdn6t97PyqrtFNaUUN+cQVlle62cbVBlbMNyipqAFwBfMP7lnu/M8U69zNnPUd76tsUa3CdNXe+NuijvMqoydj+CCcBuxCiQztttVFw3EBOgQFTTApxxq6kGA3ExeidQYlJT7zR4DpzridaH8WtL33LwO59Ql10j1ijnlhjHGnJ/t951eHQqLDZzwR5riDeOzB0B415xVYqqu2UV9mdj9XOR2uVHWu186+iuhajIapOwOcOOuJjzgSECa7X3u+bYjqTEJNGV1c7xBn1xJkNxEXriYtxPtY/W+nQHJS7g3h7JZZaV2DvHdzbKyirreBodYnnee/YNHrFdm2iZvwXqOAzKkrnGQ7Vq4vvv9zUV2N31As0zwTz3kH+0VOVzvZ0tbG1qsbTrt5tXF3jIN61n5xpu2hMMXpMsdF10hPqtXuc0XlQFxdzJvg3xRiw14bfjCT6qNYdULtV19SetT+deW2ntMKGtdpOSZGV8mpnvVuq3HXu3B/dbVBeZffsXwkxBk6UVWGMloA9kslFp0KIDm3P0VJyDx7DmHSSSpudClstFa4AtNJWS6XNjrXa+VhRXYtD09A0uGeC73csbg+ios6MefXnLGR9mqZRaav1CuC9AvyzAn3ngUDhqQqs1bWeANFd/85H51+Vq41qah3EReuJNRqIN+qJNboCe6PB9ej9PInY6BTiYwykGPX0cKXHRuuJjzGwY7eNfcZjxBr1xBj0xEZHERPtfN8YHUWs67n3GU1fBOtssb+iDYEJPN1qHQ4qqmtdgX0NFdW1lLuDe6/A3t3ehSUVnr5wZj9z72POdu6aFEtMhAaeMdF6UqP9+5WlPk3TqKqppdwVvFfZ7Cjd5e6gkUwuOhVCdGgaGjeM7c0Ll/p2U5wau4OqmlrMcR1r+sOW0ul0nvH8qUHIv9bhqBPIV9rsVNlqqbDVuh6drytraqmstnserdV2isqqqaw5E/y716+ucVBdU0tVTS02VztX2WpdaQ50OpxBvCHK+atGtJ6YaD0xhjMBfmx0FEbX8/KqmrAK2ANNHxWFOS7KtS/4/+uOaDmdTucZ5pMq9zLqEGRIjBCiQ3OgtejmQtGGqA4zY0c400dFkRAbRUJs2x042WtdQXxNrSe4dwfzziC/1hXknwn8fZ2hQwghmiIBuxCiQ9PQ0AdplhARWQz6KBL0bXuQIIQQgNzeTQjRsTk0h8yuIIQQIqxF7Bn2/+zfyZScx0JdDCFEmHOYbKTEyLzDQgghwpdO07SgZDx//nytZ8+eIZ0lZv+pIuyaIyTbFkK0DwZdFP07dQl1MYQQQoj6PFet+3WGXVGUWOA1YJArj8dUVf3Ce5lwmCVG/gkLIYQQQoj2zt+Bm08AOlVVzwduA95VFCUtYKUSQog2sHzDcvrM70PUXVH0md+H5RuWN7vOquJVXLXjKkZvGc1VO65iVfGqNiipCAaf2/Lgcvi4D/wzyvl4sB30Ez/KHBJe5Sxf0YXf/LlLi/ZH4Wdf87F/+JK3P9+jAf88ARLy/bYJLT7DrihKFHAncD2Aqqp7FEX5GZgBPB/Y4gkhRHAs37Ccu5fdTYWtAoDDJYe5e9ndAEw/f3qD66wqXsVTR56iSqsC4FjNMZ468hQAkztPboNSi0DxuS0PLoeNd0Ots59Qcdj5GqBvmPYTP8ocEvXKmWAv5s89oagc/nWi+f1R+NnXfOwfvuTtz/dowD9PgIR8v22GP2fY+wGdgVyvtBxgdEBKJIQQbWDBRws8/2TcKmwVLPhoQaPrvFz4sufL3K1Kq+LlwpeDUkYRPD635bYFZwIbt9oKZ3pr8w4WP8ocEg2U06SHP/VzPm9ufxR+9jUf+4cvefvzPdqUUO47Id9vm+HPGHb30JdSr7TTwBDvhfLz88nOzgYI6YWnQgjRkCMlR1qUDnC85niL0kX48rktKxrpD42ltyTvYPGjzCHRSHl6xZ553tT+KPzsaz72D1/y9ud7tCmh3HdCvt82ozWTD9efXqbO/ZfdF53OmTNHgnUhRNjpldKrRekAadENX6rTWLoIXz63ZXwj/aGx9JbkHSx+lDkkGinPEa+TnE3tj8LPvuZj//Alb3++R5sSyn0n5PttM/wJ2E+4HpO90pK90oUQIuw9PfVp4o1151+PN8bz9NSnG13n/u73E6uLrZMWq4vl/u73B6WMInh8bsvhT4O+3jz9+nhnemvzDhY/yhwSDZTTWguPHXA+b25/FH72NR/7hy95+/M92pRQ7jsh32+b4c+QmP1ACaBwJkgfAqwMVKGEECLY3BdELfhoAUdKjtArpRdPT326yQul3BcevVz4MsdrjpMWncb93e8PiwuSRMv43Jbui/C2LXAOGYjv5Qxsmrh4M+T9xI8yh0S9cpYbUnjsMLx7ooTePuyPws++5mP/8CVvf75HA/55AiTk+20z/LpxkqIoi4E0VVX/S1GUgcAG4BxVVY+5l8nOztZCPQ+7EEIIIYQQ7ZRnuHmr5mFXFGUD8C/gVu9gXQghhBBCCBEYfp1h98X8+fO1nj17ygwxQgghhBBCtJznDLs/Y9h94p4lRgghhBBCCOG/1kzrKIQQQgghhAgyCdiFEEIIIYQIYxKwCyGEEEIIEcaCFrDn5+eTnZ3Njz/+GKxNCCGEEEIIEfHkolMhhBBCCCHCmAyJEUIIIYQQIoxJwC6EEEIIIUQYk4BdCCGEEEKIMCYBuxBCCCGEEGFMZokRQgghhBAijMksMUIIIYQQQoQxGRIjhOiQnn12PWvXHqyTtnbtQZ59dn2j67y2Gb7Pq5v2fZ4zXQg36Sct58/+KILb13zJO9DtFsp9J9z3WwnYhRAd0pgx3bnppvc9/2zWrj3ITTe9z5gx3RtdJysN7lt15kv9+zzn66y0tiixaC+kn7ScP/ujCG5f8yXvQLdbKPedcN9vdZqmBSXj7OxsTYbECCHCmfufyz33jObVVzezYsU0xo/v2+Q67i/xGcPgnR3wymS4MKONCizaDeknLefP/iiC29d8yTvQ7RbKfScM91ud+4lcdCqE6LDGj+/LPfeM5o9//IZ77hnt0z+ZCzOcX+ZLNjofJQgTDZF+0nL+7I8iuH3Nl7wD3W6h3HfCeb8NWsDuvuj0vPPOC9YmmiUHC05SD05SD1IHbu56WLv2IK++upmFCy/l1Vc3nzUWsyHf5znPvPxmrPOx/pjH9kT6Q/DqoL31k3DoC/7sj4EWDvXQUoHua9514EvegW63UO473tv+ny01YbXfRvQY9o0bN4a6CGFB6sFJ6kHqwG3jxo2en3FXrJjGk0+OZ8WKaXXGYjbE/XPpK5Ph4Qucj95jHtsb6Q/BqYP22E9C3Rf82R+DIdT10FLB6GvuOvAl70C3Wyj3nfrbnlz7aVjttyEJ2P05gm2ro962LJvUQ3jXgb/bCud1/BGp+8SmTYV1xlqOH9+XFSumsWlTYaPb2X687pjGCzOcr7cfD1zZwrkv+LutcF7HH81tp6F+cn+f3Y32k9ZsK1Dr+COQZWtqfwznOvB3W4Fap7nvpNbUgy/fd+52i48/ATT/Pdraz+PvZ/Jlnfrb7qUraPL73d/t+LtO0C46VRTlDSC/kbf7AIdamGWkrdOW25J12nZbso60UXtYpy23JetIG7WHddpyW7KOtJEv6xxSVfUtCGLALoQQQgghhGi9iB7DLoQQQgghRHsnAbsQQgghhBBhzBDqArSWoiixwGvAIJyf5zFVVb9oZNmHgdtcL99VVfW5tillcPlaB4qi3A48CJz2Sn5cVdXvgl/KtqEoyhjgX8BT7nFfjSw3HXgI0ICvgUdUVY2I8WG+1IGiKONw9pljXskvq6r6ftALGGSKokQD9wHX4bzphBFnP1/dyPIR2RdaUg+R3B8AFEWZw5l6SAVeV1U1u5FlI7U/+FQHkd4XABRFGQjkAFeoqrqukWUiMl7w1lw9RHLMoCjKWzjHjnu7WlXV8kaWD/n3QrsP2IEnAJ2qqucripIJbFAUZbCqqnWu61UUZRJwFzDClbRVUZQcVVU/a9PSBscT+FAHLg829gXV3imKMhW4EShtZrmhwPPAUKAEWA3cC7wS7DIGm6914LK4qYOadqwH8AAwQlXVUkVRrgA+URRFUVW1wHvBSO4LtKAeXCK1PwDciTMoOeEKUnYrivKjqqp15vCL8P7gUx24RHJfAHgSsDX2ZoTHC96arAeXiI0ZVFUd58ty4fK90K6HxCiKEoXzS+gfAKqq7gF+BmY0sPhs4J+qqlapqloFLAf+u63KGiwtrINIt0lV1dsASzPL/RpYqapqkaqqDuBNIqAvuPhaB5HMAixSVbUUQFXVL4Eq4MIGlo3kvtCSeoh0M1VVPQGgqupe4BRnn12DyO4PvtZBRHP9AlkOnGxisYiMF7z5WA/CKSy+F9r7GfZ+QGcg1ystBxjdwLLuYQLey90bvKK1mZbUAcBsRVH+gPNnnRWqqi4NcvnajKqqjU0jWt8Y4H+9XucA5yiKEqeqamXgS9Z2WlAHANcpijIL5/fA/+E8q2YPTsnajqqqxcAy92tFUdzDQRr6xxTJfaEl9QAR2h8AVFXd7n6uKMoNOAOVhoZORnJ/8LUOIIL7AvAH4G7giiaWidR4wZsv9QARHDMoivIacA7OfeFZVVXXNrJoWHwvtOsz7ECa69H75//TQNdGlvVlufamJXVwHOcX9DjgJuB+17jGjqahvqADuoSkNKFRCvyA88t6MvBL4JmQlih4LgMOA9808F5H6gtN1UPE9wdFUYYpipIDvADcoqrq6QYWi+j+4GMdRGxfUBRlMrDLhxMbkRovAC2qh0iOGXYD/1BV9RJgEfCpoigjGlk2LL4X2vsZdrf6A/91Pi4XSZqtA1VVV3m9PKEoSjbwG6DBi68iXEN9obF+E3FUVf0Z59ApgHJFUZ4B3lcUJSIusHNzXZD9J+B210+ZDYn4vtBcPXSE/qCq6g5giKIoI4FViqJc28j47YjtD77UQaT2Bdfw0fnA9T6u0m4/a1NaUg+RHDOoqvqM1/NNiqJ8hnMo1D2NrBLy74X2fob9hOsx2Sst2Su9/rL1l4uEsVstqYP6jgC9A1ye9qChvqARGf3BX0eAeJyzR0QE1xCQ14EXVVX9qZHFIr4v+FgP9UVcf3BzBaSf4Qw86ov4/gDN1kF9kdIXbgP+T1XVEh+WjdR4AVpWD/VFcszQ1GcLi++F9h6w78d5xa7ilTYE2NTAspt8XK698bkOFEV5tF5SGlAYvKKFrYb6wq72Pka1JRRF+Y3rrKtbGs7ZAopDVKRgeB7YqKrqvxVFiVEUpVcDy3SEvtBsPURyf1AUpbNrzLY3K2BqYPGI7A8tqYMI7guXAFcrirJOUZR1QDrwV0VRPmlg2UiNF6AF9RDJMUMLP1tYfC+064Dd9dPu34E7wDOn6AhguaIogxVFWa0oit61+GvArYqixLq+jG5zpbVrLayDyYqiXOZaLh7nzz/Lzs41siiK0kVRlG8VRUlxJb0BXOn6JxYF3E4E9IWmNFAHo3COScTVP+bgnBWhNlRlDCRFUebjHPL3lqIoCUB/4I6O1hdaUA+R3B/MwELXdx6KoqTinI98dQfqDy2pg4jsC6qqzlZV9SJVVce5pvM7hnPKwikdJV6AFtdDJMcMDymK0hVAUZS+wBRcny1cvxfadcDu8gSgUxRlA86rum9VVfUYkITzRkLRAKqqfo5z6sP1rr83I2hO1SfwoQ5wnmlbpCjKWpwT/68F/tz2xQ0ORVHOdZ0xGAH8VlGUD11vxeGsh3gAVVV3Ao/gvJjmR2ArEBFXvvtaBzgP8m5x9YUNOH/ae6BtSxscivNeBItxBhoW198u19sdqS/4XA9EcH/AGZB8ijM4/Rr4Cue0bEvpOP3B5zogsvsCiqKMrXdmeREdK14AfKsHIjtm+AvwkWt/+Ddwv6qqX7veC8vvBZ2mReR1FUIIIYQQQkSESDjDLoQQQgghRMSSgF0IIYQQQogwJgG7EEIIIYQQYUwCdiGEEEIIIcKYBOxCCCGEEEKEMQnYhRBCCCGECGOGUBdACCFEYLlu9rIXUFRVrQh1eYQQQrSOnGEXQogI4LrV+O0AqqpWAcMkWBdCiMggAbsQQkQgVVVPh7oMQgghAkPudCqEEO2coih/Bu7BeQv6Y8BwnLfXngR8j/OW2pcB9wPXAP2BXwNjgOtxnry5VlXVk678zgVeBDTADtynqmpuG34kIYQQXuQMuxBCtHOqqv4O2AosVlV1nKqqnXAG7qiqalNVdZxr0RhVVScBrwDLgfWqql4EnMQZwKMoShLwOfCEqqqXAS8AnyiKIv8vhBAiROQLWAghOo4vXY87gXhVVTe4Xm8H+rmeXw2Uq6q6BkBV1c+AdOC8tiyoEEKIM2SWGCGE6Dgsrke713P3a6PreU8gRVGUdV7vnwQ6B710QgghGiQBuxBCCG95QL7XMBoURUkEqkJWIiGE6OBkSIwQQkQGCxCvKMpARVGea0U+/wG6KIoyBkBRFBOwFkgKQBmFEEL4QWaJEUKICKAoylRgMVAK9MA5hCUX58WkfwauAH4EZgHvAYOA/wd8BiwBYoFnVVV9wTVLzPOAzvX3rKqq/2nTDySEEMJDAnYhhBBCCCHCmAyJEUIIIYQQIoxJwC6EEEIIIUQYk4BdCCGEEEKIMCYBuxBCCCGEEGFMAnYhhBBCCCHC2P8Hhz/jbz23TSkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 931.988x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.rcParams[\"axes.spines.right\"] = False\n",
    "plt.rcParams[\"axes.spines.top\"] = False\n",
    "plt.rcParams['font.size'] = 12\n",
    "plt.rcParams['xtick.labelsize'] = 12\n",
    "plt.rcParams['ytick.labelsize'] = 12\n",
    "plt.figure(figsize=(8 * fig_aspect,2))\n",
    "a, b = plotHawkes_plot(s_plot, lambda t: mu0, alpha, w, T, 10000.0, label= r'$\\lambda_m(t)$', legend= 'accepted in counterfactual')\n",
    "a, b = plotHawkes_plot(np.array(c), lambda t: new_mu0, new_alpha, w, T, 10000.0, label= r\"$\\lambda_{m'}(t)$\", legend= 'accepted in counterfactual')\n",
    "plt.plot(mu_events, [-4 for i in range(len(mu_events))], 'x', color = 'navy',  label = r'$t \\sim \\mu_m(t)$')\n",
    "plt.plot(s, [-4 for i in range(len(s))], 'x', color = 'dodgerblue', label = r\"$t \\sim g_{m}(t)$\")\n",
    "plt.plot(all_c['mu'], [0 for i in range(len(all_c['mu']))], 'o', color = 'darkgreen', label = r\"$t \\sim \\mu_{m'}(t)$\")\n",
    "plt.plot(all_c['exp'], [0 for i in range(len(all_c['exp']))], 'o', color = 'limegreen', label = r\"$t \\sim g_{m'}(t)$ (existing $g_{m'}$)\")\n",
    "plt.plot(all_c['new'], [0 for i in range(len(all_c['new']))], 'o', color = 'orange', label = r\"$t \\sim g_{m'}(t)$ (new $g_{m'}$)\")\n",
    "plt.legend(loc='upper left', fontsize = 12, ncol=3)\n",
    "plt.xticks(np.arange(0,T +0.5,0.5))\n",
    "plt.yticks(np.arange(0,20,5))\n",
    "plt.xlabel('time')\n",
    "# plt.savefig('figs/hawkes_figs/example5box.pdf', bbox_inches = 'tight')"
   ]
  }
 ],
 "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.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}