example_diffsimgen.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example of how to generate simulated data using DiffSimGen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/ROBARTS/bkarat/venv/lib/python3.9/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html\n",
      "  from .autonotebook import tqdm as notebook_tqdm\n"
     ]
    }
   ],
   "source": [
    "from diffsimgen import diffsimrun\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Lets take a look at the function diffsimrun"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function diffsimrun in module diffsimgen.diffsimrun:\n",
      "\n",
      "diffsimrun(model, bval, bvec, SNR, output=None, numofsim=100000, noise_type='rician', delta=None, Delta=None, TE=None, parameter_distributions=[])\n",
      "    Function for generating training data (signal and microstructure parameters)\n",
      "    Current models accepted:\n",
      "    - NODDI_watson\n",
      "    - ball\n",
      "    - stick\n",
      "    - ball_stick\n",
      "    - zeppelin\n",
      "    bval: str or array (N,)\n",
      "      path to .bval file or numpy array with bvalues.\n",
      "    bvec: str or array (N,3)\n",
      "      path to .bvec file or numpy array with bvectors.\n",
      "    output: str\n",
      "      Output path/filename for saving as a .pkl. If none specified saves in \n",
      "      current working directory.\n",
      "    SNR: float or array\n",
      "      Signal to noise ratio. Example -> 40, [10,50], [10,20,30,40,50]\n",
      "      If you specify two SNR values ([10,100])\n",
      "      then a random SNR value will be drawn within that range with equal probability.\n",
      "      This can help if you want to cover the range of signals that might be present\n",
      "      on a voxel-to-voxel basis. If you have multi-shell data, then you can specify\n",
      "      the SNR at each shell ([10,20,30,40,50]), then the data will be simulated\n",
      "      using only those specified values. \n",
      "    numofsim: int\n",
      "      Number of simulations to perform (i.e. how many random microstructural\n",
      "      environments should data be simulated from) (default=100000).\n",
      "    noise_type: str\n",
      "      Type of noise to inject. One of 'rician' or 'gaussian'\n",
      "    delta: float (optional)\n",
      "      Pulse duration time in seconds (default=None)\n",
      "    Delta: float (optional)\n",
      "      Pulse seperation time in seconds (default=None)\n",
      "    TE: float (optional)\n",
      "      Echo time (default=None)\n",
      "    parameter_distributions: Array (optional)\n",
      "      Array of model parameters which will be used to simulate signal as opposed\n",
      "      to uniform sampling on the whole parameter range. Useful if you want to \n",
      "      simulate signal from realistic parameter distributions in the brain. *Note*\n",
      "      the index of the parameters in the array provided has to match the order of\n",
      "      the parameters seen in the <model>.json file in the resources folder. For\n",
      "      example, if you are providing a parameter distribution for the stick model,\n",
      "      the array will be [numofparameters,3] where the 0th index is theta, the 1st\n",
      "      index is phi, and the 2nd index is Dpar (as seen in stick_parameter_range.json)\n",
      "      numofsim samples will be drawn from the provided parameter distributions with\n",
      "      replacement.\n",
      "    Note, if you want to change fixed parameters (for example, in NODDI the free\n",
      "    diffusivity and parallel diffusivity are fixed) this can be adjusted\n",
      "    in the models.py file\n",
      "    \n",
      "    returns:\n",
      "      - Normalized noisy signal\n",
      "      - Normalized ground-truth (noiseless) signal\n",
      "      - Microstructure parameters\n",
      "      - Name of the microstructural parameters\n",
      "      - SNR used for each simulation\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(diffsimrun)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Can see our diffusion acquisition has 266 measurements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(266,)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bvals = np.loadtxt('datasets/tutorial_dwi.bval')\n",
    "bvecs = np.loadtxt('datasets/tutorial_dwi.bvec')\n",
    "bvals.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### In this example we simulate data using the NODDI watson model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Converting bvalues to s/m^2\n"
     ]
    }
   ],
   "source": [
    "signal, noiseless_signal, parameters, parameter_names, SNRarr = diffsimrun(model='NODDI_watson',bval=bvals,bvec=bvecs,SNR=[50,100],numofsim=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Here are the NODDI-defined microstructural parameters we simulated with"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['theta', 'phi', 'ODI', 'stick_fraction_within_watson', 'ball_fraction', 'total_watson_fraction']\n"
     ]
    }
   ],
   "source": [
    "print(parameter_names)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### We now have our simulated signal across 10000 different random microstructural environments defined by NODDI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1000, 266)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "signal.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### This contains the 1000 random microstructural parameters used to simulate signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1000, 6)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "parameters.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lets now check out a signal plot from two different simulated microstructural environments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "microstructural environment #1 (high free diffusion):\n",
      "total_watson_fraction = 0.08666370578231841\n",
      "ball_fraction = 0.9133362942176816\n",
      "stick_fraction_within_watson = 0.7745842594764226\n",
      "ODI = 0.7704539226117512\n",
      "microstructural environment #2 (high restricted diffusion):\n",
      "total_watson_fraction = 0.9439228492545744\n",
      "ball_fraction = 0.05607715074542562\n",
      "stick_fraction_within_watson = 0.29275266956584434\n",
      "ODI = 0.8623487757942124\n"
     ]
    }
   ],
   "source": [
    "print(f'microstructural environment #1 (high free diffusion):')\n",
    "for ii in range(4):\n",
    "    ii = ii+1\n",
    "    print(f'{parameter_names[-ii]} = {parameters[1,-ii]}')\n",
    "\n",
    "print(f'microstructural environment #2 (high restricted diffusion):')\n",
    "for ii in range(4):\n",
    "    ii = ii+1\n",
    "    print(f'{parameter_names[-ii]} = {parameters[2,-ii]}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Powder average data within each shell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "bunique = np.unique(bvals)\n",
    "s1 = np.zeros((len(bunique)))\n",
    "s2 = np.zeros((len(bunique)))\n",
    "\n",
    "for ii in range(len(bunique)):\n",
    "    bind = bvals == bunique[ii]\n",
    "    s1[ii] = np.mean(signal[1,bind])\n",
    "    s2[ii] = np.mean(signal[2,bind])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Fitting monoexponential function to data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import curve_fit\n",
    "# Function to calculate the exponential with constants a and b\n",
    "def monoexp(x, m, t):\n",
    "    return m * np.exp(-t * x)\n",
    "\n",
    "barray = np.linspace(0,6000,100)\n",
    "\n",
    "p0 = (1,0.0000003 ) # start with values near those we expect\n",
    "\n",
    "params1,cv = curve_fit(monoexp,bunique,s1,p0)\n",
    "params2,cv = curve_fit(monoexp,bunique,s2,p0)\n",
    "\n",
    "s1fit = monoexp(barray,params1[0],params1[1])\n",
    "s2fit = monoexp(barray,params2[0],params2[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### As expected, the more \"free\" environment has a lower signal, given by the fact that water can diffuse further in that environment. As well, the \"free\" signal decays much faster than the \"restricted\" signal, in line with the idea that the intracellular signal dominates at higher b-values "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAA9hAAAPYQGoP6dpAABd00lEQVR4nO3deXhN1/4/8PdJyCSDEJlIBEk0WiIxHKFKK5WgLtXBJV+JoXoVRdNbQRH01tRetK5qS9HB1IGiCKU1lCRqHhKZJKIqZokEiSTr98f5ZctJTpJzkjPn/XqePHHWXnvvdbbI+Vjrs9aSCSEEiIiIiMyEhaEbQERERKRNDG6IiIjIrDC4ISIiIrPC4IaIiIjMCoMbIiIiMisMboiIiMisMLghIiIis9LA0A3Qt9LSUvz9999wcHCATCYzdHOIiIhIDUII3L9/H56enrCwqL5vpt4FN3///Te8vLwM3QwiIiKqhStXrqBFixbV1ql3wY2DgwMAxcNxdHQ0cGuIiIhIHXl5efDy8pI+x6tT74KbsqEoR0dHBjdEREQmRp2UEiYUExERkVlhcENERERmhcENERERmZV6l3NDROorKSnB48ePDd0MIqonrKysapzmrQ4GN0RUiRACOTk5uHfvnqGbQkT1iIWFBVq1agUrK6s6XYfBDRFVUhbYuLq6ws7OjgteEpHOlS2ye+3aNXh7e9fp9w6DGyJSUlJSIgU2TZs2NXRziKgeadasGf7++28UFxejYcOGtb4OE4qJSElZjo2dnZ2BW0JE9U3ZcFRJSUmdrmPQ4ObQoUMYOHAgPD09IZPJ8PPPP9d4zoEDBxAcHAxra2v4+vpi3bp1Om8nUX3EoSgi0jdt/d4xaHBTUFCAwMBArFixQq36mZmZGDBgAJ5//nmcPn0aU6ZMwRtvvIE9e/bouKVERERkKgyac9OvXz/069dP7fqff/45WrVqhf/+978AgICAAPzxxx9YunQpwsLCVJ5TWFiIwsJC6XVeXl7dGk1ERERGzaRybuLj4xEaGqpUFhYWhvj4+CrPWbBgAZycnKQvne4InpgIfPut4jsRGb3evXtjypQphm6GzowcORKDBw82dDPqnREjRmD+/PnS6wcPHuCVV16Bo6MjZDKZTpdYiIuLQ8eOHVFaWqqze5gCkwpucnJy4ObmplTm5uaGvLw8PHz4UOU506dPR25urvR15coV3TQuJgbo1g2IjFR8j4nRzX2IqEojR46ETCbDuHHjKh2bMGECZDIZRo4cKZVt2bIFH3zwgR5bqF+ffPKJ2eclHjhwoFYBw4ABA/Dll18CAN58803MmzdP6fiFCxfwyiuvwMfHBzKZDMuWLVPrumfOnMGuXbswadIkqezrr7/G4cOHcfToUVy7dg1OTk4atVUT4eHhaNiwIdavX6+ze5gCkwpuasPa2lraAVxnO4EnJgKLFyuXLV7MHhwiA/Dy8sKmTZuU/sPz6NEjbNiwAd7e3kp1mzRpAgcHh1rdRwiB4uLiOrVV15ycnNC4ceMqjxcVFemvMUZECIGEhAT06NEDAHD48GHpz2UePHiA1q1bY+HChXB3d1f72suXL8drr70Ge3t7qSwjIwMBAQF45pln4O7urjJpVpt/FyNHjsSnn36qteuZIpMKbtzd3XH9+nWlsuvXr8PR0RG2trYGahWA1FTNyonqibQ04OTJJ19pabq/Z3BwMLy8vLBlyxapbMuWLfD29kZQUJBS3YrDUoWFhYiJiYGXl5c0I/Orr74C8KSHYPfu3ejUqROsra3xxx9/oLCwEJMmTYKrqytsbGzw7LPP4s8//6yxnVeuXMHrr7+Oxo0bo0mTJhg0aBCysrKk42VDSh9//DE8PDzQtGlTTJgwQZqqP2PGDMjl8krXDQwMlHohKg5L9e7dGxMnTsSUKVPg4uIi5SoePHgQXbt2hbW1NTw8PDBt2jSlwK13796YNGkSpk6diiZNmsDd3R1z5sxRuq9MJsMXX3yBl156CXZ2dggICEB8fDzS09PRu3dvNGrUCN27d0dGRobSedu2bUNwcDBsbGzQunVrzJ07V+neMpkMq1evxssvvww7Ozv4+flh+/btAICsrCw8//zzAABnZ+dKPXNVSUlJgRAC7dq1w61bt5Cenl7pWXbp0gUfffQR/vnPf8La2rrGawKK6cs//vgjBg4cqPTs/vvf/+LQoUOQyWTo3bs3AMDHxwcffPABIiMj4ejoiDfffBMA8Mcff6Bnz56wtbWFl5cXJk2ahIKCAul6hYWF+Pe//43mzZujUaNGkMvlOHDggFI7Bg4ciOPHj1d61vWKMBIAxNatW6utM3XqVPHMM88olQ0bNkyEhYWpfZ/c3FwBQOTm5tammaolJAgBVP5KSNDePYj05OHDhyIpKUk8fPiwTtdJTVX9zyI1VUsNVSEqKkoMGjRILFmyRPTp00cq79Onj1i6dKkYNGiQiIqKksp79eolJk+eLL1+/fXXhZeXl9iyZYvIyMgQ+/btE5s2bRJCCPH7778LAKJDhw5i7969Ij09Xdy+fVtMmjRJeHp6il27dokLFy6IqKgo4ezsLG7fvl1lO4uKikRAQIAYPXq0OHv2rEhKShLDhw8Xbdu2FYWFhdJ7cXR0FOPGjRPJyclix44dws7OTnz55ZdCCCHOnz8vAIj09HTpumVlaWlpSs+j/Pu1t7cX7733nrh48aK4ePGi+Ouvv4SdnZ0YP368SE5OFlu3bhUuLi4iNjZW6TxHR0cxZ84ckZqaKr7++mshk8nE3r17pToARPPmzcXmzZtFSkqKGDx4sPDx8REvvPCCiIuLE0lJSaJbt24iPDxcOufQoUPC0dFRrFu3TmRkZIi9e/cKHx8fMWfOHKXrtmjRQmzYsEGkpaWJSZMmCXt7e3H79m1RXFwsfvrpJwFApKSkiGvXrol79+5V+dwHDBggnJycRKNGjYSFhYVwcnIS9vb2QiaTCScnJ+Hk5KTyvJYtW4qlS5dWed0yJ0+eFABETk6OVHb79m0xduxYERISIq5duyb9XLRs2VI4OjqKjz/+WKSnp0tfjRo1EkuXLhWpqaniyJEjIigoSIwcOVK63htvvCG6d+8uDh06JNLT08VHH30krK2tRWqFf1hubm5i7dq1NbbZ2FT3+0eTz2+DBjf3798Xp06dEqdOnRIAxJIlS8SpU6fE5cuXhRBCTJs2TYwYMUKqf+nSJWFnZyfee+89kZycLFasWCEsLS1FXFyc2vfUSXAjhBBTpyr/Bo+J0e71ifREW8HNiROqg5sTJ7TUUBXKPsxv3LghrK2tRVZWlsjKyhI2Njbi5s2b1QY3KSkpAoD49ddfVV67LLj5+eefpbL8/HzRsGFDsX79eqmsqKhIeHp6isWLF1fZzm+//Va0bdtWlJaWSmWFhYXC1tZW7NmzR3ovLVu2FMXFxVKd1157TQwdOlR6HRgYKObNmye9nj59upDL5ZWeR/n3GxQUpNSWGTNmVGrLihUrhL29vSgpKZHOe/bZZ5XO69Kli4gp93sOgJg5c6b0Oj4+XgAQX331lVS2ceNGYWNjI73u06ePmD9/fqVn4+HhUeV18/PzBQCxe/duIcSTv5e7d++Kmly7dk1kZmaKsLAwERsbKzIzM8WIESPE22+/LTIzM0VmZqbK89QNbrZu3SosLS2VnqUQQkyePFn06tWr0jUHDx6sVDZmzBjx5ptvKpUdPnxYWFhYiIcPH4rLly8LS0tLcfXqVaU6ffr0EdOnT1cqCwoKUgoSTYW2ghuDTgU/fvy41KUIANHR0QCAqKgorFu3DteuXUN2drZ0vFWrVti5cyfeeecdfPLJJ2jRogVWr15d5TRwvVq0CBgyRDEU5e8PqOguJiL9aNasGQYMGIB169ZBCIEBAwbAxcWl2nNOnz4NS0tL9OrVq9p6nTt3lv6ckZGBx48fK+VrNGzYEF27dkVycjIAYNy4cfjuu++k4/n5+Thz5gzS09Mr5fs8evRIaSjh6aefhqWlpfTaw8MD586dk15HRERgzZo1mDVrFoQQ2Lhxo/R7tCqdOnVSep2cnIyQkBClPJAePXogPz8ff/31l5Sn1KFDB6XzPDw8cOPGDaWy8nXKJn+0b99eqezRo0fIy8uDo6Mjzpw5gyNHjuDDDz+U6pSUlODRo0d48OCBtEp2+es2atQIjo6Ole6tDnd3dzx+/BgJCQlYvnw5fHx8EB8fj7Vr18LHx0fj61X08OFDWFtbq70QXfmfJUCRjHz27FmlZGAhBEpLS5GZmYlLly6hpKQE/v7+SucVFhZW2irF1tYWDx48qOU7MX0GDW569+4NIUSVx1Vl+ffu3RunTp3SYavqQC5nUENkJEaPHo2JEycCgFoLhaqbt9eoUSON2jFv3jz8+9//VirLz89Hp06dVM5oadasmfTninvryGQypSm+w4YNQ0xMDE6ePImHDx/iypUrGDp0qFbbr25bKtYp+4BXVVZ2Xn5+PubOnYshQ4ZUup+NjY1G967J/PnzMX/+fAgh8ODBAyn/qqCgAGFhYVI+Vc+ePTW6bnkuLi548OABioqK1NrVuuLfRX5+Pv71r38pzbQq4+3tjbNnz8LS0hInTpxQCnoBKCUwA8CdO3eUfpbqG26cqU2pqUBGBuDrC/j5Gbo1RAZV1SSkWk5O0lh4eDiKioogk8nU6t1t3749SktLcfDgwUrraVWlTZs2sLKywpEjR9CyZUsAir25/vzzTylR2dXVFa6urkrnBQcHY/PmzXB1da3TDM4WLVqgV69eWL9+PR4+fIgXX3yx0r1qEhAQgJ9++glCCCn4OHLkCBwcHNCiRYtat00dwcHBSElJga+vb62voe5eROPGjcPrr7+Ozz77DH/99Rfmz5+PH374AXv37sWqVasAAM2bN691OwCgY8eOAICkpCTpz5oIDg5GUlJSlc8jKCgIJSUluHHjRrVBWFkPYMUE+vrEpGZLGa07d4DwcKBtW6B/f8WwVHg4cPeuoVtGZDB+fop4/8SJJ1+pqfqL+y0tLZGcnIykpKRK/8tVxcfHB1FRURg9ejR+/vlnZGZm4sCBA/j++++rPKdRo0Z466238N577yEuLg5JSUkYO3YsHjx4gDFjxlR5XkREBFxcXDBo0CAcPnxYutekSZPw119/afQ+IyIisGnTJvzwww+IiIjQ6FwAGD9+PK5cuYK3334bFy9exLZt2xAbG4vo6GhYWOj2I2L27Nn45ptvMHfuXFy4cAHJycnYtGkTZs6cqfY1WrZsCZlMhl9++QU3b95Efn6+ynpNmjSBr68vkpKS0K9fP/j6+iItLQ19+/aFr68vfH19lXrvioqKcPr0aZw+fRpFRUW4evUqTp8+jfT09Crb0qxZMwQHB+OPP/5Q/yGUExMTg6NHj2LixIk4ffo00tLSsG3bNqkH0t/fHxEREYiMjMSWLVuQmZmJY8eOYcGCBdi5c6d0nYSEBFhbWyMkJKRW7TAHDG60YfhwYN8+5bJ9+4BhwwzTHiIj4ecHBAc/+dJ3h6ama1utXLkSr776KsaPH4+nnnoKY8eOVZqGq8rChQvxyiuvYMSIEQgODkZ6ejr27NkDZ2fnKs+xs7PDoUOH4O3tjSFDhiAgIABjxozBo0ePNO7JefXVV3H79m08ePCgVqsRN2/eHLt27cKxY8cQGBiIcePGYcyYMRoFGLUVFhaGX375BXv37kWXLl3QrVs3LF26VOoFU0fz5s0xd+5cTJs2DW5ublIgoEpxcTGOHDmC5557DoBiCnzZnyv6+++/ERQUhKCgIFy7dg0ff/wxgoKC8MYbb1TbnjfeeKPWC+h16NABBw8eRGpqKnr27ImgoCDMnj0bnp6eUp21a9ciMjIS7777Ltq2bYvBgwfjzz//VFrDaePGjYiIiJByluojmagu6cUM5eXlwcnJCbm5udpZ0C81VdFjU91xDlGRCXn06BEyMzPRqlUrpbwHIqrZw4cP0bZtW2zevNkgPSe3bt1C27Ztcfz4cbRq1Urv96+r6n7/aPL5zZ6buqppkaRqujCJiMi82Nra4ptvvsGtW7cMcv+srCx89tlnJhnYaBMTiuuqTZvqj9chUY6IiExP2SrEhtC5c+dKU8zrI/bc1JW/PxAWBlRMWLS0VJRzSIqIiEivGNxow8aNQMWpo6GhinIiIiLSKw5LaYOzMxAXp9gVMD2d69wQEREZEIMbbfLzY1BDRERkYByWIiIiIrPC4IaIiIjMCoMbIiIiMisMbnQpMRH49lvFdyIyOr1795Y2uDRHI0eOrNWWDFQ3zz33HDZs2GDoZkgM9XOQlZUFmUyG06dPA1BsKNqiRYsatzTRBgY3uhITA3TrBkRGKr7HxBi6RURmb+TIkZDJZBg3blylYxMmTIBMJsPIkSOlsi1btuCDDz7QYwv165NPPsG6desM3QydOnDgAGQyGe7du6fReQMGDMCXX34JAHjzzTcxb948peOrVq1Cz5494ezsDGdnZ4SGhuLYsWM1Xnf79u24fv06/vnPf2rUHk35+Phg2bJlatXV5OegYkCiTe3atUO3bt2wZMkSrV+7IgY3upCYCCxerFy2eDF7cIj0wMvLC5s2bcLDhw+lskePHmHDhg1KmwsCip2iHRwcanUfIQSKi4vr1FZdc3JyQuPGjas8XlRUpL/GGBEhBBISEtCjRw8AwOHDh6U/lzlw4ACGDRuG33//HfHx8fDy8kLfvn1x9erVaq/96aefYtSoUdXuqP748eO6vwk1lJSUoLS0tMafA30aNWoUVq5cqfN/OwxudCE1VbNyInOm5+HZ4OBgeHl5YcuWLVLZli1b4O3tjaCgIKW6FYelCgsLERMTAy8vL1hbW8PX1xdfffUVgCc9BLt370anTp1gbW2NP/74A4WFhZg0aRJcXV1hY2ODZ599Fn/++WeN7bxy5Qpef/11NG7cGE2aNMGgQYOQlZUlHS8bSvj444/h4eGBpk2bYsKECdIH44wZMyCXyytdNzAwUOqFqDgc0bt3b0ycOBFTpkyBi4sLwsLCACh2x+7atSusra3h4eGBadOmKX349O7dG5MmTcLUqVPRpEkTuLu7Y86cOUr3lclk+OKLL/DSSy/Bzs4OAQEBiI+PR3p6Onr37o1GjRqhe/fuyKiwH9+2bdsQHBwMGxsbtG7dGnPnzlW6t0wmw+rVq/Hyyy/Dzs4Ofn5+2L59OwBFL8Pzzz8PAHB2dq7UM1eVlJQUCCHQrl073Lp1C+np6ZWe5fr16zF+/Hh07NgRTz31FFavXo3S0lLs37+/yuvevHkTv/32GwYOHFjp2axcuRL/+Mc/0KhRI3z44Yc1vnchBObMmQNvb29YW1vD09MTkyZNkv4+Ll++jHfeeQcymQwymQwAsG7dOjRu3Bjbt29Hu3btYG1tjezs7Eo/B6WlpVi8eDF8fX1hbW0Nb29vqU1le1IFBQVBJpMpbSWxevVqBAQEwMbGBk899RQ+++wzpfd57NgxBAUFwcbGBp07d8apU6cqPaMXX3wRd+7cwcGDB6t8jloh6pnc3FwBQOTm5uruJgkJQgCVvxISdHdPIi15+PChSEpKEg8fPqz7xaZOVf43MHVq3a9ZjaioKDFo0CCxZMkS0adPH6m8T58+YunSpWLQoEEiKipKKu/Vq5eYPHmy9Pr1118XXl5eYsuWLSIjI0Ps27dPbNq0SQghxO+//y4AiA4dOoi9e/eK9PR0cfv2bTFp0iTh6ekpdu3aJS5cuCCioqKEs7OzuH37dpXtLCoqEgEBAWL06NHi7NmzIikpSQwfPly0bdtWFBYWSu/F0dFRjBs3TiQnJ4sdO3YIOzs78eWXXwohhDh//rwAINLT06XrlpWlpaUpPY/y79fe3l6899574uLFi+LixYvir7/+EnZ2dmL8+PEiOTlZbN26Vbi4uIjY2Fil8xwdHcWcOXNEamqq+Prrr4VMJhN79+6V6gAQzZs3F5s3bxYpKSli8ODBwsfHR7zwwgsiLi5OJCUliW7duonw8HDpnEOHDglHR0exbt06kZGRIfbu3St8fHzEnDlzlK7bokULsWHDBpGWliYmTZok7O3txe3bt0VxcbH46aefBACRkpIirl27Ju7du1flcx8wYIBwcnISjRo1EhYWFsLJyUnY29sLmUwmnJychJOTU5Xn5uXlCRsbG7Fjx44q62zZskU0atRIlJSUKJUDEK6urmLNmjUiIyNDXL58ucb3/sMPPwhHR0exa9cucfnyZZGYmCj93d++fVu0aNFCzJs3T1y7dk1cu3ZNCCHE2rVrRcOGDUX37t3FkSNHxMWLF0VBQUGln4OpU6cKZ2dnsW7dOpGeni4OHz4sVq1aJYQQ4tixYwKA2Ldvn7h27Zr0c/zdd98JDw8P8dNPP4lLly6Jn376STRp0kSsW7dOCCHE/fv3RbNmzcTw4cPF+fPnxY4dO0Tr1q0FAHHq1Cml5yGXy5V+vsqr7vePJp/fDG50peIv9ZgY3d6PSEu0FtwYIMgv+yV+48YNYW1tLbKyskRWVpawsbERN2/erDa4SUlJEQDEr7/+qvLaZcHNzz//LJXl5+eLhg0bivXr10tlRUVFwtPTUyxevLjKdn777beibdu2orS0VCorLCwUtra2Ys+ePdJ7admypSguLpbqvPbaa2Lo0KHS68DAQDFv3jzp9fTp04VcLq/0PMq/36CgIKW2zJgxo1JbVqxYIezt7aUP6V69eolnn31W6bwuXbqImHK/1wCImTNnSq/j4+MFAPHVV19JZRs3bhQ2NjbS6z59+oj58+dXejYeHh5VXjc/P18AELt37xZCPPl7uXv3rqjJtWvXRGZmpggLCxOxsbEiMzNTjBgxQrz99tsiMzNTZGZmVnnuW2+9JVq3bl3tv4ulS5eK1q1bVyoHIKZMmaJUVtN7/+9//yv8/f1FUVGRynu1bNlSLF26VKls7dq1AoA4ffq0Unn5n4O8vDxhbW0tBTMVZWZmqgxI2rRpIzZs2KBU9sEHH4iQkBAhhBBffPGFaNq0qdLzWblypcprvfzyy2LkyJEq76+t4IYrFOvKokXAkCGKoSh/f0BF9zGRWatueFbH/x6aNWuGAQMGYN26dRBCYMCAAXBxcan2nNOnT8PS0hK9evWqtl75HZczMjLw+PFjpXyNhg0bomvXrkhOTgYAjBs3Dt999510PD8/H2fOnEF6enqlfJ9Hjx4pDds8/fTTsCy3Ka+HhwfOnTsnvY6IiMCaNWswa9YsCCGwceNGREdHV9v+Tp06Kb1OTk5GSEiINLQBAD169EB+fj7++usvKU+pQ4cOSud5eHjgxo0bSmXl67i5uQEA2rdvr1T26NEj5OXlwdHREWfOnMGRI0ekIRFAkSfy6NEjPHjwAHZ2dpWu26hRIzg6Ola6tzrc3d3x+PFjJCQkYPny5fDx8UF8fDzWrl0LHx+fKs9buHAhNm3ahAMHDsDGxqbKeg8fPqzyeMWdumt676+99hqWLVuG1q1bIzw8HP3798fAgQPRoEH1H9tWVlaV/q7KS05ORmFhIfr06VPtdcorKChARkYGxowZg7Fjx0rlxcXFcHJykq7boUMHpfcfEhKi8nq2trZ48OCB2vevDQY3uiSXM6ih+svfX7NyLRs9ejQmTpwIAFixYkWN9W1tbdW6bqNGjTRqx7x58/Dvf/9bqSw/Px+dOnXC+vXrK9Vv1qyZ9OeGDRsqHZPJZCgtLZVeDxs2DDExMTh58iQePnyIK1euYOjQoVptv7ptqVinLFhSVVZ2Xn5+PubOnYshQ4ZUul/5D0l17l2T+fPnY/78+RBC4MGDB1L+VUFBAcLCwqR8qp49eyqd9/HHH2PhwoXYt29ftUEDALi4uODu3bsqj1V87jW9dy8vL6SkpGDfvn349ddfMX78eHz00Uc4ePBgpedRnq2trVKgquq4pvLz8wEoZpBVzE0qH3yr686dO2jTpo3G52mCwY2WpKUB9+8/ee3gwG2mqJ6Ty4GpU5VnDsbE6C3gDw8PR1FREWQymZQ4W5327dujtLQUBw8eRGhoqFr3aNOmDaysrHDkyBG0bNkSgGImzJ9//iklKru6usLV1VXpvODgYGzevBmurq5wdHTU7I2V06JFC/Tq1Qvr16/Hw4cP8eKLL1a6V00CAgLw008/QQghfSgeOXIEDg4OaNGiRa3bpo7g4GCkpKTA19e31tewsrICoOj1qM64cePw+uuv47PPPsNff/2F+fPn44cffsDevXuxatUqAEDz5s2Vzlm8eDE+/PBD7Nmzp1LPiypBQUHIycnB3bt34ezsXG1ddd67ra0tBg4ciIEDB2LChAl46qmncO7cOQQHB8PKyqrG96yKn58fbG1tsX//frzxxhuVjqt6nm5ubvD09MSlS5cQERGh8roBAQH49ttv8ejRIykwTUhIUFn3/PnzePXVVzVuuyYY3GhBWprq/4ympjLAoXrOgMOzlpaW0tCQOv+79PHxQVRUFEaPHo1PP/0UgYGBuHz5Mm7cuIHXX39d5TmNGjXCW2+9hffeew9NmjSBt7c3Fi9ejAcPHmDMmDFV3isiIgIfffQRBg0ahHnz5qFFixa4fPkytmzZgqlTp2oUVERERCA2NhZFRUVYunSp2ueVGT9+PJYtW4a3334bEydOREpKCmJjYxEdHV3tdGZtmD17Nl566SV4e3vj1VdfhYWFBc6cOYPz58/jP//5j1rXaNmyJWQyGX755Rf0798ftra2sLe3r1SvSZMmaNKkCZKSkjB06FD4+voiLS0Nffv2VRlgLFq0CLNnz8aGDRvg4+ODnJwcAIC9vb3K6wOK4MbFxQVHjhzBSy+9VKf3vm7dOpSUlEAul8POzg7fffcdbG1tpSDax8cHhw4dwj//+U9YW1vXOOxaxsbGBjExMZg6dSqsrKzQo0cP3Lx5ExcuXMCYMWPg6uoKW1tbxMXFoUWLFrCxsYGTkxPmzp2LSZMmwcnJCeHh4SgsLMTx48dx9+5dREdHY/jw4Xj//fcxduxYTJ8+HVlZWfj4448r3T8rKwtXr15V+z8QtcWp4FpQvsdGnXKiekUuB0aMMMgQraOjo0Y9IytXrsSrr76K8ePH46mnnsLYsWNrXE114cKFeOWVVzBixAgEBwcjPT0de/bsqfZ/7nZ2djh06BC8vb0xZMgQBAQEYMyYMXj06JHGPTmvvvoqbt++jQcPHtRqFdrmzZtj165dOHbsGAIDAzFu3DiMGTMGM2fO1PhamgoLC8Mvv/yCvXv3okuXLujWrRuWLl0qfYCro3nz5pg7dy6mTZsGNzc3aShSleLiYhw5cgTPPfccAMUU+LI/V7Ry5UoUFRXh1VdfhYeHh/Sl6gO7jKWlJUaNGqVyuLGimt5748aNsWrVKvTo0QMdOnTAvn37sGPHDjRt2hSAYrgzKysLbdq0URrKVMesWbPw7rvvYvbs2QgICMDQoUOlHKYGDRrg008/xRdffAFPT08MGjQIAPDGG29g9erVWLt2Ldq3b49evXph3bp10tRxe3t77NixA+fOnUNQUBDef/99LFq0qNK9N27ciL59+2r0d1wbMiGE0OkdjExeXh6cnJyQm5tbp+7g8n7+GXj55crlW7cCXPmcTM2jR4+QmZmJVq1aVZs8SUSV5eTk4Omnn8bJkyd1/gFuaoqKiuDn54cNGzZUWjSxTHW/fzT5/GbPTR2lpakObABFeVqafttDRESG4+7ujq+++grZ2dmGborRyc7OxowZM6oMbLSJOTd1VNPQE4emiIjqF25Wqpqvr2+dksc1weDGWKSmAhkZgK8vs5CJiIjqgMNShnbnDhAeDrRtC/Tvr5hREh4OVLFWApG+1LN0PCIyAtr6vcPgxtCGDwf27VMu27cPGDbMMO2heq9sgTBdryBKRFRR2U71tVkcsDwOSxlSaiqwZ0/l8pISRXlaGoeoSO8sLS3RuHFjaWqonZ1dtSueEhFpQ2lpKW7evAk7O7sat5moCYObOqqwNYxmx8vtIaNSejqDGzIId3d3AKjV/j1ERLVlYWEBb2/vOv+HisFNHfn5KTpgVM2KqnELhpr21tBTVjlRRTKZDB4eHnB1dcXjx48N3RwiqiesrKy0sjI2gxstqHXnir8/EBamyLEpv0eIpSUQGspeGzI4S0vLOo99ExHpGxOK9S0xEfj2W8V3ANi4URHIlBcaqignIiIijbHnRp9iYpR3SJ46VbGxYFycInk4PZ3r3BAREdUR95bSl8REoFu3yuUJCQbZUJCIiMiUcG8pY5Saqlk5ERER1QqDG33x99esnIiIiGqFwY2+yOWKHJvyYmI4JEVERKRlTCjWp0WLgCFDFENR/v4MbIiIiHSAwY2+yeUMaoiIiHSIw1JERERkVhjcEBERkVlhcENERERmhcENERERmRUGN0RERGRWGNwQERGRWWFwQ0RERGaFwQ0RERGZFQY3REREZFYY3BAREZFZYXBDREREZoV7SxlaYiI30iQiItIi9twYUkwM0K0bEBmp+B4TY+gWERERmTwGN4aSmAgsXqxctnixopyIiIhqzeDBzYoVK+Dj4wMbGxvI5XIcO3as2vrLli1D27ZtYWtrCy8vL7zzzjt49OiRnlqrRampmpUTERGRWgwa3GzevBnR0dGIjY3FyZMnERgYiLCwMNy4cUNl/Q0bNmDatGmIjY1FcnIyvvrqK2zevBkzZszQc8u1wN9fs3IiIiJSi0GDmyVLlmDs2LEYNWoU2rVrh88//xx2dnZYs2aNyvpHjx5Fjx49MHz4cPj4+KBv374YNmxYjb09RkkuB6ZOVS6LiWFSMRERUR0ZLLgpKirCiRMnEBoa+qQxFhYIDQ1FfHy8ynO6d++OEydOSMHMpUuXsGvXLvTv37/K+xQWFiIvL0/py2gsWgQkJADffKP4vnChoVtERERk8gw2FfzWrVsoKSmBm5ubUrmbmxsuXryo8pzhw4fj1q1bePbZZyGEQHFxMcaNG1ftsNSCBQswd+5crbZdq+Ry9tYQERFpkcETijVx4MABzJ8/H5999hlOnjyJLVu2YOfOnfjggw+qPGf69OnIzc2Vvq5cuaLHFhMREZG+GaznxsXFBZaWlrh+/bpS+fXr1+Hu7q7ynFmzZmHEiBF44403AADt27dHQUEB3nzzTbz//vuwsKgcq1lbW8Pa2lr7b6AGaWnA/ftPXjs4AH5+em8GERFRvWOwnhsrKyt06tQJ+/fvl8pKS0uxf/9+hISEqDznwYMHlQIYS0tLAIAQQneN1VBammLSU6dOT778/RXlREREpFsG3X4hOjoaUVFR6Ny5M7p27Yply5ahoKAAo0aNAgBERkaiefPmWLBgAQBg4MCBWLJkCYKCgiCXy5Geno5Zs2Zh4MCBUpBjDMr32KhTTkRERNpj0OBm6NChuHnzJmbPno2cnBx07NgRcXFxUpJxdna2Uk/NzJkzIZPJMHPmTFy9ehXNmjXDwIED8eGHHxrqLRAREZGRkQljGs/Rg7y8PDg5OSE3NxeOjo46ucfJk4qhqIpOnACCg3VySyIiIrOmyee3Sc2WMhUODpqVExERkfYYdFjKXPn5KbaI4mwpIiIi/WNwoyMMZIiIiAyDwY0xSkxUdP34+3P1YiIiIg0x58bYxMQA3boBkZGK7zExhm4RERGRSWFwY0wSE4HFi5XLFi9WlBMREZFaGNwYk9RUzcqJiIioEgY3xsTfX7NyIiIiqoTBjTGRy4GpU5XLYmKYVExERKQBzpYyNosWAUOGcLYUERFRLTG4MUZyOYMaIiKiWuKwFBEREZkVBjdERERkVhjcEBERkVlhcENERERmhcENERERmRUGN0RERGRWGNwQERGRWWFwQ0RERGaFwQ0RERGZFa5QbCoSE7klAxERkRrYc2MKYmKAbt2AyEjF95gYQ7eIiIjIaDG4MXaJicDixcplixcryomIiKgSBjfGLjVVs3IiIqJ6jsGNsfP316yciIionmNwY+zkcmDqVOWymBgmFRMREVWBs6VMwaJFwJAhnC1FRESkBgY3pkIuZ1BDRESkBg5LERERkVlhcENERERmhcENERERmRXm3OhRWhpw//6T1w4OgJ+f4dpDRERkjhjc6ElamuqlaVJTGeAQERFpE4el9KR8j4065URERFQ7DG6IiIjIrHBYylQlJnJRPyIiIhXYc6MnDg6alVcrJgbo1g2IjFR8j4mpU9uIiIjMiUwIIQzdCH3Ky8uDk5MTcnNz4ejoqNd7a2W2VGKiIqCpKCGBPThERGS2NPn85rCUHmllVlRqatXlDG6IiIg4LGVyVM0nr66ciIionmFwY2rkcmDqVOWymBj22hAREf1/HJYyRYsWAUOGcLYUERGRCgxuTJVczqCGiIhIBQ5LERERkVlhcENERERmhcENERERmRUGN0RERGRWGNwQERGRWWFwQ0RERGaFwQ0RERGZFa5zY04SE7mwHxER1XvsuTEXMTGK3cIjIxXfY2IM3SIiIiKDkAkhhKEboU+abJmuD2lpwP37lcsdHDTYRTwxURHQVJSQwB4cIiIyC5p8fnNYyoDS0qrfzDs1Vc0AJzW16nIGN0REVM9wWMqAVPXYaHJcUlWEVF3kREREZKYY3JgDuRyYOlW5LCaGvTZERFQvcVjKXCxaBAwZwtlSRERU7xm852bFihXw8fGBjY0N5HI5jh07Vm39e/fuYcKECfDw8IC1tTX8/f2xa9cuPbXWyMnlwIgRDGyIiKheM2jPzebNmxEdHY3PP/8ccrkcy5YtQ1hYGFJSUuDq6lqpflFREV588UW4urrixx9/RPPmzXH58mU0btxY/40nIiIio2TQ4GbJkiUYO3YsRo0aBQD4/PPPsXPnTqxZswbTpk2rVH/NmjW4c+cOjh49ioYNGwIAfHx89NlkrXJwqNtxIiIiqsxgwU1RURFOnDiB6dOnS2UWFhYIDQ1FfHy8ynO2b9+OkJAQTJgwAdu2bUOzZs0wfPhwxMTEwNLSUuU5hYWFKCwslF7n5eVp943UgZ+fIkWmzuvcEBERkcRgwc2tW7dQUlICNzc3pXI3NzdcvHhR5TmXLl3Cb7/9hoiICOzatQvp6ekYP348Hj9+jNjYWJXnLFiwAHPnztV6+7WFAQwREZF2mdRsqdLSUri6uuLLL7+EpaUlOnXqhKtXr+Kjjz6qMriZPn06oqOjpdd5eXnw8vLSV5OrVXF1YvbWEBER1Z3BghsXFxdYWlri+vXrSuXXr1+Hu7u7ynM8PDzQsGFDpSGogIAA5OTkoKioCFZWVpXOsba2hrW1tXYbrwVVrU6s9qrE6uBGmkREVA8ZbCq4lZUVOnXqhP3790tlpaWl2L9/P0JCQlSe06NHD6Snp6O0tFQqS01NhYeHh8rAxphVtfqw2qsS14QbaRIRUT1l0HVuoqOjsWrVKnz99ddITk7GW2+9hYKCAmn2VGRkpFLC8VtvvYU7d+5g8uTJSE1Nxc6dOzF//nxMmDDBUG/BOCUmAosXK5ctXqwoJyIiMnMGzbkZOnQobt68idmzZyMnJwcdO3ZEXFyclGScnZ0NC4sn8ZeXlxf27NmDd955Bx06dEDz5s0xefJkxLBXQhk30iQionpMJoQQhm6EPmmyZbou6TTnJjFRMRRVUUICgxsiIjJJmnx+m9RsKXOiao0brc2WKttIs/zQFDfSJCKieoI9N+aMs6WIiMhMsOeGFORyBjVERFTvGHxXcCIiIiJtYnBDREREZoXBDREREZkVBjdERERkVhjcEBERkVnhbKn6htPDiYjIzLHnpj7hZppERFQPMLipL7iZJhER1RMMbuqL6jbTJCIiMiMMbuoLVbt0VldORERkohjc1Bdlm2mWx800iYjIDHG2VH2yaBEwZAhnSxERkVljcFPfcDNNIiIycxoPSyUlJWH8+PEICgqCh4cHPDw8EBQUhPHjxyMpKUkXbSQiIiJSm0Y9N7t378bgwYMRHByMQYMGwc3NDQBw/fp1/PrrrwgODsa2bdsQFhamk8YSERER1UQmhBDqVg4MDMSgQYMwb948lcfnzJmDLVu24OzZs1proLbl5eXByckJubm5cHR0NHRziIiISA2afH5rNCyVmpqKiIiIKo8PGzYMaWlpmlySiIiISKs0Cm58fHywc+fOKo/v3LkTLVu2rHOjSI8SE4Fvv+VKxUREZDY0yrmZN28ehg8fjgMHDiA0NFQp52b//v2Ii4vDhg0bdNLQ+iItDbh//8lrBwfAz09HN4uJUd6SYepUxXRxIiIiE6ZRzg0AHD16FJ9++ini4+ORk5MDAHB3d0dISAgmT56MkJAQnTRUW4w55yYtTfWCwampOghwEhMVm2dWlJDAqeJERGR0NPn81nidm+7du6N79+61bhxVrXyPjTrldVLdXlMMboiIyIRpZRG/gwcPoqCgACEhIXB2dtbGJUnXuNcUERGZKY0SihctWoRZs2ZJr4UQCA8Px/PPP4+XXnoJAQEBuHDhgtYbSXWQmgrs3q0Y8yqPe00REZGZ0ii42bx5M5555hnp9Y8//ohDhw7h8OHDuHXrFjp37oy5c+dqvZH1hYODZuXVunMHCA8H2rYF+vdX9MiEhwN37z6ps2iRIsfmm28U3xcurFW7iYiIjIlGCcXOzs44evQoAgICAACjRo1CSUkJvvnmGwBAQkICXnvtNVy5ckU3rdUCY04oBrQ4Wyo8HNi3DygpeVJmaQmEhgJxcXVuJxERkT7pLKG4uLgY1tbW0uv4+HhMmTJFeu3p6Ylbt25p1lpSopVZUampwJ49lctLShTlaWk6nF9ORERkWBoNS7Vp0waHDh0CAGRnZyM1NRXPPfecdPyvv/5C06ZNtdtC0lxGRvXH09P10w4iIiID0KjnZsKECZg4cSIOHz6MhIQEhISEoF27dtLx3377DUFBQVpvJGmoTZvqj/v66qcdREREBqBRz83YsWPx6aef4s6dO3juuefw008/KR3/+++/MXr0aK02kGrB3x8IC1Pk2JRnaakor2lIilsyEBGRCdMoofi3335Dr169YFnxQ9OEGHtCsdbcvQsMG6acexMWBmzcCFS3FhG3ZCAiIiOkyee3RsFN69atce/ePYSHh2Pw4MEIDw83uQCh3gQ3ZdLSFDk2vr7q9dhwSwYiIjJCmnx+azQsdenSJRw4cADt2rXDxx9/DDc3N7z44otYvnw5srOz69Ro0hE/P6BfP/VmR1W3JQMREZGJ0Ci4AYAOHTpg5syZOHbsGDIyMvDKK69g9+7daNu2LTp27IjZs2fj+PHjumgr6Rq3ZCAiIjOgcXBTnqenJ8aNG4ddu3bh1q1bmDVrFrKyshAeHo758+drq42kL9ySgYiIzIBGU8Hj4+Nx+/ZtvPTSS1LZN998g9jYWBQUFGDw4MFYtWoVGjRogDt37mi9saQHixYBQ4YohqL8/RnYEBGRydGo52bevHlKG2OeO3cOY8aMQWhoKKZNm4YdO3ZgwYIFsLS0RLNmzbTeWNITuRwYMYKBDRERmSSNgpvTp0+jT58+0utNmzZBLpdj1apViI6Oxqefforvv/9e640kIiIiUpdGwc3du3fh5uYmvT548CD69esnve7SpYtRb5pJRERE5k+j4MbNzQ2ZmZkAgKKiIpw8eRLdyq2Lcv/+fTRs2FC7LSTD44rFRERkQjQKbvr3749p06bh8OHDmD59Ouzs7NCzZ0/p+NmzZ9Gmpn2NyLTExCgW9ouMVHyPiTF0i4iIiKqlUXDzwQcfoEGDBujVqxdWrVqFVatWwcrKSjq+Zs0a9O3bV+uNJANJTFTeigFQvGYPDhERGTGNpoK7uLjg0KFDyM3Nhb29faU9pn744QfY29trtYGk2EHh/v0nrx0c1FtwuM6qW7GYM6mIiMhIaRTclHFyclJZ3qRJkzo1hipLS1O9QHBqqh4CHK5YTEREJqhOKxST7pXvsVGnXKu4YjEREZmgWvXcUD3CFYuJiMjEMLihmsnlDGqIiMhkcFjKyDk4aFZORERU37Hnxsj5+SlGhAwyW4qIiMgEMbgxAUYZyCQmMg+HiIiMEoelSHNctZiIiIwYgxvSDFctJiIiI8fghjRT3arFRERERoDBDWmGqxYTEZGRY3BDmuGqxUREZOSMIrhZsWIFfHx8YGNjA7lcjmPHjql13qZNmyCTyTB48GDdNpCULVoEJCQA33yj+L5woaFbREREJDF4cLN582ZER0cjNjYWJ0+eRGBgIMLCwnDjxo1qz8vKysK///1v9OzZU08tJSVyOTBiBHtsiIjI6Bg8uFmyZAnGjh2LUaNGoV27dvj8889hZ2eHNWvWVHlOSUkJIiIiMHfuXLRu3bra6xcWFiIvL0/pi4iIiMyXQYOboqIinDhxAqGhoVKZhYUFQkNDER8fX+V58+bNg6urK8aMGVPjPRYsWAAnJyfpy8vLSyttJxUSE4Fvv+W0cCIiMiiDBje3bt1CSUkJ3NzclMrd3NyQk5Oj8pw//vgDX331FVatWqXWPaZPn47c3Fzp68qVK3VuN6nAhf2IiMhIGHxYShP379/HiBEjsGrVKri4uKh1jrW1NRwdHZW+SMu4sB8RERkRg+4t5eLiAktLS1y/fl2p/Pr163B3d69UPyMjA1lZWRg4cKBUVlpaCgBo0KABUlJS0KZNG902miqrbmE/JhwTEZGeGTS4sbKyQqdOnbB//35pOndpaSn279+PiRMnVqr/1FNP4dy5c0plM2fOxP379/HJJ5/Um3yatDQj2yWcC/sREZERMfiu4NHR0YiKikLnzp3RtWtXLFu2DAUFBRg1ahQAIDIyEs2bN8eCBQtgY2ODZ555Run8xo0bA0ClcnOVlqY6ZkhNNWCAU7awX/mhKS7sR0REBmLw4Gbo0KG4efMmZs+ejZycHHTs2BFxcXFSknF2djYsLEwqNUinyvfYqFOuN4sWAUOGKKIsf38GNkREZDAyIYQwdCP0KS8vD05OTsjNzTXJ5OKTJ4FOnSqXnzgBBAfrvz1ERET6oMnnN7tEiIiIyKwYfFiKNOPgoFm5QSUmcpiKiIj0jsGNifHzU8QLRjVbSpWYGOUE46lTFXk5REREOsacG9K+xETFKsUVJSSwB4eIiGqFOTdkWNUt6kdERKRjDG5I+7ioHxERGRCDG9K+skX9yuOifkREpCdMKCbd4KJ+RERkIAxuSHfkcgY1RESkdwxuSL+49g0REekYc25If2JiFFPEIyMV32NiDN0iIiIyQwxuSD8SE5UX9QMUrxMTDdMeIiIyWwxuSD+49g0REekJgxvSD659Q0REesKEYjORlmbk+02VrX1TfmiKa98QEZEOMLgxA2lpqjtAUlONLMDh2jdERKQHDG7MQPkeG3XKDari2jecGk5ERFrGnBsyHE4NJyIiHWBwQ4bBqeFERKQjDG7MgIODZuVGgVPDiYhIR5hzYwb8/BQxgVHPlqqIU8OJiEhHGNyYCaMOZFTh1HAiItIRBjdkOJwaTkREOsDghgyr4tRwgNPDiYioTphQTMaF08OJiKiOGNyQ8eD0cCIi0gIGN2Q8OD2ciIi0gMENGQ9ODyciIi1gcEPGo2x6eHmcHk5ERBribCkzlZZmYov6leH0cCIiqiMGN2YoLU31SE5qqokEONw5nIiI6oDDUmaofI+NOuVGjVPDiYhIQwxuyHhxajgREdUCgxsyXpwaTkREtcDgxgw5OGhWbrQ4NZyIiGqBCcVmyM9P0blhkrOlyuPO4UREVAsyIYQwdCP0KS8vD05OTsjNzYWjo6Ohm0PqUDVbijOoiIjqFU0+v9lzQ8av4tTwsWOB1aufvJ46VbE+DhEREZhzQ6bkzh0gJEQ5sAE4g4qIiJQwuCHTMXw4cOyY6mOcQUVERP8fh6XINKSmAnv2VH3cxkZ/bSEiIqPGnhsyDRkZ1R//6y/g2285PEVEROy5qU9MdjNNAGjTpvrj0dFP/swEYyKieo1TwesJk99MEwDCw4F9+4CSkidlFhZAaWnlugkJnCJORGRGNPn85rBUPWEWm2lu3AiEhiqXPf206rpMMCYiqrc4LEWmw9kZiItTdEOlpwO+vorp4d26Va7LLRqIiOotBjdkevz8lMfSuEUDERGVw+CmnjCbzTRVWbQIGDKEWzQQEREABjf1htlsplmVils0xMQo9+ZwBhURUb3B2VJkfhITVefhcAYVEZHJ4mwpqt+qminFGVRERPUCgxsyP1XNlOIMKiKieoHBDZkfuVyRY1Ne2QyqxERu00BEZOaYUFyPmfR2DDVRNYOKScZERPUCE4rrKbPYjkETTDImIjJpTCimGpnFdgyaYJIxEVG9weCG6oeqkomLipiDQ0RkZowiuFmxYgV8fHxgY2MDuVyOY8eOVVl31apV6NmzJ5ydneHs7IzQ0NBq6xMBUJ1kLJcDb7wBREYqhqxiYgzTNiIi0iqDBzebN29GdHQ0YmNjcfLkSQQGBiIsLAw3btxQWf/AgQMYNmwYfv/9d8THx8PLywt9+/bF1atX9dxy02bW2zFUZdEiRY7NN98Aq1dX7q1ZvJg9OEREZsDgCcVyuRxdunTB//73PwBAaWkpvLy88Pbbb2PatGk1nl9SUgJnZ2f873//Q2RkZKXjhYWFKCwslF7n5eXBy8vL7BOKK86EKlN+RpRZz5aqybffKnpsKvrmG2DECP23h4iIqqVJQrFBp4IXFRXhxIkTmD59ulRmYWGB0NBQxMfHq3WNBw8e4PHjx2jSpInK4wsWLMDcuXO10l5TUdVMqDJlM6LqTSCjChf6IyIyWwYdlrp16xZKSkrg5uamVO7m5oacnBy1rhETEwNPT0+EhoaqPD59+nTk5uZKX1euXKlzu41dTTOezHZGlCa40B8Rkdky6UX8Fi5ciE2bNuHAgQOwsbFRWcfa2hrW1tZ6bhmZBC70R0Rklgwa3Li4uMDS0hLXr19XKr9+/Trc3d2rPffjjz/GwoULsW/fPnTo0EGXzaxX6l0ejlz+ZBG/xETlwAZQvB4yhAv9ERGZEIMOS1lZWaFTp07Yv3+/VFZaWor9+/cjJCSkyvMWL16MDz74AHFxcejcubM+mlovlOXqdOr05MvfX1FeL3ChPyIis2DwqeDR0dFYtWoVvv76ayQnJ+Ott95CQUEBRo0aBQCIjIxUSjhetGgRZs2ahTVr1sDHxwc5OTnIyclBfn6+od6C2ah3qxZXxIX+iIjMgsFzboYOHYqbN29i9uzZyMnJQceOHREXFyclGWdnZ8PC4kkMtnLlShQVFeHVV19Vuk5sbCzmzJmjz6YbrZrWqjHrtWzqoizJuPzQVNlCf2WYg0NEZPQMvs6NvtWXjTPVWeemopMnFUNRFZ04AQQHa7d9Ri0xUTEUVVSkHNiU4WabRER6ZzLr3JDu1CYJuF6uWqxKWZLxt9+qPp6ayuCGiMiIMbghiZ+f4nO7Xs2Wqk51C/2V9e6UTSEnIiKjweCGlNTbQEYVVTk4MTHAli1cC4eIyIgx54aqVe/WvVGlfC8NoNhBvCLm4RAR6RRzbkgrqtqjqmxvqnqj/EJ/VeXh7N7NYSoiIiPB4IaqVO/XvVGlqjyc8puzcpiKiMigDL6IH5FJUbXhZkWLF3PBPyIiA2JwQ6SpRYsUOTbffAPExqquwy0biIgMhsNSVCVjW/emNgsT6kxZHk5iovKQVBlOFyciMhgGN1QlY1r3pqrk5jIGS3LmdHEiIqPD4IaqpSpgMMT08JqSmA2a5LxoETBkSPXTxRcvVtRhDw4Rkc4xuCGNcHp4FThdnIjIaDC4IY1wergaOF2ciMigOFuKSNs4XZyIyKDYc0OkC+XzcDIyVM+oKpsuzqEqIiKtYs8NacRQ08Nrur6hpqdXSy4HRowA+vVTffzXXxWJx5GRiu8xMfptHxGRmeLGmaQxQ22maVTr3GgqJkZ5aviIEaoTj41pA86yXidfXxN4wERk7rhxJumUoaaHm/Tna8Xp4qmpqoMbYxiqunMHGD4c2LPnSVlYGLBxI+DsrP/2kFaZ9H8SiNTE4IbqjNPD1VR+unhVfv1VMUxVxhCzqoYPB/btUy7btw8YNgyIi9NvW0irjHYxTCItY84N1Rmnh9eCqhlVqoaqFi8GvvpKUa6P2VWpqYoem5IS5fKSEkV5Wpru20A6Y9SLYRJpEYMbIkMpvwFnQgLw4ouq673xhv6SjjMyqj+enq7b+xMRaQGHpYgMSZ2hqvJ0vY1DmzbVH/f11c19iYi0iD03VGdVTcPOzgZOnlR8cTRDDeos/gcoho4SE3UzVOXvr0getrRULre0VJQzIYOITAB7bqjOKu4enp0NvPyy4qs8JiuqofysqqIixZBURbpOOt64UZE8XH62VGioopzMhh9S0QYZSIcv0sF/mGReGNyQVqgTtDBZUU3lh6pSU2teH2fxYkWPi5WVdqaPOzsrZkWlpSlybLjOjdlwcACccQcbMBzheBK8xiEMw7ARDg6c6k/mgYv4kdadPAl06lS5/MQJIDhY/+0xeYmJyuvjlO+1UYWbclI1Cp4Lh92RfZCVPpkRJyws8aBHKBod4lR/Ml5cxI/InNQ26RjgvlWkLDUVjQ7vqVQsKy1RlKelsZeOzAITiknrqkowvn37SYIxk4xrSd2k4w8+4L5VVBmn+lM9wWEp0omKS7zfvg307Vu5HpOMa6lsqKqqpGNVVq/WXl4OmabUVKBt2+qP8x8kGSkOS5HBVfz9ePKk6npMMq6l6pKOBwwAdu6sfE75IIh5OfVT2VT/ffuUV6G2tFTMiGNgQ9pgBJvuMrghMnUVN+UEVAc35TEvp/7iVH/SFSPadJfBDRlUcvKTP3NX4jqomHQ8dapyb44qH3ygHASxN6d+4FR/0hUj2nSXOTekFzXtRlyGQ/5axLwcItIXPeRzMeeGjE7FVYyTk4H/+7/K9ZiDo0Xaysspv84OAx4iUkWdmXh6/J8rgxvSG3V+rrOzlV9zqEpLapuXc+2a8orIHLoiIlWMbNNdDkuRQVS1irEqHKrSkZiYmvNyVElIUHxnbw4RlRceXvVMPC3k3HBYioxeVQv9qcKhKh1RZ5NOVVQlIpfvFWKwQ1Q/GdFMPPbckMFUXOivbDfxir77DggIUPyZw1Q6VLEnR9Umnepgrg5R/aajmXiafH4zuCGjoe5QFYepdKhiQFIx4KkqEbmiioERc3WIqI4Y3FSDwY3xUne6OHcX17PyAQ+g2KuqNjjNnIjqgDk3ZJLUnS7OGVV6VpsFAlXhNHMi0hP23JDR4owqI1Y+KNmypXa5OqqGrpiYTGSyKuZRltHWf0A5LFUNBjemQ91hKoBJxwZXU65ObbB3h8hk1PT7Whv/AWVwUw0GN6ZF3RlVFbEnxwjUZvuHiti7Q2QSaupp10auJHNuyGzUNkC5cEE5KGJvjgFUt/2DukNXFessXqx8HfbuEJEKDG7IpKi7+J+q3h325hhQxe0f5HLAw6PuQ1fVbRFRMeBhAERUb3BYikxO+aGqqmZUqbJ1K+Dt/eQ1e3OMgDYSk1WpeG7XrsCxY09es8eHSKs4LEVUR+UDEk22cVDVm1M+4GGwYwDlh67kcu317lQMisoHNkDVPT7M5yEyC+y5IZNX26RjVTh0ZYR01btTEw5vEanN2GZLseeGTF7FfzCa9OZUxERkI6Rp7462Ah5VvTvqDm8xAKJ6puIirOUZ4vcoe27ILGmzN4dDVyagpnV2ygKOMtrs8alLfg+DICK1cZ2bajC4qZ80WRCwJnv3Ak2bPnnNgMdI1dSboo2FBtWlar0eoPK0dlU5P+wVIgLA4KZaDG7qr/K9OXXpyVGlYu8OYLrDW7peQt2oGCqfRx2qAqC6DIsxKCITx+CmGgxuqIw2h67UYQoBkD6SAo2aIYe3akudYTGgdr1EqsoYJJGBMLipBoMbqoo2h65qq+KQ1+3b+h0C08daFSZH0+EtYwyAalJVAFSxTJs9R+xdIg0xuKkGgxuqji6HrrRFnQAIUO4VUjdIYnBTS5oGQDExgBD6y/nRpdr0HBmid0lbdchgTC64WbFiBT766CPk5OQgMDAQy5cvR9euXaus/8MPP2DWrFnIysqCn58fFi1ahP79+6t1LwY3pImKQ1e3bwN9+xquPdqkKki6caP6FZ+3bgWefrrmwAnQvI4xDMvpVE0fphVzflQFQKbYK6Qtdeld0laduuQ36TMAM7Y2aolJBTebN29GZGQkPv/8c8jlcixbtgw//PADUlJS4OrqWqn+0aNH8dxzz2HBggV46aWXsGHDBixatAgnT57EM888U+P9GNxQXZlC746pqikvSd0gqTbn6bNOVed5/Z0Im+xUPPL2h2V3xYdCydEnZVc85ZXqNP0oBk1WPflQvveUHI0vPgmA7rypCJLK16HayR00Ak7bngST99p2ReOUJwHQnbGKIKn8s1ZVps55+qyj62s3+XIRtMGkghu5XI4uXbrgf//7HwCgtLQUXl5eePvttzFt2rRK9YcOHYqCggL88ssvUlm3bt3QsWNHfP755zXej8ENaVvF3p2KH1wMgEjXuiIR/khFKvxxDPJKryvWeRlbMA1PPoAWIAYyiBrL4iFHCJ4ETl9jBKJQT3qOqNayf0iA96t178ExmRWKi4qKcOLECUyfPl0qs7CwQGhoKOLj41WeEx8fj+joaKWysLAw/PzzzyrrFxYWorCwUHqdl5dX94YTlVPTUEpwcOWVO6v6n7u5DHmRfh2DXApiVL2uWHYMcmzFkEoBkDplFQOna/CoNgBSJ0hSVYfMR2lyKgD95iwZNLi5desWSkpK4ObmplTu5uaGixcvqjwnJydHZf2cnByV9RcsWIC5c+dqp8FEtaRuLknFIKjiEAcDINKWmgKgqsoqvp6ORTUGQEDNQVLFOtrqXdJWHfZS1d4jb/1PQzX7vaWmT5+u1NOTl5cHLy8vA7aIqGrqBEE1BUDq5IEwSCJt0rTnSJ062uxd0lad2vRS6TMAq20dXV87rL3+Z5oZNOemqKgIdnZ2+PHHHzF48GCpPCoqCvfu3cO2bdsqnePt7Y3o6GhMmTJFKouNjcXPP/+MM2fO1HhP5twQKaiaCWaoZF3mJZGp0DS/qaoyY6ujy2tra/FPk0so7tq1K5YvXw5AkVDs7e2NiRMnVplQ/ODBA+zYsUMq6969Ozp06MCEYiITVlNiNmD+s6XYRrbR3NqozSUeTCahGACio6MRFRWFzp07o2vXrli2bBkKCgowatQoAEBkZCSaN2+OBQsWAAAmT56MXr164b///S8GDBiATZs24fjx4/jyyy8N+TaIqI7Meo0bItIrgwc3Q4cOxc2bNzF79mzk5OSgY8eOiIuLk5KGs7OzYWFhIdXv3r07NmzYgJkzZ2LGjBnw8/PDzz//rNYaN0RERGT+DD4spW8cliIiIjI9mnx+W1R7lIiIiMjEMLghIiIis8LghoiIiMwKgxsiIiIyKwxuiIiIyKwwuCEiIiKzwuCGiIiIzAqDGyIiIjIrDG6IiIjIrBh8+wV9K1uQOS8vz8AtISIiInWVfW6rs7FCvQtu7v//7Uq9vLwM3BIiIiLS1P379+Hk5FRtnXq3t1RpaSn+/vtvODg4QCaTafXaeXl58PLywpUrV7hvVQ34rNTHZ6U+Piv18Vlphs9Lfbp6VkII3L9/H56enkobaqtS73puLCws0KJFC53ew9HRkT/8auKzUh+flfr4rNTHZ6UZPi/16eJZ1dRjU4YJxURERGRWGNwQERGRWWFwo0XW1taIjY2FtbW1oZti9Pis1MdnpT4+K/XxWWmGz0t9xvCs6l1CMREREZk39twQERGRWWFwQ0RERGaFwQ0RERGZFQY3REREZFYY3GjJihUr4OPjAxsbG8jlchw7dszQTdK5Q4cOYeDAgfD09IRMJsPPP/+sdFwIgdmzZ8PDwwO2trYIDQ1FWlqaUp07d+4gIiICjo6OaNy4McaMGYP8/HylOmfPnkXPnj1hY2MDLy8vLF68WNdvTesWLFiALl26wMHBAa6urhg8eDBSUlKU6jx69AgTJkxA06ZNYW9vj1deeQXXr19XqpOdnY0BAwbAzs4Orq6ueO+991BcXKxU58CBAwgODoa1tTV8fX2xbt06Xb89rVq5ciU6dOggLQAWEhKC3bt3S8f5nKq2cOFCyGQyTJkyRSrj81KYM2cOZDKZ0tdTTz0lHedzUnb16lX83//9H5o2bQpbW1u0b98ex48fl44b/e93QXW2adMmYWVlJdasWSMuXLggxo4dKxo3biyuX79u6Kbp1K5du8T7778vtmzZIgCIrVu3Kh1fuHChcHJyEj///LM4c+aM+Mc//iFatWolHj58KNUJDw8XgYGBIiEhQRw+fFj4+vqKYcOGScdzc3OFm5ubiIiIEOfPnxcbN24Utra24osvvtDX29SKsLAwsXbtWnH+/Hlx+vRp0b9/f+Ht7S3y8/OlOuPGjRNeXl5i//794vjx46Jbt26ie/fu0vHi4mLxzDPPiNDQUHHq1Cmxa9cu4eLiIqZPny7VuXTpkrCzsxPR0dEiKSlJLF++XFhaWoq4uDi9vt+62L59u9i5c6dITU0VKSkpYsaMGaJhw4bi/PnzQgg+p6ocO3ZM+Pj4iA4dOojJkydL5XxeCrGxseLpp58W165dk75u3rwpHedzeuLOnTuiZcuWYuTIkSIxMVFcunRJ7NmzR6Snp0t1jP33O4MbLejatauYMGGC9LqkpER4enqKBQsWGLBV+lUxuCktLRXu7u7io48+ksru3bsnrK2txcaNG4UQQiQlJQkA4s8//5Tq7N69W8hkMnH16lUhhBCfffaZcHZ2FoWFhVKdmJgY0bZtWx2/I926ceOGACAOHjwohFA8m4YNG4offvhBqpOcnCwAiPj4eCGEIpi0sLAQOTk5Up2VK1cKR0dH6flMnTpVPP3000r3Gjp0qAgLC9P1W9IpZ2dnsXr1aj6nKty/f1/4+fmJX3/9VfTq1UsKbvi8noiNjRWBgYEqj/E5KYuJiRHPPvtslcdN4fc7h6XqqKioCCdOnEBoaKhUZmFhgdDQUMTHxxuwZYaVmZmJnJwcpefi5OQEuVwuPZf4+Hg0btwYnTt3luqEhobCwsICiYmJUp3nnnsOVlZWUp2wsDCkpKTg7t27eno32pebmwsAaNKkCQDgxIkTePz4sdLzeuqpp+Dt7a30vNq3bw83NzepTlhYGPLy8nDhwgWpTvlrlNUx1Z/FkpISbNq0CQUFBQgJCeFzqsKECRMwYMCASu+Jz0tZWloaPD090bp1a0RERCA7OxsAn1NF27dvR+fOnfHaa6/B1dUVQUFBWLVqlXTcFH6/M7ipo1u3bqGkpETpBx4A3NzckJOTY6BWGV7Ze6/uueTk5MDV1VXpeIMGDdCkSROlOqquUf4epqa0tBRTpkxBjx498MwzzwBQvBcrKys0btxYqW7F51XTs6iqTl5eHh4+fKiLt6MT586dg729PaytrTFu3Dhs3boV7dq143NSYdOmTTh58iQWLFhQ6Rif1xNyuRzr1q1DXFwcVq5ciczMTPTs2RP379/nc6rg0qVLWLlyJfz8/LBnzx689dZbmDRpEr7++msApvH7vd7tCk5kaBMmTMD58+fxxx9/GLopRqtt27Y4ffo0cnNz8eOPPyIqKgoHDx40dLOMzpUrVzB58mT8+uuvsLGxMXRzjFq/fv2kP3fo0AFyuRwtW7bE999/D1tbWwO2zPiUlpaic+fOmD9/PgAgKCgI58+fx+eff46oqCgDt0497LmpIxcXF1haWlbKqr9+/Trc3d0N1CrDK3vv1T0Xd3d33LhxQ+l4cXEx7ty5o1RH1TXK38OUTJw4Eb/88gt+//13tGjRQip3d3dHUVER7t27p1S/4vOq6VlUVcfR0dGkfoFbWVnB19cXnTp1woIFCxAYGIhPPvmEz6mCEydO4MaNGwgODkaDBg3QoEEDHDx4EJ9++ikaNGgANzc3Pq8qNG7cGP7+/khPT+fPVQUeHh5o166dUllAQIA0jGcKv98Z3NSRlZUVOnXqhP3790tlpaWl2L9/P0JCQgzYMsNq1aoV3N3dlZ5LXl4eEhMTpecSEhKCe/fu4cSJE1Kd3377DaWlpZDL5VKdQ4cO4fHjx1KdX3/9FW3btoWzs7Oe3k3dCSEwceJEbN26Fb/99htatWqldLxTp05o2LCh0vNKSUlBdna20vM6d+6c0i+MX3/9FY6OjtIvopCQEKVrlNUx9Z/F0tJSFBYW8jlV0KdPH5w7dw6nT5+Wvjp37oyIiAjpz3xequXn5yMjIwMeHh78uaqgR48elZaqSE1NRcuWLQGYyO/3Oqckk9i0aZOwtrYW69atE0lJSeLNN98UjRs3VsqqN0f3798Xp06dEqdOnRIAxJIlS8SpU6fE5cuXhRCKqYKNGzcW27ZtE2fPnhWDBg1SOVUwKChIJCYmij/++EP4+fkpTRW8d++ecHNzEyNGjBDnz58XmzZtEnZ2diY3Ffytt94STk5O4sCBA0pTUR88eCDVGTdunPD29ha//fabOH78uAgJCREhISHS8bKpqH379hWnT58WcXFxolmzZiqnor733nsiOTlZrFixwuSmok6bNk0cPHhQZGZmirNnz4pp06YJmUwm9u7dK4Tgc6pJ+dlSQvB5lXn33XfFgQMHRGZmpjhy5IgIDQ0VLi4u4saNG0IIPqfyjh07Jho0aCA+/PBDkZaWJtavXy/s7OzEd999J9Ux9t/vDG60ZPny5cLb21tYWVmJrl27ioSEBEM3Sed+//13AaDSV1RUlBBCMV1w1qxZws3NTVhbW4s+ffqIlJQUpWvcvn1bDBs2TNjb2wtHR0cxatQocf/+faU6Z86cEc8++6ywtrYWzZs3FwsXLtTXW9QaVc8JgFi7dq1U5+HDh2L8+PHC2dlZ2NnZiZdffllcu3ZN6TpZWVmiX79+wtbWVri4uIh3331XPH78WKnO77//Ljp27CisrKxE69atle5hCkaPHi1atmwprKysRLNmzUSfPn2kwEYIPqeaVAxu+LwUhg4dKjw8PISVlZVo3ry5GDp0qNK6LXxOynbs2CGeeeYZYW1tLZ566inx5ZdfKh039t/vMiGEqFvfDxEREZHxYM4NERERmRUGN0RERGRWGNwQERGRWWFwQ0RERGaFwQ0RERGZFQY3REREZFYY3BAREZFZYXBDREREZoXBDZEZ6d27N6ZMmaL3+86ZMwcdO3bU+30rSklJgbu7O+7fv2/opmgkISEB9vb2cHV1RWpqaq2ukZSUhBYtWqCgoEDLrSMyPQxuiMhsTJ8+HW+//TYcHBzUPufy5cuwtbVFfn6+DltWtQsXLqB///4YNWoUevbsib59++Lq1atKde7cuYO3334bbdu2ha2tLby9vTFp0iTk5uZKddq1a4du3bphyZIl+n4LREaHwQ0RmYXs7Gz88ssvGDlypEbnbdu2Dc8//zzs7e1107BqZGVloW/fvhg3bhyWL1+O77//Xgpw7ty5I9X7+++/8ffff+Pjjz/G+fPnsW7dOsTFxWHMmDFK1xs1ahRWrlyJ4uJifb8VIqPC4IbIzBQXF2PixIlwcnKCi4sLZs2ahaq2kMvLy4OtrS12796tVL5161Y4ODjgwYMHAICYmBj4+/vDzs4OrVu3xqxZs/D48eMq26BqeGzw4MFKgUdhYSH+/e9/o3nz5mjUqBHkcjkOHDggHb98+TIGDhwIZ2dnNGrUCE8//TR27dpV5T2///57BAYGonnz5hpdY9u2bfjHP/4BABg5ciQGDx6M+fPnw83NDY0bN8a8efNQXFyM9957D02aNEGLFi2wdu1a6fysrCzIZDIpMLG1tUWXLl2QmpqKP//8E507d4a9vT369euHmzdvSufduHEDffv2xYQJEzB//nwAgKWlJb7++mv07NkTAwYMkJ7/M888g59++gkDBw5EmzZt8MILL+DDDz/Ejh07lAKZF198EXfu3MHBgwerfE5E9UEDQzeAiLTr66+/xpgxY3Ds2DEcP34cb775Jry9vTF27NhKdR0dHfHSSy9hw4YN6Nevn1S+fv16DB48GHZ2dgAABwcHrFu3Dp6enjh37hzGjh0LBwcHTJ06tdbtnDhxIpKSkrBp0yZ4enpi69atCA8Px7lz5+Dn54cJEyagqKgIhw4dQqNGjZCUlFRt78rhw4fRuXNnpbKarnHv3j388ccf+Pbbb6Wy3377DS1atMChQ4dw5MgRjBkzBkePHsVzzz2HxMREbN68Gf/617/w4osvokWLFtJ5sbGxWLZsGby9vTF69GgMHz4cDg4O+OSTT2BnZ4fXX38ds2fPxsqVKwGgyvwaCwsLfP755zU+v9zcXDg6OqJBgye/xq2srNCxY0ccPnwYffr0qfEaRGZLK3uLE5FR6NWrlwgICBClpaVSWUxMjAgICKjynK1btwp7e3tRUFAghBAiNzdX2NjYiN27d1d5zkcffSQ6deokvY6NjRWBgYFK7Zg8ebLSOYMGDRJRUVFCCCEuX74sLC0txdWrV5Xq9OnTR0yfPl0IIUT79u3FnDlzqn2/5QUGBop58+YpldV0jfXr14vOnTtLr6OiokTLli1FSUmJVNa2bVvRs2dP6XVxcbFo1KiR2LhxoxBCiMzMTAFArF69WqqzceNGAUDs379fKluwYIFo27at2u+nOjdv3hTe3t5ixowZlY69/PLLYuTIkVq5D5Gp4rAUkZnp1q0bZDKZ9DokJARpaWkoKSnB/PnzYW9vL31lZ2ejf//+aNiwIbZv3w4A+Omnn+Do6IjQ0FDpGps3b0aPHj3g7u4Oe3t7zJw5E9nZ2bVu47lz51BSUgJ/f3+l9hw8eBAZGRkAgEmTJuE///kPevTogdjYWJw9e7baaz58+BA2NjZKZTVdo/yQVJmnn34aFhZPfjW6ubmhffv20mtLS0s0bdoUN27cUDqvQ4cOSucAUDrPzc2t0jm1kZeXhwEDBqBdu3aYM2dOpeO2trbScBZRfcXghqgeGTduHE6fPi19eXp6wsrKCq+++io2bNgAANiwYQOGDh0qDXfEx8cjIiIC/fv3xy+//IJTp07h/fffR1FRUZX3sbCwqJTnUz5HJz8/H5aWljhx4oRSe5KTk/HJJ58AAN544w1cunQJI0aMwLlz59C5c2csX768ynu6uLjg7t27SmXVXaOoqAhxcXGVgpuGDRsqvZbJZCrLSktLqzyvLLisWFbxHE3dv38f4eHhcHBwwNatWyu1C1DMrGrWrFmd7kNk6hjcEJmZxMREpdcJCQnw8/ODpaUlmjRpAl9fX+mrLICJiIhAXFwcLly4gN9++w0RERHS+UePHkXLli3x/vvvo3PnzvDz88Ply5erbUOzZs1w7do16XVJSQnOnz8vvQ4KCkJJSQlu3Lih1B5fX1+4u7tL9by8vDBu3Dhs2bIF7777LlatWlXlPYOCgpCUlFSpvKprHDhwAM7OzggMDKz2vRiLvLw89O3bF1ZWVti+fXulXqoy58+fR1BQkJ5bR2RcGNwQmZns7GxER0cjJSUFGzduxPLlyzF58uRqz3nuuefg7u6OiIgItGrVCnK5XDrm5+eH7OxsbNq0CRkZGfj000+xdevWaq/3wgsvYOfOndi5cycuXryIt956C/fu3ZOO+/v7IyIiApGRkdiyZQsyMzNx7NgxLFiwADt37gQATJkyBXv27EFmZiZOnjyJ33//HQEBAVXeMywsDPHx8SgpKZHKqrvG9u3bK/XaGKuywKagoABfffUV8vLykJOTg5ycHKX3m5WVhatXryoNKRLVRwxuiMxMZGQkHj58iK5du2LChAmYPHky3nzzzWrPkclkGDZsGM6cOaPUawMA//jHP/DOO+9g4sSJ6NixI44ePYpZs2ZVe73Ro0cjKioKkZGR6NWrF1q3bo3nn39eqc7atWsRGRmJd999F23btsXgwYPx559/wtvbG4Cit2fChAkICAhAeHg4/P398dlnn1V5z379+qFBgwbYt2+fVFbdNUwpuDl58iQSExNx7tw5+Pr6wsPDQ/q6cuWKVG/jxo3o27cvWrZsacDWEhmeTFQcGCciMlErVqzA9u3bsWfPnmrrnTx5Ei+88AJu3rypMm/FFBUVFcHPzw8bNmxAjx49DN0cIoPiOjdEZDb+9a9/4d69e7h//361WzAUFxdj+fLlZhPYAIrhyBkzZjCwIQJ7boiIiMjMMOeGiIiIzAqDGyIiIjIrDG6IiIjIrDC4ISIiIrPC4IaIiIjMCoMbIiIiMisMboiIiMisMLghIiIis8LghoiIiMzK/wNRz50P/v4m2QAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax1 = fig.add_subplot(111)\n",
    "\n",
    "ax1.scatter(barray, s1fit, s=10, c='b', marker=\"s\", label='Micro-environment #1 (free)')\n",
    "ax1.scatter(barray,s2fit, s=10, c='r', marker=\"o\", label='Micro-environment #2 (restricted)')\n",
    "ax1.scatter(bunique, s1, s=20, c='b', marker=\"s\")\n",
    "ax1.scatter(bunique,s2, s=20, c='r', marker=\"o\")\n",
    "plt.legend(loc='upper right')\n",
    "plt.xlabel(\"b-values (s/mm^2)\")\n",
    "plt.ylabel(\"S/S0\")\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
  },
  "kernelspec": {
   "display_name": "Python 3.8.10 64-bit",
   "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.9.19"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}