From 8f44018baeb68b23c6c3672c647f31fce4c27baa Mon Sep 17 00:00:00 2001 From: Jiarui Li Date: Mon, 20 Oct 2025 13:47:50 +0800 Subject: [PATCH] update evaluate --- evaluate_auc.py | 49 ++++--------- evaluate_models.ipynb | 162 ++++++++++++++++++++++++++++++++---------- utils.py | 47 ++++++++++-- 3 files changed, 182 insertions(+), 76 deletions(-) diff --git a/evaluate_auc.py b/evaluate_auc.py index 635c3d7..7329ec5 100644 --- a/evaluate_auc.py +++ b/evaluate_auc.py @@ -2,12 +2,12 @@ import scipy.stats import scipy import warnings import torch -from model import DelphiConfig, Delphi +from models import TimeAwareGPT2 from tqdm import tqdm import pandas as pd import numpy as np import argparse -from utils import get_batch, get_p2i +from utils import load_model, get_batch, PatientEventDataset from pathlib import Path @@ -350,7 +350,7 @@ def evaluate_auc_pipeline( total=d100k[0].shape[0] // batch_size + 1, ): dd = [x.to(device) for x in dd] - outputs = model(*dd)[0].cpu().detach().numpy() + outputs = model(dd[0], dd[1]).cpu().detach().numpy() # Keep only the columns corresponding to the current disease chunk p100k.append(outputs[:, :, diseases_chunk].astype("float16")) # enough to store logits, but not rates p100k = np.vstack(p100k) @@ -422,13 +422,6 @@ def evaluate_auc_pipeline( def main(): parser = argparse.ArgumentParser(description="Evaluate AUC") - parser.add_argument("--input_path", type=str, help="Path to the dataset") - parser.add_argument("--output_path", type=str, help="Path to the output") - parser.add_argument("--model_ckpt_path", type=str, help="Path to the model weights") - parser.add_argument("--no_event_token_rate", type=int, help="No event token rate") - parser.add_argument( - "--health_token_replacement_prob", default=0.0, type=float, help="Health token replacement probability" - ) parser.add_argument("--dataset_subset_size", type=int, default=-1, help="Dataset subset size for evaluation") parser.add_argument("--n_bootstrap", type=int, default=1, help="Number of bootstrap samples") # Optional filtering/chunking parameters: @@ -436,10 +429,7 @@ def main(): parser.add_argument("--disease_chunk_size", type=int, default=200, help="Chunk size for processing diseases") args = parser.parse_args() - input_path = args.input_path - output_path = args.output_path - no_event_token_rate = args.no_event_token_rate - health_token_replacement_prob = args.health_token_replacement_prob + output_path = './' dataset_subset_size = args.dataset_subset_size # Create output folder if it doesn't exist. @@ -449,35 +439,26 @@ def main(): seed = 1337 # Load model checkpoint and initialize model. - ckpt_path = args.model_ckpt_path - checkpoint = torch.load(ckpt_path, map_location=device) - conf = DelphiConfig(**checkpoint["model_args"]) - model = Delphi(conf) - state_dict = checkpoint["model"] - model.load_state_dict(state_dict) + model = load_model('config_n_embd_256_n_layer_16_n_head_16.json', + 'best_model_n_embd_256_n_layer_16_n_head_16.pt', + 1270) model.eval() model = model.to(device) # Load training and validation data. - val = np.fromfile(f"{input_path}/val.bin", dtype=np.uint32).reshape(-1, 3).astype(np.int64) - val_p2i = get_p2i(val) + + val_data_path = 'ukb_real_val.bin' + + val_data_arr = np.memmap(val_data_path, dtype=np.uint32, mode='r').reshape(-1, 3) + block_length = 128 + val_dataset = PatientEventDataset(val_data_arr, block_length) if dataset_subset_size == -1: - dataset_subset_size = len(val_p2i) + dataset_subset_size = len(val_dataset) # Get a subset batch for evaluation. - d100k = get_batch( - range(dataset_subset_size), - val, - val_p2i, - select="left", - block_size=80, - device=device, - padding="random", - no_event_token_rate=no_event_token_rate, - health_token_replacement_prob=health_token_replacement_prob, - ) + d100k = get_batch(val_dataset, slice(dataset_subset_size)) # Load labels (external) to be passed in. delphi_labels = pd.read_csv("delphi_labels_chapters_colours_icd.csv") diff --git a/evaluate_models.ipynb b/evaluate_models.ipynb index 2a3845a..1b18ed9 100644 --- a/evaluate_models.ipynb +++ b/evaluate_models.ipynb @@ -7,7 +7,6 @@ "metadata": {}, "outputs": [], "source": [ - "import os\n", "import torch\n", "from models import TimeAwareGPT2\n", "from utils import load_model\n", @@ -166,8 +165,8 @@ "output_type": "stream", "text": [ "Using device: cpu\n", - "Model config: {'n_layer': 12, 'n_embd': 120, 'n_head': 12, '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_120_n_layer_12_n_head_12.pt with 2.40M parameters.\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" ] } ], @@ -178,8 +177,8 @@ "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_120_n_layer_12_n_head_12.json',\n", - " 'best_model_n_embd_120_n_layer_12_n_head_12.pt', \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)" @@ -220,7 +219,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 5, "id": "d7f2e4a1", "metadata": {}, "outputs": [ @@ -250,36 +249,29 @@ "42.0: No event\n", "=====================\n", "Generated trajectory:\n", - "44.9: E03 Other hypothyroidism\n", - "45.2: M19 Other arthrosis\n", - "45.7: H60 Otitis externa\n", - "47.1: H10 Conjunctivitis\n", - "47.6: H66 Suppurative and unspecified otitis media\n", - "48.4: M25 Other joint disorders, not elsewhere classified\n", - "49.8: M79 Other soft tissue disorders, not elsewhere classified\n", - "50.9: M65 Synovitis and tenosynovitis\n", - "51.2: D50 Iron deficiency anaemia\n", - "51.9: K57 Diverticular disease of intestine\n", - "52.6: K63 Other diseases of intestine\n", - "52.8: K66 Other disorders of peritoneum\n", - "52.9: A09 Diarrhoea and gastro-enteritis of presumed infectious origin\n", - "53.8: K20 Oesophagitis\n", - "53.8: K29 Gastritis and duodenitis\n", - "53.9: L30 Other dermatitis\n", - "54.5: M54 Dorsalgia\n", - "54.8: M15 Polyarthrosis\n", - "56.1: N30 Cystitis\n", - "56.3: F07 Personality and behavioural disorders due to brain disease, damage and dysfunction\n", - "56.4: K80 Cholelithiasis\n", - "56.4: M13 Other arthritis\n", - "56.6: B37 Candidiasis\n", - "58.0: E11 Non-insulin-dependent diabetes mellitus\n", - "58.0: K59 Other functional intestinal disorders\n", - "60.5: L03 Cellulitis\n", - "60.8: K52 Other non-infective gastro-enteritis and colitis\n", - "61.4: K31 Other diseases of stomach and duodenum\n", - "61.8: K74 Fibrosis and cirrhosis of liver\n", - "61.8: Death\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" ] } ], @@ -310,10 +302,106 @@ }, { "cell_type": "code", - "execution_count": null, + "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": [] } ], diff --git a/utils.py b/utils.py index bab1538..ef25e81 100644 --- a/utils.py +++ b/utils.py @@ -42,17 +42,22 @@ class PatientEventDataset(torch.utils.data.Dataset): """ return len(self.patient_ids) - def __getitem__(self, idx: int) -> tuple[torch.Tensor, torch.Tensor]: + def __getitem__(self, idx): """ - Retrieves, processes, and returns a single patient's event sequence. + Retrieves, processes, and returns a single patient's event sequence, + or a list of sequences if a slice is provided. Args: - idx (int): The index of the patient to retrieve. + idx (int or slice): The index or slice of the patient(s) to retrieve. Returns: - A tuple of two torch.long tensors: (event_sequence, time_sequence), - both of shape (block_length,). + If idx is an int, a tuple of two torch.long tensors: + (event_sequence, time_sequence), both of shape (block_length,). + If idx is a slice, a list of such tuples. """ + if isinstance(idx, slice): + return [self[i] for i in range(*idx.indices(len(self)))] + # 1. Retrieve and Sort patient_id = self.patient_ids[idx] records = sorted(self.patient_events[patient_id], key=lambda x: x[0]) @@ -150,3 +155,35 @@ def load_model(config_path, model_path, vocab_size, device='cpu'): print(f"Model loaded from {model_path} with {model.get_num_params():.2f}M parameters.") return model + + +def get_batch(dataset: PatientEventDataset, batch_slice: slice) -> tuple[torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor]: + """ + Retrieves a batch of data from a PatientEventDataset and prepares it for model training. + + Args: + dataset (PatientEventDataset): The dataset to retrieve data from. + batch_slice (slice): The slice defining the batch of patients to retrieve. + ignore_tokens (list, optional): A list of token IDs to be ignored in the target events. + These tokens will be replaced with -100. Defaults to None. + + Returns: + A tuple containing four tensors: + - input_events: (batch_size, sequence_length - 1) + - input_tims: (batch_size, sequence_length - 1) + - target_events: (batch_size, sequence_length - 1) + - target_times: (batch_size, sequence_length - 1) + """ + batch_data = dataset[batch_slice] + + input_events = [item[0][:-1] for item in batch_data] + input_tims = [item[1][:-1] for item in batch_data] + target_events = [item[0][1:] for item in batch_data] + target_times = [item[1][1:] for item in batch_data] + + input_events = torch.stack(input_events) + input_tims = torch.stack(input_tims) + target_events = torch.stack(target_events) + target_times = torch.stack(target_times) + + return input_events, input_tims, target_events, target_times \ No newline at end of file