minivess_2dunet.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adapted from MONAI tutorial:<br>\n",
    "https://github.com/Project-MONAI/tutorials/blob/master/3d_segmentation/spleen_segmentation_3d.ipynb<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Install libraries.\n",
    "\"\"\"!conda install -c conda-forge tqdm -y\n",
    "!conda install -c conda-forge tensorboard -y\n",
    "!pip uninstall monai -y\n",
    "!pip install tensorboard\n",
    "!pip install 'monai[all]'\n",
    "!python -c \"import matplotlib\" || pip install -q matplotlib\n",
    "!python -c 'import monai; monai.config.print_config()'\n",
    "%matplotlib inline\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import clear_output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import glob\n",
    "import os\n",
    "import tqdm\n",
    "import matplotlib\n",
    "#matplotlib.use('TkAgg')\n",
    "import matplotlib.pyplot as plt\n",
    "import nibabel as nib\n",
    "import cv2\n",
    "import random\n",
    "import math\n",
    "from skimage import metrics\n",
    "\n",
    "random.seed(42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torchvision\n",
    "import torchvision.transforms as transforms\n",
    "import torch\n",
    "from torch.nn import BCELoss, BCEWithLogitsLoss, CrossEntropyLoss\n",
    "from torch.utils.tensorboard import SummaryWriter\n",
    "from torch.optim.lr_scheduler import ReduceLROnPlateau"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "import monai\n",
    "from monai.networks.layers import Norm\n",
    "from monai.data import create_test_image_2d, list_data_collate, decollate_batch\n",
    "from monai.visualize import plot_2d_or_3d_image\n",
    "\n",
    "from monai.utils import first, set_determinism\n",
    "from monai.transforms import (\n",
    "    Activations,\n",
    "    AddChanneld,\n",
    "    AsDiscrete,\n",
    "    EnsureChannelFirstd,\n",
    "    EnsureType,\n",
    "    CenterSpatialCropd,\n",
    "    Compose,\n",
    "    CropForegroundd,\n",
    "    LoadImaged,\n",
    "    Orientationd,\n",
    "    RandAffined,\n",
    "    RandCropByPosNegLabeld,\n",
    "    RemoveSmallObjects,\n",
    "    ScaleIntensityd,\n",
    "    EnsureTyped,\n",
    ")\n",
    "from monai.networks.nets import UNet\n",
    "from monai.networks.layers import Norm\n",
    "from monai.metrics import DiceMetric\n",
    "from monai.losses import DiceLoss\n",
    "from monai.inferers import sliding_window_inference\n",
    "from monai.data import CacheDataset, DataLoader, Dataset, decollate_batch\n",
    "from monai.config import print_config\n",
    "set_determinism(seed=42) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "56 train_files, 7 val_files, 7 test_files\n"
     ]
    }
   ],
   "source": [
    "# Set directories.\n",
    "\n",
    "root_dir = '/home/charissa/minivess/data'\n",
    "\n",
    "images = glob.glob('/home/charissa/minivess/0publish/raw/'+'*.nii.gz')\n",
    "labels = glob.glob('/home/charissa/minivess/0publish/seg/'+'*.nii.gz')\n",
    "\n",
    "images = sorted(images)\n",
    "labels = sorted(labels)\n",
    "\n",
    "data_dicts = [\n",
    "    {'image': image_name, 'label': label_name}\n",
    "    for image_name, label_name in zip(images, labels)\n",
    "]\n",
    "\n",
    "# Split data for training and testing.\n",
    "random.Random(42).shuffle(data_dicts)\n",
    "split_train = int(len(images)*.8)\n",
    "split_val_test = int(len(images)*.1)\n",
    "assert(split_train+split_val_test*2 == len(images))\n",
    "\n",
    "train_files = data_dicts[:split_train]\n",
    "val_files = data_dicts[split_train:split_train+split_val_test]\n",
    "test_files = data_dicts[split_train+split_val_test:]\n",
    "\n",
    "print('{} train_files, {} val_files, {} test_files'.format(len(train_files), len(val_files), len(test_files)))\n",
    "assert((len(train_files)+len(test_files)+len(val_files)) == len(images))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdin",
     "output_type": "stream",
     "text": [
      "Input date, which will be used to save relevant plots/images for this run: (yyyymmdd): 20230127\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20230127_best_metric_model.pth\n"
     ]
    }
   ],
   "source": [
    "daten = input('Input date, which will be used to save relevant plots/images for this run: (yyyymmdd):')\n",
    "best_model_name = daten +'_best_metric_model.pth'\n",
    "print(best_model_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set transforms for training, validation, and saving validation outputs (train_transform2).\n",
    "train_transforms = Compose(\n",
    "    [\n",
    "        LoadImaged(keys=[\"image\", \"label\"]),\n",
    "        EnsureChannelFirstd(keys=[\"image\", \"label\"]),\n",
    "        ScaleIntensityd(keys=['image', 'label'], minv=0.0, maxv=1.0),\n",
    "        CropForegroundd(keys=[\"image\", \"label\"], source_key=\"image\"),\n",
    "        RandCropByPosNegLabeld(\n",
    "            keys=[\"image\", \"label\"],\n",
    "            label_key=\"label\",\n",
    "            spatial_size= (96,96,1), #(160, 160),\n",
    "            pos=3,\n",
    "            neg=1,\n",
    "            num_samples=4,\n",
    "            image_key=\"image\",\n",
    "            image_threshold=0,\n",
    "        ),\n",
    "          RandAffined(\n",
    "             keys=['image', 'label'],\n",
    "             mode=('bilinear', 'nearest'),\n",
    "             prob=1.0, spatial_size=(96,96,1),\n",
    "             rotate_range=(0, 0, np.pi/15),\n",
    "             scale_range=(0.1, 0.1)),\n",
    "        EnsureTyped(keys=[\"image\", \"label\"], data_type='tensor'),\n",
    "    ]\n",
    ")\n",
    "\n",
    "\n",
    "val_transforms = Compose(\n",
    "    [\n",
    "        LoadImaged(keys=[\"image\", \"label\"]),\n",
    "        EnsureChannelFirstd(keys=[\"image\", \"label\"]),\n",
    "        ScaleIntensityd(keys=['image', 'label'], minv=0.0, maxv=1.0),\n",
    "        CenterSpatialCropd(keys=[\"image\", \"label\"], roi_size=[-1,-1,1]),\n",
    "        EnsureTyped(keys=[\"image\", \"label\"], data_type='tensor'),\n",
    "    ]\n",
    ")\n",
    "\n",
    "test_transforms = Compose(\n",
    "    [\n",
    "        LoadImaged(keys=[\"image\", \"label\"]),\n",
    "        EnsureChannelFirstd(keys=[\"image\", \"label\"]),\n",
    "        ScaleIntensityd(keys=['image', 'label'], minv=0.0, maxv=1.0),\n",
    "        CenterSpatialCropd(keys=[\"image\", \"label\"], roi_size=[-1,-1,1]),\n",
    "        EnsureTyped(keys=[\"image\", \"label\"], data_type='tensor'),\n",
    "    ]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "image shape: torch.Size([4, 1, 96, 96, 1]), label shape: torch.Size([4, 1, 96, 96, 1])\n",
      "label values: [0 1]\n",
      "image max: 0.4471955895423889, image min: 0.006625339854508638\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAr8AAAFfCAYAAABHtaTxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAABtTElEQVR4nO29e7Bd51nm+XyWYju2Y1mSZVn3iyXL9zixiWMgTBIDCekEwgxk0t1QgQ6V1HRz7e7h0lVD90zNTHVqqCZMzRSZQGhgYAghYSYQCKm049ywI9mOZMu63++SLV9yg9iWteaPc/bn3/d6v9vHtrSlc9bzq0rlO3uv9d3X0vJ+nvW+pes6GWOMMcYY0wcuONcdMMYYY4wxZlz44dcYY4wxxvQGP/waY4wxxpje4IdfY4wxxhjTG/zwa4wxxhhjeoMffo0xxhhjTG/ww68555RStpRS3nyu+2GMMeblUUrZX0r5wSkc15VS1rzMNl72ucaQ2ee6A8Z0XXfjue6DMcYYY/qBf/k1xhhjjDG9wQ+/5pwzkMtKKf+hlPIXpZQ/KaV8q5SyuZRybSnlN0opj5ZSDpVSfhjn/WwpZdvksXtLKR8M9f5qKeVYKeVoKeXnKJmVUi4qpfxWKeVgKeVEKeUjpZRXj3vsxhgzkyilvKGUcl8p5anJ++//UUq5MBz2jsl79slSyv9WSrkA5/+Lyfv6k6WUz5VSVox5CKYH+OHXnG+8S9L/LWmupI2SPqeJfbpE0v8k6f/CsY9KeqekyyX9rKTfLqW8XpJKKW+X9K8l/aCkNZLeHNr5j5KulXTr5PdLJP3mWRiPMcb0ieck/YqkKyXdKekuSf8yHPPjkm6X9HpJPybpX0hSKeXHJP07Sf+1pAWSviLpz8bSa9MrStd157oPpueUUvZL+jlJ3y/p+7qu+6HJz9+liRvfnK7rniulvEbSNyXN7bruqSH1/H+S7um67ndKKX8g6UTXdb8x+d0aSbskrZW0R9K3Jd3Sdd2eye/vlPT/dF236myO1RhjZiKD+3jXdf8lfP7Lkv6rrut+fPLvTtKPdF33d5N//0tJ/03XdXeVUj4r6ZNd131s8rsLNHGvvr7rugOT567tum73uMZlZib+5decb5xA+R8lney67jn8LUmXSVIp5UdKKV8rpTxRSnlK0js08WuDJC2WdAh1sbxA0iWSHpyU5p6S9HeTnxtjjHmZTFrVPlNKOV5K+aak/1XP35cH8H58QBP3a0laIel3cF9+QlLRhDJnzBnDD79mWlJKuUjSpyT9lqSFXdddIelvNXGjlKRjkpbilGUon9TEg/SNXdddMfm/OV3XXXb2e26MMTOa35W0XRO/0F6uCRtDCcfwfrxc0tHJ8iFJH8R9+Yqu617ddd29Z73Xplf44ddMVy6UdJGkxySdKqX8iKQfxvefkPSzpZTrSymXSPofBl90XXda0u9pwiN8lSSVUpaUUt42tt4bY8zMZGBP+3Yp5TpJ/92QY/77UsrcUsoySb8k6c8nP/+IpN8opdwoSaWUOaWUnxxHp02/8MOvmZZ0XfctSb+oiYfcJyX9M0l/he8/K+l/l3SPpN2Svjb51dOT//9rg88npbn/ImndWDpvjDEzl3+rifvxtzTxI8OfDznm05IelLRJ0t9I+pgkdV33/0r6kKSPT96XH5H0I2e/y6Zv+IU30wtKKddr4kZ6Udd1p851f4wxxhhzbvAvv2bGUkr58cl4vnM18WvCX/vB1xhjjOk3fvg1M5kPaiIW8B5NxJ4c5j0zxhhjTI+w7cEYY4wxxvSGV/TLbynl7aWUHaWU3aWUXz9TnTLGGHPm8T3bGGNewS+/pZRZknZK+iFJhyXdL+mfdl23NTvnggsu6C64YOJ5e9asWfXz5557rjnu9OnTtXzRRRcNrevpp5+u5VKeDyE4qH9YXRmcA9Y1lXNfCqP6OazN7Jh4XNYGx8X5PnUqt73yuFe96lVDP2cbUrsWXEu2f+GFz6d25/nf/e530/6zffaZY+fxs2fPbup69tlna5lzyT3F/sZ9yL8vvvjiWuZ4s70zZ86cpq5/+Id/GNovMmrvXXLJJbX8j//4j7XMcfGcUdc15ymOeQDHy7mPeyfbbyxzvuPYz+Q1xnGx3my/cE/GcWXz0nVdjFc6bXg59+wLy0Xdxbp0TD00xpgzx3f1HT3TPT30nj172IdT5A2Sdnddt1eSSikf10SO7lEPv7rssok8AvPnz6+fP/nkk81xfFBYuXJlc/6APXv21DL/0Xv1q1+d1pX948h/6PiPIx9ypPwfdz4Yjnro4MMcH2bY/jPPPFPL/Mc5Pth9+9vfHtpntsEHjSuuuKKWT548mfZxsD6StHTp0qGfx77s3bt3aL/Y/vLly4f2d+fOnU1dHPPChQtr+fHHH69lPvyxL/PmzWvqOn78eC1zvtesWVPLTzzxxNC+S9JTTz1Vy9yH+/btq2WuHffBW97ylqaujRs31vKJE88nseP5HDvHKEk33HBDLT/88MO1zP3Ovc6Hubgnr7zy+WRLnFeybt3zUd+4X9h3Kd9vbJ/76NFHH23Oj3M+DF732X9sSO0e53X0rW99q5a5X5Ytez7OfrwHcV/MIGvYS75nX6xLdUe5a0zdM8aYM8f67u70u1die1iiNkXhYQ1JQVhK+UAp5YFSygNn+pdUY4wxU+Yl37Of1dPxa2OMmfa8kl9+p0TXdR+V9FFpwvYw+EXr6NGj9ZgoOfKXnu985zu1zF+ZMuk2yqpZXfw8swRQ+pXaX5MIf7GjjB9/Ic1+Yeb5r3nNa2qZv7Kxj1I7Tv4ylf1yzLHwV3ep/dXr6quvruWtW5//QSj7pVpqf3GlxM1f9Y4cOTK0j1xTqf21/Rvf+EYtX3fddbW8a9euWuY88Bc+qZ1//trJc1jmr9tS+0sqf92+/PLLa5m/0HLvbNmypamL88+1ZL84F/zVWmp/iWRdnMvsF8pLL21la46L+5B7j+PiOi5Z0j4r8Tj+iszri78cj/oVld9RKeB+Yxvc31K79zhHXFf299Ch558D/R/mz8N79uVl3oz52dsYYwa8kl9+j6jNz7108jNjjDHnH75nG2OMXtnD7/2S1pZSVpVSLpT0XiG9rDHGmPMK37ONMUavwPbQdd2pUsrPS/qcpFmS/qDrui0vck6VF7OoAFL+Eg2P48sqfPktWhP4QhBtD9nLazw+Rprg39/85jdrmVI9+z4q8gTlZtoOshfDrrrqqqau7KWpLCoAj4+2B7ZDOwpZvHhxLUe5+fDhw7VMGZ8vjK1YsaKWuQ4sS+280vZA2wTHRUk/StdcL87fwYMHa5mSfrSWcC5Z91T2UbS8cM64R/gCGD+P1wQtFQsWLNAwaI3g3MeIGlx/jpkvjLE92klomZBya00WjSXaXLIX2I4dOza0j2wvXp8cJ+8brJfXSmbRyfo13V98ezn3bGOMmYm8Is9v13V/K+lvz1BfjDHGnEV8zzbGGKc3NsYYY4wxPcIPv8YYY4wxpjec9VBnTWOzZ1dPKL2DMaA/PY30Z2YZx1ieagYpfs4QT2wj+lFXrVo1tE16LeltpU9Vav3E9MbSbxiTCAyI/uHM28t+McQT/ZgxoD990plXc1TmOYY6Yzvs47Zt24b2kfMV4bzG5AjD2oveVu4FjjHLXhZDgtHryro5LtabZZSTWg9wFpqO/YqeYe4L7iO2wz6yjdgX7rcsk1mWuS1eX/Tgsh3OK/sVk9CwHZZ5ftYX+ppfzvnsb+xX5us2xhgz/fEvv8YYY4wxpjf44dcYY4wxxvSGsdoepOflSMrrUVZk6CxC+ZJ2CGadYoikUefT3kDplpaAmOGNYcAokVMGZ/axRYsWNefTOsCsV7RdsE0es3///nQslMgp1a9cubKWKa9nmepivSzTwhGtJGyf33Fd41wOGCWjs/0MhkaLYbR4PjOmMaQXLRjRWsGxZJaEzEJw4MCBpi62wzViH2lfibYHSvRz5sypZc4fz8/WUcotLDFD3rDj4zpy79OmkvU9ZtHj9ZbZE7iutEvFfZidT3hOZp2S2vkf7MloqzHGGDM98S+/xhhjjDGmN/jh1xhjjDHG9Iax2h6ee+65KlPzbWzKuFL7pjWl5OyNeR4fJV5mgnvsscdqOYsEQIk1vglP2ZNWBcq67OPjjz/enE/JlcdRSuYb7OxXHBfHn2XAIvycfZda6T7aEIb1fVQmtIULF9Yyx0JLAd+sHxWhgW2yXmbB27p1ay1zT0ntunAtOHeMLBIjT9C6wDlipA5mjqPlJkYK4XrTMsM9xs+jNYXtU/pnO8xYNiqKBec4y8hHOwnXm1YaqV3XLLIKz4/XBL/LIjHwGF6ftE5J7T7MsuVxXLzuYsQZR3gwxpiZi3/5NcYYY4wxvcEPv8YYY4wxpjeM1fZQShkqy8ekCzwmk4hpA6D8Gd+SpxRNKTULiM+30aMd4xvf+EYtUyKm3EvpNVoVokw7rP0sikOUYTN7QJZA4fLLL6/lKMlnVocsWsPatWub4zive/bsqWWuHWVoRlvgnEqtpM/xsw0mzODxMdoD7RWsa+7cubVMawPnKMI9SXsBbROMBhItAFyLbL3YL9ocYvuc//Xr1w89hmOn1B/hmGkLYhv79u2r5b179zbn0/bBtbj66qtrmdaIaHPh+bw+uEdZL8cVbS5cf9oYuHd5THYNRQb9GnWMMcaY6YN/+TXGGGOMMb3BD7/GGGOMMaY3jNX20HVdlRAp90Z7AOXyefPm1TLlTwbkpxQaIxFk8ielayaToNw8ymrA49gmj2Ebsb4swH4mj0eJl3XF7wbwzXrK/pTqpVaWprWE/aIkHZMhMLkB7QkcC89nOUYiYN3sc7ZeZFR0Ds4Rk4/s3LlzaBuSdM0119TyoUOHapnrwn3AfRzrYvs8n5E3smgkUmvHoY2BFhJ+Psp+wzaZ/IPncB143cQoFDExyADaWTj2GGmE4+KYub85X1zTeH1miXN4Ducis59EuKeNMcZMf3xXN8YYY4wxvcEPv8YYY4wxpjecs2gPlFKjjJ5Jm4y+wAgRPD5Ge8gSAvDN+sy2wAQAsS+0U/CNeY4lRragLJvJ4HxLnhIzo1ZIrV2AdbGPlJQPHjxYyzEqAm0QlN5pD6Dsf/jw4eZ89pMJKE6cOFHLWRSLKMmzz5wXWjiyxAhTjYjx9a9/vZZvueWWWmZUg9gm55jrSHmf7UfbQiadZ5aCGGmE9TGiBvdYFtVgVBQM7mNeR1ldkWxPZ/YVrm88h7aVzIbAdYzXVxbphLAvtFRxTqQ2CogTXhhjzMzCv/waY4wxxpje4IdfY4wxxhjTG/zwa4wxxhhjesNYPb+nT5+uHlx69xiuSWpDNtEHGMNHDWCIpuitpD+U/kaGXMpCREWfbRaGjP1n5reYvYx+Xnpd2Ud6a+m1jCHc6FekT5d95LjYXgwvloWVol966dKltRy90OvWravl++67b2g7bINzFD2/WWYyzgX92vTl0m8stSHRuC84lzyGPtd4DvvJ0G4ZMVsc+08/LL28K1asqOUsfJ3Urgt91VwjzmP0H/M6YL94TfCYUV7mLDsi14jryP7GfsY9PgzWy/FK0ubNm4eew/7zfjIqrCHnImZENMYYM73xL7/GGGOMMaY3+OHXGGOMMcb0hrHaHqTnpU1aBSjbS7lMyVBMlEgZ9oy2BamVNhcuXFjLR44cGVoX2+a5w+oewHBPDKUUpW/KpxwL5W72hdaGKLVzzJkdg22QaDOhFM3zKX2zDVoNpDZLGGX0bLw8Joa+4vhpO+A+oFWBsnuUp/kdx0L7DPdElN25ruwnJXFaEHg8pf44Fs4x9xT7y/WVWqsE+5mFjePnHLvU2lEYLozryP7ymDhHWXixeE0P4Nhj37KxZPuA2fliXwjPzywcXDupXVeHOjPGmJmFf/k1xhhjjDG9wQ+/xhhjjDGmN4zd9jCQLSnjRymUf1PKpvRLiZUSaZQvKVlm9gZKqaMiLPCcZcuW1XIWnSL2hbI2y8wslkWBiFEs2BfK+LQXUN5mX5i9SmrlZkrcWVYyyv6xTa5RlnmO5Zili9A2wj3BKBqMFhDnm3uH53ONjh49WsvRDsL+UwZnxAHOBduPc8e/s0xonMeYoYz7kvNCaw73ISX9aC3hcZwjtsHxcr5iJrTMzjF//vxa5rhiFIsYPWIYrJfjjVFLCPvP64hjyTIzxuOy68AYY8z0xHd1Y4wxxhjTG/zwa4wxxhhjesPYbQ8DCZFSPSXG+B1l2Ux+pdwcEzAwMgElbkqpWYQHyttS+2Y638bPJNYs2kJsM4t+wGQQEVoVOH7aJthffh6lZ8q6nFd+nknqUm7PYF0sUxJnUo94fpZkgvuDtoloe4h2mmH1si6ug9SOk/uQY2GCDybJ2L17d1MXz2GbjOjBCCRz5sxpzue+yKI6MMIC1z7K9myfY3zLW95Sy3/zN39Ty9wv8VplX9jHRx99dOg5U00mwQQWhw8fHtoXXneStGjRolrmtcN6aXVgG4x0IbXX92BdY/QXY4wx0xP/8muMMcYYY3qDH36NMcYYY0xvGKvtoZRSZeosMYOURwbgcZnVICaDoExKiZx1US5nOSZzYEKC7BzKwFHiZZs8LnubnXMUoyJQoudb76yXySsoj8f5zhJTUC6nPL1gwYLmfNpJKKlfddVVtbx///5aph0kJg5hJAeuK8fIsVBqj+NavHjx0D4Sjpe2BamdS5bXrl1by5s3bx56frQHZNI7z2E5RhrJojJwjbgujF6wevXqpq7t27fXcrZGbIPzGi0v/I7rQngdcB5j/7nG3Pu0H3Fco64J7tfMYkTLTbQC0XYy2K9OdmGMMTMD//JrjDHGGGN6gx9+jTHGGGNMbxi77WEgjfLt9yilZrI4pWC+zc0ID5TdpVbupuSaRQlgG1F6Zr+yaBGUT6MMnMnaWYSGLDGD1M7f1VdfXcubNm0a2heW49v/bD9LLME1ilEo+NY835JnmXNHawpld6mV99lnjpHSN9ebESkk6bbbbhvaF0ab4P6IMjojFpBsf7KNKO9nlhu2z3qjxM52KMkfO3aslrknuPeYREXKo1hwXrk/OcexX7x2aFnhGEclicja5HzRAjEqOseotRzAsWdRXqTWgjGIyDEqeosxxpjpg3/5NcYYY4wxvcEPv8YYY4wxpjf44dcYY4wxxvSGsXp+Z8+e3WS0GhBDVGWhv+iZpVeUPsCY5Yu+Ph5H/x4/ZxsxgxSPo6eSflx6KBluSWp9oAzDxdBTbJ910fcotV5PhmRbuHBhLTPUGccbvcz0iq5ataqWv/71r9cyPafMpCW1c8w541rQn8m+nDx5sqmLnmN6g7lv6MXNMsJJrdeVe4pjYXgw+orjOZxL+rKzTIPR752F1uP53B9xvel/preVHlauKz+Pnnr2jT7pzCfL9b3iiiua7zguriXHwnWJbXD87Evmn+b+jnC/xNBlw/qbhVSMdQ32ZJYx0BhjzPTCv/waY4wxxpje4IdfY4wxxhjTG8Ye6mxgHWA2qVEZ3ih5ZlYFhtraunVrUxfrpixMuZXhlhhGivJ2rItSLmVlyvMxFNPx48eHnsOxREvCsGOkVjqnVWLlypVD+895jJIwx8/5y+wcMXRWtkb8nOPlvESrAcOocS4OHz5cy5yLURnGaAdhXdx7zNbGMUrtPF1zzTW1zDBYtOVwvmKosyw7GO0N7GPcO9yvWcg9lmkvWL58eVMXJX3ON9ug7YB7PWbk497jGDlHvKZi+DiOn3uEc0n7EOeIfYywrrgWwz7n2GObgzW27cEYY2YG/uXXGGOMMcb0Bj/8GmOMMcaY3jBW28OpU6fqW/PZ29hSKzlSuqcUmmVCi2RvtlOipURKSXlUZirKv7QNUMbN3p6P33FcmbUi2gNooaD0vmHDhlqmJL1kyZKhbUgvzIw2IBtLPJ8WDI6F0n02XkZRkFr5mfYGzjf3zqFDh2o5WgsYESTLMMcIBTHCAueV9gDuKZ6TRRKQ2rVgtAxK74899lgtx0ghbJ/zyrnkHLPeKPszIgfrGib1S21EDFqMpDZSCS0/nAu2H6Ox8Dj2hVFH9u7dW8tc+2gNWbNmTS1zvniv4Fi49jEKDf8e7LHYd2OMMdMT//JrjDHGGGN6gx9+jTHGGGNMbxir7aHruipb8g10Bu2XWon4ta99bS0/9NBDtUzJc8eOHbU8SpqkZMs3xSmlMkIArQJSK5dTOr7sssuGfk5JWmrlX0rUWWIIyv7xLXkeR+mY8j7bz5JySG2Ug8x2QAtGPJ8SP9eOb+bzHLZBa4LUzsvOnTtrmevF9ljv/Pnzm7oyawr7xTUdZZ+hJYD959oxAQTbk9r14n7nWLiPos0lJlwZwKgKXGNaKGhhkF6Y9GIYnEvWFS0ybJ9rxAgPPJ/HS+24OBdZFAsmh6G1QWr3Dtee4+Xccx1jZJeDBw/W8iCxRlxTY4wx0xP/8muMMcYYY3qDH36NMcYYY0xvGKvt4fTp01VepHw6Sm7es2dPLVPippSZSe1SK6syQkIWVYGfx7f/aVtg/ymJU9aNUjXbYd1Z5AmON0YP4N9ZVArW9Y1vfCM9nhYSzl9mdYjnU6Kn9E2JmtaSAwcO1PLll1/e1JVFdaA9gf2lVYCWFamVsmkHoT2A8vzcuXOb85lYgxI715X957ywPaldr9tuu62WGZ2D51N2l9p9SZsK54J9ifYCQnsAoyrwWsvsGNH2wDnOLBDsV5wXJmXhvHItaE/gtcZ+SdKRI0eG1sX7A+cxi5Qhtes9sG2MilBjjDFm+uBffo0xxhhjTG940YffUsqyUso9pZStpZQtpZRfmvx8Xinl86WUXZP/P/fF6jLGGHN28T3bGGNGM5Vffk9J+jdd190g6Y2S/lUp5QZJvy7p7q7r1kq6e/JvY4wx5xbfs40xZgQv6vntuu6YpGOT5W+VUrZJWiLpxyS9efKwP5L0RUm/NqquCy64oHo66VuM0N955ZVX1jI9hfSw0rsXPb/8m/5O+j6nkr1Laj1/MRTVAPqKo4+QddMTSS9xFuIpZuliXQxJNpXzOQ9Sm+mKY6RXkp/HUGccJ+ebftRjx44N7VcM4cawWBwj/bD0ozIcHUPexePYf5bZl7h3CNvPMtrRc8q5i+1s3LhxaJs8JoYjoxc6C43H8IFbtmyp5XXr1jV1sW/04zLbHuco889KrYeXdfHaYUY77gOpHT/HTP8x9zd92PH6Yj/pZWaIMpbZXvSe87vB9TFdPL9n8p5tjDEzkZf0wlspZaWk10laL2nh5E1Wko5LWpic8wFJH5gsv+yOGmOMeWm80nv2xbpk2CHGGDOtmfILb6WUyyR9StIvd13XvPLdTfzE1g07r+u6j3Zdd3vXdbf74dcYY8bDmbhnv0oXDTvEGGOmNVP65beU8ipN3ET/tOu6v5z8+EQpZVHXdcdKKYskPZrXMMGsWbNqqKIsq5rUZmKbN29eLVOWpSzJuhgKSWrtCWwny5BGGZay+6g+sy98wI8yeswiNYA2BPZ327ZttTwqPBn7OchGJbVz8fWvf72Wo+XkpptuquXNmzcPrTcLDSe1NhVCS0MWnitaS2hz4Ryzfa7D3r17azlmi2NmM1oguEYMtRXDyVEiZ5+zkHtZmLrYf84f14v9jZI+62NfWC/XLmtbauecVhF+TpsFP4+hylg398GCBQtqmfMSMzAy7B2vddoQ2AZtMbQ+xeNoZaJVg+HRuA9Yr9TO8SBkHy1N5ztn6p5tjDEzkalEeyiSPiZpW9d1/wlf/ZWk902W3yfp02e+e8YYY14KvmcbY8xopvLL7/dJ+mlJm0spmyY/+3eS/qOkT5RS3i/pgKT3nJUeGmOMeSn4nm2MMSOYSrSHr0rKzLp3vZTGTp06VeVISr/xDXLKv5SeGe2BkiUlyhUrVjR1UbqnbYIycmZViNms+NY53yzPLAjR5kAbAuumJYB9ZF9GRXvgHGVRAbIMYZK0c+fOWuYYKb1nETGkNgPX/Pnza5nSN8/nuGKWruuvv76WH3zwwVrm3FE659ijzYRzxmxzjLRB6Zx9l1o7CmVxtsO55BxzT8TvOOYsUgnXQcozqXEtaEkYlfkuy4jI9eIcc+651rHuLCIHIzxE2wStKty7vD64jrRgxPWmhYTzwigcbINlXnfx/MG4plG0hzN2zzbGmJmIM7wZY4wxxpje4IdfY4wxxhjTG15SnN9XygUXXFAlV8rglB+lVgrdvXt3LWeyMt8mP3ToUFMXv+Ob/ayXEu2aNWtqmRESpFY637dv39B6KU/Ht8Mp31I+5jmU5H/gB35gaH+lPPEAx0uJ+bbbbqvlrVu3pv2i3M7zo9xNKJHTssI15tqNqpfRG26//fahn7O/lKKjdM26aS9gX2ibYN9jP9kObQeU+rPkExFK+uxXtP8QXiM8h33MkrhEqwGhPYKRF/g590ScI84lLQlc+1EJbbj3aa9gn+fOfT4LL6/7aENgdBGuC+cr9n/YMbGdQR9j1AxjjDHTE//ya4wxxhhjeoMffo0xxhhjTG8Yq+2hlFLlRAaxp4wstZYCypGUXykDU+6Mge/ZDmVdQjmT0nWU5GlJoKzKfvHt95h0gRIvy4xwcPjw4VpmAoJRfckkbvZr165dtRwlXtoIKFFfccUVQ8+JkQh4Duc/iyqQJayQWkvAhg0bhrbP8zn2mKSCbfL8zJIwKloEx8LPaTPhvNCCILWWgMy2wD0Vo2Bke4znZIk4YqSQqSTmYL9G1ZXZHtjft73tbbXMNZXa/ZKtV4wuMqzvsW8sM2oJr8mjR4/W8q233trUtWnTphf0K+4vY4wx0xP/8muMMcYYY3qDH36NMcYYY0xvGKvtgVD2j2+DU16kfMlzKLdSko6SdiaLU1alNeKJJ54YeoyUS7RZAoQlS5Y059OewXHRwkAplxaAmHCD4+Q5lNsp8XJcUd6ndM61yJJ/cBxSa8mgVYIyONtgvXG92E6WaIFWCR4To4Zk7XBPcB+xj1K7r1g354j10r7DPRXr4rxs375dw4hJMtg3Xgfz5s2rZVoNoj2B7N+/v5YZ/SDu9wGMiBHnKGuH+3XPnj1pG/yb68oy+8h1iFEsli5dWsvZtcoy7we0WknS6tWra3kw5sw2ZYwxZnrhX36NMcYYY0xv8MOvMcYYY4zpDX74NcYYY4wxvWHsoc4GnrssA9PguGFlhiGjt3Tt2rW1fPDgwaaukydP1nIMqzWAPkD6QWN4MfaTXkv6UbNsbfE79oW+08yfGMM6sZ88jpm16LukBzV6Yzmv9AxzvDwmzmPmFeUaM0vXKF81fbOLFi2q5RtvvLGWN2/eXMuc42XLljV1RW/yAI6F/Yrj4t7jet9yyy21/PnPf76W6QWOa0//N9eF7bO96GfNsvXR28pxrVq1qpYffvjhpi56jjN/fRZabpSXmPA49jGGLeMeobeXbfJ89iuGHuP6cR/z2mHWRbYRPfXsy6DPnDdjjDHTF//ya4wxxhhjeoMffo0xxhhjTG8Yq+3h2WefrWHBKItGKZRyJMuUOSmFMhtTzD5GiZ+SKS0IlDgpHcewTrRasB3aIxj6KkqptCTEuofBemPGMGanolWA0j8lXma5iqHlOMcMF8XwT5R8Y0Y+hlTjHLOPtGZw7Gwvtsm55PiZ4YztRVmaMjj3C/cR14jhrSTp2LFjtcw9ynnNQritWbOmqYt9pgWFe2/lypXK4BrxfNbLOWZGPx4jScuXL69lWkMYXo3XRxYiUGpDuPGaYF+4dhxvPIcsWLBgaF2cY4YYlFo7DdclC2fHNiLsZ7QcGWOMmd74l19jjDHGGNMb/PBrjDHGGGN6wznL8DZK9qe0SsmTNgBKsZTx4xv7tD0wKgM/p1zKt+qvvvrqpq7rrruulimL8riFCxfW8pEjR5rzKf9S3me/OHZKyhGOn3L1iRMnapmZuTiPlIHj31mGtyz7lpRnrmOEg0OHDtUy5X32UWr3Ba0au3fvrmXaG2iN4DpKeRSPbO5ihAb2hVYHzhHXjpYZjleS7rzzzqFtcr8QWmykdl9wjWmB4DrwmCjbc79y7Wl1YHtZBkGpHT+vT14fhw8fruU4rssuu6yWOX9cC/aFaxKjPXDvsc/ZemfjjefEdowx/eJzRze9ovPftvjWM9IPc+bwL7/GGGOMMaY3+OHXGGOMMcb0hrEnuRjIjnzLO0ZFiOcMoCyaBduPdgpK9JSrKYOzL5Q4o+1g27ZttczkBEyAwH499NBDzfkcJ/vJfvENdPYxJu9gYg9GMmAbHMv8+fNref/+/U1d7DOtGhw/+xvnmPYCrhEtCVdeeWUtM8JDnCPWRVl68eLFtUxJnBaIaFsgTJhBqwXnaJTtgfuIdgzaZPh5TJDCsXCNf/qnf7qWaVWI0Kbz2c9+tpY5f+wL+x6jSNAewTFzvTkW9j1G1OB+4z7k2pNR0R54TfGayOw7PF5qx0zbAy0NrJdrH608tEcYY/rHK7U6ZHXZAnF+4F9+jTHGGGNMb/DDrzHGGGOM6Q1jtT1ccMEFVQ5moH2+DS5JJ0+erOX4FvYAysCUWOOb7ZT0+WY7k0Fs3759aBvxzXRKtG9605tq+ZZbbqnlr371q7W8bt265vz169cP7VeWqIHRC2JiCUrklGgzqZ5v5dOCILVSdJZIhHMX+8J54vlcR1o4mJSENgmpjURAuZ3nc+4og0epmnuBCSsYYYB10b4itZEUmFAhi4jBMhOaSG2iidtuu62WaTvgeDn3Uiv933jjjbW8YsWKWr7vvvtq+eGHH65lroPU7gvuaV5HWRSGaHnJ7AFZtIiY1IJzxvbZR46R80Urj9RadrgveE1xH9PmEe0YtK3EpDBm+pNJ2pakzdnGe+/8wL/8GmOMMcaY3uCHX2OMMcYY0xvGnuRiIDtS4o1B5LNA9qtXr65lyuC0B0SJkjInJVdaLVjXtddeW8t8Kz5+x2gLlMR5TpTRGbGA48+iQLDv0R5AWZmJB2gHoaScWRukdl4YMSBLNhKTZLBuytU8jpJ2FpVAauV+tn/gwIGh9bIcIxFwXXkc+0vZnhEOpNbeQYmca8y5pFROW40k3XrrrbXM6AtcV+6DuPfYPvvMupgI5JprrqnlaKGgJYOWnS9/+cu1zPXiPEYbEuvifuM+5HpHywxtDJxLjmvHjh21zP01ys7BeWW9vFY4j6MSZgzsL6OiiZjzkzP5xr6Z2ZzrveKIEOPFv/waY4wxxpje4IdfY4wxxhjTG/zwa4wxxhhjesM5C3XGjGXRq0n/Hf2Ghw4dquXrrruululHZXgqqfX+0dtK7yE9u/Qq3n777U1d/I7Z3uiDpC+Qx0itxzAL4Ub/Mn2H0ZPIvjCMGeeCbWSeWyn3X7NerkNcryxEFX22rJfe2phFj3AsDJfF9tkv9ndUn/k5Q3qxPakNT0bfLH3l9HGzHP3e9PN+7Wtfq+Xf/M3frOWf//mfH9qeJL3rXe+q5Q9+8IO1/Id/+Ie1fMcdd9Ry5m+X2jFzL7zxjW+s5S1bttRylmlPkvbu3VvLXGPuCe6DUXvvrrvuqmXeHzKPN8O/Se31wvB70cs9YFQIs2HZIKPX3Zx/nGvfpjFngriP7QE+8/iXX2OMMcYY0xv88GuMMcYYY3rDWG0Pp06dqrI+JfkoJzJkVBauirIqYSgjqZV/mRWN8vyiRYuGtsdwZlIrpTKD1lvf+tZafu9731vLH//4x5vzH3jggVqmRLtmzZpaZkiv48eP13K0PVDGp6TOemkVoPQbLRecb9opaMHILAxSK12TLPtZFpJKauefa8k2aJXIQn1J7Zjnzp1by5Te2V6Uxynpc84Ynm3hwoW1zHHt3LmzqevDH/5wLTOMGeeFMv7WrVub89/5znfWMq01e/bsqWWGU2NdtDBIrZ2DIdlo1eC6ZPaX+N2oTHAD4hxfddVVtbxhw4ZaZji6LLtgzObI7ICE53Pu2N9ovxkWZtC2h/OTM2l1sNxszkccBu3M419+jTHGGGNMb/DDrzHGGGOM6Q1jtT10XVflY8qilH6lXGKmrEkJkjJujERAmZQSL2Vwtk8LAS0IknTllVfW8k/+5E/WMm0Tf/EXf1HLX/3qV5vz2Te2mVkK5s+fX8u0b0jSnXfeWcuUi9l/yvAsR7k4y5JGGZwRMTgPUitRc71ooeAxcb0J+8boC9lYaGXJImhIrTzP/cJzWK8kLV26tJb37dtXy5wvfs7jf/RHf7Spi2P+/Oc/X8tHjhypZWZli+u9fPnyWqZV4G1ve1stM0LB29/+9lqOUSwefPDBWuZa8jiWacuJkSNoBeI1lWVui1nSeE523dMqwes5WjBoZ6GNgXuKdqlRUUfY50F0DGd4O384W1EdLCmb8x1bIM4M/uXXGGOMMcb0Bj/8GmOMMcaY3jBW24P0vLRJ6TjKspQc+ZY/A/9TphwVBYJy9+rVq4d+vnHjxlpmhAfK2FIbTYByO+uiVYLyutQmuaBEzXFRhh5lNaB0zXaYwIBtZJaL2H/aJngOZeiYSITfPfXUU7UcIzkMY5RVgRELKPVzjbNIE1KeXIFRJDhero/UzvnJkydrmePieGnz2L9/f1PXQw89VMtLliyp5XXr1tXyzTffXMuMtCG1FhRaD5j4hXPx2c9+tpYZNSS2yb3HsbCPtA3EOeIeXbVqVS1zjzDaxKg9kUXRINzr0fbA/UqbDe0sXCNGOeHxUjvmgQUlRjkx48UJLIxpsQXi5eNffo0xxhhjTG/ww68xxhhjjOkNY7c9DCR3Spx8g1tqZXjaCyhrUuKlvB2lUNoAMrmcUSDuuOOOWj527FhTF+Xiu+++u5aZ2OILX/jC0HFI+dvslIKZAIG2hYMHDzZ10SrC8zl+yr2Uq6PVgPYCnk+Zl1aJKMlzXtkmx8+xsF7K0FIrRbNfbDNGQhh2rtSOhbaFG2+8sZYp+8doD4zkwO+yRBy04uzataupi/YAns9x0WYT9w7nj4ktWO/9998/9HNaJqQ2Ogn3Am0PtHzQMhIjHnBf8DhGgbj++utrmUlJJOno0aO1nF0HbIP9jRaKLNJIZjfKEsJI7VoO9reTXBhjzlcyW5DtEMPxL7/GGGOMMaY3+OHXGGOMMcb0Bj/8GmOMMcaY3jB2z+8A+j6jD5B+VnoimfGMmbGYqSl693g+Qz4xdNYb3vCGWmaWLXoopTZrF0OafeITn6hlekjpN5ZaD2uW1Y2+5lGhr+ivpHcxZrgb1l70RZNsLlmmNzTWnY2Lfkx+HvuSZaLj56yL7bHvsR3WxZB19IbGMHk8hx5Q9pk+Y3pQ6SuW2tBhDzzwQC1/5jOfGdoX+mSlNhMa/d/XXnvt0Da4RtwrUuuZpv+X+4jeXh4T66IvnuHCOF+87uL+5HXMteTc8zqi5zeG7GOoNfaLmSEZJo7vGvC6lYaHw3Oos/EzjvBm9kQaafQ+mMo+jOefL6H5Yj+83yfwL7/GGGOMMaY3+OHXGGOMMcb0hrHaHmbPnl2tCww9FcN40ZJAyZPyK0MhrVixopYZGk1qM3hRimamKMrglL5j6CxK1Mxgxf5mMq7UyuKUdSnXU95nSCtK/VIrH2dWA8rzPD6Oi20yBFwW6oxhrKQ2lFQWkox1UR5n6KkIxzUqLNWweqV2zijDsy7K3VFGJww9xjJleErq99xzT3M+JXlaGDgWzv2GDRua87lHOP833HBDLXMf87qJmQbZf9oYfuZnfqaWt27dWsu/93u/V8vRjsF1pRWI1yfXMdoeaPXg+DmvmRUmwr3EMu8P/Jz3hsio74wx/eLlWAWyc84XO4TU7wxx/uXXGGOMMcb0Bj/8GmOMMcaY3jBW20MppUrk3/M931M/j9mwmFmL0jk/pyzKKBAxexnfWqf8yaxulDgZkYF9lFqJlhIz+8+MVVGGZxQLlilRsw0eE+Ve2juyjGm0NzDjF+XleD6zgfGNd859jPZAGZ+ZwSi3sw1aEGKEBtoQMttFFjliVAauLHMd1yhGnuD8c/64xzjHoyJqcP4ovWcRJWirkdr1pqXgIx/5SC0vWLBgaL9oMZJaqwPn8kMf+lAt33bbbbXMayVmi+N3nEtmq3v3u99dy1/84heb89evX1/L3Jecv8yKw+teavcOrTj33ntvLfNa4zzESA4xqoUxxpwJaC84FxaIrM2+WSD8y68xxhhjjOkNU374LaXMKqVsLKV8ZvLvVaWU9aWU3aWUPy+lXPhidRhjjBkPvmcbY8xwXort4ZckbZM00GM/JOm3u677eCnlI5LeL+l3R1Vw6tSpKt9SUqa1QWqjImSRCBjQnzJylHiZtIJyNaV6HkMZeP/+/U1dtAFQ1qUkzb7zjfnYPiNEUK5mm5TR+Ya+1EYPYIQFzhGjAjDQf5yjLOJCTNIxgGsX6+NYGJ2DdhBaHWjTiN+xz5xLts/oHDEiBmXxmCRkAOVu7jWpnX+uHS0Qt9xyy9DjY8IM2hhos2G9tOjEsWTRKmhBoWXnHe94Ry3HaA+HDh2qZc5/PG7YMVxHqbVBrF69upZpQfjCF75Qyxyj1EZ/yO4J3Mccb7Qq0KbEMbJMa0Zmn5HaOR7sySyBzHnMK75nnwvGLQX3Te415w/n2gKRMaovM+UamdIvv6WUpZL+iaTfn/y7SHqrpE9OHvJHkt59FvpnjDHmJeJ7tjHG5EzV9vBhSb8qafDT1nxJT3VdN/ip7LCkJUPOUynlA6WUB0opDzg9qDHGjIUP6wzcs5/V08MOMcaYac2L2h5KKe+U9GjXdQ+WUt78Uhvouu6jkj4qSZdddll30003SWqtApRrpfZN70z+5DF8M5tveUutRM26KD1feumltUyJd/PmzU1db3/722uZyQUeeuihoX2MySCyt+wZLYKfU4aPc0QLBiVxStdMHjIqgQOtBuwj547zHSV9fkc7CCVuRjWghSPWRen6vvvuq2XaBqJ0PiCOkRYUWhJomWFSkrh3uBfYT7bDtVu4cGHaF1pL2M6aNWtqmZI8E2bE7zJ7ACM8MBpJtJZwLMeOHRvaF+5pRnSIdgxaCjguWgTYPq0RUmtNYZKMaEEZwL3O9ZXafcHzM5sN7wEcuyStXLmylgfXxKZNm4b26XzjTN6zLy/zxv6LxfkqBRtzNhllJzifroOZYhOaiuf3+yT9aCnlHZIu1oR/7HckXVFKmT35S8JSSUdG1GGMMWY8+J5tjDEjeFHbQ9d1v9F13dKu61ZKeq+kL3Rd988l3SPpJyYPe5+kT5+1XhpjjJkSvmcbY8xoXkmSi1+T9PFSyv8saaOkj73YCd/97ne1fft2Sa0UGeVmSpuU+yml8q1vyu6LFy9u6mJCAMqktCRQ+r755ptrOUqvGzZsqOUdO3bUMqXbe+65Jz2fb/ZT/qUFhJYAjjfaHmhp4PxlSR8olcfEEvyOcj0jJPDzGJ2D/aRVhOtIGGmD45Ckr3zlK7VMC0YW7YB9ifNNiZ7ncB/wmDgv3Fe0xtDCkEVroGwutevKeWWkDNoRGIFEau0s3C+0pvAYfh7Hxf3KOWN0DY6dUT+i1YDz9/3f//21zOvzs5/9bC2/5S1vac7/4R/+4Vr++7//+1rmdUvLDG0WcVxsk3uHlhlGI+G+5/pK0t69e2t5ENFiVBKVacJLvmf3geks3Zp+kO3Rc22HmM4WiJf08Nt13RclfXGyvFfSG858l4wxxpwJfM82xpgX4gxvxhhjjDGmN/jh1xhjjDHG9IZX4vl9yXRdVz2S9EAyjJPUevkyDyw9u/yc4YukNuMZPYFvetObavlLX/rS0PaYRU1qPZ30tvI4+kFjWChmdWNf2CYzY3Eeos+WfkXOH73EmTc1+jbpZcyysmVh2mI/6T9eu3ZtLdNDyXrp55TataRPNwt9xc+jJ5N9uf3222uZvlG2wTHGv6O/dFj7bC9mS2Pf2CbXm+sYPagMo8Y22Ud6hjmPcb35HcOjcY3ok+XeoxdYasf8yCOPDB0LyzHrIeeF+/X48eO1zFBr119/fS3Ha53h4XhNcl44d/yc10dscxAeLY7dnH3OltdxuvkTjRnGVPbx2fQFT+fryL/8GmOMMcaY3uCHX2OMMcYY0xvKOFMOX3zxxd0gBBQlTloApFYWprzPvlKyZHamKH1T4mZ2p9e//vW1vHPnzlpmGC6GoYrnMxTTjTfeWMu0WcTsZZSyt23bVssMYxbl1wEMTxXbYTmbL0q20bbAv2nB4LpwTaIFgH3mHLHPnEuGLWPfpXYuuBZsg2G82K+4lzkuSvfcI6P2Py0ZWXg1WllYF8ORSbk9gmHEaEGI/aK9gf3P1msQnkt6YUY8XntZ+C6u4/Lly2v51ltvbY5jmMCNGzfWMrPojYJ9YztXXnllLT/88MO1zH0QLTPZGvOa4HrRjhGzMQ6by4ceekjf/va381SJM5DLy7zujnLXue6GMeYMMRUbxHS2M5D13d36ZvfE0Hu2f/k1xhhjjDG9wQ+/xhhjjDGmN4z19eULL7ywSpuUOymxSq1ETFmXkQx27drV1DvseKmViympb9q0qZYpt1L+jHUxMxclV0q/JNo5CMdIKTezKsRIBJT0aSOg1YJ18Q3/CNthXyjDZ9EKYt84Ls43P+f5sS72c1T2tWHEujh+ljlfjBgQz6dcznlhH48dO1bL3Mej+st2KL2zvWhV4BpxjbkPuQ5cu5hpj2PhGtGOQWsE+xhtKtdee+3Q7xitge3HqB133HFHLe/bt6+WGbmBY8muFam1WnBd2JelS5fWMqNgRJvG7t27a3lgVxqnRcwYY84GM8XS8ErxL7/GGGOMMaY3+OHXGGOMMcb0hrHaHp5++ulqV8iiGkjtW9uUcu+9995a5vmUkUe9Zc5oC4xEQFl0z549tcw35qU8YgIjGfAtdcrY8TvK1QcOHKhljoXHR0mesjhleLZJ6X1U8g3+TRsA55hlyuOxHUrf/JxzSWtJlOSZBIFtcu2YgIHSNaNDSO280CrA8XIscY4HkUkkafv27bVM+Zt1cU9FeZ97LItYwHMYBULKrwnOMfvFJBvcRxGOmZEQaBug3ef+++9vzuf8cy6YvIOJU6JtglFT+B0TofB8zkuM0BCTeQygfYnt0e7Ea1Bq53uwXrY9GGPMzMC//BpjjDHGmN7gh19jjDHGGNMbxmp7OH36dJUmKTfzbXKplYXnzJlTy5Qis8gL8S37LOIB26CFIEukMep8ts8ID1EmjYkPBjBKAKVvzkuU0dkXziVtBAsXLqxlyvZR3uec0TrAt+yZCCTOMe0glJ7ZR8rgnC+uaRwL6x0VrWIA7RCSdPDgwaHH0dpBm0i0B9CqkMF1YTn2l+PkutIqwX5Fawq/W7BgQS1nkR+ypCBSHpWCET2yuY9RR7he7AvbZzSUuN6bN28eWhctDbRwsC4mtxkF54VlRpSI0TVolRicY9uDMcbMDPzLrzHGGGOM6Q1++DXGGGOMMb1hrLaH5557rtoCKL1TnpdaWZZyOSV5Ri8YZVXg35Rs+Tn7QomXCQCkVgpmHynLUvoeFVWBb9bTDsF+UYaOb8lnCSAoqdMGwDGyXqm1KnCOOF7K0DFSB+eJ1gFGwWC9lJ7jHGXWAR7H/o9KNkLLDGVt9p91xaQkrJtzkSXM4Brx3DgWWh24j7lGUWJnP1mmbYPnMFIF96c0tb3P/U2bRbQosc0sOgiJ9h22n0UdYV8437QtxO/YZ65dlpQk2jH43WC+bHswxpiZgX/5NcYYY4wxvcEPv8YYY4wxpjf44dcYY4wxxvSGsXp+u66rnj/68GL4JPoKGW6L2ZnoU6VvMIbxYkYo+iOzNtiv6LOlD5ChkOht5efR95iFFKPnl75CZrZixisp9y+zj6xrlJ+Ufl6GkGP/ly1bVssxLBTnMst+xjlmvdGjTe8l/aFZeDS2NypsGtee/WJfsgxhUrteWdg1+tNjFjz6dOnB5bwwC168JjI/LfvMOWKYttiXzLvK/ckxHjp0aGjbEXpuY9i57Hyuf7Z32H/6uHk/kFr/NMscC8/nuwbRLz7Mex698sYYY6Yn/uXXGGOMMcb0Bj/8GmOMMcaY3jBW20MppUqrlGij9E0pmBI7LQnLly+vZcrNLEutdYDf0QZASfzSSy+t5Sir8jhmzWIb/JwSa2T//v21zFBhWX9jiChKzJSFaXugdJ5J/VIbVopQ6o8hrgjrZp8pfXO9Dhw4MLSP8RzuA1pL7rjjjqH92rJlS1MX1y+zgHC9OXdSLqPzfK4LLTbRapBlIeTe597h2KXW0sH14z7gXE5Voudx3AdZOLdomeC4OH9ZNsQ4rixkIfvFcXGOR9lcuC+ycHbcq6NC7g2uSe5BY4wx0xf/8muMMcYYY3qDH36NMcYYY0xvGKvtgVCKjJnUli5dWsuMPrB27dpapvw56i17RlKgFMs2+MY8pd+bb765qYuyNNukjL1kyZJajhIx38Bnli9me2P7nBdaEKQ2gxXnKMsyxr4wC5vUSv+UdjnHlLQ5DqmVtbNoDbt3765l2gY4d7H/Wb07duyo5eytfqm1qWRv/HOMkcWLF9cy99FUMr/FSCHsC60d7Bel+riPKcuz7swOwblnBBIpv/Z4PttnmXstwr7w+sgsFFIb4YJ1c4+wX+xvtEtldpAsGkmW0S62MyBm7TPGGDM98S+/xhhjjDGmN/jh1xhjjDHG9Iax2x4GsiPf/o9vUdNeQEmeUiiD0vNt8JjkgnIxrQ60HVAKpTwdIxHQHkHJlXIt640yKeVuSrSZ1M++MzFC7HMmd1PeZjkmqaDEm0W0iPYEQumaUjLLlKGzZArxb/aF45pqghPK8CwzqQfnPkr6tETQXkGrA8ujImrceOONtbxr165a5tzzfLYX26EFhudndooYFYF10xLBSAy0eWT2k3h+tt+yY6TWHsEx8rpnf0clSMmuCe4LjpF7alRSksEcj0rwYYwxZvrgX36NMcYYY0xv8MOvMcYYY4zpDWO3PQzkREYMiAkcaE+gLH306NEX1CO1cnGU5yn/0l6RSf2UXqP0TCmWERMov9IqwbpiO3wDne1wXDw+yrJTsUpQBqf0GyMRUArm/DHyRPb2fWyH9gaeQ+mbcx/f/udc0DaSRUugVSH2i/uK59P2wagZMdnHqlWrannTpk1Dx5JFAKA1I9bN8Wd2khgphGvG9brzzjtreevWrbXM/cV9K7VzQQsIE6/EvTuA8yXlURVoj+DejVEU2Df2i2XC+WaUlPgdI3UcPHhwaB9pRYnRVDiX8XoxxhgzvfEvv8YYY4wxpjf44dcYY4wxxvQGP/waY4wxxpjeMFbP74UXXqiVK1e+4HOGfpKkp556qpbpfWSZ3tYsRJLU+j7plaRPlz5Tehijh5TtM/QYfZv0ucbsYVdffXUtc4z8nL7LzLMr5ZmuOC8sc15iiCj+TX8j26SfMoYnm4o3l6GzmE0r7gf2+dixY0P7snr16lpmqKy49vTZcu3YX3q0o/84W1eew7o4jzGk1759+4bWxTZ4/qi9x3B669evr2XOC/sY/cdx/Qdwv2b+2+iDj3UPoDeYe3XhwoXNcfyOey8LW8b9ETMVcl75fgD7T782vcj0AkutB3jRokWSpJ07d8oYY8z0x7/8GmOMMcaY3uCHX2OMMcYY0xtKDKl0Nrngggu6mCFKeqHEO2fOnKHfscxQRjw+hk2jTHvgwIFapnRP6TlKqYSSLW0LbJN9HJUVjf2i9EsZOQv3JLWyeraGlIizrGBSKwUzdBjPjzI+Yd2U22lvoE2E9gS2LbU2BtbFeWW9/HzJkiVNXQxtF+X6AZzjaC2hTYayOOVy7j3OEcchtdacLGwc92Scb1oduN9oYWDYsizTntTuF/aLdXHvsi/RMpGNKyPu6alkTWN/M5uI1K5X1kfef3h+vIa49wf1Hj16VE8//fSLD3IGcXmZ191R7jrX3TDGmJfM+u5ufbN7Yug927/8GmOMMcaY3uCHX2OMMcYY0xvGGu2h67oqNVJ+jG/ZU+I9cuRILVMG5+eU6hk5QWqzO1GuZl2UPClXRxmX0vmyZctqmXJzJttLbfQCfkcZm32ZirVBajPi8S13vr1PqwEtG/E4fke5mFENoiTPvl1zzTW1zHWhXM31ZrY1qZ0X2hh4HKMCxMx3JMvcx/mmJSBG56D0vWLFilqmBSKzLcR+sS6SrXG0zNAOwzYZVYHXFOcx7sMsagrh3LGPMfMbr5EsiyA/j1nZOOe0IbB92ldobYgZ+Thn7CfXhfOYZbGTWlvVIAJLtFQZY4yZnviXX2OMMcYY0xv88GuMMcYYY3rDWG0Ps2fP1rx58yS1tgFGApDyJAKMKkBJnbJmjGRAuZ+yJWVV2hYoCVO6jX259dZba3n37t21fPz48VqOtolsXDfddFMt7927d2j7TBgh5dELGK2CkQjYr2hb4N+MpJBFe4hRFWgD2LNnTy1zvjN5P779zzmjPYJRIbJkEI8++mhTF+eYZbbB/RGtCuw/5fIs8QrXN0ZY4Dgpz2eRF6Ikz37yO85Rlogk7kOewz1GyxCtNExCE+tim7yms6ghcVy0CTGaCo+jdYmf0/IhtfcH2n8437w/0CbCfS+1lorBXI6KeGKMMWb64F9+jTHGGGNMb/DDrzHGGGOM6Q1jtT1cdNFFWrNmjaRWyhxYIQZQZqU0Sel5qm/JL1++vJYp/VOiphRK6XblypVNXbQ3/PVf/3UtU0amHSG+2U65nrLqjh07aplvs2dv4kutDSKzakxF9pdaeZ/fcb4pu8d+cc55fpZEYFQUC8rPnMssUQHrGhXVgLYDrhetIZwjqZ1Lts/oA1kUjijvZ5EnaBnh57Fe9oXRUGgZ4t4ZNcdZxANaEGh74PpyH0htBBNG+uBcch1j8g8exwguXBeON1sT6fmoDFJroWC9tLZwf0abCvfxICJFtFQZY4yZnviXX2OMMcYY0xv88GuMMcYYY3rDWG0PTz/9tHbu3CmpfUs/QpmV0iSlWEIplm+ZS60sm70BzzJlWZ4rtZI8+085lP0dFX0gSzTAPo56u5xjZr9oh6A8T9k92h7YT/aLx3FNYrB/ytrsczYunh9tKoy8kdkeeA5tA3G+2H+uF9eIaxel72yNGBUhk/dHRfpgXziWUWvPvyndZxYKtvHmN7+5qYuWI47lvvvuq2WOnfsowogLGzZsqOVobxgQo3twXmlvyKI1cI5iNBaOmWsxsFpFOC7aJCTp0KFDtTywPWQRVowxxkwv/MuvMcYYY4zpDVN6+C2lXFFK+WQpZXspZVsp5c5SyrxSyudLKbsm/3/ui9dkjDHmbON7tjHG5Ez1l9/fkfR3XdddJ+m1krZJ+nVJd3ddt1bS3ZN/G2OMOff4nm2MMQklC9VUDyhljqRNklZ3OLiUskPSm7uuO1ZKWSTpi13XrRtV10UXXdQtXrxYUhvqi6GrpDYUU5aRif6+UZmx6BGMGbwG0I9JD2L0bfJ8emvp+2Rd0SPIugc+wlimb5Njj15JjpMeWobxYh/pc33qqaeUwUxqWdi0mG2O/uPDhw8P7TPr5bjiHHOOOC6G8eI5o0Kd0UP65JNPDj2HeyqG3OP6DfZt7Bcz8mUh72I72X7jnqYXV2r9qfRFZ95zrhfDhsU2sxB09PLSEz7Kq885yjy/8Xxe62wz82XzGuS5UrsuvD7Yr1WrVtUyQ6MxDKLU+ugH5UcffVTPPPNMu2HPQ87kPfvyMq+7o9x1VvtrjDFng/Xd3fpm98TQe/ZUfvldJekxSf+5lLKxlPL7pZRLJS3sum7wRthxSQuHnVxK+UAp5YFSygOj/uE0xhhzRjhj9+xn9fSwQ4wxZlozlYff2ZJeL+l3u657naTvKMhlk78uDP0Jueu6j3Zdd3vXdbfHXy+NMcaccc7YPftVumjYIcYYM62ZSqizw5IOd123fvLvT2riRnqilLIIEtqjaQ2TnDp1qoYio9waQ0wxtNHq1atrmXIxQ5pRho+/LlMmpVxOiZvyJyXhGIYryyxGqwItDDE0GyXyLAMX7SCZJB3/ptyeSeqU4UdljqOUnEnyo7J08RzOfRYmKmbBY12U29kvzjH3zhve8IamrgceeGBom5TEaSGItgdmT6OdIwtbxjWJdo4sBFxmO4rhxfgfjmyf7Vx55ZW1TJtHnHvaIDKbD9vnHMcMb+wLQwOyv1nfI+wzrz1e37zW4hzTqsD5ps1n27ZttUxbz6ishaOul/OUM3bPNsaYmciL3tW7rjsu6VApZeANu0vSVkl/Jel9k5+9T9Knz0oPjTHGTBnfs40xZjRTTXLxC5L+tJRyoaS9kn5WEw/OnyilvF/SAUnvOTtdNMYY8xLxPdsYYxKm9PDbdd0mSbcP+eolvQZcSqmSM9/gjlBmpUxK+ZR2CMq9mzZtauqi9YCyLCXuffv21TKjCkSPMmVVyq+MZECJOb6NntkzKNGyTUqvlPol6bHHHhv6HcfIt98pr0c7B49jNjBK/bt3767lGFGD8jG/4zpyvJSRKXVL7Vqyriw7Hut98MEHm7r4HS0rXC9K58wqFttk9IEotw/gPHAcUrsXWO8oq0TWF46F53Ad5859PoQrrTRSnqmQVgfOMa0Oo/rIunjd8vxom+CcZWXuV7bPcUjtvGRWHO43jjHu6WFr9GKRcc4nztQ92xhjZiLTzsxmjDHGGGPMy8UPv8YYY4wxpjdM1fN7Rjh9+nS1C1C+jEkT+HY2pVxCKZeyaHybPEuMwcgLJEuEEaEsS9h+lJtpA6CEOhVJnzYLqbU6xMgAA5YtW1bLlPTjnLLNe++9t5azuYifZ0k6KOlzjSjJx75w7TN5nu1xH8WoBpSys4QbXCOug9RG0cj2GPcB5zHuQ/aFdWW2hxjdg3NOGZ92DNbLfc9xSK3lh5E7sr7QssG1k9poHUeOHKllXtOjkn9wzrPENZzXUdFYMmsKr4OtW7fWMqNexGuV9qfB9eU45cYYMzPwL7/GGGOMMaY3+OHXGGOMMcb0hrHaHmbPnl3fgqcF4Id+6Iea4+65555a5lvzW7ZsqeUTJ07UMuXOKMsyqsQHP/jBWv6TP/mTWqZEumrVqlo+fvx4UxclfbbDZAyUjmPyDsrtfOucEjnPp/QbrQZZcoL58+fX8ooVK4a2wUgRsU3K+JzXLMFHrJvScGbtoE1jwYIFTV2ZjYEWCPaLfRllx6CMndkWMiuL1I6LZfaX1ooYzYRrz34xkcaoZArcC+xnZufgXDCRh9SuC/d4ltSEcxzHxXmlvSFGdRgQbS6M2rJjx46h7WeRSuJ60WrBMiOVZFFHIpzLQQSV6RTtwRhjTI5/+TXGGGOMMb3BD7/GGGOMMaY3+OHXGGOMMcb0hrF6fruuq54/+gO3bdvWHMfvGJKM/kr67xgeix5EqfXp/vEf/3Et00NK3yn9oDFLF9t54oknapk+QnoFo+eXIaPoieTn9GrSHxlDRPF8ehf5Of2ZnC/6d6XW08nxZ/VG3ybDZTGMFteCvmjOS/Rdsh16OhmKinPM42M4ONbNc+iTZb0xyxf/5ljo0+Vc0lsb/eJZdkLuIx4TfdX8m/PPa4V9ya4bqZ0XrgXb5z7knoi+V46Z+41rx75H/zHHT58u1yXLaBdDj3H92U+ewzLHFeeI1+RgXuz5NcaYmYF/+TXGGGOMMb3BD7/GGGOMMaY3jNX2cOrUqSoZM2QRpc/BccPKmSxM+TNm6SKUMqP8OoAh1KKFgu1kUjBl6JhZi7Iwx8/jKN1Sno6Sa2ZpoL2Ac8GxcB5jXZxvyv48JloVMqsG+8K5Y/sM+Sa1EjdtJ4cPH67lzA4xKgMX68pCcnF/SHnGtiyEG8tx7TkvDOnFucsyAEpTCz22f//+WuZ6x3Fx/lnmetO+Mgj1Jb0wGyOvXY6R/WcWOO6J2E4MozaAc7R06dKhbUutdSHbF9xfXPtoUeK+HMxLlm3SGGPM9MK//BpjjDHGmN7gh19jjDHGGNMbxmp7KKVUmTXKjIRSLKMy0CrADGuUmw8dOtTURfmV9TJyA+VSthFlWErfhFIq5eaYcYz1UZbN+pu9lR+/Y/vMfpa1EevKsqRx7ikDc46kdv4y2wnPodQdJXn2jRaQLFrC1q1bazlGa+BcMKsd54Vjf/3rX9+c/8ADD9Qyo0VkkScY9SJaMNiXm2++uZYZ6YRWiRjtgfuKc8z2ub+yLHRSbk9g+7w+sugYUrtfskgKXNM4L5n9h2PMonbEDG+0VHBf0arBvnCM0aayZMmSWh6sXby3GGOMmZ74l19jjDHGGNMb/PBrjDHGGGN6w1htD7Nmzap2BUZbiMkJKGfyTXFCWZhyKxNWSK38SRmf8i2lUJajjE6rBaVjyt1MgBCh9E0pl9Ivx3XZZZfVcpSLeRz7yePYBmVojkOSTp48ObSPlMTZRrQ2UBb/qZ/6qVreuHFjLVMypu2Bcxfbf9Ob3lTLn/vc52qZc58l+5Datece436hvH306NHmfO4XRotgvXPmzKllznGMqMG/uUZZ8pBoVcjWlXYUnj9//nxNBa7xLbfcUsvcxzt37qxlro/UWm547SxcuHBofzmnUmtV4J6mVYFzxKgfce9ktiTeTzivrHfNmjXNORzXgQMHJDnJhTHGzBT8y68xxhhjjOkNfvg1xhhjjDG9Yay2h9mzZ1c5NpNupTzRAi0QlOqzt/djXZmlgDL2DTfcUMuU7aVW1qX0T/mUb41HGZZvqhO2TxmecxTlZtb9ne98Z2i9tICw7divLFpElighzjHl4Pvuu29om1migxh5gvI+Ze1ly5bVMi0MtHDEpAdZIhRaAnhMjGRAiZ6RAQj7v2jRolqOkQi4x3ft2lXLtFPwnFHzwvWLUQoGcC64P6XWZsM2Gd2C88JylP6593it8nNeK3FcWTQWRpFgHzn2aAXifmO9PJ/1Xn311UP7GPs/iBTiJBfGGDMz8C+/xhhjjDGmN/jh1xhjjDHG9Iax2h6ee+65Ki1Tro1yIr9j4H7K+JShKaVGWZZSKG0EmaWAdcVIEydOnIhDekF/MxlYaiVXSveUlY8cOTK03tgXSvSUuBcvXlzLHCPf3qfNQWojb3BeuS6UmBnhQGrHNXgzXmrfrGebjAQQ7Rxsn7aTLNEB5y7aSmhtofTPuacdYZSkz36yTVozsggisf+ZvM/P47xwL2VRB2gFGpV4hG3SAkJ7QGZ/iZaXzFrC6zaLzBLPIWyf5RhFg7BujoXWGNp3uHbROsR5GezduCbGGGOmJ/7l1xhjjDHG9AY//BpjjDHGmN7gh19jjDHGGNMbxur5PXXqVPXNvva1r62f7927tzmOIaYGYYak1tPJkGBZSC6pzSy2atWqWqbHj95QelajVzL6FQdkmaWuvfba5m/6IOk3XLlyZS0fP368lumHjGG86Aemz5d9YZm+yzhHHD/nMmbQGhDDeHEus9BVrJfjitnm2E+eQ18y55FtM4RXbId+XB5H/2vMPsb+08PLfcD1Yn9jpkHCvmTZx6JnmPNE/zbPj/s1q4t/04vMtaPfm8dwr8Q+c7+xLp4fffCcC/af7XMdmIlt69atTV1xnMPgnuDaRR8198igj5k/2RhjzPTCv/waY4wxxpje4IdfY4wxxhjTG8Zqe7jooouqxE9ZkXKxlEuulFUpfTPEUwxHRFmWYcT4OS0E/Pyaa65p6jp8+HAtZyHBKJ/GzHW0GzD0F6VfSquclyjpci6ycFm0EFDSZxin+DdtIpTUmYks2j+YsY1zyf5zjlhvnCOuN+eL65qFyYuSOsObsc8MG0fpO9o8aLPhfmWbXEfOEfdn7AutGlxHtnH55Zc353NfcZw8LssCF+FxLM+bN6+Ws6yJ0fLCvUNrxs6dO2uZcx+vde5Lrj33e7ZfYmi7LCQar1WOl+OK4fs4f4O6Ypg3Y4wx0xP/8muMMcYYY3qDH36NMcYYY0xvGKvt4ZlnnqnWAUq8lLGlVsanXEwpctGiRbV88ODBWo7Z4ihLZ/YElildR6tBrHtYf3kO7RhSO5ZMlmVfKPfGtinfUu7lvLI9SsSMLhHb5HGU1EfNcRbVgVIyo2hwjkZlQuPb/+wjx5uto9RG9/ja175Wy7QNXH311bWcRbeQ2qxwtC1wLtjfaA+IEv2wcyirL1u2rDmO9gxeB4y+wHHRahBtKhwzsxZyH3Id2ceYHZBzwcx9t9xySy3TShNtD1mkEM4rry9Gf4nZ3mhn4fXBteN8cb2jxYnnDK5jzo8xxpjpi3/5NcYYY4wxvcEPv8YYY4wxpjeM1fZwwQUXVCme0nWUZSlZUgZntAVKxJR7oyzLdrK3tSmxUlaNySCy5AyU7tk+JV2plWn5ljxtHzyG/WJEhQjlYvYli5oR7QGUpVevXl3LjFhA20JMTEGrCOV+JubI1iGuCRN2HDt2bOj5tHZwXFGW3rJlSy1z7bJIHdGaQOmb/eQc0xrCNYryPueS+50WBI6RNpPY/2y/sw1+zjmScssR7Q1ZFIpoeeFxjPDAOWJCmmgt4ZjZPsdL+wzHHqN7sJ/sFy0k7D+vzxh1hHshzp8xxpjpje/qxhhjjDGmN/jh1xhjjDHG9Iax2h5Onz5dZfFRb05HGwTPH1bOAtdLrWTJemkvoPRKuTgmgzh69GgtU9amPE5LQBwHJV+2/7rXva6WabugbWAU7DPbz+b4iiuuSPu1bdu2Wl6xYkUtUy6mVC+18jHlbsrYlL459ijvc/04FkYloNxNm0VMVMB14VowegCJ0ncWHeSmm26q5f379w+td1RSEo6Rc8/2KOFLeSQH9jlLkMJ1kPKEIVn/s6QiUnt9ZcliaAeJ1yfb5JizJBscY4yoka0X9wGTavC6GWVxGpzjaA/GGDMz8C+/xhhjjDGmN/jh1xhjjDHG9Iax2h66rqvyb/aWudRKoZQ2+WY9z6GMT3k8HpdZICi38nhGOJBauZgyNiVWythR+mYQfkqotBowMUYmz0t5tAh+nkU4YBSE2E8mU+A569atGzoOSXrwwQdrmXI1+89zGMUhyvv8m3PJtaOMTak7zjfHT4mba09rRbSZUAqnXL5v375a5j5ghICYgIHWkHnz5tUyJX0miWDkBqmdM44rS/LBOYrRPThPtFNwjtkGxxitQDyHe5cWFO5P9ktq7Rlc+yxxSmaHkNp15XxzH3IfMxHHI4880tTF+8ANN9wgybYHY4yZKfiXX2OMMcYY0xv88GuMMcYYY3qDH36NMcYYY0xvGKvnt5RSPYb06dJnKrUexal4QFlmtrFYF32T/JwZoLIscvEcegLpNaRPlF5HqQ1LNch0J7XhzTIPK4+XWu8jvYiPP/740Lrom4z9Yl0MY8b+0m/NstR6XTn+zD+dZfCT2mxznG/OMbPAcY04dqmdPx6XZUWLnmG2yfnL6mJ/Y4Y37mOuF+eIfti43qwv8+zSl8y+RC8z67rssstqOWZfG0Cve8xamGUX5LpyjugLltq55H7jfF133XW1zLnbvHlzUxf3JT3D7POGDRtqmX7vuHduvvnmWh5cL870ZowxMwPfzY0xxhhjTG/ww68xxhhjjOkNY7U9zJo1q9odKGvGcFeUzletWlXLzAbG8EUx6xTJ2snCe2XnxnMo91K6pVTOrGZSa7ugVYA2gCwrGcNNSa11gbI05yLLZsV+SNKCBQuGnsM5ZvtHjhxpzqd0n0n/nCPWy7BfUjt+2jEo71Oe59hjRj2OmfPCcgwpRrKMZ7TscB0onceQXlxjtk/bAecxWhC4Lgw3RqsAz6fUH20u7D/XiGtM+wutBjHcF+eca5T1JYZd4xpzXljvQw89VMtZBkGpvW9wXJx7zis/5zUgtRaHwdzb9mCMMTMD382NMcYYY0xv8MOvMcYYY4zpDecs2gPlZsqPUivxHj16tJb5BjwlWtoW4lvyfAOcUi6Po5RKS0CU0SndU6KllEupnZELpFampRQdIwMM4BvrUUbn+GM2sQGcR2YPi1ENKL0vWrSolrMIBTEKBuVgjiWzo9AOEo/hXHJcmbUjs5/EfnG9r7nmmlreuHHj0OOl1vbA72jV2L9//9A+MjNhHAv3JDOhcV0o4Ut5JrZ3vetdtfzJT36ylnlNxX1MqwTbYfu0B3Du3vnOdzZ1fepTn6plrj3HyPbj3uMe47qy/7RG8HjOndTaTngdci24X7jvly5d2tS1e/fuWh5cR7HvxhhjpidT+uW3lPIrpZQtpZRHSil/Vkq5uJSyqpSyvpSyu5Ty56WUC1+8JmOMMWcb37ONMSbnRR9+SylLJP2ipNu7rrtJ0ixJ75X0IUm/3XXdGklPSnr/2eyoMcaYF8f3bGOMGc1UbQ+zJb26lPKspEskHZP0Vkn/bPL7P5L0HyT97qhKTp8+XaVhSpmUi6X2zWtKlpQyKQOzrhgVIUtGwcQS2Zv4UQbn+VkChizQv9RKtidOnKhlvg3Pcyi7R2sD26QszfOzZAwxScWWLVtqmUlCbrrpplqmPWD79u3N+ewb7Ryc45UrV9bynXfeWcv33XdfUxfnj0kI2H9G/eA6xPnOIkxEG8CwY6TWEsC6uF9ZZh9ZjudzHzByxmtf+9pajlEV2A7neOfOnUNG0loQli9f3nzHurP5o+2B1+C2bduaumgdoJUps69Eiw/b5LWaWVuyspRHR+Hcc3+zDe5Vqb2/DGxVMbrEec4ZuWeb85PPHd30ose8bfGtZ70fxkxXXvSX367rjkj6LUkHNXED/YakByU91XXd4F+Pw5KWDK/BGGPMuPA92xhjRjMV28NcST8maZWkxZIulfT2qTZQSvlAKeWBUsoD8YUkY4wxZ5Yzec9+Vk+/+AnGGDPNmIrt4Qcl7eu67jFJKqX8paTvk3RFKWX25C8JSyUdGXZy13UflfRRSZo1a1Y3kFBpNYhB+Ck5Uj7mw3OUPF/sc+mFb9APyOwM0WrA5AKUrjMZPj7sZ33LoiJQho8JMyiDs52TJ0/WcpbYYu3atU1d73//89a/TO5+3/veV8u/8Au/0JzP+eNa3nzzzbVMqf1LX/pSLTPygiTt27dPw+Ce4DwyikVcX8rwnAtK929961trOdomvvKVr9RyJnlzjUZZK2jNefTRR4cex4ge0R5AawzXm1EJ+Dll/3h9cb9xXJwj7i9GADlw4EBTF9eV88frOyaLIdxvtFCw/+zjqOQ4jNjAtWfEGNpq2N6oaCqDtRh1bznPOGP37MvLvGkzaNNCa4QtEMa0TCXaw0FJbyylXFIm/qW+S9JWSfdI+onJY94n6dNnp4vGGGNeAr5nG2PMCKbi+V0v6ZOSvi5p8+Q5H5X0a5L+dSllt6T5kj52FvtpjDFmCviebYwxo5lStIeu6/69pH8fPt4r6Q1nvEfGGGNeEb5nG2NMzlgzvHVdV72E9HDGbFj0CzKbFj24Uw1Plvn0Mp8v/ZHRZzsq+9sAjit6ftn/7HxCDyt9i5K0evXqWmbGLvojOXfMbsewZZK0Z8+eWn7Pe95Ty8ycF7PwEc5Z1md6O+nVZOgpqQ1vduutt9byQw89VMv0pjLMG32eEYbPozeUYzx8+HBzTha2jvuN3lzOUcwGloXG4x6h7/Sqq65qzqdP+Iknnqhl+owzb228vjgWzlnmX+YcRR889zt9wpnnN16fPJ97l/PH+eIx8dpmuDKez/a5d7IsdlJ7vQz82n5h1xhjZgZTyvBmjDHGGGPMTMAPv8YYY4wxpjeM3fYwkENHWQAohb7mNa+pZdoO+PmoEESXXHJJLVMKpYQZQ0ENiJnnKB+zz5S+GT4pys2E8i0l5jlz5gw9JtbF4x5//PFa5rxQRmdIqiiPsy6uC8NrffGLX6xlhheL7dNqQUmdMjjrjRnKuJa0Z9DawfninoiyNCV2Hsf+c41jyLko8Q/Iwp4xc16cY0rsJMtUGLPNcV3ZT64x9zrLnG+p3Uvc+zyHY8ysGdILbQwDMgvIqEyF7AvHyLlgv+I6ZGPJsvDRohMtM7RBDPoyjUKdGWOMGYF/+TXGGGOMMb3BD7/GGGOMMaY3jNX2QChLRkmf2ZWYsYyf83zKnbEuZmWj/EnJk9IzJVpK5VIbYYHyL/vIN/FpJ5Ba+TjL2EV5nxYCjleSNm/ePHQsUS4fwLHEKBYc15YtW2qZ9gDKwpTgpXbOOC+Z3J1FG5CkxYsX1zJlbcroWda/KIPTesD5ZlY01hvtN5xzyvCUvxmd4oorrqjluA+nElWB+zjaC3g+9yv3CMfIued1E//mmDkXmbUjSv88n2vM+eJcjMqaOOo6HsC1Z0QGqV1Xziuve0au4NzF8fIaGWS1m0qEFmOMMec//uXXGGOMMcb0Bj/8GmOMMcaY3jB228NAWqVEzEgEUivrkuyt81HSNf9et25dLTNCwYkTJ2o5k7rjcewLpWDKsvF8nsM+Z/2nhYLzFetmvZS7+cY6z49v6F9zzTVD+79169ahfYkJAZiAIYNtZkkepNZCQhl8+fLltcwxMpFGlNRpSciiUNAqEPcdpXPuF8K14zGM/BD7Rkk/S9oQk4pkkUq4X9kmj4nXF9eY++K2226r5c985jO1PEruZz85xlWrVtUy7QmxL4cOHaplrn2W5IL7KCYlufnmm2uZ1yrnlVE3aBfinMR+DhKxTGWfG2OMOf/xL7/GGGOMMaY3+OHXGGOMMcb0hrHaHi666CItW7ZMUiudX3XVVc1xfBud0jelUL4lTll0zZo1TV2URillsl4mqeCb4WxPaqXkTMbPjo91MxkGoxJQlqUkPiq5QGYj4DmcU0Z3kFrp+qabbqplyv5MckG5WJJuueWWof3fs2fP0P6yX6NsKoTRKjhflLTjfHPMtInwjX+OPdo5uC8pyfM4rhGjYIyyGtAewDm+9dZba/nLX/5yej6h1YL2glGJJThPjz32WC1/6UtfqmVeH5zHUXVxjbm/d+zYUcsxukcWVYLXJNvPbBaSdPDgwVpm1BCSJRWJ18SwPZL11Zjznc8d3dT8/bbFt56TfhhzvuBffo0xxhhjTG/ww68xxhhjjOkNY7U9vOpVr6pyJO0IUaKkZEyZlBI35VpKr1FWjW/ND6C8Thk+i+IgtW/WU4amxMw3+WNUBf5N6wAjDvAYyugxUUH2BnxmE2EkgGih2LBhQy1TSuZxtJlESZ9v1jP6AqVv9p/zGqNYfO/3fm8tb9y4sZazpCAkrhdlakriWZKJuHd27txZy1kUDe6JbE9KbfQDRg1gn2ntiOfTBkB7w65du4aez/WOdg4mX2H7lPoHEQ6kdl5o3Yl1c+/t27evlmlTidcE2+ce43XLvcNjskQYUnt9sX3ub84xE2RIbXSQwXWYRfwwxhgzvfAvv8YYY4wxpjf44dcYY4wxxvQGP/waY4wxxpjeUKJP8mwya9asjh5JfN78TV8gQ1TRb0jPK8MX0cMZz8k8tDyGvsUI/an0R7K/0RNJ6D3kPNCbyrrob6SHUmrnjH3mfNG3SW8ms2pJ0pve9Kah7WeZxJhtTZIeeeSRWl6/fn0t03dKPyzP379/f1PXokWLaplzzL5w7a+//vq0riwEHDPabd++vZZjRj7COWa99KaO8jLTM0z/MutlXdHPyj3OffHkk0/WMn2q9LxGvzivg8wvz8+5jrFfK1asGPodfc3cO9GDT185j+M+nOq1So85fc1cV2Yq5DpyTuL5g3vKjh079A//8A95ursZyOVlXndHuetcd8OMIIYxe6k47JmZqazv7tY3uyeG3rP9y68xxhhjjOkNfvg1xhhjjDG9YayhzqTnpUZKqTHEVBamKJOlM+lZauXiLDQSz2EIM0rKsX22SbmVtoto52CIsMzqQCmXoZUGmfGGnU9LAcfCYygXR+mZcjPDizE815EjR2o5Svrvfve7a5kSNedvwYIFtcxwaHG9KJdzXrI1Zpau22+/vamL7dDqwRBmo/ZOBiVxjnfu3LlDP5fy7H7ZfomSPr+jhYJ95v6kbYH2lwgz8mXzzXK0FXG92OfMwjAqS1oW8o/tZxkApfZewXtKdt2xj5yH+PfAipOF2DPGGDO98C+/xhhjjDGmN/jh1xhjjDHG9Iax2h5Onz7dWBoGRCmTkQkoa1P6ZZYrZouLsiqzxx04cKCW+Zb6nj17apkSKaV6qZV4KcVSomUmrmjnILQO8A30zPIR540yOi0klI5p8+A8Rh544IFapu2D52Rv30vSvffeW8uUxSnDcy5Yb4wccfLkyVqmZYXjZZlZzTh2qd1HLGdSP9uT2jXifuP5tHYwYkC0UHD9OMccC8+J9gJ+x+sgZkwbQKsDrw+ptSfQgsHIC7wm2d9oD6AFg+cz8gT7Evcx9xL3COul3YDXBK+1COvlOrL/3Hu09UjtvA7sEdlcGzOdYbQIR34wfcF3c2OMMcYY0xv88GuMMcYYY3rD2KM9DORbvhkfJUcmLqANgdJvFuyecqnUSv+UYimvU+KlbYI2BylPAEHbQia1x75xXDyOx2QJOuLflHIZiYBjp7wdbQucC0rEtIYw8gSjKEh5YhBK9ydOnKhl2gsYqUJqx8yxsEzbRBZdQmrXkmvMfcTPo/2GEj2TjGR7jzL8KNsC2+F8Z5EbpHYtM5sN9wEtILQKxHO4Fixzf3Me4vWVwb7QNhGtJRwn9/RVV11Vy5x7Wj5i9BaOmZEcOH72P+49wjkaXB+jkqAYc66gVeGVJrwwpi/4l19jjDHGGNMb/PBrjDHGGGN6wzmzPVA6nj9/fnPM8ePHa5nyMaMv0GqQvbEu5RELKJlSSqXcG9/YZ4QHyvBZYokYFJ+SL8/h2/C0Q7Ac37JfuHBhLVOizdqnzSRK17SdsE1aHRhRIvaF88J2aAnIElNwH0Q4X9wHWUSJaDXI1puwX/Ftfta3cuXKWub8MXkGI0/ESCHRtjKAlhPOV5wXtsl54bXDteNYRkUtIbQH3HjjjbV8//3313K0hmR2DnLw4MFajhYZzgvnglYLRoHIEtLEujhHtO/Q6sB9TEuU1M7FYO/zejBmJuLID6Yv+JdfY4wxxhjTG/zwa4wxxhhjesPYbQ8DaZS2hShfZm+gM2IApfZR0jmPo72BcjEl1ixCQqyLZFEoYr9oQ2DdtCow8gOPiXPEvlD6p1zL9ih1cx5iXfyO59BOEt/YZ3QORlzYv39/LTOqAeebb+VL7VxmlgTaEdiv+DY+x8/15t7j3McEKYwk8dBDD9Uy5fUsqUlMKsJ1yWR8zj1tFpK0ffv2Wua4MqsDpX5aM6TWQsG55D7I9ku0AmW2B9aVfS61NojsOB7DNYnrldlcNm/ePLRezmO0ArHuQYITJ7kw5zuZVcFRIIxp8d3cGGOMMcb0Bj/8GmOMMcaY3uCHX2OMMcYY0xvG7vkd+Obo9YweVGZqmkpGqVHhqujro9+PYZV4Dv2s9HPGdhj2iP3PfMERehIznys9kKP8rPyOnmF6UxmGK2b84nH0VGYZ7aLXkr5T+nnXrVs3tH16oXm89MIwagO4X+in5RzFtV+0aNHQ/nOO6Xkd5R3Psq/Rc0xvKn29sZ/XXXddLXMfsa6Bz3QAfbvcO5w/huuiZzfWxb3Pulhm/3kdRO85r9XoBx7AsY8KR8f547i4pzPPrtTOH/cL15XXN9uOIRI55sHejfveGGPM9MS//BpjjDHGmN7gh19jjDHGGNMbxm57GEApMsrVlJWzsGc8hnJrtCoQSrSEIcUol9ICILVSLCVQts+QXjGrGMdM6wGlXNY7KnMdpWS2n4WxInG+KXHTwsE2mfErZrqiVYEhtjiuw4cP1zLDocW+MNsdpXuuXSZ9X3XVVU1dixcvruXdu3cPbXP58uW1zExkUr5Hua7M9Me5Y6Y7qV2jbdu2DW2DMKua1Er/3Idxjw7gHNFyEuviHqN9JbOWRGsDbRBZODd+HjO8sU22k1lxeH3EvcPzeQ7nmNYQXh8xNN2wuqLNwhhjzPTEv/waY4wxxpje4IdfY4wxxhjTG0r2hvZZaayUxyR9R9LJFzt2BnOl+jv+Po9d6vf4Z8LYV3Rdt+Bcd2KcTN6zD2hmrN/LxWPvL30e/0wYe3rPHuvDrySVUh7ouu72sTZ6HtHn8fd57FK/x9/nsc8E+rx+Hns/xy71e/wzfey2PRhjjDHGmN7gh19jjDHGGNMbzsXD70fPQZvnE30ef5/HLvV7/H0e+0ygz+vnsfeXPo9/Ro997J5fY4wxxhhjzhW2PRhjjDHGmN7gh19jjDHGGNMbxvrwW0p5eyllRylldynl18fZ9rgppSwrpdxTStlaStlSSvmlyc/nlVI+X0rZNfn/c1+srulKKWVWKWVjKeUzk3+vKqWsn1z/Py+lXPhidUxXSilXlFI+WUrZXkrZVkq5sy9rX0r5lck9/0gp5c9KKRf3ae1nEr5n9+ueLfX3vt3ne7bUv/v22B5+SymzJP2fkn5E0g2S/mkp5YZxtX8OOCXp33Rdd4OkN0r6V5Pj/XVJd3ddt1bS3ZN/z1R+SdI2/P0hSb/ddd0aSU9Kev856dV4+B1Jf9d13XWSXquJeZjxa19KWSLpFyXd3nXdTZJmSXqv+rX2MwLfs3t5z5b6e9/u5T1b6ud9e5y//L5B0u6u6/Z2XfeMpI9L+rExtj9Wuq471nXd1yfL39LEhbREE2P+o8nD/kjSu89JB88ypZSlkv6JpN+f/LtIequkT04eMpPHPkfSD0j6mCR1XfdM13VPqSdrL2m2pFeXUmZLukTSMfVk7WcYvmf36J4t9fe+7Xu2pJ7dt8f58LtE0iH8fXjysxlPKWWlpNdJWi9pYdd1xya/Oi5p4bnq11nmw5J+VdLpyb/nS3qq67pTk3/P5PVfJekxSf95Uj78/VLKperB2nddd0TSb0k6qImb5zckPaj+rP1Mwvfsft2zpf7et3t7z5b6ed/2C29nmVLKZZI+JemXu677Jr/rJuLMzbhYc6WUd0p6tOu6B891X84RsyW9XtLvdl33OknfUZDLZvDaz9XEryWrJC2WdKmkt5/TThnzEujjPVvq/X27t/dsqZ/37XE+/B6RtAx/L538bMZSSnmVJm6if9p13V9OfnyilLJo8vtFkh49V/07i3yfpB8tpezXhFT6Vk34qa6YlFSkmb3+hyUd7rpu/eTfn9TEjbUPa/+DkvZ1XfdY13XPSvpLTeyHvqz9TML37An6cN1K/b5v9/meLfXwvj3Oh9/7Ja2dfHvwQk2Yqf9qjO2PlUmv1Mckbeu67j/hq7+S9L7J8vskfXrcfTvbdF33G13XLe26bqUm1vkLXdf9c0n3SPqJycNm5Nglqeu645IOlVLWTX50l6St6sHaa0I2e2Mp5ZLJa2Aw9l6s/QzD9+wJ+nDd9vq+3fN7ttTD+/ZYM7yVUt6hCU/RLEl/0HXd/zK2xsdMKeX7JX1F0mY975/6d5rwkH1C0nJJByS9p+u6J85JJ8dAKeXNkv5t13XvLKWs1sQvCvMkbZT0U13XPX0Ou3fWKKXcqomXRi6UtFfSz2riPzZn/NqXUv5HSf+tJt6e3yjp5zThFevF2s8kfM/u3z1b6ud9u8/3bKl/922nNzbGGGOMMb3BL7wZY4wxxpje4IdfY4wxxhjTG/zwa4wxxhhjeoMffo0xxhhjTG/ww68xxhhjjOkNfvg1xhhjjDG9wQ+/xhhjjDGmN/z/qYfyWEtlXZcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "torch.Size([4, 1, 96, 96, 1])\n"
     ]
    }
   ],
   "source": [
    "# See example of image, label.\n",
    "check_ds = Dataset(data=train_files, transform=train_transforms)\n",
    "check_loader = DataLoader(check_ds, batch_size=1)\n",
    "check_data = first(check_loader)\n",
    "image, label = (check_data[\"image\"], check_data[\"label\"])\n",
    "label = label.int()\n",
    "print(f\"image shape: {image.shape}, label shape: {label.shape}\")\n",
    "print(f'label values: {np.unique(label)}')\n",
    "print(f'image max: {image.max()}, image min: {image.min()}')\n",
    "plt.figure(\"check\", (12, 6))\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.title(\"image\")\n",
    "plt.imshow(image[0,0,...], cmap=\"gray\")\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.title(\"label\")\n",
    "plt.imshow(label[0,0,...])\n",
    "plt.show()\n",
    "\n",
    "print(image.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create Datasets and DataLoaders.\n",
    "\n",
    "# Dataset and DataLoader for training.\n",
    "train_ds = monai.data.Dataset(data=train_files, transform=train_transforms)\n",
    "train_loader = DataLoader(\n",
    "    train_ds,\n",
    "    batch_size=2,\n",
    "    shuffle=True,\n",
    "    num_workers=4,\n",
    "    collate_fn=list_data_collate,\n",
    "    pin_memory=torch.cuda.is_available(),\n",
    ")\n",
    "\n",
    "# Dataset and DataLoader for validation.\n",
    "val_ds = monai.data.Dataset(data=val_files, transform=val_transforms)\n",
    "val_loader = DataLoader(val_ds, batch_size=1, num_workers=4, collate_fn=list_data_collate)\n",
    "\n",
    "# Dataset and DataLoader for testing.\n",
    "test_ds = monai.data.Dataset(data=test_files, transform=test_transforms)\n",
    "test_loader = DataLoader(test_ds, batch_size=1, num_workers=4, collate_fn=list_data_collate)\n",
    "\n",
    "\n",
    "# Metrics.\n",
    "dice_metric = DiceMetric(include_background=False, reduction=\"mean\", get_not_nans=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "device: cuda:0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "monai.networks.nets.unet UNet.__init__:dimensions: Argument `dimensions` has been deprecated since version 0.6. Please use `spatial_dims` instead.\n"
     ]
    }
   ],
   "source": [
    "# Create UNet, DiceLoss and Adam optimizer.\n",
    "device = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu')\n",
    "print('device:', device)\n",
    "\n",
    "model = monai.networks.nets.UNet(\n",
    "    dimensions=2,\n",
    "    in_channels=1,\n",
    "    out_channels=1,\n",
    "    channels=(16, 32, 64, 128, 256),\n",
    "    strides=(2, 2, 2, 2),\n",
    "    num_res_units=2,\n",
    ").to(device)\n",
    "\n",
    "loss_function = monai.losses.DiceLoss(sigmoid=True)\n",
    "optimizer = torch.optim.Adam(model.parameters(), 1e-3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Typical PyTorch training loop.\n",
    "\n",
    "epochs = 200\n",
    "val_interval = 2\n",
    "best_metric = -1\n",
    "best_metric_epoch = -1\n",
    "epoch_loss_values = list()\n",
    "metric_values = list()\n",
    "\n",
    "for epoch in range(epochs):\n",
    "    print(\"-\" * 10)\n",
    "    print(f\"epoch {epoch + 1}/{epochs}\")\n",
    "    model.train()\n",
    "    epoch_loss = 0\n",
    "    step = 0\n",
    "    for batch_data in train_loader:\n",
    "        step += 1\n",
    "        inputs, labels = batch_data[\"image\"].to(device), batch_data[\"label\"].to(device)\n",
    "        inputs = inputs[...,0]\n",
    "        labels = labels[...,0]\n",
    "        optimizer.zero_grad()\n",
    "        outputs = model(inputs)\n",
    "        loss = loss_function(outputs, labels)\n",
    "        loss.backward()\n",
    "        optimizer.step()\n",
    "        epoch_loss += loss.item()\n",
    "        epoch_len = len(train_ds) // train_loader.batch_size\n",
    "        #print(f\"{step}/{epoch_len}, train_loss: {loss.item():.4f}\")\n",
    "    epoch_loss /= step\n",
    "    epoch_loss_values.append(epoch_loss)\n",
    "    print(f\"epoch {epoch + 1} average loss: {epoch_loss:.4f}\")\n",
    "\n",
    "    if (epoch + 1) % val_interval == 0:\n",
    "        model.eval()\n",
    "        with torch.no_grad():\n",
    "            val_images = None\n",
    "            val_labels = None\n",
    "            val_outputs = None\n",
    "            for val_data in val_loader:\n",
    "                val_images, val_labels = val_data[\"image\"].to(device), val_data[\"label\"].to(device)\n",
    "                val_images = val_images[...,0]\n",
    "                val_labels = val_labels[...,0]\n",
    "                roi_size = (96, 96)\n",
    "                sw_batch_size = 4\n",
    "                val_outputs = sliding_window_inference(val_images, roi_size, sw_batch_size, model)\n",
    "                val_outputs = [post_trans(i) for i in decollate_batch(val_outputs)]\n",
    "                # Compute metric for current iteration\n",
    "                dice_metric(y_pred=val_outputs, y=val_labels)\n",
    "            # Aggregate the final mean dice resul\n",
    "            metric = dice_metric.aggregate().item()\n",
    "            # Reset the status for next validation round\n",
    "            dice_metric.reset()\n",
    "            \n",
    "            metric_values.append(metric)\n",
    "            if metric > best_metric:\n",
    "                best_metric = metric\n",
    "                best_metric_epoch = epoch + 1\n",
    "                torch.save(model.state_dict(), best_model_name)\n",
    "                print(\"saved new best metric model\")\n",
    "            print(\n",
    "                \"current epoch: {} current mean dice: {:.4f} best mean dice: {:.4f} at epoch {}\".format(\n",
    "                    epoch + 1, metric, best_metric, best_metric_epoch\n",
    "                )\n",
    "            )\n",
    "print(f\"train completed, best_metric: {best_metric:.4f} at epoch: {best_metric_epoch}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot loss and mean Dice.\n",
    "plt.figure(\"train\", (12, 6))\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.title(\"Epoch Average Loss\")\n",
    "x = [i + 1 for i in range(len(epoch_loss_values))]\n",
    "y = epoch_loss_values\n",
    "plt.xlabel(\"epoch\")\n",
    "plt.plot(x, y)\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.title(\"Val Mean Dice\")\n",
    "x = [val_interval * (i + 1) for i in range(len(metric_values))]\n",
    "y = metric_values\n",
    "plt.xlabel(\"epoch\")\n",
    "plt.plot(x, y)\n",
    "plt.savefig(daten + best_model_name[:-4] + '_diceloss.jpg')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write epoch loss values and dice values to file.\n",
    "export_prefix = daten + '_unet'\n",
    "export_diceloss = export_prefix + '_diceloss.txt'\n",
    "export_dicemetric = export_prefix + '_dicemetric.txt'\n",
    "\n",
    "with open(export_diceloss, 'a') as myfile1:\n",
    "    myfile1.write('\\n'.join(str(item) for item in epoch_loss_values)+'\\n')\n",
    "with open(export_dicemetric, 'a') as myfile2:\n",
    "    myfile2.write('\\n'.join(str(item) for item in metric_values)+'\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# testing\n",
    "best_model_name = '/home/charissa/minivess/code/20230120_best_metric_model.pth'\n",
    "post_trans = Compose([EnsureType(data_type='tensor'), AsDiscrete(threshold=0.6)]) \n",
    "model.load_state_dict(torch.load(best_model_name))\n",
    "test_dice=0\n",
    "test_hausdorff=0\n",
    "\n",
    "model.eval()\n",
    "with torch.no_grad():\n",
    "    for test_data in test_loader:\n",
    "        test_images, test_labels = test_data[\"image\"].to(device), test_data[\"label\"].to(device)\n",
    "        test_images = test_images[...,0]\n",
    "        test_labels = test_labels[...,0]\n",
    "        \n",
    "        # define sliding window size and batch size for windows inference\n",
    "        roi_size = (96, 96)\n",
    "        sw_batch_size = 4\n",
    "\n",
    "        test_outputs = sliding_window_inference(test_images, roi_size, sw_batch_size, model)\n",
    "        test_outputs = post_trans(test_outputs)\n",
    "        test_labels = decollate_batch(test_labels)\n",
    "        \n",
    "        # compute metric for current iteration\n",
    "        dice_metric(y_pred=test_outputs[0], y=test_labels[0])       \n",
    "        \n",
    "        # visualize outputs\n",
    "        plt.figure(\"check\", (18, 6))\n",
    "        plt.subplot(1, 3, 1)\n",
    "        plt.title(f\"image\")\n",
    "        plt.imshow(test_data[\"image\"][0, 0, :, :], cmap=\"gray\")\n",
    "        #print('val_data-image:' ,np.unique(val_data[\"image\"][0, 0, :, :]))\n",
    "        plt.subplot(1, 3, 2)\n",
    "        plt.title(f\"label\")\n",
    "        plt.imshow(test_data[\"label\"][0, 0, :, :], cmap='gray')\n",
    "        print('test_data-label:' ,np.unique(test_data[\"label\"][0, 0, :, :]))\n",
    "        plt.subplot(1, 3, 3)\n",
    "        plt.title(f\"output\")\n",
    "        plt.imshow(test_outputs[0].detach().cpu()[0, :, :], cmap='gray')\n",
    "        print('test_outputs shape:', test_outputs[0].shape, test_outputs[0].dtype)\n",
    "        #print('test_outputs:', np.unique(test_outputs[0].detach().cpu()))\n",
    "        plt.show()\n",
    "        clear_output(wait=True)\n",
    "        \n",
    "        if False:\n",
    "            for test_output in test_outputs:\n",
    "                saver(test_output)\n",
    "                \n",
    "                \n",
    "    # aggregate the final mean dice result\n",
    "    test_dice += dice_metric.aggregate().item()\n",
    "    \n",
    "        \n",
    "    \n",
    "    # reset the status\n",
    "    dice_metric.reset()\n",
    "    \n",
    "    print(\"evaluation metric:\", test_dice)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:jupyterhub] *",
   "language": "python",
   "name": "conda-env-jupyterhub-py"
  },
  "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}