notebook_fatwater_TEs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Noise Propagation in CSE-MRI\n",
    "\n",
    "    Author: Diego Hernando, PhD, University of Wisconsin-Madison\n",
    "    Date created: November 25, 2021\n",
    "    Date last modified: December 8, 2021\n",
    "    \n",
    "## Fat quantification signal model (magnitude signal shown): \n",
    "\n",
    "\n",
    "<font size=\"4\">\n",
    "\n",
    "<br>\n",
    "\n",
    "$s(TE) = \\left| \\left( A_{\\hbox{water}} + A_{\\hbox{fat}} \\sum_{p=1}^P \\rho_p e^{i 2 \\pi TE f_p} \\right) e^{-  TE \\, R_2^* } \\right|$\n",
    "\n",
    "**Known parameters**: TE values (say, 6 echo times), $\\rho_p$ and $f_p$ (fat NMR peak relative amplitudes and frequencies)\n",
    "\n",
    "\n",
    "**Unknown parameters**: $A_{\\hbox{water}}$, $A_{\\hbox{fat}}, R_2^*$\n",
    "\n",
    "**Fat quantification**: FF = $ A_{\\hbox{fat}} / (A_{\\hbox{fat}} + A_{\\hbox{water}})$\n",
    "    \n",
    "**Note**: The signal fitting code below measures FF directly with the formulation:\n",
    "    \n",
    "$s(TE) =  \\left| A \\left( (1-FF) + FF\\sum_{p=1}^P \\rho_p e^{i 2 \\pi TE f_p} \\right) e^{-  TE \\, R_2^* } \\right|$    \n",
    "\n",
    "where $A$ is the overall signal amplitude\n",
    "\n",
    "</font>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Load some python packages\n",
    "##\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.legend import Legend\n",
    "from matplotlib import lines\n",
    "from IPython.display import display, Image\n",
    "from dataclasses import dataclass\n",
    "from scipy.optimize import curve_fit\n",
    "import scipy.special as special\n",
    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
    "import ipywidgets as widgets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Declare a couple helpful classes and functions\n",
    "\n",
    "@dataclass\n",
    "class trueParams: # Tissue parameters to simulate chemical shift encoded signal\n",
    "    A: float\n",
    "    ff: float\n",
    "    r2star: float\n",
    "    fieldmap: float = 0.0\n",
    "        \n",
    "@dataclass\n",
    "class imDataParams: # Acquisition parameters to simulate chemical shift encoded signal\n",
    "    FieldStrength: float\n",
    "    TE: list\n",
    "        \n",
    "\n",
    "@dataclass\n",
    "class algoParams: # Specifies signal model for chemical shift encoded signal\n",
    "    fatFreqs: list\n",
    "    fatAmps: list\n",
    "    USE_R2: bool    \n",
    "        \n",
    "\n",
    "# This function synthetizes some CSE signal given tissue, acquisition, and signal model        \n",
    "def create_signal(imP, trueP, algoP):\n",
    "    sw = trueP.A*(1-trueP.ff*0.01)*np.ones(len(imP.TE))\n",
    "    sf = 0*imP.TE\n",
    "    \n",
    "    for n in range(0, len(algoP.fatFreqs)):\n",
    "        sf = sf + algoP.fatAmps[n]*trueP.A*trueP.ff*0.01*np.exp(1j*2*np.pi*imP.TE*algoP.fatFreqs[n]*imP.FieldStrength*42.576)\n",
    "        #sf = sf + imP.TE\n",
    "    s = sw + sf\n",
    "    \n",
    "    if algoP.USE_R2:\n",
    "        s = s*np.exp(-imP.TE*trueP.r2star + 1j*2*np.pi*imP.TE*trueP.fieldmap)\n",
    "    else:\n",
    "        s = s*np.exp(1j*2*np.pi*imP.TE*trueP.fieldmap)\n",
    "    \n",
    "    return s\n",
    "\n",
    "# Wrapper for create_signal function above     \n",
    "def create_signal_r2star(TE, A, ff, r2star):\n",
    "    \n",
    "    global myImDataParams\n",
    "    global myAlgoParams\n",
    "    \n",
    "    fieldmap = 0\n",
    "    myTrueParams = trueParams(A, ff, r2star,fieldmap)\n",
    "\n",
    "    s = create_signal(myImDataParams, myTrueParams, myAlgoParams)\n",
    "    \n",
    "    s2 = np.abs(s)\n",
    "       \n",
    "    return s2\n",
    "\n",
    "\n",
    "# Wrapper for create_signal function above, without including R2* in the model     \n",
    "def create_signal_nor2star(TE, A, ff):\n",
    "    \n",
    "    global myImDataParams\n",
    "    global myAlgoParams\n",
    "    \n",
    "    r2star = 0\n",
    "    fieldmap = 0\n",
    "    myTrueParams = trueParams(A, ff, r2star,fieldmap)\n",
    "\n",
    "    s = create_signal(myImDataParams, myTrueParams, myAlgoParams)\n",
    "    \n",
    "    s2 = np.abs(s)\n",
    "    \n",
    "    return s2\n",
    "\n",
    "\n",
    "# Make the fonts large\n",
    "font = {'weight' : 'normal',\n",
    "        'size'   : 26}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqgAAAHoCAYAAACFCB2tAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABss0lEQVR4nO3deVxU9f4/8NcMm4DsiywigqgobkmIK+6pLZpLamWut67Z1bppu5WlXevW9aZ+va1qlllqeW9pi6WJ4g4qoOKGbArIvsrOnN8f/DgxLDNnmJkzA7yej8c8mjPncz7zZiTn5eecz+coBEEQQERERERkJpSmLoCIiIiIqCEGVCIiIiIyKwyoRERERGRWGFCJiIiIyKwwoBIRERGRWWFAJSIiIiKzwoBKRERERGaFAZWIiIiIzIqlqQsg01GpVMjIyICDgwMUCoWpyyEiIqJ2ThAElJSUwMfHB0ply+OkDKgdWEZGBvz8/ExdBhEREXUwt27dQteuXVvcz4DagTk4OACo+yVxdHQ0cTVERETU3hUXF8PPz0/MIC1hQO3A6k/rOzo6MqASERGRbLRdWshJUkRERERkVhhQiYiIiMistOuAWlJSgqioKHz44YeYP38+QkJCYGlpCYVCAYVCge7du+v9HmlpaXjnnXcwdOhQ+Pj4wMbGBn5+fhg9ejQ2btyInJwcs+ybiIiIyFwpBEEQTF2EMfTu3Rs3btyAph/P398fKSkprX6PzZs346WXXkJ5eXmLbVxdXfHRRx9h9uzZZtN3veLiYjg5OaGoqIjXoBIREZHRSc0e7XYE9fr16xrDqb5ef/11rFixQi1ABgUFYfTo0QgKChJfy8/Px5w5c7B9+3az6JuIiIjI3LXbgFrP3t4ew4cPx9/+9jds374dkydP1rvPffv2Yd26deJ23759ERMTgxs3biAyMhI3btxATEwM+vbtK7Z56qmncPbsWZP2TURERNQWtNtT/F9//TUGDx6M3r17q92pYOHChdixYweA1p3ir66uRnBwMJKSkgAAXbt2RVxcHFxdXZu0zc/Px4ABA5Ceng4AGDVqFI4dO2aSvpvDU/xEREQkpw5/iv/xxx9Hnz59NN5GqzV27dolBkgA2LBhQ7MBEqi7RnTDhg3idlRUlMYQacy+iYiIiNqKdhtQjWXv3r3icx8fH0yfPl1j+xkzZsDHx6fZ4+Xsm4iIiKitYEDVQXl5OQ4dOiRuT548GZaWmm/GZWlpiUmTJonbP/74o+x9m5IgCLhTegeHkw5j85nNOHnrpKlLIiIiIjPHW53qICEhAZWVleL2iBEjJB03YsQIcaZ9Wloa8vPzm5y6N2bfcopKjcK5zHNIyEkQHwUVBeL+lcNWYrjfcJPVR0REROaPAVUHCQkJats9e/aUdFzjdgkJCRg5cqRsfctpzdE1+CP5jxb3J+QktLiPiIiICOApfp00nvHfrVs3Scf5+/tr7MfYfcspxCNE434GVCIiItKGI6g6KC4uVtt2dnaWdJyTk5PadklJiax916usrFS7jKDxexpCX4++GvenFqWitKoUna07G/y9iYiIqH3gCKoOSktL1bZtbW0lHde4XXMh0ph911u/fj2cnJzEh5+fn6T30IW2gAoAV3OvGvx9iYiIqP1gQNVBdXW12ra2Wfb1rKysNPZj7L7rvfLKKygqKhIft27dkvQeupASUC9nXzb4+xIREVH7wVP8OrC3t1fbrqiogJ2dndbjysvLNfZj7L7r2djYwMbGRmuf+nC3c4eHnQdyynJabMPrUImIiEgTjqDqoHNn9esmy8rKJB3XuJ2Dg4OsfcstxFPLRKlcBlQiIiJqGQOqDjw8PNS2MzMzJR3XuJ27u7usfcutr7vm0/wcQSUiIiJNGFB1EBwcrLadmpoq6bjG7Rr3Y+y+5abtOtTkgmSUVUsbISYiIqKOhwFVByEh6qeuz58/L+m4hu2sra0RFBQka99y0xZQBQicyU9EREQtYkDVgZ+fHwIDA8Xto0ePSjquYbuRI0fCwsJC1r7lpu0aVICn+YmIiKhlDKg6mj59uvg8MjISaWlpGtunpaWphciZM2eapG85edh5wM3WTWMbBlQiIiJqCQOqjhYtWgSlsu5jU6lUWLt2rcb2b7/9NlQqFYC6mfqzZ882Sd9yUigUWk/zM6ASERFRSxhQdRQSEoJ58+aJ259//jk+++yzZtt++umn2Lp1q7i9atUqjbPsjdm33LQF1Ms5XKyfiIiImqcQBEEwdRHGsG7dOqxbt67J69XV1eKoI4BmF65/4oknWgyGAJCdnY2hQ4ciOTlZfO2hhx7C3Llz4ePjg4yMDHz77bfYv3+/uD8sLAyRkZFaF983Zt+NFRcXw8nJCUVFRXB0dNTpWG02ndmEZ399tsX9SoUSpa+UwtZK2i1diYiIqO2Tmj3a7Z2kampqUFlZqbVdc2003S4UADw9PfHzzz9j0qRJ4nWi+/fvVwuNDQ0YMAAHDhyQFCCN2becQjw0T5RSCSpcz7uOgV4DZaqIiIiI2gqe4m+l4OBgXLx4EcuWLWtyF6h6bm5uWL16NaKjo+Hp6WkWfctF2yl+gNehEhERUfPa7Sl+OZWXlyMyMhKpqakoKCiAu7s7evTogYiICFha6jdIbcy+jXmKXxAEuP7TFYUVhS22WT1qNdaO0zwRjIiIiNqPDn+KX062traYMmVKm+vbmOpn8p+8dbLFNpwoRURERM3hKX4ymr7uXGqKiIiIdMeASkaj7Y5SifmJqKzRPpGNiIiIOhYGVDIabROlaoVa3Mi/IVM1RERE1FYwoJLRSJnJfzmb16ESERGROgZUMhpfB184WDtobMPrUImIiKgxBlQymvqZ/Jok5DKgEhERkToGVDIqbXeU4ggqERERNcaASkalbQT1et51VNdqvrUsERERdSwMqGRU2gJqjaqGM/mJiIhIDQMqGZWUmfw8zU9EREQNMaCSUfk5+cHeyl5jGwZUIiIiaogBlYxKqVBqn8nPgEpEREQNMKCS0TGgEhERkS4YUMnotAXUa3nXUKOqkakaIiIiMncMqGR02gJqVW0VbubflKkaIiIiMncMqGR02hbrB3ian4iIiP7EgEpG5+/sD1tLW41tGFCJiIioHgMqGZ1SoUQfjz4a21zOuSxTNURERGTuGFBJFpzJT0RERFIxoJIs+rprDqhXc6+iVlUrUzVERERkzhhQSRYhnponSlXWViK5MFmmaoiIiMicMaCSLLSd4gd4mp+IiIjqMKCSLAKcA2BjYaOxzeVsTpQiIiIiBlSSiYXSAsHuwRrbJORyBJWIiIgYUElGnMlPREREUjCgkmy03VHqSs4VqASVTNUQERGRuWJAJdloG0EtrylHcgFn8hMREXV0DKgkG21LTQHApexLMlRCRERE5owBlWTTw6WH1pn8DKhERETEgEqysVBaaD3NfymHAZWIiKijY0AlWfXz7KdxP0dQiYiIiAGVZKUtoF7LvYbq2mqZqiEiIiJzxIBKstK21FS1qho38m/IVA0RERGZIwZUkpW2EVSAp/mJiIg6OgZUklU3p27obN1ZYxsGVCIioo6NAZVkpVAoOFGKiIiINGJAJdn182BAJSIiopYxoJLstI2gJuYnory6XKZqiIiIyNwwoJLstN3yVICAq7lXZaqGiIiIzA0DKsmOM/mJiIhIEwZUkl0X+y5ws3XT2IYBlYiIqONiQCXZSZrJn8OASkRE1FExoJJJcKkpIiIiagkDKpmEtoCaVpSG4spimaohIiIic8KASiYhZaLU5ezLMlRCRERE5oYBlUwixEPzUlMAcDmHAZWIiKgjYkAlk3CxdYGPg4/GNrwOlYiIqGNiQCWT4UQpIiIiag4DKplMPw8GVCIiImqKAZVMRtsIatbdLOTczZGpGiIiIjIXDKhkMpJm8nOiFBERUYfDgEom08ejj9Y2PM1PRETU8TCgksl0tu6MAOcAjW24FioREVHHw4BKJqV1Jn8OR1CJiIg6GgZUMikpS00JgiBTNURERGQOGFDJpLQF1MKKQmSUZMhUDREREZkDBlQyKSkz+TlRioiIqGNhQCWT6u3WGxYKC41tLmZflKkaIiIiMgcMqGRSNpY26OnWU2MbBlQiIqKOhQGVTG5AlwEa98dnxctUCREREZkDBlQdKRQKvR7du3dvse+UlJRW9RkcHCzfB2AEAzw1B9SEnARU11bLVA0RERGZGgOqzNzd3U1dgtnRNoJaVVuF63nXZaqGiIiITM3SGJ1WV1cjLS0N2dnZuHv3LgDA3t4enp6e6NatG6ysrIzxtrKYNGmSTu1TUlJw7do1cXvevHmSj42IiICtra3Wdn5+fjrVZG60BVSg7jrUEM8QGaohIiIiUzNIQM3KysIvv/yCP/74A9HR0bhx40aLi6srFAr06tULYWFhGDt2LKZMmYIuXboYogxZ/Prrrzq1nzZtmhhQra2t8cQTT0g+dseOHRovCWgvujl1g6ONI4ori1tsE58Vj7n95spYFREREZlKqwNqRUUFdu/ejR07duDYsWNqgVTTnX8EQcC1a9dw7do17Ny5EwqFAiNHjsTChQsxZ84cSSOGbcWdO3fw888/i9sPP/ww3NzcTFiReVIoFOjv2R8nbp1osQ0nShEREXUcOl+DWlBQgDfeeAPdunXD4sWLcfToUahUKgiCID60adhWpVIhKioKS5YsgZ+fH15//XXk5eW16ocxN19++SVqamrE7SVLlpiwGvPGmfxERERUT/IIanl5Od5//31s2LABJSUlANRHSn19fTFkyBAMHDgQQUFB6Nq1K1xcXGBnZwdBEFBeXo78/Hzcvn0biYmJiIuLw9mzZ5GRkSH2U1BQgH/84x/YtGkT/v73v+OFF16Avb29gX9k+Wzbtk183q1bN0yYMMGE1Zg3bQH1VvEtFFYUwrmTszwFERERkclICqh79+7FypUrkZ6eLoZJhUKB4cOHY9asWXjggQfQs6fmxdZbcv36dfzyyy/Yu3cvTp06BUEQUFJSgrVr12Lr1q3YsGEDHnnkkVb1bUrHjx9Xmxy1aNEiKJVcNKElkiZKZV3EKP9RMlRDREREpiQpMc2ZM0cMpy4uLli1ahWuXbuG48eP47nnnmt1OAWAXr164dlnn8Xx48dx48YNvPDCC3B1dYUgCEhPT8ejjz7a6r5NqeHoqUKhwKJFi0xYjfnr59lPaxue5iciIuoYJA/pubu74/3330daWhr++c9/IigoyODFBAYG4r333kNaWho++OADeHh4SLqm1dyUlpZiz5494vaECRPg7++vcz8vvvgi+vfvD2dnZ9jY2MDLywthYWF47rnncPToUUOWbHKONo7o7txdYxsGVCIioo5BUkBdvXo1bt68iZUrV8pyTaidnR2ef/55JCUlYfXq1UZ/P0P79ttvxfVfgdZPjtq7dy8uXbqEoqIiVFVVISsrCzExMdi4cSPGjBmD8PBwXLp0yVBlm5zWiVLZDKhEREQdgUJoi0OUZm748OE4deoUAMDV1RUZGRmwsbHRelxKSgoCAgLEbTc3N/To0QMODg4oLS3FzZs3kZubq3aMra0t9uzZgwcffFBr/5WVlaisrBS3i4uL4efnh6KiIjg6Okr98Yzm9T9ex7qodS3u72zdGUUvF0Gp4LW8REREbVFxcTGcnJy0Zg9+0xvY1atXxXAK1N05Sko4rRcaGootW7YgKSkJubm5OHPmDA4dOoTTp08jJycHMTExmDVrlti+vLwcc+bMQWxsrNa+169fDycnJ/Fhbneg0jaCWlpVipTCFHmKISIiIpNhQDWwrVu3qm3rcnq/e/fuiImJwbJly9RGUhsKDQ3F3r17sWnTJvG1srIyrFixQmv/r7zyCoqKisTHrVu3JNcmh/5d+mttw+tQiYiI2j8GVAOqrq7GV199JW7fe++9GDBA+/JJrbF8+XIsXrxY3I6KikJMTIzGY2xsbODo6Kj2MCdBrkHoZNlJYxsGVCIiovaPAdWADhw4gKysLHH7L3/5i1Hf77XXXlPb/vXXX436fsZmqbREiEeIxjYXsy/KVA0RERGZiuQ7SR07dkyvN7K2toaLiwvc3d3b7f3oG659amdnZ/Q1XAMDA+Hv74/U1FQAdde/tnUDugzAucxzLe7nCCoREVH7JzmgjhkzBgqFwiBv2rVrV4waNQpPPfUUIiIiDNKnqWVmZuKXX34Rt2fNmiXLKXRvb28xoDae4d8WaZsodSPvBsqqy2BnZSdTRURERCQ3nU/xC4Kg9+PWrVv45ptvMHbsWIwbNw6ZmZnG+Nlk9cUXX6C2tlbcbu3ap7oqKysTn9va2srynsakLaAKEHA5+7JM1RAREZEp6BRQDb1kqiAIiIyMxODBg5GcnGzQvuW2fft28XlQUJAsI8OVlZVITEwUt728vIz+nsbW31P7TH5eh0pERNS+ST7F3zCAtUZFRQVKSkqQnJyMCxcu4MyZMwDq7lOflZWFhx9+GNHR0bC2ttbrfUzh2LFjuHHjhrgt1+jpvn371EZQR44cKcv7GpOHvQe8OnvhTumdFtvwOlQiIqL2TXJAXbBggUHf+Pbt23j11Vexc+dOAMClS5ewfft2/PWvfzXo+8ih4eQoCwsLg39WzcnOzsbLL78sbtvb22Py5MlGf185DOgygAGViIioAzPZMlNdu3bFl19+iVdeeUV8bePGjaYqp9VKSkqwd+9ecfv++++Ht7e3zv2cOnUKS5cuxbVr17S2vXjxIsaOHYu0tDTxtZUrV7ab1REGeGq+DjU+K97gl5sQERGR+ZA8gmosb7/9Nnbv3o2kpCRcu3YNGRkZ8PHxMXVZkn3zzTdqp9lbe3q/srISn3zyCT755BMMHDgQ48aNw4ABA+Dl5QUHBweUlpYiMTERBw8exE8//QSVSiUeO3HixCZrorZl2iZK5ZXnIbM0Ez4Obef3hIiIiKQzeUC1sLDA3Llz8Y9//AMAcObMGUyfPt3EVUnX8PR+ly5d8MADD+jdZ1xcHOLi4iS1XbBgAbZs2dImr91tidRbnjKgEhERtU9mcSepsLAw8Xl2drYJK9FNQkKCONkLqAuLlpaty/wBAQGYN28egoKCtLa1tLTE1KlTcfjwYXzxxRewt7dv1Xuaqz7ufWChsNDYJu6OtABPREREbY/JR1ABqF07WVBQYMJKdNO3b1+DXQvp7++Pr776CkDdgvvx8fHIyclBXl4eCgoK0KlTJ7i4uKBnz564995728Wapy2xsbRBH48+uJR9qcU2sVmx8hVEREREsjKLgFpZWSk+b0+nqlvL3d0d48aNM3UZJjXIa5DmgHonVr5iiIiISFZmcYo/PT1dfO7q6mrCSshcDOoySOP+a7nXcLfqrjzFEBERkazMIqAeOXJEfB4YGGjCSshcDPIapHG/AEHjCCsRERG1XSYPqOnp6di9ezcAwMrKCuHh4SauiMzBQK+BWtvwND8REVH7ZNKAWlBQgOnTp6OiogIKhQKTJk2CjY2NKUsiM+Fu546ujl01tmFAJSIiap9MElBv3ryJ9957D8HBwTh37pz4+ksvvWSKcshMaTvNz5n8RERE7ZPkWfz6ziqvrKxEaWkpUlJSUFpaCgAQBAEKhQIKhQLPPPMMhg8frtd7UPsyqMsgHLh+oMX98VnxqFXVwkKpec1UIiIialskB9TIyEgoFAq93qzhmqH1wVQQBDz66KP497//rVff1P5oG0Etqy5DYn4ierv3lqcgIiIikoXJrkEVBAHdu3fH559/jq+//hoWFhwFI3XaAirA61CJiIjaI8kjqN26ddNrBNXKygrOzs5wd3dHaGgoIiIiMH78eCiVJl9IgMxUgEsAHKwdUFJV0mKb2DuxmNNvjoxVERERkbFJDqgpKSlGLIOoKaVCiYFeA3E87XiLbThRioiIqP3h8CWZNW13lOIpfiIiovaHAZXMmrbrUO+U3sGd0jvyFENERESyMJuAmp2djQ8++MDUZZCZkTJRKu5OnPELISIiItmYNKCqVCr89NNPmDFjBvz8/LhQPzUR4hkCC4XmFR54mp+IiKh9kTxJypASExOxbds27NixA3fu1J2erV+0n6ihTpad0MejDy5lX2qxDSdKERERtS+yBdTy8nLs3bsX27ZtQ1RUFIA/F+5nMCVNBnkN0hxQOYJKRETUrhg9oEZHR2Pr1q349ttvUVJSt55l49FSLy8vPPLII3j00UeNXQ61QYO6DMJO7Gxx/7Xca7hbdRf21vYyVkVERETGYpSAmpeXh6+++grbtm3D5cuXAajf5hQA3NzcMGPGDMydOxejR4/mKCq1SNtEKQECLmVfQnjXcHkKIiIiIqMyWEAVBAG//fYbtm7dih9//BHV1dXi60Ddafz6kdOffvoJEydO5O1NSZKBXgO1tonLimNAJSIiaif0DqjJycnYvn07vvjiC6SnpwNoOloaHh6OuLg4VFRUAAAmT56s79tSB+Ju546ujl1xu/h2i214HSoREVH70aqAWllZie+//x7btm1DZGQkBEFoEkp9fX3xxBNPYOHChejVqxc8PDzEgEqkq0FegzQG1POZ52WshoiIiIxJp4B64cIFbN26Fbt27UJRUREA9dFSW1tbPPzww1iwYAEmTpzI60rJYAZ1GYQD1w+0uD8+Kx41qhpYKk2ychoREREZkORv88GDByMuru6OPY1HS4cNG4aFCxdizpw5cHR0NGyFRADu8b5H4/7ymnJcybmC/l36y1QRERERGYvkgBobG6u27efnJ57CDwoKMnRdRGpCvUO1tjmfeZ4BlYiIqB3Q6VanCoUCCoUC999/P86cOYN169YxnJIsujl1g6utq8Y25zLPyVQNERERGZNOAbXeL7/8gm7dumHq1Kn4/vvvxSWliIxFoVBoHUVlQCUiImofJAfUhQsXwt7eXpyxX1NTg59++gmzZ8+Gt7c3li9fjujoaGPWSh2ctoAaeycWtapamaohIiIiY5EcULdt24bMzEx89tlnGDZsGACIYTU/Px//+c9/MHToUISEhOD9999HRkaG0Yqmjmmw92CN+8uqy3At75pM1RAREZGx6HSK397eHkuWLMGJEyeQkJCA559/Hp6engD+DKtXr17Fyy+/DH9/f0yZMgW7d+9GZWWlUYqnjiXUR9pEKSIiImrbWnUNKgAEBwfjgw8+wO3bt/Hdd9/h/vvvh4WFhRhUa2tr8dtvv+Gxxx6Dl5cXiouLDVk3dUABzgFw7uSssc25DF6HSkRE1Na1OqDWs7S0xIwZM3DgwAGkpqZi7dq16NGjB4A/R1WLiopQU1MjHrNmzRpcuXJF37emDkahUGg9zc+JUkRERG2f3gG1IR8fH7z22mu4ceMG/vjjDzz22GOwtbUFUBdW6+8stXbtWvTr1w8DBw7Eu+++i+TkZEOWQe3YYC/NAfXCnQtQCSqZqiEiIiJjMGhAbWjMmDHYuXMnMjMzsWXLFoSGhoojqkBdYL106RJee+01BAUFITw8HB9++KGxyqF2Qtt1qKVVpbiRd0OmaoiIiMgYjBZQ6zk6OuLpp59GdHQ0YmNj8be//Q0uLi4AoBZWo6OjsXLlSmOXQ22ctlP8ACdKERERtXVGD6gNDRgwAJs2bUJGRgZ27dqFCRMmyPn21A4EuQbBwdpBYxteh0pERNS2yRpQ61lbW2Pu3Ln47bffkJSUhNdffx1+fn6mKIXaGKVCyYlSRERE7ZxJAmpD/v7+eOutt5CSkoJffvnF1OVQG6AtoJ7PPM+JUkRERG2YyQNqPYVCgfvuu8/UZVAboO2Wp8WVxUgqSJKpGiIiIjI0swmoRFJxohQREVH7JimgZmdnG7uOFmVlZZnsvck89XLrBXsre41teEcpIiKitktSQA0KCsLbb7+N0tJSY9cjKi0txZo1a9CzZ0/Z3pPaBgulBQZ5DdLYhhOliIiI2i5JAbW0tBRvvfUW/P398eabbxp1RDUrKwtvvPEG/P39sXbtWty9e9do70Vtl7brUM9nnhfX2SUiIqK2RVJAvffeeyEIAgoKCrBu3Tp069YNc+fOxU8//YTq6mq9i6iqqsKPP/6IRx55BN26dcM777yDgoICCIKAsLAwvfun9kfbdagFFQVIKUyRpxgiIiIyKEkB9cyZM/j444/h4eEBQRBQVVWFvXv3YurUqfD09MSsWbOwZcsWREdHo7y8XGt/ZWVlOHPmDDZt2oTp06fDw8MD06dPx759+1BdXQ1BEODh4YGPP/4Yp06d0vuHpPZH2y1PASA6I1qGSoiIiMjQFIIO50Hv3r2Lf/3rX/jwww9RWFj4ZycKhdpzb29v+Pr6wtXVFba2thAEAeXl5cjPz8ft27eRlZWldvq14XNnZ2c8//zz+Pvf/w57e80TYUg/xcXFcHJyQlFRERwdHU1djk5qVDVwetcJZdVlLbZZNWwV3r/vfRmrIiIiIk2kZg+dAmq9u3fv4vPPP8dHH32E69evN99xg9DaUEtvFxwcjKeffhqLFy9mMJVJWw6oADBq+ygcTzve4v4I/wgcXXhUxoqIiIhIE6nZo1XroNrb2+PZZ5/F1atXcfLkSTz33HNNZtsLgtDso6Hg4GCsXLkSp0+fRkJCApYvX85wSpKF+Wi+PvlcxjnUqmplqoaIiIgMxVLfDoYOHYqhQ4diw4YNyMjIQHR0NBISEpCcnIycnBxxFr69vT08PT0REBCAfv36ISwsDF26dNH7B6COa4jvEI3771bfxZXcK+jn2U+mioiIiMgQ9A6oDfn4+GDatGmYNm2aIbslapa2EVQAiE6PZkAlIiJqY3irU2qzAl0C4WrrqrHN2fSzMlVDREREhsKASm2WQqHQOorKpaaIiIjaHgZUatO0XYcanxWPipoKmaohIiIiQ2BApTZN2whqtaoacXfiZKqGiIiIDIEBldq0MF8JE6V4mp+IiKhNYUClNs2rsxf8HP00tuFEKSIioraFAZXaPG3XoXIElYiIqG1hQKU2T9t1qFdzr6KookimaoiIiEhfDKjU5mkbQQWAc5nnZKiEiIiIDIEBldq8UJ9QKKDQ2CY6naf5iYiI2goGVB2lpKRAoVDo/AgODtbpfdLS0vDOO+9g6NCh8PHxgY2NDfz8/DB69Ghs3LgROTk5RvoJ2x5HG0cEu2v+fM9mcKIUERFRW8GAaoY2b96M4OBgrF69GmfOnEFmZiaqqqpw+/ZtHDt2DM899xyCg4OxZ88eU5dqNrQtN8URVCIiorbD0tQFtHURERGwtbXV2s7PT/NSSPVef/11rFu3Tu21oKAg+Pr6Ij09HYmJiQCA/Px8zJkzB3fv3sWiRYt0L7ydGeIzBF/Gfdni/lvFt5BZkglvB28ZqyIiIqLWYEDV044dO9C9e3eD9LVv3z61cNq3b198+eWXCA0NFV87d+4c5s+fj4SEBADAU089hZCQEAwZon2iUHsmZcH+07dPY3qf6TJUQ0RERPrgKX4zUV1djRdeeEHc7tq1K6KiotTCKQCEhoYiKioKvr6+AICamhqsWrVK1lrN0cAuA2FtYa2xzenbp2WqhoiIiPTBgGomdu3ahaSkJHF7w4YNcHV1bbatq6srNmzYIG5HRUXh2LFjRq/RnNlY2mCw92CNbU7dPiVTNURERKQPSaf45Qw/ERERsr2XOdm7d6/43MfHB9Onaz4VPWPGDPj4+CAjI0M8vqN+dvWGdR2mcZQ0JiMG1bXVsLKwkrEqIiIi0pWkgDpmzBgoFJrXmTQEhUKBmpoao7+PuSkvL8ehQ4fE7cmTJ8PSUvMfjaWlJSZNmoTt27cDAH788Uds3rzZqHWau6Fdh2rcX15TjviseIT6hGpsR0RERKZlNqf4BUGAIAimLsMkEhISUFlZKW6PGDFC0nEN26WlpSE/P9/gtbUlw7oO09qGp/mJiIjMn+SAWh8gjfVoq1588UX0798fzs7OsLGxgZeXF8LCwvDcc8/h6NGjkvqon5Ffr2fPnpKOa9yucT8djZ+TH3wdfDW2YUAlIiIyf5ICqkqlkuVRW1tr7J/X4Pbu3YtLly6hqKgIVVVVyMrKQkxMDDZu3IgxY8YgPDwcly5d0thHSkqK2na3bt0kvbe/v7/Gfjoibaf5OZOfiIjI/JnNKf62ys3NDUOGDMH48eMRHh4Od3d3tf1nz57FkCFDcODAgRb7KC4uVtt2dnaW9N5OTk5q2yUlJRrbV1ZWori4WO3R3mg7zZ9UkITsu9kyVUNEREStwYDaCqGhodiyZQuSkpKQm5uLM2fO4NChQzh9+jRycnIQExODWbNmie3Ly8sxZ84cxMbGNttfaWmp2raUO1M1105bQF2/fj2cnJzEh9S7W7Ulw/wkXId6i6f5iYiIzBkDqo66d++OmJgYLFu2DAEBAc22CQ0Nxd69e7Fp0ybxtbKyMqxYsaLZ9tXV1Wrb2mbw17OyUl8uqXE/jb3yyisoKioSH7du3ZL0Pm3JYO/BsFJqXkaKp/mJiIjMGwOqES1fvhyLFy8Wt6OiohATE9Oknb29vdp2RUWFpP7Ly8s19tOYjY0NHB0d1R7tTSfLTrjH+x6NbThRioiIyLwxoBrZa6+9prb966+/NmnTuXNnte2ysjJJfTdu5+DgoGN17ZO261CjM6JRo+p46+0SERG1FdLOJetIEAQkJCTg9u3bKCwsVFvjU5v58+cboySTCQwMhL+/P1JTUwEAV69ebdLGw8NDbTszM7PJZKvmZGZmqm1LOaYjGNp1KDae2dji/rLqMlzMuqh1pJWIiIhMw6ABNTY2Fv/617/w448/Npn4I1V7C6gA4O3tLQbU3NzcJvuDg4PVtlNTU9G/f3+t/db32VI/HZXUBfsZUImIiMyTwU7xv/feexgyZAh27dqFkpISjQvyt7fF+rVpeCq+uRn6ISEhatvnz5+X1G/DdtbW1ggKCmplhe1LN6du8O7srbENr0MlIiIyXwYJqB9//DFeeeUV1NTUiEHTzc0NnTp1AgAoFAr4+/vD1dUVCoVCPE6hUEChUMDPzw/+/v6SF6hvSyorK5GYmChue3l5NWnj5+eHwMBAcVvqHagaths5ciQsLCz0qLT9UCgUXLCfiIioDdM7oObn5+OFF14AUBcMHnjgAVy/fh05OTkYPXq02C45ORm5ubkoKirCwYMH8cgjj4j7BgwYgLi4OCQnJ+tbjtnZt2+f2gjqyJEjm203ffp08XlkZCTS0tI09puWlqYWUGfOnKlnpe2LttP8ifmJXLCfiIjITOkdULdt24a7d+9CoVAgIiICP/zwg8ZTzfb29pg4cSJ2796N3377DY6Ojvj5558xYcIE1NS0r5nV2dnZePnll8Vte3t7TJ48udm2ixYtglJZ98ehUqmwdu1ajX2//fbbUKlUAOpWAZg9e7aBqm4fpCzYfyLthAyVEBERka70DqhHjhwRn7/77rtiyJJi/Pjx+O9//wsAOHfuHNasWaNvOUZ16tQpLF26FNeuXdPa9uLFixg7dqzaSOjKlSvh5ubWbPuQkBDMmzdP3P7888/x2WefNdv2008/xdatW8XtVatWcQZ/I6HeoVoX7D+edlymaoiIiEgXCkHP2Un+/v64desW3N3dkZ2tfsp0ypQpOHjwIBQKBWpqatSuP21o1qxZ2LdvH1xdXZGVlWW211JGRkZi7NixAICBAwdi3LhxGDBgALy8vODg4IDS0lIkJibi4MGD+Omnn8QRTgCYOHEiDhw4AGtr6xb7z87OxtChQ9UudXjooYcwd+5c+Pj4ICMjA99++y32798v7g8LC0NkZCTs7Ox0/nmKi4vh5OSEoqKidrlo//CtwzVOhgrzCcPZJ8/KWBEREVHHJjV76L3MVH5+PhQKBXr16tVkX8OgWV5e3mKImjZtGvbt24eCggIcPXoU48aN07cso4uLi0NcXJyktgsWLMCWLVs0hlMA8PT0xM8//4xJkyaJI6/79+9XC6QNDRgwAAcOHGhVOO0IRnYbqTGgns88j7tVd2FvrfkOXERERCQvvU/xV1VVAWj+NpsN72yUk5PTYh/+/v7i85SUFH1LMpqAgADMmzdP0nJOlpaWmDp1Kg4fPowvvvhC621I6wUHB+PixYtYtmxZkztM1XNzc8Pq1asRHR0NT09PnX6GjmRUt1Ea99cKtZzNT0REZIb0HkF1cnJCXl5eswvzu7q6is+TkpLUgmhD9SEX0BxkTc3f3x9fffUVgLoF9+Pj45GTk4O8vDwUFBSgU6dOcHFxQc+ePXHvvfc2u+apFI6OjtiyZQs++OADREZGIjU1FQUFBXB3d0ePHj0QEREBS0uj3ASsXRnuN1xrm+NpxzE+cLwM1RAREZFUeqecgIAA5ObmIiMjo8m+vn37is+PHTsmXr/Z2JkzZ8TnrQ11cnN3dzf6pQi2traYMmWKUd+jPXOzc0Nfj75IyElosc3xW5woRUREZG70PsVff0vOtLQ0lJSUqO0bMWKE+Hzr1q1N9gN117D+5z//Ebebu5aVqLW0neY/desUalTta3kzIiKitk7vgBoRESE+/+OPP9T2DRo0CH369AEApKen47777sOxY8dQXl6O8vJyHDx4EKNHj0ZmZiaAulPbDfsj0tfIbs3fGKHe3eq7iL0TK08xREREJIneAXXSpEni2qffffddk/3r168Xb3969uxZjB07Fp07d0bnzp1x//33IyGh7vSrQqHAqlWrOCOdDErbCCoARKVGyVAJERERSaV3QO3SpQvmz5+PwYMHIzk5We22ngAwdepUvPrqq2JIFQRBfNRvA3W36nz11Vf1LYdITTenbujq2FVjG16HSkREZF70DqhA3e1Oo6Ojcfz48WZHQNetW4f//ve/GDx4sNrrgiCgR48e+M9//oM9e/a0uJA/UWspFAqtp/mPpx2HnverICIiIgOSba2iadOmYdq0acjNzUVaWhqqq6vh4+MDPz8/uUqgDmpUt1H49tK3Le7PvpuNG/k30MuNE/SIiIjMgeyLabq7u/O+8SQrbSOoQN0oKgMqERGReTDIKX4icxbiEQInGyeNbaLSOFGKiIjIXDCgUrtnobTAiG4jNLbhTH4iIiLzwYBKHcJIP82n+W8W3ERGSdO7oREREZH8DH4NamVlJS5cuIDLly+jsLAQZWVlOs2QfuONNwxdEpGk61AjUyLxWP/HZKiGiIiINDFYQM3Ozsbrr7+Ob775Bnfv3m11PwyoZAxhvmHoZNkJFTUVLbY5knyEAZWIiMgMGCSgnjt3DpMnT0Z+fr5e60lyHVQylk6WnTCs6zAcSTnSYpvI1Ej5CiIiIqIW6R1Qi4uL8cADDyAvL08tYHp5eaFr166wt7fX9y2IDGJs97EaA2pifiJuF9/WeucpIiIiMi69A+rmzZuRnZ0NhUIBQRCwZMkSvPjii+jZs6ch6iMymDHdx2htE5kSiXkD5hm/GCIiImqR3rP4Dxw4ID5fvXo1PvvsM4ZTMktDfIfA1tJWY5vIlEh5iiEiIqIW6R1QExMTAQD29vZYvXq13gURGYuNpQ2G+w3X2EbTJQBEREQkD70DamlpKRQKBQYMGABra2tD1ERkNGO7j9W4P6kgCWlFaTJVQ0RERM3RO6B6e3sDAKysrPQuhsjYpF6HSkRERKajd0Dt168fBEFASkqKAcohMq4w3zDYWdlpbMOASkREZFp6B9QnnngCAJCWlob4+Hi9CyIyJmsLa4zwG6GxDa9DJSIiMi29A+qMGTMwdOhQCIKAv//971CpVIaoi8hotF2HmlKYgpTCFHmKISIioib0DqgWFhb49ttv4efnh8jISMycOROFhYUGKI3IOHgdKhERkXlTCPrcmxTAsWPHAABZWVl45plnkJeXBwcHB8ydOxcjRoyAl5cXbGxsJPcXERGhTzmkg+LiYjg5OaGoqAiOjo6mLkc21bXVcHnPBXer77bY5vH+j2PnjJ0yVkVERNT+Sc0eegdUpVKpdovT+u4aviaVQqFATU2NPuWQDjpqQAWAyTsn4+DNgy3u72LfBZkrM1v1e0xERETNk5o99D7FX69hMG3Nl7ogCNAzKxNJNiFwgsb9WXezcDH7okzVEBERUUMGCaj1wbI+ZLbmQSQnbQEVAH6/+bsMlRAREVFjegdUlUplsEdtba0hfiYirQZ0GQAPOw+NbQ4lH5KpGiIiImrIYKf4idoSpUKpdRT1aMpRVNZUylQRERER1WNApQ5rYuBEjfvLa8px8tZJmaohIiKiegyo1GFJug41idehEhERyY0BlTosPyc/9HbrrbENAyoREZH8GFCpQ9N2mv9cxjnkl+fLVA0REREBgKW+Hbz99tt6F2FjYwMnJyd4eHjgnnvuQWBgoN59EkkxscdE/F/0/7W4X4CAP5L/wKy+s2SsioiIqGPTO6CuWbPG4Hfb8fLywpIlS7B8+XJ4eGheCohIH6P9R8NCYYFaoeUlzn6/+TsDKhERkYwMtlB/SwvvS3294SMzMxPvvPMO+vXrh4MHW74dJZG+nDo5IbxruMY2vyX9xptJEBERyUjvEdSIiAhxBPXSpUvIz88Xg6aVlRV69OgBZ2dnCIKAoqIi3Lx5E9XV1QDqbovq5uaGkJAQVFRUICsrCykpKWLfOTk5ePjhh/HHH39g2LBh+pZK1KyJgRM1LieVUpiCa3nXEOweLGNVREREHZfeI6iRkZH4448/MHz4cDGczpgxA4cOHUJpaSkSEhJw8uRJnDp1CgkJCSgtLcWhQ4cwffp0CIKA/Px8jBgxAidPnkRSUhKysrKwdu1adOrUCQqFApWVlVi4cCFUKpUhfl6iJrRNlAKAn2/8LEMlREREBBjoFP+aNWvw7rvvwtraGnv27MF3332HcePGwcrKqklbKysrjBs3Dt9//z12794NS0tLrF+/HmvWrAEAeHh44LXXXsOxY8dga2sLAEhMTMT3339viFKJmgjvGg7nTs4a2zCgEhERyUfvgBoXF4d33nkHALB27VrMmiV9MskjjzyCtWvXQhAEvPPOO4iLixP3hYaG4vXXXxe3f/nlF31LJWqWpdISk3pM0tjmWOoxlFaVylQRERFRx6Z3QP3888+hUqlga2uLZ555Rufjn3nmGdjZ2UEQBHz++edq+5588knx+tYzZ87oWypRi+7veb/G/dWqahxOOixTNURERB2b3gH1yJEjUCgU6Nevn3hKXhd2dnbo168fBEFAZGSk2j5XV1f06tULgiAgOztb31KJWjQ5aLLWNjzNT0REJA+9A2p6ejqAuqDZWvXH3r59u8m+Ll26AACKiopa3T+RNp72ngjzCdPY5ufEn7ncFBERkQz0Dqg1NTUAoLY8lK7qj63vqyELCwsAaHbCFZEhTQmaonH/7eLbuJR9SaZqiIiIOi69A6qvry8EQUBqaipiYmJ0Pj46OhopKSlQKBTw8fFpsj8/v+4+6G5ubvqWSqSRtutQAeCXRE7WIyIiMja9A+q4cePE50uXLkVpqfSZziUlJVi6dGmzfQGASqXC9evXoVAo0L17d31LJdLoXp974W7nrrENr0MlIiIyPr0D6pIlS8SZ9hcuXMCoUaNw9uxZrcedPXsWERERiI2NBVB3V6knn3xSrc2ZM2dQXl4OABg8eLC+pRJpZKG00DpZ6njacRRV8HpoIiIiY9L7VqehoaFYtmwZtmzZAoVCgbi4OAwbNgyDBg3C+PHj0adPHzg7O0OhUKCgoABXrlzB4cOHxWAK1IXTZcuWNQmhe/fuFZ+PHTtW31KJtLo/6H7sjN/Z4v5aoRYHbx7E7JDZMlZFRETUsegdUAFg48aNyM3Nxe7du6FQKCAIAmJjY9VCaGOCIIht586di40bNzZpU1hYiGnTpkGpVGLSJM0LqRMZwn097oNSoYRKaPnWuj9c+4EBlYiIyIgUggHXzfnoo4/w6quviktC1Z/6b6z+LR0dHfHuu++qXYdK8ikuLoaTkxOKiorg6Oho6nLMxohtI3Dy1skW9zvZOCHnhRxYWXBlCSIiIl1IzR56X4Pa0NNPP43bt2/jk08+wUMPPQQPDw8IgqD28PDwwIMPPohPPvkE6enpDKdkdqb2mqpxf1FlEY6mHpWpGiIioo7HoCOozamqqkJBQQEAwMXFBdbW1sZ8O9IBR1Cbdy33GoK3BGts80zYM/i/+/9PpoqIiIjaB5OMoDbH2toaXbp0QZcuXRhOqU3o7d4bwe6aA+oP137gXaWIiIiMxOgBlagtmtZ7msb9t4tv43zmeZmqISIi6lgYUIma8XDww1rb/O/q/4xeBxERUUfEgErUjCG+Q+DV2Utjmx+u/SBTNURERB2LpHVQ3377bbXtN954o8V9+mrYN5GpKBVKTO01FZ+e/7TFNhezL+Jm/k30cO0hY2VERETtn6RZ/EqlUm1N09ra2hb36ath32RcnMWv2S83fsH9u+7X2OZf9/0Lzw97XqaKiIiI2jajzOJvKcs2Xuu0tQ8iczIuYBw6W3fW2GbflX0yVUNERNRxSDrFHxER0eIoqaZ9RG2ZjaUNpgRNwd6EvS22OXHrBNKL0+Hr6CtjZURERO2bpIAaGRnZqn3tVWVlJU6cOIHIyEicP38eCQkJyMnJQWVlJZycnNC1a1eEh4dj+vTpuO+++yQF+JSUFAQEBOhcS+/evXH16tXW/BgkwcPBD2sMqADw/ZXvsSJ8hUwVERERtX+SAirVycrKwt///nccOHAAJSUlzbbJzc1Fbm4uYmNj8cknnyAkJARbt25FeHi4zNWSITzY60HYWNigsrayxTZ7Lu9hQCUiIjIgBlQd3Lp1C998802T1729vdG1a1c4ODjgzp07uHr1KlQqFQDg8uXLGDlyJHbv3o0ZM2ZIfq+IiAjY2tpqbefn5yf9ByCdOdo4YnLQZI1LSp24dQK3i2+jq2NXGSsjIiJqvxhQWyk8PBwLFy7E5MmT0b17d7V9d+7cwTvvvIMtW7ZAEATU1NTg0UcfRXx8PHr37i2p/x07djTpl0xjdshsrWuefpfwHZ4b+pw8BREREbVzXKhfB0qlEtOmTcO5c+dw+vRpLF26tNkQ6eXlhc2bN2Pjxo3ia1VVVXjttddkrJYM5aFeD8HGwkZjmz2X98hUDRERUfsnaR1UfWVkZGDr1q2IiopCRkYGLC0t4ePjgwkTJuDxxx9Hly5djF2CyYSHh+Ps2bMAgE6dOiEvLw92dnZN2jWeJJWcnGz0EVSugyrd9N3Ttd7aNO25NPg58ZILIiKilhhlHdTo6Gj89ttv+O2335CYmCjpmC1btiAoKAhr1qzB4cOHkZCQgPj4eBw8eBAvvPACevfujS+++EKXMtqUqVOnis8rKiqQkpJiumKo1Wb3na21zXcJ38lQCRERUfsnOaCqVCpMmjQJU6ZMwZQpU5CWlqb1mE2bNmHFihWoqKgQF+JXKBRqyy4VFxdjyZIl2LFjRyvKN39ubm5q28XFxSaqhPTxYK8H0cmyk8Y2uy7tkqkaIiKi9k1yQI2OjkZhYSEEQUCvXr0wbtw4je1v3LiBVatWAfgzlAqCAHt7e3Tp0kXcrv/vihUrkJmZqd9PY4Yaj5h6enqaphDSi4ONA+7vqfm2pzEZMbiayzVpiYiI9CU5oJ4+fRpAXdh85JFHtLZ/++23UVNTI2736tULR44cQXFxMTIyMpCbm4s1a9ZAqVRCoVCgtLQU//73v1vxI5gvQRDw3Xd/nvb19vaWvBj/iy++iP79+8PZ2Rk2Njbw8vJCWFgYnnvuORw9etRYJZMGUk7z74zfKUMlRERE7ZvkgHr+/HnxecPrKptTXFyM77//XhwddXFxwR9//IHRo0eLbZydnfHGG2/g/fffhyAIEAQB3377bSt+BPO1a9cu3Lx5U9x+/PHHJd8Wdu/evbh06RKKiopQVVWFrKwsxMTEYOPGjRgzZgzCw8Nx6dIlY5VOzXio90PobN1ZY5ud8TuhElQyVURERNQ+SQ6oN27cAADY2Njgnnvu0dj2999/R0VFBYC6Edfnn38e3t7ezbZdsWIFgoKCAADp6elISkqSWpJZu337Np599llx29nZGa+88ork493c3DBkyBCMHz8e4eHhcHd3V9t/9uxZDBkyBAcOHJDcZ2VlJYqLi9UeJJ2dlR1m9pmpsU1qUSqOpx2XqSIiIqL2SXJATU1NhUKhQM+ePWFhYaGxbVRUFACIE6Pmz5/fcgFKJWbO/PNL/+LFi1JLMltlZWWYMWMG8vLyxNc++eQTuLq6ajwuNDQUW7ZsQVJSEnJzc3HmzBkcOnQIp0+fRk5ODmJiYjBr1iyxfXl5OebMmYPY2FhJda1fvx5OTk7ig3eh0t0TA57Q2uaruK9kqISIiKj9khxQi4qKADSdld6cM2fOiM979eqFrl013wKy4YhsW58oVVNTg7lz5yI6Olp87ZlnnsHs2ZqvX+zevTtiYmKwbNmyFq9TDQ0Nxd69e7Fp0ybxtbKyMqxYIe0+8K+88gqKiorEx61btyQdR38a030MfB18NbbZm7AXFTUVMlVERETU/kgOqJWVlQCg9RpKlUqF+Ph4ceZ+eHi41r4bLtRfUlIitSSzo1Kp8MQTT2D//v3ia7Nnz1a7o5QhLF++HIsXLxa3o6KiEBMTo/U4GxsbODo6qj1INxZKCzze/3GNbYoqi3DguvRLL4iIiEid5IBaH2Zyc3M1trt06RLKy8vF7UGDBmntu+ElAw1n/rclKpUKCxcuVJvoNXPmTHz99ddaL4lojca3Tf31118N/h7UvHkD5mltsz12uwyVEBERtU+SA6qHhwcEQcD169fF0dTmNL7+dMiQIVr7zs/PF5937qx5lrQ5UqlUWLJkCb766s9rD6dPn45vvvkGlpaWRnnPwMBA+Pv7i9tXr3L9Tbn079IfA7sM1Njm18Rfcbv4tkwVERERtS+SA2r9SGhVVZXaKezG/vvf/4rPO3XqhLCwMK19N7wrVVtbyF6lUuEvf/mL2u1aH374YezevRtWVlZGfe+GKyNoG9kmw9I2WUolqLD9AkdRiYiIWkNyQG1456jVq1c3e61oTEwMjhw5Il5/OnHiREkh7ezZs+Lz+iWn2oL6cLp9+59B5OGHH8aePXuMHk6BuglS9WxtbY3+fvSnxwc8Dkul5tHxrRe2ck1UIiKiVpAcUOfOnQs7OzsAdWuijh49GocPH0Z5eTnKysrwww8/YPr06QD+PL2/YMECSX3X3xnJwsICffv21ekHMBVTh9PKykokJiaK215eXkZ/T/qTV2cvPNTrIY1tUotScSjpkEwVERERtR86TZJ64403xPAZGxuL++67D507d4aDgwNmzJiB9PR0AHUz/UNCQsTAqsmxY8dw69YtKBQKDBw4sE2MBAqC0CScTp8+XbZwCgD79u1TG0EdOXKkLO9Lf/rL4L9obfP5+c9lqISIiKh9kRxQAWDVqlWYNm0aBEEQb2Pa8FH/mo2NDXbs2CGpz23btonPG15GYK4EQcCTTz6pFk5nzJghyzWn9bKzs/Hyyy+L2/b29pg8ebIs701/mtRjEro6al7j939X/4ecuzkyVURERNQ+6BRQlUol9u7di5deeqnZ2emCIKB79+747bfftN4OFQBSUlKwa9cucXvq1Km6lCM7QRCwdOlSbN26VXxt5syZeofTU6dOYenSpbh27ZrWthcvXsTYsWPVJpatXLlS0g0UyLAslBZYPGixxjbVqmp8EfuFPAURERG1Ewqh/py9jrKzs7F//35cvXoVxcXFcHd3x7BhwzBp0iTJYe3YsWP4/fffAdRdf7pmzZrWlCKbPXv2YM6cOeK2QqHAuHHjdFpKauXKlZg4caLaa5GRkRg7diwAYODAgRg3bhwGDBgALy8vODg4oLS0FImJiTh48CB++uknqFR/TryZOHEiDhw4AGtra51/nuLiYjg5OaGoqIiL9rdSamEqAjYGQEDL/xt1d+6OxOWJsFAafj1cIiKitkRq9mj1Ip2enp5YsmRJaw8HAERERCAiIkKvPuTU8JpPoG5E9fDhwzr1MXfuXI374+LiEBcXJ6mvBQsWYMuWLa0Kp2QY/s7+uK/HfTh482CLbVIKU3Dg+gFMC54mY2VERERtl06n+Mk4AgICMG/ePElLbFlaWmLq1Kk4fPgwvvjiC9jb28tQIWny5OAntbbZfHazDJUQERG1D60+xU/GkZubi/j4eOTk5CAvLw8FBQXo1KkTXFxc0LNnT9x7770GW+mAp/gNo7q2GoGbArXeOerS05cQ4hkiU1VERETmx+in+Mk43N3d28RqBvQnKwsrLLt3GV7941WN7Taf3YyPH/xYpqqIiIjaLp7iJzKAJ0OfhI2FjcY2X8V/hYLyApkqIiIiarsYUIkMwN3OHY/1f0xjm7LqMnx67lOZKiIiImq7GFCJDGT5kOVa23x45kNU1FTIUA0REVHbxYBKZCD3eN+DUd1GaWxzp/QOF+4nIiLSggGVyIBWhK/Q2uafJ/6JGlWNDNUQERG1TQyoRAY0PXg6glw1r2ebXJiMvZf3ylQRERFR28OASmRAFkoLvDj8Ra3t3j3xLrgEMRERUfMYUIkMbP7A+fDu7K2xTXxWPA5cPyBTRURERG0LAyqRgdlY2uD5Yc9rbff6kdehElQyVERERNS2MKASGcFfQ/8Kl04uGtvEZcXhu4TvZKqIiIio7WBAJTICBxsH/G3I37S2e+PIG5zRT0RE1AgDKpGRPBv+LBysHTS2uZZ3DTvjd8pUERERUdvAgEpkJG52bpKuRX3jyBsoqy6ToSIiIqK2gQGVyIieH/Y8XG1dNba5VXwLH5z8QKaKiIiIzB8DKpEROdo44qURL2lt996J93C7+LYMFREREZk/BlQiI/vbkL/Bq7OXxjZl1WV49fCrMlVERERk3hhQiYzMzsoOa8eu1druq/ivcDztuAwVERERmTcGVCIZLBq0CAO7DNTa7qn9T6GyplKGioiIiMwXAyqRDCyUFvj3pH9rbXcl9wreO/GeDBURERGZLwZUIpmMDRiL6cHTtbZ7J+odXM29KkNFRERE5okBlUhG7098H50sO2lsU1Vbhfn/nY/q2mqZqiIiIjIvDKhEMurh2gNvRLyhtV10RjTePvq2DBURERGZHwZUIpmtGr4K/Tz7aW33j+P/QFRqlAwVERERmRcGVCKZWVlY4bOHPoMCCo3tVIIKj+97HNl3s2WqjIiIyDwwoBKZwNCuQ7F8yHKt7W4V38LsvbN5PSoREXUoDKhEJvLuhHfRx72P1nZHU49i5W8rZaiIiIjIPDCgEpmIrZUtds3cBWsLa61tN5/djI2nN8pQFRERkekxoBKZ0CCvQVg/fr2kts8dfA4743cauSIiIiLTY0AlMrHnhj6HB3o+IKntwv8txH+v/NfIFREREZkWAyqRiSkVSuycsRO93HppbVsr1OKRvY/gq7ivZKiMiIjINBhQicyAcydn/G/O/+Bg7aC1ba1Qi/n/m48NpzZAEAQZqiMiIpIXAyqRmejj0Qe7Zu6CUiHtf8uVv63E/P/NR1l1mZErIyIikhcDKpEZebDXg/j8oc8lt98ZvxNDPx+K2DuxxiuKiIhIZgyoRGZm0T2L8P7E9yW3v5h9EWGfheH1P17naCoREbULDKhEZmjV8FV4I+INye1rVDVYF7UOPTf3xGfnPuOdp4iIqE1jQCUyU2+NfQv/uu9fOh2TUZKBpw48hYCNAVgftR7Zd7ONVB0REZHxKAROA+6wiouL4eTkhKKiIjg6Opq6HGrBtgvb8OT+J6ESVDofq1QoMT5gPGb2mYmJPSYi0CXQCBUSERFJIzV7MKB2YAyobcfhpMOY+/1c5Jbl6tVPgHMAhvsNx2DvwbjH6x708+wHdzt3KBQKA1VKRETUMgZU0ooBtW1JK0rDzD0zEZMRY9B+HW0cEegSiECXQAQ4B8C7sze8Hbzh3dkbXp294O3gDScbJ4ZYIiLSGwMqacWA2vZU1lTi7aNv470T76FWqJXtfTtZdoJXZy/4OPjAx8EHvg6+4n99HX3FbXtre9lqIiKitocBlbRiQG27zmeex+IfFiMuK87UpahxsnFCgEsAglyD0MOlh/jfHq490NWxq+SbEBARUfvEgEpaMaC2bTWqGnwZ9yXejHwTt4tvm7ocrWwtbdHXoy/6d+mP/p790c+zH/p79odXZy9ePkBE1EEwoJJWDKjtQ3l1OT459wk2ndmE5MJkU5ejMzdbNwz2HowwnzCE+YZhiO8Q+Dj4mLosIiIyAgZU0ooBtX2pVdVi//X9+PTcp/g96XfUqGpMXVKr+Tj4IMynLqwO9xuOcN9w2FrZmrosIiLSEwMqacWA2n7lluXi+4Tv8dONnxCZEomSqhJTl6QXK6UV7vW5F6O6jcIo/1EY4TcCLrYupi6LiIh0xIBKWjGgdgzVtdWIyYjBsdRjuHDnAi7cuYDreddNXZZeFFCgn2c/jOo2CmMDxmJs97Fws3MzdVlERKQFAyppxYDacZVUluBG/g0kFSThZv5NJBUkIakwCZklmcgszUR+eb6pS9SJAgrc430PxgeMx/iA8RjlPwp2VnamLouIiBphQCWtGFCpJZU1lci6myUG1vr/ZpRkIL0kHenF6cgoyUBeeZ6pS22WtYU1hnUdhgmBEzA+YDzCfMNgqbQ0dVlERB0eAyppxYBK+qqoqUBmSaYYWpMLk3Ez/yYSCxJxM/8mbhXfMnWJAAAHaweM6T4GEwInYELgBPRx78OlrYiITIABlbRiQCVjq6ipQHJBMq7mXsXF7Iu4lH0JF7Mv4kbeDVnvhNWYj4NPXVgNqAus3g7eJquFiKgjYUAlrRhQyVQqaipwLfcaYu/EIiYjBmczziL2TiyqaqtMUk+IRwgmBE7AxMCJiPCPgIONg0nqICJq7xhQSSsGVDInVbVViM+KR3R6NKIzonHy1klcy7smex2WSksM7ToUEwMnYkLgBIT5hMHKwkr2OoiI2iMGVNKKAZXMXfbdbBxPO46o1ChEpUXhwp0LUAkqWWuov361PrAGuwfz+lUiolZiQCWtGFCprSmpLMGp26dwLPUYjqQcwZnbZ2S/lrX++tWJgRMxPmA8r18lItIBAyppxYBKbV1JZQmOpR7D4eTDOJR0CBezL8peQ4hHiDi6yutXiYg0Y0AlrRhQqb3JvpuNP5L/wOGkwziUfAgphSmyvr+l0hKh3qFqt2TlHa6IiP7EgEpaMaBSe5dUkIRDSYdwKOkQDicfNskdsvp69MWobqMwsttIjOo2Cv7O/rLXQERkLhhQSSsGVOpIVIIKsXdi8fvN33Eo+RCiUqNQWVspex1+jn4Y2W0khnYdiiG+QzDIaxA6WXaSvQ4iIlNgQCWtGFCpIyuvLsfJWyfxe9LvOJR0COczz0OA/H8dWiotMbDLQAzxHYIhvkMQ7huO3u69oVQoZa+FiMjYGFBJKwZUoj/lleXhj+Q/6i4JSD6EpIIkk9XiYO2AUJ9QDOoyCAO9BmKQ1yD0ce8DG0sbk9VERGQIDKjtQFpaGr766ivs378faWlpyMvLg6enJwIDAzFjxgw89thj8PDwaHX/DKhELTOH61cbslRaoo97Hwz0GoiBXeoefT36wsfBh+uyElGbwYDaxm3evBkvvfQSysvLW2zj6uqKjz76CLNnz27VezCgEkmjElS4kHlBHF011fWrzXGwdkCwe7D46OPeB8Huwejh2gPWFtamLo+ISA0Dahv2+uuvY926dWqvBQUFwdfXF+np6UhMTFTbt23bNixatEjn92FAJWqd8upynLh1QhxhNdX1q5pYKCzg7+yPQJdABDgHIMA5oO65S91/3WzdOPJKRLJjQG2j9u3bh5kzZ4rbffv2xZdffonQ0FDxtXPnzmH+/PlISEgAAFhaWuLEiRMYMmSITu/FgEpkGLlluWq3ZD2feV72O1zpqrN1Z/g7+cPX0Rc+Dj7wdWj0X0dfeNp7wlJpaepSiagdYUBtg6qrqxEcHIykpLrJGV27dkVcXBxcXV2btM3Pz8eAAQOQnp4OABg1ahSOHTum0/sxoBIZx92quzh9+zSi0qJwPO04Tt0+hbLqMlOXpTMFFHC1dYW7nTvc7NzgbucOd1t3tW1XW1c42jjCwdoBDjYO4nN7a3uuREBETTCgtkE7duzAwoULxe09e/bgkUceabH9nj17MGfOHHH76NGjiIiIkPx+DKhE8qiurcaFOxcQlRqFM+lncDb9LFKLUk1dllEpoEBn685wsHGAg7UDbK1s0cmyE2wsbNDJspPao+FrNpY2sFRawkJhAUulZd1zZYPnDV5vbp+F0gIWCgsoFUooFUpYKBs81+H11vTBSyaItJOaPXjuxozs3btXfO7j44Pp06drbD9jxgz4+PggIyNDPF6XgEpE8rCysBLXOa13p/QOotOjcTb9LM5mnMXZ9LMorCg0XZEGJkBASVUJSqpKTF2KrFobbjW9bog+1F6HmdZloNd1aauAgv+wMFMcQTUT5eXlcHFxQWVl3czgxYsXY+vWrVqPW7x4MbZv3w4A6NatG1JTpY/KcASVyHwIgoDE/EScST+DC5kXEJcVh7isOOSW5Zq6NKIOT4GmIba5YGuKdsZ+z+9mf4f7etzXpF1rcQS1jUlISBDDKQCMGDFC0nEjRowQA2paWhry8/ObvWaViMybQqFAT7ee6OnWE/MGzANQF1ozSzMReycWcXfixNB6Pe86VILKxBUTdRzNrdIheXyvjQ8D1qpMM+GTAdVM1M/Ir9ezZ09JxzVul5CQgJEjRxqsLiIyHYVCAR8HH/g4+OD+nveLr1fUVOBG3g1czb2KK7lXcDX3qvgor2l57WQioraCAdVMpKSkqG1369ZN0nH+/v5N+mkpoFZWVqqN0hYXF+tWJBGZhU6WndC/S3/079Jf7XWVoMKtolu4mnsVNwtuIqkgCcmFyUguSEZSQRKKKotMVDERkW4YUM1E47Do7Ows6TgnJye17ZKSlickrF+/Hm+99ZbOtRFR26BUKOHv7A9/Z/9m9xeUF4iBNaUwBekl6cgoyfjzv8XpZnOHLCLq2BhQzURpaanatq2traTjGrfTFFBfeeUVPP/88+J2cXEx/Pz8dKiSiNoyF1sXuNi6YLD34Gb3C4KAgooCpBfXBdacshzkluUitywXeWV5yC3PVd8uy0W1qlrmn4KI5GSqVQ4YUM1EdbX6X/KWltL+aKysrDT205CNjQ1sbGx0L46IOgSFom5hfldb1yaXDzRHEARU1FSgpKoExZXFKKn8//+tKmnyvLK2EpU1laioqUBFbUXdf2sq/nytfru2ErWqWtQKtahR1aBGVYNaVYPnjV439zt2EVHrMKCaCXt7e7XtiooK2NnZaT2uvFx9QkTjfoiIjEWhUMDWyha2VrbwtPc0SQ2CIKBWqEWtqhbVqmqoBBVUggq1qto/nwu1srxuivcUX4e0Y2Wvq9HrRFIxoJqJzp07q22XlZVJCqhlZeq3T3RwcDBoXURE5kyhUMBSUXcnKRvwDJG5q/8HhbkE/cahWepyUqZoZ6ra+nn2a/KaHBhQzYSHh4fadmZmJtzd3bUel5mZqbYt5RgiIiJTqP8HBZE2SlMXQHWCg4PVtqXeEapxu8b9EBEREbU1DKhmIiQkRG37/Pnzko5r2M7a2hpBQUEGrYuIiIhIbgyoZsLPzw+BgYHi9tGjRyUd17DdyJEjYWFhYfDaiIiIiOTEgGpGpk+fLj6PjIxEWlqaxvZpaWlqAXXmzJlGq42IiIhILgyoZmTRokVQKuv+SFQqFdauXaux/dtvvw2Vqm4GYufOnTF79myj10hERERkbAyoZiQkJATz5s0Ttz///HN89tlnzbb99NNPsXXrVnF71apVnMFPRERE7YJCaG7RKzKZ7OxsDB06FMnJyeJrDz30EObOnQsfHx9kZGTg22+/xf79+8X9YWFhiIyMlLRuakPFxcVwcnJCUVERHB0dDfYzEBERETVHavZgQDVDV69exaRJk7RegwoAAwYMwO+//w5PT93v4sKASkRERHKSmj14it8MBQcH4+LFi1i2bFmTO0zVc3Nzw+rVqxEdHd2qcEpERERkrjiCaubKy8sRGRmJ1NRUFBQUwN3dHT169EBERAQsLfW7GwdHUImIiEhOUrMH7zdm5mxtbTFlyhRTl0FEREQkGwbUDqx+8Ly4uNjElRAREVFHUJ85tJ3AZ0DtwEpKSgDU3cWKiIiISC4lJSVwcnJqcT+vQe3AVCoVMjIy4ODgAIVCYZT3KC4uhp+fH27dusXrXKnD4O89dTT8nSepBEFASUkJfHx8xJsTNYcjqB2YUqlE165dZXkvR0dH/qVFHQ5/76mj4e88SaFp5LQel5kiIiIiIrPCgEpEREREZoUBlYzKxsYGb775JmxsbExdCpFs+HtPHQ1/58nQOEmKiIiIiMwKR1CJiIiIyKwwoBIRERGRWWFAJSIiIiKzwoBKRERERGaFAZUMLi0tDe+88w6GDh0KHx8f2NjYwM/PD6NHj8bGjRuRk5Nj6hKJ9JaSkgKFQqHzIzg42NSlE6kpKSlBVFQUPvzwQ8yfPx8hISGwtLQUf2e7d++u93vwe4F0xVn8ZFCbN2/GSy+9hPLy8hbbuLq64qOPPsLs2bNlrIzIsFJSUhAQEKDzcb1798bVq1eNUBGR7nr37o0bN25AUxTw9/dHSkpKq9+D3wvUGrzVKRnM66+/jnXr1qm9FhQUBF9fX6SnpyMxMREAkJ+fjzlz5uDu3btYtGiRKUolMriIiAjY2tpqbefn5ydDNUTSXL9+3aj983uBWosjqGQQ+/btw8yZM8Xtvn374ssvv0RoaKj42rlz5zB//nwkJCQAACwtLXHixAkMGTJE9nqJ9NV4BDU5Odkgp0KJ5KRQKAAA9vb2GDhwIAYPHozQ0FDs3r0bv/76K4DWj6Dye4H0wYBKequurkZwcDCSkpIAAF27dkVcXBxcXV2btM3Pz8eAAQOQnp4OABg1ahSOHTsma71EhsCASu3B119/jcGDB6N3795QKv+clrJw4ULs2LEDQOsCKr8XSF+cJEV627Vrl/iXEABs2LCh2b+EgLrrjDZs2CBuR0VF8S8iIiITefzxx9GnTx+1cGoI/F4gfTGgkt727t0rPvfx8cH06dM1tp8xYwZ8fHyaPZ6IiNo+fi+QvhhQSS/l5eU4dOiQuD158mRYWmqee2dpaYlJkyaJ2z/++KPR6iMiInnxe4EMgQGV9JKQkIDKykpxe8SIEZKOa9guLS0N+fn5Bq+NiIjkx+8FMgQGVNJL/czLej179pR0XON2jfshamtefPFF9O/fH87OzrCxsYGXlxfCwsLw3HPP4ejRo6Yuj0g2/F4gQ2BAJb00ntnZrVs3Scf5+/tr7Ieordm7dy8uXbqEoqIiVFVVISsrCzExMdi4cSPGjBmD8PBwXLp0ydRlEhkdvxfIEBhQSS/FxcVq287OzpKOc3JyUtsuKSkxVElEJuHm5oYhQ4Zg/PjxCA8Ph7u7u9r+s2fPYsiQIThw4ICJKiSSB78XyBAYUEkvpaWlattS7qTTXDv+RURtUWhoKLZs2YKkpCTk5ubizJkzOHToEE6fPo2cnBzExMRg1qxZYvvy8nLMmTMHsbGxpiuayMj4vUCGwIBKeqmurlbb1jZTs56VlZXGfojMXffu3RETE4Nly5apLdjfUGhoKPbu3YtNmzaJr5WVlWHFihVylUkkO34vkCEwoJJe7O3t1bYrKiokHVdeXq6xH6L2ZPny5Vi8eLG4HRUVhZiYGBNWRGQ8/F4gQ2BAJb107txZbbusrEzScY3bOTg4GKwmInP02muvqW3X3+ecqL3h9wIZAgMq6cXDw0NtOzMzU9Jxjds1nlBC1N4EBgaqzVK+evWqCashMh5+L5AhMKCSXoKDg9W2U1NTJR3XuF3jfojaI29vb/F5bm6uCSshMh5+L5AhMKCSXkJCQtS2z58/L+m4hu2sra0RFBRk0LqIzFHDU5hSZzYTtTX8XiBDYEAlvfj5+SEwMFDclnrHnIbtRo4cCQsLC4PXRmROKisrkZiYKG57eXmZsBoi4+H3AhkCAyrpbfr06eLzyMhIpKWlaWyflpam9hfRzJkzjVYbkbnYt2+f2gjqyJEjTVgNkXHxe4H0xYBKelu0aBGUyrpfJZVKhbVr12ps//bbb0OlUgGom+05e/Zso9dIZErZ2dl4+eWXxW17e3tMnjzZhBURGRe/F0hfDKikt5CQEMybN0/c/vzzz/HZZ5812/bTTz/F1q1bxe1Vq1Zxpia1OadOncLSpUtx7do1rW0vXryIsWPHqo0grVy5Em5ubsYskcik+L1A+lIIgiCYughq+7KzszF06FAkJyeLrz300EOYO3cufHx8kJGRgW+//Rb79+8X94eFhSEyMhJ2dnamKJmo1SIjIzF27FgAwMCBAzFu3DgMGDAAXl5ecHBwQGlpKRITE3Hw4EH89NNP4sgQAEycOBEHDhyAtbW1qconEq1btw7r1q1r8np1dbXa762NjU2TNk888USLoRPg9wLpR9r9x4i08PT0xM8//4xJkyaJI0X79+9X+4unoQEDBuDAgQP8S4javLi4OMTFxUlqu2DBAmzZsoXhlMxGTU0NKisrtbZrro22W5Hye4H0wVP8ZDDBwcG4ePEili1b1uROIvXc3NywevVqREdHw9PTU+YKiQwjICAA8+bNk7QMjqWlJaZOnYrDhw/jiy++4O0bqUPh9wK1Fk/xk1GUl5cjMjISqampKCgogLu7O3r06IGIiAhYWnLgntqP3NxcxMfHIycnB3l5eSgoKECnTp3g4uKCnj174t577+Wap0Tg9wLphgGViIiIiMwKT/ETERERkVlhQCUiIiIis8KASkRERERmhQGViIiIiMwKAyoRERERmRUGVCIiIiIyKwyoRERERGRWGFCJiIiIyKwwoBIRERGRWWFAJSIiIiKzwoBKRG1WSkoKFAqF+EhJSTF1SW3CwoULxc9s4cKFpi6nTYuNjYWlpSUUCgVmzpxp6nJklZ+fDxcXFygUCvTv3x81NTWmLonaEQZUImpizZo1asGvNQ+GxdbT97Nv6REZGWnqH61dEQQBzzzzDGpra2FpaYn169ebuiRZubq64uWXXwYAXLp0Cf/3f/9n4oqoPWFAJSIiaoU9e/bg5MmTAIAnnngCvXr1MnFF8lu+fDk8PDwAAG+99RYKCwtNWxC1G5amLoCIzN+kSZN0PsbW1tYIlXQMUj7vixcvIiMjQ6djXF1d9aqL/lRbW4s33ngDAKBUKsWRxI7Gzs4Ozz77LFavXo3CwkL885//xD/+8Q9Tl0XtgEIQBMHURRCReVmzZg3eeustcdtc/5pISUlBQECAuJ2cnIzu3bubriAZLVy4EDt27BC3zfXPqL3auXMnnnjiCQDAgw8+iP3795u4ItPJz8+Hj48PKisr0blzZ9y6dQvOzs6mLovaOJ7iJyIi0tGGDRvE50uXLjVhJabn6uqKWbNmAQBKS0vx6aefmrgiag8YUImIiHRw8uRJXLhwAQDg4eGByZMnm7gi01uwYIH4/D//+Q9H9ElvDKhEZFKHDx/Gc889h8GDB8PLywvW1tbo3LkzgoKCMG3aNGzYsEHnFQHKy8uxfft2TJgwAX5+frCxsYGnpyfGjBmDjRs3oqKiQuc6L126hFdffRVDhgwR6/Tw8MCgQYPw/PPP48yZMzr3aSpSl5kaM2aM2G7NmjXi6wcPHsSjjz6Knj17ws7ODg4ODggLC8MHH3yA8vLyJv1UVlbi008/xYQJE+Dh4QErKyt06dIFU6ZMwXfffdeqnyEhIQFvvvkmRowYAV9fX9jY2MDFxQV9+/bF008/jaNHj7aqXykaXloxY8YMWFhYaD2me/fu4mf5xRdfAKi7jvW7777DtGnTEBgYiE6dOsHZ2RmjRo3CZ599htra2ib9lJSU4F//+hdGjhwJFxcXWFtbw8fHBzNnzsThw4d1+jmOHDmCv/71r7jnnnvg6uoKKysrdOrUCR4eHhg8eDAee+wxbNq0CampqVr7Gjt2LNzc3AAAqampRv38qYMQiIgaefPNNwUA4sMYYmJihPDwcLX3aemhUCiEffv2NekjOTlZrV1ycrIQFxcn9OnTR2N/QUFBws2bNyXVeffuXWHJkiWCUqnUWueMGTOEvLw8Q39UzVqwYEGr/4waHrtgwYIW240ePVps9+abbwqFhYXCrFmzNH4GAwYMEHJzc8U+Ll68KPTu3VvjMbNnzxaqq6sl1V5QUCAsWLBA0p/HlClThJycHJ0+G21qa2sFNzc38T1++OEHScf5+/uLx2zfvl1IT09X+3ybe4wfP14oKysT+zh69Kjg7e2t8Zjnn39eay3Z2dnChAkTJP2/V/8oLS3V2u+jjz4qtl+6dKmkz4WoJZzFT0Sy++677zB//ny10TaFQoHevXvDy8sL1dXVSE9PF0dOBUFAUVGR1n6vXbuGuXPnikvd9OzZE76+vigrK0NcXBwqKysBAImJibjvvvsQHx8POzu7FvsrLi7GpEmTcPr0abXXe/fuDW9vb+Tl5eHy5ctQqVQAgH379uHKlSs4dOgQfHx8dPlIzF5tbS1mzJiBP/74AwDg4+ODoKAgVFdXIy4uDmVlZQCA+Ph4TJkyBWfOnMGNGzcwevRo5OfnAwB69eoFX19fFBYWIj4+Xhwh3LNnD7p164b3339fYw1paWmYPHkyrly5Ir6mVCoRHByMLl26oLy8HJcuXUJpaSkA4JdffsHw4cNx7NgxeHl5GeRzOH36NPLy8sT3joiI0LmP0tJS3Hfffbh8+TIAICAgAP7+/igrK0NsbCyqqqoA1J1deOKJJ/Ddd98hKioKkyZNQkVFBRQKBUJCQuDp6Yns7GxcunRJ7HvDhg0IDAzEM8880+x7V1ZWYsKECYiPjxdfs7KyQq9eveDh4QGlUonCwkIkJiaiuLhYbCNIOGU/duxYfPPNNwCAn3/+WefPhUiNqRMyEZkfY46gnj59WrC2thb7trOzE9566y0hOzu7Sdvs7Gzh008/FQYPHixs3769yf7GI6ju7u4CAGHmzJnCjRs31NoWFhYKixcvVmu/bt06jbU+9thjau0feOCBJv1mZGQI8+fPV2s3ZswYoba2VvcPRwdyj6DWjxr27dtXOHLkiFq70tJS4a9//ataPTt37hT69+8vABAeeuihJp9bSkqKMGzYMLG9paWlkJSU1GItFRUVwqBBg8T2nTp1EtatW9dkxLqqqkrYunWr4OjoqDYSaag/j3Xr1on9BgcHSz6u4Qhq/Wc5fPhw4fz582rtcnNzhWnTpql9lr/++qs4crp48WIhIyND7ZiLFy8KvXr1Ets7OjoKJSUlzdaxadMmsZ2VlZXw3nvvCcXFxc22vXLlivDee+8JPXr0aLG/huLi4tTqbvxnTqQLBlQiasJYAbW6ulro0aOH2K+zs7Nw5swZScc29wXZOKACEJ5++mmN/TQ8tRkQENBiu99//73JaWhNIefvf/+7Wvtt27ZJ+rlaS+6ACkDo06ePkJ+f32L7kSNHim3r/xHy2GOPCSqVqtn22dnZakHyrbfearHvl19+WWzXuXNnrb83Fy5cEOzs7MRjvv/+e43tpWoYHmfPni35uIYBFYAwevRooaKiotm2lZWVQmBgYJPP8qWXXmqx/4SEBMHCwkI8ZseOHc22Gz9+vNjm1VdflVR7bW1ti3+GDVVVVan943PXrl2S+idqDidJEZFWut5Ss6WJN99++y1u3rwpbm/ZsgVDhgyRVEPnzp21tunevTs+/PBDjW1eeOEF8XlycrLaYvcNbdy4UXzu4eGBTz75BEply39lvvfee+jdu7e4ra2Otujjjz+Gi4tLi/uffvpp8XlVVRWcnZ3x0UcfQaFQNNvew8MDs2fPFrePHz/ebLvi4mL85z//Ebf/+c9/av29GTRokNri+Zs2bdLYXqqEhATxeVBQUKv6sLCwwPbt22FjY9PsfmtrayxZskTcrqqqQq9evbBu3boW++zTpw/Gjh0rbrf0Wd66dUt8PmrUKEn1KpXKFv8MG7KyskK3bt3E7YaXHhDpigGViGTz9ddfi8+Dg4Px2GOPGbT/J598EtbW1hrbjBw5Ui1o1l8H2FB5eTl+/fVXcXvJkiVaFx63srLCihUrxO34+Hi1MN7WBQcHa73ecujQoWrbc+fOhaOjo+RjGoa/hvbt2ydeD+ns7KwW3jRp+A+lkydPitfJtpYgCGoz2n19fVvVz8SJE9VuMNGcxp/l4sWLYWmpedqIlM+y4R3e6pfKMqSuXbuKz3VdfYOoIU6SIiKtdL3Vaf/+/Zu8VlNTgxMnTojbjzzyiN51NTZixAitbezs7ODq6orc3FwAaPbe4dHR0aipqRG3H3roIUnvP23aNLXJKSdPnkSPHj0kHWvuhg0bprVN44lIUo7x9vYWnxcUFDTbpuGSRREREVr/EVLPz88Pzs7OKCwsFCdzSampJUVFReIEJqD1t4415WcZHh6OuLg4AMBbb70FNzc3zJ8/H506ddLavxQNR9izs7MN0id1TAyoRKRVw9HE1rp9+zZKSkrE7bCwML37bEzqTG17e3sxoDY3qnbjxg217YEDB0rq19fXF25ubuIs78TEREnHtQVSPtvGKyLoekxLI5wNZ5xfuHBBp4XxG655m5OTI/m45ty9e1dtu+FopC5M+VmuWLECO3bsQGVlJSorK/HXv/4Vq1atwsSJEzF69GgMHz4cgwYN0jpaK6WGxp8XkS4YUIlIFvVLDdXz9PQ0+Hu0dE2fJkIzy+c0HH2ytbWFvb295P7c3d3FgNrSKFZbJHXUUt9jmlP/eQJ111A2vI5SF82Nluujud8dKeT4LFuqLSQkBHv27MG8efPEfzCWlJRg37592LdvHwDAwcEBY8eOxWOPPYaZM2fqFFZb+5kQNcZrUIlIFo3v3tSaMCmX+vVSAd2DQcOfqzV3rKKmDDUSV79ebWs1/odKc3fNagumTp2KxMREvPLKK/Dz82uyv6SkBD/++CPmzp2LPn364NSpU5L7bviZ6PIPO6LGGFCJSBaNZ39LWXjfVBpOiCotLdVpVKjh4ubaJlaRNA0/x9deew1C3RKJOj803dZVCicnJ7V/gDQ+K9CWeHp64h//+AfS0tJw7do1bN26FYsWLWoyeSsxMRETJkyQPKGq4WdijLMk1HEwoBKRLBpO4ADq7vpkrhp+sdbW1kq6FzlQN/Kanp7ebD/Ueg2vv7x+/brJ6lAoFPD39xe3b9++bbJaDKlXr15YvHgxtm3bhqSkJMTGxmL+/Pni/rKyMrUluzRp+Jl0797d0KVSB8KASkSycHZ2Rp8+fcTtyMhI0xWjRWhoqNr2yZMnJR0XHR2N6upqcfvee+81aF0d1fDhw8Xnhw8fFm+Ragr9+vUTn7enSXANDRw4EDt27MDjjz8uvnbkyBG1S1+aU1VVhbS0NHG7udU8iKRiQCUi2TScfb1v3z5kZWWZsJqWBQYGql2b99VXX0k6bseOHeJzW1tbyTchIM3uv/9+8Xl+fj527dplsloarj7RcHWB9mjmzJni8+rqaq2XNFy5ckXtH2jGWKmDOg4GVCKSzfLly2FhYQEA4hI35jrr96mnnhKf//rrr/j99981to+Pj8cXX3whbj/xxBNNlgqi1hk9erRa2H/xxRdNdnq94d2abty40eauQ9Xl/7fS0lK1bU13EQOAM2fOiM/9/PzazRrAZBoMqEQkm4CAACxfvlzc/uGHH/D4449rnDAlCAL279/f4q0bjeXpp59Wu4Z0zpw5al/ADV27dg0PPviguLi/nZ0dXnzxRVnq7Cj+/e9/w8rKCgBw584djB49GtHR0VqPS0tLw+rVq7Fy5UqD1BEWFgYPDw8Adb+bDW8i0Bb0798f27ZtaxI+G8vJycH69evF7WHDhmldzP/IkSPi84aj3kStwXVQiUgrXRZGr7d06VI8/PDDTV7/5z//iejoaPGuUt988w1+++03PPbYYxg1ahS6dOmC6upq3L59G9HR0fjhhx9w+/ZtbN++HSNHjtT3R5HMzc0N27dvx0MPPQSVSoWCggKMGDECjz/+OB544AF4eXkhLy8Phw4dwrZt29SWlNq4cSNHjwxs+PDh2LJlizjqnpSUhCFDhmDixImYMmUKgoOD4eDggLt37yIrKwvx8fE4evQoYmJiAAALFiwwSB1KpRIzZ87Exx9/DAD46aefMH36dIP0LYfLly9jyZIleOaZZzBx4kSEh4ejT58+cHNzg1KpRFZWFk6dOoUvvvhCbXR4zZo1GvutqalRO8swe/ZsY/0I1EEwoBKRVgcPHtT5mJZCrZWVFX777TfMmzcP//3vfwHULcS+efNmbN68Wa86De3+++/H7t27MW/ePFRWVqK2thZffvklvvzyy2bbK5VKbN68GX/5y19krrRjePLJJ+Hu7o4FCxaIi8z//vvvWi+/MLQFCxaIAfV///sfPv7441bfeclUKioqsH//fuzfv19jOwsLC2zatAn33XefxnZHjhwRb6jg5+endikEUWvwFD8Ryc7Ozg779u3D//73P60z3d3c3PDUU0+Z7Atv1qxZuHjxIh555BHxFHNjSqUSkyZNQkxMDJYtWyZzhR3L9OnTcePGDbzwwgtwd3fX2NbGxgbjxo3Dli1bsGHDBoPVMHToUAwePBhA3T+ufvnlF4P1bWxbtmzB/fffD0dHR43trKysMHXqVMm/0w0nEj799NNQKBR610odm0Iw1xkKRNRhZGZm4sSJE7hz5w4KCwtha2sLb29vhISEoH///lAqzePf0iUlJTh69Chu3bqFgoICODo6wtfXF6NHj4arq6upy+twBEFAXFwcLl68iLy8PJSWlsLe3h6enp7o1asX+vfvr/W6ydb6+uuvMW/ePADAAw88gAMHDhjlfYxFpVLhypUruHbtGm7fvo3S0lIoFAo4OzujV69eCA0NlXyjiYKCAvj4+KCiogL29vZIS0vj/w+kNwZUIiIiHdXW1qJv3764fv06lEolEhIS0Lt3b1OXZRLr16/Hq6++CgB46aWX8O6775q4ImoPGFCJiIhaYffu3Zg7dy4AYNGiRdi2bZuJK5JfeXk5unfvjuzsbDg5OSEpKYmjp2QQ5nHejIiIqI2ZPXu2eJerr776yqS3YTWVzZs3Izs7GwDwxhtvMJySwXAElYiIqJViY2Nx7733ora2FjNmzMD3339v6pJkk5+fjx49eqCwsBAhISGIjY1tc6sZkPliQCUiIiIis8JT/ERERERkVhhQiYiIiMisMKASERERkVlhQCUiIiIis8KASkRERERmhQGViIiIiMwKAyoRERERmRUGVCIiIiIyKwyoRERERGRWGFCJiIiIyKwwoBIRERGRWfl/V3COYR1eZDMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Generate some reasonable parameters and signal\n",
    "field_strength = 1.5 # Field strength in Tesla\n",
    "\n",
    "global myAlgoParams \n",
    "global myImDataParams\n",
    "global t\n",
    "global s1\n",
    "\n",
    "# Set some parameters for plotting a continuous-looking CSE signal\n",
    "N = 1001\n",
    "tmax = 13*0.001\n",
    "t = np.linspace(0.0,tmax,N)\n",
    "\n",
    "# Set some true FF and R2* (can try different values but for this \n",
    "# exercise we will end up fixing them at 30.0 and 100.0, respectively)\n",
    "true_FF = 30.0\n",
    "true_R2star = 667\n",
    "\n",
    "myTrueParams = trueParams(100., true_FF, true_R2star, 0.)\n",
    "\n",
    "fFreq = [-3.80, -3.40, -2.60, -1.94, -0.39, 0.60];\n",
    "fAmps = [0.087, 0.693, 0.128, 0.004, 0.039, 0.048];\n",
    "myAlgoParams = algoParams(fFreq, fAmps, True )\n",
    "\n",
    "myImDataParams = imDataParams(field_strength, t) # Field strength is given here\n",
    "s1 = create_signal(myImDataParams, myTrueParams , myAlgoParams);\n",
    "\n",
    "plt.figure(figsize=(7,5))\n",
    "plt.plot(t*1000,abs(s1),'green',linewidth=6.0);\n",
    "plt.xlabel('Echo Time (ms)');\n",
    "plt.ylabel('Signal (AU)');\n",
    "plt.rc('font', **font)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.stats as stats\n",
    "\n",
    "# Plot the least-squares signal fit given some echo time combination\n",
    "def showFitNoise(FF, R2star, TE1,dTE, SNR, Npixels):\n",
    "    global myAlgoParams \n",
    "    global myImDataParams\n",
    "    global t\n",
    "    global s1\n",
    "    \n",
    "    font = {'weight' : 'normal',\n",
    "        'size'   : 16}\n",
    "    \n",
    "    A = 100.\n",
    "    myTrueParams = trueParams(A, FF, R2star, 0.)\n",
    "\n",
    "    \n",
    "    # Generate continuous-looking noiseless true signal\n",
    "    N = 1001\n",
    "    tmax = 20*0.001\n",
    "    t = np.linspace(0.0,tmax,N)\n",
    "    myImDataParams = imDataParams(1.5, t)\n",
    "    s1 = create_signal(myImDataParams, myTrueParams, myAlgoParams );\n",
    "\n",
    "    \n",
    "    TE1 = TE1*0.001\n",
    "    dTE = dTE*0.001\n",
    "    nTE = 6\n",
    "    TE = np.linspace(TE1,TE1 + (nTE-1)*dTE,nTE)\n",
    "    myImDataParams = imDataParams(1.5, TE)   \n",
    "    s = create_signal(myImDataParams, myTrueParams , myAlgoParams);\n",
    "\n",
    "    \n",
    "    sigma = A/SNR;\n",
    "    #Npixels = 1;\n",
    "    \n",
    "    ffest = np.zeros(Npixels)\n",
    "\n",
    "    r2est = np.zeros(Npixels)\n",
    "    \n",
    "    sfit3 = np.zeros((N, Npixels))\n",
    "    snall = np.zeros((nTE, Npixels))\n",
    "    for k in range(Npixels):\n",
    "        sn = s +  stats.norm.rvs(0, sigma, size=(nTE)) + stats.norm.rvs(0, sigma, size=(nTE))*1j\n",
    "        snall[:,k] = np.abs(sn) \n",
    "\n",
    "        p0 = np.array([90.,40.,50.])\n",
    "        \n",
    "        try:\n",
    "            popt, pcov = curve_fit(create_signal_r2star, TE, abs(sn), p0);\n",
    "        except RuntimeError:\n",
    "            popt = p0\n",
    "\n",
    "\n",
    "        sfit2 = create_signal_r2star( TE, popt[0], popt[1], popt[2])\n",
    "        \n",
    "        myEstParams = trueParams(popt[0], popt[1], popt[2], 0.)\n",
    "        mySimParams = imDataParams(1.5, t)\n",
    "        sfit3[:,k] = np.abs(create_signal(mySimParams, myEstParams , myAlgoParams));\n",
    "        ffest[k] = popt[1]\n",
    "        r2est[k] = popt[2]\n",
    "        \n",
    "    \n",
    "    # Display results\n",
    "    fig, axs = plt.subplots(1, 2)\n",
    "    fig.set_size_inches(15, 7)\n",
    "    fig.subplots_adjust(hspace=0.3,wspace=0.3)\n",
    "    p2 = axs[0].plot(t*1000,abs(s1),'gray',linewidth=6.0, label=\"Noiseless signal (continuous)\");\n",
    "    p1 = axs[0].plot(t*1000,abs(sfit3),'#DD010120',linewidth=2.0, label=\"Fitted signal (continuous)\");\n",
    "    p3 = axs[0].plot(TE*1000,abs(snall),'or',markersize=6, label=\"Noisy signal (sampled)\");\n",
    "    p4 = axs[0].plot(TE*1000,abs(s),'oc',markersize=12,label=\"Noiseless signal (sampled)\");\n",
    "\n",
    "    axs[0].set_xlabel('Echo Time (ms)');\n",
    "    axs[0].set_ylabel('Signal (AU)');\n",
    "    axs[0].set_title('Magnitude fitting: \\nMeasured and estimated signals' );\n",
    "\n",
    "    axs[0].legend([\n",
    "           lines.Line2D([0], [0], linewidth=6.0, ls='-', c='gray'),\n",
    "        lines.Line2D([0], [0],marker='o',markersize=12,  c='c',ls=''),\n",
    "        lines.Line2D([0], [0], marker='o',markersize=6,  c='r',ls=''),      \n",
    "        lines.Line2D([0], [0],  linewidth=1.0,ls='-', c='#DD010160')],\n",
    "          ['Noiseless signal (continuous)',\n",
    "           'Noiseless signal (sampled)',\n",
    "           'Noisy signal (sampled)', \n",
    "           'Fitted signal (continuous)'])\n",
    "    \n",
    "    axs[1].boxplot([r2est,ffest]);\n",
    "    axs[1].set_ylim([-10, 100])\n",
    "    axs[1].set_ylabel('R2* and FF estimates');\n",
    "    axs[1].set_title('True R2* = ' + '{:.1f}'.format(R2star) + ', Estimated R2* =' + '{:.1f}'.format(np.mean(r2est)) + '$\\pm$' + '{:.1f}'.format(np.sqrt(np.var(r2est))) + ',\\n True FF ='  + '{:.1f}'.format(FF) + ', Estimated FF =' + '{:.1f}'.format(np.mean(ffest)) + '$\\pm$' + '{:.1f}'.format(np.sqrt(np.var(ffest)))  );\n",
    "    axs[1].set_xticklabels(['R2*', 'FF'])\n",
    "    plt.show()\n",
    "    plt.rc('font', **font)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "fa56ec22298a4885b6f92056430a4249",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=25, description='FF', max=50), IntSlider(value=60, description='R2star',…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Interactive plot to demonstrate noise propagation in \n",
    "## magnitude based CSE-MRI estimates of R2* and PDFF\n",
    "y = interactive(showFitNoise,FF=(0,50,1),\n",
    "                R2star=(20,100,10),\n",
    "                TE1=(0.4,4.0,0.01),\n",
    "                dTE=(0.6,4.0,0.01), \n",
    "                SNR=(2,50,1), \n",
    "                Npixels=fixed(100))\n",
    "display(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Effective number of signal averages (NSA)\n",
    "\n",
    "How can we quantify the noise propagation from our acquired echo images to our resulting output images obtained from CSE-based fitting? One reasonable approach is to compare the noise variance ($\\sigma^2$) at the input (in the echo images) vs the noise variance at the output (say, in the fat-only image). This has the advantage that we can compare to a simple averaging process, where the variance will decrease proportionally to the number of averages $N$ ($\\sigma_{ave}^2 = \\sigma^2/N$). This way, if we know how the variance in the fat-only image ($\\sigma^2_{fat}$) relates to the variance in the echo images ($\\sigma^2$), we can calculate the effective number of signal averages (NSA) that we are obtaining in the water-fat separation process (i.e., what is the equivalent noise performance that we would obtain through simple averaging). \n",
    "\n",
    "$$NSA = \\frac{\\sigma^2}{\\sigma^2_{fat}}$$\n",
    "\n",
    "where the value of NSA will be between $0$ (in a terrible case where we are not able to obtain any reasonable fat images) and $N_{echoes}$ (in an ideal case where we optimally use our echoes to generate high-SNR water-fat separated images). As we show below, NSA is highly dependent on the choice of echo times, and also depends heavily on how we fit our signals (magnitude vs complex fitting, with the latter giving higher NSA). Note that $NSA = N_{echoes}$ is achievable if we do not include $R_2^*$ in the model (as in the original IDEAL echo times), but generally not achievable in the presence of $R_2^*$ decay. \n",
    "\n",
    "In the following cells, we calculate NSA for both magnitude and complex fitting, for a range of echo time combinations indexed by the initial echo time $TE_1$ and the echo spacing $\\Delta TE$. \n",
    "\n",
    "In each case, we calculate a theoretical NSA using the Cramer-Rao Lower Bound (CRLB). The CRLB provides a measure (actually a bound) of the noise propagation achievable in the presence of a given signal model, sampling parameters, and noise distribution. Note that we could also estimate noise propagation in each case using simulations, as in the cells above. However, this would be computationally expensive and potentially algorithm-dependent. In contrast, the CRLB is a good tool to rapidly predict NSA in an algorithm-independent manner.\n",
    "\n",
    "We will not go over the details of the CRLB in this description, but a few good references include:\n",
    "\n",
    "    Scharf and MacWhorter, Geometry of the Cramer-Rao bound, 1992\n",
    "\n",
    "    Pineda, Reeder, Wen, and Pelc, Cramér-Rao bounds for three-point decomposition of water and fat, 2005\n",
    "\n",
    "    Hernando, Liang, and Kellman, Chemical shift–based water/fat separation: A comparison of signal models, 2010"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Water-fat separation signal model (complex signal): \n",
    "\n",
    "\n",
    "<font size=\"4\">\n",
    "\n",
    "<br>\n",
    "\n",
    "$s(TE) = \\left( A_{\\hbox{water}} + A_{\\hbox{fat}} \\sum_{p=1}^P \\rho_p e^{i 2 \\pi TE f_p} \\right) e^{-  TE \\, R_2^* } e^{i 2 \\pi f_B TE + i \\phi } $\n",
    "\n",
    "**Known parameters**: TE values (say, 6 echo times), $\\rho_p$ and $f_p$ (fat NMR peak relative amplitudes and frequencies)\n",
    "\n",
    "\n",
    "**Unknown parameters**: $A_{\\hbox{water}}$, $A_{\\hbox{fat}}, R_2^*, f_B, \\phi$\n",
    "\n",
    "**Fat quantification**: FF = $100 \\cdot A_{\\hbox{fat}} / (A_{\\hbox{fat}} + A_{\\hbox{water}})$\n",
    "\n",
    "</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "def CRLB_CSE_complex(fm,a,fF,relAmps,r2,phi,t,sigma):\n",
    "    N = np.size(t)\n",
    "    \n",
    "    FS, T = np.meshgrid(fF,t)\n",
    "    RA, T = np.meshgrid(relAmps,t)\n",
    "    \n",
    "    D = np.zeros((N,5)) + 1j*np.zeros((N,5)) \n",
    "        \n",
    "    # Derivative wrt a[0]\n",
    "    D[:,0] = np.exp(1j*(phi + 2*np.pi*fm*t) - r2*t)\n",
    "    \n",
    "    # Derivative wrt a[1]\n",
    "    D[:,1] = np.sum(RA*np.exp(1j*2*np.pi*fm*T)*np.exp(1j*phi+ 1j*2*np.pi*FS*T),1)\n",
    "    \n",
    "    # Derivative wrt r2\n",
    "    D[:,2] = -t*np.exp(1j*(phi + 2*np.pi*fm*t) - r2*t)*(a[0] + a[1]*np.sum(RA*np.exp(1j*2*np.pi*FS*T),1))\n",
    "    \n",
    "    # Derivative wrt phi\n",
    "    D[:,3] = 1j*np.exp(1j*(phi + 2*np.pi*fm*t) - r2*t)*(a[0] + a[1]*np.sum(RA*np.exp(1j*2*np.pi*FS*T),1))\n",
    "\n",
    "    # Derivative wrt fm\n",
    "    D[:,4] = 1j*2*np.pi*t*np.exp(1j*(phi + 2*np.pi*fm*t) - r2*t)*(a[0] + a[1]*np.sum(RA*np.exp(1j*2*np.pi*FS*T),1))\n",
    "\n",
    "    D = np.matrix(D)\n",
    "    \n",
    "    C = np.linalg.inv(np.real(np.matmul(D.H,D)))*(sigma**2)\n",
    "\n",
    "    return C \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Water-fat separation signal model (magnitude signal): \n",
    "\n",
    "\n",
    "<font size=\"4\">\n",
    "\n",
    "<br>\n",
    "\n",
    "$s(TE) = \\left| \\left( A_{\\hbox{water}} + A_{\\hbox{fat}} \\sum_{p=1}^P \\rho_p e^{i 2 \\pi TE f_p} \\right) e^{-  TE \\, R_2^* } \\right|$\n",
    "\n",
    "**Known parameters**: TE values (say, 6 echo times), $\\rho_p$ and $f_p$ (fat NMR peak relative amplitudes and frequencies)\n",
    "\n",
    "\n",
    "**Unknown parameters**: $A_{\\hbox{water}}$, $A_{\\hbox{fat}}, R_2^*$\n",
    "\n",
    "**Fat quantification**: FF = $100 \\cdot A_{\\hbox{fat}} / (A_{\\hbox{fat}} + A_{\\hbox{water}})$\n",
    "\n",
    "</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "def CRLB_CSE_magn(a,fF,relAmps,r2,t,sigma):\n",
    "    N = np.size(t)\n",
    "    \n",
    "    FS, T = np.meshgrid(fF,t)\n",
    "    RA, T = np.meshgrid(relAmps,t)\n",
    "    \n",
    "    D = np.zeros((N,5)) + 1j*np.zeros((N,5)) \n",
    "\n",
    "    # Signal magnitude\n",
    "    shat = np.abs(a[0]*np.exp(-r2*t) + a[1]*np.sum(RA*np.exp(1j*2*np.pi*FS*T)*np.exp(-r2*T),1))\n",
    "\n",
    "    # Matrix A: partial derivatives\n",
    "    k1 = np.real(np.sum(RA*np.exp(1j*2*np.pi*FS*T),1))\n",
    "    k2 = np.imag(np.sum(RA*np.exp(1j*2*np.pi*FS*T),1))\n",
    "    J = np.zeros((4,N))\n",
    "    ds1 = (a[0]*np.exp(-2*r2*t) + a[1]*k1*np.exp(-2*r2*t))/shat\n",
    "    ds2 = (a[0]*k1*np.exp(-2*r2*t) + a[1]*(k1**2 + k2**2)*np.exp(-2*r2*t))/shat\n",
    "    ds3 = -t*shat\n",
    "    A = np.matrix([ds1, ds2, ds3]).T\n",
    "    \n",
    "    # Matrix V: diagonal with noise variances\n",
    "    V = sigma**2*np.eye(N);\n",
    "    \n",
    "    # Matrix ER: The rician \"touch\"\n",
    "    rdiag = np.zeros((N))\n",
    "    for kn in range(N):\n",
    "        trueamp = shat[kn]\n",
    "        ms = np.linspace(0,1.5*trueamp+8*sigma,1000)\n",
    "        dm = ms[1] - ms[0]\n",
    "        zs = ms*trueamp/(sigma**2)\n",
    "        \n",
    "        if trueamp/sigma < 8:\n",
    "            pric = ms/(sigma**2)*np.exp(-((ms**2 + trueamp**2)/(2*sigma**2)))*np.i0(zs);\n",
    "            rdiag[kn] = -trueamp**2/sigma**2 + 1/(sigma**2)*np.sum(dm*pric*(ms**2)*special.iv(1,zs)**2/np.i0(zs)**2);\n",
    "        else:\n",
    "            mu2 = np.sqrt(trueamp**2+sigma**2);\n",
    "            pric = (2*np.pi*sigma**2)**(-0.5) * np.exp(-0.5*(ms-mu2)**2/sigma**2);\n",
    "            rdiag[kn] = 1;\n",
    "    ER = np.diag(rdiag); \n",
    "    \n",
    "\n",
    "    \n",
    "    # Finally, invert Fisher matrix to get CRLB\n",
    "    F = np.matmul(A.H, np.matmul( np.linalg.inv(V), np.matmul(ER,A)));\n",
    "    C = np.linalg.inv(F);\n",
    "    return C        \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the NSA for fat estimation across various echo time combinations\n",
    "# for some choice of FF, R2*,  Number of Echoes, and SNR\n",
    "def showFatNSA(FF, R2star, Nechoes, SNR):\n",
    "\n",
    "    fm = 0.\n",
    "    a = [100-FF,FF]\n",
    "    fF = np.zeros((1,6))\n",
    "    fF[:] = [-3.80, -3.40, -2.60, -1.94, -0.39, 0.60];\n",
    "    fF = fF*42.576*1.5\n",
    "    relAmps = [0.087, 0.693, 0.128, 0.004, 0.039, 0.048];\n",
    "    r2 = R2star\n",
    "    phi = 0.\n",
    "    sigma = 100/SNR\n",
    "\n",
    "    nTE1 = 41;\n",
    "    ndTE = 25;\n",
    "    TE1s = np.linspace(0,2.5,nTE1)\n",
    "    dTEs = np.linspace(0.8,2.8,ndTE)\n",
    "    nTE = Nechoes\n",
    "\n",
    "\n",
    "    f_nsa_magn = np.zeros((nTE1,ndTE))\n",
    "    f_nsa_complex = np.zeros((nTE1,ndTE))\n",
    "\n",
    "    for k1 in range(np.size(TE1s)):\n",
    "        for k2 in range(np.size(dTEs)):\n",
    "            TE1 = TE1s[k1]*0.001\n",
    "            dTE = dTEs[k2]*0.001\n",
    "            TE = np.linspace(TE1,TE1 + (nTE-1)*dTE,nTE)\n",
    "\n",
    "            C = CRLB_CSE_magn(a,fF,relAmps,r2,TE,sigma);\n",
    "            f_nsa_magn[k1,k2] = sigma**2/C[1,1]\n",
    "\n",
    "            C2 = CRLB_CSE_complex(fm, a, fF, relAmps, r2, phi, TE, sigma)\n",
    "            f_nsa_complex[k1,k2] = sigma**2/C2[1,1]\n",
    "\n",
    "\n",
    "    # Display the bias and the standard deviation in R2* estimates \n",
    "    # as a function of TE combination             \n",
    "    plt.figure(figsize=(15,7))\n",
    "    plt.subplot(1,2,1)\n",
    "    plt.imshow(f_nsa_magn,interpolation='bicubic', extent = [dTEs[0], dTEs[-1], TE1s[-1], TE1s[0] ], vmin=0, vmax=nTE/3*2);\n",
    "    plt.xlabel('$\\Delta$TE (ms)');\n",
    "    plt.ylabel('TE$_1$ (ms)');\n",
    "    plt.title('Fat NSA (magnitude fitting)' );\n",
    "    plt.jet()\n",
    "    plt.colorbar()\n",
    "    plt.rc('font', **font)\n",
    "    plt.subplot(1,2,2)\n",
    "    plt.imshow(f_nsa_complex,interpolation='bicubic',extent = [dTEs[0], dTEs[-1], TE1s[-1], TE1s[0] ], vmin=0, vmax=nTE/3*2);\n",
    "    plt.xlabel('$\\Delta$TE (ms)');\n",
    "    plt.ylabel('TE$_1$ (ms)');\n",
    "    plt.title('Fat NSA (complex fitting)' );\n",
    "    plt.jet()\n",
    "    plt.colorbar()\n",
    "    plt.rc('font', **font)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cd44162a2bde46b09403f073ac14c65e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=25, description='FF', max=50), IntSlider(value=50, description='R2star',…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Create interactive slide bar to calculate NSA for magnitude and complex fitting\n",
    "y = interactive(showFatNSA,\n",
    "                FF=(0,50,1),\n",
    "                R2star=(0,100,10),\n",
    "                Nechoes=(3,12,1), \n",
    "                SNR = fixed(20))\n",
    "display(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}