FigB_fast_mri_noise_experiment_complete-R2_R4.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "4aaf82bd-c34d-446c-9506-ad38b82037bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "#os.environ[\"CUDA_DEVICE_ORDER\"] = \"PCI_BUS_ID\" # so the IDs match nvidia-smi\n",
    "os.environ[\"CUDA_VISIBLE_DEVICES\"] = \"-1\" # \"0, 1\" for multiple"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ab44ab40-57f4-42de-8157-4aed9a834918",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-11-14 15:46:13.954289: I tensorflow/core/platform/cpu_feature_guard.cc:210] This TensorFlow binary is optimized to use available CPU instructions in performance-critical operations.\n",
      "To enable the following instructions: AVX2 AVX512F FMA, in other operations, rebuild TensorFlow with the appropriate compiler flags.\n",
      "2024-11-14 15:46:14.647350: W tensorflow/compiler/tf2tensorrt/utils/py_utils.cc:38] TF-TRT Warning: Could not find TensorRT\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import sigpy as sp\n",
    "import sigpy.mri as mr\n",
    "import sigpy.plot as pl\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.io import loadmat, savemat\n",
    "from scipy import fft, fftpack\n",
    "from random import randrange\n",
    "import re\n",
    "\n",
    "import h5py\n",
    "from scipy.ndimage import zoom\n",
    "\n",
    "#from skimage.metrics import structural_similarity as ssim\n",
    "#from skimage.metrics import normalized_root_mse as nrmse\n",
    "#from skimage.metrics import peak_signal_noise_ratio as psnr\n",
    "\n",
    "# functions for metrics adapted from .py file provided by Efrat\n",
    "from metrics import calc_NRMSE as nrmse\n",
    "from metrics import calc_SSIM as ssim\n",
    "\n",
    "from ulf_recon_fns import ulfl1recon, applyMask, poissonDiscSigpy, find_lamda_mask, find_iter_mask, coil_combine, awgn\n",
    "#import display_fns\n",
    "\n",
    "import automap_fns #import automapRecon\n",
    "import unrolling_fns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "6a76a2a0-b114-4b5a-aaae-bceeef46bae2",
   "metadata": {},
   "outputs": [],
   "source": [
    "#font = {'family' : 'normal',\n",
    "#        'weight' : 'bold',\n",
    "#        'size'   : 22}\n",
    "\n",
    "#plt.rc('font', **font)\n",
    "plt.rc('font', size=9)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f818cac8-9ec5-4814-8739-9a2d0e39b33b",
   "metadata": {},
   "source": [
    "### Set mask parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "861ed2b3-ffc7-4171-8098-59f1680527ef",
   "metadata": {},
   "source": [
    "## Importing fastMRI data and downsampling"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25676454-d49b-4296-8cc9-aa56c8060c56",
   "metadata": {},
   "source": [
    "The downsampled FastMRI data was saved externally to save processing time, it can be regenerated by uncommenting the following block of code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d4392e05-2fd8-41e2-bbb7-fc9078a84f59",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# loading previously downsampled fast MRI data\n",
    "#inputfile = 'data/fastmri_75_64_25_train_1.npz'\n",
    "inputfile = 'data/fastmri/fastmri_75_64_25_val_all.npz'\n",
    "data = np.load(inputfile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "64a0805e-a2ac-43ac-8059-8db7f580dd2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "ksp_fs_sc_all = data['ksp_fs_sc_all']\n",
    "ksp_us_sc_all = data['ksp_us_sc_all']\n",
    "img_gt_sc_all = data['img_gt_sc_all']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "0c94d5b5-c929-4860-a522-ea867296d7a8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1, 75, 64, 25, 1378)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ksp_fs_sc_all.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f4cf1ead-c25d-46d2-9d74-d9262f2d2be7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(75, 64, 25, 1378)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "img_gt_sc_all.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "5cd02c97-d11b-4a61-84a1-3a4038176bd8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# comment out unneeded code block below\n",
    "\n",
    "#quick run to generate images at high/low SNR\n",
    "noise_dB = np.array([60,70,80]) \n",
    "n_ex_cs = 3\n",
    "n_ex_test = 3\n",
    "\n",
    "#test run to generate manuscript plots\n",
    "# noise_dB = np.arange(50,92,2)\n",
    "# n_ex_cs = 3\n",
    "# n_ex_test = 3\n",
    "\n",
    "# #long run to generate manuscript plots\n",
    "# noise_dB = np.arange(50,92,2)\n",
    "# n_ex_cs = 20\n",
    "# n_ex_test = 100\n",
    "\n",
    "noise_dB_zero = noise_dB - 40"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "f245945c-597d-4f89-a5e4-3c55ad77b6f3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([20, 30, 40])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "noise_dB_zero"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c094e8ef-8233-44f9-99bf-67750b4f07af",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "#loss_fn = 'perceptual'\n",
    "loss_fn = 'mse'\n",
    "# defining mask\n",
    "R=2\n",
    "calib_size=10\n",
    "in_seed = 0\n",
    "imSize = img_gt_sc_all[:,:,:,0].shape\n",
    "mask = poissonDiscSigpy(imSize,R, in_seed,calib_size)\n",
    "\n",
    "mask_R2 = mask\n",
    "\n",
    "mps = np.ones(ksp_fs_sc_all[:,:,:,:,0].shape,dtype=complex) #single coil"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ed57f3b3-c309-413b-adde-3970ad0f833b",
   "metadata": {},
   "outputs": [],
   "source": [
    "if R==2:\n",
    "    unroll_model_pth = 'trained_models/Unrolling/DisCorNet_mse_R2.pth'\n",
    "    model_real_dir = 'trained_models/AUTOMAP/train_25x75_fmri-hcp_r2_real_noise_bkgd_20240318/checkpoint/model.keras' #fast mri model\n",
    "    model_imag_dir = 'trained_models/AUTOMAP/train_25x75_fmri-hcp_r2_imag_noise_bkgd_20240318/checkpoint/model.keras'\n",
    "    \n",
    "elif R==4:\n",
    "    unroll_model_pth = 'trained_models/Unrolling/DisCorNet_mse_R4.pth'\n",
    "    model_real_dir = 'trained_models/AUTOMAP/train_25x75_fmri-hcp_r4_real_noise_bkgd_20240328/checkpoint/model.keras' #fast mri model\n",
    "    model_imag_dir = 'trained_models/AUTOMAP/train_25x75_fmri-hcp_r4_imag_noise_bkgd_20240328/checkpoint/model.keras'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "dd744ca5-bd93-4d4d-a93f-e55b35e37a03",
   "metadata": {},
   "outputs": [],
   "source": [
    "#mdic = {\"mask\":mask}\n",
    "#savemat(\"mask_R2_PD.mat\",mdic)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "338a191e-f33e-4b68-9569-130887771f1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "lamda_opt_vals = np.zeros((noise_dB.size,n_ex_cs))\n",
    "iter_opt_vals = np.zeros((noise_dB.size,n_ex_cs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "cb9608af-dbf8-42d8-87a0-112472a4f4db",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/david/repos/Low_Field_MRI_Recon/metrics.py:39: RuntimeWarning: invalid value encountered in divide\n",
      "  im2_mag_uint8 = (I_pred * 255 / np.max(I_pred)).astype('uint8')\n",
      "/home/david/repos/Low_Field_MRI_Recon/metrics.py:39: RuntimeWarning: invalid value encountered in cast\n",
      "  im2_mag_uint8 = (I_pred * 255 / np.max(I_pred)).astype('uint8')\n"
     ]
    }
   ],
   "source": [
    "for j, noise_val in enumerate(noise_dB):\n",
    "    for k in range(n_ex_cs):\n",
    "        kspTestScaled = ksp_fs_sc_all[:,:,:,:,k]\n",
    "        volumeGTscaled = img_gt_sc_all[:,:,:,k]\n",
    "        \n",
    "        kspInput = awgn(kspTestScaled,noise_val,L=1)\n",
    "        uskspace = applyMask(kspInput,mask)\n",
    "\n",
    "        lamda_opt_vals[j,k] = find_lamda_mask(uskspace, volumeGTscaled, mps,calib_size, show_plot=False)\n",
    "        iter_opt_vals[j,k] = find_iter_mask(uskspace,volumeGTscaled, lamda_opt_vals[j,k], mps, show_plot=False)\n",
    "\n",
    "lamda_opt = lamda_opt_vals.mean(axis=1)\n",
    "iter_opt = iter_opt_vals.mean(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "e8bea416-8079-4195-8e91-e05cca92d6ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "nrmse_cs_vals = np.zeros((noise_dB.size,n_ex_test))\n",
    "nrmse_ifft_vals = np.zeros((noise_dB.size,n_ex_test))\n",
    "nrmse_amap_vals = np.zeros((noise_dB.size,n_ex_test))\n",
    "nrmse_unroll_vals = np.zeros((noise_dB.size,n_ex_test))\n",
    "\n",
    "ssim_cs_vals = np.zeros((noise_dB.size,n_ex_test))\n",
    "ssim_ifft_vals = np.zeros((noise_dB.size,n_ex_test))\n",
    "ssim_amap_vals = np.zeros((noise_dB.size,n_ex_test))\n",
    "ssim_unroll_vals = np.zeros((noise_dB.size,n_ex_test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "c0ea5bde-7d9f-45d0-af2e-8cca3d608a1d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1, 75, 64, 25)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-11-14 15:53:25.382293: E external/local_xla/xla/stream_executor/cuda/cuda_driver.cc:282] failed call to cuInit: CUDA_ERROR_NO_DEVICE: no CUDA-capable device is detected\n",
      "2024-11-14 15:53:25.382328: I external/local_xla/xla/stream_executor/cuda/cuda_diagnostics.cc:134] retrieving CUDA diagnostic information for host: autobeast\n",
      "2024-11-14 15:53:25.382334: I external/local_xla/xla/stream_executor/cuda/cuda_diagnostics.cc:141] hostname: autobeast\n",
      "2024-11-14 15:53:25.382391: I external/local_xla/xla/stream_executor/cuda/cuda_diagnostics.cc:165] libcuda reported version is: 550.54.14\n",
      "2024-11-14 15:53:25.382409: I external/local_xla/xla/stream_executor/cuda/cuda_diagnostics.cc:169] kernel reported version is: 550.54.14\n",
      "2024-11-14 15:53:25.382414: I external/local_xla/xla/stream_executor/cuda/cuda_diagnostics.cc:248] kernel version seems to match DSO: 550.54.14\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1, 75, 64, 25)\n",
      "(1, 75, 64, 25)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/david/repos/Low_Field_MRI_Recon/metrics.py:38: RuntimeWarning: invalid value encountered in divide\n",
      "  im1_mag_uint8 = (I_true * 255 / np.max(I_true)).astype('uint8')\n",
      "/home/david/repos/Low_Field_MRI_Recon/metrics.py:38: RuntimeWarning: invalid value encountered in cast\n",
      "  im1_mag_uint8 = (I_true * 255 / np.max(I_true)).astype('uint8')\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1, 75, 64, 25)\n",
      "(1, 75, 64, 25)\n",
      "(1, 75, 64, 25)\n",
      "(1, 75, 64, 25)\n",
      "(1, 75, 64, 25)\n",
      "(1, 75, 64, 25)\n"
     ]
    }
   ],
   "source": [
    "for j, noise_val in enumerate(noise_dB):\n",
    "    for k in range(n_ex_test):\n",
    "        kspTestScaled = ksp_fs_sc_all[:,:,:,:,k+n_ex_cs]\n",
    "        volumeGTscaled = img_gt_sc_all[:,:,:,k+n_ex_cs]\n",
    "        \n",
    "        kspInput = awgn(kspTestScaled,noise_val,L=1)\n",
    "        uskspace = applyMask(kspInput,mask)\n",
    "\n",
    "    \n",
    "        volifftrecon_mc = sp.ifft(uskspace,axes=[1,2,3])\n",
    "        volifftrecon = coil_combine(volifftrecon_mc,mps)# np.sqrt(np.sum(np.abs(volifftrecon_mc)**2,axis=0)).reshape(imSize)\n",
    "\n",
    "        volcsrecon = ulfl1recon(uskspace,mask,lamda_opt[j],iter_opt[j],mps)\n",
    "\n",
    "        amap_scale = 1\n",
    "        unroll_scale = 1#1/8\n",
    "        volunrollrecon = unrolling_fns.unrollingRecon(uskspace*unroll_scale,mask,unroll_model_pth)/unroll_scale\n",
    "        #volunrollrecon = np.roll(volunrollrecon,(1,1),axis=(0,2))\n",
    "        volamaprecon, volamaprecon_mc = automap_fns.automapRecon(uskspace*amap_scale,mps,model_real_dir,model_imag_dir)\n",
    "        volamaprecon = volamaprecon/amap_scale\n",
    "        volamaprecon_mc = volamaprecon_mc/amap_scale\n",
    "\n",
    "        mask_metrics = np.ma.getmask(np.ma.masked_less(abs(volumeGTscaled),0.05)) # a mask in the image domain for calculating image metrics (not for masking kspace)\n",
    "\n",
    "        volumeGTscaled[mask_metrics]=0\n",
    "        volcsrecon[mask_metrics]=0\n",
    "        volifftrecon[mask_metrics]=0\n",
    "        volamaprecon[mask_metrics]=0\n",
    "        volunrollrecon[mask_metrics]=0\n",
    "\n",
    "        nrmse_cs_vals[j,k] = nrmse(np.abs(volcsrecon[:,:,1:]),np.abs(volumeGTscaled[:,:,1:]))\n",
    "        nrmse_ifft_vals[j,k] = nrmse(np.abs(volifftrecon[:,:,1:]),np.abs(volumeGTscaled[:,:,1:]))\n",
    "        nrmse_amap_vals[j,k] = nrmse(np.abs(volamaprecon[:,:,1:]),np.abs(volumeGTscaled[:,:,1:]))\n",
    "        nrmse_unroll_vals[j,k] = nrmse(np.abs(volunrollrecon[:,:,1:]),np.abs(volumeGTscaled[:,:,1:]))\n",
    "\n",
    "        ssim_cs_vals[j,k] = ssim(np.abs(volcsrecon[:,:,1:]),np.abs(volumeGTscaled[:,:,1:]))\n",
    "        ssim_ifft_vals[j,k] = ssim(np.abs(volifftrecon[:,:,1:]),np.abs(volumeGTscaled[:,:,1:]))\n",
    "        ssim_amap_vals[j,k] = ssim(np.abs(volamaprecon[:,:,1:]),np.abs(volumeGTscaled[:,:,1:]))\n",
    "        ssim_unroll_vals[j,k] = ssim(np.abs(volunrollrecon[:,:,1:]),np.abs(volumeGTscaled[:,:,1:]))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "b762d5cd-6ae4-40de-9830-16198e5304cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "nrmse_cs = nrmse_cs_vals.mean(axis=1)\n",
    "nrmse_ifft = nrmse_ifft_vals.mean(axis=1)\n",
    "nrmse_amap = nrmse_amap_vals.mean(axis=1)\n",
    "nrmse_unroll = nrmse_unroll_vals.mean(axis=1)\n",
    "\n",
    "nrmse_cs_sd = nrmse_cs_vals.std(axis=1)\n",
    "nrmse_ifft_sd = nrmse_ifft_vals.std(axis=1)\n",
    "nrmse_amap_sd = nrmse_amap_vals.std(axis=1)\n",
    "nrmse_unroll_sd = nrmse_unroll_vals.std(axis=1)\n",
    "\n",
    "ssim_cs = ssim_cs_vals.mean(axis=1)\n",
    "ssim_ifft = ssim_ifft_vals.mean(axis=1)\n",
    "ssim_amap = ssim_amap_vals.mean(axis=1)\n",
    "ssim_unroll = ssim_unroll_vals.mean(axis=1)\n",
    "\n",
    "ssim_cs_sd = ssim_cs_vals.std(axis=1)\n",
    "ssim_ifft_sd = ssim_ifft_vals.std(axis=1)\n",
    "ssim_amap_sd = ssim_amap_vals.std(axis=1)\n",
    "ssim_unroll_sd = ssim_unroll_vals.std(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "46c5be34-8140-4eb7-9594-4f675338ee43",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhAAAAH9CAYAAABLMkVIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAsTElEQVR4nO3df2xVd/3H8RcbcLP23t57S3HC7m2zZLBhBC6GGA1xzDg3F12GZDKz1CGxDAItRpxxBTrC1t7qNMFl8E1NhGj8w43GbM6YTMZ/dCG4TO9iSoQxlZYfHRuUc1ugl177+f5helntz3d7b8tpn4/kJPee+z6n73s5953X/cWZ5ZxzAgAAMLhtqhsAAAD+Q4AAAABmBAgAAGBGgAAAAGYECAAAYEaAAAAAZgQIAABgRoAAAABmBAgAAGBGgAAAAGamALFv3z6tXLlSgUBAa9asGbE2nU7rySefVElJie6880698MILE+kTwDTGbAH8Z7aleOHChdq1a5eOHDmis2fPjlhbU1Ojy5cvq62tTRcvXtSDDz6oiooKPfXUUxNqGMD0w2wB/McUINauXStJSqVSIz7Jr127pldeeUVvv/22IpGIIpGIampqdODAgRGf5JlMRplMJne9r69Ply9f1rx58zRr1ixLq4BvOOfU1dWlhQsX6rbbZuaniswWIP8KPVtMAWKsTp48qRs3biiRSOTWJRIJJZPJEbdrbGzUnj17CtEScMtrb29XLBab6jZuacwWwK5Qs6UgAaK7u1vFxcWaPfvm7iORiLq6ukbcrra2Vtu3b89d9zxP5eXlam9vV0lJSSFaBaZcOp1WPB5XKBSa6lZuecwWYOwKPVsKEiCCwaCuXbumbDabe6J7njfqnQgEAgoEAoPWl5SU8CTHtMdb6aNjtgB2hZotBfnA9d5779WcOXP03nvv5dalUiktXbq0EH8OwAzBbAFuHaYAkc1m1dPTo2w2q76+PvX09OjGjRuD6oqKivTEE0+orq5Onufp/fff18svv6yqqqq8NQ5g+mC2AP5jChD19fW644471NDQoD/+8Y+644479NBDD0mSHnnkkQFfZNq3b5/C4bBisZhWrVql733ve/zMCsCQmC2A/8xyzrmpbmI46XRa4XBYnufxOSWmLY7zycdjjpmg0Mf5zPzROQAAmBACBAAAMCNAAAAAMwIEAAAwI0AAAAAzAgQAADAjQAAAADMCBAAAMCNAAAAAMwIEAAAwI0AAAAAzAgQAADAjQAAAADMCBAAAMCNAAAAAMwIEAAAwI0AAAAAzAgQAADAjQAAAADMCBAAAMCNAAAAAMwIEAAAwI0AAAAAzAgQAADAjQAAAADMCBAAAMCNAAAAAMwIEAAAwI0AAAAAzAgQAADAjQAAAADMCBAAAMDMHiN7eXlVXVysajaq0tFQ1NTXKZrND1p47d05r1qzRvHnzVFZWpnXr1umjjz6acNMAph9mC+Av5gBRX1+vlpYWnThxQq2trTp69KiSyeSQtVu3bpUknTlzRv/617/U09Ojbdu2TaxjANMSswXwF3OAOHjwoHbt2qUFCxZowYIF2rlzpw4cODBk7T//+U+tW7dOwWBQoVBITzzxhP7+978Pu+9MJqN0Oj1gATAzMFsAfzEFiM7OTp09e1aJRCK3LpFIqK2tTZ7nDarfvn27mpub5Xmerly5ot/97nd69NFHh91/Y2OjwuFwbonH45b2APgUswXwH1OA6O7uliRFIpHcuv7LXV1dg+pXrVqlixcv5j7T7OzsVG1t7bD7r62tled5uaW9vd3SHgCfYrYA/mMKEMFgUJIGvCLovxwKhQbU9vX16atf/apWrVql7u5udXd3a9WqVXrooYeG3X8gEFBJScmABcD0x2wB/McUIKLRqGKxmFKpVG5dKpVSPB5XOBweUHv58mWdOXNG27ZtU1FRkYqKilRTU6Pjx4/r448/zkvzAKYHZgvgP+YvUW7YsEENDQ3q6OhQR0eHksmkqqqqBtWVlZXpnnvu0f79+9XT06Oenh7t379fsVhMZWVleWkewPTBbAH8ZbZ1g7q6Ol26dElLliyRJFVWVmrHjh2SpM2bN0uSmpqaJEl/+MMf9IMf/EB33XWX+vr6tGLFCr3xxhv56h3ANMJsAfxllnPOTXUTw0mn0wqHw/I8j88sMW1xnE8+HnPMBIU+zvmvrAEAgBkBAgAAmBEgAACAGQECAACYESAAAIAZAQIAAJgRIAAAgBkBAgAAmBEgAACAGQECAACYESAAAIAZAQIAAJgRIAAAgBkBAgAAmBEgAACAGQECAACYESAAAIAZAQIAAJgRIAAAgBkBAgAAmBEgAACAGQECAACYESAAAIAZAQIAAJgRIAAAgBkBAgAAmBEgAACAGQECAACYESAAAIAZAQIAAJgRIAAAgBkBAgAAmJkDRG9vr6qrqxWNRlVaWqqamhpls9lh69944w0lEgkVFxdr4cKFampqmlDDAKYnZgvgL+YAUV9fr5aWFp04cUKtra06evSoksnkkLVvvvmmtmzZol/84hdKp9NqbW3VAw88MNGeAUxDzBbAZ5xRLBZzzc3NueuHDh1y5eXlQ9auXLnS/fKXvxzzvnt6epznebmlvb3dSXKe51nbBHzD8zyOc8dsAfKt0LPF9A5EZ2enzp49q0QikVuXSCTU1tYmz/MG1F69elXvvvuuzp07p8WLF+vTn/60vvWtb+nChQvD7r+xsVHhcDi3xONxS3sAfIrZAviPKUB0d3dLkiKRSG5d/+Wurq4BtZ2dnXLO6fXXX9dbb72l06dPKxAIqLKyctj919bWyvO83NLe3m5pD4BPMVsA/5ltKQ4Gg5Ikz/NUVlaWuyxJoVBoyNpt27apoqJCkrRnzx4tWrRIV69eVXFx8aD9BwIBBQIB410A4HfMFsB/TO9ARKNRxWIxpVKp3LpUKqV4PK5wODygNhKJqLy8fMj9OOfsnQKYtpgtgP+Yf4WxYcMGNTQ0qKOjQx0dHUomk6qqqhqy9umnn9bLL7+sc+fO6fr163r++ef1la98JfcKAgD6MVsAfzF9hCFJdXV1unTpkpYsWSJJqqys1I4dOyRJmzdvlqTc77GfffZZXb58WcuXL5ckffnLX9Zvf/vbvDQOYHphtgD+Msvdwu/5pdNphcNheZ6nkpKSqW4HKAiO88nHY46ZoNDHOf+VNQAAMCNAAAAAMwIEAAAwI0AAAAAzAgQAADAjQAAAADMCBAAAMCNAAAAAMwIEAAAwI0AAAAAzAgQAADAjQAAAADMCBAAAMCNAAAAAMwIEAAAwI0AAAAAzAgQAADAjQAAAADMCBAAAMCNAAAAAMwIEAAAwI0AAAAAzAgQAADAjQAAAADMCBAAAMCNAAAAAMwIEAAAwI0AAAAAzAgQAADAjQAAAADMCBAAAMCNAAAAAM3OA6O3tVXV1taLRqEpLS1VTU6NsNjviNtevX9c999yjSCQy3j4BTHPMFsBfzAGivr5eLS0tOnHihFpbW3X06FElk8kRt3nuuedUUVEx7iYBTH/MFsBfzAHi4MGD2rVrlxYsWKAFCxZo586dOnDgwLD17777rt588039+Mc/HnXfmUxG6XR6wAJgZmC2AP5iChCdnZ06e/asEolEbl0ikVBbW5s8zxtUn81mtXHjRu3fv19z584ddf+NjY0Kh8O5JR6PW9oD4FPMFsB/TAGiu7tbkgZ83th/uaura1D9z372M61YsUL333//mPZfW1srz/NyS3t7u6U9AD7FbAH8Z7alOBgMSpI8z1NZWVnusiSFQqEBtadPn1ZTU5P+9re/jXn/gUBAgUDA0hKAaYDZAviP6R2IaDSqWCymVCqVW5dKpRSPxxUOhwfUtrS06MMPP9TixYtVVlamxx57TOl0WmVlZTp+/HhemgcwPTBbAP8xvQMhSRs2bFBDQ4NWrVolSUomk6qqqhpUt27dOj344IO568eOHVNVVZVSqZQ+9alPTaBlANMRswXwF3OAqKur06VLl7RkyRJJUmVlpXbs2CFJ2rx5sySpqalJRUVFKioqym03f/58zZo1S7FYLB99A5hmmC2Av8xyzrmpbmI46XRa4XBYnueppKRkqtsBCoLjfPLxmGMmKPRxzn9lDQAAzAgQAADAjAABAADMCBAAAMCMAAEAAMwIEAAAwIwAAQAAzAgQAADAjAABAADMCBAAAMCMAAEAAMwIEAAAwIwAAQAAzAgQAADAjAABAADMCBAAAMCMAAEAAMwIEAAAwIwAAQAAzAgQAADAjAABAADMCBAAAMCMAAEAAMwIEAAAwIwAAQAAzAgQAADAjAABAADMCBAAAMCMAAEAAMwIEAAAwIwAAQAAzAgQAADAzBwgent7VV1drWg0qtLSUtXU1CibzQ6qy2Qy2rhxo+6++26FQiHdd999OnjwYF6aBjD9MFsAfzEHiPr6erW0tOjEiRNqbW3V0aNHlUwmB9Vls1ktWLBAR44cUTqd1q9//Wv98Ic/1OHDh/PSOIDphdkC+Mss55yzbBCPx7V37149/vjjkqTm5mY988wzOnPmzKjbrl27Vp/97Gf1/PPPD3l7JpNRJpPJXU+n04rH4/I8TyUlJZY2Ad9Ip9MKh8Mz/jhntgD5VejZYnoHorOzU2fPnlUikcitSyQSamtrk+d5I27b09Ojv/zlL1q2bNmwNY2NjQqHw7klHo9b2gPgU8wWwH9MAaK7u1uSFIlEcuv6L3d1dQ27nXNOVVVVWrRokdauXTtsXW1trTzPyy3t7e2W9gD4FLMF8J/ZluJgMChJ8jxPZWVlucuSFAqFhtzGOactW7bo5MmTOnLkiG67bfjMEggEFAgELC0BmAaYLYD/mN6BiEajisViSqVSuXWpVErxeFzhcHhQvXNOW7du1fHjx3X48OEhawCA2QL4j/lXGBs2bFBDQ4M6OjrU0dGhZDKpqqqqIWurq6v19ttv66233lI0Gp1wswCmL2YL4C+mjzAkqa6uTpcuXdKSJUskSZWVldqxY4ckafPmzZKkpqYmnTlzRv/3f/+nQCCgioqK3PaVlZVqamrKR+8AphFmC+Av5p9xTiZ+3oaZgON88vGYYya4pX7GCQAAIBEgAADAOBAgAACAGQECAACYESAAAIAZAQIAAJgRIAAAgBkBAgAAmBEgAACAGQECAACYESAAAIAZAQIAAJgRIAAAgBkBAgAAmBEgAACAGQECAACYESAAAIAZAQIAAJgRIAAAgBkBAgAAmBEgAACAGQECAACYESAAAIAZAQIAAJgRIAAAgBkBAgAAmBEgAACAGQECAACYESAAAIAZAQIAAJgRIAAAgBkBAgAAmJkDRG9vr6qrqxWNRlVaWqqamhpls9kJ1wKY2ZgtgL+YA0R9fb1aWlp04sQJtba26ujRo0omkxOuBTCzMVsAf5nlnHOWDeLxuPbu3avHH39cktTc3KxnnnlGZ86cmVCtJGUyGWUymdx1z/NUXl6u9vZ2lZSUWNoEfCOdTisej+vKlSsKh8NT3c6UYbYA+VXw2eIMLl++7CS5999/P7fu1KlTTpK7cuXKuGv77d6920liYZmRywcffGB5Ok4rzBYWlsIthZots2XQ3d0tSYpEIrl1/Ze7uroGJBxLbb/a2lpt3749d/3KlSuqqKhQW1ubr16Z9ac+P7668Wvvfu1buvlquLS0dKpbmTLMlrHx83Hu19792rdU+NliChDBYDDXVFlZWe6yJIVCoXHX9gsEAgoEAoPWh8Nh3/3DSVJJSYkv+5b827tf+5ak226buT+KYrbY+Pk492vvfu1bKtxsMe01Go0qFosplUrl1qVSKcXj8UEp3lILYGZjtgD+Y44lGzZsUENDgzo6OtTR0aFkMqmqqqoJ1wKY2ZgtgL+YPsKQpLq6Ol26dElLliyRJFVWVmrHjh2SpM2bN0uSmpqaRq0di0AgoN27dw/51uOtzK99S/7t3a99S/7uPZ+YLaPza9+Sf3v3a99S4Xs3/4wTAABg5n5rCwAAjBsBAgAAmBEgAACAGQECAACYESAAAIDZpAeIfJ6ydzJP6TvWv5XJZLRx40bdfffdCoVCuu+++3Tw4MEBNQ888IACgYCCwWBuOX/+fEH6tvQuSd/97nc1d+7cAb0dO3ZsXPuazL4/2W8wGNScOXO0bNmyMd+vfNq3b59WrlypQCCgNWvWjFibTqf15JNPqqSkRHfeeadeeOEF0+24idnCbClE38yWERTkDBsjeO6559zy5cvd+fPn3fnz593y5cvdnj17xlVr2ddk9d3d3e3q6urc6dOnXV9fnzt27JiLRCLuz3/+c65m9erVbu/evQXpcyK9O+fc+vXr3fe///287Gsy+/5fS5cudQ0NDbnro92vfPr973/vXnvtNbd161b32GOPjVj71FNPuYcffth1dna6kydPung87n7zm9+M+XbcxGxhthSi7//FbLlp0gNELBZzzc3NueuHDh1y5eXl46q17GuiJvK3vvnNb7q6urrc9cl+klt6H+3J4IfH/Pjx4+722293586dy62bzCd5v927d4/4JL969aqbO3eue+edd3LrXnzxRXf//feP6XYMxGxhtowVsyU/s2VSA0Q+T9k7nlP6Tkbf/+v69evurrvuGnCwrl692s2bN89Fo1GXSCQK+orS2vv69etdNBp10WjUfeYzn3E///nP3X/+859x7Wsy+/6kp59+2n3jG98YsG6k+1Uooz3J//rXvzpJrre3N7fu8OHDLhKJjOl23MRs+S9mS/77/iRmy0CT+h2I0U7Da6m17Guixvu3nHOqqqrSokWLtHbt2tz6xsZGffDBB/rwww/1k5/8RDU1NXrttdfy2vN4e9+2bZtOnjypjz76SAcOHNBLL72kl156aVz7msy++129elWvvPLKoPMijHS/pkp3d7eKi4s1e/bN/1E+Eonk7t9ot+MmZst/MVvy33c/ZstgkxogPnka3n5jOWXvULWWfU1m3/2cc9qyZYtOnjyp119/fcDpVL/4xS8qHA5rzpw5evjhh7Vp0ya9+uqree15vL1/7nOf0/z583X77bfrC1/4gp599tlcb7f6Yy5Jzc3NKioq0te//vUB60e6X1MlGAzq2rVrA7685Xle7v6NdjtuYrb8F7Ml/333Y7YMNqkBIp+n7J3MU/pa/5ZzTlu3btXx48d1+PDhUfsp1LnapYmf+viTvd3Kj3m/X/3qV1q/fv2AZD2UQj7mY3Xvvfdqzpw5eu+993LrUqmUli5dOqbbcROzZWjMlsGYLXmcLaYPPPKgrq7OrVixwl24cMFduHDBrVixYthvv45Wa9nXZPa9ZcsWt2zZMvfxxx8Puq2zs9P96U9/clevXnXZbNYdOXLEhcNhd+jQoYL0be391VdfdZ7nub6+PvfOO++4iooK9+KLL45rX5PZt3PO/eMf/3CzZs1yp06dMt+vfOrt7XXXr193O3fudI8++qi7fv26y2QyQ9Z+5zvfcY888oi7cuWKO3XqlCsvLx/wufVot+MmZguzpRB9O8dsGc6kB4gbN264LVu2uEgk4iKRiKuurs59kWPTpk1u06ZNY6ody+1T0fe///1vJ8kFAgFXXFycW/pvv3jxovv85z/vQqGQC4VCbunSpe7AgQMF6dnau3POfelLX3LhcNgVFxe7xYsXu5/+9KcDvhB0Kz7m/X70ox8N+y3i0e5XPu3evdtJGrCsXr3aOefc1772tQE/AfM8z3372992wWDQzZ8/f9AQG+123MRsYbYUom/nmC3D4XTeAADAbOo/rAEAAL5DgAAAAGYECAAAYEaAAAAAZgQIAABgRoAAAABmBAgAAGBGgAAAAGYECAAAYEaAAAAAZgQIAABgRoAAAABmBAgAAGBGgAAAAGYECAAAYEaAAAAAZgQIAABgZgoQ+/bt08qVKxUIBLRmzZoRa9PptJ588kmVlJTozjvv1AsvvDCRPgFMY8wWwH9mW4oXLlyoXbt26ciRIzp79uyItTU1Nbp8+bLa2tp08eJFPfjgg6qoqNBTTz01oYYBTD/MFsB/TAFi7dq1kqRUKjXik/zatWt65ZVX9PbbbysSiSgSiaimpkYHDhwY8UmeyWSUyWRy1/v6+nT58mXNmzdPs2bNsrQK+IZzTl1dXVq4cKFuu21mfqrIbAHyr9CzxRQgxurkyZO6ceOGEolEbl0ikVAymRxxu8bGRu3Zs6cQLQG3vPb2dsVisalu45bGbAHsCjVbChIguru7VVxcrNmzb+4+Eomoq6trxO1qa2u1ffv23HXP81ReXq729naVlJQUolVgyqXTacXjcYVCoalu5ZbHbAHGrtCzpSABIhgM6tq1a8pms7knuud5o96JQCCgQCAwaH1JSQlPckx7vJU+OmYLYFeo2VKQD1zvvfdezZkzR++9915uXSqV0tKlSwvx5wDMEMwW4NZhChDZbFY9PT3KZrPq6+tTT0+Pbty4MaiuqKhITzzxhOrq6uR5nt5//329/PLLqqqqylvjAKYPZgvgP6YAUV9frzvuuEMNDQ364x//qDvuuEMPPfSQJOmRRx4Z8EWmffv2KRwOKxaLadWqVfre977Hz6wADInZAvjPLOecm+omhpNOpxUOh+V5Hp9TYtriOJ98POaYCQp9nM/MH50DAIAJIUAAAAAzAgQAADAjQAAAADMCBAAAMCNAAAAAMwIEAAAwI0AAAAAzAgQAADAjQAAAADMCBAAAMCNAAAAAMwIEAAAwI0AAAAAzAgQAADAjQAAAADMCBAAAMCNAAAAAMwIEAAAwI0AAAAAzAgQAADAjQAAAADMCBAAAMCNAAAAAMwIEAAAwI0AAAAAzAgQAADAjQAAAADMCBAAAMCNAAAAAMwIEAAAwI0AAAAAzc4Do7e1VdXW1otGoSktLVVNTo2w2O2TtuXPntGbNGs2bN09lZWVat26dPvroowk3DWD6YbYA/mIOEPX19WppadGJEyfU2tqqo0ePKplMDlm7detWSdKZM2f0r3/9Sz09Pdq2bdvEOgYwLTFbAH8xB4iDBw9q165dWrBggRYsWKCdO3fqwIEDQ9b+85//1Lp16xQMBhUKhfTEE0/o73//+7D7zmQySqfTAxYAMwOzBfAXU4Do7OzU2bNnlUgkcusSiYTa2trked6g+u3bt6u5uVme5+nKlSv63e9+p0cffXTY/Tc2NiocDueWeDxuaQ+ATzFbAP8xBYju7m5JUiQSya3rv9zV1TWoftWqVbp48WLuM83Ozk7V1tYOu//a2lp5npdb2tvbLe0B8ClmC+A/pgARDAYlacArgv7LoVBoQG1fX5+++tWvatWqVeru7lZ3d7dWrVqlhx56aNj9BwIBlZSUDFgATH/MFsB/TAEiGo0qFosplUrl1qVSKcXjcYXD4QG1ly9f1pkzZ7Rt2zYVFRWpqKhINTU1On78uD7++OO8NA9gemC2AP5j/hLlhg0b1NDQoI6ODnV0dCiZTKqqqmpQXVlZme655x7t379fPT096unp0f79+xWLxVRWVpaX5gFMH8wWwF9mWzeoq6vTpUuXtGTJEklSZWWlduzYIUnavHmzJKmpqUmS9Ic//EE/+MEPdNddd6mvr08rVqzQG2+8ka/eAUwjzBbAX2Y559xUNzGcdDqtcDgsz/P4zBLTFsf55OMxx0xQ6OOc/8oaAACYESAAAIAZAQIAAJgRIAAAgBkBAgAAmBEgAACAGQECAACYESAAAIAZAQIAAJgRIAAAgBkBAgAAmBEgAACAGQECAACYESAAAIAZAQIAAJgRIAAAgBkBAgAAmBEgAACAGQECAACYESAAAIAZAQIAAJgRIAAAgBkBAgAAmBEgAACAGQECAACYESAAAIAZAQIAAJgRIAAAgBkBAgAAmBEgAACAGQECAACYESAAAICZOUD09vaqurpa0WhUpaWlqqmpUTabHbb+jTfeUCKRUHFxsRYuXKimpqYJNQxgemK2AP5iDhD19fVqaWnRiRMn1NraqqNHjyqZTA5Z++abb2rLli36xS9+oXQ6rdbWVj3wwAMT7RnANMRsAXzGGcViMdfc3Jy7fujQIVdeXj5k7cqVK90vf/nLMe+7p6fHeZ6XW9rb250k53metU3ANzzP4zh3zBYg3wo9W0zvQHR2durs2bNKJBK5dYlEQm1tbfI8b0Dt1atX9e677+rcuXNavHixPv3pT+tb3/qWLly4MOz+GxsbFQ6Hc0s8Hre0B8CnmC2A/5gCRHd3tyQpEonk1vVf7urqGlDb2dkp55xef/11vfXWWzp9+rQCgYAqKyuH3X9tba08z8st7e3tlvYA+BSzBfCf2ZbiYDAoSfI8T2VlZbnLkhQKhYas3bZtmyoqKiRJe/bs0aJFi3T16lUVFxcP2n8gEFAgEDDeBQB+x2wB/Mf0DkQ0GlUsFlMqlcqtS6VSisfjCofDA2ojkYjKy8uH3I9zzt4pgGmL2QL4j/lXGBs2bFBDQ4M6OjrU0dGhZDKpqqqqIWuffvppvfzyyzp37pyuX7+u559/Xl/5yldyryAAoB+zBfAX00cYklRXV6dLly5pyZIlkqTKykrt2LFDkrR582ZJyv0e+9lnn9Xly5e1fPlySdKXv/xl/fa3v81L4wCmF2YL4C+z3C38nl86nVY4HJbneSopKZnqdoCC4DiffDzmmAkKfZzzX1kDAAAzAgQAADAjQAAAADMCBAAAMCNAAAAAMwIEAAAwI0AAAAAzAgQAADAjQAAAADMCBAAAMCNAAAAAMwIEAAAwI0AAAAAzAgQAADAjQAAAADMCBAAAMCNAAAAAMwIEAAAwI0AAAAAzAgQAADAjQAAAADMCBAAAMCNAAAAAMwIEAAAwI0AAAAAzAgQAADAjQAAAADMCBAAAMCNAAAAAMwIEAAAwI0AAAAAzAgQAADAzB4je3l5VV1crGo2qtLRUNTU1ymazI25z/fp13XPPPYpEIuPtE8A0x2wB/MUcIOrr69XS0qITJ06otbVVR48eVTKZHHGb5557ThUVFeNuEsD0x2wB/MUcIA4ePKhdu3ZpwYIFWrBggXbu3KkDBw4MW//uu+/qzTff1I9//ONR953JZJROpwcsAGYGZgvgL6YA0dnZqbNnzyqRSOTWJRIJtbW1yfO8QfXZbFYbN27U/v37NXfu3FH339jYqHA4nFvi8bilPQA+xWwB/McUILq7uyVpwOeN/Ze7uroG1f/sZz/TihUrdP/9949p/7W1tfI8L7e0t7db2gPgU8wWwH9mW4qDwaAkyfM8lZWV5S5LUigUGlB7+vRpNTU16W9/+9uY9x8IBBQIBCwtAZgGmC2A/5jegYhGo4rFYkqlUrl1qVRK8Xhc4XB4QG1LS4s+/PBDLV68WGVlZXrssceUTqdVVlam48eP56V5ANMDswXwH9M7EJK0YcMGNTQ0aNWqVZKkZDKpqqqqQXXr1q3Tgw8+mLt+7NgxVVVVKZVK6VOf+tQEWgYwHTFbAH8xB4i6ujpdunRJS5YskSRVVlZqx44dkqTNmzdLkpqamlRUVKSioqLcdvPnz9esWbMUi8Xy0TeAaYbZAvjLLOecm+omhpNOpxUOh+V5nkpKSqa6HaAgOM4nH485ZoJCH+f8V9YAAMCMAAEAAMwIEAAAwIwAAQAAzAgQAADAjAABAADMCBAAAMCMAAEAAMwIEAAAwIwAAQAAzAgQAADAjAABAADMCBAAAMCMAAEAAMwIEAAAwIwAAQAAzAgQAADAjAABAADMCBAAAMCMAAEAAMwIEAAAwIwAAQAAzAgQAADAjAABAADMCBAAAMCMAAEAAMwIEAAAwIwAAQAAzAgQAADAjAABAADMCBAAAMCMAAEAAMzMAaK3t1fV1dWKRqMqLS1VTU2NstnsoLpMJqONGzfq7rvvVigU0n333aeDBw/mpWkA0w+zBfAXc4Cor69XS0uLTpw4odbWVh09elTJZHJQXTab1YIFC3TkyBGl02n9+te/1g9/+EMdPnw4L40DmF6YLYC/zHLOOcsG8Xhce/fu1eOPPy5Jam5u1jPPPKMzZ86Muu3atWv12c9+Vs8///yQt2cyGWUymdz1dDqteDwuz/NUUlJiaRPwjXQ6rXA4POOPc2YLkF+Fni2mdyA6Ozt19uxZJRKJ3LpEIqG2tjZ5njfitj09PfrLX/6iZcuWDVvT2NiocDicW+LxuKU9AD7FbAH8xxQguru7JUmRSCS3rv9yV1fXsNs551RVVaVFixZp7dq1w9bV1tbK87zc0t7ebmkPgE8xWwD/mW0pDgaDkiTP81RWVpa7LEmhUGjIbZxz2rJli06ePKkjR47ottuGzyyBQECBQMDSEoBpgNkC+I/pHYhoNKpYLKZUKpVbl0qlFI/HFQ6HB9U757R161YdP35chw8fHrIGAJgtgP+Yf4WxYcMGNTQ0qKOjQx0dHUomk6qqqhqytrq6Wm+//bbeeustRaPRCTcLYPpitgD+YvoIQ5Lq6up06dIlLVmyRJJUWVmpHTt2SJI2b94sSWpqatKZM2f0f//3fwoEAqqoqMhtX1lZqaampnz0DmAaYbYA/mL+Gedk4udtmAk4zicfjzlmglvqZ5wAAAASAQIAAIwDAQIAAJgRIAAAgBkBAgAAmBEgAACAGQECAACYESAAAIAZAQIAAJgRIAAAgBkBAgAAmBEgAACAGQECAACYESAAAIAZAQIAAJgRIAAAgBkBAgAAmBEgAACAGQECAACYESAAAIAZAQIAAJgRIAAAgBkBAgAAmBEgAACAGQECAACYESAAAIAZAQIAAJgRIAAAgBkBAgAAmBEgAACAGQECAACYESAAAICZOUD09vaqurpa0WhUpaWlqqmpUTabnXAtgJmN2QL4izlA1NfXq6WlRSdOnFBra6uOHj2qZDI54VoAMxuzBfCXWc45Z9kgHo9r7969evzxxyVJzc3NeuaZZ3TmzJkJ1UpSJpNRJpPJXfc8T+Xl5Wpvb1dJSYmlTcA30um04vG4rly5onA4PNXtTBlmC5BfBZ8tzuDy5ctOknv//fdz606dOuUkuStXroy7tt/u3budJBaWGbl88MEHlqfjtMJsYWEp3FKo2TJbBt3d3ZKkSCSSW9d/uaura0DCsdT2q62t1fbt23PXr1y5ooqKCrW1tfnqlVl/6vPjqxu/9u7XvqWbr4ZLS0unupUpw2wZGz8f537t3a99S4WfLaYAEQwGc02VlZXlLktSKBQad22/QCCgQCAwaH04HPbdP5wklZSU+LJvyb+9+7VvSbrttpn7oyhmi42fj3O/9u7XvqXCzRbTXqPRqGKxmFKpVG5dKpVSPB4flOIttQBmNmYL4D/mWLJhwwY1NDSoo6NDHR0dSiaTqqqqmnAtgJmN2QL4i+kjDEmqq6vTpUuXtGTJEklSZWWlduzYIUnavHmzJKmpqWnU2rEIBALavXv3kG893sr82rfk39792rfk797zidkyOr/2Lfm3d7/2LRW+d/PPOAEAAGbut7YAAMC4ESAAAIAZAQIAAJgRIAAAgBkBAgAAmE16gMjnKXsn85S+Y/1bmUxGGzdu1N13361QKKT77rtPBw8eHFDzwAMPKBAIKBgM5pbz588XpG9L75L03e9+V3Pnzh3Q27Fjx8a1r8ns+5P9BoNBzZkzR8uWLRvz/cqnffv2aeXKlQoEAlqzZs2Itel0Wk8++aRKSkp055136oUXXjDdjpuYLcyWQvTNbBlBQc6wMYLnnnvOLV++3J0/f96dP3/eLV++3O3Zs2dctZZ9TVbf3d3drq6uzp0+fdr19fW5Y8eOuUgk4v785z/nalavXu327t1bkD4n0rtzzq1fv959//vfz8u+JrPv/7V06VLX0NCQuz7a/cqn3//+9+61115zW7dudY899tiItU899ZR7+OGHXWdnpzt58qSLx+PuN7/5zZhvx03MFmZLIfr+X8yWmyY9QMRiMdfc3Jy7fujQIVdeXj6uWsu+Jmoif+ub3/ymq6ury12f7Ce5pffRngx+eMyPHz/ubr/9dnfu3Lncusl8kvfbvXv3iE/yq1evurlz57p33nknt+7FF190999//5hux0DMFmbLWDFb8jNbJjVA5POUveM5pe9k9P2/rl+/7u66664BB+vq1avdvHnzXDQadYlEoqCvKK29r1+/3kWjUReNRt1nPvMZ9/Of/9z95z//Gde+JrPvT3r66afdN77xjQHrRrpfhTLak/yvf/2rk+R6e3tz6w4fPuwikciYbsdNzJb/Yrbkv+9PYrYMNKnfgRjtNLyWWsu+Jmq8f8s5p6qqKi1atEhr167NrW9sbNQHH3ygDz/8UD/5yU9UU1Oj1157La89j7f3bdu26eTJk/roo4904MABvfTSS3rppZfGta/J7Lvf1atX9corrww6L8JI92uqdHd3q7i4WLNn3/wf5SORSO7+jXY7bmK2/BezJf9992O2DDapAeKTp+HtN5ZT9g5Va9nXZPbdzzmnLVu26OTJk3r99dcHnE71i1/8osLhsObMmaOHH35YmzZt0quvvprXnsfb++c+9znNnz9ft99+u77whS/o2WefzfV2qz/mktTc3KyioiJ9/etfH7B+pPs1VYLBoK5duzbgy1ue5+Xu32i34yZmy38xW/Lfdz9my2CTGiDyecreyTylr/VvOee0detWHT9+XIcPHx61n0Kdq12a+KmPP9nbrfyY9/vVr36l9evXD0jWQynkYz5W9957r+bMmaP33nsvty6VSmnp0qVjuh03MVuGxmwZjNmSx9li+sAjD+rq6tyKFSvchQsX3IULF9yKFSuG/fbraLWWfU1m31u2bHHLli1zH3/88aDbOjs73Z/+9Cd39epVl81m3ZEjR1w4HHaHDh0qSN/W3l999VXneZ7r6+tz77zzjquoqHAvvvjiuPY1mX0759w//vEPN2vWLHfq1Cnz/cqn3t5ed/36dbdz50736KOPuuvXr7tMJjNk7Xe+8x33yCOPuCtXrrhTp0658vLyAZ9bj3Y7bmK2MFsK0bdzzJbhTHqAuHHjhtuyZYuLRCIuEom46urq3Bc5Nm3a5DZt2jSm2rHcPhV9//vf/3aSXCAQcMXFxbml//aLFy+6z3/+8y4UCrlQKOSWLl3qDhw4UJCerb0759yXvvQlFw6HXXFxsVu8eLH76U9/OuALQbfiY97vRz/60bDfIh7tfuXT7t27naQBy+rVq51zzn3ta18b8BMwz/Pct7/9bRcMBt38+fMHDbHRbsdNzBZmSyH6do7ZMhxO5w0AAMym/sMaAADgOwQIAABgRoAAAABmBAgAAGBGgAAAAGYECAAAYEaAAAAAZgQIAABgRoAAAABmBAgAAGBGgAAAAGb/D6PKmsGUABYcAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 600x600 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ((ax1, ax3),(ax2, ax4)) = plt.subplots(2,2)\n",
    "fig.set_size_inches(6,6)\n",
    "fig.subplots_adjust(wspace=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "283084f6-0e87-481b-83e3-0f2c1d9bfa1b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.axis.YTick at 0x7f62a0622740>,\n",
       " <matplotlib.axis.YTick at 0x7f62a06214b0>,\n",
       " <matplotlib.axis.YTick at 0x7f62a05a73d0>]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "line1 = ax1.plot(noise_dB_zero,nrmse_ifft, '-')\n",
    "line2 = ax1.plot(noise_dB_zero,nrmse_cs, '-')\n",
    "line3 = ax1.plot(noise_dB_zero,nrmse_amap, '-')\n",
    "line4 = ax1.plot(noise_dB_zero,nrmse_unroll, '-')\n",
    "#ax1.set(title='Comparison of Reconstruction Methods, R='+str(R),\n",
    "ax1.set(ylabel='NRMSE',\n",
    "xlabel='Input SNR (dB)')\n",
    "#ax1.set_xscale('log')\n",
    "ax1.legend(['Zero-filled', 'CS', 'AUTOMAP', 'Unrolled'])\n",
    "#ax1.legend([line1, line2, line3], ['label1', 'label2', 'label3'])\n",
    "#plt.ylim(0,0.085)\n",
    "#plt.xlim(55,90)\n",
    "fill1 = ax1.fill_between(noise_dB_zero,nrmse_ifft-nrmse_ifft_sd,nrmse_ifft+nrmse_ifft_sd,alpha=0.25)\n",
    "fill2 = ax1.fill_between(noise_dB_zero,nrmse_cs-nrmse_cs_sd,nrmse_cs+nrmse_cs_sd,alpha=0.25)\n",
    "fill3 = ax1.fill_between(noise_dB_zero,nrmse_amap-nrmse_amap_sd,nrmse_amap+nrmse_amap_sd,alpha=0.25)\n",
    "fill4 = ax1.fill_between(noise_dB_zero,nrmse_unroll-nrmse_unroll_sd,nrmse_unroll+nrmse_unroll_sd,alpha=0.25)\n",
    "\n",
    "step=10\n",
    "ax1.set_xticks(np.arange(noise_dB_zero[0],noise_dB_zero[-1]+step,step))\n",
    "ax1.set_yticks([0,0.05,0.1])\n",
    "\n",
    "line1 = ax3.plot(noise_dB_zero,ssim_ifft, '-')\n",
    "line2 = ax3.plot(noise_dB_zero,ssim_cs, '-')\n",
    "line3 = ax3.plot(noise_dB_zero,ssim_amap, '-')\n",
    "line4 = ax3.plot(noise_dB_zero,ssim_unroll, '-')\n",
    "#ax3.set(title='Comparison of Reconstruction Methods, R='+str(R),\n",
    "ax3.set(ylabel='SSIM',\n",
    "xlabel='Input SNR (dB)')\n",
    "#ax1.set_xscale('log')\n",
    "ax3.legend(['Zero-filled', 'CS', 'AUTOMAP', 'Unrolled AI'])\n",
    "#ax1.legend([line1, line2, line3], ['label1', 'label2', 'label3'])\n",
    "#plt.ylim(0,0.085)\n",
    "#plt.xlim(55,90)\n",
    "\n",
    "fill1 = ax3.fill_between(noise_dB_zero,ssim_ifft-ssim_ifft_sd,ssim_ifft+ssim_ifft_sd,alpha=0.25)\n",
    "fill2 = ax3.fill_between(noise_dB_zero,ssim_cs-ssim_cs_sd,ssim_cs+ssim_cs_sd,alpha=0.25)\n",
    "fill3 = ax3.fill_between(noise_dB_zero,ssim_amap-ssim_amap_sd,ssim_amap+ssim_amap_sd,alpha=0.25)\n",
    "fill4 = ax3.fill_between(noise_dB_zero,ssim_unroll-ssim_unroll_sd,ssim_unroll+ssim_unroll_sd,alpha=0.25)\n",
    "\n",
    "step=10\n",
    "ax3.set_xticks(np.arange(noise_dB_zero[0],noise_dB_zero[-1]+step,step))\n",
    "ax3.set_yticks([0.8,0.9,1])\n",
    "\n",
    "#plt.show()\n",
    "\n",
    "#fig.savefig('export_figs/fastmri_noise_curves_R'+str(R)+'.pdf',bbox_inches='tight',dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "99d79d88-6e17-4e5b-a9a7-24fa87a00828",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAisAAAH9CAYAAADbITyHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAACa40lEQVR4nOzdeXxU5aH/8c+Zfcss2feFHdlBUHBDFNRqcalaBRdstVQLFq33VgH3vXpvvdefXrVFcKFVrK1araJUrbKKRUBAAYEAIfs2+37O74+TTBKSQAIhG8/79ZpXZs6cOefJZHLyzbNKiqIoCIIgCIIg9FKani6AIAiCIAjCkYiwIgiCIAhCrybCiiAIgiAIvZoIK4IgCIIg9GoirAiCIAiC0KuJsCIIgiAIQq8mwoogCIIgCL2aCCuCIAiCIPRqIqwIgiAIgtCribAiCIIgCEKv1qNhJRqNMm/ePFwuF8nJycyfP59YLNbmvv/v//0/Tj31VIxGI5dddlmr5z0eD7NmzcJut5ORkcHDDz98gksvCIIgCEJ36NGw8sgjj7B69Wp27NjB9u3b+fLLL3nsscfa3Dc7O5vFixdzyy23tPn8/Pnzqa2t5cCBA3z55Zf84Q9/4NVXXz2RxRcEQRAEoRtIPbmQYV5eHr///e+58sorAXjrrbe466672L9/f7uveeCBB9i8eTPvvPNOYlsgEMDlcrFmzRpOPfVUAJ566inef/99/vWvf7U6RjgcJhwOJx7LskxtbS0pKSlIktRF350g9C6KouD1esnOzkajOXlbgP/f//t/LFu2jG+//ZaLLrqoxbXkcB6Ph1/+8pe8//77mM1m5s2bx7333tvhc8myTGlpKUlJSeLaIvRb3XFt0Z2Qo3ZAXV0dJSUljB07NrFt7NixHDhwALfbjcPh6PCxdu7cSSQSaXWs9mppHn/8cR588MFjLbog9GkHDx4kNze3p4vRYxpraVetWkVJSckR921eY1tZWcn5559PQUEBN9xwQ4fOVVpaSl5eXlcUWxB6vRN5bemxsOLz+QBwOp2JbY33vV5vp8KKz+fDarWi0zV9O06nE6/X2+b+99xzD3feeWfisdvtJj8/n4MHD2K32zvxXQhC3+HxeMjLyyMpKamni9KjrrjiCgA2b958xLASCAR44403WLNmDU6nE6fTyfz581myZEm7YeXwWtvGimtxbRH6s+64tvRYWLHZbIAaFFJTUxP3gU5/wzabjUAgQCwWSwQWt9vd7nGMRiNGo7HVdrvdLi4oQr8nmiM6prM1ttB+ra24tggngxN5bemxhmuXy0Vubi6bN29ObNu8eTN5eXmdqlUBGDp0KHq9ni1btrQ41qhRo7qquIIgnGQ6W2MLaq2t2+1O3A4ePNgdRRWEfq9He9nddNNNPProo5SXl1NeXs5jjz3GzTff3Oa+sViMUChELBZDlmVCoRCRSAQAi8XCT3/6U+69917cbje7d+/m2WefbfdYgiAIR9O8xrbRkWpsQa21baxFEbUpgtB1eqwZCODee++lpqaG4cOHA3DdddexcOFCAH75y18C8MILLwDqMOfm1atms5lzzjmHzz//HFB7+M+dO5fc3NxEr/2OdoITBEE4XPMa2wkTJgCixlY4ecixGLGKCqKlpcTKy9GYzSSdf36PladHhy73Fh6PB4fDgdvtFv8JtSMejxONRnu6GEIH6PV6tFptq+3ic66KxWLEYjEeeeQRtm7dyooVK9BoNBgMhlb73nDDDVRXV/PnP/85MRro4Ycf7vA/QuI9F3obJRpF9geIlpcRLS0lWl5OrKKCWEUlsapKYtXVxGpqidfVQTyeeJ2hsJCBH33Y5jG743PeozUrQt/g8/koKSlB5Nq+QZIkcnNzE53YhZaOVEt70UUXcdZZZyVqeEWNrdAXKPE4SjSKEo0S93qJlpURKysjWl6hBpGqKmLV1cRraojV1iJ7PNCR67kkgSMJKSUZ0yk9W6MoalYQ//0cSTweZ/fu3VgsFtLS0sRIkl5OURSqqqoIBAIMHjy4RQ2L+Jx3P/GeC8dKURQ1gESiKNEIciSCXF9PrKycaEVDbUhlFbGahpqQWrU2RPb7O3YCrRbFkYTisiMnO4il2AkkW3C7DNS6tFTZodwWoy7mxRPxkJuUy9PnPN3moUTNitDjotEoiqKQlpaG2Wzu6eIIHZCWlkZxcTHRaLTN5iBBEHqOEoslakEab3I4TLymRm2SKW9WE1JbS6yujnhdHfHaWpRmc/gc8RwGPYozCcVpJ+6yE3KZ8bhM1Lt0VDskKuwy5cYw9bKX+rCb+vAB/NFmIScO1DXcGlQFq7r0fegsEVaEDhE1Kn2H+FkJQvdSZLllAIk03o8gh0LEKiubbs1qQeJ1dcRqa4nX10M7i/i2OpfFhOK0I9tthJ1m/E4TbpeeWoeWSgeUJsWo1Aeoj6lBxBttY+JDT9vH1mv0uIwunCYnTqN6cxgdpFvSmZAx4djfoC4gwoogCIIgtOHwGpDDb3G/v6FPSHUigMTqaonV1BKrUx8r7o71D1EkIMmG7LARdVgI2o14HHrqHFqqnRLldplSa4wqrZ/6qAdP5CAKbRy3jVYgraRVw0ezEJK4NWxzGV2YdWYMGj0GRcGgKBgVBb0iY9YacWSOOf439DiIsCIIXeQPf/gD9913H16vly+//JLrrruOJ598kksuuYRly5bxzDPPJCZBLCws5JlnnuGyyy7r9Hnq6+txuVzs27ePwsLCLv0eBKG/U2QZJRJpHT4izR9HiHs8xGtqidWqwSNWW0O0trFWpB65rh462D9E0WpQ7DbiDhtBhwl/kh63XUuNQ0OlQ6HMFqPEGqJG8eGJ1RJXqts+UKjlQ42kwWFwHDWI2PQ2JElCp9FhkPQYIBFIDIqMQZYxyHEMkTiS7Gt9XkPPNyeLsCL0ecuXL2fu3LmttgeDQaZOnco///nPE16GaDTK7bffzscff8xZZ50FwPbt20/4eQVBaEmJRIi73YcFkaZwIsfiyF4v0doaIjXVDTUhNcTr64jX1qHU16PUeyDU8f4hcYeNiMNMwG7Am6Sjzq6hyqZQ7pA5ZItwyBygLu4jpgTaP1Ck6a6EhN1gP2IAcRqdJBmS0Ejq3K5aSYtBa0Cv0WNAgxEFg6ygl+MYFQWDHEMTiULM37GRQAAhD3gOQaAaqnbC2Fkde90JIMKK0OfNnj2b2bNnt9j25ZdfcsEFF7B48eJOHy8ajaLX6zv1mvLyckKhkJgwTBB6gCLLahNMdQ2x6io8e3cRralu6Jiq9g1R6twobjfUe5Bi8aMfFJDNJmINTTJ+h4H6JA01SVBpkzmUFOWgNUyZwUcEL9D+Mgw0645i09tahY7D+4nYDXa0mqbaDAkJg9aAUWtEr9VjQIdBkTEC+ngcgxxHJ8cgHIKYB+SOfX9q2cLgKVVDSeJrw/1ws84tyQNEWBH6DkVRCEY78YtwnMx6bac7jB46dIirrrqKxx57jHPPPReAVatWsXDhQnbt2kVOTg6PP/44M2fOBGDOnDlotVq8Xi8fffQRjz76KHPmzOE3v/kNf//73wG49NJL+a//+i+sVmur833zzTeceeaZAOTm5pKRkcGePXs61dRzpPKFw2EWLFjAm2++icPhYNGiRZ16PwShv4p7POqombo6YjW11K1aif+zz5G8bTfPNF5JFAlkq5mww0zAYcSdpKM+SaIqSabMFueQLUyxNYhHH6Xd3qiHseqsrTqmHl4b4jA60Gta/yNk0BowaAzNvuoxyGBQ4hgUBb0cg2gIQkGI1UM80roARyLHwV8F7pLWocR/lFE+llRw5kHuRLVGpoc68IuwInRKMBrnlPtWdtv5djx0ARZDxz+mkUiEK6+8kunTp7NgwQIAtm7dylVXXcXbb7/N1KlTWbt2LRdffDFfffUVQ4cOBeDPf/4zf/vb33jjjTcIhULMmzeP4uJitm3bhqIoXHnlldxxxx289NJLrc45btw4tm/fTlFRESUlJTidzk59j0cr36OPPsq6devYtm0bFouFWbN67r8bQehpcjCYmOBMDkeI7N+Pe+VHBDdsgFgcCYhZjfhTrfiS9NTbtdQkQYUtTqktwgFbiFJrhLg2Chx9Vm6T1nTU5hin0YlB23oGZFBH2Oi1eowaoxpGGppqjEgYZBm9LCPFQ2oNRygA0VqIhzveVNNIUSBUrwYQ96GWocRbBvIRRhsZrGDPBXs22HPAkaN+TcoGvQkMFsgRo4EEocvMmzePcDjcIlS8+OKLzJkzh2nTpgFw5plncskll7BixQruvfdeAGbMmMEFF1wAgMlkYvny5XzxxRekpKQA8NhjjzFt2jReeOEFNJquXf/zaOVbvnw5jz/+ONnZ2QDcf//9fPhh29NeC0J/pESjxGpqiFXXIPv9KLJMcPNmvCtXEv7++8R+/oI03pwQ4uPBAWTNkWpEJAwaw1E7pjqMDsy69ueX0krapqaZxpqRxhtaDHIcTTwKsZB6i3gb7oc711TTXDTQFEIODyXRI/SJ0eibwsjhocRo77Eak44SYUXoFLNey46HLujW83XUSy+9xF//+le+/vrrFhPYFRcX8+mnn7J06dLEtlgs1mKmxfz8/MT9qqoqIpFIi5E2AwYMIBwOU11dzSeffJLo0FtQUHDcHWmPVr7S0lIKCgoSzzW/Lwj9VfN+KHF3PShqrYr/yy/xfvIJscpKdT+NRHTMUD6eqOe1lO9RULDqrGTZstsNII3DdI/UxKyRNK0DSPPHkh6t3BBEog1hJOyBWFANI/HjWEstHgVfRVNTTfNQEqw9wgslsKU3BJLDQoklFTQduJ7qjKC3qLUp+oaboXXzd3cTYUXoFEmSOtUs013Wr1/PHXfcwXvvvddqOG9eXh6//vWveeKJJ9p9ffPakrS0NAwGA8XFxWRkZABqoDAajaSmprbZofd4HK182dnZ7N+/n9NOOw2AAwcOdNm5BaG3Ufuh1BCvq0Vp6Agbq6rCu2oVvn/9CyUYBNTJ0WKTx7HvtDye1/6L0lgNAJOzJzNr2Cys+vb/wEpIatOM1ohBY2hx36BVH+s1eohFmsJHLATRuqaakWNpqmlOUSBQ06xDa7NQ4qsARW7/tSZHs0DSLJTYs6Cd5qhWdEbQm5vCSGMw0fa+6zuIsCL0A+Xl5fzkJz/hoYce4rzzzmv1/Ny5c7nwwgu54IILOPvss4nFYmzatAmn08nw4cNb7a/RaJg1axaLFi1ixYoVKIrCwoULuf7667u8Cagj5bv22mt54oknOPPMM7FYLDz00ENdXgZB6ElyMKgOH66uRg6rnUcVRSG8ezfelSsJ/vvfiWAgpycTO3si4QnDeVfexPu+d5BjMnaDnRtH3Mi49HHq8N2GmpBEM81hjwGIx5qaaKIhiHkh2iycHCkwdFTY1zKQeErVjq7eUvU87dEZDwskzUKJsROLlGoNh9WSWEBv7bWhpD19q7SC0IaXXnqJ0tJS7r//fu6///5Wz/t8Pv785z+zePFivvvuOzQaDWPHjuXpp9telAvgf/7nf7jzzjs55ZRTAJg5cyb/9V//dULKP27cuCOWb/HixVRWVjJy5EjsdjuLFy/mgw8+OCFlEYTuokSj6oRrVdUtFt9TYjECGzbg/fhjIsXFie3xoUXEzpmIPHQAB5Rq/uB/mwOhMgAmZk7kuuHXkWRIItOaSa4tV23ikeVm/UXC6kiaRC1J6MidTjsjHmnoN1LaOpSEj9B3RtJAUmbbocSS0rl+JFr9YbUkZvV+YzDr48Sqy4iVUY8kFAqxb98+ioqKMJlMPV0coQPa+5mJz3n3E+95S4osE6+vV6enb+iH0iju9eL77DN8//ynulYOgF6PPHEUkTPHo2SlEZMUPlC28G7Vp8SVODa9jetOuY5JmZPQSlqKHEW44nE1KDQ213SVxuG/bc1H4quEtqa+b2ROburM2vyWlAGaTtYZNIaS5v1KejiUiFWXBUEQhD4v7vWqAaVZP5RGkZISvB9/TGDtWpSo2ilV43Agnz2R4KQRYLOgSBKHrGH+WPUO+7z7ARiXPo4bTrkBh9GBRWdhoMGFqWYfRDo2BX6bFEWtCWlrPhJPGchH6DSrt7QcYZO4ZanPdZZG19RkY2ioKdFbQdfBPin9jAgrgiAIQpdrqx9KI0WWCX37Ld6VKwk1G02nLyxAnjoZ34h80GlRJIm408bK2Bb+uu9dYnIMs87M7OGzmZw1GUmSSEFLQTiE1run44WLhtQ+I22FkiOFHY1OnXvEnt06lJgcxzb8V6NrarJp3qfkJA0l7RFhRRAEQegSiX4o1dXIvtZ/9OVwGP/q1erQ4zK1vwmShHnCBDh3MvXZScgNzSmyy065Q2bJrlf5of4HAEamjGTOyDkkm5KRIj7yYgoZUjvNH3KsYfhvG6EkUHOE70ICa1pT35HmocSa1rHhv23RaJv1KWk2CkdnPLbjnWREWBEEQRCOWaIfSnW12tekja4bsZoafP/8J77PP090ppXMZmznnIPmnCnUWuNE4hFAQXEkEctM4bOqtbz177eIxCMYtUauGXoNZ+eejRT1o687wCBdErbDJ2yLhWHbX2DfF+AtB+UIE68Z7U0hxJHdrB9J1vEFCI22qcmm+SgcvejzdzxEWBEEQRA6Le71qgGltnU/lEbhPXvwrlxJYONGdWQOoEtPJ2nGDHSTJ1KteAhEAxAHxWFDzkylWvLz8rZn+b5WnZl2WPIwfjbyZ6RqzVC3D3tcZoAlC/3hNRwHN8BXL6m1KY20xmZzkBwWSoxJx/cGNA8lzZtxRCg5IURYEQRBEDpEDoXUCdtqqpFDbY+0UeJxAl9/jXflSiJ7mvqRGIcPJ2nGDPRjRlIbrqU+VK7un2RFzkpDsZj4ouQL3tj5BuF4GIPWwFWDr+LczNPQ+KsgXEqmMZlca0rLmWd9FWpIObhBfWxJhQk3QsbIhuG/xzk3kqQ5rE9JQzjRmXr9FPX9iQgrgiAIQruO1g+lkez34/v8c7yrVhGvbZgSXqfDevrpakgpyKcuVMchTzGyIqPYzMiZaZBkpS5Ux9JN/8e26m0ADHIO4ufDZpMRV6B2L1pJQ5E1G5e+2WRo8Shs/ytsXaHOJitpYcRlMPoaNUx0VmMoad7J1WARoaSXEGFFEARBaKF5PxTZ7UaR259DJFpWhveTT/B/+SVKRB31o7HbSZo2Ddu0aWgdDvxRH4fc+4jEIyhWM3JmKthtKIrCutK1/Om7PxGIBdBpdFwxYCYzXCPRBNwAmLUGBlmyMTWfRr70G9jwgtpZFiBzFJx2KzjzDy9ea5JGbapp1afELEJJLybCiiAIggA064dSV4cSbX92V0VRCO/YgWflSkJbtiS26/PySLrgAqynnYZkMBCOhynzHcQf8aOYjchZueBQ+4q4w25e3fEq31R+A0CRvYCfF15CtmRIzPqaYkiiwJyBtrEpx18NG/8I+1erj80uOPXnUHRO66AhSU01JYfP7CpCSZ8jwoogCMJJrCP9UBL7RiIE1q7F+/HHRA811GpIEuaxY0m64AKMw4YhSRJxJU5NoIK6UJ0aUopywNk0s+lX5V/x+o7X8UV9aCUtl+ZO5SLXSLRIoChISOSZ08gwOhtOHIMd78KWP6sz00oaGHYJjJ3dekVgkwOSB6jbRSjpN0RYEfqV1atX8+ijj7J+/XoURaGgoIDZs2ezYMEC/vnPf3L//fezc+dONBoNgwYN4sEHH+RHP/pRTxdbELqVEosRq6klXlNN3Os76v7x+nq8//wnvs8+Q/Z6AZCMRqxnn03S9OnoG1YnV1CoD9dRHawmpteiFGSjuOyJ0OCNeHn9u9fZWL4RgDxLJj/PPZ98U2riXHqNjkGWrKZhyeXfwvr/A3fDauNpw+H0W9VA0pykAVcBOHKP560ReikRVoTOURSIBrrvfHpLh/87ev/997n22mt5+OGHee2110hNTeX777/niSeeYN26dVx11VW8/vrr/PjHPyYajbJhw4YTsoqyIPRGiiwTd7uJN8yHcqR+KI0ixcV4Vq4ksGEDxNXhydrUVJKmT8d21llorE21Gv6oj8pAFWGtjJKbjpLcckbXbyq/4ZXtr+CJeNCg4eKMifw4bRK6ZkOQ7TpL07DkYB18/TLs/Ux90miHCTfBoPNaj/AxWCBtWOtaFqHfEGFF6JxoAB7L7r7zLSzt0AVIURRuv/12fvvb37JgwYLE9mHDhrFs2TL+8pe/kJGRwWWXXQaAVqvlnHPOOUGFFoTeI+7zEauqOmo/lEaKLBPctAnvxx8T3rkzsd04eDBJF1yAefx4JG1TwIjEI1QGK/ERQslORUlxtggp/qifP33/J9aVrgMg25jCz/NmUGTJaHHeTGMyuaYUJEWG7/4O37zW8I+RBEMuhPE3tD03iiMHnIUg/vHo10RYEfqF3bt3s2/fPq699to2n58wYQKlpaXceuutXHrppUyaNInk5ORuLqUgdI/O9ENJvCYYxPfFF3g/+YR4VZW6UavFMmkSSTNmYBzQstklrsSpCdZQF/cgp6egpOa0CgzfVn3L0u1LqQ/XIyFxYdp4LsuYjL7ZSsNaSUORJVMdllz5HWz4P6jdqz6ZMlht8kkd0rrAOqO63ezs8Psi9F0irAido7eotR3deb4OqGq4uObk5LT5fFFREWvWrOH3v/89N998M2VlZUybNo0XX3yRAYddhAWhL+psP5RGscpKvJ98gu+LL1BCIQA0Viu2adOwnXceOper5XlQcIfdVEfriKbaUVIHgrZlSAnGgrz53Rt8UfolABkGJz/Pm8Ega8ta2cSw5GgQvvof+OET9QmDDcbfCINntL0Wjy0dkgeCVvwJO1mIn7TQOZLUK9uFU1PVDnqHDh1i4MCBbe4zfvx4XnvtNQD27NnD3Llzue6661i7dm23lVMQupKiKMTr6zvVD6XxdeGdO/GuXEnwm2/UvmiAPjubpAsuwDJ5Mhpj6/VxAjE/laEagilWlLTCViEF4Lvqbby87WVqwvUAnJ8ylp9knYFR03LBwRRDEgWmNLS7P4FNr0CkIWANOl/tm2JytC64Rgcpg8CW1qHvU+g/RFgR+oUhQ4ZQWFjIG2+8waJFi466/8CBA/n1r3/NrFmzuqF0gtC14j4f8epqYrW1HeqH0kiJRvFv2KAOPd6/P7HdNHo0STNmYBo5suVU9g2icoTKUA0epx6lIBd0rWs7wtEgb+14nU/L1b4pqQY7P8+dwVBb69E5+eZ0MnzV8OnvoGa3utFVCKffBumntF14s1Nt9hGrFB8zRVGIyQrxhltMVpAbvypN2+PNtsXiCgadxKD041xL6TiJsCL0C5Ik8eyzz3Lttddit9uZNWsWKSkp7Nq1iyeffJIbbriB7du3c9lll5GdnU15eTl/+MMfmDJlSk8XXRA6RA6H1YBSU4McDHXqtXGPB9+nn+L99FNktzozrGQwYD3jDHUq/Oy2O83LSpyacB01NgU5Lx10bfzJkGV2V2xiyfd/ojJSD8DU5FFclXUW5uazztIwLFmbhO2bP8POfwCKOknb2OvUeVPaavKRNGqQcbTdxHsyaB4i4opCPK5+jckyskzLr4cFkua3Dla8tWIxtPFz6WYirAj9xiWXXMKHH37II488wr333gtAfn4+119/PS6Xi5UrV/Lwww/j8XhwOBxMnz6dJ598sodLLQjtU2IxYrW1ajNPJ/qhNIocPIj344/xr1sH0SgAWpcL2/nnY5s6Fa3N1u5r3VEPlZYY0RwX6Nv4U6EoRLwV/G3PX/m48msUwKW3cVPudEYmFbTaPUlrYlD59+g2vQIhNTBRdI46A62lnc7uBiukDe2VTc9HoyiH1VA0r8WQDwsUSutg0by2QznGkNGfiLAi9CtnnnkmH330UZvPvfvuu91cGkHovGPth5J4vSwT2roVz8qVhHfsSGw3DBig9kc59VSktmpIGgTjQSrMEQJZVjDoW++gKBCsY2/FFpbsf5+ycB0AZ7hO4drss7FoTa1ekhf0k7H5j0iVDeVx5Klr+WSNbrsQkgT27B4ZkiwfFjBa1VA0q9mIyzJxmcNCR9M2oeuIsCIIgtALxH1+4tVVne6H0kgOhfCvXo3344+JVVSoGzUaLKeeStIFF2AYOLDN/iiNYkqMKlOEOpcejG2MwlMUCNUT9ZTyXukXfFj5NTIKDp2FG3PPZ6y99ag6fSzCsB/+hWnXSlBkdQXjMdfC8JmgbSMIgbpP6uBjHpKsKAr+SJx4XG0maavWojFYxOKtaz1ELUbvJMKKIAhCD5IDAcI//NDpfiiNYjU16tDjf/0LJaDOLi1ZLNimTiXp/PPRpaQc+fzI1JllqhwgG9tpbgm6wVfOAe9B/njwY0pC1QCc5hzK7OypTVPjN1IUMsu+JXfbu0hBteaFgikw8RawHmEkz3EMSY7FZap8YcrdIUJRudOvF3o3EVYEQRB6kBKNdjqoKIpC5Icf8K5cSeDf/wZZ/eOsy8ggacYMrGeeicbUujnmcF6rRHkSRI3t1HIE3eCrIBb184/Kr/l7xQbiyNi0Zm7IncapjsGtXmLyVjBo27uYK79XNyRlwWm/hJwJ7RfkOIYkByNxyj0hqrxh0fTSRRRFIRyT8Ydj+MIx4rJCMBrn9AFHDr4nkggrgiAIfYQSixHYuBHvxx8T2bs3sd14yinYL7gA0+jRSB3o4xGxGSmzx/DrZaCNkR4hD/gqIBrkUKiGPx5cyf5gJQDj7QO5Ifc87LqWTUWaWJjsnZ+QufdfSHIctAYYdRWM/Il6vz3HOCS5zh+h3BOiPhDt1OtOFoqiEInL+MNx/OGYGjwisYb78UQQ8Ueanlf3UR8fHvwGpFr59K6pPfPNIMKKIAhCrxf3+fB9/jm+VauI1zU0q+j1WCdPJmnGDAx5eR06jmwzU+mUqNUG294h7FVDSiSArMisrNrE3yrWEVPiWLVGZmefy2nOoS37vigKzrJvKdj2LobGJp/ciTBpLiRltl+YYxiSHJcVqrxhyj0hgpF4h1/Xl0VishoqwjH8kaag4Y80hI1mYaTx+cb9Y8dZ06SVJKxGLTaTjvyUjs0mfqKIsCIIgtBLRUtL1aHHa9agRCIAaBwOks47D9u556K12zt0HCnJSp1LT7nGi6y00Z8j4gNvBUT8AJSH61hy8GP2BMoAGJ1UyI2556vr9zRj9FWR/+3fcDY2+dgyYNIvIO+0Ixeok0OSQ9E45e4QVb4wsXjfa+qJxuWmGo1IDF+oMVjE8DVsa6rdaAgbDduix/n9aiSwGnVYDTqsRi1Wow6bUdd6m0GH1aTDatAmnjfqNEiShMWgZUyes2vejGMkwoogCEIvoigKoW3b8K5cSejbbxPb9fn5JF1wAdbTTkPSt9PH5DAamxV/mpVDkptIPACH/92L+NWalLA6h4usKPyzZjNvl60hosQwaQxcm30OZ7pOaVGbIsUjZO3+lKzdn6KRYygaHdLIn6jNProj9JWRJLDngLOgQ0OS6wNNTT09PUonFpdbN5tEmjWnHLataXucSPz4OvxKElgNjSFD2xAyDgsfh21r3N+k1xxxFFhfIcKKIAhCLyCHwwTWrlWnwi9tWCxUkjCPG0fSBRdgHDq0w390NBYLscxkDmo9eCNVrXeIBBpCijexqSri5uWDH7PTfwiAU2z53JR7PimGlrU3jvId5H/7N0yBGnVD9jikSb88enNOB4ckn8imnpgsE2inGeXwIHL443DsOAMHYGkWNFoHD/Vx8+BhawgeJr0WTT8IHMdDhBVBEIQeFK2spP4vf8H32WfIPrWGQzKZsJ19NknTp6NLT+/wsTQWM1JmOmWGANXBQ3D43/poUG3uCXsSmxRF4fPab1lR9iVhOYpRo+fqrDOZmjy6RTgyBGrJ//YdXOXb1NdZkpEm/gIKzlD/9T8SWwYkDzjikORQNE6FJ0Sl99ibehRF4eviOr45WI8vHD0siMSOe0izhDr1fCJcGA5vVmn2nLFlk4rZIALH8RBhRehXfvazn7F06VJ27NjB8OHDAVi2bBnPPPMMmzdvbrHvnDlzcDqdTJw4kblz5wLqxS4QCGC1NrWlv/jii8yePZuKigruuecePvjgA7xeL0VFRdx6663Mmzcvse8DDzzAgw8+yKJFi3jkkUcS2zdu3MikSZMYM2ZMq3K0VWaAzz//nHPPPRer1YokSSQnJ3P99dfz0EMPoenmWT2FE6Purbcof+BBiKupQpuWRtL06djOPhuN2XyUVzfRmE1os7OoNkYpC5QgBw/7oxwNga8SQvUtNtdGvCwt+YTtvgMADLHm8LPcGaQbm1Y8luIxMvd8TtauT9DGoyiSFmX4TDRjrwX9UTpdavXqkGRraru7uANRyj0h6gKRY27qURSFbw+5eWPjQfZV+4+6v1mvbarFaKv/RjuPLQYtGo0IHD1BhBWh3/B6vaxYsYLk5GSWLFnC008/3aHXzZ49m9mzZwNQXFxMUVERJSUlOJ3OxD719fWcccYZjBs3jq+++ors7GxWr17NnDlz2L9/P0899VRi36FDh/Lqq6+2CBVLly5l2LBhnS6zw+Ggvr4egG+//Zbzzz+fgoICbrnlls68NUIvZR49BuJxjEOHkjRjBubx4zs09LiRxmhAn5ODx6alxFdC2B9uuUMsrDb3hNw0TwKKorCmbgd/Lv0XQTmCXtLyk8wzOD91XIv//u2Vu8j/9q+Yfeqw5UjaMAyT5yO5Wq/90/qbc6nNPm0MSY7LCtUNE7gFjrOpZ2e5lze/PsB3ZWqTlkmvYdqwDDLtxqa+G6am8GEx6NCKwNHniLAidIqiKARj7Qx7PAHMOnOH2+nffPNNrFYrjz76KIsWLeLxxx9H38GOiEfzzDPPoNPpeOONN9Bq1Xkpzj33XF5//XWmTp3KrbfeyoAB6nTjw4YNo6KiglWrVjFjxgxCoRBvvfUWd9xxBytWrDjmMo8aNYqzzjqLrVu3dsn3JPQ809AhFL61AjnQud8pyaBHn51DxGVlj68Ej9vTcodYBPyVEKzj8OqK+qifV0pWscW7D4ABlkx+njuDLFPTYoL6YD35294juXQzAFFjErEJN2EeNP3oTT6SBpKL1LV9DtPY1FPlDR/3KJfiGj8rNh7km4P1apm1EtNPyeTSMdnYzV3ze38y0WkldBoJvVbTcF+DXiuh02ow6Xq+JleEFaFTgrEgp/3pKMMSu9CGWRuwHK2qucGSJUuYPXs211xzDQsWLODvf/87V1xxRZeUY+XKlVx99dWJoNLorLPOIjs7m1WrVvGLX/wisf2mm27i5ZdfZsaMGfztb39j0qRJZGe3vnh3psxbtmzhiy++ECtF9zOG/HxC3+/s0L6SXo8+JxslxckhfxlVtftb7hCPgK8KgrWtQoqiKGyo38ny0s/wx8PoJC2XZpzOhWkT0ErqHyNJjpO+90tyvl+JNh5GQaJu4FSSTr0Zs8nBURmskDYMDC1/Z93BKBWeELX+Y2/qaVRWH+Stf5ewbq/awVcjwdSh6VwxLocUW+cmluuvJEkNb3qtpkXoaB5G9JqGUNJwv7c3b4mwIvQLO3bsYP369bzwwgvYbDYuv/xylixZ0mVhpbq6us2wAZCdnU1VVcsRFz/96U+5++67qaurY+nSpfziF7/A19B5sjNldrvdOJ1OJEkiPT2d+fPnM2fOnC75noS+Q9Lr0GdloU1PpzJURWnNduJKs+aTeFTtk9JGSAHwxAK8duhT/u3+AYB8Uxo3519ArqmpL4mteg8FW9/G4i0HwOcqxHfqTWRkjT967WbjkGRXYaLmRW5s6vGE8IePf1RPjS/M25sO8a9dlTTOdTZlYApXTsgly9Hx/j190eG1HnqtplXY0CXCiRpM+hsRVoROMevMbJi1oVvP1xFLlixhzJgxjBkzBoAbb7yRCy+8kEOHDqHX64lGW0/JHY1GO9xMlJqaSmnjcNLDlJaWkpbWck0Th8PBj370I5588kk2b97MzJkz+dOf/tThMufk5CSO09hnRTj5SDot+sxMdJmZuKNeDtbuIBRvto5QPKY29wRq1VWN2/Bv9w+8WvJPvPEgWjRckjGJi9MnopPUWkJdyEve9r+TWvI1AFGDldIRPyZp6I/JNHagNkVngrQh0FDzEo7FqXCHqfSGjrupB9RamXc3H+KTHRWJGVnH57u4+tRcClI6Nqlcb6KRQKdtqO1oo9aj8XHz5/vDPCnHS4QVoVMkSepws0x3iUajvPbaa/h8PjIz1em9FUUhHo+zbNkyzjnnHEpKSojFYuh0TR/5PXv2MGXKlA6dY/r06bz11lvcf//9LZqCVq9eTWlpKeedd16r19x0001Mnz6defPmYTC0XBvlaGVetGhRp98Hof+QtBp0mZnoMzMJEWWv+wc8kWb9UuJxCFSBv7rdkOKLhfhT6Wesr1ebmHJMKdycewEFloah0IpM+r415Hz3IbpYCAWJqsLTqRlxKYWuIZiPtJ5PI1sGpAwEjRZPKEq5u2uaegD84RgffFvGP74tS8xxckqWnZ9OzGNIRtLxn6CLNK/RaK+JJVEbotWIzr3HSIQVoc9777338Hg8bN68ucUInueff56XX36Zu+66i4yMDBYtWsR9992HXq9n+fLlbN++nUsvvbRD57jjjjt4/fXXmT17Nk899RSZmZmsXbuWG2+8kdtvv52BAwe2es20adP45JNPGDVqVKfLvHDhwk6/D0LfJ2k16DIy0GdmEtPCAd8hKgOVTTskQkoNKO03rWzx7GNZySrcMT8SEj9KO5WZGaeh16iXfGttMQVb38bqVieA8zty2T/mJxjTRzLEkpHow9KuhiHJsjmFan+YCrcPXzh23N8/qDUzK7eV897W0kTz0YBUKz+dmMeoHMcJrWU4vNbDoFO/Ht7EImo9up8IK0Kft2TJEq699tpWQ4Nvv/12nnrqKdasWcM//vEP7rrrLgYMGEA0GmXkyJF89NFH5ObmdugcLpeLNWvWsHDhQiZMmIDX66WwsJA777yT+fPnt/kaSZLarHHpSJk/++wzMZfKSUTSSOjS09FnZYFeT2WgklJ/KTG5IQDIcQjUqDUpcvuhIBAP80bpF6yu2w5AptHFz/NmMNCSBYA24idvx/uk7VebcmN6MyXDf0RV4WTyLBlkGl1HL6zZRcQ1iAq/TGV1HZFY18yDH4vLfLqzkr9tOkR9UG22zXGa+empeZxa6Op0KJAkEuHi8FoPva7tJhhR69F7SYrS0ysu9DyPx4PD4cDtdmPv4MJgJ4tQKMS+ffsoKirCZDrCmh9Cr9Hez0x8zrtfR95zuXGBQoMBd9jNQe/BpukBZLkhpFQdMaQAbPceYGnJJ9RGvUjA9NTxXJE5BYNGB4pM6v4N5O34AF00AEBV3kRKRlyCZHYx0JJF0tH6h0kafNY8ymQXNV3U1KN+iwqrf6jmL/8uocqnzhOTnmTkygm5nDEwtVOjVCQJMuwmsp0mDNr+sSZOX9Ad1xZRsyIIgtCDNAYDwViQg3W7cIfd6kZZVkf2+CqPGlJC8Qhvla/msxp1/p00g4Of581giFXtpG2pL6Fg61+w1amz1AbsWewf/RN8KQNI0pkZaMlKNA+1RVYU3LKJUn0+Ho8eiBz/N43aR2tjcR0rvj7IoXo1nDnNei4fn8O0oemdHtFiM+ooSrNiM4o/a/2R+KkKgiD0IH/Uz46aHeoDRVFH9vgr1eHIR7HTV8LLJZ9QFVFDzrkpo7kq80xMWgPaaJCc7/5B+r61SCjEdUYODbuQiqIzQaMl0+gi15Tabu1DNK5QF4hQoUnDZ8qBWNfUUjROjf/mxoPsbZga32rUMnNMDheMyMCo0x7lCC3ptBK5LjOZdpOoSenHRFgRBEHoQTE5poaUYJ1akxI/es1FRI7xdvkaVlV/gwKk6JO4KXc6pyTlg6KQcmAjeTv+jj6szu1TkzOOgyNmEjU70EgaiswZJBvaHlETjMap9Yepi+gI2AcSN3Rdtf6uCi9vbGyaGt+o03DxqCwuHp2FxdD5P0cpNgMFKZZOBxyh7xFhRRAEoSfFQlC9S13HpwP2+Mv448GPqYjUAXCWawTXZJ+NWWvE7CmjYOvbJNXsBSBoS2f/6J/gTRsMgFlrYKAlu9WwZAUFTzBGrT9CIBInYk4nlFwAmq4JAftr/Kz4+iCbDtQDDVPjD89g5tgcHMcwNb5Jr6Eo1YrT0oHh1UK/IMKKIAhCT4qFOxRUonKMdyrW81HVv1FQcOis3JR7PqPtRWiiIXK+e5eMvV8iKTJxrYHSodOpGHgOSkN/lGR9EoWHDUuOyWpTT60/QiyuoGh0BJ1DiDVbJ+h4lNUHeWtTCev2dM3U+BoJsp1mcpzmXj89vNC1RFgRBEHo5YoDFSw5+DGHwuof/cnOYVybPRWb1kjyoW/I2/YuhpA6aVxt1igOjryMiKVpGHKeOa3FsGS1qSeCOxhNjOqJGZ0E7QNQOjIZ3FG0NTX+5IEpXDU+lyznsU2NbzfrGJBqw2wQTT4nIxFWBEEQeqmYEuf9iq94v/IrZBTsOgvX50xjgmMQJm8l+d/+FUfVLgBC1hT2j7oCT8bwxOv1Gl1iWPLhTT2NFElDKKmAqCXjuMvraZwa/7uKxFT74/OdXHVqHoXHODW+XiuRn2IhPUlMnXAyE2FFEAShFyoJVvPHgys5EFIXyTzVMZjrcs7FiZasHR+Q+cPnaJQ4skZH2ZDzKBs0DUXb1P+jcViyhJYqXzjR1NNcXG8l6BiE3ME1uNoTiMT4YGsZ/9hWRiiqTo0/PCuJaybmH9fU+Ol2I/nJFvT9cGE+oXNEWBGEY7Bs2TKeeeYZNm/eDEBhYSHPPPMMl112WaePVV9fj8vlYt++fRQWFnZpOYW+J67IfFT1Ne9UrCeuyFi1Jq7LOZdJ9sG4KraT/+07GINq59r6jFM4MOpywtaUFsfINLpI0bqo9EZxBwKtJ3CTJMKWbMK23MQqycciHIuzcnsF72051KVT41sMWorSrNhNne98K/RPIqwI/UJbYaG4uJiioiLq6uparL8jCL1VWaiWPx5cyb5gBQBjkwZwY+55pEdC5H+1BGfFdwCEzS4OjLqc+swRLcKGJEmkaVKJBQzsCwfaPIesNRF0DCTeztDljmhvavyrT81j4jFMjd9Iq1HnTMlyiDlThJZEWBGEwyiKgizLLVZXFoQTSVZkPqnezF/L1xBV4pg1BmZlT+UM+yCyf/iMrN3/RCPHkCUt5YPOpWzI+ci6po6wcVkhFtXgUJLxyVqg7UUOI+Z0QknHPiRZlhXW7FGnxq/0qiOY0mzq1PhnDurc1PiHS7aqc6aY9OL3TmhNhBWhUxRFQQkGu+18ktncJf9hTZ06lcmTJ7Np0ybWrl3L4MGDeeWVVxIrIhcWFjJ37lzeffddtmzZwldffYXZbOZXv/oVGzduxOVyMX/+fBYsWNCh861atYqFCxeya9cucnJyePzxx5k5cyYA4XCYBQsW8Oabb+JwOFi0aNFxf39C33XIX8aTe/7C7kApACNs+dyUO52iuhLyP38Kk18dAeROG8KBUVcQSkpPvDYSl/GEouhiZjJ1qWjaWS1Z0egI2gcc85BkRVH4uriON7toavzmjHoNhSlWkq1izhShfSKsCJ2iBIPsHD+h2843dNO/kSyWLjnWa6+9xgcffMCIESO47bbbmD9/Pp9//nni+WXLlvHee+8xaNAg4vE4o0ePZubMmbz77rvs2rWLCy+8kPT0dGbNmnXE82zdupWrrrqKt99+m6lTp7J27VouvvhivvrqK4YOHcqjjz7KunXr2LZtGxaL5ajHE/qvj/Z9xL1rFhOKhzFq9Pw062ymm3Io2PwXXGXfAhAx2Tkw8jLqsseAJKGgEIjIeENRQlGZDH0KKXpHu+c4niHJXT01fnOSBFkOE7kui1jtWDgqEVaEk8Z1113HmDFjALjxxhu58MILWzx/6623MnToUADWr19PWVkZjzzyCAaDgdGjRzNv3jyWLVt21HDx4osvMmfOHKZNmwbAmWeeySWXXMKKFSu49957Wb58OY8//jjZ2dkA3H///Xz44Ydd/e0KfUC2LZtIPMoway4/y5nGyAPfkL1rOdp4BEXSUDHgbA4NnYGsNxGXFXyhKN5QlGhcQSfpKDRkY9G2PaT3eIck76rw8ubGg+woU+dvMeo0/GhUFhePysJ6nIsFJpl0DEizHtMU+8LJSXxShE6RzGaGbvp3t56vI/R6PdFoy4XfGh/r9eqIgszMzMRzVqsVn8/XYv/8/PzE/ZKSErKzszEYmv4bHTBgAK+//vpRy1JcXMynn37K0qVLE9tisVhi6fTS0lIKCgoSzzW/L5xcRqeN5vmzf4dr16cUrXkJs68SAG/KAPaP/glBexaRuIzXH8EXjiE3zLBm1ZrJNaSjk9q+hMf1toYhyZ2fm+TwqfF1Gonpp2Rw6TFOjd+cXiuRn2wh3S7mTBE6R4QVoVMkSeqyZpmuVFBQwL59+1ps27NnD6mpqVitHZuMSqNpanfPzc2ltLSUaDSaCDvFxcXk5uYe9Th5eXn8+te/5oknnmjz+ezsbPbv389pp50GwIEDBzpUPqEf8pQxcf1SDDvVmrWo0cbBETOpyZ1AICrj8YQIRlp2lk3RO8nQJbfdl0uSCFtzCFtzOj0kucwd5C//LmFts6nxzxmSzhXjc0g9hqnxD5eWZKQgRcyZIhwbEVaEfuG6667j/vvvZ/r06YwdO5YDBw7wwAMPMHv27GM63qRJk8jIyOC+++7jgQceYPfu3Tz77LP87ne/O+pr586dy4UXXsgFF1zA2WefTSwWY9OmTTidToYPH861117LE088wZlnnonFYuGhhx46pjIK/cD2v2LY+SEKEpVFZ3Bg2IW4ZT3e+mBiBthGGklDjiENu9bW5qGOdUhyjS/MX785xOc7m02NPyCFKyfkkn2MU+M3ZzZoKUq1HnetjHByE2FF6BduvPFGvF4vs2fPpqSkhNTUVK688koeeOCBYzqeXq/n/fffZ968eWRmZuJyubjzzjs71Bl23Lhx/PnPf2bx4sV89913aDQaxo4dy9NPPw3A4sWLqaysZOTIkdjtdhYvXswHH3xwTOUU+rhJvyBycD3fpQygzJiJzxdDliOtdjNpDOQaMjBq2u4kG7GkE7J1bkhyW1Pjj8tzcvXEY58avzmNBDkuM9kOseigcPwkRWk1t+FJx+Px4HA4cLvdiX4FgioUCrFv3z6KioowmUQ7c1/Q3s9MfM67X0fe85LKvazd8n67x3DobGTr09oclixr9IQcA4g1W6TwaNqaGn9Ypjo1/tDMY58orjmnRU9RqlXMmXKS6I5ri6hZEQRB6EGK3PZ2SZLUYcm6tocld3ZIcjgW5+PtFby3pRRfOAZAUaqVa45zavzmDDoNhSkWUrqgj4sgNCfCiiAIQi+jk3TkGTLaHJasDkkuJGpJb+OVrcXiMp/trOSv3xyiPqCOkMt2mrj61DwmFbbTUbeTJAky7SbyksWcKcKJIcKKIAhCL3KkYcmdGZLc1tT4qTYDV07I46zjnBq/OZtRR1GaFdtxzr0iCEciPl2CIAi9RLvDkjsxJLmtqfEdZj1XjMvh3GHpXTZ0WKeVyHNZyLAbxaKDwgknworQIaIfdt8hflZ9z5GGJcs6M0H7gKMOSVYUhW2lHt7ceIA9Vc2mxh+dzYwRmV3a2TXVZqAgxYpBJ+ZMEbqHCCvCEen1eiRJoqqqirS0NPEfVC+nKApVVVVIkpSYzE7o3QwaPQOMOW0OS+7okOQTOTV+cya9hgGpNhwW8dkSulenPsUTJ05k48aNice33HILf/jDHxKPs7OzKS0t7brSCT1Oq9WSm5tLSUkJxcXFPV0coQMkSSI3NxetVgwb7QtMOlOroNLRIcnq1PglbDpQB6hT459/SgaXjsnGaem6VYw1EmQ7zeQ4xZwpQs/oVFj5/vvvWzx+++23W4QVr9fbNaUSehWbzcbgwYNbrb0j9E56vV4ElT4sanQRcgxA0bRfe9E4Nf66PTUoqN1Ypg5J44rxuV0yNX5zDrM6Z4rZID5TQs/pVFg5vAng8LZx0UTQf2m1WvEHUBBOoI4MSa7xhfnbN4f4rNnU+KcPSOaqCXldMjV+cwadRH6ylbQkMWeK0POOqzFThBNBEITjd7QhyZ5glHe3lPLJjvLE1Phj85xcfWoeRanHPzV+c5IE6UlG8pMt6MSig0Iv0amwEovF+Pvf/56oUYlGoy0ex+PxI71cEARBOJzBhj95RJtDkgORGB98W8Y/vm05Nf5PJ+YxLLPrpzW3GtVFB5NMogOt0Lt0Kqykp6dz++23Jx6npqa2eJye3rEZFQVBEIQGGl2roNLe1Pg/PTWP0bldMzV+c1qNRK7LTJbDJGrMhV6pU2FFjAYRBKEnbN269aj7jB49uhtKcmJ1x9T4h0uxGShIsWDUiT5pQu/VJQPw6+rqcLk6vuqnIAhCZ4wdOzbxh7qtSe8kSerTzdDtT42fy5mD0k7IejtGvYaiFCsua9cNcRaEE6VTYeXVV18lPT2dCy+8EIBvvvmGmTNnUlpaypAhQ/j73//OoEGDTkhBBUE4eZ1yyinIsszPfvYzrr/+ejIyMnq6SF1CURQ+/b6S//pkJyV1TVPjXz4uh2ldODV+c5IE2Q4zOS6zWHRQ6DM69Zvw3//936SlpSUe33bbbZxzzjls2bKFs88+m7vvvrvLCygIgrBt2zZeeeUV9u7dy6hRo5g5cybvvPMOsVisp4t2XF78Yi93rNhMSV0Qq0HLNRPzeOanY7lgROYJCSpJJh2jcx3kp4jVkYW+RVI6sZBIcnIyVVVVaLVaKioqyMnJoaSkhMzMTOrr6xk+fDhlZWUnsrwnhMfjweFw4Ha7sdu7voe9IPQG/eVzHg6H+etf/8rSpUvZunUrn3zyCaNGjerpYrXpaO95pTfEBb//gnOGpHPJ6K6dGr85vVYiP9lCuv3oqzULQmd1x7XlmH8zNmzYQFFREZmZmQA4HA78fn+XFUwQBKEt4XAYt9uN2+3GbDaj0/XdJc7Sk0ysXHB2YuHBEyEtyUhBiuWE1NQIQnfp1Kd37NixvPzyywAsX76c6dOnJ547cOCA6GQrCMIJ889//pPZs2dTVFTEmjVreOyxx9i7dy/Dhw/v6aIdlxM1Csdi0DIix86gdJsIKkKf16l/SR5//HEuuugi5s+fj8PhYN26dYnn3nzzTaZMmdLlBRQEQSgsLMRisXDTTTfxxBNP4HA4gJbrkfXlpq2upNVI5LjMZIs5U4R+pFN9VgB8Ph87d+5kyJAhJCUlJbbv3LmTpKQksrOzu7yQJ1p/acsXhCPpy59zjaapZqCtNcp669Dljrzn7kCUHWWeLjmfy6qnMMWKSS/mTBG6T3dcWzpdN2iz2ZgwYUKLoAIwdOjQTgWVaDTKvHnzcLlcJCcnM3/+/HZ79h9t3zlz5mAwGLDZbIlb81ofQRD6tn379iVue/fubXFr3HYyM+g0DMmwMSzTLoKK0C91qhnoiiuuOOo+f/3rXzt0rEceeYTVq1ezY8cOAC666CIee+wx7rvvvmPa97bbbuOZZ57p0LkFQehbCgoKuvR40WiUO+64g+XLlyNJErNnz+b3v/99m5119+zZw7x581i/fj0Wi4Vf//rX/Od//meXludYSRJk2k3kJYuhyEL/1qmw8s477zB8+HAuu+wyDIbjm/Xw5Zdf5ve//z1ZWVkALFq0iLvuuqvNsNKZfTsiHA4TDocTjz2erqmCFQThxPjf//3fo+7TfJ2yo+noP0vxeJyZM2dy2WWX8d5777F3716mT59Obm4us2bN6tw30cWSTDqKUq0nbLizIPQmneqz8tZbb7Fs2TI2btzIT3/6U2666SbGjx/f6ZPW1dWRnJzM7t27EzPe7t69myFDhlBfX5/oPNfRfefMmcN7770HQFZWFj/72c+44447WrRzN/fAAw/w4IMPttreF9vyBaGj+nKflXPPPfeIz0uSxKefftrh4+Xl5fH73/+eK6+8ElCvbXfddRf79+9vsd+OHTsYPXo0gUAg8Q/agw8+yGeffcbnn3/e6rht/SOUl5fXpX1WdI1zpiQZRQdaoVfodfOsXHXVVVx11VWUlZXxyiuvMGvWLIxGIz/72c+4+eabsVqtHTqOz+cDwOl0JrY13vd6vS3CSkf2vf3223nqqadITk5m48aNXH311Wg0Gu644442z3/PPfdw5513Jh43XlAEQeidPvvssy47Vl1dHSUlJYwdOzaxbezYsRw4cAC3293i+iPLMtByPSJZlttdWPHxxx9v8x+hrpKWZCA/2YpBJ4YiCyeXY/rEZ2Vlcffdd7N9+3YuueQSfvOb37Bhw4YOv95mswFqTUajxvuHd9ztyL7jx48nLS0NrVbL6aefzt13382bb77Z7vmNRiN2u73FTRCEvqW+vp5vvvmGYDDYqdcd7R+g5oYOHUphYSH33Xcf4XCY7du38/LLL7fbdHzPPfckJqxzu90cPHiwU2Vrj9mg5ZQsO4PSk0RQEU5Kx/Sp/+GHH7jnnnsoLCzkyy+/5KWXXurUHCsul4vc3Fw2b96c2LZ582by8vJa/FfT2X0btdf8IwhC3/TUU0/x9ttvJx6vWrWKvLw8JkyYQH5+Pl9//XWHj9WZf5b0ej3vvvsu33zzDTk5OcyePZubbrqJlJSUNo/d1f8IaSTIdZkZnePAYdEf17EEoS/r1F/1l19+mbPOOotp06ah1Wr59NNP+eKLL/jZz36GydS5NSduuukmHn30UcrLyykvL+exxx7j5ptvPqZ9V6xYgcfjQVEUvv76a5544gl+8pOfdKo8giD0XkuXLuWUU05JPF6wYAFz587F4/Ewb948Fi9e3OFjdfYfoBEjRvDxxx9TXV3N5s2bCYfDnHPOOcf1/XSEw6xnTJ6TvGQLGjHSRzjJdaqDrUajYdiwYZx//vntrsfx3//93x06VjQaZcGCBfzpT38C4LrrrksMHfzlL38JwAsvvHDUfQHOPvtstm7dSiwWIycnh5///OfcddddHa5h6csdDwWho/ry59zpdFJXV4ckSRw4cICBAwdSVVWF0+kkGAxSUFBAZWVlh49333338f777/OPf/wDgB/96EdcdtllbY4w3Lp1KwMHDkSv1/P+++8zd+5c/vnPfzJ69OijnudYJoUz6CQKUqyk2owd/n4EoSf1ug62N9xwA5IktWrXPRZ6vZ7nnnuO5557rtVzjSGlI/sCfPHFF8ddHkEQei+9Xk8kEsFoNLJhwwaGDRuW6GdiNBoJhUKdOt69995LTU1NYl2h6667joULFwK0+mdpxYoV/N///R+hUIgxY8bwzjvvdCiodJYkQYbdRJ7LjE6s5SMILXR6uv3+qC//xykIHdWXP+czZ85k3Lhx/OIXv2DOnDmMHTuWp556ClCX+rj44ov54YcferiUrXW0ZuVAbYCiNCs2MWeK0Af1yun22xMKhXj66ae76nCCIAgJTz/9NH/605/Iy8vj0KFD/Md//Efiuddff52zzz67B0t3fKxGLSNz7CKoCMIRdPq34/PPP+ff//43Q4YM4cc//jHxeJznn3+eRx99lJSUFO66664TUU5BEE5iQ4YMYffu3dTU1CRG4tTX17Nv3z5uu+22dkcG9gWiyUcQjq5TYeWZZ55h0aJFjBgxgu3bt/Ob3/yGTz/9lEAgwHPPPdehtYMEQRA666mnnmLAgAGJUX6rVq3i8ssvx+/3k5KSwkcffcSECRN6uJSCIJwonYr0zz//PKtWreKrr77iww8/5NFHH+X8889n06ZN/OQnPxFTPwuCcEIcPnT517/+dYuhy4sWLerB0gmCcKJ1qoNtYweaRiaTCZ/P1+4w5r6iL3c8FISO6suf864eutxd+vJ7Lggd1es62B6eaywWS58PKoIg9H6NQ5eBLhm6LAhC39KppBEIBFqssuzxeFqturxp06auKZkgCEKDyZMn89hjj/GLX/yCl156iQsvvDDx3O7du0lPT+/B0gmCcKJ1KqwsWbLkRJVDEAShXU8//TQXX3wxDz/8MMOGDWP58uWJ5/r60GVBEI6uU2Hl3HPPPVHlEARBaFdbQ5cb/eY3v8FgMPRQyQRB6A6dCiuFhYUtRvw09mFp3CZJErFYrAuLJwiC0KSt1Y4b+670VbIsJ/rjCH2HXq9Hq9X2dDFOGp0KK3V1da22xeNxXnnlFR577DEGDRrUZQUTBEHo7yKRCPv27UOW5Z4uinAMnE4nmZmZYtqObtCpsHL4LJHvvfceCxcuRJZlXnzxRTEpnCAIQgcpikJZWRlarZa8vLwOrxIv9DxFUQgEAonh8llZWT1cov7vmMYdr169mt/+9rfs37+f+++/n5///OfiF00QBKETYrEYgUCA7OxsLBZLTxdH6CSz2QxAZWUl6enpoknoBOtUwti2bRs//vGPmTlzJjNnzuSHH37glltuEUFFEAShk+LxOIDoHNyHNYbMaDTawyXp/zpVszJ27FiSk5O5/fbbMZvNvPTSS632uf3227uscIIgCP2d6O/Qd4mfXffpVFg588wzkSSJf/3rX20+L0mSCCuCIAiCIHSpTrXffP7553z22Wft3j799NMTVU5BEAThJPWHP/yBrKwsbDYb33zzDSNGjOD9998HYNmyZYwdOzaxb2FhIe+8884xnae+vh5JkiguLj7+QgtdSnQ2EQRBENq0fPlybDZbq5tWq+W8887rljJEo1Fuv/12VqxYgc/nY9y4cWzfvp1LLrmkW84v9A4irAiCIAhtmj17Nj6fr8Xtww8/xGg0snjx4k4f71g6opaXlxMKhRg1alSnXyv0HyKsCIIg9AKKohCIxLrt1jgDeWccOnSIq666isceeyyx/MqqVauYNGkSTqeTESNG8N577yX2nzNnDj//+c+5+uqrsdvtvPDCC3i9Xn7xi1+QlZVFVlYWv/zlL/H7/W2e75tvvmHYsGEA5ObmMnDgQKBzTT1HKl84HObWW28lOTmZoqIi/vKXv3T6PRG6xzHNsyIIgiB0rWA0zin3rey28+146AIsho7/CYhEIlx55ZVMnz6dBQsWALB161auuuoq3n77baZOncratWu5+OKL+eqrrxg6dCgAf/7zn/nb3/7GG2+8QSgUYt68eRQXF7Nt2zYUReHKK6/kjjvuaHN0aWOTT1FRESUlJZ1eWuFo5Xv00UdZt24d27Ztw2KxMGvWrE4dX+g+omZFEARBOKp58+YRDodbhIoXX3yROXPmMG3aNDQaDWeeeSaXXHIJK1asSOwzY8YMLrjgAjQaDSaTieXLl/P444+TkpJCamoqjz32GK+++uoJWXLgaOVbvnw5CxcuJDs7G6fTyf3339/lZRC6hqhZEQRB6AXMei07HrqgW8/XUS+99BJ//etf+frrrxMztwIUFxfz6aefsnTp0sS2WCyG3W5PPM7Pz0/cr6qqIhKJUFhYmNg2YMAAwuEw1dXVfPLJJ8ydOxeAgoICtm/ffizfWofLV1paSkFBQeK55veF3kWEFUEQhF5AkqRONct0l/Xr13PHHXfw3nvvtQgZAHl5efz617/miSeeaPf1zWc4T0tLw2AwUFxcTEZGBqAGCqPRSGpqKrNnz2b27NldVvajlS87O5v9+/dz2mmnAXDgwIEuO7fQtUQzkCAIgtCm8vJyfvKTn/DQQw+1OVR57ty5LF26lM8++4x4PE44HGbdunV89913bR5Po9Ewa9YsFi1aRG1tLTU1NSxcuJDrr7/+hCzbcrTyXXvttTzxxBOUlpZSX1/PQw891OVlELqGCCuCIAhCm1566SVKS0u5//7725xvZdy4cfz5z39m8eLFpKWlkZOTw7333ks4HG73mP/zP/9DYWEhp5xyCiNGjGDQoEH893//9wkp/9HKt3jxYk499VRGjhzJ2LFjueyyy05IOYTjJynHMn6tn/F4PDgcDtxud4u2VkHoT8TnvPsd6T0PhULs27ePoqIiTCZTD5VQOB7iZ6jqjmuLqFkRBEEQBKFXE2FFEARBEIReTYQVQRAEQRB6NRFWBEEQBEHo1URY6aBwLI47ECUa7/pZFgVBEARBaF/vm4GolwpFZHaUeQAw6jVYDTqsRi02ow6LQYdBJ3KfIAiCIJwIIqwcg3BUJhyNUNtsoVCDToPVqMVq0KkBxqjFqOv4dNaCIAiCILRNhJUuEonJRGIydf5oYptBp06fbTPqsBp1WAxaTJ1Yj0MQBEEQBBFWOi4WRhMLIuvMR9+3QSSmEIlFqQ80BRi9tnmA0WI16kSAEQRBEIQjEB0tOioWxFa9BWvNNvSBSpDjx3SYaFzBHYxyqD7Irgof3xyoZ2NxLTtKPeyv8VPtCxOMHNuxBUEQToTVq1dz0UUX4XK5cDqdjBkzht/97ndEIhE+/PBDJk2ahMPhwOVyMXHiRP7xj3/0dJGFfkbUrHRApTfE6+tLSQlBkcNHcsSHyVtM1JRC1JxG3HB80wvHGgKMO9hUA6PVSIk+MNaGWhizXoskScf77QiCIHTY+++/z7XXXsvDDz/Ma6+9RmpqKt9//z1PPPEE69at46qrruL111/nxz/+MdFolA0bNpyQRQmFk5sIKx2w5aCb//2iJPHYboAiu0yRvYoiexWFTgMZ6elgS0PRGrvknHFZwROM4QnGEtu0GgmLQZsIL1aD2g9GBBhB6AcUBaKB7juf3gJHuXYoisLtt9/Ob3/7WxYsWJDYPmzYMJYtW8Zf/vIXMjIyEgsAarVazjnnnBNYaOFkJcJKB7gsei4Ylsz2klpKfeCJwJZq9aaKoJVKyLOVUODSk5fqIC89hfwUK06LocvKEZcVvKEY3lBTgNFIJDrvqqOQdFhFgBGEvicagMeyu+98C0vBYD3iLrt372bfvn1ce+21bT4/YcIESktLufXWW7n00kuZNGkSycnJJ6K0wklOhJUOOLUwmcyLMtm2/TusUhI1AQt7PbDPA8Ue2OsBfxSKvVDsjcKBakBNMg6TjvwUKwUpFvKTLRSkWMl2mtB1UTWprJAIMBWoy55rJLAY1OHTiZFIei0ajQgwgiB0XFVVFQA5OTltPl9UVMSaNWv4/e9/z80330xZWRnTpk3jxRdfZMCAAd1ZVKGfE2GlI/w1pL43lzRjOodShiOZUjjVZuf8PBtaSYuiQFVQDS/Nb6V+cIdifHvIzbeH3InDaTUSuU4z+SkWCpIbgkyKBbtJ3yXFlRXwhWP4wjEqGwKMJIHFoG0xEsli0KEVAUYQege9Ra3t6M7zHUVqaioAhw4dYuDAgW3uM378eF577TUA9uzZw9y5c7nuuutYu3Zt15VVOOmJsNIROz/AdGgDE4Dxe9+nyl5AScoIDqSOwGDJxqWzk24xkW6B0zKbXhaKNdS2NISXvR4o9kgEYgr7awPsrw3wJYm2JJwWPQUNtS/5yWpNTLbT3CWBQlHAH47jD8ep8jYFGLNemxhCbTXqsIoAIwg9Q5KO2izT3YYMGUJhYSFvvPEGixYtOur+AwcO5Ne//jWzZs3qhtIJJxMRVjpi4DTqT7+D+La3SPGVkO4pJt1TzPh9H1CdlEdJygjK0sZhshXg0NrQSGoTj0kHw1zqrZGiKFQEYa9Pxx6/mb1eDfvrQlR4wtQHotQH3GwpaaqF0Wslcl2NTUgNX5Ot2EzH/6NTFAhE4gQicaq8EUC9Xpr0WqyG5gFGi04revcLwslGkiSeffZZrr32Wux2O7NmzSIlJYVdu3bx5JNPcsMNN7B9+3Yuu+wysrOzKS8v5w9/+ANTpkzp6aIL/YwIKx0g27P5JHsoUc3N2CNecmp2kFu9nTTvAVK9B0n1HoTij6ix5VCaOor69ImYkgZi0rTuXCtJkGmBTEuMKXgBiBnseHTZ7A2Y2F8b4kCtn/01AQ7UBgjHZPZV+9lX7W9xnGSroaEWRg0w+SlWsuym4+6XoigQjMQJRuJU+yKJ7Sa9JhFebA39YfQiwAhCv3fJJZfw4Ycf8sgjj3DvvfcCkJ+fz/XXX4/L5WLlypU8/PDDeDweHA4H06dP58knn+zhUgv9jaQoitLThehpjb9kbrcbu731nClfl3/NTStvwiQZOMU8gFHmweQZMrBEvOTUNgQXz34kmt7KOmsWFWljCWRMwZA0oEOjcxSNlqixce6WJGRFocobZn9NgP21fg40BJjKhmacwxm0GvKSzeQnW1vUxFiNJyaTigUd+5ajfc6Frnek9zwUCrFv3z6KioowmUw9VELheIifoao7ri2iZqUD6sP1pJtTqQxWsynwPZsC35OitzPeNoQxOSM4kDcFTdBNTvV35NZsJ81djMtfhstfBsUf4rZkUps+gVDGGcj2wnbPI8lxDMFKDMFK4jozUXMamdZUMuzJTCpqGg4YiMQ4UBvgQI3a7+VAbYCDDbUwe6r87KlqWQuTajOQ39CRtyBZ7cybYTehOc7hze0t6GhrMZRaLOgoCIIgHB9Rs0LHUmFt/T7e//YV1tZ9x9fu3YTlptlmB1gyOd05jPFJgzFghKCX5PJtZFZuI61+DxpFTuzrM6dTnzGRUMYZhJMKjjopE5JE1OAkak4jZnS1u78sK1R4QonwojYj+Vs05TRn1GnIS24KLwXJVvKSzVgMXZ9fDTop0XlXLOjYc0TNSvcTNSv9m/gZqrrj2iLCCh17o93uA+za+wkAYTnKN+49rKv/nu3e/cgNzT9aNIyyFzLFOZwx9iL0Gh1S2I+9bBvJpVtJrt6FRmla9ydgTsObcRq+jNMJJRUdNbjIGj1RcypRc3qHF1T0hZtqYRr7whysCxCNt/1jT08yJpqQCpKt5KdYSEsyHnctzOH0WgmTXotGkpAkmn1VO/VJ0OI5AI2mabtGAhpf12xfqeG5xFcat7fc92Scc0aEle4nwkr/Jn6GKhFWukmHwoq/ml2H1kLEB9Fg0/aonw31O1lX/z37g5WJ7WaNgYnOIUxxDWeQJRuNJKGNBnGU7yC5dAuOyu/RyE0z0UbMqfgyT6c+bRJB+8CjBpe43kbEnE7UlAKaztVSxGWFck+IAzV+dQh1Q1+YWn/btTBmvTbRF6axKSkv2dLna0cOD0mgBpxWIanZY01D8DlSoGoMSW0FKo0kNWxrHahaBLeGY3UlEVa6nwgr/Zv4GapEWOkmHQorYTe76napD+IxNbSEverXuNokdChUw7q671hfv5PaqDfx2lS9ndNdw5jiGk6mUR3HrImGcFZ+h6t0K46K79DGm4KCbEklknM6/szT8DkGE44rRGNt/5gUSUPMlEzEnH7cCyp6Q9FmTUjqraSdWhgJyLCbGpqQmoZWp9qMYqr/LtThWidJwqTXMCDN1u6xRFjpfiKs9G/iZ6gSYaWbdOSN9kf9fF/7PXKz/icJ0SCEfWpwifiR5Tg7/SWsq/uOr90/EJKbgkiROYMpruFMdA7BrlNnkNTEIjgqv8dVugVn+Q608WajfczJkD+ZeMEZhF1DCcuS2rE1LhOKxok1CxKy1kTEnEbUnNplCyrGZJmy+lBDiPGrX2sD1Aeibe5vMWgTE9o17wsjOtmeeBaDljF5znafF2Gl+4mw0r+Jn6FKhJVu0tE3WlZkPGEPdeE66sP1xJo14zTtJEPUr4aXsJdwxMtmz17W1X3HtsP6t4xMKmCyazjj7APQa9SOrVI8iqNyJ67SrTjLt6GLhZqObXJC/ulQcCZkjgKNlrisEIrFicRkwlGZUEwmHI8T0tqJmNKImZJB6vrhxJ5gVO3M22xYdUl9kLjcRi2MBFkNtTD5ydbE/DDJVoOohelCIqz0PiKs9G/iZ6gSYaWbdOSNVmQZJCnxx1VRFHxRH3WhOurCdUTibff3IB5taDLy4faV81XdDtbWfdeqf8upjsFMcQ1nsDUnUbUvxWPYq3fhKt2Kq2wbuubLxxvtDcHlDMgcDdqW6wrFZIVwLE5Y1uLXu/DpUwhgItJOc1JXiMVlDtUHE01I+xuGVnuCbdfCGLQaDDoNOo2ETiuh02jQaqTWjxvuq1+lhn2aHqv7a1rcP/w4Oo2Ettn9pmO1PH/TMZvO3dWdi08UEVZ6HxFW+jfxM1SJsNJNOvJGx91uwrt3o0myo3XY0drtaCxNC4EFogHqwnXUheoIxoJtHgNINBmVuvexrvLfrK/7nppm/VtS9Elq/xbncLJMTXOrSHKcpOofcJVuIblsG7qIr+mYBivkNQSX7HGtgkvz/aKWdAL6ZIJxDYFIjGBUna22vdFBXaE+EGkRXg7U+CmtDxHvIx89jUTLINMsGGmbPda2Ck3NwlDz55oFInWflo8TYemwMKZtJ5g1nttu0jG+wNVubZUIK91PhJX+TfwMVSKsdJOOhpXQ9ztbbJMMerQOB1q7Gl4kgzq9figWoj5cT12oDl/U19bhVLKMHPGyq+pb1lV8zdd13xNs1r+l0JzBZNcwTnMOTfRvUV8XJ6lmL6ll3+Iq+xZtqGktIfRmyDsN8qdAzgTQtdF3RdKAJRlsGWBW526JxGR1mv1onEAkRqDhfuwEhZhoXKbWHyEmK8TiMjFZIS4rLR/HFWKy3LCt6X78sMfq/abXqfvITdsbnmt1/MMeNx6zjZasPiPXZWb1b6e1+ZwIK91PhJX+TfwMVWIG215OiUSJVVUTq1JXTtZYzGgdDvR2OxlJaWRaM4nGo4kaF0/E0/IAGg0ak4NheWcyLO9MZkf8bC7/inVlG9hW/wPFwQqKgxW8WfoFI5MKmeIaxlj7QAwaHd60wXjTBrNv1GXYaovJLN+B49BmNMFa2Pu5etOZIHeiWuOScyroG36ZFBn81epNZwRbOgZbBgaLGQcta2XCsTihiEwg2hBguijE6LUaMuy985dbbgwxcvOA0xCO4u2EIVluCFct92t8XatgdHgAO+y18cOP02J7y+eaE+s1CSfKz372M5YuXcqOHTsYPnw4AMuWLeOZZ55h8+bNLfadM2cOTqeTiRMnMnfuXEBtOg8EAlitTStLv/jii8yePZuKigruuecePvjgA7xeL0VFRdx6663Mmzcvse8DDzzAgw8+yKJFi3jkkUcS2zdu3MikSZMYM2ZMq3K0VWaAzz//nHPPPRer1YokSSQnJ3P99dfz0EMPodGI36HeSISVDgjv20fJr+ZhKCjAOHQoxiFD0KWktNpPDgSRA0GiZeVIGgmNzYbWbifF4STNlUZcieMOu6kL1eGOuFuNLDIYrEzKP5dJ+efiiXj4qmQN68rWsc9XwlbvPrZ692HWGJjgGMxk1zCGWnPRSBp8KQP4IWUAnPIjXPWlZJXvwHLo30j+Kij+Ur1pjWpNS8EZaoAxNNTUxMJQf1C9mRxqbYs1NTF3i1GnTpffVogJNqzYHIjECUXVr211sO1rNBoJg0bCQO+/aCmKQlxRg49Bq+GUbFFj0lcpinLkJuQuZtaZO9zB3ev1smLFCpKTk1myZAlPP/10h143e/ZsZs+eDUBxcTFFRUWUlJTgdDoT+9TX13PGGWcwbtw4vvrqK7Kzs1m9ejVz5sxh//79PPXUU4l9hw4dyquvvtoiVCxdupRhw4Z1uswOh4P6+noAvv32W84//3wKCgq45ZZbOvS9Cd1LhJUOCHz9NZG9e4ns3Yvvs88A0KamYhwyBOOQIZiGDEGXnd3iF1+RFeIeL3GPF0oOIel1aJOSsNsduBy54BxwxJFFdoOd8wdcxPkDLqLMV8a60rWsK11LTbiO1XXbWV23nWR9Eqc7hzLZNZwcUwpIGupcudS5cpGGzSDLX0Na2TYMBzaArxwOrFVvGj3kjFeDS94kMDTMzRFyq7favWpgsWWAqe0/fo0hxmlpuT3U0Acm0PC1sSamP4SY3kiSJHSShE6jdrB1WVuv9C30DcFYkNP+dFq3nW/DrA1Y9Jaj7wi8+eabWK1WHn30URYtWsTjjz+OXt9O37hOeuaZZ9DpdLzxxhtoteo/Seeeey6vv/46U6dO5dZbb2XAgAEADBs2jIqKClatWsWMGTMIhUK89dZb3HHHHaxYseKYyzxq1CjOOusstm7d2iXfk9D1RFjpAPv06UgGA75PPyO8axeR/fuJV1cTqK4msHYtABqbLRFejEOGYCgoQNI1vb1KNEastg5q69T9TUbMdjs2h4MCVy4+JUh9qL7NkUVZtiyuGPITLht8ObvrdrOudB0bKzZSG/Xyj6qv+UfV1xSY05nsVPu3OPRWFAlKbSmUDj4Hy9DpZAXqcZZuQbN/LXgOwcEN6k2jg6yxUDBF7aRrsoMcA2+5etNbICkDrOmgO/ofQpNeXffH1WyboiiEG/rEqCFGbVIKx2RkWUEBRM8pQei9lixZwuzZs7nmmmtYsGABf//737niiiu65NgrV67k6quvTgSVRmeddRbZ2dmsWrWKX/ziF4ntN910Ey+//DIzZszgb3/7G5MmTSI7O/u4yrxlyxa++OILnnzyyS75noSuJ8JKB2idTpKmTkWfmQWAHAwS3rOH8K5dhHfuJLJnD7LPR3DTJoKbNgEgGQwYBg3CNGQIxqFDMQwciMbY1NlVDoWRQ1XEKqtAAr3FSqbTQY69iKBRoj7ipj5c36JaWCNpGJo8lKHJQ5k9fDabqzazrnQd31Z/y/5gJfuDlawo+5IRSQVMdg5jnGMgRo2egBxhj8mCduAZpAy7kPSgB/PBr2D/Gqg/AIe+Vm/S/4PMMVB4hhpczE6IBqB2H9QVqxPU2TLUzrmdGM4rSVKbIaY5RVFQFFAAueG+3JBgGu8rDfvJCtDGNoWGYzTcl5Vmxz18W7NzNn9t47k5rByNgarxtSJkCV3NrDOzYdaGbj1fR+zYsYP169fzwgsvYLPZuPzyy1myZEmXhZXq6uo2wwZAdnY2VVVVLbb99Kc/5e6776auro6lS5fyi1/8Ap+v5UCGjpTZ7XbjdDqRJIn09HTmz5/PnDlzuuR7ErqeCCvHQGM2Yx45EvPIkQAosRiR4mLCO3cS2rWLyO7dyH4/4R07CO/Yob5Iq8VQWNii9kVra2h+UUD2+5H9fqKHStFoNaQmJZFhTyFqNeLWhKkL1+GP+hNl0Gv1TMycyMTMiXgiHjaWbWRt2Vr2uffxrbeYb73FmEqNTHAMZopzKEOtuYBMZcRNpRZsg84m/ZSZuPy1aA6sU4NL3T4o+0a9rX8eMkaqTUX5k9WAEqhRb1oD2NLU4GKw0hWkhunjAbT0jXlNGjUFlyOHpsaQdXhokhtSWluhSTl8G4eFpmb7GXS9v4+N0D5JkjrcLNOdlixZwpgxYxgzZgwAN954IxdeeCGHDh1Cr9cTjbaeRykajXa4mSg1NZXS0tI2nystLSUtLa3FNofDwY9+9COefPJJNm/ezMyZM/nTn/7U4TLn5OQkjtPYZ0Xo/URY6SitFkmnRYnFWz0l6XQYBw3COGgQ9osvRpFloqWlhHfuVGtfdu0iXltLZM8eInv24P3wQwD02dmJDrvGoUMTnXaVuEy83k28Xh2S7NDrSXbYka0uPKY49bK/xcgiu8HOeQXncV7BeZT5y1hXuo71ZeupDlazpnYba2q34TI4OD1lJFPsg8kxOvHFgvhiQXQaLalDzid91E8w+iph/1o1uNT8AOVb1duGFyD9FLXGJX+K2p/FfUi9GZMaOuWmgfbk/Dg1LTjYt0KWIBxNNBrltddew+fzkZmZCTR06o7HWbZsGeeccw4lJSXEYjF0zZq99+zZw5QpUzp0junTp/PWW29x//33t2gKWr16NaWlpZx33nmtXnPTTTcxffp05s2bh8HQsnn6aGVetGhRp98HoeednH9djoHWZsMyYQJyKKTWgvh86tdAACXeclSPpNFgyM3FkJtL0nnnqb8o1dWEd+0i1BBgYmVlREtLiZaWNnXaTUlJBJfmnXaVaJRYdQ1U12AD7GYTSlIKPpOCWx/DE/clRhZlWbO4YvAVXDboMn6o/4G1pWvZWL6RuoibD8vW8GHZGvJtuUxJHcNpjkE4kCkP11IersWus5Ax7CIcI69E8lWoweXAGqjaCZXb1dtXL0HaULXGpeAM9RsOe9VaGUtKw9wtzm78yQiCcKK89957eDweNm/e3GIEz/PPP8/LL7/MXXfdRUZGBosWLeK+++5Dr9ezfPlytm/fzqWXXtqhc9xxxx28/vrrzJ49m6eeeorMzEzWrl3LjTfeyO23387AgQNbvWbatGl88sknjBo1qtNlXrhwYaffB6HnibDSSRqTCY3JBI21IIqCEgwSbwwv/gBywA/N+jJIkoQuLQ1dWhrWM9Q/8HGPJ1Hrkui0W1NDYN06AuvWqeey2TAOHpwIMI2dduVgCIIhrIBVghyrFb9JwmOS8egixJQ4GknDENcQhriGMHvYbLZUbWFd2Tq2Vm3lgK+EA74S3kRiRMoIpqSNY1xSIZ54BI+/FINGT5rBTtopM9GPvAL8VU01LpXfqeGlaid8/TKkDG4ILlNAjoOvUp2YzpauBpe2JqUTBKFPWLJkCddee22rocG33347Tz31FGvWrOEf//gHd911FwMGDCAajTJy5Eg++ugjcnNzO3QOl8vFmjVrWLhwIRMmTMDr9VJYWMidd97J/Pnz23yNJElt1rh0pMyfffaZmEulDxIz2NL1s+8psowcCCT6ocg+H3Io1CLAHE4OhYjs2ZOoeYns2YMSaTkqqLHTbuNwacOgQS067QKg1RAya/EYZdyGGFFDy19Kb8TLxvKNrCtdxx73nsR2o9bIhIwJTM6YyHBrDppIACI+krVm0oyOphl0AzXQ2MelYrs6wVyj5AFNNS6OXLUTrsnZ0Ck3BcQFokeJGWy7n5jBtn8TP0OVmG6/m3THG63E4y3Di9+PHG5n8UOaddptGHEUbui024JWq05U11Dz0qLTboOwVsZnVHAb44QsWmjWrlzhr2Bd2TrWlq6lOlid2O4yujgt6zSmZE0m15QCES+mWJQ0SUeqPgmd1NCuHKyDA+vV4FK+tWVwcRY0BRdnvrpekS1dvRmTjv2NFI6ZCCvdT4SV/k38DFUirHSTnrqIK9Eost9P3OdPBBmljZ71QLuddg+nz85uGnE0dCi61NTEc+F4GL8ujtck4zdrwGoGjQZFUfih/gfWla7jq/KvCMSaVnfOS8pjcvZkTs86Hac+CU00SDJa0hUN1uYfnZBbnbdl/xoo3QxKs47I9lw1tBSeAa4iMNqaZsrVGjo1DFo4diKsdD8RVvo38TNUibDSTXrTRVwOh5tqYBoDTBsjkABi1dWJ8BLatYtYG8P/tMnJTSOOhgxBn52NpNEQk6N4Y358BhmfSUFJsoHFRFSOsrVqK2tL17K1aivxhtAhIXFKyilMyZ7C+PTxGHVGLJKODI2RZFlBE/ZCvCFohX1QsgGK10DpJnWSuUZJWU01LimD1KYhrVHt26IzqRPP6UxNj7VG0XzURXrT5/xkIcJK/yZ+hioRVrpJb7+Iy8FgIrjEfT6UQACljenr415vU6fdnTuJ7N8PcsuRShqrtUXNi6GgAFkr4Yv68MUD+MwSss2MYrPgkyJsLFfnb9lT37J/y/iM8UzJmsLwlOHoNXpSTSmkac2YYyG1eSjsVSf/iASgpGECukP/huaz85pdTf1ZLCnqXC6J+w03nUmtfdEZm91MLQPOSTpkurN6++e8PxJhpX8TP0OVCCvdpK9dxBVFQQkEiPv9yI1NSMFAqw68He60O3Cg2ml36FB0A4sIauL4ol68mihxmwnFZqVC42Nd5QbWl66nMliZeL3T6OT0rNOZnD2ZvKQ87AY7aZY0XHo7UsgNoXo1vESD6u3Q12qNy6GN6iKKR6O3tg4wrUJNsrow4+EhpjHciKYmoO99zvsDEVb6N/EzVImw0k36w0U8MQIpMYTarw5xbr5PLEZk//4W/V5addrVaBKrSxuGDEYuyiVgAm/UR8ykR7ZZ+IFK1tZtYmP5Rvyxptfn2nKZkj2F07JOI92STqo5lXRLOgatAaINNS6hegjWQ8SnTuMfqIFgbdPsuM1vHQkzAJJGXQqgrZoZSzJY0sCerc7/0jzEnGRNTf3hc97XiLDSv4mfoUqElW7SXy/iSizWou9L3OdDiTR14E102m3WdNRWp11ddjamIUPQDC4iXJhFIElPmBhRi4Et8gHWerewpebbFv1bhqcMZ3LWZCZkTCDTmkm6JR2H0dFwYkVtJgp71GaheKzha0Tt9yJH1earaKCNEHNYsAnVtxyFdCR6sxpgzG2EGluGGmjsWeoSAv2wqam/fs57MxFW+jfxM1SJsNJNTqaLuBKJqM1H/gCy39cwAqmpA2yi0+7u3YR37iTaTqdd/eBByAPyCBWmE0514DXE+Cq+l3X+b9nt25fY16A1MD59PFOypzA2bSxZtixSzCnoNUdZNyQebQovifsRtbNuItjE1KaltmplDq+tiQaPfL5GkkadG+bw5iZrGiRlqoHGkacGHr25ZS1NL3cyfc57CxFW+jfxM1R1x7Wl7/+7KHSKZDCgMxjA1bT+cWIEks+H1p6EPiO9aaZdrzcRXBIz7dbWEt/wFWz4Cg1gsVqxDSzkkgG5nF8wmYqsc1nLD6wJbqcyUsP6svWsL1uPw+DgtKzTOCP7DMakjSHZnIxBa0Cv0aPX6JGa9yvR6tVbRzTWzMjR9kNOyAO+cnWG3fZqaoK1ai1NsFa91RzhnDpjU5AxJ6vDsG0ZDYEmR705ctQ5ZU6ipiZB6ArLli3jmWeeYfPmzQAUFhbyzDPPcNlll3X6WPX19bhcLvbt20dhYWGXllPoPiKsCGiMRnUm3ORkoGkJgcbmI31WJtZTJ6DICnI4TOSHHwg1LhPwww8ofj/xrdth63ZMQIFBT0FhHlcWDuJgwWg+zahhjbIbd8TNx/s/5uP9H5Njy+HUjFNJNidjN9hJ0ifhMrlIMaVgM9jQa9UAY9AYEvf1Gj0GrQGNdNgffa2u4800cvywQNMs5ESD4KtQA423DHxVzYJMs4AT8av9aTyl6u1ImtfSWBtWqk5qaG5y5KkT5lnT1FqaftDUJPQ/bQWF4uJiioqKqKura7H+Tl/26quvcuONN/L8889z6623tnhOkiS++eYbxo4d2zOFE0RYEVqTJAnJYkFjsUDD8uyKojQ0Hfkx5OZimXgqcjCIEm3otNtsnSPZ54Nde9Hv2ssAYIBGw025GZTnZ/JVTpCPMis5xCEO+Q61eX6DxkCSISlxsxvsLR47DA412JhTSDGlYDfY0Wv16DS6FqFGp2nj463Rqjd9e1W2I5vuynJDkIm0DDghtxpSfOXgLW8KOM1raoK1apNVqF691e5p53yoo5UStTNpYMtU56NJagw0DaHGJJpuhL5PURRkWW6xwnJvsGTJEpKTk1myZEmrsCL0PBFWhA6RJAmtzYrWZk1sU+Jx5EAA48ABWCdNVDvxBoLEysrUmpeGpqN4TQ3aA2XkHIDLUW++NCv7cwzU2KEySaHCGqXEGqE6ScZjCVMTilATOlI7TJMjhRuXyYXL6CLZlEyqOZVUcypJhiQ11GjVmhudRteyCaqRRgOatvqj5EHmyJabFKVls1M8ooYY9yG1lsZb3tQM5a9qCjaNnYx9Dc9XHOEbNdrVQGNNh6QMNczYs9V1mBz56hIHtjQxTLuPaqzR7C6S2dz25/4YTJ06lcmTJ7Np0ybWrl3L4MGDeeWVVxKrIhcWFjJ37lzeffddtmzZwldffYXZbOZXv/oVGzduxOVyMX/+fBYsWNCh861atYqFCxeya9cucnJyePzxx5k5cyYA4XCYBQsW8Oabb+JwOFi0aNFRj7d7926++OIL3nnnHS6//HK2bNnCmDFjjvn9ELqeCCvCMZO0WrRJSWiTmtb6aRyBZD1tUkNHXj/RsvKm4dINnXZtVX5GVPnbPK6i0RC1mwjbjfiTDHiTdNQnSdTYoDJJpjwpziFrmCpdkJgSIyKrwabD4UZrIEnfMtg4jA5cRhdOk5NkUzIpphQ13FhScRgciaaoVk1QiTdDaph919C0LSkTstq54DUGm7AH3KXgLgFvaVOo8VaoocZfrQYbOaruG/ZA7d72vzlbJty1s0Pvg9C7KMEgO8dP6LbzDd30bySLpcuO99prr/HBBx8wYsQIbrvtNubPn8/nn3+eeH7ZsmW89957DBo0iHg8zujRo5k5cybvvvsuu3bt4sILLyQ9PZ1Zs2Yd8Txbt27lqquu4u2332bq1KmsXbuWiy++mK+++oqhQ4fy6KOPsm7dOrZt24bFYjnq8QBefvllxo0bx6WXXspZZ53FkiVL+N///d/jfUuELiTCitClJJ0OrcOB1uGgsXusPCKCdfLpiT4w0fJyQju+I3rgAPH6emJ1dcTr6ojX1yN7PEiyjKHej6HeTxKQ2d7JjAYUh4Oow0LQbsSfpMOdpKPWrqHWLlFhlym3xnBLQTwxnzpXjBIjEo9QEz+2cOMwOnAYHDhNzkSNTYo5hRRzCmnmNNIsadgN9raboJpr7EBssKo1JLlH+CMVizbU0hwEzyH15i1TA02itqZabW6ypHToexKErnbdddclaiNuvPFGLrzwwhbP33rrrQwdOhSA9evXU1ZWxiOPPILBYGD06NHMmzePZcuWHTVcvPjii8yZM4dp06YBcOaZZ3LJJZewYsUK7r33XpYvX87jjz9OdnY2APfffz8ffvhhu8eLx+O88sor/Pa3vwXghhtu4D//8z956qmnMB6+qr3QY0RYEU44jcGAptkIJNPQoVhPO03t8xKJokSjKNEIxGLIwSCxqipi1dXEqmvUENMQZJqHGiUYhHAEqbIGQ2UNBsABZLdxfikpCclpR3KmE7dbCDjN+J1G3E49NQ4NVUkKlZYobtmPN+bHG/PhiXjwRrzE5GMPNw6jA4fRgdPoxGVykWxUg01jc1RjuLHoj/LfrU4Pzlz1diSxhpoaoU+SzGaGbvp3t56vI/R6PdHDFlhtfKzXN43Yy8xs+rfCarXi8/lavCY/Pz9xv6SkhOzsbAyGpprIAQMG8Prrrx+1PMXFxXz66acsXbo0sS0WiyWGzJaWllJQUJB4rvn9tvzjH/+guro6EZKuuuoq5s+fz9/+9jeuueaao5ZH6B4irAg9QmMyoTnKvASKLDcEmaYbDV/jXi+xikqilRVqqKmubhVo4nV1EI+jeL0oXi8cVDv0WhpuacCgxpPpdGidTrROJxqXC03KAEhJJppix+sy4nYZqHZqqNGHqIt5cUfc1IfrcYfdeCKe4w43doMdp9GZuCWb1aaoFFMKqRY12KSYU0g2JWPWtfNHRmcAXWrbzwm9XmPH9t6moKCAffv2tdi2Z88eUlNTsVqt7byqNU2zofu5ubmUlpYSjUYTgae4uJjc3KMEciAvL49f//rXPPHEE20+n52dzf79+znttNMAOHDgwBGPt2TJEmRZTvSvATWMLVmyRISVXkSEFaHXkjQaJKMR2quKbbq2qJ0To1F1ht6YGmjkcJh4TQ3RykpiFZXEqiqJVdcQq6kmXlffoumJWEwNPNXVrU5jarhloP43qnW50LlcaNNS0acORJeRjjYjHSkvjWCanXqnllrFT224lppgDXWhOvUWrsMdduOOuHGH3XgjXqJylEg8QnWwmupg63O3xag1JmpsnEa1j01jB+JsWzYXD7i40++1ILTnuuuu4/7772f69OmMHTuWAwcO8MADDzB79uxjPuakSZPIyMjgvvvu44EHHmD37t08++yz/O53vzvqa+fOncuFF17IBRdcwNlnn00sFmPTpk04nU6GDx/OtddeyxNPPMGZZ56JxWLhoYceavdYFRUVfPDBB7z66quJZiWAzZs386Mf/Yji4mIxN0svIcKK0C9IkoRkMECzamUAcnM5vB5CUZREDY0SiyEHAkQrKoiVlxOtrCJWWUmsqop4TQ2xmppELY0SCqEEg8SCQWJtzOzbvCzpdjtZycnoUlPRpqaiz8hAlzkQfVYWuoHZGHJy0TgdBGIBKoOV1ARrqA5WUxuspSakBpzacC31IbX2pj5cjyfiISpHCcfDVAYqqQxUtjp3ljVLhBWhS9144414vV5mz55NSUkJqampXHnllTzwwAPHfEy9Xs/777/PvHnzyMzMxOVyceedd3aoM+y4ceP485//zOLFi/nuu+/QaDSMHTuWp59+GoDFixdTWVnJyJEjsdvtLF68mA8++KDNY73yyivk5+dzzTXXtKj5ufDCCxk/fjwvv/zyEcOO0H3EdPuIaciFI2tsgoq73UTLyoiWlxOrqCDWGGxqatRgU1tLvL4e4vGOHVivR+dyqYEmLRVdegb6jAz0WZnos7PRZ2ejS09H09C3QFEUArEAtcFatdYmUENVqCoRcOpD9SSbk1l42sI2Tyc+591PTLffv4mfoUpMty8IvYCk1yPp9WgsFvRZWUfcV45EiFVWEi0vJ1pWRqy8oqGmppJYVXUi2MheL0Sj6nOVrWtImtNYrWhTU9GlpqLLSMeQnkFudhZFWVnoM09Bl5WOLiUFSSd+nQVB6J/E1U0QupDGYMCQm4vhKB0F5UhEraEpLVWDTXm5GmwaQk28Wg02SuO6TX4/0f37j3BiDVqXC+PgwRQsW9r+foIgCH2QCCuC0AM0BgPG/HyMzYZzHk5RFGSfT62pKWsINGVlav+aiorEEO94bS3IstoU1U/WaREEQWhOhBVB6KUkSUrMEGwcOLDd/ZR4nHhtLdGKSnW+GkEQhH5GhBVB6OMkrRZdWhq6hkUnhb5FjHHou2RZ7ukinDREWBEEQegBer0eSZKoqqoiLS2tyxYVFE48RVGIRCJUVVWh0WhazMQrnBgirAiCIPQArVZLbm4uJSUlFBcX93RxhGNgsVjIz89vMUeLcGKIsCIIgtBDbDYbgwcPbrX2jtD7abVadDqdqBHrJiKsCIIg9CCtVotWq+3pYghCr9ZjdVfRaJR58+bhcrlITk5m/vz5xGKxY9q3M8cSBOHk1pnrxaFDh7jssstISUkhNTWVq6++mqqqqm4usSAIPRZWHnnkEVavXs2OHTvYvn07X375JY899tgx7duZYwmCcHLrzPXiV7/6FQD79+9n3759hEIhbr/99u4sriAI9ODaQHl5efz+97/nyiuvBOCtt97irrvuYn8bs3Qebd/OHAsgHA4TDocTj91uN/n5+Rw8eFCsmSL0Wx6Ph7y8POrr63E4HD1dnB7TmevF6NGjufvuuxML7C1fvpzHH3+cbdu2tXlscW0RTkbdcm1RekBtba0CKLt3705s27VrlwIo9fX1ndq3M8dqdP/99yuAuInbSXk7ePBgF/9G9x2dvV4sXbpUueyyy5T6+nqlrq5Oufjii5W777673eOLa4u4ncy3E3lt6ZEOtj6fDwBns6nBG+97vd4Wyexo+yoNFUMdOVaje+65hzvvvDPxWJZlamtrSUlJabdnd2NyFP8hdQ3xfna9o72niqLg9XrJzs7ugdL1Dp259gCcccYZ/OEPf8DlcgEwefJk7rnnnnaPL64tPU+8n12vN1xbeiSs2Gw2QK0iTU1NTdwHSEpK6tS+jTMIduRYjYxGI0ajscU2ZwfXVLHb7eIXoAuJ97PrHek9PZmbf6Bz1x5Zlpk+fTpXX301n3zyCQAPPPAAM2bMYP369W0eX1xbeg/xfna9nry29EgHW5fLRW5uLps3b05s27x5M3l5ea2+4aPt25ljCYJwcuvM9aK2tpb9+/dz++23Y7FYsFgszJ8/nw0bNlBdXd3NJReEk1uPjQa66aabePTRRykvL6e8vJzHHnuMm2+++Zj27cyxBEE4uXX0epGamsqgQYN47rnnCIVChEIhnnvuOXJzcxO1MoIgdI8emxTu3nvvpaamhuHDhwNw3XXXsXDhQgB++ctfAvDCCy8cdd+OPN8VjEYj999/f6sqXuHYiPez64n3tGM6c+159913ueOOO8jJyUGWZcaNG8d7773XpeURP7euJd7Prtcb3tMeG7osCIIgCILQEWL1JUEQBEEQejURVgRBEARB6NVEWBEEQRAEoVcTYUUQBEEQhF5NhBVBEARBEHo1EVYOEw6HueWWWygqKiIpKYlhw4bx8ssvJ573eDzMmjULu91ORkYGDz/8cA+Wtm+YP38+eXl52O12cnJyWLBgAZFIBBDv5/EKBoMMGjSoxSyp4j3tncS1peuJa8uJ0+uuLSds1aE+yufzKffee6/yww8/KLIsK+vWrVOcTqeycuVKRVEU5YYbblAuuOACpa6uTtm5c6eSl5envPLKKz1c6t5tx44dis/nUxRFUaqqqpSpU6cqDz/8sKIo4v08XnfddZcybdo0xeFwJLaJ97R3EteWrieuLSdOb7u2iLDSAZdffrly7733Kn6/XzEYDMrGjRsTz/3ud79Tzj777B4sXd9SWVmpTJs2TbnhhhvE+3mcvv76a2XkyJHKypUrExcU8Z72LeLa0nXEtaXr9MZri2gGOopQKMRXX33F6NGj2blzJ5FIhLFjxyaeHzt2LFu3bu25AvYRTzzxBDabjfT0dLZs2cL8+fPF+3kcYrEYt9xyC8899xwGgyGxXbynfYe4tnQNcW3pWr312iLCyhEoisLNN9/M4MGDueKKK/D5fFitVnS6plUKnE4nXq+3B0vZN9x99934fD527NjBL3/5SzIzM8X7eRyeeuopxo0bx9lnn91iu3hP+wZxbek64trStXrrtUWElXYoisJtt93Gzp07eeedd9BoNNhsNgKBALFYLLGf2+1utbS80L7hw4czZswY5syZI97PY/TDDz/wwgsv8NRTT7V6TrynvZ+4tpwY4tpy/HrztUWElTYoisKvfvUrNmzYwMcff5xYOn7o0KHo9Xq2bNmS2Hfz5s2MGjWqp4raJ0WjUXbv3i3ez2O0evVqKioqGDJkCKmpqVx66aV4PB5SU1PxeDziPe3FxLXlxBLXluPTq68t3dIzpo+57bbblNGjRyvV1dWtnrv++uuViy66SKmvr1d27dql5Ofnix7mR+D1epWXX35ZqaurU2RZVrZu3aoMHz5cueWWWxRFEe/nsfD7/crBgwcTtxUrVih2u105ePCgEg6HxXvai4lrS9cR15au15uvLSKsHKa4uFgBFKPRqFit1sRt7ty5iqIoitvtVq655hrFZrMpaWlpyoMPPtjDJe7dfD6fcv755yvJycmK1WpVioqKlLvuukvx+/2Kooj3syt89tlnLYYXive0dxLXlq4lri0nXm+6tkiKoijdU4cjCIIgCILQeaLPiiAIgiAIvZoIK4IgCIIg9GoirAiCIAiC0KuJsCIIgiAIQq8mwoogCIIgCL2aCCuCIAiCIPRqIqwIgiAIgtCribAi9Dv79+9nyJAhhMPhNp/fvHkzkiR1+Hher5eBAwdSXV3dVUUUBKEPEteWniPCSj8xdepUnnnmmW4515w5c1iwYMER9ykrK2PWrFlkZmaSlJTEgAEDuOOOO1ocQ5Ik3n///RavczqdfP755wB8/vnnSJKEzWYjKSmJgoIC7rnnHmRZPuK577vvPubPn4/RaOzw92MwGLDZbNhsNgYOHMhzzz2XeD4pKYkbbriBRx99tEPHE4T+RFxbmohrS88RYUU4Ia6//npMJhPff/89brebTz75hLFjx7bYJyUlhYULFx7xAuFwOPD5fHi9Xj766COWLVvGsmXL2t2/pqaGv/71r8yePbtT5b3tttvw+Xz4fD7eeOMNfvvb37JmzZrE8zfeeCNLly4lEAh06riCIHQtcW05OYmw0g99/vnnOJ1O/vjHP5KXl0dKSgr/+Z//mXh+2bJljB07loULF5KSkkJ+fj7PP/984vnD/7upr69HkiSKi4v53//9X5YvX87zzz+PzWZjxIgRbZZh/fr13HTTTTidTjQaDQMHDuTGG29ssc9Pf/pTAoEAr7/+eoe+r+HDh3PmmWfy73//u919Vq5cyfDhw0lOTm5R/quvvhqn08mwYcP44osvjnieiRMncsopp7B9+/bEtsLCQlJSUvjXv/7VobIKQn8kri3i2tJTRFjpp7xeLzt27GD37t2sXr2a5557LlEFCrBt2zYkSaKsrIw333yTu++++6i/aAC33347s2fPTvy30PyXrrkzzjiDBQsW8Oqrr7Jr164299Hr9Tz88MPcd9997bYBN/ftt9/yxRdfMGTIkHb32bx5M8OGDWtV5vr6eoqLi/n000959dVX2329oiisW7eOnTt3cvrpp7d47pRTTmHz5s1HLacg9Gfi2tKyzOLa0j1EWOmnFEXhkUcewWQyMXz4cKZMmdLivwar1coDDzyAwWBg8uTJzJ49+4i/aJ311ltv8eMf/5hnnnmGESNGUFBQwJ/+9KdW+11zzTUkJyfzf//3f20ex+1243Q6MZvNjB49mh/96Efceuut7Z63rq4Ou92eeByPx3nzzTd55JFHcDqdZGdn8x//8R+tXvd///d/OJ1ObDYbU6ZM4cYbb2TUqFEt9rHb7dTV1XX0LRCEfklcW1Ti2tK9RFjpp+x2OxaLJfHYarXi9XoTj7Ozs9Hr9YnHBQUFHDp0qEvP/8ADD7Bp0ybq6uq4/fbbueGGG/juu+9a7CdJEk888QSPPvpoi/I1cjgc1NfX4/P5eOmll1i9ejXBYLDd87pcLjweT+JxdXU1kUiEgoKCxLbm9xvdeuut1NfX4/f7KSkpYevWrSxcuLDFPh6PB5fL1eH3QBD6I3FtUYlrS/cSYeUkVVpaSjQaTTw+cOAAOTk5ANhsthadvcrKylq8VqPp3MfGZrPxm9/8BofDwY4dO1o9P2PGDMaMGcNTTz3V7jG0Wi233HILp5xyCg888EC7+40dO5bvv/8+8Tg1NRW9Xs/+/fsT2w4cOHDE8ubk5HDllVe2Gk2wY8eOVh35BEFoSVxb2ieuLcdOhJWTlN/v5+GHHyYSibBhwwaWL1+e6OU+fvx4Vq5cSVlZGV6vlwcffLDFazMyMti7dy+KorR7/P/4j/9g8+bNRCIRIpEIf/zjH/H7/UyYMKHN/Z944gmeeeYZQqHQEct977338sILL7T7n9qMGTP47rvvElWqWq2Wq6++mvvuu4/6+npKS0uPeOECqKys5O23325RVbt//36qq6s5++yzj/haQTjZiWtL+8S15diJsHKSGjlyJLFYjKysLK688koeffRRzj33XACuu+46zjnnHIYNG8bYsWO5+OKLW7z25ptv5tChQyQnJzN69Og2jx8Oh7nmmmtISUkhMzOTpUuX8u6771JYWNjm/qeeeioXXXTRUTvDnXrqqZx99tntzkuQmprK5ZdfzvLlyxPbnn32WWw2GwUFBUybNo3rr7++1esaRyDYbDZGjx5NXl4ezz77bOL5V199lTlz5mC1Wo9YPkE42YlrS0vi2tJFlF7q2WefVSZMmKAYDAbl0ksvPeK+brdbufbaa5WkpCQlPT1deeihh7qnkH3U0qVLlTFjxvR0MU6Yffv2KYMHD1ZCoVCXHM/j8SgDBw5UKisru+R4Qs8S15YTR1xbOkdcWzpO19NhqT3Z2dksXryYVatWUVJScsR958+fT21tLQcOHKCyspLzzz+fgoICbrjhhm4qrdCbFBYWtjuk8VgkJSXxww8/dNnxhJ4lri3CsRLXlp7Ta8PKFVdcAahj2490QQkEArzxxhusWbMGp9OJ0+lk/vz5LFmypN0LSjgcblElKMsytbW1pKSkdGpdh74qGAwSj8db9GwX+j9FUfB6vWRnZ3e6I2N/Iq4tJ464tpycuuPa0mvDSkft3LmTSCTSoif12LFjeeyxx9p9zeOPP96qY9fJyOFw9HQRhB5w8OBBcnNze7oYvZ64thw7cW05OZ3Ia0ufDys+nw+r1YpO1/StOJ3ONsfVN7rnnnu48847E4/dbjf5+fkcPHiwxaQ/gtCfeDwe8vLySEpK6umi9Ani2iIIHdMd15Y+H1Yax+3HYrHERcXtdh/xTTMajW2ummm328UFRej3TobmiK4gri2C0Dkn8trS5xuuhw4dil6vZ8uWLYltmzdvbjWdsSAIQmeIa4sg9B69NqzEYjFCoRCxWAxZlgn9//buP6auu/7j+Ktb25vBvdx7GVVx3JIlax3GtrdGjYa4zqibiy6bTe2MQSuRtoRCjbPGQceaTbjoNKnLVsMfttH4h27EbM6YzK7/gWlwmV5jILbrVKA/2A8K5wLSO7Cf7x8N1/IFCrSXe9+nPB/JSe49530P73Nz+s7r/uq5dEnvvfferLqCggI98sgjam5ulud5euONN/Tss8+qpqYmD10DsI7ZAviP2bDS0tKi2267Ta2trfr973+v2267Tffdd58k6YEHHpjxJbfnnntO4XBYZWVlqqys1Le+9S1+WghgTswWwH9WOXeN/9d4hUilUgqHw/I8j8+VcdPiPM89nnOsBLk4z82+swIAACARVgAAgHGEFQAAYBphBQAAmEZYAQAAphFWAACAaYQVAABgGmEFAACYRlgBAACmEVYAAIBphBUAAGAaYQUAAJhGWAEAAKYRVgAAgGmEFQAAYBphBQAAmEZYAQAAphFWAACAaYQVAABgGmEFAACYRlgBAACmEVYAAIBphBUAAGAaYQUAAJhGWAEAAKYRVgAAgGmEFQAAYBphBQAAmEZYAQAAphFWAACAaYQVAABgGmEFAACYZjasTE5Oqr6+XtFoVMXFxWpoaNDU1NSctefOndPDDz+s22+/XSUlJdq5c6feeeedHHcMwA+YLYD/mA0rLS0t6urqUm9vr3p6etTZ2alEIjFn7b59+yRJfX19+te//qVLly5p//79uWwXgE8wWwD/MRtWjh07pscff1ylpaUqLS3VwYMHdfTo0Tlr//nPf2rnzp0KBoMKhUJ65JFH9Pe//33efafTaaVSqRkLgJWB2QL4j8mwMjw8rLNnzyoej2fWxeNx9ff3y/O8WfWPPvqoOjo65HmeRkZG9Otf/1oPPvjgvPtva2tTOBzOLLFYbDkOA4AxzBbAn0yGlbGxMUlSJBLJrJu+PTo6Oqu+srJSb7/9duYz6OHhYTU2Ns67/8bGRnmel1kGBgay2j8Am5gtgD+ZDCvBYFCSZrzSmb4dCoVm1F6+fFmf//znVVlZqbGxMY2NjamyslL33XffvPsPBAIqKiqasQC4+TFbAH8yGVai0ajKysqUTCYz65LJpGKxmMLh8Izaixcvqq+vT/v371dBQYEKCgrU0NCg7u5uvfvuuznuHIBlzBbAn0yGFUmqrq5Wa2urBgcHNTg4qEQioZqamll1JSUluuuuu3TkyBFdunRJly5d0pEjR1RWVqaSkpI8dA7AMmYL4D+r893AfJqbmzU0NKSKigpJUlVVlZqamiRJtbW1kqT29nZJ0u9+9zt95zvf0R133KHLly9r69atevnll/PTOADTmC2A/6xyzrl8N5FvqVRK4XBYnufxGTNuWpznucdzjpUgF+e52Y+BAAAAJMIKAAAwjrACAABMI6wAAADTCCsAAMA0wgoAADCNsAIAAEwjrAAAANMIKwAAwDTCCgAAMI2wAgAATCOsAAAA0wgrAADANMIKAAAwjbACAABMI6wAAADTCCsAAMA0wgoAADCNsAIAAEwjrAAAANMIKwAAwDTCCgAAMI2wAgAATCOsAAAA0wgrAADANMIKAAAwjbACAABMI6wAAADTCCsAAMA0wgoAADCNsAIAAEwjrAAAANNMh5XJyUnV19crGo2quLhYDQ0Nmpqamrf+5ZdfVjweV2FhoT74wQ+qvb09h90C8AtmC+AvpsNKS0uLurq61Nvbq56eHnV2diqRSMxZ+8orr6iurk4//elPlUql1NPTo3vvvTe3DQPwBWYL4C+rnHMu303MJxaL6fDhw9qxY4ckqaOjQwcOHFBfX9+s2o9//OPavXu39uzZs+B+0+m00ul05n4qlVIsFpPneSoqKsreAQCGpFIphcNhznMxW4BsysVsMfvOyvDwsM6ePat4PJ5ZF4/H1d/fL8/zZtSOj4/r9ddf17lz57Rx40Z94AMf0Fe+8hVduHBhzn23tbUpHA5nllgstpyHAsAQZgvgP2bDytjYmCQpEolk1k3fHh0dnVE7PDws55xeeuklvfrqqzpz5owCgYCqqqrm3HdjY6M8z8ssAwMDy3IMAOxhtgD+szrfDcwnGAxKkjzPU0lJSea2JIVCoTlr9+/fr/LycknSk08+qQ0bNmh8fFyFhYUz6gOBgAKBwLL2D8AmZgvgP2bfWYlGoyorK1MymcysSyaTisViCofDM2ojkYjWr18/534MfyUHQB4wWwD/MRtWJKm6ulqtra0aHBzU4OCgEomEampq5qzds2ePnn32WZ07d04TExN66qmn9NnPfjbzyggApjFbAH8x+zGQJDU3N2toaEgVFRWSpKqqKjU1NUmSamtrJSnz/x089thjunjxorZs2SJJ+sxnPqNf/epXeegagHXMFsBfTP90OVf4SSdWAs7z3OM5x0qwon+6DAAAIBFWAACAcYQVAABgGmEFAACYRlgBAACmEVYAAIBphBUAAGAaYQUAAJhGWAEAAKYRVgAAgGmEFQAAYBphBQAAmEZYAQAAphFWAACAaYQVAABgGmEFAACYRlgBAACmEVYAAIBphBUAAGAaYQUAAJhGWAEAAKYRVgAAgGmEFQAAYBphBQAAmEZYAQAAphFWAACAaYQVAABgGmEFAACYRlgBAACmEVYAAIBphBUAAGAaYQUAAJhmOqxMTk6qvr5e0WhUxcXFamho0NTU1DUfMzExobvuukuRSCQ3TQLwHWYL4C+mw0pLS4u6urrU29urnp4edXZ2KpFIXPMxTzzxhMrLy3PUIQA/YrYA/mI6rBw7dkyPP/64SktLVVpaqoMHD+ro0aPz1r/++ut65ZVX9P3vf/+a+02n00qlUjMWACsHswXwF7NhZXh4WGfPnlU8Hs+si8fj6u/vl+d5s+qnpqa0e/duHTlyRGvXrr3mvtva2hQOhzNLLBbLdvsAjGK2AP5jNqyMjY1J0ozPh6dvj46Ozqr/8Y9/rK1bt+qee+5ZcN+NjY3yPC+zDAwMZKVnAPYxWwD/WZ3vBuYTDAYlSZ7nqaSkJHNbkkKh0IzaM2fOqL29XX/9618Xte9AIKBAIJDFbgH4BbMF8B+z76xEo1GVlZUpmUxm1iWTScViMYXD4Rm1XV1deuutt7Rx40aVlJTooYceUiqVUklJibq7u3PcOQDLmC2A/5h9Z0WSqqur1draqsrKSklSIpFQTU3NrLqdO3fqc5/7XOb+yZMnVVNTo2Qyqfe973056xeAPzBbAH8xHVaam5s1NDSkiooKSVJVVZWampokSbW1tZKk9vZ2FRQUqKCgIPO4devWadWqVSorK8t90wDMY7YA/rLKOefy3US+pVIphcNheZ6noqKifLcDLAvO89zjOcdKkIvz3Ox3VgAAACTCCgAAMI6wAgAATCOsAAAA0wgrAADANMIKAAAwjbACAABMI6wAAADTCCsAAMA0wgoAADCNsAIAAEwjrAAAANMIKwAAwDTCCgAAMI2wAgAATCOsAAAA0wgrAADANMIKAAAwjbACAABMI6wAAADTCCsAAMA0wgoAADCNsAIAAEwjrAAAANMIKwAAwDTCCgAAMI2wAgAATCOsAAAA0wgrAADANMIKAAAwjbACAABMI6wAAADTTIeVyclJ1dfXKxqNqri4WA0NDZqamppVl06ntXv3bt15550KhUK6++67dezYsTx0DMAPmC2Av5gOKy0tLerq6lJvb696enrU2dmpRCIxq25qakqlpaU6ceKEUqmUfvGLX+i73/2ujh8/noeuAVjHbAH8ZZVzzuW7ifnEYjEdPnxYO3bskCR1dHTowIED6uvrW/Cx27dv10c+8hE99dRTs7al02ml0+nM/VQqpVgsJs/zVFRUlL0DAAxJpVIKh8Oc52K2ANmUi9li9p2V4eFhnT17VvF4PLMuHo+rv79fnudd87GXLl3Sn//8Z23evHnO7W1tbQqHw5klFotls3UAhjFbAP8xG1bGxsYkSZFIJLNu+vbo6Oi8j3POqaamRhs2bND27dvnrGlsbJTneZllYGAga30DsI3ZAvjP6nw3MJ9gMChJ8jxPJSUlmduSFAqF5nyMc051dXU6deqUTpw4oVtumTuLBQIBBQKBZegagHXMFsB/zL6zEo1GVVZWpmQymVmXTCYVi8UUDodn1TvntG/fPnV3d+v48eNz1gAAswXwH7NhRZKqq6vV2tqqwcFBDQ4OKpFIqKamZs7a+vp6/elPf9Krr76qaDSa404B+AmzBfAXsx8DSVJzc7OGhoZUUVEhSaqqqlJTU5Mkqba2VpLU3t6uvr4+/exnP1MgEFB5eXnm8VVVVWpvb8994wBMY7YA/mL6p8u5wk86sRJwnucezzlWghX902UAAACJsAIAAIwjrAAAANMIKwAAwDTCCgAAMI2wAgAATCOsAAAA0wgrAADANMIKAAAwjbACAABMI6wAAADTCCsAAMA0wgoAADCNsAIAAEwjrAAAANMIKwAAwDTCCgAAMI2wAgAATCOsAAAA0wgrAADANMIKAAAwjbACAABMI6wAAADTCCsAAMA0wgoAADCNsAIAAEwjrAAAANMIKwAAwDTCCgAAMI2wAgAATCOsAAAA0wgrAADANLNhZXJyUvX19YpGoyouLlZDQ4OmpqZuuBbAysZsAfzHbFhpaWlRV1eXent71dPTo87OTiUSiRuuBbCyMVsA/1nlnHP5bmIusVhMhw8f1o4dOyRJHR0dOnDggPr6+m6oVpLS6bTS6XTmvud5Wr9+vQYGBlRUVLQMRwPkXyqVUiwW08jIiMLhcL7byRtmC5BdOZktzqCLFy86Se6NN97IrDt9+rST5EZGRq67dtqhQ4ecJBaWFbm8+eaby/MP1weYLSwsy7cs52xZLYPGxsYkSZFIJLNu+vbo6OiM5LaU2mmNjY169NFHM/dHRkZUXl6u/v5+373inE60fnvl5te+Jf/2Pv0qv7i4ON+t5A2zZfH8ep77tW/Jv73nYraYDCvBYFDSlSegpKQkc1uSQqHQdddOCwQCCgQCs9aHw2FfnSBXKyoq8mXvfu1b8m/vt9xi9qtqy47ZsnR+Pc/92rfk396Xc7aYnFrRaFRlZWVKJpOZdclkUrFYbNark6XUAljZmC2AP5kMK5JUXV2t1tZWDQ4OanBwUIlEQjU1NTdcC2BlY7YA/mPyYyBJam5u1tDQkCoqKiRJVVVVampqkiTV1tZKktrb2xesXYxAIKBDhw7N+fatdX7t3a99S/7t3a99ZxuzZXH82rtf+5b823su+jb702UAAADJ8MdAAAAAEmEFAAAYR1gBAACmEVYAAIBphBUAAGDaTRtWsnkZ+FxfJn6xfy+dTmv37t268847FQqFdPfdd+vYsWMzau69914FAgEFg8HMcv78+bz2LUnf/OY3tXbt2hl9nTx58rr2lever+45GAxqzZo12rx586KPLZuee+45fexjH1MgENDDDz98zdpUKqWvfe1rKioq0vvf/3794Ac/WNJ2XMFsuYLZkv3emS3XsGxXHcqzJ554wm3ZssWdP3/enT9/3m3ZssU9+eST11W7lH3lsvexsTHX3Nzszpw54y5fvuxOnjzpIpGI++Mf/5ip2bZtmzt8+PCy9Xo9fTvn3K5du9y3v/3trOwrG27k723atMm1trZm7i90bNn029/+1r344otu37597qGHHrpm7Te+8Q13//33u+HhYXfq1CkXi8XcL3/5y0VvxxXMliuYLYvDbMnObLlpw0pZWZnr6OjI3H/hhRfc+vXrr6t2KfvKhhv5e1/+8pddc3Nz5n4uB8pS+l7oH51fnvPu7m536623unPnzmXW5XKgTDt06NA1B8r4+Lhbu3ate+211zLrnn76aXfPPfcsajv+h9lyBbNlcZgt2ZktN2VYyeZl4K/nMvG56v3/m5iYcHfccceMfxjbtm1zt99+u4tGoy4ejy/bK+Wl9r1r1y4XjUZdNBp1H/7wh91PfvIT99///ve69pXr3q+2Z88e96UvfWnGumsd23JZaKD85S9/cZLc5ORkZt3x48ddJBJZ1HZcwWxhtixn71djtsx0U35nZaFLuy+ldin7yobr/XvOOdXU1GjDhg3avn17Zn1bW5vefPNNvfXWW/rhD3+ohoYGvfjii3nve//+/Tp16pTeeecdHT16VM8884yeeeaZ69pXrnufNj4+rt/85jezrhVzrWPLl7GxMRUWFmr16v9dYSMSiWSOb6HtuILZwmxZzt6nMVtmuynDytWXdp+2mMvAz1W7lH1lw/X8Peec6urqdOrUKb300kszLtP9qU99SuFwWGvWrNH999+vvXv36vnnn8973x/96Ee1bt063XrrrfrkJz+pxx57LNOXH55zSero6FBBQYG++MUvzlh/rWPLl2AwqP/85z8zvtjneV7m+BbajiuYLcyW5ex9GrNltpsyrGTzMvC5vkz8Uv+ec0779u1Td3e3jh8/vmBPVw+bbLrR5+nqvqw/59N+/vOfa9euXTNeMcxluZ7zpfjQhz6kNWvW6G9/+1tmXTKZ1KZNmxa1HVcwW+bHbJmN2ZLF2bKkD418pLm52W3dutVduHDBXbhwwW3dunXeb2AvVLuUfeW697q6Ord582b37rvvzto2PDzs/vCHP7jx8XE3NTXlTpw44cLhsHvhhRfy3vfzzz/vPM9zly9fdq+99porLy93Tz/99HXtK9e9O+fcP/7xD7dq1Sp3+vTpWdsWOrZsmpycdBMTE+7gwYPuwQcfdBMTEy6dTs9Z+/Wvf9098MADbmRkxJ0+fdqtX79+xvcMFtqOK5gtzJbl6t05Zst8btqw8t5777m6ujoXiURcJBJx9fX1mS/47N271+3du3dRtYvZnq/e//3vfztJLhAIuMLCwswyvf3tt992n/jEJ1woFHKhUMht2rTJHT16NO99O+fcpz/9aRcOh11hYaHbuHGj+9GPfjTji2JWn/Np3/ve9+b9NvtCx5ZNhw4dcpJmLNu2bXPOOfeFL3xhxs8ePc9zX/3qV10wGHTr1q2bNTAX2o4rmC3MluXq3Tlmy3xWOefcDbwDBAAAsKzy/4EXAADANRBWAACAaYQVAABgGmEFAACYRlgBAACmEVYAAIBphBUAAGAaYQUAAJhGWAEAAKYRVgAAgGmEFQAAYNr/AcjUhRgJwmi6AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 600x600 with 4 Axes>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e598e7c0-0568-4295-b7cb-6f0ce368d35a",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Low SNR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "66119f54-77f1-423e-bfeb-58c90affa07a",
   "metadata": {},
   "outputs": [],
   "source": [
    "test_no = 125\n",
    "kspTestScaled = ksp_fs_sc_all[:,:,:,:,test_no]\n",
    "volumeGTscaled = img_gt_sc_all[:,:,:,test_no]\n",
    "mask_metrics = np.ma.getmask(np.ma.masked_less(abs(volumeGTscaled),0.05)) # a mask in the image domain for calculating image metrics (not for masking kspace)\n",
    "\n",
    "l = 0\n",
    "noise_val = noise_dB[l]\n",
    "lamda_val = lamda_opt[l]\n",
    "iter_val = iter_opt[l]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "8b1fb15d-d3ad-4ecf-b1b7-90fa1ba37d67",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "60"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "noise_val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "4fdeff84-bbd0-4936-a7f8-a6365482b855",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1, 75, 64, 25)\n"
     ]
    }
   ],
   "source": [
    "kspInput = awgn(kspTestScaled,noise_val,L=1)\n",
    "uskspace = applyMask(kspInput,mask)\n",
    "\n",
    "\n",
    "volcsrecon = ulfl1recon(uskspace,mask,lamda_val,iter_val,mps)\n",
    "    \n",
    "volifftrecon_mc = sp.ifft(uskspace,axes=[1,2,3])\n",
    "volifftrecon = coil_combine(volifftrecon_mc,mps)# np.sqrt(np.sum(np.abs(volifftrecon_mc)**2,axis=0)).reshape(imSize)\n",
    "\n",
    "amap_scale = 1\n",
    "unroll_scale = 1#0.25#1/8\n",
    "volunrollrecon = unrolling_fns.unrollingRecon(uskspace*unroll_scale,mask,unroll_model_pth)/unroll_scale\n",
    "#volunrollrecon = np.roll(volunrollrecon,(-1,-1),axis=(0,2))\n",
    "volamaprecon, volamaprecon_mc = automap_fns.automapRecon(uskspace*amap_scale,mps,model_real_dir,model_imag_dir)\n",
    "volamaprecon = volamaprecon/amap_scale\n",
    "volamaprecon_mc = volamaprecon_mc/amap_scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "dc9a3832-04f8-46f2-b511-39b60d70a9e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # calculating SNR of input kspace (fully sampled)\n",
    "# fs_recon = sp.ifft(kspInput,axes=[1,2,3])\n",
    "# im_slice = fs_recon[0,:,32,:]\n",
    "# fig=plt.figure()\n",
    "# plt.imshow((abs(im_slice)))\n",
    "# noise = np.std(np.abs(im_slice[0:8,10:20]))\n",
    "# signal = np.mean(np.abs(im_slice[20:40,0:10]))\n",
    "# snr=signal/noise\n",
    "# print('SNR')\n",
    "# snr"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f93965b-45c9-4e31-85a5-e011984bba85",
   "metadata": {},
   "source": [
    "### Background Masked"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "d8ad5671-5d8a-4190-b7b6-130f380b995f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'trained_models/AUTOMAP/train_25x75_fmri-hcp_r2_real_noise_bkgd_20240318/checkpoint/model.keras'"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_real_dir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "736ff93b-cd11-4f9f-a252-a31fcd8a88f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "volumeGTscaledmasked = np.copy(volumeGTscaled[:,:,1:])\n",
    "volumeGTscaledmasked[mask_metrics[:,:,1:]]=0\n",
    "vol1 = np.copy(volifftrecon[:,:,1:])#np.ma.masked_array(volifftrecon, mask = mask_metrics)\n",
    "vol1[mask_metrics[:,:,1:]]=0\n",
    "vol1title = 'Zero-filled IFFT\\nInp. SNR='+str(noise_val)+', R = '+str(R)\n",
    "vol2 = np.copy(volcsrecon[:,:,1:])\n",
    "vol2[mask_metrics[:,:,1:]]=0\n",
    "vol2title = 'Compressed Sensing\\nInp. SNR='+str(noise_val)+', R = '+str(R)\n",
    "vol3 = np.copy(volamaprecon[:,:,1:])\n",
    "vol3[mask_metrics[:,:,1:]]=0\n",
    "vol3title = 'AUTOMAP\\nInp. SNR='+str(noise_val)+', R = '+str(R)\n",
    "vol4 = np.copy(volunrollrecon[:,:,1:])\n",
    "vol4[mask_metrics[:,:,1:]]=0\n",
    "vol4title = 'Unrolling\\nInp. SNR='+str(noise_val)+', R = '+str(R)\n",
    "\n",
    "vol1_low_R2 = vol1\n",
    "vol2_low_R2 = vol2\n",
    "vol3_low_R2 = vol3\n",
    "vol4_low_R2 = vol4\n",
    "\n",
    "#slc=10\n",
    "#fig_comp = display_fns.compareFourPlot(volumeGTscaledmasked,vol1,vol1title,vol2,vol2title,vol3,vol3title,vol4,vol4title,slc, scale_err=3)\n",
    "#fig_comp.savefig('export_figs/fastmri_imgs_masked_snr'+str(noise_val)+'dB.pdf',bbox_inches='tight',dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92b09322-96ee-4454-aec3-e522e9b544d8",
   "metadata": {},
   "source": [
    "## High SNR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "6f6f43b5-3ad3-4e8e-907e-d7163156a28a",
   "metadata": {},
   "outputs": [],
   "source": [
    "#test_no=125\n",
    "kspTestScaled = ksp_fs_sc_all[:,:,:,:,test_no]\n",
    "volumeGTscaled = img_gt_sc_all[:,:,:,test_no]\n",
    "mask_metrics = np.ma.getmask(np.ma.masked_less(abs(volumeGTscaled),0.05)) # a mask in the image domain for calculating image metrics (not for masking kspace)\n",
    "\n",
    "l=-1\n",
    "noise_val = noise_dB[l]\n",
    "lamda_val = lamda_opt[l]\n",
    "iter_val = iter_opt[l]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "d88729a5-950a-4957-a6e9-dcdab66dbb2e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "80"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "noise_val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "eacadaf5-1fa8-430e-8c52-7d19ac04fd09",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1, 75, 64, 25)\n"
     ]
    }
   ],
   "source": [
    "kspInput = awgn(kspTestScaled,noise_val,L=1)\n",
    "uskspace = applyMask(kspInput,mask)\n",
    "\n",
    "volcsrecon = ulfl1recon(uskspace,mask,lamda_val,iter_val,mps)\n",
    "    \n",
    "volifftrecon_mc = sp.ifft(uskspace,axes=[1,2,3])\n",
    "volifftrecon = coil_combine(volifftrecon_mc,mps)# np.sqrt(np.sum(np.abs(volifftrecon_mc)**2,axis=0)).reshape(imSize)\n",
    "\n",
    "amap_scale = 1\n",
    "unroll_scale = 1#0.25#1/8\n",
    "volunrollrecon = unrolling_fns.unrollingRecon(uskspace*unroll_scale,mask,unroll_model_pth)/unroll_scale\n",
    "#volunrollrecon = np.roll(volunrollrecon,(-1,-1),axis=(0,2))\n",
    "volamaprecon, volamaprecon_mc = automap_fns.automapRecon(uskspace*amap_scale,mps,model_real_dir,model_imag_dir)\n",
    "volamaprecon = volamaprecon/amap_scale\n",
    "volamaprecon_mc = volamaprecon_mc/amap_scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "7911cd70-92fb-44a5-9bcb-45a288af797d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculating SNR of input kspace (fully sampled)\n",
    "# fs_recon = sp.ifft(kspInput,axes=[1,2,3])\n",
    "# im_slice = fs_recon[0,:,32,:]\n",
    "# fig=plt.figure()\n",
    "# plt.imshow((abs(im_slice)))\n",
    "# noise = np.std(np.abs(im_slice[0:8,10:20]))\n",
    "# signal = np.mean(np.abs(im_slice[20:40,0:10]))\n",
    "# snr=signal/noise\n",
    "# print('SNR')\n",
    "# snr"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "063a3189-21af-4290-ba73-43ac8803f0ef",
   "metadata": {},
   "source": [
    "### Background Masked"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "8d480b98-211a-4485-8e9b-a95cc3bbac81",
   "metadata": {},
   "outputs": [],
   "source": [
    "volumeGTscaledmasked = np.copy(volumeGTscaled[:,:,1:])\n",
    "volumeGTscaledmasked[mask_metrics[:,:,1:]]=0\n",
    "vol1 = np.copy(volifftrecon[:,:,1:])#np.ma.masked_array(volifftrecon, mask = mask_metrics)\n",
    "vol1[mask_metrics[:,:,1:]]=0\n",
    "vol1title = 'Zero-filled IFFT\\nInp. SNR='+str(noise_val)+', R = '+str(R)\n",
    "vol2 = np.copy(volcsrecon[:,:,1:])\n",
    "vol2[mask_metrics[:,:,1:]]=0\n",
    "vol2title = 'Compressed Sensing\\nInp. SNR='+str(noise_val)+', R = '+str(R)\n",
    "vol3 = np.copy(volamaprecon[:,:,1:])\n",
    "vol3[mask_metrics[:,:,1:]]=0\n",
    "vol3title = 'AUTOMAP\\nInp. SNR='+str(noise_val)+', R = '+str(R)\n",
    "vol4 = np.copy(volunrollrecon[:,:,1:])\n",
    "vol4[mask_metrics[:,:,1:]]=0\n",
    "vol4title = 'Unrolling\\nInp. SNR='+str(noise_val)+', R = '+str(R)\n",
    "\n",
    "vol1_high_R2 = vol1\n",
    "vol2_high_R2 = vol2\n",
    "vol3_high_R2 = vol3\n",
    "vol4_high_R2 = vol4\n",
    "\n",
    "#slc=10\n",
    "#fig_comp = display_fns.compareFourPlot(volumeGTscaledmasked,vol1,vol1title,vol2,vol2title,vol3,vol3title,vol4,vol4title,slc, scale_err=3)\n",
    "\n",
    "#fig_comp.savefig('export_figs/fastmri_imgs_masked_snr'+str(noise_val)+'dB.pdf',bbox_inches='tight',dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "bbaec820-99a7-4c37-a904-6c674aaa80ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "GT_R2 = volumeGTscaledmasked"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "c910358c-0403-4c20-81f0-222306101b7d",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "#loss_fn = 'perceptual'\n",
    "loss_fn = 'mse'\n",
    "# defining mask\n",
    "R=4\n",
    "calib_size=10\n",
    "in_seed = 0\n",
    "imSize = img_gt_sc_all[:,:,:,0].shape\n",
    "mask = poissonDiscSigpy(imSize,R, in_seed,calib_size)\n",
    "mask_R4 = mask\n",
    "\n",
    "mps = np.ones(ksp_fs_sc_all[:,:,:,:,0].shape,dtype=complex) #single coil"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "31a1e3a4-f8d4-4544-95df-826f6efbfaf7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# if R==2:\n",
    "#     if loss_fn=='perceptual':\n",
    "#         unroll_model_pth = 'trained_models/Unrolling/DisCorNet_perceptual_R2.pth'\n",
    "#         model_real_dir = 'keras_amap/experiments/train_25x75_fmri-hcp_r2_real_noise_bkgd_20240523_percepmse/checkpoint/model.keras' #evenly weight perceptual and mse loss\n",
    "#         model_imag_dir = 'keras_amap/experiments/train_25x75_fmri-hcp_r2_imag_noise_bkgd_20240523_percepmse/checkpoint/model.keras'\n",
    "#     elif loss_fn == 'mse':\n",
    "#         unroll_model_pth = 'trained_models/Unrolling/DisCorNet_mse_R2.pth'\n",
    "#         model_real_dir = 'keras_amap/experiments/train_25x75_fmri-hcp_r2_real_noise_bkgd_20240318/checkpoint/model.keras' #fast mri model\n",
    "#         model_imag_dir = 'keras_amap/experiments/train_25x75_fmri-hcp_r2_imag_noise_bkgd_20240318/checkpoint/model.keras'\n",
    "#     #model_real_dir = 'keras_amap/experiments/train_25x75_fmri-hcp_r2_real_noise_bkgd_20230911/checkpoint/' #fast mri model\n",
    "#     #model_imag_dir = 'keras_amap/experiments/train_25x75_fmri-hcp_r2_imag_noise_bkgd_20230911/checkpoint/'\n",
    "\n",
    "    \n",
    "# elif R==4:\n",
    "#     if loss_fn == 'perceptual':\n",
    "#         unroll_model_pth = 'trained_models/Unrolling/DisCorNet_perceptual_R4.pth'\n",
    "#     elif loss_fn == 'mse':\n",
    "#         unroll_model_pth = 'trained_models/Unrolling/DisCorNet_mse_R4.pth'\n",
    "#     #AUTOMAP model\n",
    "#     #model_real_dir = 'keras_amap/experiments/train_25x75_fmri-hcp_r4_real_noise_bkgd_20240304/checkpoint/' #fast mri model\n",
    "#     #model_imag_dir = 'keras_amap/experiments/train_25x75_fmri-hcp_r4_imag_noise_bkgd_20240304/checkpoint/'\n",
    "#     model_real_dir = 'keras_amap/experiments/train_25x75_fmri-hcp_r4_real_noise_bkgd_20240328/checkpoint/model.keras' #fast mri model\n",
    "#     model_imag_dir = 'keras_amap/experiments/train_25x75_fmri-hcp_r4_imag_noise_bkgd_20240328/checkpoint/model.keras'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "3e9c1074-1889-459a-a72c-cccc6b223e41",
   "metadata": {},
   "outputs": [],
   "source": [
    "if R==2:\n",
    "    unroll_model_pth = 'trained_models/Unrolling/DisCorNet_mse_R2.pth'\n",
    "    model_real_dir = 'trained_models/AUTOMAP/train_25x75_fmri-hcp_r2_real_noise_bkgd_20240318/checkpoint/model.keras' #fast mri model\n",
    "    model_imag_dir = 'trained_models/AUTOMAP/train_25x75_fmri-hcp_r2_imag_noise_bkgd_20240318/checkpoint/model.keras'\n",
    "    \n",
    "elif R==4:\n",
    "    unroll_model_pth = 'trained_models/Unrolling/DisCorNet_mse_R4.pth'\n",
    "    model_real_dir = 'trained_models/AUTOMAP/train_25x75_fmri-hcp_r4_real_noise_bkgd_20240328/checkpoint/model.keras' #fast mri model\n",
    "    model_imag_dir = 'trained_models/AUTOMAP/train_25x75_fmri-hcp_r4_imag_noise_bkgd_20240328/checkpoint/model.keras'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "eea0694a-7f96-4d15-9576-bb888cb6fd98",
   "metadata": {},
   "outputs": [],
   "source": [
    "#mdic = {\"mask\":mask}\n",
    "#savemat(\"mask_R2_PD.mat\",mdic)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "94a2fd6c-a346-4d59-a203-458942db32fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "lamda_opt_vals = np.zeros((noise_dB.size,n_ex_cs))\n",
    "iter_opt_vals = np.zeros((noise_dB.size,n_ex_cs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "5f52701d-0eed-4d35-be78-70ab1928ec57",
   "metadata": {},
   "outputs": [],
   "source": [
    "for j, noise_val in enumerate(noise_dB):\n",
    "    for k in range(n_ex_cs):\n",
    "        kspTestScaled = ksp_fs_sc_all[:,:,:,:,k]\n",
    "        volumeGTscaled = img_gt_sc_all[:,:,:,k]\n",
    "        \n",
    "        kspInput = awgn(kspTestScaled,noise_val,L=1)\n",
    "        uskspace = applyMask(kspInput,mask)\n",
    "\n",
    "        lamda_opt_vals[j,k] = find_lamda_mask(uskspace, volumeGTscaled, mps,calib_size, show_plot=False)\n",
    "        iter_opt_vals[j,k] = find_iter_mask(uskspace,volumeGTscaled, lamda_opt_vals[j,k], mps, show_plot=False)\n",
    "\n",
    "lamda_opt = lamda_opt_vals.mean(axis=1)\n",
    "iter_opt = iter_opt_vals.mean(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "5a44fa6f-af31-4405-bd15-b4fe36570e48",
   "metadata": {},
   "outputs": [],
   "source": [
    "nrmse_cs_vals = np.zeros((noise_dB.size,n_ex_test))\n",
    "nrmse_ifft_vals = np.zeros((noise_dB.size,n_ex_test))\n",
    "nrmse_amap_vals = np.zeros((noise_dB.size,n_ex_test))\n",
    "nrmse_unroll_vals = np.zeros((noise_dB.size,n_ex_test))\n",
    "\n",
    "ssim_cs_vals = np.zeros((noise_dB.size,n_ex_test))\n",
    "ssim_ifft_vals = np.zeros((noise_dB.size,n_ex_test))\n",
    "ssim_amap_vals = np.zeros((noise_dB.size,n_ex_test))\n",
    "ssim_unroll_vals = np.zeros((noise_dB.size,n_ex_test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "6fee2302-1c08-494c-a2a5-86c1d204bb27",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1, 75, 64, 25)\n",
      "(1, 75, 64, 25)\n",
      "(1, 75, 64, 25)\n",
      "(1, 75, 64, 25)\n",
      "(1, 75, 64, 25)\n",
      "(1, 75, 64, 25)\n",
      "(1, 75, 64, 25)\n",
      "(1, 75, 64, 25)\n",
      "(1, 75, 64, 25)\n"
     ]
    }
   ],
   "source": [
    "for j, noise_val in enumerate(noise_dB):\n",
    "    for k in range(n_ex_test):\n",
    "        kspTestScaled = ksp_fs_sc_all[:,:,:,:,k+n_ex_cs]\n",
    "        volumeGTscaled = img_gt_sc_all[:,:,:,k+n_ex_cs]\n",
    "        \n",
    "        kspInput = awgn(kspTestScaled,noise_val,L=1)\n",
    "        uskspace = applyMask(kspInput,mask)\n",
    "\n",
    "    \n",
    "        volifftrecon_mc = sp.ifft(uskspace,axes=[1,2,3])\n",
    "        volifftrecon = coil_combine(volifftrecon_mc,mps)# np.sqrt(np.sum(np.abs(volifftrecon_mc)**2,axis=0)).reshape(imSize)\n",
    "\n",
    "        volcsrecon = ulfl1recon(uskspace,mask,lamda_opt[j],iter_opt[j],mps)\n",
    "\n",
    "        amap_scale = 1\n",
    "        unroll_scale = 1#1/8\n",
    "        volunrollrecon = unrolling_fns.unrollingRecon(uskspace*unroll_scale,mask,unroll_model_pth)/unroll_scale\n",
    "        #volunrollrecon = np.roll(volunrollrecon,(1,1),axis=(0,2))\n",
    "        volamaprecon, volamaprecon_mc = automap_fns.automapRecon(uskspace*amap_scale,mps,model_real_dir,model_imag_dir)\n",
    "        volamaprecon = volamaprecon/amap_scale\n",
    "        volamaprecon_mc = volamaprecon_mc/amap_scale\n",
    "\n",
    "        mask_metrics = np.ma.getmask(np.ma.masked_less(abs(volumeGTscaled),0.05)) # a mask in the image domain for calculating image metrics (not for masking kspace)\n",
    "\n",
    "        volumeGTscaled[mask_metrics]=0\n",
    "        volcsrecon[mask_metrics]=0\n",
    "        volifftrecon[mask_metrics]=0\n",
    "        volamaprecon[mask_metrics]=0\n",
    "        volunrollrecon[mask_metrics]=0\n",
    "\n",
    "        nrmse_cs_vals[j,k] = nrmse(np.abs(volcsrecon[:,:,1:]),np.abs(volumeGTscaled[:,:,1:]))\n",
    "        nrmse_ifft_vals[j,k] = nrmse(np.abs(volifftrecon[:,:,1:]),np.abs(volumeGTscaled[:,:,1:]))\n",
    "        nrmse_amap_vals[j,k] = nrmse(np.abs(volamaprecon[:,:,1:]),np.abs(volumeGTscaled[:,:,1:]))\n",
    "        nrmse_unroll_vals[j,k] = nrmse(np.abs(volunrollrecon[:,:,1:]),np.abs(volumeGTscaled[:,:,1:]))\n",
    "\n",
    "        ssim_cs_vals[j,k] = ssim(np.abs(volcsrecon[:,:,1:]),np.abs(volumeGTscaled[:,:,1:]))\n",
    "        ssim_ifft_vals[j,k] = ssim(np.abs(volifftrecon[:,:,1:]),np.abs(volumeGTscaled[:,:,1:]))\n",
    "        ssim_amap_vals[j,k] = ssim(np.abs(volamaprecon[:,:,1:]),np.abs(volumeGTscaled[:,:,1:]))\n",
    "        ssim_unroll_vals[j,k] = ssim(np.abs(volunrollrecon[:,:,1:]),np.abs(volumeGTscaled[:,:,1:]))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "ae0fcd8a-94e7-4809-aa33-c0ff72dc3c59",
   "metadata": {},
   "outputs": [],
   "source": [
    "nrmse_cs = nrmse_cs_vals.mean(axis=1)\n",
    "nrmse_ifft = nrmse_ifft_vals.mean(axis=1)\n",
    "nrmse_amap = nrmse_amap_vals.mean(axis=1)\n",
    "nrmse_unroll = nrmse_unroll_vals.mean(axis=1)\n",
    "\n",
    "nrmse_cs_sd = nrmse_cs_vals.std(axis=1)\n",
    "nrmse_ifft_sd = nrmse_ifft_vals.std(axis=1)\n",
    "nrmse_amap_sd = nrmse_amap_vals.std(axis=1)\n",
    "nrmse_unroll_sd = nrmse_unroll_vals.std(axis=1)\n",
    "\n",
    "ssim_cs = ssim_cs_vals.mean(axis=1)\n",
    "ssim_ifft = ssim_ifft_vals.mean(axis=1)\n",
    "ssim_amap = ssim_amap_vals.mean(axis=1)\n",
    "ssim_unroll = ssim_unroll_vals.mean(axis=1)\n",
    "\n",
    "ssim_cs_sd = ssim_cs_vals.std(axis=1)\n",
    "ssim_ifft_sd = ssim_ifft_vals.std(axis=1)\n",
    "ssim_amap_sd = ssim_amap_vals.std(axis=1)\n",
    "ssim_unroll_sd = ssim_unroll_vals.std(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "52d8cf35-6f6a-4e27-be4b-3ad32dfbb6af",
   "metadata": {},
   "outputs": [],
   "source": [
    "#ax2.clear()\n",
    "#ax4.clear()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "698b8574-0f19-41a5-8133-b04b5632bcdb",
   "metadata": {},
   "outputs": [
    {
     "ename": "FileNotFoundError",
     "evalue": "[Errno 2] No such file or directory: 'export_figs/fastmri_noise_curves_all.pdf'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mFileNotFoundError\u001b[0m                         Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[43], line 62\u001b[0m\n\u001b[1;32m     59\u001b[0m fig\n\u001b[1;32m     60\u001b[0m \u001b[38;5;66;03m#plt.show()\u001b[39;00m\n\u001b[0;32m---> 62\u001b[0m \u001b[43mfig\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43msavefig\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mexport_figs/fastmri_noise_curves_all.pdf\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43mbbox_inches\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mtight\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43mdpi\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m300\u001b[39;49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/repos/lowfieldrecon/LowFieldReconPrivate/.venvkerasandtf/lib/python3.10/site-packages/matplotlib/figure.py:3390\u001b[0m, in \u001b[0;36mFigure.savefig\u001b[0;34m(self, fname, transparent, **kwargs)\u001b[0m\n\u001b[1;32m   3388\u001b[0m     \u001b[38;5;28;01mfor\u001b[39;00m ax \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39maxes:\n\u001b[1;32m   3389\u001b[0m         _recursively_make_axes_transparent(stack, ax)\n\u001b[0;32m-> 3390\u001b[0m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcanvas\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mprint_figure\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfname\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/repos/lowfieldrecon/LowFieldReconPrivate/.venvkerasandtf/lib/python3.10/site-packages/matplotlib/backend_bases.py:2193\u001b[0m, in \u001b[0;36mFigureCanvasBase.print_figure\u001b[0;34m(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)\u001b[0m\n\u001b[1;32m   2189\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[1;32m   2190\u001b[0m     \u001b[38;5;66;03m# _get_renderer may change the figure dpi (as vector formats\u001b[39;00m\n\u001b[1;32m   2191\u001b[0m     \u001b[38;5;66;03m# force the figure dpi to 72), so we need to set it again here.\u001b[39;00m\n\u001b[1;32m   2192\u001b[0m     \u001b[38;5;28;01mwith\u001b[39;00m cbook\u001b[38;5;241m.\u001b[39m_setattr_cm(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mfigure, dpi\u001b[38;5;241m=\u001b[39mdpi):\n\u001b[0;32m-> 2193\u001b[0m         result \u001b[38;5;241m=\u001b[39m \u001b[43mprint_method\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m   2194\u001b[0m \u001b[43m            \u001b[49m\u001b[43mfilename\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   2195\u001b[0m \u001b[43m            \u001b[49m\u001b[43mfacecolor\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mfacecolor\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   2196\u001b[0m \u001b[43m            \u001b[49m\u001b[43medgecolor\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43medgecolor\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   2197\u001b[0m \u001b[43m            \u001b[49m\u001b[43morientation\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43morientation\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   2198\u001b[0m \u001b[43m            \u001b[49m\u001b[43mbbox_inches_restore\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m_bbox_inches_restore\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   2199\u001b[0m \u001b[43m            \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m   2200\u001b[0m \u001b[38;5;28;01mfinally\u001b[39;00m:\n\u001b[1;32m   2201\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m bbox_inches \u001b[38;5;129;01mand\u001b[39;00m restore_bbox:\n",
      "File \u001b[0;32m~/repos/lowfieldrecon/LowFieldReconPrivate/.venvkerasandtf/lib/python3.10/site-packages/matplotlib/backend_bases.py:2043\u001b[0m, in \u001b[0;36mFigureCanvasBase._switch_canvas_and_return_print_method.<locals>.<lambda>\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m   2039\u001b[0m     optional_kws \u001b[38;5;241m=\u001b[39m {  \u001b[38;5;66;03m# Passed by print_figure for other renderers.\u001b[39;00m\n\u001b[1;32m   2040\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mdpi\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mfacecolor\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124medgecolor\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124morientation\u001b[39m\u001b[38;5;124m\"\u001b[39m,\n\u001b[1;32m   2041\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mbbox_inches_restore\u001b[39m\u001b[38;5;124m\"\u001b[39m}\n\u001b[1;32m   2042\u001b[0m     skip \u001b[38;5;241m=\u001b[39m optional_kws \u001b[38;5;241m-\u001b[39m {\u001b[38;5;241m*\u001b[39minspect\u001b[38;5;241m.\u001b[39msignature(meth)\u001b[38;5;241m.\u001b[39mparameters}\n\u001b[0;32m-> 2043\u001b[0m     print_method \u001b[38;5;241m=\u001b[39m functools\u001b[38;5;241m.\u001b[39mwraps(meth)(\u001b[38;5;28;01mlambda\u001b[39;00m \u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs: \u001b[43mmeth\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m   2044\u001b[0m \u001b[43m        \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43m{\u001b[49m\u001b[43mk\u001b[49m\u001b[43m:\u001b[49m\u001b[43m \u001b[49m\u001b[43mv\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mfor\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mk\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mv\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;129;43;01min\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mkwargs\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mitems\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mif\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mk\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;129;43;01mnot\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[38;5;129;43;01min\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mskip\u001b[49m\u001b[43m}\u001b[49m\u001b[43m)\u001b[49m)\n\u001b[1;32m   2045\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:  \u001b[38;5;66;03m# Let third-parties do as they see fit.\u001b[39;00m\n\u001b[1;32m   2046\u001b[0m     print_method \u001b[38;5;241m=\u001b[39m meth\n",
      "File \u001b[0;32m~/repos/lowfieldrecon/LowFieldReconPrivate/.venvkerasandtf/lib/python3.10/site-packages/matplotlib/backends/backend_pdf.py:2800\u001b[0m, in \u001b[0;36mFigureCanvasPdf.print_pdf\u001b[0;34m(self, filename, bbox_inches_restore, metadata)\u001b[0m\n\u001b[1;32m   2798\u001b[0m     file \u001b[38;5;241m=\u001b[39m filename\u001b[38;5;241m.\u001b[39m_ensure_file()\n\u001b[1;32m   2799\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m-> 2800\u001b[0m     file \u001b[38;5;241m=\u001b[39m \u001b[43mPdfFile\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfilename\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmetadata\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mmetadata\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m   2801\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[1;32m   2802\u001b[0m     file\u001b[38;5;241m.\u001b[39mnewPage(width, height)\n",
      "File \u001b[0;32m~/repos/lowfieldrecon/LowFieldReconPrivate/.venvkerasandtf/lib/python3.10/site-packages/matplotlib/backends/backend_pdf.py:688\u001b[0m, in \u001b[0;36mPdfFile.__init__\u001b[0;34m(self, filename, metadata)\u001b[0m\n\u001b[1;32m    686\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39moriginal_file_like \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mNone\u001b[39;00m\n\u001b[1;32m    687\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mtell_base \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m0\u001b[39m\n\u001b[0;32m--> 688\u001b[0m fh, opened \u001b[38;5;241m=\u001b[39m \u001b[43mcbook\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mto_filehandle\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfilename\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mwb\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mreturn_opened\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mTrue\u001b[39;49;00m\u001b[43m)\u001b[49m\n\u001b[1;32m    689\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m opened:\n\u001b[1;32m    690\u001b[0m     \u001b[38;5;28;01mtry\u001b[39;00m:\n",
      "File \u001b[0;32m~/repos/lowfieldrecon/LowFieldReconPrivate/.venvkerasandtf/lib/python3.10/site-packages/matplotlib/cbook.py:483\u001b[0m, in \u001b[0;36mto_filehandle\u001b[0;34m(fname, flag, return_opened, encoding)\u001b[0m\n\u001b[1;32m    481\u001b[0m         fh \u001b[38;5;241m=\u001b[39m bz2\u001b[38;5;241m.\u001b[39mBZ2File(fname, flag)\n\u001b[1;32m    482\u001b[0m     \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m--> 483\u001b[0m         fh \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mopen\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mfname\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mflag\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mencoding\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mencoding\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    484\u001b[0m     opened \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mTrue\u001b[39;00m\n\u001b[1;32m    485\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m \u001b[38;5;28mhasattr\u001b[39m(fname, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mseek\u001b[39m\u001b[38;5;124m'\u001b[39m):\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: 'export_figs/fastmri_noise_curves_all.pdf'"
     ]
    }
   ],
   "source": [
    "line1 = ax2.plot(noise_dB_zero,nrmse_ifft, '-')\n",
    "line2 = ax2.plot(noise_dB_zero,nrmse_cs, '-')\n",
    "line3 = ax2.plot(noise_dB_zero,nrmse_amap, '-')\n",
    "line4 = ax2.plot(noise_dB_zero,nrmse_unroll, '-')\n",
    "#ax2.set(title='Comparison of Reconstruction Methods, R='+str(R),\n",
    "ax2.set(ylabel='NRMSE',\n",
    "xlabel='Input SNR (dB)')\n",
    "#ax2.set_xscale('log')\n",
    "ax2.legend(['Zero-filled', 'CS', 'AUTOMAP', 'Unrolled'])\n",
    "#ax2.legend([line1, line2, line3], ['label1', 'label2', 'label3'])\n",
    "#plt.ylim(0,0.085)\n",
    "#plt.xlim(55,90)\n",
    "fill1 = ax2.fill_between(noise_dB_zero,nrmse_ifft-nrmse_ifft_sd,nrmse_ifft+nrmse_ifft_sd,alpha=0.25)\n",
    "fill2 = ax2.fill_between(noise_dB_zero,nrmse_cs-nrmse_cs_sd,nrmse_cs+nrmse_cs_sd,alpha=0.25)\n",
    "fill3 = ax2.fill_between(noise_dB_zero,nrmse_amap-nrmse_amap_sd,nrmse_amap+nrmse_amap_sd,alpha=0.25)\n",
    "fill4 = ax2.fill_between(noise_dB_zero,nrmse_unroll-nrmse_unroll_sd,nrmse_unroll+nrmse_unroll_sd,alpha=0.25)\n",
    "\n",
    "step=10\n",
    "ax2.set_xticks(np.arange(noise_dB_zero[0],noise_dB_zero[-1]+step,step))\n",
    "ax2.set_yticks([0,0.05,0.1])\n",
    "\n",
    "line1 = ax4.plot(noise_dB_zero,ssim_ifft, '-')\n",
    "line2 = ax4.plot(noise_dB_zero,ssim_cs, '-')\n",
    "line3 = ax4.plot(noise_dB_zero,ssim_amap, '-')\n",
    "line4 = ax4.plot(noise_dB_zero,ssim_unroll, '-')\n",
    "#ax4.set(title='Comparison of Reconstruction Methods, R='+str(R),\n",
    "ax4.set(ylabel='SSIM',\n",
    "xlabel='Input SNR (dB)')\n",
    "#ax2.set_xscale('log')\n",
    "ax4.legend(['Zero-filled', 'CS', 'AUTOMAP', 'Unrolled AI'])\n",
    "#ax2.legend([line1, line2, line3], ['label1', 'label2', 'label3'])\n",
    "#plt.ylim(0,0.085)\n",
    "#plt.xlim(55,90)\n",
    "\n",
    "fill1 = ax4.fill_between(noise_dB_zero,ssim_ifft-ssim_ifft_sd,ssim_ifft+ssim_ifft_sd,alpha=0.25)\n",
    "fill2 = ax4.fill_between(noise_dB_zero,ssim_cs-ssim_cs_sd,ssim_cs+ssim_cs_sd,alpha=0.25)\n",
    "fill3 = ax4.fill_between(noise_dB_zero,ssim_amap-ssim_amap_sd,ssim_amap+ssim_amap_sd,alpha=0.25)\n",
    "fill4 = ax4.fill_between(noise_dB_zero,ssim_unroll-ssim_unroll_sd,ssim_unroll+ssim_unroll_sd,alpha=0.25)\n",
    "\n",
    "step=10\n",
    "ax4.set_xticks(np.arange(noise_dB_zero[0],noise_dB_zero[-1]+step,step))\n",
    "ax4.set_yticks([0.8,0.9,1])\n",
    "\n",
    "ax1.set_title('(a) R = 2 Subsampling',x=0,y=1.04)\n",
    "ax2.set_title('(b) R = 4 Subsampling',x=0,y=1.04)\n",
    "\n",
    "for ax in [ax1,ax2]:\n",
    "    ax.vlines(23,0,0.1,color='tab:gray',linestyle='--')\n",
    "    ax.set_ylim(0,0.1)\n",
    "    #ax.set_xlim(10,50)\n",
    "\n",
    "for ax in [ax3,ax4]:\n",
    "    ax.vlines(23,0,1,color='tab:gray',linestyle='--')\n",
    "    ax.set_ylim(0.7,1)\n",
    "    #ax.set_xlim(10,50)\n",
    "\n",
    "fig.tight_layout()\n",
    "\n",
    "fig\n",
    "#plt.show()\n",
    "\n",
    "fig.savefig('export_figs/fastmri_noise_curves_all.pdf',bbox_inches='tight',dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a146114-7e9a-4eaa-8f5b-123e85f1b208",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Low SNR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93e1176e-f943-4ae4-b657-74e5b3cb33ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "test_no = 125\n",
    "kspTestScaled = ksp_fs_sc_all[:,:,:,:,test_no]\n",
    "volumeGTscaled = img_gt_sc_all[:,:,:,test_no]\n",
    "mask_metrics = np.ma.getmask(np.ma.masked_less(abs(volumeGTscaled),0.05)) # a mask in the image domain for calculating image metrics (not for masking kspace)\n",
    "\n",
    "l = 0\n",
    "noise_val = noise_dB[l]\n",
    "lamda_val = lamda_opt[l]\n",
    "iter_val = iter_opt[l]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3df5538e-a755-4cd8-890a-490d804f5648",
   "metadata": {},
   "outputs": [],
   "source": [
    "noise_val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83949ba8-94e9-4ec0-9e30-dc08da57a195",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "kspInput = awgn(kspTestScaled,noise_val,L=1)\n",
    "uskspace = applyMask(kspInput,mask)\n",
    "\n",
    "\n",
    "volcsrecon = ulfl1recon(uskspace,mask,lamda_val,iter_val,mps)\n",
    "    \n",
    "volifftrecon_mc = sp.ifft(uskspace,axes=[1,2,3])\n",
    "volifftrecon = coil_combine(volifftrecon_mc,mps)# np.sqrt(np.sum(np.abs(volifftrecon_mc)**2,axis=0)).reshape(imSize)\n",
    "\n",
    "amap_scale = 1\n",
    "unroll_scale = 1#0.25#1/8\n",
    "volunrollrecon = unrolling_fns.unrollingRecon(uskspace*unroll_scale,mask,unroll_model_pth)/unroll_scale\n",
    "#volunrollrecon = np.roll(volunrollrecon,(-1,-1),axis=(0,2))\n",
    "volamaprecon, volamaprecon_mc = automap_fns.automapRecon(uskspace*amap_scale,mps,model_real_dir,model_imag_dir)\n",
    "volamaprecon = volamaprecon/amap_scale\n",
    "volamaprecon_mc = volamaprecon_mc/amap_scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "413a4b82-096d-456b-aefb-bee3d65b0690",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # calculating SNR of input kspace (fully sampled)\n",
    "# fs_recon = sp.ifft(kspInput,axes=[1,2,3])\n",
    "# im_slice = fs_recon[0,:,32,:]\n",
    "# fig=plt.figure()\n",
    "# plt.imshow((abs(im_slice)))\n",
    "# noise = np.std(np.abs(im_slice[0:8,10:20]))\n",
    "# signal = np.mean(np.abs(im_slice[20:40,0:10]))\n",
    "# snr=signal/noise\n",
    "# print('SNR')\n",
    "# snr"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15523123-e3ac-456c-9d97-39bdc0526b87",
   "metadata": {},
   "source": [
    "### Background Masked"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0579a994-8e54-4caf-a452-712d52887a0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_real_dir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "951a5cfb-ac23-413c-843c-485f63877f8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "volumeGTscaledmasked = np.copy(volumeGTscaled[:,:,1:])\n",
    "volumeGTscaledmasked[mask_metrics[:,:,1:]]=0\n",
    "vol1 = np.copy(volifftrecon[:,:,1:])#np.ma.masked_array(volifftrecon, mask = mask_metrics)\n",
    "vol1[mask_metrics[:,:,1:]]=0\n",
    "vol1title = 'Zero-filled IFFT\\nInp. SNR='+str(noise_val)+', R = '+str(R)\n",
    "vol2 = np.copy(volcsrecon[:,:,1:])\n",
    "vol2[mask_metrics[:,:,1:]]=0\n",
    "vol2title = 'Compressed Sensing\\nInp. SNR='+str(noise_val)+', R = '+str(R)\n",
    "vol3 = np.copy(volamaprecon[:,:,1:])\n",
    "vol3[mask_metrics[:,:,1:]]=0\n",
    "vol3title = 'AUTOMAP\\nInp. SNR='+str(noise_val)+', R = '+str(R)\n",
    "vol4 = np.copy(volunrollrecon[:,:,1:])\n",
    "vol4[mask_metrics[:,:,1:]]=0\n",
    "vol4title = 'Unrolling\\nInp. SNR='+str(noise_val)+', R = '+str(R)\n",
    "\n",
    "vol1_low_R4 = vol1\n",
    "vol2_low_R4 = vol2\n",
    "vol3_low_R4 = vol3\n",
    "vol4_low_R4 = vol4\n",
    "\n",
    "#slc=10\n",
    "#fig_comp = display_fns.compareFourPlot(volumeGTscaledmasked,vol1,vol1title,vol2,vol2title,vol3,vol3title,vol4,vol4title,slc, scale_err=3)\n",
    "#fig_comp.savefig('export_figs/fastmri_imgs_masked_snr'+str(noise_val)+'dB.pdf',bbox_inches='tight',dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e981321c-2d83-4c0f-95d7-c0e08135b49f",
   "metadata": {},
   "source": [
    "## High SNR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06ed8a1d-c488-4aaa-903d-4746e1d7858c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#test_no=125\n",
    "kspTestScaled = ksp_fs_sc_all[:,:,:,:,test_no]\n",
    "volumeGTscaled = img_gt_sc_all[:,:,:,test_no]\n",
    "mask_metrics = np.ma.getmask(np.ma.masked_less(abs(volumeGTscaled),0.05)) # a mask in the image domain for calculating image metrics (not for masking kspace)\n",
    "\n",
    "l=-1\n",
    "noise_val = noise_dB[l]\n",
    "lamda_val = lamda_opt[l]\n",
    "iter_val = iter_opt[l]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "364db9e8-693f-4055-a748-d78cee421c33",
   "metadata": {},
   "outputs": [],
   "source": [
    "noise_val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd865459-851a-491e-8242-bd0bacf8f3e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "kspInput = awgn(kspTestScaled,noise_val,L=1)\n",
    "uskspace = applyMask(kspInput,mask)\n",
    "\n",
    "volcsrecon = ulfl1recon(uskspace,mask,lamda_val,iter_val,mps)\n",
    "    \n",
    "volifftrecon_mc = sp.ifft(uskspace,axes=[1,2,3])\n",
    "volifftrecon = coil_combine(volifftrecon_mc,mps)# np.sqrt(np.sum(np.abs(volifftrecon_mc)**2,axis=0)).reshape(imSize)\n",
    "\n",
    "amap_scale = 1\n",
    "unroll_scale = 1#0.25#1/8\n",
    "volunrollrecon = unrolling_fns.unrollingRecon(uskspace*unroll_scale,mask,unroll_model_pth)/unroll_scale\n",
    "#volunrollrecon = np.roll(volunrollrecon,(-1,-1),axis=(0,2))\n",
    "volamaprecon, volamaprecon_mc = automap_fns.automapRecon(uskspace*amap_scale,mps,model_real_dir,model_imag_dir)\n",
    "volamaprecon = volamaprecon/amap_scale\n",
    "volamaprecon_mc = volamaprecon_mc/amap_scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0fcc6a98-0a25-498f-a4e1-634f0fbd9dc1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculating SNR of input kspace (fully sampled)\n",
    "# fs_recon = sp.ifft(kspInput,axes=[1,2,3])\n",
    "# im_slice = fs_recon[0,:,32,:]\n",
    "# fig=plt.figure()\n",
    "# plt.imshow((abs(im_slice)))\n",
    "# noise = np.std(np.abs(im_slice[0:8,10:20]))\n",
    "# signal = np.mean(np.abs(im_slice[20:40,0:10]))\n",
    "# snr=signal/noise\n",
    "# print('SNR')\n",
    "# snr"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34606b09-0b5a-41d0-8415-90ada9e130af",
   "metadata": {},
   "source": [
    "### Background Masked"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "abe82658-e504-4ab9-92c7-dc63247c6ba4",
   "metadata": {},
   "outputs": [],
   "source": [
    "volumeGTscaledmasked = np.copy(volumeGTscaled[:,:,1:])\n",
    "volumeGTscaledmasked[mask_metrics[:,:,1:]]=0\n",
    "vol1 = np.copy(volifftrecon[:,:,1:])#np.ma.masked_array(volifftrecon, mask = mask_metrics)\n",
    "vol1[mask_metrics[:,:,1:]]=0\n",
    "vol1title = 'Zero-filled IFFT\\nInp. SNR='+str(noise_val)+', R = '+str(R)\n",
    "vol2 = np.copy(volcsrecon[:,:,1:])\n",
    "vol2[mask_metrics[:,:,1:]]=0\n",
    "vol2title = 'Compressed Sensing\\nInp. SNR='+str(noise_val)+', R = '+str(R)\n",
    "vol3 = np.copy(volamaprecon[:,:,1:])\n",
    "vol3[mask_metrics[:,:,1:]]=0\n",
    "vol3title = 'AUTOMAP\\nInp. SNR='+str(noise_val)+', R = '+str(R)\n",
    "vol4 = np.copy(volunrollrecon[:,:,1:])\n",
    "vol4[mask_metrics[:,:,1:]]=0\n",
    "vol4title = 'Unrolling\\nInp. SNR='+str(noise_val)+', R = '+str(R)\n",
    "\n",
    "vol1_high_R4 = vol1\n",
    "vol2_high_R4 = vol2\n",
    "vol3_high_R4 = vol3\n",
    "vol4_high_R4 = vol4\n",
    "\n",
    "#slc=10\n",
    "#fig_comp = display_fns.compareFourPlot(volumeGTscaledmasked,vol1,vol1title,vol2,vol2title,vol3,vol3title,vol4,vol4title,slc, scale_err=3)\n",
    "\n",
    "#fig_comp.savefig('export_figs/fastmri_imgs_masked_snr'+str(noise_val)+'dB.pdf',bbox_inches='tight',dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c7980f3-a284-484b-8371-d2aa3ac2d2ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "GT_R4 = volumeGTscaledmasked"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f321aa4-3e6e-4e61-a397-13ac76d93766",
   "metadata": {},
   "source": [
    "## Comparing images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a7b2124-a633-4940-9ac0-532eab9a89fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "slc=10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e850bdb-0d08-4e37-9ef9-9b46a4cb1c2d",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    " vol1title = 'Zero-filled'\n",
    "# vol2 = np.copy(volcsrecon[:,:,1:])\n",
    " #vol2[mask_metrics[:,:,1:]]=0\n",
    " vol2title = 'Comp. Sensing'\n",
    "# vol3 = np.copy(volamaprecon[:,:,1:])\n",
    "# #vol3[mask_metrics[:,:,1:]]=0\n",
    " vol3title = 'AUTOMAP'\n",
    "# vol4 = np.copy(volunrollrecon[:,:,1:])\n",
    "# #vol4[mask_metrics[:,:,1:]]=0\n",
    " vol4title = 'Unrolled AI'\n",
    "\n",
    "font_lw = 500\n",
    "\n",
    "#    GT_R2 = np.swapaxes(GT_R2,0,1)\n",
    "#    vol1 = np.swapaxes(vol1,0,1)\n",
    "#    vol2 = np.swapaxes(vol2,0,1)\n",
    "#    vol3 = np.swapaxes(vol3,0,1)\n",
    "#    vol4 = np.swapaxes(vol4,0,1)\n",
    "\n",
    "# RO 64, 2.5 mm\n",
    "# PE1 75, 3.5mm\n",
    "# PE2 15, 8 mm\n",
    "asp_ratio = (GT_R2.shape[1]*3.5)/(GT_R2.shape[0]*2.5)\n",
    "\n",
    "#plt.style.use('default')\n",
    "plt.style.use('dark_background')\n",
    "\n",
    "figimgs, axs = plt.subplots(4, 5,figsize = (8.5,5))\n",
    "\n",
    "figimgs.subplots_adjust(hspace=1.4)\n",
    "figimgs.subplots_adjust(wspace=-0.6)\n",
    "\n",
    "max_GT_R2 = np.amax(np.abs(GT_R2[:,:,slc]))\n",
    "\n",
    "im = axs[0, 0].imshow(np.abs(GT_R2[:,:,slc].T), cmap= 'gray', vmin = 0, vmax = 1.05*max_GT_R2)\n",
    "\n",
    "axs[0, 0].set_xticks([])\n",
    "axs[0, 0].set_yticks([])\n",
    "axs[0, 0].set_aspect(asp_ratio)\n",
    "axs[0, 0].set_xlabel('PE1')\n",
    "axs[0, 0].set_ylabel('RO')\n",
    "im.set_clim(0, (1.05*max_GT_R2))\n",
    "\n",
    "im = axs[0, 1].imshow(np.abs(vol1_high_R2[:,:,slc].T), cmap = 'gray', vmin = 0,vmax = 1.05*max_GT_R2)\n",
    "axs[0, 1].set_title(vol1title,fontsize=10, fontweight=font_lw)\n",
    "axs[0, 1].set_xticks([])\n",
    "axs[0, 1].set_yticks([])\n",
    "axs[0, 1].set_aspect(asp_ratio)\n",
    "im.set_clim(0, (1.05*max_GT_R2))\n",
    "axs[0, 1].set_xlabel('SSIM = '+str(round(ssim(np.abs(vol1_high_R2[:,:,:]),np.abs(GT_R2[:,:,:])),2))\n",
    "                     +'\\n NRMSE = '+str(round(nrmse(np.abs(vol1_high_R2[:,:,:]),np.abs(GT_R2[:,:,:])),3)),fontsize=7)\n",
    "\n",
    "axs[0, 1].set_ylabel('SNR '+str(noise_dB_zero[-1])+' dB',fontsize=10, fontweight=font_lw)\n",
    "\n",
    "im = axs[1, 1].imshow(np.abs(vol1_low_R2[:,:,slc].T), cmap = 'gray', vmin = 0, vmax = 1.05*max_GT_R2)\n",
    "#axs[1, 1].set_title(vol1title + ' Diff')\n",
    "axs[1, 1].set_xticks([])\n",
    "axs[1, 1].set_yticks([])\n",
    "axs[1, 1].set_aspect(asp_ratio)\n",
    "im.set_clim(0, (1.05*max_GT_R2))\n",
    "axs[1, 1].set_xlabel('SSIM ='+str(round(ssim(np.abs(vol1_low_R2[:,:,:]),np.abs(GT_R2[:,:,:])),2))\n",
    "                     +'\\n NRMSE ='+str(round(nrmse(np.abs(vol1_low_R2[:,:,:]),np.abs(GT_R2[:,:,:])),3)),fontsize=7)\n",
    "axs[1, 1].set_ylabel('SNR '+str(noise_dB_zero[0])+' dB',fontsize=10, fontweight=font_lw)\n",
    "\n",
    "\n",
    "im = axs[0, 2].imshow(np.abs(vol2_high_R2[:,:,slc].T), cmap = 'gray', vmin = 0, vmax = 1.05*max_GT_R2)\n",
    "axs[0, 2].set_title(vol2title, fontsize=10, fontweight=font_lw)\n",
    "axs[0, 2].set_xticks([])\n",
    "axs[0, 2].set_yticks([])\n",
    "axs[0, 2].set_aspect(asp_ratio)\n",
    "im.set_clim(0, (1.05*max_GT_R2))\n",
    "axs[0, 2].set_xlabel('SSIM = '+str(round(ssim(np.abs(vol2_high_R2[:,:,:]),np.abs(GT_R2[:,:,:])),2))\n",
    "                     +'\\n NRMSE = '+str(round(nrmse(np.abs(vol2_high_R2[:,:,:]),np.abs(GT_R2[:,:,:])),3)),fontsize=7)\n",
    "\n",
    "im = axs[1, 2].imshow(np.abs(vol2_low_R2[:,:,slc].T), cmap = 'gray', vmin = 0, vmax = 1.05*max_GT_R2)\n",
    "#axs[1, 2].set_title(vol2title + ' Diff')\n",
    "axs[1, 2].set_xticks([])\n",
    "axs[1, 2].set_yticks([])\n",
    "axs[1, 2].set_aspect(asp_ratio)\n",
    "im.set_clim(0, (1.05*max_GT_R2))\n",
    "axs[1, 2].set_xlabel('SSIM = '+str(round(ssim(np.abs(vol2_low_R2[:,:,:]),np.abs(GT_R2[:,:,:])),2))\n",
    "                     +'\\n NRMSE = '+str(round(nrmse(np.abs(vol2_low_R2[:,:,:]),np.abs(GT_R2[:,:,:])),3)),fontsize=7)\n",
    "\n",
    "im = axs[0, 3].imshow(np.abs(vol3_high_R2[:,:,slc].T), cmap = 'gray', vmin = 0, vmax = 1.05*max_GT_R2)\n",
    "axs[0, 3].set_title(vol3title, fontsize=10, fontweight=font_lw)\n",
    "axs[0, 3].set_xticks([])\n",
    "axs[0, 3].set_yticks([])\n",
    "axs[0, 3].set_aspect(asp_ratio)\n",
    "im.set_clim(0, (1.05*max_GT_R2))\n",
    "axs[0, 3].set_xlabel('SSIM = '+str(round(ssim(np.abs(vol3_high_R2[:,:,:]),np.abs(GT_R2[:,:,:])),2))\n",
    "                     +'\\n NRMSE = '+str(round(nrmse(np.abs(vol3_high_R2[:,:,:]),np.abs(GT_R2[:,:,:])),3)),fontsize=7)\n",
    "\n",
    "im = axs[1, 3].imshow(np.abs(vol3_low_R2[:,:,slc].T), cmap = 'gray', vmin = 0, vmax = 1.05*max_GT_R2)\n",
    "#axs[1, 3].set_title(vol3title + ' Diff')\n",
    "axs[1, 3].set_xticks([])\n",
    "axs[1, 3].set_yticks([])\n",
    "axs[1, 3].set_aspect(asp_ratio)\n",
    "im.set_clim(0, (1.05*max_GT_R2))\n",
    "axs[1, 3].set_xlabel('SSIM = '+str(round(ssim(np.abs(vol3_low_R2[:,:,:]),np.abs(GT_R2[:,:,:])),2))\n",
    "                     +'\\n NRMSE = '+str(round(nrmse(np.abs(vol3_low_R2[:,:,:]),np.abs(GT_R2[:,:,:])),3)),fontsize=7)\n",
    "\n",
    "im = axs[0, 4].imshow(np.abs(vol4_high_R2[:,:,slc].T), cmap = 'gray', vmin = 0, vmax = 1.05*max_GT_R2)\n",
    "axs[0, 4].set_title(vol4title, fontsize=10, fontweight=font_lw)\n",
    "axs[0, 4].set_xticks([])\n",
    "axs[0, 4].set_yticks([])\n",
    "axs[0, 4].set_aspect(asp_ratio)\n",
    "im.set_clim(0, (1.05*max_GT_R2))\n",
    "axs[0, 4].set_xlabel('SSIM = '+str(round(ssim(np.abs(vol4_high_R2[:,:,:]),np.abs(GT_R2[:,:,:])),2))\n",
    "                     +'\\n NRMSE = '+str(round(nrmse(np.abs(vol4_high_R2[:,:,:]),np.abs(GT_R2[:,:,:])),3)),fontsize=7)\n",
    "\n",
    "im = axs[1, 4].imshow(np.abs(vol4_low_R2[:,:,slc].T), cmap = 'gray', vmin = 0, vmax = 1.05*max_GT_R2)\n",
    "#axs[1, 4].set_title(vol4title + ' Diff')\n",
    "axs[1, 4].set_xticks([])\n",
    "axs[1, 4].set_yticks([])\n",
    "axs[1, 4].set_aspect(asp_ratio)\n",
    "im.set_clim(0, (1.05*max_GT_R2))\n",
    "axs[1, 4].set_xlabel('SSIM = '+str(round(ssim(np.abs(vol4_low_R2[:,:,:]),np.abs(GT_R2[:,:,:])),2))\n",
    "                     +'\\n NRMSE = '+str(round(nrmse(np.abs(vol4_low_R2[:,:,:]),np.abs(GT_R2[:,:,:])),3)),fontsize=7)\n",
    "\n",
    "\n",
    "im = axs[1,0].imshow(np.abs(mask_R2.T),cmap='gray')\n",
    "axs[1, 0].set_xticks([])\n",
    "axs[1, 0].set_yticks([])\n",
    "axs[1, 0].set_xlabel('PE1')\n",
    "axs[1, 0].set_ylabel('PE2')\n",
    "axs[1, 0].set_title('Mask',fontsize=10, fontweight=font_lw)\n",
    "#axs[1, 0].axis('off')\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "max_GT_R4 = np.amax(np.abs(GT_R4[:,:,slc]))\n",
    "\n",
    "im = axs[2, 0].imshow(np.abs(GT_R4[:,:,slc].T), cmap= 'gray', vmin = 0, vmax = 1.05*max_GT_R4)\n",
    "\n",
    "axs[2, 0].set_xlabel('PE1')\n",
    "axs[2, 0].set_ylabel('RO')\n",
    "im.set_clim(0, (1.05*max_GT_R4))\n",
    "\n",
    "im = axs[2, 1].imshow(np.abs(vol1_high_R4[:,:,slc].T), cmap = 'gray', vmin = 0,vmax = 1.05*max_GT_R4)\n",
    "axs[2, 1].set_title(vol1title,fontsize=10,fontweight=font_lw)\n",
    "im.set_clim(0, (1.05*max_GT_R4))\n",
    "axs[2, 1].set_xlabel('SSIM = '+str(round(ssim(np.abs(vol1_high_R4[:,:,:]),np.abs(GT_R4[:,:,:])),2))\n",
    "                     +'\\n NRMSE = '+str(round(nrmse(np.abs(vol1_high_R4[:,:,:]),np.abs(GT_R4[:,:,:])),3)),fontsize=7)\n",
    "axs[2, 1].set_ylabel('SNR '+str(noise_dB_zero[-1])+' dB',fontsize=10, fontweight=font_lw)\n",
    "\n",
    "im = axs[3, 1].imshow(np.abs(vol1_low_R4[:,:,slc].T), cmap = 'gray', vmin = 0, vmax = 1.05*max_GT_R4)\n",
    "#axs[3, 1].set_title(vol1title + ' Diff')\n",
    "im.set_clim(0, (1.05*max_GT_R4))\n",
    "axs[3, 1].set_xlabel('SSIM ='+str(round(ssim(np.abs(vol1_low_R4[:,:,:]),np.abs(GT_R4[:,:,:])),2))\n",
    "                     +'\\n NRMSE ='+str(round(nrmse(np.abs(vol1_low_R4[:,:,:]),np.abs(GT_R4[:,:,:])),3)),fontsize=7)\n",
    "axs[3, 1].set_ylabel('SNR '+str(noise_dB_zero[0])+' dB',fontsize=10, fontweight=font_lw)\n",
    "\n",
    "\n",
    "im = axs[2, 2].imshow(np.abs(vol2_high_R4[:,:,slc].T), cmap = 'gray', vmin = 0, vmax = 1.05*max_GT_R4)\n",
    "axs[2, 2].set_title(vol2title,fontsize=10,fontweight=font_lw)\n",
    "im.set_clim(0, (1.05*max_GT_R4))\n",
    "axs[2, 2].set_xlabel('SSIM = '+str(round(ssim(np.abs(vol2_high_R4[:,:,:]),np.abs(GT_R4[:,:,:])),2))\n",
    "                     +'\\n NRMSE = '+str(round(nrmse(np.abs(vol2_high_R4[:,:,:]),np.abs(GT_R4[:,:,:])),3)),fontsize=7)\n",
    "\n",
    "im = axs[3, 2].imshow(np.abs(vol2_low_R4[:,:,slc].T), cmap = 'gray', vmin = 0, vmax = 1.05*max_GT_R4)\n",
    "#axs[3, 2].set_title(vol2title + ' Diff')\n",
    "im.set_clim(0, (1.05*max_GT_R4))\n",
    "axs[3, 2].set_xlabel('SSIM = '+str(round(ssim(np.abs(vol2_low_R4[:,:,:]),np.abs(GT_R4[:,:,:])),2))\n",
    "                     +'\\n NRMSE = '+str(round(nrmse(np.abs(vol2_low_R4[:,:,:]),np.abs(GT_R4[:,:,:])),3)),fontsize=7)\n",
    "\n",
    "im = axs[2, 3].imshow(np.abs(vol3_high_R4[:,:,slc].T), cmap = 'gray', vmin = 0, vmax = 1.05*max_GT_R4)\n",
    "axs[2, 3].set_title(vol3title,fontsize=10,fontweight=font_lw)\n",
    "\n",
    "im.set_clim(0, (1.05*max_GT_R4))\n",
    "axs[2, 3].set_xlabel('SSIM = '+str(round(ssim(np.abs(vol3_high_R4[:,:,:]),np.abs(GT_R4[:,:,:])),2))\n",
    "                     +'\\n NRMSE = '+str(round(nrmse(np.abs(vol3_high_R4[:,:,:]),np.abs(GT_R4[:,:,:])),3)),fontsize=7)\n",
    "\n",
    "im = axs[3, 3].imshow(np.abs(vol3_low_R4[:,:,slc].T), cmap = 'gray', vmin = 0, vmax = 1.05*max_GT_R4)\n",
    "#axs[3, 3].set_title(vol3title + ' Diff')\n",
    "\n",
    "im.set_clim(0, (1.05*max_GT_R4))\n",
    "axs[3, 3].set_xlabel('SSIM = '+str(round(ssim(np.abs(vol3_low_R4[:,:,:]),np.abs(GT_R4[:,:,:])),2))\n",
    "                     +'\\n NRMSE = '+str(round(nrmse(np.abs(vol3_low_R4[:,:,:]),np.abs(GT_R4[:,:,:])),3)),fontsize=7)\n",
    "\n",
    "im = axs[2, 4].imshow(np.abs(vol4_high_R4[:,:,slc].T), cmap = 'gray', vmin = 0, vmax = 1.05*max_GT_R4)\n",
    "axs[2, 4].set_title(vol4title,fontsize=10,fontweight=font_lw)\n",
    "\n",
    "im.set_clim(0, (1.05*max_GT_R4))\n",
    "axs[2, 4].set_xlabel('SSIM = '+str(round(ssim(np.abs(vol4_high_R4[:,:,:]),np.abs(GT_R4[:,:,:])),2))\n",
    "                     +'\\n NRMSE = '+str(round(nrmse(np.abs(vol4_high_R4[:,:,:]),np.abs(GT_R4[:,:,:])),3)),fontsize=7)\n",
    "\n",
    "im = axs[3, 4].imshow(np.abs(vol4_low_R4[:,:,slc].T), cmap = 'gray', vmin = 0, vmax = 1.05*max_GT_R4)\n",
    "#axs[3, 4].set_title(vol4title + ' Diff')\n",
    "\n",
    "im.set_clim(0, (1.05*max_GT_R4))\n",
    "axs[3, 4].set_xlabel('SSIM = '+str(round(ssim(np.abs(vol4_low_R4[:,:,:]),np.abs(GT_R4[:,:,:])),2))\n",
    "                     +'\\n NRMSE = '+str(round(nrmse(np.abs(vol4_low_R4[:,:,:]),np.abs(GT_R4[:,:,:])),3)),fontsize=7)\n",
    "\n",
    "\n",
    "im = axs[3,0].imshow(np.abs(mask_R4.T),cmap='gray')\n",
    "axs[3, 0].set_xlabel('PE1')\n",
    "axs[3, 0].set_ylabel('PE2')\n",
    "axs[3, 0].set_title('Mask',fontsize=10,fontweight=font_lw)\n",
    "\n",
    "\n",
    "\n",
    "img_scale = 1.6\n",
    "\n",
    "for j in range(4):\n",
    "    for k in range(5):\n",
    "        l, b, w, h = axs[j,k].get_position().bounds\n",
    "        axs[j,k].set_position([l,b,w*img_scale,h*img_scale])\n",
    "        if j>1:\n",
    "            axs[j,k].set_position([l,b-0.15,w*img_scale,h*img_scale])\n",
    "\n",
    "col0shift=-0.05\n",
    "col_delta=-0.013\n",
    "\n",
    "axs[0, 0].set_title('Reference',fontsize=10, fontweight=font_lw)\n",
    "axs[2, 0].set_title('Reference',fontsize=10,fontweight=font_lw)\n",
    "\n",
    "l, b, w, h = axs[0,0].get_position().bounds\n",
    "axs[0, 0].set_position([l+col0shift+col_delta,b,w,h])\n",
    "\n",
    "l, b, w, h = axs[1,0].get_position().bounds\n",
    "axs[1, 0].set_position([l+col0shift,b,w,h])\n",
    "\n",
    "l, b, w2, h2 = axs[1,0].get_position().bounds\n",
    "axs[1,0].set_position([l,b,w2/img_scale,h2/img_scale])\n",
    "\n",
    "l, b, w, h = axs[2, 0].get_position().bounds\n",
    "axs[2, 0].set_position([l+col0shift+col_delta,b,w,h])\n",
    "\n",
    "l, b, w2, h2 = axs[3,0].get_position().bounds\n",
    "axs[3,0].set_position([l,b,w2/img_scale,h2/img_scale])\n",
    "\n",
    "l, b, w, h = axs[3,0].get_position().bounds\n",
    "axs[3, 0].set_position([l+col0shift,b,w,h])\n",
    "\n",
    "axs[1, 1].set_title('(a) $R = 2$ undersampling',x=-1.8,y=2.8,fontweight=font_lw)\n",
    "axs[3, 1].set_title('(b) $R = 4$ undersampling',x=-1.8,y=2.8,fontweight=font_lw)\n",
    "\n",
    "for ax in axs.flatten():\n",
    "\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "    ax.set_aspect(asp_ratio)\n",
    "    ax.tick_params(color='k')\n",
    "    if ax not in [axs[1,0],axs[3,0]]:\n",
    "        for spine in ax.spines.values():\n",
    "            spine.set_edgecolor('k')\n",
    "        ax.set_xlim([10,65])\n",
    "        ax.set_ylim([59,5])\n",
    "\n",
    "for ax in [axs[1,0],axs[3,0]]:\n",
    "    #for spine in ax.spines.values():\n",
    "    #    spine.set_edgecolor('w')\n",
    "    ax.set_xticks([0,75])\n",
    "    ax.set_yticks([0,25])\n",
    "    #ax.set_xlim([])\n",
    "    #ax.set_ylim([])\n",
    "\n",
    "figimgs.savefig('export_figs/fastmri_imgs_compare_snr_R.pdf',bbox_inches='tight',dpi=300)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venvkerasandtf",
   "language": "python",
   "name": ".venvkerasandtf"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}