{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "08bd379a", "metadata": {}, "outputs": [], "source": [ "import torch\n", "from models import TimeAwareGPT2\n", "from utils import load_model\n", "from tqdm import tqdm\n", "import pandas as pd\n", "import numpy as np\n", "import textwrap\n", "import matplotlib.pyplot as plt\n", "\n", "plt.rcParams['figure.facecolor'] = 'white'\n", "plt.rcParams.update({'axes.grid': True,\n", " 'grid.linestyle': ':',\n", " 'axes.spines.bottom': False,\n", " 'axes.spines.left': False,\n", " 'axes.spines.right': False,\n", " 'axes.spines.top': False})\n", "plt.rcParams['figure.dpi'] = 72\n", "plt.rcParams['pdf.fonttype'] = 42\n", "\n", "#Green\n", "light_male = '#BAEBE3'\n", "normal_male = '#0FB8A1'\n", "dark_male = '#00574A'\n", "\n", "\n", "#Purple\n", "light_female = '#DEC7FF'\n", "normal_female = '#8520F1'\n", "dark_female = '#7A00BF'\n", "\n", " \n", "delphi_labels = pd.read_csv('delphi_labels_chapters_colours_icd.csv')" ] }, { "cell_type": "code", "execution_count": 2, "id": "1d8375ab", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
nameICD-10 Chapter (short)
46A41 Other septicaemiaI. Infectious Diseases
95B01 Varicella [chickenpox]I. Infectious Diseases
1168C25 Malignant neoplasm of pancreasII. Neoplasms
1188C50 Malignant neoplasm of breastII. Neoplasms
374G30 Alzheimer's diseaseVI. Nervous System Diseases
214E10 Insulin-dependent diabetes mellitusIV. Metabolic Diseases
305F32 Depressive episodeV. Mental Disorders
505I21 Acute myocardial infarctionIX. Circulatory Diseases
603J45 AsthmaX. Respiratory Diseases
1269DeathDeath
\n", "
" ], "text/plain": [ " name ICD-10 Chapter (short)\n", "46 A41 Other septicaemia I. Infectious Diseases\n", "95 B01 Varicella [chickenpox] I. Infectious Diseases\n", "1168 C25 Malignant neoplasm of pancreas II. Neoplasms\n", "1188 C50 Malignant neoplasm of breast II. Neoplasms\n", "374 G30 Alzheimer's disease VI. Nervous System Diseases\n", "214 E10 Insulin-dependent diabetes mellitus IV. Metabolic Diseases\n", "305 F32 Depressive episode V. Mental Disorders\n", "505 I21 Acute myocardial infarction IX. Circulatory Diseases\n", "603 J45 Asthma X. Respiratory Diseases\n", "1269 Death Death" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Delphi is capable of predicting the disease risk for 1,256 diseases from ICD-10 plus death. \n", "# For illustrative purposes, some of the plots will focus on a subset of 10 selected diseases - the same subset in used in the Delphi paper. \n", "\n", "diseases_of_interest = [46, 95, 1168, 1188, 374, 214, 305, 505, 603, 1269]\n", "delphi_labels.iloc[diseases_of_interest][['name', 'ICD-10 Chapter (short)']]" ] }, { "cell_type": "code", "execution_count": 3, "id": "f8d86352", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Using device: cpu\n", "Model config: {'n_layer': 16, 'n_embd': 256, 'n_head': 16, 'max_epoch': 200, 'batch_size': 128, 'lr_initial': 0.0006, 'lr_final': 6e-05, 'weight_decay': 0.2, 'warmup_epochs': 10, 'early_stopping_patience': 10, 'pdrop': 0.0, 'token_pdrop': 0.0, 'betas': [0.9, 0.99]}\n", "Model loaded from best_model_n_embd_256_n_layer_16_n_head_16.pt with 13.29M parameters.\n" ] } ], "source": [ "seed = 1337\n", "\n", "torch.manual_seed(seed)\n", "torch.cuda.manual_seed(seed)\n", "device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\n", "print(f'Using device: {device}')\n", "model = load_model('config_n_embd_256_n_layer_16_n_head_16.json',\n", " 'best_model_n_embd_256_n_layer_16_n_head_16.pt', \n", " 1270)\n", "model.eval()\n", "model = model.to(device)" ] }, { "cell_type": "code", "execution_count": 4, "id": "a8658fb6", "metadata": {}, "outputs": [], "source": [ "# Let's try to use the loaded model to extrapolate a partial health trajectory.\n", "\n", "example_health_trajectory = [\n", " ('Male', 0),\n", " ('B01 Varicella [chickenpox]',2),\n", " ('L20 Atopic dermatitis',3),\n", " ('No event', 5),\n", " ('No event', 10),\n", " ('No event', 15),\n", " ('No event', 20),\n", " ('G43 Migraine', 20),\n", " ('E73 Lactose intolerance',21),\n", " ('B27 Infectious mononucleosis',22),\n", " ('No event', 25),\n", " ('J11 Influenza, virus not identified',28),\n", " ('No event', 30),\n", " ('No event', 35),\n", " ('No event', 40),\n", " ('Smoking low', 41),\n", " ('BMI mid', 41),\n", " ('Alcohol low', 41),\n", " ('No event', 42),\n", "]\n", "example_health_trajectory = [(a, b * 365.25) for a,b in example_health_trajectory] " ] }, { "cell_type": "code", "execution_count": 5, "id": "d7f2e4a1", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Input trajectory:\n", "0.0: Male\n", "2.0: B01 Varicella [chickenpox]\n", "3.0: L20 Atopic dermatitis\n", "5.0: No event\n", "10.0: No event\n", "15.0: No event\n", "20.0: No event\n", "20.0: G43 Migraine\n", "21.0: E73 Lactose intolerance\n", "22.0: B27 Infectious mononucleosis\n", "25.0: No event\n", "28.0: J11 Influenza, virus not identified\n", "30.0: No event\n", "35.0: No event\n", "40.0: No event\n", "41.0: Smoking low\n", "41.0: BMI mid\n", "41.0: Alcohol low\n", "42.0: No event\n", "=====================\n", "Generated trajectory:\n", "43.9: B95 Streptococcus and staphylococcus as the cause of diseases classified to other chapters\n", "44.6: J30 Vasomotor and allergic rhinitis\n", "44.8: L03 Cellulitis\n", "46.6: L30 Other dermatitis\n", "46.8: K56 Paralytic ileus and intestinal obstruction without hernia\n", "47.8: K76 Other diseases of liver\n", "48.8: N20 Calculus of kidney and ureter\n", "48.8: L24 Irritant contact dermatitis\n", "49.5: I10 Essential primary hypertension\n", "49.5: K59 Other functional intestinal disorders\n", "50.6: B96 Other bacterial agents as the cause of diseases classified to other chapters\n", "51.4: E14 Unspecified diabetes mellitus\n", "51.6: E55 Vitamin d deficiency\n", "51.7: A41 Other septicaemia\n", "52.0: F41 Other anxiety disorders\n", "52.1: F32 Depressive episode\n", "52.3: M19 Other arthrosis\n", "53.5: G96 Other disorders of central nervous system\n", "53.7: E87 Other disorders of fluid, electrolyte and acid-base balance\n", "53.8: J18 Pneumonia, organism unspecified\n", "53.8: E66 Obesity\n", "53.8: E11 Non-insulin-dependent diabetes mellitus\n", "54.4: Death\n" ] } ], "source": [ "max_new_tokens = 100\n", "\n", "name_to_token_id = {row[1]['name']: row[1]['index'] for row in delphi_labels.iterrows()}\n", "\n", "events = [name_to_token_id[event[0]] for event in example_health_trajectory]\n", "events = torch.tensor(events, device=device).unsqueeze(0)\n", "ages = [event[1] for event in example_health_trajectory]\n", "ages = torch.tensor(ages, device=device).unsqueeze(0)\n", "\n", "res = []\n", "with torch.no_grad():\n", " y,b,_ = model.generate(events, ages, max_new_tokens, termination_tokens=[1269])\n", " # Convert model outputs to readable format\n", " events_data = zip(y.cpu().numpy().flatten(), b.cpu().numpy().flatten()/365.)\n", " \n", " print('Input trajectory:')\n", " for i, (event_id, age_years) in enumerate(events_data):\n", " if i == len(example_health_trajectory):\n", " print('=====================')\n", " print('Generated trajectory:')\n", " event_name = delphi_labels.loc[event_id, 'name']\n", " print(f\"{age_years:2.1f}: {event_name}\")" ] }, { "cell_type": "code", "execution_count": 6, "id": "0b937d75", "metadata": {}, "outputs": [], "source": [ "from utils import PatientEventDataset, get_batch\n", "train_data_path = 'ukb_real_train.bin'\n", "val_data_path = 'ukb_real_val.bin'\n", "train_data_arr = np.memmap(train_data_path, dtype=np.uint32, mode='r').reshape(-1, 3)\n", "val_data_arr = np.memmap(val_data_path, dtype=np.uint32, mode='r').reshape(-1, 3)\n", "block_length = 128\n", "train_dataset = PatientEventDataset(train_data_arr, block_length)\n", "val_dataset = PatientEventDataset(val_data_arr, block_length)\n", "dataset_subset_size = 2048" ] }, { "cell_type": "code", "execution_count": 7, "id": "6a3bb2dc", "metadata": {}, "outputs": [], "source": [ "input_events, input_times, target_events, target_times = get_batch(val_dataset, slice(256))\n", "with torch.no_grad():\n", " p = model(input_events.to(device), input_times.to(device))\n", " p = p.cpu().detach().numpy().squeeze()\n", "t = (target_times-input_times).cpu().numpy().squeeze()\n", "target_events = target_events.cpu().numpy().squeeze()\n", "ignored_token_ids = [0, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]\n", "mask = ~np.isin(target_events, ignored_token_ids)\n", "p = p[mask]\n", "t = t[mask]" ] }, { "cell_type": "code", "execution_count": 8, "id": "1c32bd6e", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ4AAAELCAYAAAAofGgWAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAALEwAACxMBAJqcGAAAUUhJREFUeJztXQe0FEXWvgRBsiigmMCAJAmKIAgICoqrK6irmAVFZMW0KJhFwfgvSlJ0XVZXXTMuimFRJCOCgiggCKIgBiTnR3jAm/98BXesV1R1mp6e6Zn6znnnTeiurr7TdevmWyKRSCTIwsLCwgdK+jnYwsLCArCMw8LCwjcs47CwsPANyzgsLCx8wzIOCwuL3GMc33//faanYGFhETfGAezatSvTU8gaOlha/AFLj8zRokQc4jh27NhBBx54IOU7QAfA0iI1emzbto0KCgqoQoUKVL58+f3eq8eVKlWKtmzZIhbmwQcfLP74+w0bNojXVatWFf9xPOaF/yVKlKA9e/ZQUVERpRu7d+8W/0uXLp3SOCtXrqRp06bRJZdcQocffnjyXlWkdpWIEAPeFgksHVKnBxb7smXLxLlY2DVr1qTff/89+f6YY45JMhM+TsaKFSuSzEr+HgwE58f5NwLTuPfee2nAgAHJezUxjlioKhYWYQGSAC9u/N+8eXOx9/hePU7F5s2btd/HmWms37yNbr/zHrrgr3fTnK1V6J+z1zseHwuJw8IiLEAdYckA/ytXrpxkAniP79XjVFSuXFlIHHGXML78dRu9MW8Trdy6m7YUFlFR+7vpn9+XhfxEB5YuQTecopc2YmPjKCwspDJlyhT7DPrmr7/+mtRz8wH8U+GB9Qs86EceeSQdcMABlCvAcwGoz0bYNg7YKdauXZv8/Nhjj02eB/Ee/1WULFkyEtsGA3MEYI/xgu/W7KT7xq+iGhVKUaNDD6RDK5be+1dh7//KZUvSEUccEW8bh26hgGlUqlSJateuHWghxRH8IOKh9Mtw1q1bJ2gGHT5XEPR3x6KXGYT6Xv18zZo1xT4vKChIHq9jGkCVKlWEwRJGVd28eRPAM4w/tp34BRgFMw2v+H3LLnp48mqqWLiemhZ8TX/tfK32OCejcyxsHDrODUnjkEMOyRumkQpAI9Aq16QzPBded3UscDAA00J3OobVFkBWZwr22UN0gLFUxzQAWcgHA9q+fTsFhV+msXnnHnpo4moq3LSG9kwYRud16mg81un+YiFxmB4OyzS8Ixdp5YdpyJ4U9px4PQb/8V5VZ0p5VAvc5sYqVyrwYnHYtSdBj05ZQys376AKnz1LTzw6kI466ijj8U73FwvGYWHhFTp7hepJkVUNhtsx8lj83u9ur4NJKgkbRYkEDZ2xlhas3kn92h5GrS591tU25HR/lnFkGfBjhbGT5SNMUoPqSWFVQ4bbMbqxS2l+J9grZGZQrVo18ZvifBj0ozSYynht7kaa8tM2uqbpQdSu9v73r4OOTrGycfg1BkaJCy64gJo1a0YNGzakf/7zn/SPf/yD+vXrl/z+pZdeoptvvlm8fvXVV6lFixbUtGlT6tWrV5KjV6xYke644w5q0qQJzZgxgwYOHEjNmzenE088kW644YbkTjhr1ixxLv5wDXwPYBy8xzmNGzem559/PlIaeLEfpOu5kJ8NndQgqxqHHnqoVk3xcow69oYNG0TgmG4ceF0wDv7Ddbtx40bauXNnKEwDjEknKYCZmdTRcT9spdc+/5EOWTiKLmlY2dN14H3T0SlWEocb4xjwwQJauGJzqNdscHhlevD8hq7Hvfjii8JlBQMXFu6ECROodevWNGjQIPH9W2+9Rffddx9999134vX06dPFj9K7d2967bXX6JprrhEP5amnnkpPPfXU3ms3aED9+/cXr6+++mr68MMP6fzzz6cePXoIpoDx77777uQcXnjhBWHFB2PBA4rvzz777Eg8KF7sB1E9F05Sg8lzIsPpGDWuY/v27fvZFfiasloD+oQZ8YBNwo9NZP6qHTTs00W089Oh9NDQx0KzdcWCcWRzqMnw4cPp3XffFa9/+eUXsYiw08ycOZPq1KlDixYtEgt5xIgR9NVXXwnmwg9ejRo1xGuIvH/5y1+SY06aNIn+/ve/i4du/fr1Qppp27atEIFbtWoljrniiisEQwHGjRtH8+bNo3feeUe837RpEy1ZsiQSxuHFfhDVc2EyYoYBHhuSxsaNG/fzUCFXpVy5cslcFQ5jDwMYF65RXNuPVLe2YDc9MXUVbZnwHI188lE6ptbRvnNfYs043IxQXiSDdGDy5Mk0fvx4oV7gwWrfvr14aC677DJ6++23qV69enThhRcmd6pu3brR448/vt84eChYX8b5kEZmz54tLN4PPfSQqxsVYz/99NPUqVMnihpe7AdRPhdeJIugRlb8V0PNK1asmNwAdLktYUlWTq5RQL0uPCiPT1tDhXtK0CvPP0O1Dynn+5qO3/sazaIYsLNjp8EDBckCUgYAZjFmzBh64403BBMBOnToICSC1atXi/eQJJYvX77fmMwkYFTbunVrUoo46KCDhH77xRdfiPdvvvlm8hwwjOeeey6ZVo0aJm4PWljwYj8IA6nEYQS5FpjAqlWrxH8ej12nJaSYDjANHUMJExjbr8t25FfrafHaQrqt1SG+mYac7RtriSNbcc455whjaP369alu3brUsmXLJNHx2cKFC4UxlO0WjzzyiLA9wEgGOwfUl1q1ahUbEwyiZ8+ewvB52GGHJVUbYOTIkcKoit2gXbt2wq4BXH/99fTTTz/RySefLB7e6tWr03vvvRcZHcLa5dMVh+EXUAlUQ6gqUVStWjW5abjltqjA74ZNJ10YNfNHeu2fL1OPm/tQm1rBJEA3z17Gc1VgNBw7diwtXbpU2At0IhIMfmXLlt3vPCzOfAKyMiEag0ZPPPGE0KOHDRvm+fy40gxSBHZ/BqQbMEc8FwCeDdMxQfDbb78lmQUABgFPhtv4y5cv9xSXgZoZbjYEv+DxZi/5jW7tew+1urofDe7alEqVDGYM5XycyFUVxN5jB4T+LhOpT58+wtB32223ifd4kPHDcD2DuLljo8RHH30kaAppBMVW7r//fsoHmEK+ZXes6RgnmFQbVUzHe7fxt+2L0/CCsJkGzwnh5HcO+D867sK/0cALGgdmGpA23KS1tKkqcFHCNQl9nzFnzhyht+Ohv/HGG4X7EKL4tddem6yNwOI3A/o6iq0wsUEg3BiOz1QwTabQtWtX8SczWD80AM3gzeEFh/eygVFeiPLDzTQHcLwspHLFKfVzHI/z1HwSt2urn2MMqHXI7MWzAw8DPpPdoTiei/JgAcMWhHPU2AmeE87DWPCC4TWuC5URY/NYqH4FOxTe4z8+wzON9+XLlxfjQ93AOHwMz5dpLUP+zeTvTJ97GUv+vHD3bhr8+Qaqev6d1L9jDapQOlFszfi5NuxrYIJOmdRpYxyQNNTsOhgPzzrrLPG6Y8eOwhsB4oOhQEe/6qqrtGOBAKmWRLOIDljUWLRYYGGVOcRYauATLwz+HAsff14ib+X4CvzHe5zLwPzBFLC4YLBm2wk8XeXKlRPHQzXRMW4wojDC0b0Cc3t9/hb6ZlUh3dS8KtU5xF+ZgSCIdDXC/w3dCYBksWDBAsFA8GfCCSecIB5E+UeVI+VyMXkr7LR6gB98lY4mmHYbt8+xoDiGAbtwmJ4W07WDMCdID9i0mC54HsGA2A2LP91mxZLVhg0b9otclRHVRvfTryvojseepp2n/ZXOOrY8/alucYk9CFB7BBGvThJHpMYD/DhQRwD8hwchaHsEPCyoMZHNwWHZAq7HEUWRY1PYd7aEqTu5kWU3rC52hlW2ZcuWCWkkCFQjfyr4+vtf6Kre/aigfhe6olElurG5+3ryYyDOGokDUY8ImYaejsCp7t27ezpPJ1VA50VhGrXISi4jjApg2RAQ5td1aqrQFYQO6ljyeGosBmwleM8BepA2EIeTSGGzwvXYG5QKPlteQHcPGEKHnncb3d+5ETWqHm5lNzdDb9oYBy78pz/9iebOnSsClB577DGRj4EfAV4VJGpxjIMToKrobgJiVC5Vs/ICpkM2l//zEvbtJ0zdiclgMeN7uEC9hJjrxuL54HyV6UHd4FgNqF2yO1YGfg+vHpVUmcbuogT9e84GGrNoC7W4qi/dc3oNqlHxgNA9NWCaWR3H4QVuFt58QRwYhxf4kTic4jNgo4CBEnYGL5KLOhaYAuxuKiOB/YLjOPhzSBpB1ZOwgNyT/5u2hr5bW0jn161E151clQ4oVSJQzVEngGEiXcIJWR8gARtHlBbqbAbokAu08BOm7hQ/gYXMBmMne4ppLD5PPl/13vDnMBaqKOOzSHJQ7NhdRK/P20jXvTqXpr/0KPVtfQj1an5wkmnwPMOSATj3xgmx8HHGQCiKBLlEB7YvsJHUpGroVB+2U7Ck4TXoSx0LYImDJRjs2F4T98qWLRtK2T8T9hQlaMLSrfTq3E20etVK2j1+GD396ABqdkxFSie4cLKpwnlsVBWdOzYfwUVt40oLXVsCv/kl8jmQvrA7QrR2snHojKtqe0c5VBzBX7DFyeeoag7A95BKxXETvlqxnV6cs4GWb9xF9aqVpRKf/4tu7dnNWCM0rBaQXsPOYyFxWMQfOiYRpJaHeg5UFaecFJNBVP5MdZHCvnHccccVm4tO6tilGEQxTqru5XXbdtPQGevo69930GEVS9PdbatR66PLU4lzok8vcPo9SsfBxpFv3hMT4hw9q2MSQWp5yOdAVUHSn9/r8mv+D2OzHLcByQ4MwI2J7VIYR6pMAx6Tx6eupZ82FlLPZlXpT3UqUpnS3syQ6QiEjH3NUYv4Q2fkDFLLQz4H+SVu5+iuq34GiUV1P8KLwowArljEDKUbr3yzkRat3Um3tDyETq26jfre/rf9cm2iAtQ1J9rGwsaB3SCKqMdsB++KcaWFKZAr6Hk6ejjZM5w+03Wn50ZWcvvHdPZxHTh5jZAy/lJrV7JrvFPfEz82DtjF/DR+citLEAvZNwa8LRLEnQ663iRucDKgqvQwHasrNKRrA4nj4VFghsQZ2yZUq1ZN3EsqndiA1QXIbF1Hx1Y9gHqecjD9Y8TTvpiGG3BvKArlp7Shm3ST9aqKKVfFIn4wleQLI/fF77HqvOA1AcNQIzt18RsMMJVUmQbsGgjq2pNI0F1tq1OZUiXo1ltvLcY0/AT8IahNlTpwf5zh6xVu3qFYSBwWuYEgXhQ/BlR5l/RTzMe0E2MRYqdGoBeYBOwNMrMrDCGG4+WvN4raoH1PO4iOqKxnEKZwdnhxVEYHlUSuXsbwq27FXuJArortbLYXoEOcaRGkSpeTAVWmBxa03PEdRX28qEJORYY5XgaBULVr1xZMJEx88cs2eve7zXR6tW305t/v8pzvwpCZBufWmCQgv0zOTeLIesYB5FPNDSfEvf5I0IroOA6GOvV4mR7qLutVhVCZGewWDNQVkSUM7tJWtWrVlH+H1Vt305AZ6+iIEhtp4VuDRNMuryqJqQYIFrtO2ggCt2LKsVBV8q1EYC7TIcyK6KnSgz0rkE6w6MBEZLuISZ0qU6ZM8hyOFsV/WeJxwtbCInpin12j3sYv6OKB/gyhbvcN1zIYZyoZs269fErHwTgalnU57sgFxpEueqiZrm59QRCbwQtd9daYbCpuIfJblNB1HRav3Ul/n7aG1m7bI4yhpx19A4UBWd3yUmndC5yC4LKecVjkNoLGdrhJDk5jyUwD4FKAPA8Oh4cUIbuO1X4rK1asEIFSACfdmVCUSNB7320WxtCqZfbQ3zsdQXWrhVMNjIt3h505nVLIOZokow0h6h5A9GFuiz4oUYBrjlrkHsJoooRng2uccqVzpzgR1YjKUOtvgIHIc8O4kGhUcX7p0qWujZg27dhDQz5fR7NXbKcmFbfSig+GUK3zn6FsLrfgZrx2ZRzokD5kyBBq1qxZxiz6cfYkhIlco0OqDatBDzk+gaUAwMSIdLEdsAmweC9LH/K4cMeamEPCgWl8u2oHDZq+VjCPy4/ZRRNeeIoGDhgQevRvWEZz0AI0c5PaSnopMIwSgEhfRvgt/0WJOHsSwkTcvSphuGdl8Dk6mpgCwNTjoWrAYyN/Bsahui/90j6RSNA7CzbRveNXUdlSJeipc2pSyZ9nhRoR6geyt8hNgvOiNrrmqtx9991CDLrooouKpR+jo1hUwI8YVbWlbAY/zLlEiyA2Dj4H7kucAxWabRKy2mJSfeRrAnwunnOuYZoKihIJGjl7A32weAu1rVVeJK2VPyC9kQ9OpQMR3wFJwo+r1q0FpKuqwt3RZ8+enfwMP8rEiRMpKuRCubwwkIt08OueVQv5YPdGgJY8BtQKhIp7YRqyHQOSdKpMY09RgobNXEcTlxZQl3qV6PL6ZUNlGqZiQU77PxiryjTASLgYMzJ/VQnLTW10ZRyTJk2iTMK6Yy1kqLYHMAIucSc3g+LCwnL5O9UFi74+8lhBs2BLliwpXMOFe/bmnXzx63a6uslBIiL01lv6isbgblXDvSKMYkEAGCS7rXUlEFMOOUdCEgyksHMACxcupBdeeCHVeVtYhGIXceqLAiaBRYY/NBhSXbA8hhvUfsYqRA5LYRE9OHEVffnrdrqx+cHUrvp2EQ0Km0ZYTANQmQZsF0Hd2KCBTtrwIt262jjAMNAU+tFHHxU9UiD2nHTSSTR//nyKCojJD7MDVlzBuQn5TgvZxoF8EqYHPlfDBOTAMBXcjpRtHKbIT4yB700SyaYde+jBiatp2YZC6nNaNWp/TAUaM2aMsANGJS2HXXPUzcbhKnGAWOi8xsEtmFjUbsEgvVJzEU69SuMOPy0hOXcF0ocbPTj2yFThiseCe9QkfYBhm5jGyq276K5xK+nnTbvo/nbVBdMAunTpEqmKHbXHzZU94cdB31GeFDrOu4luYSPuBWzCQq7SIWggmEoP1f0K75Nq7IT0wN3Z1HNN9NUxs1179kaCvjl/E5UqSTSwQw06quzO5D3EHSkbRwcPHkydO3emH3/8kVq3bi12hXfeeYeiRC56E4IgV+kQNBBMpgfbMmSoujsYxhFHHOE6lhvmrtxOz325nn7dvJta16pIPU6uQkVb1lKfPvfSk08+6diPJAyobRnSsamkXMgHxJ4yZQotXrxYTK5u3br0zTffhDlHizxHkGrnXovxMNiDoIsbgbfFi0dl/bbd9MKcDTTlp22idcHgLsfTacdUoW+//TZZIzTdTAPwaipgbw8D9w8pDOfjnp0yYN2yY12Nowg1f//995OceurUqXTTTTdFahy1DZlyoyFT2IFgTA+4Xk0NoWX1BFDzT7CzovaE00JBbMZH32+hV+duFC7XSxpWoYsbVqZqVasIVejTTz8VvVajtGnARoN4FXY764yj8LjIDBHncBIg4MRscWxKndxmzZpFvXv3pg8++IDmzJlD99xzD3344YeREQlxHEcffXRsK3uHibhXOU8XPbCrygwB8RlyVW+2aei6sblh6YZCenrmOlqyrpBOrnmg6Nl6xL4Sf34rh4cJ3BOuz54gHeOAdIGqYm4eJTATmcGA4bhVO3NVVZo3b07Dhw+ns88+Wzyw48ePdyybng7kWnJXUFg66OmB51LXX5aDwRA1yQvNlMmKRSbbRHbuLhKGz9ELN1PFMiWpX5tqdHqt8knD59q1a432kigASUcu1KMzyDqVCwRNMH9mqH6fMyPjOP/884tNBj8EvCkIBgOgvkSVVu+3FqNF/tXlUFsvqF4SLBTEc+j6pHDJQN6956/aIaSMFVt2U8djK9B1zapS5bJ/LKaVK1cKm8ZTTz3lWjDIK/z2nQXTSCU8HudyoR6VUaTEOPr27UvZAhDJT4n4XAXvMNlOC79MIKg7VqaHOgbsFyrwnWrLQAtJ5LVA5SgqVZZGTF9Bn/ywlQ6tWJoe6VCDmtYsbk9ipjFgwAA68sgjPbdgcIPKNHBPUMHcmAkkJRhB2dahgu0aavIeaIdgOXyvBr7hPaS4QO7Ydu3aJV9DL4StA2jRooVIsbfIPYRVjcsvE0i1LoduDN1ujPlsKSpDnyzZSuu376YN2/fQ9sRmWrFhCW3YsUe8L0oQXVS/Ml3RpAodqOnbisU2YF9qvI5p6FylQaCTsmG/UOuIOqkjCHWX83h0NDE1nEo5juPtt9+mfv36Ufv27cUPcsstt9CgQYPo4osvdjvVIs+qcYXZO8UvE1PHUBcYduYft5ejfi/Pp60793530IGl6OAKu6lKmZJ0ZJUDqOqBpUQa/PGHmEP6TzvtNMd5hME0TFDvCVKJkxovqyMmQLLQSStuLnFXxoEcFUgbLGXAkNKxY0fLOHIMYez6QWMyuG2CKd3dCxNTx4BKwh4PjPPxskIaOnkJ1a9ZmR7vUo+KNq0ij43gRQNqqDPZ5s3apTANqF0iS1eSQjid3pSHA3sP0xaMCXk/ulYUKlxJh4nIqgmMS1FX2w4rcSfuAB3SRYtUq3GF2TvFaytHlR48BgC9HgbPMuXK08h522nwpOXUqeFhNOqvrahmeSrGNJwKI8GmgWJWYfUrSRcgfYC5we6iVjODeuUUp8JRt2A4JluJCten8JxzzqFOnTrR5ZdfLt6/9dZbyRR7i9yBumOn0vsk1d4pzMS+X7uDhs9cRyfVKqTBl1X3rXJt2lFET32xmeb8solu7VCH/tahDpUsWQIXKOaWNdkJwDSQGv/II49oja1eVQYd3Aoc+wWkBbifEcMCycPJ4yLXWFXhVdp0DQADRo8eTZ999pl43bZtW7rwwgspSoBbZpuYmAnkSwDY7j1FNPzTRTRi6k+iCfP2XUU07LKm1KXpEY70ANOAWoFdEynuD09eTZt2FtGTlzSl85vsbWPAQH0OWYpgt6qcgo/UCkjYCEBMdZGXUeJEwgC8KZC4MK6ftHoEf4F2JvXFLaVeXNvtInfddZeoN4pkN/yBaeCzKJGrWaFB6JDrtFi+roC6Pj+Dhk9eRuc1qknT7+5AJx99ED3w3rf0+6btRnqwpAGmMW/lDur3yUraU0T0n+4nC6ahpu3LYfucx4LFLdO3adOmwnuCz1Jl1tU8Fgv2A3hM/I7L5QTUcosykA3vBlfGgTh8FWPHjvU6TwsLT8DifPPLn+lPw6bRktVbhYQx/PKT6OAKZWhw16a0uyhB/UbNoyL4SzVgu8iP6wvp4SmrqUaF0vRa9ybU4vjDkkwFYQX4jwQviPVqg2qoSDhG1zc11eJJa9asCTX+BnNFWDgYAJiB1/nJcSEmW6UXO4eRcTz33HPUqFEjkRXbuHHj5B/0YPyPCshVscjtIj3rtu6kHi99SXePnk+NDq9En/ztdDrrhKrJc2pXq0D3n9eAPvthLb0y4yfjQlq5ZZco31fxgJI08MwadEj5vWK7W48UXkzYadF8DAtK3Y3dmjC7ATaQMCOgcU9MTzCP2rVru2bmQrJCVCjT1XQ8bCRuMCpEV1xxhTCCIqntiSee0AaVRAWboxFfOniJD3nsowU0dcla6nFyVbqgQWU6sGg7LVtWvM3B5S2OovHfraLHxy6iNnWq0fE1KhWjx7aiUvTQlHW0u4jo8bMOpeoVD0jGg8i2BYwH7wMzE/YgYYPq3r07PfDAA8n0+2xCCY0xFVITS0ughRNj4ghZzt8BIKngDy5ZmUZeJCNPxtFMAwTJ9jDrKMAPRpxooWaksi1B9twsWb6CFv+8kmpXLZN8yGVxGa5duFlXb9lBnYZMpSOrlqfRvU8jKtorKRQWlaDLnv+cvl+1lR7pWIPqVz9QLAjYJdSsWc6UhbrCbRRw3Lhx44S47yXru6RS58IvgnhUnHJZcE94JsAUdJsLh+DLafgM2EhAC/V+3DJkY1HA0nZp/4MOcaOF2jkNXgu2NfCufkT1g+iYg8sWkwh0MSU1Kh1Ij1/UiOb/tomenrBE0GJH4W668bU5tOD3LXRn22qCaQBYZKqKAoYlZ85iEeE/vCvHH3+851IRRSn+Bl5UARVOOSuYv2yzkYH7BdPg+1UBaUN3P/jcSeqKRWRV3BZLuhBHOsjxIRCH2QUqxwvoYkggBehiSk4/tgqd1+AQGjHpR2pZuwq9MftXmvr9Ghr457rUrOofNT+Z2ci7O+/GzFAQp4Hv3RYyssJTtXHIcIqxCOq2daqXiucmiGLhFM/hyR3r5TMLCxM4ohMitSk6VY4c1b2X7SXXnFieDi5fkq575Wv6YN4quuOsE+iaNsfvF7HKuy0Duy7GEIbUfVmucgyECVEW6ylMMdZDp8a6lQE0wSl62NXGgd4QqPwlA16VefPmUVSwpQNzp3RgKhm4sr0EsRoPTFhFXZsdQY/9pYmxsrhqYwFjwUK66qqr6LHHHou8Yn9QgFawrUAS09VHlQPAnCQktv3IOT2wfeA/xmAJzC1fpbSTO/bZZ58V1c1l9yvELLcMwbARR29CPtLBC1PwE46ujicn0DWpWY4m3NqCah78hz3ET+btv//9b+Eh9FtKMBMoUaJE0t4AOwWMuNycSz6GAcaAiFeohWqODdtKINngT9iJduwQ9AHTgQHZrWygo8QBjoWLZoM71npVst+rElZavtt4MjNhOnAhHxPT0nWn91KwF+NgFw47VNwvnHJLdN3q5ZBxtV8uG0qdFI2Uao5C3OFSgbVq1Sr23csvv0zdunWjqBADjzHlOx3CSst3G0+WWHhBuzEtvEZwF/qkclg5LyKI5WDI2AyxU8P2wbs7q8iZZhx+wJ4jBt8XM0sIA27PEVQXN8bhahwdOHAg3XjjjeLCEOtQixQVz6NErjYiCkKHbKVFWGn5fsZjesiLAf+R6Ca7En/55Re65pprilX8xn/sxNjJIapjF8Z/+Twckw2BYFs81BblvB0wOXXOcrkBuMPdAHe4G1wZB5oxHXfccSLhp02bNiKiNOpObhbZj6B1OEwh6V7Hg1Sg6vEiO3ZfnAikDDCNf/3rX3TiiSca7SFYdBDr42bHOuCAA/brn2uqX+LU5hJ0gY3Di5riKY4DP8qXX34pmAd+hOXLl+dMf0yLcBGkDoeTmuFlPJNEwOoNxPQXXngh2UeE40V03ekhcbhFdUJ1KSoq2s84GRb8SpRqmLmTtCcbitVz/DJ7V4mjZcuWopjPxx9/LEoIgtjoIRslcrVDu1/kYrd6r9W+dAAtYKOQVRp1AWEHZaYhi+2cVaru8FwIxwRIODvTxDRUgyh3oPPSgoHdqDoGAObK9UfUe8P9c76LH7hKHGjAhEImzG3RnAltIKNEri2WoMhFOqTSNxb0wPFqvVIkrMEu0LBhQ+O5csMmtS6Hm1jvBDAqtrv4gZr/gvmDeYDJQRLyMh5sEzqm4eQ5wlxZ8uL2kF6YSEkvhHj44YepZ8+e4v2SJUuMJdXz0ZsQJXKxkE8qthGmhxxlCu9Jnz59hDThBJ13AdKGHAKPOWHx+sG6fUVwsG78FP/RpROAeaBeqBq5CikBY6vSkk7N0TFAXZUwMA81hyglxnHttdeKgJMZM2aI92gbd//991OUyFZPQtTIZq9KKtCFlwehB8RxGEKffvppsfuaFgA+13kXZJWA56S2O3UqbAxwy0lEdwYN9dZBVsNwzxhbfRZ0xXzUJENmkCZ4VRddVRVEjqJA8RtvvJEkaK7teha5EZaOBYFoZ+zeHBHKNgJ5HHUXxu7NJfVU4DPYSHg+8P4URhzXwX1v4WbW5dVwrQ1dKgJLT6zqMC3AAPEZpA9IVay2eVUXXRkHLgBRibkWGEmqZdRkoA7CV199JUQl7BQW+Y1UIlDZGCqHkXPPWDnyFAtftquYmAaOlRccdvlU+rUGAebHcSa6DRvf6yKJVeYrR90CMA7LEeCmbOTAjAPt7uBVQRDNlVdeSdOnT6eXXnrJdWAwgj//+c+0cOFC4VdnvQr65+zZs0Xy3LBhw+jss88Wf3/961+NY1nXb/7QwU8EKuiB5xJ/Z555ptHlyOoDMxH+DpsiF/uFJCEvGtWomO6+KhUqVNCqCLg+1B4d0+BQdFaN4MSA1KFjvoDuM5W5eIUr4zjrrLPEIp85c6a4KBa7l8rK4GYTJkwo1koBWbZgItOmTRPRqHDvNm/enIYOHWoMYYeFHO4iJhzXTQSg48kEZeakFrzB8ThP/ZzdmxhD1hdlt6csGmb62nw+JyWp11bnFOV9h3VtLGYei4OSOFZBFdPBMGC0R0ImGxDZLck2DJ4TGIPMePE5xvv555+T88H3aGiExQepBdflc3SxD0E+l7/jfBs26qr3x+eothquZib3kgXNsLZw72AkcnwHaMFtMfnauD8wHfW+eXz8fk45UZ4K+eBBhaiGC0OCAE4//XTHcyD6qFZlMB8wIgBtJGFwBfNAgBnHjOTDrmphBh5eVOLinqemEgLIKQHTGDFihLBByAsViwpJml5tcbKEg+siTsOtcE8lD4lnbsBuj+fdVFXdZEthxgGGLKtcTCswJJZS8DnoiP+yJCZvPnzffso1uDIOFO2BcRQ+cd6NMAE3xqEDOB8H4yCBbsGCBdS/f3/q3bu38ZwTTjjBWI8jiixR0zUycW3eVU1uvjDn5Pe+U702Hly4+fHQq/q37hrYIRFGjv8qPfCsoBSgWnXMBHnxYZFBzXZrbFRCksLSAdyr6jUBs8K8ZGMvGAxHwuK+cR6YCiQ31WYBmrC6hrF5/nzfkFhwvBcG4so43nvvPdEiIQyDKJgFx4Dgv5NbSFZVvNaCtMgfgygWCJiGCbJBULZryOAAK4AXmdfI1QIfEa5BAMapGkQh4chBavJ9qrEeOpuF2pdXV/3cK/1d4zggIYTVD6JVq1bC7sERqVBNLCy8hp3DpvH+++/7Sphjd6QayAU7na5MoS7uIWqU2efJhG0PQWgyo3CiD85x6l+jyzpGQ3nZPhVaHAeIiczYDh06FJM6EHruBDAb9GWZO3euaFqNMm2nnnqqECvRfxZjtmjRwnWCrKpY5G58hpewc06Nh3oSRHKRK3yDacAYyYtFjfWAJAx7Hs7JRMxS4b7qXFArsHGDHiw1qfRhmmK+iPNg9QMuZp26Z2ou7jfs37XmKIr26BBlIR9bASz7K4Clqo44BX3B0HfJJZcIpoEsbRM9dD1c8MeNqBlq3xaAjYu8W3OhH13afro6zusASQnX4BBzk8tYrgDmp3l00KC7rG/IBBsHHrI4LpZ8YBxeHzhd0WA1nNsEeA9wPs5xoodbQpef8nkAmA5SLPgesTh1hYIZbFNJF1QpQqapjnEwjVMpEB3rviogSjYtlkyBH45soYUfo2aqWbAq09DRA9eGmmGSECBpQKfnvi1qQJgKfCerMEg4M6GEIYKTgeulmrsCh4LMOGQa6u6BCzOHWQuWkfV52rBxZLlQlLfZsX5qafjNgoVN47XXXvNND7V3i65nLMeIQJrAXEz1LvieuJKY21y2G2xxuDYWPK7j5ub1W9LP5JrnMPpU6p04wfNdsE4YpH2dRW7CrxThNazZiyGUgcXKJf9YOlDrc3Ccg8nlCC8GFhrGgpGRA7twHKQanaRRal9krBztaSpsjPHVamPqWG5Zz3BMgElwaDxgUstABz6W58nHhRV74mrjmD9/vvgR8ePgUOhMMJiifmMU4DiOODchyuWGTGHrzxDHu3Tpsp8hVAc8kwgZ513cSRTX2VgwZ13+hpNK4oQSJUqIfia4h6AZtFxHRMdoZCatqmRyQyaMITNJzEm2zZg8LqFKHL169aLBgwfTGWecId5PnjyZbrjhBvr8888pKmR7I6J8pkOQOqNu4vjo0aM9lctDaLiaA8KiuE7iUKWjsMX4RCIhroPAtKDMx0nykOcKyJIE0wFMQY3LUG0rYEqQSFL53VwZB4jJTANo37592qPmVBtHWAFocUemA5OighemAagPPoeByyK8zCzgTVHL4+lULXVB8muO0Ew4COnM3GEo9fvc8vydbCU8V3Y1M1MEE8X8EJ0NKVA+FsxYdT9DWlHd4WrNjpQYB3zBKB149dVXi/evvvpqseKvUSCOXdrTgVylA2wa//vf/4R060cVwq6JZmG8kfFCUlPq+T8zDT7eFAyls5Hw97JrdrNGJXGyZeiAe4DawLkjOjexXLBYnSv/lwso6+4Ldhs5KU+2zeCeZAlJrl8SmHG8+OKL9OCDD9JFF10kuBeiPtF3MyrYXJXcZhxOhlA3VyJ7PLDw5M9NLQBkaUQeT5fTweNj4XG8icrESpUqlXLvWdir2N4AO4xOmpG7s5kKDsHeg+9gRNUxW9yDzDgggbB3SdUg3OqgeK5yroaXjxo1SkTyWVikAojkiEA2GUKdivrgoQfTwedIS5eZAJfKk2M0WE3xUiRIDSTj0n2ywbFmzZqhlBDEuOwpMdmwnLxVPFe5dojOe4Q/3AerI/L9q+N78ZC5xnE8/vjjnj5Lp43DIjeBRYN2oibvibqQ5PdOhk2uqSnv3qymyDEeas4HJ4ip6g4WG1QQ+XorVqzYL9DMrZCxCq71CakFqoJOzYGXxClwTKUD1CcTXdQYF75/rqvKKpGXOBujxDF27Fihd6Jy9K233pr8HBNLJYglV7wJmUAu0sFpZ1M9DPJ7nMeVu3S00cWY6HR/NgoyE2D3ZRAU7pNAwEB00ogaPapWTdMBtghmKDoXKt8n/3GQmy62xmTT4e9CKR0It84pp5wi0pibNWtWjEsOGTKEokS+eBPygQ5QL1AYqm/fvikFmOEhR49T7NgA/ssuRtMikReILrcFr3X5KH4S2go1TAMZueq4uorlXsPNGXyfUNfYq+JUeDgs97lrAFg2ZKbCYhxmZfW4gi3ncaWFbAh1C+7y4lXBTszRkX4T53RBYU7ARsr5LRt8Fi7G78XVyoNCDdpS6RL1s+Gqc2SaaVivSm54VbABXX/99Z6ZhrwwTMxAztPwmzjn1IRZB6gVOiNjuXLlXF2wQXvNwiSA9YfrgWkwTXTh86o0mo6M2Nil1duQ8+wNOfcCfoi5HqaX401uWHlB4DvQxE/PU3UMwK02qa6dgPwakgQzD7WuR1Do7tuUmwJJC8yEXdNAOjJiZZT2u+PB/6vL0ksXbAWweENlAnJMggkmN6w6FlyiWKR+GKmOKXHNClOKPa4DyOclpMhSjMHd3lgacAIHc8nX40Awk33C1ASb41Nk1zSYsxe3cypwdcdeccUVwjCDiyOxrUGDBjRo0CCKErnYpT0I5L4n2Qi13iceZlSx95sPYnKb6lyPkBLk3V2dgwp1DJwPzyH+w5uC7G91kUGiUc9j8BgMt5B0gMO/+XqwX6ACOQeyyX10+X441wbAf5wDSQNMjQuA63JXmH5udAld4kAfFUgYqI2AGqJPPPGE8LL069ePokI2L5Yokc10UHdyqCUIIUdbT78FfEweEdUmwTELWLg4XlYZTH1j1TG8GDq5O7x8HwklzsOPwRT2Hvl4ZiSsZphsGWqujU59UfNY1BT8sFSX0l5uEn9ok3DzzTeLByJqt2CWm2EiQzbTQd6R8XDfcccdSUNoEEOdKRScGYrccZ6Zh7p41b6xujG8NFViw6jMzNasWePpXK/JbpyxqjI/tVAQig+Z1BeogZBC1PB0OZQ9LNXFdQvDrlG7dm1xMTRhWr58eaQ2DhhH3Yqc5AtAh2ylhaxeYKdEWgJ7T1Tx2w1OYjWHSHP7Qt5FTdCpRzwfJ5giLPk+SnsMgvSTIcsRqiaACco0kWmO+4T3BhIKVC/TcUG8T4G8KnLHJ54gPosqetR6VeLjVQnDBeilRibHX3AAFRYzRHOd18FJNGfbBgPHIMBRlxWrm+fSgDU3dNAV59FBjlXhqFc8F6zumIoahe2edZU46tSpQ3feeSd99913yQlFGXJuc1XiARhC77//fl+ShQ5eiuuoO6hcMhALCx4O/LnVN1VzNxCJKqfdu91LiYAqO1QKdf7suXEbkzdxZMOCcYFx6ALLdPk7qf42Mlw5ABoqvfnmm9SjRw/hjr3uuuvosssui1RdyYVQ61ymg58aoWHUMZUzYGVbiGoXYQMiZ5862UzcjIhq7Mfq1asdWyk6AeHnsttVvq5bQ2v2oLgFnYWlkoQSADZlyhThnoWudfHFF9MDDzwg3EjpRjaEvWcDsrGvCh6fSy+9VGRMew0jd4NXsdqNHn5bA5h6v2B3lxdqCcWrIgeIgZnJBY914PJ+fH+663ppli1D7qsCiQZje6nklTaJAxP66KOPRPGen376SVjLr7zySpo2bRqde+65wgZhkb/AokHSWpjSUFiJWE71PLxKO2A+6u6eUCQN9GoBZCbllBUrjwcmIifF8XXVto8YwxQIyRW+1BaSXip5pY1xwMaBmqOI2zjttNOSn0PimDp1KkUBcPBs2mUzBTYGZhstMqVCudEjSPsGNX4E0oAJJfY1bcZxMLTqAsTYmA2Ph064l0PVYRzFsXx9kxrF6gxsjRyFyl4mGbriRzpJTo4b8Rq+76qqwFqb6V4q4LTZ6knIR68KbBoDBw6kkSNHZj09UvUmmHJEykjqAODkYeGWB7CLOKXSYyxZylClhZUrV4rveQz5GKhTCDdHkKDcF4a/N6ltpiCylGuOgquNGDGCFixYUEzcQi1Si/xDmIbQdMJLdq0XyIZY2eZQKHWUx+4vg0sNsjThp4CxSbUCY1DreajHwHDKXhdcX7ZxmNQ2XQ6MF7XO1R2L6ubgdJ988gm1a9dOtMJTCZVOwIaSrd6EqMFVnjIJFK72U08jE/TgXRQGR/xPNT8DCwgRm6bnfotkX8B8IBE4NafWAWNjsZsCtdibIkM+Rs2RwXt4ePjeTUFg8ue6cQOrKieddBJ9/fXX1LhxY5o3b56wZKPS+cyZMykK2G712eVVcYvUdELYQUgmepi8I0HnCgQp4OPGKPCHBS6rJ4jnYBVMlhhUzw6fzzYJ0AJR3Wo+kxowprNlBLFxeC7kA8PNt99+K4JkoKtFBduQKfOQmUUqTCPdNSJMiWx+a7WqizQdKF26tEhqU/NIZEYie0U4AQ6SB+JF4GWR6QnmgLKBWPiyhCXfO9PbS4sIN7iqKmj3CE6LpkydO3cWafWIJI0S2ZqfkQ+5KrBpwP2eaoJdOrqmm+iBRcCRmIAssrtB5351QqWAajuXA1CZGtQcE53APJA3hv8qPbnqunqfKn3C+h1cJQ6UewNg3wgzNj9XskJzmQ6yITRV24pf12iq9JAXjJ+MUD8L6fB9dUB1agQYgFsBKi/Xcuq14qXsoZfq76EyDjSadsLtt98e6IIW8cGTTz4ZmiHUqTR/OhB0gagLUledXGVOcu1TgG0KplJ/gDwneZ5qHphJwpTpKVcA4/YOPGev1d9DYxxsKV68eDHNmjVLqCkAGui0aNEi0MUs4oVhw4aFOl5YEaFerxVkgejO414lamGdClIWrQy8B/OQF7acwwKJRE44UwO9ZGOpE8NT6ck1R+X56KStMH4HV68KanAg5Jx1OTCU8847L7KoUcDmqkTnVUEiYyYrjfnxvGTCy6Sb3zYPhl+oM/jDXJlpuEVy+mF4Mi2iMES72jjg1pJb2+F1qo12LbITEHdvueUW+u9//5uRrnFhPvDpag+g2+ULCgqKlfUDuE4o17NhGwiCKLH5wj4iSyFcrFjXdMnv/UWhFroyDhjHoJpceOGF4j1KCHbv3p2iRJz7icSFDrIhNFOtJv0mpZnoEZXrd5vmOm62DQakD/kYLhmIaNR169Y5zl13XfU3S7da6Mo47rvvPlGkGNmwALJkERQWJSzjSD8dnn/++YxHhPo1aJro4ZcB+cW2fbs9FrlaMR0qgxfvlyk2SW0aDakE0oypVQLfX9T5ZLYhU54nuYUt0oeRVOb1fBM9vEocQebq1BgpCNRixgjiQrKaCrliu+7+2F0e1TqJtu28RVYBu5vsxktVpA9DRQhDxPai4wedq6kxkhfAjoH4DjlDVle7A8Zpt4rt6v1F3bQsext1SCHnmdK5sw2gQ1i0AMO44IILkkVkwojmTEd0aFB6uNXYDDJXMBvQq4SHYDhd4hgMolz0hwFXr5p8xqn6KuR5qvcX5rORMxJHpjNCc40ObAgN0iwp6ujQVOlhUkdMc3VykcoqSlUldV6GnIAme1Z4XDAQOeeEg7zkpDYcoxvbxBzcaBG2SmpkHLhxp8no0nzTBWscDZcOr7/+ekrNkrIlOtSNHm7qiFqxy+l4XfOjg/f1ekXSJ7cnAOBuxXunIjrcfV7+XJY05D6yYCymSFAvtEiHl8k1chQFiWHVRV0OXBitIN2a6oYNyzjCpcNdd92VNrddlNGhbvRw8q6oqgkWPwyVpuPVUPTCwkKxIPE91A+VsZiK5Xj5XKWlXL7Q5CVyokU6vEyuNo7333+fevfuLSQQiFY33ngjjRkzJqWLWkQPXe+NXIdTBzO1gA0kBG4pqTuepSnuxbJhw4ZkkSD5O935XoroOKl2qXZiS0cnN1fGgYtAyoC4BK6G1+nWXS3CBWwaSBMIq1O5jDC6oAcdA54EBEvpzpMjOnWNmXixy/EP2I2hBmCTRP0ZFTgHKkpCY1TlKmG6RlAyY+FxuQsb7CQY0ykOg893azBloqff80OJ40BLhNtuu42mT58uuFXr1q1p6NChoi5AVECF6LJly1K+A3QA/NBCjggNO7grDN056Bg474cffhDnIaNUbZ7kdUynuAzduanMd5mH+A8wniALm58NbPBRRM26elXAIDKtmmQy6SqbEIQOyGZOV0RoGLpz0DGczsNO7jSmahDG4lINnHwuxpKPVQ3AgKlTnHw9tfObCSzByOPrXqvX4mdDjTwNO2rWM+NA5CbsGkhsQ+lA1B2F3QN9QqNClge3ZjUdYJ9KF8Jwv4ZRN0N1p8rBU+qYJolBZ+AEdO0Kyu/78yJ9+I00hbtVHlM+T75fk4coKpe46xbWs2dP0d6PU5dRtBi9ZKOELR3or3SgU/vBVO0O8udh6M7qGIAXewfOO/LII0XRGtVtKgM2Bbc8Dx5PtWvA3qEeK99/gUMQGR8nSz+AnGkut2yEXYXjO1Q3McMpYI2fDdNvEoYtypfEgQuphXui7FZv4Q9s00AWM/IegsBL855UCt2q8LODy4Bhkf8Y6o6rRmE67chqQyTEUsiFe1RpoGbNmkapR9cOEq9RUUwO7OIFj3wVleE7SRx+C/xEFsfBwI3++OOPSQK88847xQrBWmQPZENoUKYBBI07YAQNKvM7PhaymtTlFoRm+t5UV4OlEDAVdX579uzRjqUeJzMvSDF4r+aiqNIFNmeoT2pX+3TSNVTGgS5uqHS+aNEi4W4CoeCSjRI25NwbHSZNmhSKIdS0K3vRn1PZ3fyOj4V79NFHa5mH0zW97MiAGtWphnuXKlVKey3dcVxjAxINmBNLNibpAolwkHYwF7knjOm+nJ4N3XxShSvjqFWrFo0fP15wKcRxRNnFjWGT3LzRAdJGGDDtyl5CylPZ3fyO7zWwzU0CMtkq3JLg9hjsTern8NbIY0HaYMlGboKkNn3y4hHy8myo8wnDZujKOPBDnnPOOXTppZfSmWeemfIFLcIFdjLo4mFLZaZd2203T9Wq73d8L8FQbhKQac5eP1OhjqcyNwSu4U8Xs6GqMDJDCCrNpcPT4hoAhsl++OGHwpMyZ84c+vOf/0yXXXYZtWnThqKC7VavL1zDNo233nprv3TtTCJd9T5NNg6nZ8NrO0hTEWIvn7lVCTO1jtTN5bfffit2vHyM0724FXkK+zdxdcfiIl27dqXRo0eLHrJwF6E5U1hAXAjqQsCGYhHMEJpNTENXKyId42MBYDG4FbDR5WnoXJO6OfNngOyCru7h3vg4U20NnpsKzoWR5+t0L5n6TTz5VadMmSJ2tY8//phOOeUUevvtt13PgcsJ0snChQuFjscu3D59+tDs2bPp5JNPFn07EBcCxmHhD6BhpmuEZgosssOFiQWkSzNn6CI9/Yj7qboyy5cvL1QSSBFc+QtrQa694TRfXX5NVGULUg45R3FiSB2DBg3yzOWgd0+YMCFZHR2AqgMmgsLHiEZFo6fmzZu7jmW9KvvTQaZrvsGvAVa2m6hNnt3ODcOVWd5nrIvT8abvol4jjowD1tfrrruO+vfv73tg7AJqa7yZM2fSWWedJV537NiRZsyYIcTscePGidD2Bx98cL8ELnwOCzT/eByIw/OTTTQs1cD7I9cnwPE4T/0c8f34Y9ee+jmg1ofM5LXhnkOkJKC7tjqnKO87rGurn8vX5jmxJ4KlAMxBVllM1+YoTXmu+AwJYnxtbGysknA2LI/Dc1m/fr04B9/jefVy3/y5eg1sxOp9s/EUGbM4xstvgWPwOQodY3zZ9uP3t+BrODW6crRx4GQYRsMCfNeo6QEgQAnv4e5FRarHHnvMZsC62DTQzwbJUvkOLAhUvkdwYpAK+LKdQAYWK+iMalv4j/dYuHhGYR/AtYBffvlF/A7Lly/3XSQYDEO9hm4OkIwwvp8QcRyLc3CubuxIVRWk0d98883CHSurKbBR+AWYBZccxH9dzQNdsWK4s1TpJarWf6ZrRHltPAQIwnvuuedEfoaOFmHPye99R/1bgFmwyG6ihw7YkWXJCjstb1jwgMjuT+ziuCaeW47E5Q5t8jF+GBfsMvL5antTSAvy95iTKQpYPg9rBGtKvbd0eSNdGcc333wj/svqCjj1xIkTfV+sVatWovEP7CUIKvPSEY77quQz4HGCIRRJUOnIFJbdmxyMlEnDm1cEoYVTTINqv9PZ87wc4+QGVa8PmqvtIoPEXOB4MAmvDavTzjgQxhwE4KToADd37lzq1KmTUEVOPfVUsTu0bduWmjZtarveewTbhdIheurSvtNZACbTcPNMyIvWhBIejnHyyJi62PPYcr6Mn98AjCMqr4unptP33nuvcK+OHTtWuFdh1OzRo4fjeRCjIFWogAvWD6CqRN1sJhsAeqsp4emArgZFOgvAZANMngkvHhQ/XpYCw7H8J3t4GGyo1AWppXJvYcM1AAzqBCQGTgXGQkbpwKgAVSXf0vhh07jyyiuLVbcGQIewaaEW7QXSLeaGhbDp4SXAyk8QVgWXY8OkfTqeDcfruR0A6y9sEijmI04oXdomnUUUEQprfrqhis5xsnHIkCt/mYKr5GP9BFjpSg0W7KOXXAhINz7UDm6+pJsTOwgwFgycpuNSoUs6VBdXxoELIpGKOSNiMVKp9eAXkHC4EGs+ALkKpojQdPWXiUq8DRtMDyyOpUuXJj+Xe6z6jQT1UwRnmUuxI9kegsUrR7iq4/Lx6nF+aeHnXtPKOAYPHkydO3cWxXzgmoX4jGI+USKfGjK1bNnS+F0+0cELmB5qyruT7cFvJKjfokYFmsZMumup53uZuxdapHKvoTIOxGsgV2Xx4sXi4nXr1o3EZ59P7thff/1V/KAI07fwD9Um4Cfl3c2e4LeoUQWl45uT61f+Tj4+LPtSOgsXuzKOUaNGiXocDRs2pEceeUTkm6DCeZAAMAuzTePFF1+0jMOnrcKUSOZkJGQbhSnV3XS8zk5QsWJFEXaA303+XC43COjsIfiT7R9yiUAv9+zVdiHPJUx11LUeB7JXkfr+2WefiT6yffv2pYEDB9IXX3xBUQHRc2p16FyA32ZJoAOQi7TQ2Src9HKVHl51+lR1/22GeQJeizybjvV6ffVcZpZ+aZE2dyx7UD766CPRKgGtBPkHiwq5mh2L3cZPajzokKu0cLJVeKWHU7sC9TpejvM7zwLDuLrPU5mD7tygtEgb40CB4l69eol6HOeee67wcERppIONI1eNglD//NTTULMb89lWoaNHlE2cZfAYpnF1n6cyB925QWmR1tKBKODTqFEjqlOnjgiRnT9/Pp199tkUpXE0V0oHQj0Bghh83crD5ZuNQ0cPv6X9gsY3bDPM0zRuKmUIvcw/FVqkhXEAMIjCxgHOBZds1IbRXKk5yjaNF154QRjz4sI40l1DNCiipodKh20pJAcGYSRO30dNC1evCgyh8KxcdNFF4v21115Ll1xySaS9Y3PNEBqEaWQK6TayxZUONWvW1CaoeaFPEGNptv0OrhIH4jaQ4co1D8DZkNmKuI6oALtK3Iv8/Pzzz8J1l0qNUI6gjZIWXquEZwJR0kOlQ8WKFffrbu+VPjqaAk50dvsdon42XCUO1ICQC+lggjCYRgkuTBJnoONYHOmQziCiONFDpUPlypX3ixL1Sh8/fVzczsnUs2GUOG655RYxQeyUKCqMmhB4/+mnn4o6GmiXEJVxFAWT4xi7APUEu1L9+vVDGS9TcRzZauOImh7ZbOMojJgWRsbx8ssvO57YrVs3igJx9ar4De7ygnzwqviBpUfmaOFq44Ca8sMPP4jXxx9/vK/6jvnqVUkH0wDsQjHTQ96NAZ1koPtOBydJwo9UAEQlqWUN40DMPyp/IYcCdSFwGBYEvCqPPvpopIlucWMcKEOAorhhe08s49DTA8+mLpVd5/2Qv9N5JpxKKfrxfHi5VpyfDaNFpV+/fqJ/BAjx1VdfiVgOpNaj1gHyVaJE3Iyj6H2SDperXMHa4g96OKWyI6zfKc3dTylFtzDuApeU+nQi6mfDeCX0Uxk5ciRVqlQp+RksySjR/7///Y+iRBwWC6SxL7/8Mq3XsIxDTw+1BJ8cao1n1vSdW2lA+RwvYeIVHOaRbm9U1M9Gab8JVVy+PSrAOMpiYhxsGhbRQ9cfVrYtyCnr6ndOY+lsHE5VxMu7zCOXYGQcDRo0oFdeeUUsCBmvvvoq1atXj6IEN8bJ1lJ/6TCE6sCtALOVFlFDpoda8s+pHKDbIg7SuzXoteL6bBiNo1gQCDOHsaVZs2bJDukwwrz77ruRBoFls3EUOwoKOkdRWNgaR4vD0iMLvSoMdGxbsGBBUgrp0KEDRY1sZhz5slCc3JBRBojBYM9Vs5gOOnq4BVgBTvcjZ77ienzNww47zHWO8vmYm5NLF/CSDexGf9AlyvKTnrJjMw1T79hM2jRge4maiYIOQNS0cEqwijL5CouD+/sAaDqN0ngqPbwkkckuU/V+5OpeKnDNwxyYh86dK19HnocO8MY5uYhN9EceFL5D6YusaMiUDcimPi5sCEUYfCbokAlaOLkh011pSgY3LJevraOHl4pbDN39+JmDF3euaR46uLmITfRnh0W63b6M/GqRliKQnRiVITSb4NaoOaokOKgKckaqHCrgZb5eqoq7zR9zcIJa5Zyhm4fp/CD0V6+RblhVxQcgDq5evTry7OBMqyrZauPga+nokW82jg0bNojPubJ6uhELxmGNo3thvQjFYemRhSHn2QIYITNt03jvvfcyOgcLi2xD1jOObDCEolCzhYVFzFQV2BaijpZEhuvFF1+cVYZQ0AGwkaN7YemROVpYxmEAyAJDKNeDzAbYhVIclh6Zo0Us3LGZaEIE11Y2MQ0grs2YvHhd/PQjCUKPdHl+tgWYdzrmG/TZCNKrF7CMQ7FpTJkyha666irKRsSRcXiJLDUd43auV3qkK7p1W8B5p2O+QZ4NNUoWtXa8ztUaRxVDaKtWrTI9lZyCl8hSPz1X0zWHMMctCKE3bRTRuH579cqwjGNfGHE+RoRGAS+RjX56rqZrDmGOWyGE3rRRRIP67dUbO+MoSr+nu+w7xDQkTGUzMtUeIVttHH7okes2jsKAz0ZQG0csGEcmvCrZCOtFKA5Lj8zRIhaqSjqMgrBp/OMf/6A4AXSIo4E0XbD0yBwt8pJxsCEU3eniBLtQisPS4w9YxpHmXBVkmHbv3t0aQi0sUkAs4jjCBFKw33///axqnmxhETfEwji6c+dOKlu2LOU7QAfA0mIvLD0yR4usV1WAVBvNwKYxaNAgijtsQ6bisPTIwk5u2YRUhCI2hKLVQ9wBOsRAQIwMlh6Zo0UsbBwovRa0Sc3111+fM4bQoHTIVVh6ZI4WOV86MIqo06hgS+UVh6VH5mgRC4kjFeQK0/CLKAsIpwIvhXtNxX+D9DDWhVg7hY3Lx8KVz0WLD5YaHzmFbQcN6Vbpw31sdf1sw7p/P2PnnFcFNo2RI0fSwIEDKV8t51E2SUoF6jzlR1FuYKRrkITvUW0ejMTrs6GmkWOMmjVr0u+//65NjXdqzHT44YcL5qFrwMRNlZy+8zpfU/Mm9Xf18mzoxsN98P2bxo6tcdSrtZgNod26daN8tpxH2SQpFZiaJKkNjHTA9xDP/XgSdGnkkCBMqfFeGjMVaBos6dLrTXNwm69pX1d/Vy/Phm48+f5NY8eWcXgBwm1vuummnDGEpoJMNelJdZ4yeN6mueN7v1KULo0caocpNd5LY6YKGpVJl15vmoPbfE3qWKrp+gz5/v2MnVPGUdYBcxV+DGD5ZOPwYxDMVRvHdo/PRl7ZOGxDpr2wXoTisPT4A7Yhkw/AptGnT59MT8PCIu8QC3esTs9jQyhsGvmCIO7HXIalR+ZoEQvGodotoF3169cv7wyhuWy/CQJLj8zRIralA9nvnk+wpfKKw9IjjxsyzZo1iz788EPhTn344Ye1x8DSqxIk35iGnI9gF8peWHpkjhZpM46uWLGCTj75ZFE4B8lmDBgz27ZtS7fddpt4j07wDz30ENWoUUNEsOkA6QI2DSSs5TNsNmhxWHpkjhZpYxzwc0+YMIFatmyZ/GzOnDm0detWmjZtmkg+g7ThBb/++qswhN5zzz3pmq6FhUU2qCqQNPAnY+bMmckCwR07dqQZM2ZQly5dRF4JRC3kDehqjkKFeeaZZ0RcPYJw2BCEc2QuW7p0aW3hVhwP1Ub9nMN0MYacliyH78rSEsbI5LVZj1XzOvja6pyivO+wrq1+Ll9bnhNgoofp2vx5GNcu7fO+w7q26bfgdYF5pfoM8jWc1J5IbRxoeoQkH6BKlSq0YMECatGihfgzAUk5DzzwAJ1wwgn7fRdEn+Pix7rxTNfXHS9f2zSm18+d3suvf/75Z8e5h0kP01hu9EgnLXAN+T3UV92YQeH0bOjo8b3Ls+T0ve47P7RQf4soaRGJV6V9+/Y0fvx4wQlHjBhB1atXp65du9Lo0aOFCnLrrbem8/IWFhZxjxxFQ2fYPQAwE9n+YWFhER+kjXFA/4QdY+7cudSpUyf64osvkl4WeFWgQzmpKBYWFtmLWASAWVhYZBdineRmYWGRGcSScSD+48EHHxTelnzHvHnz6IILLqBFixZRvmPcuHH0+OOP0y233EL5ju+++44GDx5MN998c1p6ymYV4wgz2jTu8EqLxo0bC8aRy/BKi7PPPlsECXJ8Rz7Ton79+qJQD4r2pCM9I6sYR5jRpnGHpUUwWgwdOjRna876pcW1115LHTp0SNZHzakkt3REm+YCvNICUhdEdATsQH3LxT6qXmmBBfPll1+Kz7CwcjER8kCPtNi0aZNgKD/99BNdddVVuc04woo2zVXoaFGrVi16/fXXKd+go0X//v2pd+/elG/YqKEFGAj+0oWsUlV0ACFY1ML/gw46iPIVlhZ/wNIis7TIesZho03/gKXFH7C0yCwtsopx2GjTP2Bp8QcsLbKPFjZy1MLCIt4Sh4WFRTxgGYeFhYVvWMZhYWHhG5ZxWFhY+IZlHBYWFr5hGYeFhYVvWMYREuA/b9q0afLviSeeiCTU+Nlnn/V9HjKLn3zySdfjKlasSOkCktG2bdtGmQCyqxcuXOj7O7WW7uzZsylfYRlHSChXrhx98803yb+77747axlHNiDujCPfYRlHGoEMxbp169LixYvF+8svv5xGjhyZ3M1RQ6Fhw4Yi9XnNmjXi8x9//JHOOeccatasmYgE5AI9q1atogsvvJCaNGki/j7//HPBnHA8JBw04QYGDRpEzZs3F3U6kC3LePTRR0Wp+zZt2iTnowLl9hG+3KhRI7r//vuTnyNlG3NEhCK+GzNmjPgcSWVgAIz77ruPhg0bJmqknH766WJeJ554okj3ljF8+HBRV+KMM84Qf8Abb7whxsbxd911l3Z+tWvXFvfE82DaFBQU0HXXXSciJk866aTk/FCbAlnUwCeffCLmBLq9//77gl6YH+jH0H2HTQAh3KAn6I/6FjJQJKd79+6CXsjWxrlM/+eff14cM3nyZCGhXHzxxVSvXj268sor49+BDpGjFqmjZMmSiSZNmiT/3nzzTfH5uHHjEi1btky88cYbiU6dOiWPB+lfffVV8XrAgAGJm266Sbw+88wzE99//714PXPmzMQZZ5whXnft2jUxZMgQ8Xr37t2JjRs3JpYtW5Zo2LBhcsxPPvkk0bNnz0RRUVFiz549ifPOOy8xZcqUxOzZsxMnnnhioqCgILFp06bEcccdlxg0aNB+93D++ecnXn75ZfH6mWeeSVSoUEG83rVrlzgPWLNmjTgf18D1TzrpJPE5rnfssccm1q5dm3jyyScTjzzySHKumzdv3u9atWrVEmMBv/32W+Koo45KrF69WlwL9/zuu+9qzxk+fLh4PWLEiESPHj3E63vuuSfxn//8R7zesGFDok6dOomtW7eK+23QoEFi4sSJiRNOOCHxww8/iGO6deuWGDVqlPZ3VL9r1KhRYvLkyeL1Aw88kLjtttvE63bt2iVmzJiRuOyyy5L3+vzzzycefvhh8XrHjh2JZs2aJZYuXZqYNGlSonLlyolffvlF0AnPw7Rp0xJxRtan1cdNVVGBOgmjRo2im266SeQXMNBB69JLLxWvUS/hoosuEjs7dr1LLrkkedzOnTvF/4kTJ9Irr7yStKcgI1Ld/VCXA3/YdQGMt2TJEtqyZYvYLcuXLy8+79y5s/Yepk+fTv/973/F66uvvjq584PP3XvvvTR16lQx799++01IQJAADjnkEPr666/Fe1wX77HjQgJAXgWqk2H3dgLqaGBHRs8dADsyrqWrbAY6AZDI0JuH7xuSAttt0NUMzatQBQsSHiSNIUOG0HHHHUd+JUaog+3atRPvUSBI/m169eolegRB0uJ5oJTjO++8kzwf9C9TpoyQho488kjxOeiBOhmQ/uIKyzjSDIiyqP+IRYuFzg+PCm7Rh5RoHQPyAixwlM7DAy1DVifcoCt+89prrwlV6quvvhIdzcAwsDgBNAJ/6aWXaOXKlYJZAFioWPgfffSREONvv/120fs3DHChIjBPLp2H+wbDg1qoYv78+YKZQTUKG6eddhpNmjSJ7rjjDpFkhnk8/fTTIvlMBlQVucCSPPe4wto40gzsdNj5UGwHpdy4HiaYBO9M+A67T+XKlemYY44REgqAB5GlFNgYnnvuOfEaujR2s0qVKglpgoEH9sUXXxSSBgDJYPXq1WIhw+i3fft2cfwHH3ygnWvr1q3pzTffTDILBq6FSmNgGlgoy5cvT34HSebjjz8WUgMvGHx/6KGHUs+ePQVjQSUqFfLcsRtPmTKF1q5dK+4N9g7e5b0A18WCZbsBJCCex1NPPSXejx07VmSSqtd2mhekOtTtZBvNf/7zn2Lz6tGjB5177rlC6gAjwDzwG/FvjKpssL/kJDKtK+WqjeOuu+5KLFq0KFGvXr2kjt+nT59E//79xWvYD/AeNgro9NDvAejEsIU0btw4Ub9+fWH/AFauXJno3LmzsFVg/M8//1x8fvnll4sx+vbtK94PHTpUHIM/6NKs10MPh+7funVrcY7OxoFr4xyce9999yVtHLBF8Ofdu3cX9wT7BqNXr17ifhkvvfSSmFPTpk0Tbdq0EeOqgK0Cdof27duL96+//roYH+fdeeedWhrLdpFZs2YJOwOwbdu2xA033CDOh00Dth3YYDp06JAYM2aMOIbtPNu3b0989tlngraYH9OHoX739ddfJ0499VRh6+jSpUti/fr14jhcG3MA8JvC1gH7BewtfB+4N9iiYOPAnBiwZ/373/9OxBk2rT5DgFeFJYM4A5ITvByQkurUqZPp6VhEBKuqWAQG4h2OP/54oUZZppFfsBKHhYWFb1iJw8LCwjcs47CwsPANyzgsLCx8wzIOCwsL37CMw8LCwjcs47CwsCC/+H9oOVWKrcFjqAAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from scipy.special import logsumexp\n", "\n", "# Calculate expected waiting times from model predictions using competing exponentials theory\n", "# In Delphi's framework, each possible event has an exponential distribution with rate λᵢ = exp(logits[i])\n", "# The expected time until any event occurs is 1/sum(λᵢ) = 1/exp(logsumexp(logits))\n", "# logsumexp provides numerical stability vs. calculating exp(logits) directly\n", "\n", "# Let's see how the predicted waiting times compare to the observed waiting times\n", "\n", "plt.figure(figsize=(4, 4))\n", "# Calculate expected time to next token (inverse of hazard rate)\n", "expected_t = 1/np.exp(logsumexp(p, axis=-1))\n", "\n", "# Define bin width for logarithmic binning\n", "delta_log_t = 0.1\n", "log_range = np.arange(1.75, 4, delta_log_t)\n", "\n", "# Calculate average observed time for each logarithmic bin\n", "observed_t = []\n", "for i in log_range:\n", " # Create mask for current bin and valid times\n", " bin_mask = (expected_t > 10**i) & (expected_t <= 10**(i+delta_log_t)) & (t > 0)\n", " # Calculate mean for this bin\n", " bin_mean = t[bin_mask].mean() if bin_mask.sum() > 0 else np.nan\n", " observed_t.append(bin_mean)\n", "plt.axes().set_aspect('equal')\n", "plt.scatter(expected_t, t+0.5, marker=\".\", c='lightgrey', rasterized=True)\n", "plt.xlabel('Expected days to next token')\n", "plt.ylabel('Observed days to next token')\n", "plt.plot(10**(np.arange(1.75,4,delta_log_t)+delta_log_t/2.),observed_t, label='average')\n", "plt.yscale('log')\n", "plt.xscale('log')\n", "plt.legend()\n", "plt.xlim(1,2e3)\n", "plt.ylim(1,2e3)\n", "plt.plot([0,1],[0,1], transform = plt.gca().transAxes, c='k' , ls=(0, (5, 5)), linewidth=0.7)\n", "\n", "plt.gca().tick_params(length=1.15, width=0.3, labelsize=8, grid_alpha=1, grid_linewidth=0.45, grid_linestyle=':')\n", "plt.gca().tick_params(length=1.15, width=0.3, labelsize=8, grid_alpha=0.0, grid_linewidth=0.35, which='minor')" ] }, { "cell_type": "code", "execution_count": null, "id": "4e5e3c8f", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "base", "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.12.9" } }, "nbformat": 4, "nbformat_minor": 5 }