{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "97c8997c-be02-49aa-bd82-3675b0689d1d",
   "metadata": {},
   "source": [
    "# Mode Identification\n",
    "\n",
    "PBjam is meant to perform (semi) automated peakbagging of solar-like oscillators. The idea being that it should be relatively easy and painless to get mode frequencies of stars with a precision that is suitable for stellar modeling. Let's see if we've succeeded...\n",
    "\n",
    "The code can be roughly divided into two main sections, the mode identification stage and the detailed peakbagging stage. This tutorial only covers the mode ID, whereas the peakbagging will be covered in another notebook. The mode ID stage fits a model consisting of a known set of radial orders and angular degrees to the spectrum, so that if a peak is found in the spectrum and the model labels it l=0 for example, then roughly $\\Delta \\nu$ further along the spectrum we should expect to see another l=0 mode. If we don't then we need to vary the parameters of the model, by for example increasing or decreasing $\\Delta \\nu$.\n",
    "\n",
    "This is all done in a Bayesian framework, so we have a prior and we sample the posterior distribution which is the product of the prior and the likelihood of the data given the model. This prior is built from thousands of previous measurements of the model parameters of other stars, so that if we provide an initial estimate of $\\Delta \\nu$ and $\\nu_{\\mathrm{max}}$ for example, then PBjam knows what range of other mode parameters are important, and crucially also any correlation there might be between them. Of course this all hinges on the solar-like oscillators having a predictable pattern of modes, where it's just the scale of the pattern that remains the unknown quantity ($\\Delta \\nu$, $\\nu_{\\mathrm{max}}$ etc.), and so naturally the models in PBjam wouldn't work on other types of oscillators.\n",
    "\n",
    "Once a rough pattern is established in the mode ID stage, the mode parameters are passed to the detailed peakbagging stage, where the model is relaxed considerably such that the mode frequencies for example are all left as free parameters. This allows PBjam to capture more detailed structure in the spectrum, which could occur due, for example, to acoustic glitches, or simply to correct for inaccuracies in the mode ID model.\n",
    "\n",
    "You can read more about the idea behind PBjam in [Nielsen et al. (2021)](https://ui.adsabs.harvard.edu/abs/2021AJ....161...62N/abstract) and the current model setup in Nielsen et al. (2025, in prep). \n",
    "\n",
    "We can start with setting the target name and loading the data. You can use a custom time series or power density spectrum, but the `IO` class does the tedious work for you, so for now we'll just use this.  \n",
    "\n",
    "**Note:** The IO.psd class normalizes the PSD consistently (Parseval). If you are using a custom PSD, be sure it is Parseval normalized, otherwise the built-in prior on mode heights won't be very helpful. The easiest way to do this is to start with the custom time series and use the IO.psd class to compute the psd, otherwise you will have to normalize the spectrum so the integral of power is equivalent to variance in the light curve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "1ad76671",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pbjam import IO\n",
    "\n",
    "tgt = 'KIC5184732'\n",
    "\n",
    "psd = IO.psd(tgt, lk_kwargs={'exptime': 60, 'mission':'Kepler', 'author':'Kepler'})\n",
    " \n",
    "psd()\n",
    "\n",
    "f = psd.freq\n",
    "\n",
    "s = psd.powerdensity"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1df1390-e9cf-4cd1-9a39-1b8c13c5de39",
   "metadata": {},
   "source": [
    "For targets with very long time series (~years), and high SNR of the oscillations, it's usually not necessary to use the spectrum at full resolution. The mode ID process just takes longer. So below we've reduced the resolution by a factor of 10, which is still enough to get a good estimate of the background and also the mode ID.\n",
    "\n",
    "**Warning:** If you're working with sub-giants or red giants, some modes can be very narrow and you may therefore need the full resolution.\n",
    "\n",
    "If you are using a shorter time series, adjust the downsampling factor accordingly. Just be sure you can still make out the modes in the spectrum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "2005dcaa-2d78-4b8d-b546-701f3fb49cd4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'Frequency [$\\\\mu$Hz]')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABVsAAAMVCAYAAABp2yRaAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAz8JJREFUeJzs3XmYHGW5N+Cnk5BAQhJBRQRU4LiB4CeCKLIIgijKpogCKhx3z/HIQcGjoEBA2WRTAZFNJIiCyCpLEghJyB6yEEjIvu97Zs3s/f0RZphJZu2pmV7mvq8rF0N3ddVT3dXVVb96631T6XQ6HQAAAAAAdEqvbBcAAAAAAFAIhK0AAAAAAAkQtgIAAAAAJEDYCgAAAACQAGErAAAAAEAChK0AAAAAAAkQtgIAAAAAJEDYCgAAAACQAGErAAAAAEAChK0AAAAAAAko2LC1uLg4brjhhvj0pz8d73znO6Nfv36x3377xQknnBBDhgyJrVu3Nvu6kpKSGDJkSBx66KGx++67x+DBg+MTn/hE3HzzzVFVVdW9KwEAAAAA5I1UOp1OZ7uIpI0aNSrOPffcWLduXURE9O3bN/r3798kYJ0xY0Z87GMfa/K6ZcuWxfHHHx9Lly6NiIj+/ftHbW1tVFZWRkTEYYcdFiNHjow99tijO1YDAAAAAMgjBdeydfz48fGlL30p1q1bF1/5ylfilVdeiYqKitiyZUuUlZXFlClT4le/+lUMHjy4yetqamritNNOi6VLl8a73/3ueOGFF6KsrCzKy8vj4YcfjoEDB8aMGTPim9/8ZpbWDAAAAADIZQXVsrW8vDwOPfTQWLx4cfzkJz+JP/7xj+1+7X333Rff+973IiJiwoQJcdRRRzV5/h//+Eecd955ERHx4osvxoknnphc4QAAAABA3iuolq0PPvhgLF68OPbee+/43e9+16HXPvDAAxERccIJJ+wUtEZEnHPOOXHAAQdERMTQoUM7XywAAAAAUFAKKmytD0HPPvvs2HXXXdv9uvLy8hg/fnxERJxyyinNTpNKpeILX/hCRESMGDGik5UCAAAAAIWmYMLWysrKmDp1akREHH744bF8+fL4wQ9+EO95z3uib9++8a53vStOO+20ePbZZ3d67Zw5c6Kuri4iIg455JAWl1H/3Nq1a2Pz5s1dsBYAAAAAQL7qk+0CkrJ06dKoqqqKiGjos7WkpCT69u0bAwYMiPXr18czzzwTzzzzTHzve9+Lu+++O1KpVERErF69umE+++67b4vLaPzc6tWrY88999xpmsrKyqisrGz4/7q6uti8eXO8/e1vb1geAAAAANA+6XQ6SkpKYp999olevXK77WjBhK1btmxp+Pu3v/1tvO1tb4tHH300zjjjjNhll11i+fLlcckll8Sjjz4a9957bxx00EHxs5/9LCIiSkpKGl7bv3//FpfR+LnGr2nsuuuui6uuuqqzqwMAAAAANLJixYrYb7/9sl1GqwombK3vBqD+7/vuuy/OPPPMhsfe+973xsMPPxzz58+PmTNnxrXXXhsXXnhh9OmT7Ftw6aWXNoS4ERFFRUXx3ve+N1asWBGDBg1KdFkAAAAAUOiKi4vjPe95TwwcODDbpbSpYMLWxm/2Bz7wgSZBa71evXrFJZdcEt/61rdi06ZNMW3atPjkJz/Z5LXl5eUtLqPxcy19uP369Yt+/frt9PigQYOErQAAAACQoXzoojO3OznogMb9qX74wx9ucbqDDz644e9ly5ZFRMQ+++zT8NiqVatafG3j5xq/BgAAAACgYMLWPffcs9XBreql0+mGv+vT8IMOOqihc91Zs2a1+Nr65/bee+9mB8cCAAAAAHqugglbIyJOPvnkiIiYM2dOi9O88cYbDX8fcMABEbF94Kujjz46IiKGDRvW7OvS6XQMHz68yXIAAAAAAOoVVNj67W9/OyIiFi5cGE8++eROz9fV1cVNN90UEdu7Hfj4xz/e8NwFF1wQERGjRo2KyZMn7/TaRx99NBYvXhwREeeff37SpQMAAAAAea6gwtZjjz02vvrVr0ZExPe+97147LHHoqamJiIili9fHueee2689tprERFxzTXXNHQdELE9bD300EMjnU7HWWedFSNHjoyI7QHto48+Gt///vcjIuKUU06JE088sTtXCwAAAADIA6l0405MC0BZWVl88YtfjJdffjkiIvr16xf9+/ePLVu2NExz5ZVXxpAhQ3Z67dKlS+OEE06IpUuXRsT27gXq6uqioqIiIiIOO+ywGDlyZOyxxx7trqe4uDgGDx4cRUVFMWjQoMxXDAAAAAB6oHzK1wqqZWtExIABA2LUqFFxzz33xHHHHRcDBgyI0tLS2HfffeOcc86J8ePHNxu0RkTsv//+8dprr8UVV1wRhxxySKRSqdhll13i8MMPj5tuuikmTZrUoaAVAAAAAOg5Cq5la67Jp+QdAAAAAHJNPuVrBdeyFQAAAAAgG4StAAAAAAAJELYCAAAAACRA2AoAAAAAkABhKwAAAABAAoStAAAAAAAJELYCAAAAACRA2AoAAAAAkABhKwAAAABAAoStAAAAAAAJELYCAAAAACRA2AoAAAAAkABhKwAAAABAAoStAAAAAAAJELYCAAAAACRA2AoAAAAAkABhKwAAAABAAoStAAAAAAAJELYCAAAAACRA2AoAAAAAkABhKwAAAABAAoStAAAAAAAJELYCAAAAACRA2AoAAAAAkABhKwAAAABAAoStAAAAAAAJELYCAAAAACRA2AoAAAAAkABhKwAAAABAAoStAAAAAAAJELYCAAAAACRA2AoAAAAAkABhKwAAAABAAoStAAAAAAAJELYCAAAAACRA2AoAAAAAkABhKwAAAABAAoStAAAAAAAJELYCAAAAACRA2AoAAAAAkABhKwAAAABAAoStAAAAAAAJELYCAAAAACRA2AoAAAAAkABhKwAAAABAAoStAAAAAAAJELYCAAAAACRA2AoAAAAAkABhKwAAAABAAoStAAAAAAAJELYCAAAAACRA2AoAAAAAkABhKwAAAABAAoStAAAAAAAJELYCAAAAACRA2AoAAAAAkABhKwAAAABAAoStAAAAAAAJELYCAAAAACRA2AoAAAAAkABhKwAAAABAAoStAAAAAAAJELYCAAAAACRA2AoAAAAAkABhKwAAAABAAoStAAAAAAAJELYCAAAAACRA2AoAAAAAkABhKwAAAABAAoStAAAAAAAJELYCAAAAACRA2AoAAAAAkABhKwAAAABAAoStAAAAAAAJELYCAAAAACRA2AoAAAAAkABhKwAAAABAAoStAAAAAAAJELYCAAAAACRA2AoAAAAAkABhKwAAAABAAoStAAAAAAAJELYCAAAAACRA2AoAAAAAkABhKwAAAABAAoStAAAAAAAJELYCAAAAACRA2AoAAAAAkABhKwAAAABAAoStAAAAAAAJELYCAAAAACRA2AoAAAAAkABhKwAAAABAAoStAAAAAAAJELYCAAAAACRA2AoAAAAAkICCClv/+te/RiqVavPfiy++2OI8Fi1aFD/84Q/jgAMOiF133TXe+c53xuc///l47LHHunFNAAAAAIB80yfbBXSFXr16xTvf+c4Wn+/Xr1+zjz/33HNx9tlnR3l5eUREDBo0KDZv3hwjRoyIESNGxLe//e247777IpVKdUndAAAAAED+KqiWrfXe8573xNq1a1v8d+yxx+70miVLlsTXvva1KC8vj6OPPjrmzZsXRUVFUVRUFFdccUVERNx///1x4403dvfqAAAAAAB5oCDD1kxcccUVUVZWFnvvvXc888wz8cEPfjAiInbfffe46qqr4gc/+EFERFxzzTWxZcuWbJYKAAAAAOQgYWtElJWVNfTJ+l//9V/xtre9badpLr300oiIKC4ujieffLIbqwMAAAAA8oGwNSLGjRsX27Zti4iIU045pdlp9t9//zjooIMiImLEiBHdVhsAAAAAkB8KMmzdsGFDHH744bH77rvHbrvtFgceeGB885vfjNGjRzc7/axZsxr+PuSQQ1qcb/1zs2fPTrReAAAAACD/FWTYWl5eHtOnT4++fftGXV1dLFmyJB566KE44YQT4jvf+U7U1NQ0mX716tUREbHHHnvEbrvt1uJ899133ybTN6eysjKKi4ub/AMAAAAACl9Bha377LNPXHnllTFz5syoqKiIzZs3R3l5eYwfPz5OOumkiIi4//7746c//WmT15WUlERERP/+/Vudf/3z9dM357rrrovBgwc3/HvPe97TmVUCAAAAAPJEQYWtJ598cgwZMiQ++tGPRr9+/SIionfv3vHpT386hg8fHmeccUZERPzpT3+KBQsWdEkNl156aRQVFTX8W7FiRZcsBwAAAADILQUVtramV69ecdNNN0VERF1dXfz73/9ueG7gwIERsb37gdbUP18/fXP69esXgwYNavIPAAAAACh8PSZsjYh4//vfH+94xzsiImLx4sUNj++zzz4REbFly5bYtm1bi69ftWpVk+kBAAAAAOr1qLC1JYccckjD37NmzWpxuvrnPvKRj3R5TQAAAABAfulRYeuiRYti48aNERFxwAEHNDx+zDHHxG677RYREcOGDWv2tcuWLYs5c+ZExPa+YQEAAAAAGiuYsDWdTrf5/M9//vOI2N5/66mnntrw3IABA+Kss86KiIg777wzioqKdnr9DTfcEBHb+2s988wzE6oaAAAAACgUBRO2Llu2LI488si46667YvHixQ3ha11dXUyaNClOOeWUeOKJJyIi4oc//GF86EMfavL6q6++OgYMGBBr1qyJ0047LRYsWBAREWVlZXH11VfHn//854iI+PWvfx177LFHN64ZAAAAAJAPUum2moTmiaVLlzbpGqBfv34xcODAKCkpicrKyobHv/3tb8fdd98dffr02Wkezz33XJx99tlRXl4eERGDBw+O0tLSqK2tbXjtfffdF6lUqt11FRcXx+DBg6OoqCgGDRqU6eoBAAAAQI+UT/lawYSt27Zti/vuuy8mTpwYr776amzYsCG2bNkSu+66a+y3337x6U9/Or7zne/E0Ucf3ep8Fi1aFDfccEO88MILsWbNmhg4cGAcdthh8cMf/rChq4GOyKeNAQAAAAByTT7lawUTtuaqfNoYAAAAACDX5FO+VjB9tgIAAAAAZJOwFQAAAAAgAcJWAAAAAIAECFsBAAAAABIgbAUAAAAASICwFQAAAAAgAcJWAAAAAIAECFsBAAAAABIgbAUAAAAASICwFQAAAAAgAcJWAAAAAIAECFsBAAAAABIgbAUAAAAASICwFQAAAAAgAcJWAAAAAIAECFsBAAAAABIgbAUAAAAASICwFQAAAAAgAcJWAAAAAIAECFsBAAAAABIgbAUAAAAASICwFQAAAAAgAcJWAAAAAIAECFsBAAAAABIgbAUAAAAASICwFQAAAAAgAcJWAAAAAIAECFsBAAAAABIgbAUAAAAASICwFQAAAAAgAcJWAAAAAIAECFsBAAAAABIgbAUAAAAASICwFQAAAAAgAcJWAAAAAIAECFsBAAAAABIgbAUAAAAASICwFQAAAAAgAcJWAAAAAIAECFsBAAAAABIgbAUAAAAASICwFQAAAAAgAcJWAAAAAIAECFsBAAAAABIgbAUAAAAASICwFQAAAAAgAcJWAAAAAIAECFsBAAAAABIgbAUAAAAASICwFQAAAAAgAcJWAAAAAIAECFsBAAAAABIgbAUAAAAASICwFQAAAAAgAcJWAAAAAIAECFsBAAAAABIgbAUAAAAASICwFQAAAAAgAcJWAAAAAIAECFsBAAAAABIgbAUAAAAASICwFQAAAAAgAcJWAAAAAIAECFsBAAAAABIgbAUAAAAASICwFQAAAAAgAcJWAAAAAIAECFsBAAAAABIgbAUAAAAASICwFQAAAAAgAcJWAAAAAIAECFsBAAAAABIgbAUAAAAASICwFQAAAAAgAcJWAAAAAIAECFsBAAAAABIgbAUAAAAASICwFQAAAAAgAcJWAAAAAIAECFsBAAAAABIgbAUAAAAASICwFQAAAAAgAcJWAAAAAIAECFsBAAAAABIgbAUAAAAASICwFQAAAAAgAcJWAAAAAIAECFsBAAAAABIgbAUAAAAASICwFQAAAAAgAcJWAAAAAIAECFsBAAAAABIgbAUAAAAASICwFQAAAAAgAcJWAAAAAIAECFsBAAAAABIgbAUAAAAASEDBh63XX399pFKphn+tKSkpiSFDhsShhx4au+++ewwePDg+8YlPxM033xxVVVXdVDEAAAAAkI9S6XQ6ne0iusq8efPiYx/7WFRUVDQ81tLqLlu2LI4//vhYunRpRET0798/amtro7KyMiIiDjvssBg5cmTsscceHaqhuLg4Bg8eHEVFRTFo0KDMVgQAAAAAeqh8ytcKtmVrXV1dfOc734mKioo46qijWp22pqYmTjvttFi6dGm8+93vjhdeeCHKysqivLw8Hn744Rg4cGDMmDEjvvnNb3ZT9QAAAABAvinYsPW2226LCRMmxDe+8Y04+eSTW532gQceiNdffz0iIh577LE46aSTIiKiV69e8fWvfz3uuuuuiIh47rnnYuTIkV1bOAAAAACQlwoybF2yZEn86le/ire//e1x6623tjn9Aw88EBERJ5xwQrOtYM8555w44IADIiJi6NChyRYLAAAAABSEggxbv//970dZWVnccsst8c53vrPVacvLy2P8+PEREXHKKac0O00qlYovfOELERExYsSIZIsFAAAAAApCwYWt99xzT4wcOTJOOumkOP/889ucfs6cOVFXVxcREYccckiL09U/t3bt2ti8eXMyxQIAAAAABaNPtgtI0qpVq+LnP/957Lbbbg39rLZl9erVDX/vu+++LU7X+LnVq1fHnnvu2ex0lZWVUVlZ2fD/xcXF7aoDAAAAAMhvBdWy9Yc//GEUFRXFkCFD4sADD2zXa0pKShr+7t+/f4vTNX6u8Wt2dN1118XgwYMb/r3nPe9pVx0AAAAAQH4rmLD1b3/7Wzz77LPxsY99LH72s59lrY5LL700ioqKGv6tWLEia7UAAAAAAN2nILoRWLduXVx00UXRu3fvuOeee6JPn/av1sCBAxv+Li8vb3G6xs81fs2O+vXrF/369Wv38gEAAACAwlAQLVt/+ctfxqZNm+IHP/hBfPjDH47S0tIm/6qqqhqm3fGxffbZp+G5VatWtbiMxs81fg0AAAAAQESBhK1LliyJiIg777wzBg4cuNO/6667rmHa+sf+7//+LyIiDjrooOjVa/vbMGvWrBaXUf/c3nvv3eLgWAAAAABAz1UQYWtn9O/fP44++uiIiBg2bFiz06TT6Rg+fHhERJx88sndVhsAAAAAkD8KImwdPXp0pNPpFv9deeWVDdPWP/b73/++4bELLrggIiJGjRoVkydP3mn+jz76aCxevDgiIs4///yuXRkAAAAAIC8VRNjaWRdccEEceuihkU6n46yzzoqRI0dGRERdXV08+uij8f3vfz8iIk455ZQ48cQTs1kqAAAAAJCj+mS7gFzQp0+fePrpp+OEE06IpUuXxkknnRT9+/ePurq6qKioiIiIww47LB566KEsVwoAAAAA5CotW9+0//77x2uvvRZXXHFFHHLIIZFKpWKXXXaJww8/PG666aaYNGlS7LHHHtkuEwAAAADIUal0Op3OdhGFrLi4OAYPHhxFRUUxaNCgbJcDAAAAAHkln/I1LVsBAAAAABIgbAUAAAAASICwFQAAAAAgAcJWAAAAAIAECFsBAAAAABIgbAUAAAAASECfTF70ne98J+k6IiJi8ODBceutt3bJvAEAAAAAulIqnU6nO/qiXr16RSqVigxe2nIhqVS8613vitWrVyc2z1xQXFwcgwcPjqKiohg0aFC2ywEAAACAvJJP+VpGLVsjInbdddf42te+llghDzzwQGLzAgAAAADobhmHrYMHD477778/sUKErQAAAABAPjNAFgAAAABAAjJq2XrhhRfG4MGDEy2kK+YJAAAAANBdMhogi/bLpw58AQAAACDX5FO+phsBAAAAAIAEZDxA1o6WL18eERHvfve7Y5dddmlz+lWrVkVtbW28973vTaoEAAAAAICsSaxl6/777x8HHHBAHH300bFmzZo2pz/iiCPiwAMPTGrxAAAAAABZlWg3Aul0OqZOnRqf+tSn4rXXXmvX9AAAAAAAhSDRsHXQoEFx0EEHxYoVK+KYY46JZ555JsnZAwAAAADkrETD1v79+8eECRPis5/9bJSWlsaXv/zluPXWW5NcBAAAAABATko0bI2IGDx4cAwfPjy++93vRm1tbVxyySXxox/9KGpra5NeFAAAAABAzkg8bI2I6N27d9xzzz1x3XXXRUTEPffcE6ecckoUFRV1xeIAAAAAALKuS8LWer/4xS/ikUceiV133TVGjhwZRx11VCxZsqQrFwkAAAAAkBVdGrZGRHz1q1+NUaNGxV577RVz586NT33qUzFu3LiuXiwAAAAAQLfq8rA1IuLII4+MSZMmxUc+8pHYsGFDfO5zn4vNmzd3x6IBAAAAALpFt4StERHve9/7YsKECfG5z30uKisro7q6ursWDQAAAADQ5fokNaMrr7wydt9991anGThwYDz33HNx4YUXxnPPPZfUogEAAAAAsi6VTqfT2S6ikBUXF8fgwYOjqKgoBg0alO1yAAAAACCv5FO+1m3dCAAAAAAAFDJhKwAAAABAAjLus3Xo0KGJFHD++ecnMh8AAAAAgGzKuM/WXr16RSqV6tzCU6moqanp1DxyXT71KQEAAAAAuSaf8rWMW7ZGRBhbCwAAAABgu4zD1iVLljT7eDqdjiOPPDI2bdoUixcvzrgwAAAAAIB8knHY+r73va/F53r37t3mNAAAAAAAhaRXtgsAAAAAACgEwlYAAAAAgAQIWwEAAAAAEiBsBQAAAABIgLAVAAAAACABwlYAAAAAgAQIWwEAAAAAEtAn0xf27t2709OkUqmoqanJtAQAAAAAgJyRcdiaTqeTrAMAAAAAIK9lHLZeeeWVSdYBAAAAAJDXUmlNVLtUcXFxDB48OIqKimLQoEHZLgcAAAAA8ko+5WsGyAIAAAAASICwFQAAAAAgARn32bqjM844Iz7+8Y/HYYcdFh//+Mdjv/32a3Ha4uLiqKqqine84x1JLR4AAAAAIKsS67O1V69ekUqlGv7/7W9/exx22GEN4ethhx0WH/jAByJie9i65557xo9+9KO4/fbbk1h8zsqnPiUAAAAAINfkU76WWNh68cUXx4wZM+LVV1+NrVu3vrWARgHs7rvvHh/96EejT58+MWbMmNhrr71i7dq1SSw+Z+XTxgAAAAAAuSaf8rXEwtbGFi9eHDNmzIjp06fH9OnTY8aMGbF+/frtC3wzfE2n08JWAAAAAKBV+ZSvJdZna2MHHnhgHHjggXHWWWc1PLZ48eJ4+OGH46abboqioqL49Kc/HZdffnlXLB4AAAAAoNv16q4FHXjggXHZZZfF9OnTY5999om6urr4zGc+012LBwAAAADoUt0Wttbbf//948Ybb4xJkybFtdde292LBwAAAADoEt0etkZEnHbaaRER8fDDD2dj8QAAAAAAictK2DpgwIDYZZddYsWKFdlYPAAAAABA4hILW4cMGRJPP/10rFq1qs1ply9fHtXV1fGud70rqcUDAAAAAGRVn6RmdPXVV0cqlYqIiHe+853x8Y9/PD7+8Y/H4YcfHh//+Mfjfe97X0REVFVVxSWXXBIRERdccEFSiwcAAAAAyKpUOp1OJzGjz3zmMzFz5swoLi5+a+Zvhq8REYMHD453vetdsXLlyigvL4/TTjst/vnPf0bfvn2TWHzOKi4ujsGDB0dRUVEMGjQo2+UAAAAAQF7Jp3wtsbC13qJFi2LGjBkxffr0mDFjRsyYMSPWr1+/84JTqdh1113j0EMPjcMOOyw+9rGPxWGHHRYf/ehHY9ddd02ypKzKp40BAAAAAHJNPuVriYetzVm9enVD+Fr/3+XLl79VRKMWsL17946qqqquLqnb5NPGAAAAAAC5Jp/ytYzD1pdeeimOP/746NUrszG2Nm/e3CR8nT59eixcuDDS6XTU1tZmNM9clE8bAwAAAADkmnzK1zIOW3v16hV77rlnfPGLX4zTTz89TjnllBgwYECniikrK4uZM2fGpz/96U7NJ5fk08YAAAAAALkmn/K1jMPWj33sY/Haa69tn0kqFX379o0TTzwxTj/99Dj99NNj7733TrTQfJVPGwMAAAAA5Jp8ytc61Wfr8uXL46mnnoonn3wyxo4dGzU1NZFKpSKVSsURRxwRZ555Zpxxxhlx0EEHJVlzXsmnjQEAAAAAck0+5WuJDZC1devWePbZZ+PJJ5+M4cOHR2lpacPAV//xH/8RZ555Zpx++ulx9NFHNxkQq9Dl08YAAAAAALkmn/K1xMLWxqqqquKll16KJ554Ip555plYs2bN9oWlUvGOd7wjTjvttDj99NPj5JNPjl133TXpxeeUfNoYAAAAACDX5FO+1iVh646mTJkSTz75ZDz11FMxZ86c7QtOpWK33XaLz33uc3HGGWfEGWecEXvssUdXl9Lt8mljAAAAAIBck0/5WreErY0tXLiwoZ/XiRMnRl1dXaRSqbjyyivjiiuu6M5SukU+bQwAAAAAkGvyKV/r1d0LfP/73x8XX3xxjB07NtauXRv33ntvnHbaadG/f//uLgUAAAAAIDHd3rK1p8mn5B0AAAAAck0+5Wvd3rIVAAAAAKAQdShs3bp1a0yZMiWWLFnS4jRLliyJoUOHdrowAAAAAIB80u6wdciQIfGud70rjjrqqHj/+98fRx99dMyePXun6SZMmBDf/va3Ey0SAAAAACDXtStsfeqpp+Lqq6+Oww8/PG644Ya48MIL44033ogjjzwynnrqqa6uEQAAAAAg5/Vpz0S33HJLfPKTn4zx48dHKpWKiIhLLrkkvvzlL8fZZ58df/nLX+Kb3/xmlxYKAAAAAJDL2tWydc6cOXHuuec2BK0REfvuu2+MGTMmPvvZz8Z//ud/xj333NNlRQIAAAAA5Lp2tWytrKyM3XbbbafHd9ttt/j3v/8dZ599dvzoRz+KioqK2HPPPRMvEgAAAAAg17WrZesBBxwQU6dObfa5XXbZJf71r3/FV77ylbjooovizjvvTLRAAAAAAIB80K6w9bOf/Ww8/vjjUVVV1ezzffr0iUceeSTOPffcmDBhQqIFAgAAAADkg3Z1I/Ctb30r1qxZE9OmTYujjjqq2Wl69eoVDz74YLzrXe+K6dOnJ1okAAAAAECuS6XT6XS2iyhkxcXFMXjw4CgqKopBgwZluxwAAAAAyCv5lK+1qxsBAAAAAABaJ2wFAABI2F/GLYkhT88ONxICQM/Srj5bO6q2tjYWLFgQW7Zsierq6lanPe6447qiBAAAgKy5+pk3IiLitP+3Txz+vj2yXA0A0F0SDVtXrlwZl112WTz++OOxbdu2NqdPpVJRU1OTZAkAAAA5o7zK+Q4A9CSJha2LFy+Oo48+OtavX9/uW2XcUgMAAAAAFIrE+my97LLLYt26dfGOd7wj7rvvvli5cmVUV1dHXV1dq/8AAAAAAApBYi1bX3zxxUilUvHwww/HCSeckNRsAQAAAADyQmItWysqKmK33XYTtAIAAAAAPVJiYesBBxygD1YAAIBGUpHKdgkAQDdKLGz9+te/HhUVFTFy5MikZgkAAAAAkDcSC1svvvji+H//7//FD37wg1iyZElSswUAAMhb6XD3HwD0JIkNkLXbbrvFiy++GN///vfj0EMPja9+9avxiU98IgYOHNjq684///ykSgAAAAAAyJrEwtaIiKVLl8a6deuivLw8HnzwwXjwwQdbnT6VSglbAQAAAICCkFjY+tprr8Xxxx8fZWVlERHRt2/feMc73hF9+iSa5wIAAOQNA2QBQM+SWBJ65ZVXRmlpaRx44IFxzz33xGc+85no1SuxLmEBAAAAAHJaYmnohAkTIpVKxSOPPBInnHBCVoLW6dOnx1VXXRWnn356fPjDH463v/3tscsuu8Tb3/72OProo+Oaa66JzZs3tzqPdevWxcUXXxwf+tCHYrfddos999wzjj322Lj33nsjnda5PQAAAADQvFQ6oQRx4MCBkUqlori4OInZZeR//ud/4o477mj4/1133TV22WWXKCkpaXjsHe94Rzz99NNx1FFH7fT6adOmxec///nYtGlTRETsvvvuUVFRETU1NRER8fnPfz6efvrp6Nu3b7trKi4ujsGDB0dRUVEMGjQo01UDAADyyP6/fDYiIv723U/GMR94R5arAYD8lk/5WmLNT9///vdHdXV11NbWJjXLDjvyyCPjxhtvjIkTJ8aWLVti27ZtUVxcHCUlJfHAAw/EO9/5zti4cWOceeaZUVRU1OS1RUVFceqpp8amTZviwx/+cLzyyitRUlISZWVlcfvtt8cuu+wSw4cPj4suuig7KwcAAOSdlC5bAaBHSSxsPf/886OysjKefvrppGaZUQ2XXHJJfOpTn4q3ve1tDY/vvvvucf7558ff/va3iIhYv359PPPMM01ee9NNN8XatWtjt912i+eeey6OOOKIiNg+0NePf/zjuOqqqyIi4u6774758+d3zwoBAAB5TU9kANCzJBa2XnjhhfHZz342fvjDH8bEiROTmm2iPvWpTzX8vXLlyibPDR06NCIizjnnnDjggAN2eu1PfvKT2H333aO2tjYeeuihri0UAAAAAMg7fZKa0TXXXBNHHXVUTJ8+PY455pg45phj4sgjj4yBAwe2+rorrrgiqRLaNHbs2Ia//+M//qPh73nz5sXy5csjIuKUU05p9rW77757HHvssfH888/HiBEjGlq6AgAAAABEJBi2DhkyJFJvdkiUTqdj7NixMW7cuDZf19Vha2VlZaxZsyaeeeaZhmW9//3vj9NOO61hmlmzZjX8fcghh7Q4r0MOOSSef/75eOONN1pdXmVlZcP/Z3PAMAAAAACg+yQWth533HENYWsu2HXXXZuEnvWOPvro+Pvf/x79+vVreGz16tUNf++7774tzrP+ueLi4igtLY3dd999p2muu+46rV4BAICIMEAWAPQ0iYWto0ePTmpWidh7772joqIiSktLo6ysLCIiTjjhhPjd734X733ve5tMW1JS0vB3//79W5xn4+dKSkqaDVsvvfTS+NnPftbw/8XFxfGe97wn4/UAAAAAAPJDYgNk5ZqlS5fG2rVro7S0NNatWxc33XRTvPrqq3HkkUd2adcF/fr1i0GDBjX5BwAAAAAUvoINWxvba6+94uKLL45hw4ZFKpWK3/zmN/HMM880PN94EK/y8vIW59P4ubYG/gIAAAAAepbEuhHY0ezZs2Pq1Kmxfv36iNgeeB5xxBHxkY98pKsW2aYjjzwyjjnmmHj55Zfj7rvvjlNPPTUiIvbZZ5+GaVatWtVia9RVq1ZFRMSgQYOa7UIAAAAAAOi5Eg9bn3nmmbjsssti9uzZzT7/kY98JH7729/G6aefnvSi26V+kKuFCxc2PHbIIYc0/D1r1qw46KCDmn3trFmzIiLi4IMP7sIKAQAAAIB8lGg3AldffXWcccYZMWvWrEin09G7d+/Ya6+9Yq+99orevXtHOp2OWbNmxZe//OUYMmRIkotut8WLF0dE024APvjBDzYMmjVs2LBmX1dWVhZjx46NiIiTTz65i6sEAAAAAPJNYmHrsGHDYsiQIZFOp+O4446LESNGRElJSaxZsybWrFkTpaWlMWLEiDj++OMjnU7Hb37zmxg+fHhSi4/a2tpIp9OtTjNy5MiYMmVKREQcf/zxDY+nUqk4//zzIyLi4YcfjqVLl+702jvuuCNKS0ujd+/e8Y1vfCOxugEAAACAwpBY2HrLLbdERMTZZ58do0aNipNOOin69evX8Hzfvn3jpJNOipEjR8bZZ58d6XS64TVJWLFiRRx22GFx1113xeLFi5sErytWrIjrr78+zjjjjEin07HnnnvGT3/60yavv+SSS2LvvfeO8vLy+NKXvhTTpk2LiIiqqqq488474/LLL4+IiB/84AfxwQ9+MLG6AQAAAIDCkEq31Ry0nfbcc88oKiqK5cuXN/SL2pKVK1fGe9/73njb294WmzdvTmLxsXTp0jjggAMa/r9v374xaNCg2LZtW5SVlTU8fsABB8Rjjz0Whx122E7zmDZtWnz+85+PTZs2RcT2rgYqKiqiuro6IrZ3H/D00083CZHbUlxcHIMHD46ioqIWB94CAAAKy/6/fDYiIv73xA/ETz+nsQYAdEY+5WuJtWytqqqKt73tbW0GrRER++23X+yxxx4NIWYS9tlnn3j00Ufjxz/+cRxxxBHxjne8I4qLi6Ouri7e+973xmmnnRb33ntvzJ49u9mgNSLi8MMPj9mzZ8dPf/rT+MAHPhDV1dUxYMCAOOaYY+Kee+6J559/vkNBKwAA0LP9YeSCbJcAAHSjxFq2fvSjH4158+ZFSUlJ9O3bt9VpKysrY9CgQfHhD384Zs6cmcTic1Y+Je8AAEAy6lu2RkQsvf5LWawEAPJfPuVribVsPe+886K6ujqGDh3a5rQPPvhgVFdXx3nnnZfU4gEAAAAAsiqxsPXiiy+OY445Ji688MJ44IEHWpxu6NChceGFF8axxx4bF198cVKLBwAAAADIqj5Jzei6666L4447Ll5//fX4zne+E1deeWUcf/zxDX24rlq1KsaMGRPLly+PwYMHx7HHHhvXXntts/O64oorkioLAAAAAKBbJNZna69evSKVSkVERP0s6/+/XkuP76i2tjaJknJCPvUpAQCdUV5VE9OXbY1PHbhn9Omd2M0zAHlJn60AkJx8ytcSa9l63HHHtRmiAgCF63sPTI0JizbFTz77/rj45A9luxwAAIBul1jYOnr06KRmBQDkoQmLNkVExD+mLBe2AgAAPZJ7/AAAAAAAEiBsBQAAAABIQGLdCOxo9uzZMXXq1Fi/fn1EROy1115xxBFHxEc+8pGuWiQAkAOSGXoTAAAg/yQetj7zzDNx2WWXxezZs5t9/iMf+Uj89re/jdNPPz3pRQMAAAAAZE2i3QhcffXVccYZZ8SsWbMinU5H7969Y6+99oq99torevfuHel0OmbNmhVf/vKXY8iQIUkuGgAAAAAgqxILW4cNGxZDhgyJdDodxx13XIwYMSJKSkpizZo1sWbNmigtLY0RI0bE8ccfH+l0On7zm9/E8OHDk1o8AAAAAEBWJRa23nLLLRERcfbZZ8eoUaPipJNOin79+jU837dv3zjppJNi5MiRcfbZZ0c6nW54DQAAAABAvkssbJ06dWqkUqm45ZZbIpVKtThdKpWKm2++OSIiXnnllaQWDwAAAACQVYmFrVVVVfG2t70t9t133zan3W+//WKPPfaI6urqpBYPAAAAAJBViYWtBx54YJSWlkZVVVWb01ZWVkZpaWkceOCBSS0eAAAAACCrEgtbzzvvvKiuro6hQ4e2Oe2DDz4Y1dXVcd555yW1eAAgR6SzXQAAAECWJBa2XnzxxXHMMcfEhRdeGA888ECL0w0dOjQuvPDCOPbYY+Piiy9OavEAAAAAAFnVJ6kZXXfddXHcccfF66+/Ht/5znfiyiuvjOOPP76hD9dVq1bFmDFjYvny5TF48OA49thj49prr212XldccUVSZQEAAAAAdItUOp1O5G6/Xr16RSqVioiI+lnW/3+9lh7fUW1tbRIl5YTi4uIYPHhwFBUVxaBBg7JdDgB0mf1/+WxEROw5oG9Mv/xzWa4GILvq94kREUuv/1IWKwGA/JdP+VpiLVuPO+64NkNUAKDwORoAAAB6qsTC1tGjRyc1KwAgjxkgCwAA6KkSGyALAAAAAKAnE7YCAAAAACQgo7D16aefjhEjRiRaSFfMEwAAAACgu2TUZ+uZZ54Z7373u2PVqlWJFdIV8wQAul86rddWAACgZ8q4G4GuOJFycgYAAAAA5KuMWrZGRGzbti2GDh2aZC0AQAFIpVLZLgEAACArMg5bi4uL49vf/naStQAAAAAA5K2Mwtb3vve9XdJqZa+99kp8ngBA99ItEAAA0FNlFLYuXbo04TIAAAAAAPJbxgNkAQAAAADwFmErAAAAAEAChK0AQKL02Ar0dPquBoCeS9gKAAAAAJAAYSsAkKhUtgsAAADIEmErAAAAAEAChK0AQKL0VAgAAPRUwlYAAIAEGR8LAHouYSsAAAAAQAKErQAAAAAACRC2AgAAAAAkoE9SM9q6dWs8+eSTMWbMmFi0aFFs3rw5IiLe/va3x3/8x3/E8ccfH2eeeWYMGjQoqUUCAAAAAOSMRMLWG264Ia6//vooLi5ueCz9Zq/wqVQqxo0bFw888EBcdNFFcdlll8Ull1ySxGIBAAByjvGxAKDn6nTY+q1vfSv+/ve/N4SrvXv3jgMPPDD23HPPiIjYvHlzLF68OGpra2Pr1q3xi1/8ImbPnh33339/ZxcNAAAAAJAzOtVn61133RUPPfRQpNPpOOyww+LRRx+NrVu3xrx582LixIkxceLEmDdvXmzdujX++c9/xmGHHRbpdDqGDh0a9957b1LrAADkkLQmXUAnbKuqjeWbyrNdBgBARjIOW6urq+PXv/51pFKpOPfcc2PSpElx1llnxYABA3aadsCAAfHVr341Jk2aFOecc06k0+n41a9+FTU1NZ0qHgAAKCwn3TImjrtxVMxaVZTtUgAAOizjsPXpp5+OTZs2xQEHHBD33Xdf7LLLLm2+Zpdddom//OUvccABB8TGjRvj3//+d6aLBwAACtCqrdsiImL47LVZrgQAoOMyDltHjRoVqVQq/ud//id23XXXdr9u1113jR//+MeRTqdj5MiRmS4eAAAgJ6X1pwIAPVbGYeuMGTMiIuJzn/tch1/7+c9/vsk8AAAAAADyXcZh6/LlyyOVSsXBBx/c4dcefPDB0atXr1i+fHmmiwcAAAAAyCkZh63FxcUxcODASKVSHX5tKpWKQYMGRXFxcaaLBwAAAADIKRmHraWlpbHbbrtlvOB+/fpFWVlZxq8HAAAAAMglGYetSXT6ruN4ACg8ft+Bns5eEAB6rozDVgAAAAAA3tKnMy9et25d9O7dO6PXptPpjPp7BQAAAADIRZ0KW90mCAAAAACwXcZh65VXXplkHQAAAAAAeU3YCgAAkCA3AAJAz2WALAAAAACABHSqz9b2qKqqimHDhsW8efOiX79+8fGPfzyOOeaYrl4sAAAAAEC3yjhsLSkpiSeeeCIiIr7+9a9Hv379dppm6tSpcdZZZ8XKlSubPP7JT34yHn/88dh7770zXTwAkKPcPQskwa34AEA+yrgbgZEjR8Z//ud/xu9///tmg9b169fHF7/4xVi5cmWk0+km/yZPnhynn356pwoHAAAAAMglGYetY8eOjYiI8847r9nnb7jhhti4cWNERFxwwQUxfvz4mDlzZvz0pz+NdDod06ZNi3/961+ZLh4AAChgqVS2K8hcWht/AOixMu5GYMqUKZFKpeILX/hCs88/9NBDkUql4rTTTov777+/4fGbb745Nm/eHA888EA89thj8dWvfjXTEgAAAAAAckbGLVvXrFkTffr0iYMPPnin52bPnh3r16+PiIgLL7xwp+f/93//NyIiZsyYkeniAYAclceN0QAAADol47B13bp1MWjQoOjVa+dZTJkyJSIi+vbtG8ccc8xOzx9yyCGRSqVi9erVmS4eAMhRbp4FAAB6qozD1tra2iguLm72uWnTpkVExEEHHRR9+/bd6fk+ffrEHnvsEdu2bct08QAAAAAAOSXjsHWvvfaKmpqaWLRo0U7PTZw4MVKpVHziE59o8fWlpaUxYMCATBcPAACQk9Ka+ANAj5Vx2Prxj388IiLuvvvuJo8vWLAgXn311YiI+MxnPtPsa5ctWxZVVVWx3377Zbp4AACggAksAYB8lHHYeu6550Y6nY5bb701brzxxpg3b16MHDkyzj777Ein0zFgwIA47bTTmn3tyy+/HBHb+24FAAqMgAQAAOihMg5bzz777DjuuOOipqYmfvnLX8bBBx8cJ598crz++uuRSqXiZz/7WQwcOLDZ1z7yyCORSqWaHTwLAAAglcp2BQAAHZdx2BoR8dRTT8Wpp54a6XS64V9ExPe+97244oormn3NggULYtiwYRER8cUvfrEziwcAcpGABAAA6KH6dObFgwcPjqeffjoWLlzY0E/rJz7xiXjf+97X4mt22WWXeOqpp2KXXXaJAw88sDOLBwAAAADIGZ0KW+u9//3vj/e///3tmnb//feP/fffP4nFAgC5SJ+tQBuWbyqPR6Yuj28ffUC8Y/d+2S4HACAxiYStdXV1MXfu3CguLo4999wzPvjBDyYxWwAAoACd9ecJsaGkMl5dsTUe+t6nsl0OAEBiOtVna3V1dfziF7+IPffcMw499NA4+uij46CDDop3vvOdcc011zT04QoAAFBvQ0llRES8smRLi9M4lQAA8lGnWraeeeaZMWzYsJ1C1U2bNsUVV1wRCxYsiL/+9a+dWQQAAAAAQF7IOGx99NFH4/nnn4+I7X22nn322bHffvvF0qVL46GHHorVq1fHgw8+GN/+9rfjM5/5TGIFAwAA5DKtcgGg58o4bP3b3/4WEREnn3xyPPXUU9Gv31sd2//qV7+Kz372szFjxox46KGHhK0AAECHpFLZrgAAoOMy7rN1+vTpkUql4tZbb20StEZEDBo0KG644YZIp9MxY8aMThcJAAAAAJDrMg5bN27cGLvuumscdNBBzT5/xBFHNEwHAAAAAFDoMg5bKysrY/DgwS0+X/9cZWVlposAAPKQrgoBAICeKuOwFQAAoKvk8yBTaZedAKDHErYCAAAAACSgT2devG7duujdu3eLz6dSqVanSaVSUVNT05kSAAAAAAByQqfC1nQ+39sDAADkrFQq2xUAAHRcxmHrlVdemWQdAAAAAAB5TdgKAABkRaEOJOUGQADouQyQBQAAAACQAGErAJAofboDAAA9lbAVAADIilS0PAqW6zYAQD4StgIAAAAAJEDYCgAAkCCNcgGg5xK2AgAAOSfVcg8DAAA5S9gKAAB5YuaKrXHziHlRUV2b7VISkdYGFAAoMAUVtm7atCnuv//++OY3vxkHH3xwDBgwIPr16xf77bdfnHnmmfHEE0+0OY+SkpIYMmRIHHroobH77rvH4MGD4xOf+ETcfPPNUVVV1Q1rAQAAzTvjjvFx20sL487Ri7JdCgAAzeiT7QKStPfee0dNTU3D/++6666xyy67xKpVq2LVqlXx1FNPxSmnnBL/+te/on///ju9ftmyZXH88cfH0qVLIyKif//+UVlZGVOnTo2pU6fGQw89FCNHjow99tiju1YJAAB2smB9SbZLAACgGQXVsrWmpiaOPPLI+NOf/hSLFi2Kbdu2RWlpaSxZsiS++93vRkTE888/Hz/84Q+bfe1pp50WS5cujXe/+93xwgsvRFlZWZSXl8fDDz8cAwcOjBkzZsQ3v/nN7l4tAMgrbgqGrpcukC9aKlrumDWf1zGdz8UDAJ1SUGHrSy+9FJMnT47/+q//igMPPLDh8f333z/uvffehpD1b3/7W6xYsaLJax944IF4/fXXIyLisccei5NOOikiInr16hVf//rX46677oqIiOeeey5GjhzZHasDAAAAAOSRggpbTzjhhFafr2/dGhExderUJs898MADDfM46qijdnrtOeecEwcccEBERAwdOrSzpQIAQMY0nAQAyE0FFba2Zdddd234u7b2rRFcy8vLY/z48RERccoppzT72lQqFV/4whciImLEiBFdWCUAAAAAkI96VNg6evTohr8PPfTQhr/nzJkTdXV1ERFxyCGHtPj6+ufWrl0bmzdv7poiAQCgh0i30stzquXuXAEAclafbBfQXbZu3RrXXXddREQce+yx8aEPfajhudWrVzf8ve+++7Y4j8bPrV69Ovbcc8+dpqmsrIzKysqG/y8uLu5U3QAAQH7RywMA9Fw9omVrXV1dfOtb34o1a9bErrvuGrfffnuT50tKShr+7t+/f4vzafxc49c0dt1118XgwYMb/r3nPe/pZPUAAAAAQD7oEWHr//7v/8YzzzwTERF33HFHfPSjH+2yZV166aVRVFTU8G/FihVdtiwAAHqm1m6/zyep0FcAAFBYCr4bgUsuuaShJeutt94a3/nOd3aaZuDAgQ1/l5eXtzivxs81fk1j/fr1i379+mVaLgDkPaOkA0mwLwEA8lFBt2z9v//7v7j55psjIuKmm26Kiy66qNnp9tlnn4a/V61a1eL8Gj/X+DUAANCdBJEAALmpYMPWn//853HjjTdGRMTvfve7uPjii1uc9qCDDopevba/FbNmzWpxuvrn9t5772YHxwIAABCGA0DPVZBh6yWXXBI33XRTRGwPWn/+85+3On3//v3j6KOPjoiIYcOGNTtNOp2O4cOHR0TEySefnGC1AFBYUrpgBNqpUPqeBQCoV3Bh6yWXXNKk64C2gtZ6F1xwQUREjBo1KiZPnrzT848++mgsXrw4IiLOP//8hKoFAACa48INAJCPCipsbdxH6y233NJq1wE7uuCCC+LQQw+NdDodZ511VowcOTIiIurq6uLRRx+N73//+xERccopp8SJJ56YfPEAUCDcPgtdz9cMACA3FUzYunz58oY+Wnv16hU33HBD7L333i3+q+9moF6fPn3i6aefjv333z9WrVoVJ510UgwYMCAGDBgQX/va16K4uDgOO+yweOihh7KxegAAAABAjuuT7QKSUldX1+TvdevWtTp9aWnpTo/tv//+8dprr8VNN90Ujz/+eCxZsiR22WWX+MhHPhLnnntu/OQnP4m+ffsmXjsAAHREobQgT0XLfQXk9Trmc+0AQKcUTNi6//77RzqBI7KBAwfGVVddFVdddVUCVQEAAAAAPUXBdCMAAADkl7QmoABAgRG2AgAAAAAkQNgKAAB5R4tQAIBcJGwFAAByTqrlsbNynu4RAKDnErYCAIkSMgAAAD2VsBUAAMiKVORx81UAgGYIWwEAgJyTLqBG8ulCWhkAoFXCVgAAyDOFkt31lG5Hbn1xQbZLAAC6ibAVAAAgQTuG4X8cKWwFgJ5C2AoAAHmmZ7QHBQDIP8JWAAAg56SMnQUA5CFhKwAAAABAAoStAECiCmXgHqDrpULzVQCgsAhbAQAyUFlTm+0S6MHSPeCqRj6vYh6XDgB0krAVAKCDhs1aGx/69bB4cNKybJcCeS0tlgQACoywFQCgg/77oWkREXH5k7OyXAk9lYgSACA3CVsBAAAAABIgbAUA6CCtCsk2w0oBAOQmYSsAZNm64or4+l0T49nX1mS7FCBP9ITAP5XHiXJPGMAMAGiesBUAsuzqZ96IyUs2x4//Pj3bpQA5rBADvJQ2ugBAgRG2AkCWFZVXZ7uERBVeHAS54dLHX2/4uwBz1530hHUEAAqPsBUAoIOEQGTDw6+syHYJiUu7PAMAFBhhKwAAAABAAoStAJBlWnbln3weuAfoevbqANBzCVsBADpINwJkm00QACA3CVsBIMuMxg0AAFAYhK0AAJBn0j2gebXuOgCAfCRsBYAs02crAABAYRC2AgDJkh0DCcjnxrv5XDsA0DnCVgAAAACABAhbAQAAAAAS0CfbBQAABcagNpCourp0VNXWZbsMOkBf3ADQcwlbAYBkyRggUWfcMT5mrS7KdhkAALSDsBUAAHLY66t2DloNwJTjfD4A0GPpsxUAAMg5KV2SAAB5SNgKAFmmhRrAzvJ535jHpQMAnSRsBQASZWAYAACgpxK2AkCWuVUW6CgXNQAAcpOwFQAAIEH53AUCANA5wlYAyLJCOylPhaa60NUKbb8BAFAohK0AAAAJ0s0DAPRcwlYAIFFCBgAAoKcStgJAlhkgC2Bn9o0AQD4StgJAlul7EeionrDfyOd1zOfaAYDOEbYCAAAAACRA2AoAAJAgDVsBoOcStgIAAAAAJEDYCgAAAACQAGErAADkmbQb1XNa2ghZANBjCVsBgETJGIAkpFLZrgAAoOOErQAAPdicNcXxjXsnxYzlW7JdCh3QEy5q5PM65nPtAEDnCFsBIMuclJNN37x3coxfuCm+/KcJ2S4FAADynrAVAKAH21RWle0SAACgYAhbASDL9EsIAABQGIStAABERERJRXW2S6CdCqX3Ed2oAACFRtgKAFkmbCBXzF5dnO0SoCDYrwNAzyVsBQAgIgREeaVAPivdqAAAhUbYCgAkqkAyICDL8jmITdsTAkCPJWwFAKCgPT59ZZxxx/hYW1SR7VLoAC2tAYB8JGwFAKCg/eyfM2Pmiq3x22ffyHYpAAAUOGErAAA9QmllTbZLoIfQKhcAei5hKwAA5JlC6RO0tVCyuKK6+wqBiJi+fEv8afTCqK0rjO8XANkhbAUAICIKJ8BridaG+WXoxGXZLiFjNrX89JU/TYjfDZsXj7yyItulAJDHhK0AkGWFHnBBriikb1qhBMepVLYrgJ0tXF+a7RIAyGPCVgAAgASlCyUNBwA6TNgKAFmWisJq2iVkACh86XQ6xi3YGBtKKrNdCgDkFGErAAA9ggsBkJzhs9fGN++bHMfc8FK2SwGAnCJsBYAs02cr0FH2GrmtJ3w+o+dtiIiIypq6LFcCALlF2AoAJCplxBugnTQ2BgAKjbAVAEiUW7XzWIF/dIW0afqe5bae8PG4rgYAzRO2ApDz6uoK+6y10AbIAgAA6KmErQDktMqa2jjpljHx44emZ7uULqPPVqCjCqW7jgJZDQCABsJWAHLa2PkbY/HGsnj29TXZLgXIc4V0YUM3ArnO5wMAPZWwFYCc5nQ1//jM8pfPDmg/zZIBoDnCVgAAeoRCagxaQKtSkAppWwMAOkbYCgAAZMWOoeSDE5dmpQ4AgKQIWwGARLmxlFyltWHuu/yp2dkuAQCgU4StAECi5FlAT9cT9oMpV9YAoFnCVgAAyDOF0kpXYAcAFBphKwA5LV0oiUIresAqQk5I94j2huQC+3UA6LmErQAA9AiFFIAV0KqQpzRKBoDmCVuhCyzZWBY/fmh6zFpVlO1SALpdIQVaPY3Pju5WqNucVtQA0HMJW6ELfPevr8Szr6+J024fl+1SgDygz0LoHuIvSI7fLgBonrAVusDijWURUbitNYBkFdq+wgk4OavAvmsAAOQeYSsAAOSbQrtKU2B8PADQcwlbAchpzlcBCpeW8ABAoRG2AgCJ0qIrfxX6oD6Fvn75oriiOtsldLmesB9MhaQcAJojbAUAgDyTz1neVU+/ke0SAAC6jLAVAIAeoSe0NswHM5ZvafjbZwIAFBphKwBAB5T0gFugC5Vcj+6iywoA6LmErQAAHXD67eOzXUKX0cqQ7mAzKwwGNwOA5glbAchpPSH86QGrWFCWbCzLdgnQI/aN+cznAwA9l7AVAADoNhpEAgCFTNgKAFkmeIDukdbcMCc0/hTcil4Ylm0qi8ueeD2WbdLyHwD6ZLsAAGidcARIhr0JJKdxTv6NeyfHyi3b4uX5G2LcLz6btZoAIBdo2QoAWSYAAjoqn0e7bxzS9aTGxiUV1dkuocus3LKtyX8BoCcrqLC1vLw8nn/++fjtb38bX/nKV+J973tfpFKpSKVSMWTIkHbNY926dXHxxRfHhz70odhtt91izz33jGOPPTbuvfdet54BAOQxh3K5oSd8DM1ta5vLqrq/kARtLquKMfM3RF1dT/gEASBzBdWNwJQpU+KLX/xixq+fNm1afP7zn49NmzZFRMTuu+8eJSUlMW7cuBg3blz861//iqeffjr69u2bVMkAANBhgmO62xf/MDbWFlfEtV8+NM775HuzXQ4A5KyCatkaEbHHHnvEiSeeGD//+c/jH//4R+y9997tel1RUVGceuqpsWnTpvjwhz8cr7zySpSUlERZWVncfvvtscsuu8Tw4cPjoosu6toVAEiYFii5z/gw5Ap7C0hGc908pPJ8b7+2uCIiIobPXhsREakCHt2sgFcNgG5QUC1bjz322Ni8eXOTx375y1+267U33XRTrF27Nnbbbbd47rnn4oADDoiIiL59+8aPf/zjKC4ujssuuyzuvvvuuOiii+KDH/xg4vUDJK2ssiZOvvXl+PR/vD1uPPv/ZbscWiDggu7hu0Y25XM/uz2NluMAdEZBtWzt3bt3xq8dOnRoREScc845DUFrYz/5yU9i9913j9ra2njooYcyXg5Ad3p65upYtXVbPDptZbZLyZgTHiAxdig5RwtCAKDQFFTYmql58+bF8uXLIyLilFNOaXaa3XffPY499tiIiBgxYkS31QYAQMdVVNfGbSMXxBuri7NdCq0o1Py7ufXK924EAID2EbZGxKxZsxr+PuSQQ1qcrv65N954o8trAqAHKdCwAbLpjlEL4+YX5scX/zj2rQcLqBlloYaUAAD5rqD6bM3U6tWrG/7ed999W5yu/rni4uIoLS2N3XfffadpKisro7KysuH/i4u1pgAA6G6zVhXt/KCEkm5iSwOAnkvL1ogoKSlp+Lt///4tTtf4ucavaey6666LwYMHN/x7z3vek1yhAB0kV8gThdPYDugmdu8AALlJ2JqwSy+9NIqKihr+rVixItslAQC0S7rQr9AUUDcChaLAt7iC5usEAM3TjUBEDBw4sOHv8vLyGDRoULPTlZeXN/uaxvr16xf9+vVLtkCAHqxHnIj3iJUE2FltXTruHbs4vnfsgdkuJVGFdOHimmffiP59nTYCQHtp2RoR++yzT8Pfq1atanG6+ucGDRrUbH+tAABAx/z22TnZLoEWrNxSHveMXRJ/GLkg26UAQN4QtkbEIYcc0vD3rFmzWpyu/rmDDz64y2sCAOhKizaUxvriimyXQYYKqeVkIWru08nH2+4ra+qyXQIA5B1ha0R88IMfjPe+970RETFs2LBmpykrK4uxY8dGRMTJJ5/cbbUBACRtXXFFnHjzmDjy2pFNHhffAe2VMrojADRL2BoRqVQqzj///IiIePjhh2Pp0qU7TXPHHXdEaWlp9O7dO77xjW90c4UAmUmLTnLK+uKKKKmoznYZEHPXlmS7BChoGh4DQM9VcGHrli1bYuPGjQ3/6uq23/pSXl7e5PHS0tImr7vkkkti7733jvLy8vjSl74U06ZNi4iIqqqquPPOO+Pyyy+PiIgf/OAH8cEPfrB7VwqgByuUE9bNZVVx5LUj49AhI7JdCmiPBrRLT91X5GOXDwDkjoILWw877LB45zvf2fBvxYoVERFx4403Nnn8f/7nf5q8bvDgwfHMM8/E29/+9njjjTfiiCOOaBgI67//+7+jqqoqTj755Lj11luzsVoA5LnZq4uyXQIUpLLKmvjhg1PjqVdbHuQUyExK6ggAHVZwYWtnHH744TF79uz46U9/Gh/4wAeiuro6BgwYEMccc0zcc8898fzzz0e/fv2yXSYAAG+6fdTCGD57Xfzvw69muxRopEBuy2iFHBYAmtcn2wUkrbn+VjviXe96V9xyyy1xyy23JFMQAECOKaSQZNmmsmyXAAAADQoubAWgcFz3/JwYO39jtsuAniMPG+OVVdZmu4SsyOf+rNP5XHw7FcoqFtB1GQDoNroRAChg+Xyyt3xTedw1ZnG8saY426VAwUkVUIRSl887Otp098uLsl0CAECHCFsByEkVNT2ztRrQMR3KWgWzOaEjgy5d+9zcLqyEziicSzYAkCxhK3Sxx6evjImLNmW7DCCHpfPx3m3yWiH12dpTbavO3wtSPaIbgWwXkJDW9hX2IwDQPH22Qhf72T9nRkTE0uu/lOVKAKCHK6B0aEtZVbZLAACgGVq2ApCTCicS2W6CFu7QJXpsy/BC20kWmB7QeBcAaIGwFQC6wZ2jWx7kpZAGKyI/9IQtrtCzrl4F1EoXAKCQCFsBClihhw2Fose2zCPn5OO2WOgtCOvqml9BWSvdobWLgR0Z6CzfFPp+BYCuJWylU6pr66K2hZMAgKT1hEFVoFsUbkbSoFBWccQb65p9XMvW3FYov1c2MwDoOGFrHpmwcGN88Q9j49UVW7NdSkREVNbUxieueTFO+cPL2S4FMjZnTXHcMmJelFXWZLsUADJQIJlWizaVVTb7eC8hWN4RXOYPnxUAndEn2wXQfufdOzkiIr557+SYddXns1xNxPy1pbG1vDq2lldnuxTI2Cl/GBsREcUVNTHk9I9kuZrkTV5sUCZgZ/oJhq5V4NcAYsqSzXH3y4uzXQYA5CQtW/NQqRZ4kLjZq4uyXUKXeOa1NdkuIVGF3oKtkDw9c3W2SyADvmP5w2dFNn3tronZLgEAcpawFQBI3IX/mBHLNpVluwx6gHwc1KsjhKr5yeeWfwqln10Ask/YCtBDnHnH+Kiprct2GfQgG0qa72uS7Cuk/gh7aj7SQ1ebHFCI297KLeXxyWtHZrsMAAqEsJWMFXpLEig0r67YGlOWbM52GZ1irwOwnVZ4dIdCujDTmptHzI/1LhACkBBhK0APUpdH5+Y95QQvomnLvDHzN8SLb6zLXjH0CIX09cqj3Ro9SCE3Siik/Ue9OhcvoFNKK2viz2MW6UIK3iRsBYAccsFfpsT3hk6NLWVV2S4F8k5lTW22S+g2W8qr4/gbR9lX5JFUgVxFFEsCO7rm2Tlx/fNz4/O/fznbpUBOELYC9CD53tJm6aayqMun5rnt1Nz5d3FFdfcXQo9RKKHPjn7y9xnZLiFxre3xlm4qj3vGLu62WuiAwvupgpySTqdja7mLTbli8uJNERFRUW18CIgQtgIZenXF1rjm2TeitLIm26XQg5x485gY8u/Z2S6jW6QK8kbN/LZgXUkUlQvBc06jUGtED+yCo9btzznpqVdXZ7uERDR3YaYQf50KcZ0K3UWPvBofu/qFhpAPIJf0yXYBQH46847xERFRU5eOK0/7SJaroScZOnFZXH3GIdkuI1Gyktw3Z01xnPKHsdG3d2Ffp7Yt5iGfWU56ZOqKbJcABa3+gsafxyyKTx749ixXA9BUYZ8xAF1uwbrSbJdAB+RXkNKz25kU6F3eeWvcgo0REVFVWxi3xxXS9pXv3aMAAIWjvKpGFxcIWwEAyF/5dREpA22sYKGvfiFZsbk82yUA0MUOHTIiPnb1C8Zf6OGErQA9iJNyulNtXTqenLEqVm7Jr4BhbVFFvLZya7bL6FIF1LC1x0u/GcauK66If76yIiqqa7NcUdt66m/ROXdPyrv9Ybrgr2ZsV6iDBvYEPjtyTe2bg/nOX1uS5UrIJn22koh0Ou2Hrody+yZ0Xk1dYX6P/jZ5efx75upIpSKWXPelbJfTbp+6bmRERAy/6Lj40N4Ds1xN1yikn+xMvz2FcuxSn4V96Y/jYmNpZSzaWBqXnnJQdouiRdOWbYn99uif7TLoARauL4l+fXrHe/a0vQF0Ny1bgawqr6qJ2auLekzLCWjOovWF2ffxuAUbIiJ/b/N+dcWWbJfQ7fL0o2qiph396n7vgVfiK3dOiLoCuNBRvwYbSysjImLU3PXZK6adWoq4C+VYoFDWg/y1pawqTrrl5Tj2d6OyXUqXy/9LZkAhErYCWXXqbePiS38cFy/Oye7JYU85L3ICSHfaUp7ffVWlnMLlpYOvHB6rtm5rdZoX56yPGcu3xsINuX+ho629dj7u1lsquQCy70in0/H1uydlu4zE5OP2RcTKLa3vAwHoWsJWEuFAjEwt3lAWERFPz1yd1TqmLtvS0CoIckEB3N0M3WLHi0hVNXVx/7gl7Z4+3+3YnU8+r14hfDbbqmtjypLNLT5fCF1XAACtE7aSFa+t3BrLN+XXAAEUnudfX9Pk/69/fm6XLq+6ti7GL9wY26qyN3hJ/p/GAsnouYFP4zwvH96FtvLHAsgnGxTQqpBn8mFfAED+ELZm0c/++Wp8f+jUgriK35E1WLG5PE6/fXwcd2Ph9yFEbvuvh6Y3+f+tGd7yXFJRHd/56yvxxIyVrU73u2Fz4xv3To6f/GN6q9NBhNZPZEc+HpN0tOL8W8OeIw83vw57YMLSbJcABWVLeVW2SwDYibA1S6pr6+Lx6avihTfWxfLNPauF54L1JdkuARJ115jF8dLc9fHTR2a2Ot0DE5ZFRGS9f1ryX0lFdfzysddi/MKN2S6lsBVw3t0VWf7MFVtjyNOzo6ib++rtaEDXOFAuhGsa+RiQt6SugNalJdOW9byB96ArTV++NdslAOykT7YLoDAGA6DnytfzoiRblrii3jUKIQTpKre+sCAefmVFPPzKilh6/ZeyXU5EbA+AyR8tjgbfiXmeccf4iIgoq6yJG8/+f52YU9fK05+tdiv09aN75etxXoc55gAgQcJWErG9VUX7jlJ6zEEbOe3Kp2dnu4Ts8P3LSR39WHLxjoiqmrpsl5A4596Zmb++NNsltKrpcUj+f8o77j9yuaVrcUV1XPfc3FjWQr/9OVw65L31xRWxS+9esceAvtkuBaDg6UYgSxxMQu7p8paUOXBOv+Oo1TSvri6d9ZaarW0uWv3SWS31CdwTjk/ybT/YVniaT5/ZzcPnxT+mLG/x+Z7QjUAh8Cnln9LKmjjy2pFx2G9eyHYp7KCiOnsD5wJdR9hKxhwPU2hs09T71l8mx6FDRsTSjWXZLoUsahxI5ltA15N0eICsRi8ohAsXO26bubylLmmhRWu9XK69p7Lvy3/pdDqWt/HdIztuHjEvPnz5sJi8eFO2SyEDpZU1UatPSFogbM0TuXxLWETHDo5zfFWyZtmmsrhj1MIo1u8h7bRwfUn86onXY/XWbdkupeCMX7j9oPex6Su7ZXnN5T2thUAFkA+Rs/LwR7oDBxaF+N0ppOOqXD/eZbtC/B5BNtz20sKIiPjNs29kuRI6amNpZRxy5fD40h/HZrsUcpQ+W/PEo1O754Sf7Dn51pejsqYuFm8oi5u/lszAIhtLK2Ph+tL45AF7tnjLaK7Ihepy/C3ayRm3j4+yqtqYtaoonvqfY7JdTuLy7OPoFPFC7kv1qC2ycKVjh5atWaskOYW0/yikdSlkPqf84hoGJG/0vA0RETF3bUmWKyFXadmaJR29JeevCY6c3hUy/RFfuaU87h+/JMqrapItKA9Vvjm4zCtLNyc2z2NvGBXn3D0pRs5Zn9g8yVzSJ/VlVdv7eHp9VVG7X+OAOzd0tvVWvl0Y6ElqavNjoLCWNqF83Ed0uBuBPImKKqq3X0xrq9p8/Mxaks6Pr0+PUkjbV6t2WM9ZHTi2AoAdCVvJWBIn+1/647i46t9vxHXPze38zNjJtjc7XB81T9jKdj3mpCnPlFb2zAtOpZU1MXTi0lhfXJHtUprV0Z+5WauK4sOXD4vbRi7oknqS1FN3Bfl0neKCv0yJU28bF/9s8+6mdKv/m0/yJQjv6fLpe9Qea4sq4vEZq5o8NnlJco0fco3uOiBZz72+Jh6ctCzbZZBjhK1Z0tHfuFz8SUzid7po2/b+Sccv3Nj5mUEnTVy0KS9HBM31LiJo3cL1zd9+1No+Nhdvac9kO/z1E6/HFU/Njq/dNbELKup+V/17dtTUpePmF+Znu5SM5eLxRpJ27EYgl9WHPXPWFLc6XV0BtQbNl8+mNYWwDo01tzoFtopxx6iF2S6hS6WjaSOZQttGIdv++6HpcfmTswysSxPCVhLR6ZYIuZcb0AOVVtbET/4xo8vmLxOlXuMTnalLt2SvkCyr7+JkaY6Oktze7+zW8qquLaQLFFLLps5cwC7Ei1X5/Mnmc+35avLiTXHd83Oisib/LjbTMq3EofttKa8qqOOrllRU18a6HL0rLZcIW/NEIX1pC2dN8kdXnksW0KYZEREvvLEu2yWQQ7pj+65rYRmtLbpQ8qFC2X188tqR2S4hMYW2T99RKgrrmCoiv0KVtt77ugL4bPJtDb5+96S4a8ziuG/ckna/pkB+gmLYrDXx479Pj7IC78pnx+9dvm2jPdnoeevjm/dOjpVbcvOidL3iiupYrFVnj3Hs70bFJ68dGctztLFEruiT7QIAck1tXTp69yqUU4mmHGBD+7U31K4f4DCftLQvyKfgrl5Hal69dVvTlq3Jl9Ptdswn8zlMzuPS815Lt782tz0Vysf0o79Nz3YJWbH9My2EvV/h+8/7X4mIiP/712vx9+9/KsvVtOzcuydlu4Ru19JvbU/4HdtQUhkREWMWbIhvvf19Wa4md2nZmiW5+iXM9AA9V9eH7XKxf0eSk06nY97akqhrqYkkOSfd5O+WDtZa/jwLpWUrJKGlr0pzD68vqSy4Y5ZCWp18DPt3lM9hdy7w/iXPO5r/6sOtXDV7det9i/cUO3ZN5Hi9ZxO20uD65+fGJ68dmZWdea7th2Ys3xJn3D4uXllauCOR0rJ8O9CvS0d8/vcvx1X/np3tUhKVX59C5npyRp5r+/4dpSIVc9YUx4rNbd8mtXrrtoK4sJVnu7+IiNhc1sE+cxutYyGcCHXmM9tWVRszlm/JnYt1OVIGreuqr831z8+NY24Y1fHvNDtp7feoEL5mM5b33P7uyV35dg5J1xK2ZkkuXrn/85hFsb6kMu5+eVG2S8m6r989KWauLIqz/9z9I2QXwolfvsuVc86OemDisjancRCQgzK4DakQQr2I3D/hW19SEaf8YWwc+7tRbf5uX/7krG6qKhmFtCtYU9T8IA0tfUty8RgsSR1Zu2/dNzm+/KcJ8dCU5V1WT0fk6+9vYwWwCk105/r8ecyiWLV1W/ylA/3H0raduxrJTh1JuuTRmdkuoVs5P8wfBfD1IiHCVnaSjR/gXBsNuCqL/e8VwgFQvuvoiXh7N9/GAdktL8yP028fF+VVhT0oQ0TE2qKKmLas4y0QCvm70Dj0zmg1c2uXWbAWrX+rD8MX56xvddpNBdISq4C/dg0a71sK4cJFZy6iTX1z3/zIK7kRthZ6EE775Mp2sGhDabZLAPJIrmUaXa1nrW3HCVuzpJBDBHJPD9vvd1p3fD//OHJBvLayKB6durLrF9ZId+96qmrq4lPXjYyz7pwQr63c2s1L7xpJtw7OZHa+0t1v4frWT7rzr9V4vtWbnEJf87zbFBvJ59rrFcI6NHbLiPnZLiFr/j1zdbZL6BK5EmZ3Rs8LtXrW+na3ovLquPWF+S0OFNgR+Xc8SFcRtrKTTHYPnd2n+Pl4y47HDul0Onf6UmvGxMWbcrq+XFddm3+jmHfEr554veHvqUvzu3+tdDod/3n/lDjn7kmJHkg5KMtdVY2+n73a+KEqlN1gT9geG69jIZyvJ/GJ5cLHXlReHef/ZUq2y6CRom3V8ezra7JdRpdYsK6kzWmWJBC85IJ0pJvs63Lh+07HFMJvVS677InX4w8jF8Qpfxjb6Xkt3dR2P/+FwnbZOmFrluTyb1xdF/8C94QTuSSde8+kOOnWMVHTiVCuq/eDI95Y28VL6F49eROtqqmLV1dsjdqEkqNHp73VcrfjP8i59UFsq66N0fM2xOQlm2PV1m0dem1lTW2T/0+38HeTaVrrs7VAjm5y/ffgiRmrGv1f6+95V/92dkRFdW38Y8ryVrfTHCq3S7T2FemOVb95xLxO9Tv598ntv61/x+9RZ1qtdfV3srXZ//GlBW22IM8LBfTd6syxZ6774YPT2pzmqVcLs2VreVVt2xNBASqrrGm2kdCUNwfF3lbd+e/GBS4a8iZha57IxZOiHCwpp2V6AjNp8eZYvKEs5q5t+wp8Z1TV1MU5d0+MG4fPjYiIreVVccuIee26ql+0rbpLa+tuXXV7VS7kY21thj//18w4847xccsL87p02bkesiXtn610F6FFZGHIpc/xtpcWxKWPvx4n3zKmxWnWFjc/qFRXqqiujUdeWR6rO3ixImldvaku2VgWt720MK5+5o2M53FZo7sC2pLU6kxctCk+cc3IGD47OxdQt5QXRr/HhSSHdmuJ21hame0Sus2O+zwDkNETrS2qiI9cOTy+dlfnB8CetaooHnllRbPPdbRBBoVL2JoluXxSmklphdD3T083bPbamLR4c9wxalFERPzysdfjjy8tjC8mcDtFvtkxFOzKbhJyrYVifSuOP41e1GXLqKiujc///uX4xb9ea3W6XNtNNldPTW1duwbUK2olRGjp96C1/WpubTWZy7XtvzPS6XTOfDBjF2yMiIiyVlov/fKx9od5SbntpQXxi8deT+Q2vaR0xX6mrLJ7Bz5MYqTxdDrivHsnxcbSyna1+KNljonJdcUV+d9IIkd+brvN3LUlBd/1WFd75rXt5zhTMxi0d0en3jYukflQ2IStZCzTHzmHoLmpcWBUVF4dU5e1/3aKQuu0vfE2+t0HpsbJv3+5gA5wsvcNrM/VRs5ZH/PXlcYjU5u/IpwPUqlUpNPp+MyNo+NT143s8PaRa0FyEgpwlZqorWv9M86lbgTas0cubSEQ7MrVGDN/Q0Rk/26IQgvDklqbJD57fbgXjnQ6nZWLMnS9dDpi6cayrO+L6Zh/Tdv5Lqmbhs+LO0YtzEI1QFuErVmSy4ei7W11m+Q6FFDjphZ19iSmO8/jR89f330Ly0GNQ5OX5q6PhetL49UVW7NWT01tXYyauz62FsgtlrkUSnVGZU1drNq6LTaXVXXJbdGt99ma+OJ406xVRc0+vqW89ZPSnNqsO7GBFFoQ2azGdy8UwPrmyj71jyMXxBHXvBjLe9DgIIVkxwvnU5ZsjhfnrMtSNXRWa/u2pZvK4vibRscRv32hGyuis9YXN+36YvXWbXH7qIVx4/B57brLCuhewtY80Z0nA5ksKUeO8wtaZ7aBpG7VLa6ojp/989V4+c3WSW2Zs6Y4/jR6YVQk0Nl4d8q17fn+8Uvj2399Jb78pwndtsxceA9yoIQmkqyn8fc5V4KSHc1aVRRf+P3LMWpu+y6+ZHbrcu6t+4rNmQVF2focX19ZlGgfmz99ZGZiA+TlqkJbu2df6/xo8Um8J7e8MD82l1XF9cPmJDC37lFXl45XV2xt93FKTW1dLN7Q+iBeObhby0ghD6JUU1sXxRXt6+5jWwG8DztukxMWbYqIiOra3NtYV2/dFmuL2u5TvLlTmx8MnRpn/3lCj2lhn2/nVxSeQru7NWnC1izJ9oHYqq3bWu4nMAu1dfSLOmHhxhZbH7XmjdXFsS4Lg4JEZBpi59bBwq0vzI/Hp6+K83ccZbGFj++UP4yN3w2bF3eNWdzmvHOqpV4zb3s2P4pnXt9+It2ewcpyWf1HnMRb+Y8py2PUvK5vgd3SRY4kN9eWtq1Ew910usMH5d97YGrMXVsS3/7rKwlWkvsyvTiVrX3EabePix8+OC3mrClObJ7z13XtgIzZ1rRf7uzVkaTGAWC21ynT5Wdj4LT7xi2JM+8YHz9oZz+13xs6NT5785h46tVVXVxZ99vx9661XWE+fW1mrSqKTTsMhvX3Kcvb/fohT89uc5rq2rq46t+z46W5udkSeNGG0nhw4rKG/8/VC2oV1bXx6etfik9dNzJqOtw9UzpGvLEuXlm6JRa2cUEkH1TV1MV590xq8tiO38nc/BR7tlw6nSX7hK15IskD5wcnLYujr38prvp386PkZtKCsjt39qu3bovz7p0cp942buc60um4Y9TCZlteLt9UHl/849j45LUju6y20sqauHds28FiJrJ98hQRsWpLZidCr2cQjDcsc+u2OPr6l+KuMV03YNOOmvsOJBF858IPcC5sR43dO3ZxLGrhoLi1Wi99/PX49v3dGwK2tg105n3N5KUd3Za+ed/k+PDlwzo0+nJJBwfQKITbsSMi/j1zdUav6+6WrS/NXdck8Fna6GJMZ/c1beXNizaUttq6b/mmspxoWdRSBd16t1Ar28XG0sqYv64ko4vHOyqrzJ0WTq19FVp77yct3twF1bTurxOWRkS0+46d0fM2NHldc7K/5SejV05dBc/MzBVb49TbxsXhv32xyeNz1rT/gtJj03fuJ3NHD09ZHvePXxrf+evUDtfYHU69bVw8NLn9AXO2bCp7q7usig7eFj9u4caky8mq52etaWiBXK+1wRBnr+787wh0VAH8THQpYWu2ZPFI7Prntt/e1dKBYnvPF5P8bnXki7qylcBv5Jz1cePweTu3vIyIWd3wI3TNs2/Eb5/t3O1zjd+KJq1vOjXXZJRVde8IyxERNw6bG6u2bovrnp/bbcts7juQzfc/0VaUCc4rU43X57fPzokTbx6TtVq62raq2hg6cWncNGJ+k8c707Kuti4da9pxi11j4xduP2B/flZyt5t31Etz1+Vca/3mPPt6Zrdkd3e2+J2/To3/ffjVZp/ryoPfiuraOPHmMfHZm8dEZU3zAd+W8ur41ZPZHVjnnrFLGoKxHXXnb2tLm/xLc9fFEb99MU6+9eU49bZxnR4dPJcuduRq1yjdJV/Xf8e7zFrbj7S1i8mVff34Rd0TwK3amp275lqSTqfz/m6ojvrWfTuf++Wzyg6GzUl0NVZcUR13jl6UcXdKQFPC1jxWW5eO+8Yt6XCLiLZukWzp8GjHViodPYxKp9NRtK262ROPuWtbvsK84wFbr1bKX7ml5R+H7rjwMnZBywd19esxZv6GdvevlxuHqm+pD2x2lPR7m06nY8byLVFWWRM1WWgd1VVLzPT25MYve3rm6rw9CKpf/1zbrturtbpb+mh/8+wbccVTrd+C2FJI0tLJ6vceeCUmL+n+FmBtauOD/c5fp8Z/PzS9e2rJQGcHl8h2wNJdrQtKGvVz2Fpryn9MWdEd5WQkF/ZBfxjZdPTojSXtb3XeFZIMx+p/tq959o24/MlZDY//e+bqFo8jCskT0wuji4HWWra2tbXkSNYavVtch/YX2J4pc+liR8T27jFauhi3o2WbcieU7cjPmL4iI5L+NbviyVlxw7C5ccYd4xOdL/RUwtYsafyj/LW7Jsb3Huj47bCPvLIifvPMG3HqbeNaDRl31PbV6OaWtTwOGTI8Ji9u/iC5uh0nqf/1t+nx/64aEY9Obf8J2C0vzI9PXjuySUfp+dxcva4uHRf8ZUr88MFp7bqdt/HJT1e3EtjxbW28uEmLN8VtIxd06fIbe3z6qvjynybEWXdOSGxwr3ovz98QP24j8Gnuve6qt789a9d4mgv/MSOO/d2orimmi3V0G861k5dWb41t4bkX3mi7/7aOblujWmitl7SyLhgU5PlZaxPtl3He2pIOd3fQnBffWBcf/PXzMXTi0ozn0Z3hQlvfpTz+mewW3fnbmo29WLZbFabT6Sivqol7xi6JByctaziG+8k/ZmS1rkzVr097TV2WgxfDMlAI+5FC6AohEzfvcDdNaz5z4+iuK6SDcuuoL/fseFyc9K7+yVe3d6O0uVF3DkDmhK3dpP7Ad0NJZczboRXnhpLKeHHO+lY7Am9uX9p4MIxjbhgVD3egs/fW7by0Xzz2epRX1bbYKumw37zQZng47M3WnCPbObJ1RMQfRy6I9SWV8ceXGgd9mV1pb+t465YR8+IX/3oto5OU9rwmHU3rK9rW0f4QW/fc62vitZVbm32us8ea59w9KW5+of0Hbp1V3z9Way2eM3X+X6a0eZtwc+91rgV/3eXulxfFjcO7rwuHfJFKde4gt7KmrmG/0eJghZnPPivaW29nb5euN3Xp5vj871+OE24a3el5/ffft/+2tdUKuS0987Q+/2QrGF+6sSxGztl+AWanC5ydXE4utfJKR9NuNWrqOtdqPNt+9eSsOPiK4fHqiq3tmr4rt6+VW8rjudfXdEufyK1d7G6z4UaypWSs5VVo//elXVNmcYVLK2ti+aamjW4K4Zg1FRHXPjcnvviHsbGtgxd/c2dv2LXy/1PuvGxfXITWCFu7yflv9iPziWtejM///uVYumnnlqj1u4rlm8rjz2MWRWll61fRdzyAuKW9YVgbv0CZ7rO667apxuvdsR1s6yv+x5cWxiNTV8T8dR0bwfLO0Yvik9eOjBWby9sMNRvXu2JzeZzyh7HxxIyWO95v79q9tnJr/PdD0+P02zO/7WPHZXUkoE269WmTeXfZnFuW7duB6z316qoOtQSPiDjjjvGttsDp6Kpd+9zcuGPUop0O5DPR0e0kRz6Gt+w4MEE7vqEtrfH/u2pEQfUvNm3Zlhg9r/0X0pJQ3x3LxtLcaIHRnQf8bS2q8XetorrrBk4qhBPaLu+ztdHfx980Or77wNSYsHBj4nfp7Lg/uvTx1+PSx9vfb26Sm2+u/IYm5e9vDizU1h0+5VU18eU/jY8R7bijIVPH3DAq/vuh6fHUzK455p65YmvcOHxubKuqbbXbrra6XcmVACTplq0lFdXxwISlsb4kd/poPerakXHcjaNi4fq3zl+yOT7hptLKuPvlRbEhga5R7n55cbyxpjie7qLtPZcVwu9rV6utS8eX/zQh/utv07JdSsYqqmt3aohH4RC2dpMZO1wNf3X5lp2mqT84PeUPL8f1z8+NQ64cHqfeNrbF4CTTVhFdtfPuyquoqRb+7sjBRHuPtzrab98Nw+bG+pLKOPZ3o2LF5tZvj21c7uVPzYo5a4rjp4/MbHn6dg6i0/gAqzm51OIlLzTftLVZFdW1DSdiSaqurYv/ffjV+Pm/Xost5Tu3BPzGvZOafd3MFVvjn6+0HNB25Hva+GSpooWBcLrCis3lMX7hxhiWxcGcmpP0Pq5+5NyWvttdfa6a5MnwWXdOiF881r5gJ9cG/ktKd57ctrWoxnv8TO56qalNx1Ovroo1Re3/TWuPXPktynYOtOMxYUSyNW0ur4p/TFke/5iyPLZk4XbQHb8LdXURW8tz46JIV3p4yoqYsXxrtyxr8uKu6argjDvGxx2jFsXtoxa0eoF0eZ70Hd9aYJyJXz7+elz59Ow49+5JkU6nY9aqoiiuqM7qb1nJm41zxi3oni6G2vJfD02Pa5+bG98fOjWxeTZ382cP7SGiiWz/lnXGtqraTgfys1cXxasrtmZ18NfO+sqfJsTnf/9yu7ody0U7fg1XbC6P028fF0+92vMukDRH2JpD6neYjfvJm7WquNngJJ1O73QQ1N4dbpsDZGW44+6uHX7j+ndsPdHammX7NzmdblpvR2+Jae20Non3fsf3J8nWYp05IMrGwVQHstb4w8gF7Q5ZmluVltavttFMy5pp5d7aICNJDSqWjYO4dDodx/5uVHzj3snxh27sJ7g9Gr8fqUi16/1pa5LVW7fFv6a33Lq9o0bPWx9Tl7Z9En7FU7PiuBtHtbu/0yFPz87bgdnaJYFtvbkwvqyyJm4cPrfDA1m2uawOfDkz6Xv3L28OrnLSzWM6/NrW5MqtrY3rSKcjfv/i/HhgwtKuWVYzq5xKde0xSePfj2y0Mt1x+zzv3knxsatf6PY6ultVK92BJa0zH2t7+p+du6ak1eOvVW30vZ0b3/SI3gmnrS++GYgs2lAWY+ZviFNvGxefuyXZ/WQisvgBTHlzAM/2drvRHjtui5MWb+qSrsbyTa78pmbiE9e8GJ+45sWGwHXxhtJYurFjg7XlYtjc0YsMb7zZLeRj05I7F8imy5+aFa+tLGr3AH2FTtjajZoMyNDM8y0dEM9cufNJ2ieuaTpoVMtz7bhMD8y7cn/X+Ee28e9tbTbvk8lA07e27QPA9v6IJn0y1dGra5kMupbUvLvCLx57rd3Tju2GlgTZunqf9LerPevx/aG5eyvQju9HEu/PiTePiWUtdtHQsSWsL66I/7z/lfjqnye2Oe3QictixeZt7T64++uEpfGt+yZ3qJ6WJLU9J9l9SRIhSXP7uRuHz4s7Ri2KU28b1+n5N1lWs4++9X509q0ZM3/7fq25oLbQWhMt2VgWv39xQVz5dOf6621Jc7/jzbfwTW6Pm+0T0Lp0uskartyS3KB43aWjg2JFdOx4ZUkHQ4UdZRqy/HnMojj4iuHxfBt910d0rlVotrfBei3/TnS+wPq7b9YVV+ZEtwkrGn3PujKEK66o7pb1bW0J59zd/N1d9Z57PX9bOrYmBzazxNR3lzht2ZbYVlUbn715TBx/0+io7saLVl1hfTta624qrYynZ66Oym68a7Cz6urS7WqJXFrRsd/NQids7UaNd5BX/fuNnZ5vKTd8YsaqnX7UNpZWNgw41VFdeaI0f11J/PdD02L+upK49rk58bW7JrZ7p9naD3eqhZPIHV/S+H93bA2Yab9NSR1QlFbWxPg3bxuOiDYHFNu+7Ob/3mm6zhTWjImLWm412ZzHW+l3trPaE6hUVNfGGbePS2wgp9FvjvbedMTqRGbdpnQ6HX8es6jDn0GmXp6/Id5YXdz2hI105UH2i3Ny7zaadDod/565ukmfSjtulpnuV7e10p9mR9/m9hzg7agj16ua62uc1nX0u9VebfbZmvV7OZqXK3U1fv+KOzhYZRJSqeT7Ou/sbjnJcKYrx8O6/MlZ8X//mhnbqmpjxOy1HQ5E2+tn/5wZB18xPOaubf47PGP51k61RBo6cWnGr+2M65/ffpx0yaMtd2EVsf24siv74+8ufxm3JNH55XLWdd+4JV3eCGXS4k3x0SEjmu0Puqa2LtHjw8ZdoDTeElsbULrerS82P47J8Nlr41v3TY71xbnT525H7HTRP5c3yFb8Y4fujbY06mamsoPd+eWjc+6eFBf+Y0bcPKL7Bp/urO8NnRqfuObFmLS46flpUudDhUrY2o3a2h92tnViu7sRaGs+GS5/1Nz1ce7dk+K519fG1++aGHe/vDimLNk+WvRnbxrd5uvbU39lTW0Mnbis4f9be8+O+92oJv/f2pe/tYODpPqNPO+eSfHdB9q+taClg9vW3p62Dm66esfX2i3tmWh6u3bbHp++KmauLIo7Ri1KtI4mNWX4zZi6dHOcfOuYmLBoY9sTx/YDweufnxvf/usrGS2vLY3f20UbSuP8v0yJL/5xbAvT7rzOFdW1cfKtL8fF/2z9RK2QTF6yOX7yjxlx7j1vtaRIp9sXxnfmq/e5W1+Oacu2tHv6TH5C8vQ4PSJyZwCWet1ZTUf2R5ns/9s7945+Brlyy2PjKrr6NvvmZn/983N3ujspyTKam9UfXlwQ5/9lSre0Gkp30SddVlkTD05aFv+cujIuuH9K/ODBafHzR3e+E+Wfr6xouJU5U0/M2N7f3L1jWw7rLt4hsOzId62lz7ussiZueWF+zFlT3KWtntrqbii9Q+vkjsqV7/riTrYgjmjf/jBXfo6OueGl2FZV22X1/P7NEPPhV1bEqyu2xrG/eymGzVoT5VU1ceS1I+OC+zt37Nr4N+V7jc6Z6r9b/5y6Ij58+bCM5//DB6fF2AUbm2301N3S6XSrt803dz6YTqejri4d33vglbjiqVk5s911VOOwvieGcwveHGvl2dfavsMgV7w0d/sguPePb/0CVq5cVM8VwtYcku7k8W97Txg602drSTN9R9abvGRzbHrzKmTjAX0Wbyhr18FOa/XXl/zn0YvjX41aErT2mk3tGBTi9ZVF8eO/T28ysNWOB4jz1r3Vmq28qib+MWV5RldE29u3UCYBQkdeUluXjvPumRSXPfHWD92GksqdTho6q/F61HSimcvjM9ruYLuqHScks1YVdfhza6ll8aINpc32o9qcr989KeavK43z7pncruSt5VvKO+aKp2a1edv3kg2tfy+b26xGz1sfC9aXxmOd6Gc0myFZXV06vvbnifGDDvSp1NLBcHesxf/8fXq7p61rEv5u//vif86MM25vegt740EAs/FZtLXI0sqanAtS2yOdbv2kYfbqorj++bkZHVyXV9XEvWMXx/IM9g/ddeB73fNzumS+FdW1MWtVUaLbRDbuWthRW31edlRbJ6y3vjg/Xp6/oVsGHqxLJ79vSafTUdtonvVh6rM73A4/bdnm+L/HXouv3dV2dypJ1bVwfWn8+snXY/XWzFrLjZq3Pq59bk7U1NbFLS/Mjz+OXBCn/GFsfOjXw+LPY5q/iNzZt7etFpDpyPyOsIjcCR+7UpM77Zp5ftqyzbFgXff2K7qmqCLGLtjQZccnjbeJ7z3wSqzYvC1+9Lfp8fL8DbG5rCpenp9c11qNz7/qf8f+71+vJTIuQXvuLuxqfxm/NI5vR2OkxurS6Zi1uihenLM+hk5cljMXNTojqSOUpH9z5q8riT+NXhgVrdyB1ln5eqxL+wlbu9FDk5e1+nxnW1c0N2J5JlrbcVfV1EVdF92i0tpc63fEU5Y2bUHZkVKaO2Y87fZx8exra+KHf2u5n8jGH8s1z86JSx9/Pc7684T2L7gT2tONwI4nIPU2NTqQaLzq05ZtiQmLNsXfJ791C8cdoxZ2ttSdNO4DcU2GJyDt1Xjtn3ltdXztrolNgtX560ri1NvGxZHXjsx4vvVmLN8SJ948Jk5q54AIbZ3QXP3MGxkMlta2oROXxdgFO7em7ci3t7ltLtd/ZKtr6+KJGStbHEV98cbSmLJ0c4x4Y127b7cb0K/PTo911wFuR97vxpPWv+6x6St36ve7udv/dtSVB5etmbWqKA65cnhc9MirWVl+Z2wsrYxJO4wQ3ng7eWL6qvjzmEXx4w4E6PVueH5u/PbZOXHy77fvd1oadOmt/2n+8eKK6rhz9KI2Bzvr6GCTq7dui7vGLG51nttf2/Kct1XVxi0vzN9pMLEL/jIlTr1tXDw6Nbnualpq2TpnTdNbxqcv3xLLNnW+ZVxHa+qsxu/yjvNt6bud6L69C3aP//3Q9HbVmNQFy+aMfLNlT2Nf+uO4OOmWMfG3Scvjrx0YZK3xSfa3738l7n55e2OC13fY/utv+9/p9e1eUvPabtlauC3Obhu5IP4xZefBh5O0tqgizrpzYnzu1pe7dDnNqUunuyzEabxNVFQ3vnDbJYtrtOBkZ1e0rTqmLduc1bDr5hHzMnpddW12D8RLK2viwYlLE+uKoSPdlVTW1LZ4/JL0R3nyrS/H74bNi9tfSv4cuV6On1I1q82aC/R3I1PC1m50xVOtD77QWtjaVQfhmSys/gAt+VYLLT/X0o74n6+saLWF4cYmgWPLaz6/0dXTf7XSB1f9wFGNW8ImrfG6Nh0xufk36Fv3TYlfPTFrp8e/cW/TVo31fRw118q0qw822hNMdep2tUaz/5+/z4gpSzbHb559q5XV9A7cit3iMt78b32roDU7DVCX4XzTEX8anewPeVL9rN2XcF9n9bqyH7h7xy6Jnz4yM05sYRT1phcw0lFcUR2/+NdrMWFh8908tHyRo+0Du1tfmJ9RP6qNteetWl9cEcf+7qX4U6OLJq2V1lyr5JrauiYHsO3p+iVTrb1vd728PbB76tXVbc4nqQuM3WVDJ1rS1Ie4jU9uW9NC7hq/fmJW3DBsbpxxx/iMa2nu42tvH2ut/RbcPmpB/HHkgp0GE5v8ZgvGli5YZ3IBuPE2WNroGOLrjVpDLt1YFl/504T4zI2jOzz/lpbVlXLpQlhdOvnLUc93Q4vcTLyxJrN+mZt7f5Ju7dwZ6ej6VvEbSirjjyMXNDPg7/Y+cbvqot/NL3Ssn8SW3oXGge2O37/lbVzQaq/Wjs9beq62LvkQp6K6NoY8PbvJeALdmamkYntjh6TMXVsSZ905MUbO2fkCSnfJ5P3b8ecuG/v9y5+cFZc/NbvNQco64vvtvNvstNvGxbG/GxXT39wW2mpdnoSZK7fu9FhbF6zbK5d+t9ur7TEDaEzYmkO6uE/zZjX3Q91WGV3Voqu1+f51wtJ4dcXWnR6/5rk5cflTOweN9Zrc/ttkYK2my2rcwq1xn7D5YFwLIVHjbgsmL9kcH/z183H/+CU7HTxX19YlMsr6ju9p01uVOzav9gRMVTV18bNHXo0nZ6xqtt6tjTpbzzTbe+a1twKf+vXriq2/s6MSd5Ubhr3Voqa+xVdnctL6l3ZluD963vYD5/IWWgs3rv/1VUVx6wvz45GpK+K8e1vuciHTwYz/MHJB6xO0w46Lbm5b+f3IBbFi87YY8cZbg4vtOPhAS+o/iv+8/5U49nejYuSbA5St7sDFhI5+nmfdOaEhIMt0UyitrGn1wljEW/2IJ3VQnE07boPNvW8XPfxqTF68c//ZjV9b33f05ja62WlPNzzt8dLc9g9419ZgYs1tKjW1dXHy7zNpOfbW3H7b6MJccUVNHP6bF+Lpmas71F9y60tq30Z+8q0vJ3aLf+MlZuPEZ3vLui6YcQLzHDlnXdwwbG6X3aXVXrl450iS10Hbsy7//dC0uOWF+XHBX6bs9NyERZviBw+2fNdZprqrpXoS72V1bV2c8oex8V/N3H1328gF8anrRsbqZgL6rvj+3fPy4vjrhKVNz1W7eefy5T8lf1fhiDe65iJORXVtTFy0qdU+sjPppqOr+xhvjxffPNZsTxeB26pq49UVW1s9TiyvqonZ7RxMdP667X2dPt3MBfn6ZWwpq4oRs9d2af/kx+4wLkymcuHz7LimNe+YKzTerHOhBXm2CVtzSGsb4uI2+lbsiMYnUYf/9sWdrqS39YWof7q7R9I9s4WWOK31f9e4xsbVZnqM3d27iiat8Doxn9dXFUVdOnbqEL66ti6OueGlRALmz948piFg/ePIBXHokBE7TTNq3vr4/YvzE9npPvLK8nh8xqq46JFX2x4gLMMjwv99+NWdHutM7Zl+Yzpa/4ZWWlNmWn9z70VjlTW1cd+4JbHwzU7fc9db7+VPH3k1Vm1pvSVRS7u57toX7LiffbCZ72ptM7eU/frJli9CNVYfBNVftGnv65rMo4Nvxrbq2ljTwu1n7d0+xy1ouV+4mtq6qK6ti2//9ZWYuHhT/OKxnQfQ6S5JHV82Pin7++TlzQZ426pr4+sJtjRpSXO/Sy3toR6b3rTP7db2ZaPmtd7XX3MnJfPXlWa0z2ntc9lUVhUX/mNGk37MH5++sqGV3bRlW1oMYmvr0jFp8abYlmGLvB+10qVRR2TU93siS96uaFt1PDq1a2/Tbo+F63fuL/O7D0yNO0cvimdeb1/fyV11jtjcdzgd7R+Uqr11dSRU3vHOj67uRuCVpdu/R/Na6Nc0yf4/63W2pXpLdvw8k3jrXlm6OeauLWm2VffNL8yPdcWVcUszrXS7IsBZ1sxFy8breOHDMxJfZpNl5VmfFv/78PZBVW8a3kpXAZmsUp5lVt+4d1Kcecf4eOSV9v8eZHqeUv+qL/9pfPzgwWlx98ttd23UHp3d9q59bm6LfTfn4sc5a1VRPDhxaYufQ9stW996v77w+5fjrDsn5uydKd1B2JpDsnGRfXNZVRx9/UtNHmurjEUbSuOqf89u9rafzsj02KDxPnDqDidALQYlHVjY0k1l8bthc2NzWVW3jOLbWOMqu+KAf/66klhXnEwn8Us2lsUVb7Yybu7gL2J7v2S/f3FBDJ/d+k63tf5p63Wk5VUSx2j1S87GxbmO1t/SgBqdne+OGt/W++fRi+M3z7zR0JftU6+uSvSWr6Q0XueaunTs0vutn8GWRn5ucXts9A0du3BjXPCXKTl1G2h77Lhu7eke4/InZ8W5d09quCMgk69Ec5vea83cqtWcl+dviB/9rfl+T9PpdBx/0+g46rq3fteyNRDGko1libUw7NXoaO2yJ15v+2C3cZ+tLXYq0Hlt1ZHsoFaJzapd/RY39rN/zozrn58b26pq46w7J8RZd05o9hbnP41aGP+/vfsOj6Jq2wB+b3olkNBr6EiTKl2NoML7ih0rUpQilg+UIjb6K6CAHVBQioI0QaVICyAEpITeSxqQEFp6L3u+P8IuW2Z2Z3dns5tw/65Lr7DTzpYzc+aZc57z0o/7jWbRdsU1Q61DJqfnIaegCP9cuIken+/Afome01JibmYb9RhWi1wv4Yy8Qskefr3m7JbN331d4lx3JikD/9t4BukuSlGi9m8lt6AYD8/ahVF2BsIcKY87TtxTmr2r1Gh3Gmb9+r/fjkqO0pAKrCrNRy/n6+0X0fur3cjIs1wPDINQzs4jWrZCrcCW0yW9PxftjZddx1rPVqmlpt+3SzoMmhRMqxWy7bcjl0te/82WYKuSIlj46OLv5O3++5Ttk5FKHsvB7Q8npMqOwHHHHp9PfBuFT/88jfUyndmsldjwu9HdUzDYSm5hnsp5G+1lrd4/8/0+LNobj7eWqTu8x7BhlnA7G+HjN5qtk5ptfuHPKyyZECe/qNisl6vhCdLwomZLO+TPY0mYuysG7aZuQ5pBA9xaI6Qs+O83UZKv23thWSFzMb1o0vPI3hl763+4CcN/iUZ2fpFNDQxVnojfOZ6S305WfpHVYbquFHszC//G3Lb5ZuD3w1eNglyGuYIPGwRWj15OxcgVx8yHfLmgZ8L7q45hvEHPRsPzgBCAt+fdfzf9ZDO+3q5s6P/640lGv8FP/ziFfy7cxAdr1O1F6eyPzJ5m3i/7E/Bv7G0ciCsJvDjU29vg/T353V5F5bHUEz+noBhXU3MlA6w7z98o1ZQdKdnGZXDkPGR6U2bpcyoo0hoFWI0Cr6qcCoXk31IkUoTbf1wV70lMH8wqsenkNWTm373u5xYUI7+oWN+z9lB8ij4XZHru3fVcfSu12iTdhlx5TOvx1dQcdJ4eiU7/i8TAnw/iSkquqjn67DF9k/RkUd1m7EDXGTsk673ppGeW/OebPViwJw6TN1ieY0ENkmkEUPKAX9H2Cn5ZW88k43JKDv5QkAPbbP/CsYCp1PtLSsvFznM3XBZgcGanFtO3ZDT/gp3v1zCw9tfxJLy1zPwho9RoAUff55fbL+BccqbkSBqjYzuhfaLmT8MtAlkWPiN7Pj+znK0uv8IA83fH4Mnv7M8Db8rer81ZX7cq7SaZsrk4m40Zw/R552SunWb1SlNyX/nh2pO4fDvH4uf1U1Sc4g5B5QWDrW5kyb8JkkOebHHSZNZpe1ir97pZ5tU+Qejq7sXrmbLDfOSS5b+38ji+kgiSGKURMKj8agyx+cjGnjGWyH3vhic0pRdUNYJ8jnw6ckPYMw2C045cuLacvo62U7YZvSb1dRpOnKZme1DJ99By4ha0m7rN6D27k0dm/4OXF+zHxet3g+BFCnptGw6rNWX4GbsyB63ht3MjIw9rjyRixaEr+klwDMtZrBXw8jS+DH65XdnEGTP+Pic5KZBcLyp7OVJXPlpnfo76xiSPrCOnQod6tmrsP/72s/J5QKV2p4EG0fEpGLzoECJm7bL9gA4wOoebvNnTSelGDVtLbMntZvogNK+wGH8cTVSvt56Q+VtCsYp3P+6W20yjAfovPIBec/7BltPJ6Df/X8n1ioq12HPxJnIK5CfzdCa52exNxZikq9p3ZyKcTAuTkJa2lTKpCTLzSsp4PtmxNrSOYf5gZz3wkvo1a4VQbbSRPYwnm1G/vnWdsQODFx/CdhdNSjRnm32zvzs6esPe+yXT82fsTfOUKVITXsoNWbaVtdF8pTpBlh0HM+zhG3n2OsZaaMM6i6ViSy0rLNbq2+NSPxvTy6A7BOsUDde35fqtYFWphwzuEHi2lVs8ELgjJbsA7yy/OwpC9sGsxAsv/LAfvx28jMGLD8qmi8otKMbUDWcw4+9zuO2iEWeuwGCrm0lz8Gao73fSPRVt4aqKv+pOI/rRL+Unu7B0s7X9jPkNuMed+v7nsUTVcrfo/KNiLqlec2SGF8j+Q167qdusr+REqTnSwV7DIUaO/sQKTBqAUhfYI5fT8O2dwJKavblsKbvSXMtWh2QoP6RNDCdRe32JsplA5RjdpLlJ26FIa/ibK/nb9KGLYc9WW0ndiAiUNPDVOo/m5FvP/SiX/mD5AfNJskwfhqyyMa+ibgIyQ/a8VanJDaxJzy1Es0//llz24+4YvLfymOzwSakJFh1l/byiMU4jYLLBf7+JwjvLj+JQfIrVY3mYHMvS72v72RtGZZu19QJGrTyGIUsP2ZciznTYoszfSrZVSun3VZqd5YX+f3fpck5K1TWdWVvP47WfDsqmvnAGVU4/bnIet4WPV9m5rZH8jmz5zJ38/QjhaBoBeYYz2pem73eWXo8qw/O9vQ+KTM+fSkdH/KDS/Y71dDWldwK2K9hq8AbeWBJt1svfXuk5hZiz9Txibmbh5NV0vPTjvzguc82SK/fZaxlINbnnL9aWPJB46ItdsrmWtSa5lN3tIaQcS6U0fQtKgqZSn6uaH8ULP9x9eLrrvPMelrrTt+dIByHdqJKYm9myv/kig6FOh+JT8OIP8vWmPCk7rZJ7hDtUOsMyxN7MwmaVcp5YM3n9GVyXmTRFx9YneLqnKyNXHMOei7cM9uP4J10aM9k6Ix+ekouYI4EiJflW1Sa3a92QTmcOAbFFRp7rewgZvo+NBhOE2DsRRUp2AXp/tRu7DCa3sTZEt8jJub0sMXziqrQK21La7PwidJkeKTnczx5FCgppzxBRHaleyM/OlR8ONmjRIf3f9k48BwDT5XrbWXi7C/fEIq9QuqfNZ5vOYd3RRKdMqCLH2rvXaIArKdZ7OuuGoQshcC45A7kF5sFzW9IIyNEFBm210cJEQtavUXYd0iwFg9y+Svs+M/KcdI88Sw9ff91fEogtzd+mJYfiUvD9zktW2zBybYXo+BQcu5KG00nppdIOskWxVvmkTpbaQqXzuxLQagX2xdwyeMUyWx+O2RoMMw962M+0zRdpYUSCEo7mIXX02JFnrxv1OrTGuGdrSdkLi7U29egyTcPiyrylayQCle6eR1XNNDY6SWm5eGnBfnyz4xL6fLUH/X7Yh/2xKXh+/j7J9eXaSX2+3mP22o2MPNzMzEdiWi5yCosxbeMZiS2NuaKDlOE7EkLIjjo1ZMt9tzPfktLRDwfjjB+Cz1XxQY3ROV/mvWbkFWL7mev6iadNXc/IQ99vo7DKhly41pj+VuXv5802VLZ/g+vRm78ewYE4+XpTnjDY6mbUPMH8fvgq5u2y4+RgUIZHZv9Tqr0xrD1VsdTYkloy7JdoyQvR3yeTjfI42kNJEERNSoaTKOHKQJeO2g/DrX0V9gSFTG9qhCgZlr14X7zN+7LKoPxSn02SypPR6ZhexJ/63vacS9/tuGTUQxYAfjso38sLAL7ZoSwvqqOMh0Sav3YrKx+/HbS/oSJV/65n5ONGZr5qyeBdMQGvbkIDpUpr6FaORBDSVLbEkGcBYdTQ//NYIl6Y/y9uZDpWr6ydd6Iu3lI0OZeu1+quCzfR+6s9eOp78xEqZsFWqz2OrB5WsZOJJemJLt/OMcuXbTVnq5WCFhRpJa/R0vksXXvtKirW2jyxVmlS+umsPnwVX2w5j3VHEy3vT2aHz8//F09/vxf//SZKcdqV0qLVCvOAoYN1xVk3/nmFWszfHYNXFhwwOJblg42zkBN81/kbmLn5nGpBSTXed1Jarn40xRsOjJ6Jjk/BfRM2Y4nC9teW08mqpzN6Y0k0vttxSfa8b/jdmU5kpVv09Pd70X7adv0DNmvMzp93ciMu3KPuSD051upGpoIOBIfiU6wGqHMLiktyYlu4B7Tn96hmGhudrjN26PNAFxRr9Q+AC4sFPlp30iyVlC3XYsOJfwuLtJKjXRfvi0eCQV5nW+dtFkJg0d44hI/fiGkbrAdzrRm6NFr2IfgfBtcYS1/FF1uM03sUC4GVhy7jwvVM7L10C2kSIyalPlYlX7e9vZvjFObSVsLwnC/XRhr480EMWRqNL7dfkKwX0zedxcnEdIwziWWoPZF3UlquxdFNtpC6vjl7Yj13wGCrm1HzCdXo1ccxc/M5o3ygSvIJuvaGxvJVydZG5LX0PMRKNLhGrz4uO5mTUs4aumH0CRje1ErefEqz9DsqUvCo15GhQfZu+sexJPtnzLXym7WnTFI3NVJDRU8lKp98Q05+UbHs00u1Wfqs7BnOkSczhF2KLg9ewm3zGXXtZem37kjPy7v7d3gXDjF8B2evZeDnvXEuK4spR/KuyrH0+zQdSi9FySVi5IpjOBifgv85YbZ0Q3O2GQei5IqvO9+uO1JyY3LhuvmNuNk5zM4HTKYpWBQRJXnFH/xiJ7rO2GH0He06fxNxt7Lx637piVQstaNvZeWjxcTNGP6L+WSbUr8pV3eiNB3yqcb5RcrEP08p7j2nlNQkUboguiO+3eEeE7vq2BJcsfT9lUY7eN3RRHy+2TjIsEFm9mepHsSmrwxadAjzdsUYBThsZZqz1ZF7kqz8InSdsQMd/7fd7n3ojFp5DAVFWkz8y/rEZbsv3MTwXw4jYtYu1Xv9yZ3nAOPvo8fnOzFw0UH9v3X3C6fvtIH+Oq5sNIrp165BSSeYaU64dkkFtKx9fEquJ/3m/2t2LTT18bqTeGvZEby9/Kjy3nQKlHZv6OUHLqPL9B24YTBK096rhKXgnmHbxfQ3/vy8fbIjDoQQGPDzQUxeXxJkXRhl3J605wGFpfzLo1YeMzi2/D4S04xjE78fvooPfj+Jx77cjVcXHpCdzNnUXAcnGrc1v64qZD6Xo3c6PczbFYNWk7Ya9YYFpDsfDFlyCI0//lu/rqWRJ0KYn99NRzEtO5CArjN26H8vOrsv3DTqzWz62UjNMaCB5Waru42SURODrW5Grd/ajnN3h+ss2ZeAMauP4+TVdIxeZT05+KaTybI5AJ3tr2OWG4mWGtJywU9n9eR05kX8RkYeTiWmW23wyzUku0zfIbuNO/Rs1TH9nf1xLAlzd13Cv7HW83ldMpgooDTySin91Oxp3O++eAtNPvlbUf10N7a8X2s9Xu2htIe5rphKimuWJ01uvVK4ITf87fZfeMDCmuXDNonc2zpKJomSuw5IvWwaPHMV3buy9Gsyfe9HrqRa3GfUpVuSr9ubF37nuZIh8Om5hSgsulvSD9eeRMSsXWY3bTqWzg9/HE1EYbHAVonvPCHF/IGM1L5c0fNbx1n1f8m/CVgVrU5uQZ1zyebBVmujNJS8Oz/vu7cRageI7ZEtkeNa7jdSLAS2nk7G1VR1JzV0xDWZUSxS1zm5dvro1cex4UQShBCYst7+nmtCOJZGwDCFiqNBTyWb387KhxDCKN+z2g9Lb9swAa3hudZsBnkFBTufnGkW3HVmjlS18rxK+Unm+qCz9s4Dgt0XbsqeV+35Km0N4Agh8OHaE/haYsJlW3xpsL2939lz8+SHV9/Kuvs7NH2L0QmpGL9Wugd8fpHWKKWeqcsm190YiQnZAPtSotnyTRy5bNy+0QVjjQJ8Eh+row//LJUxPbcQyXfOz2o+xFG6p9lbTR7eG7x/3YguXeD7lQUH8G3kRbSbts0ogL7m8FW8smA/0nIK8MqCA3ji2yhk5BViX8wtZOcXYeZm4/Reup7rUm0Fwzp9wmRydsl8urCQ8m/reXT433azoHt5wWCrm1Gj8d7n6z14ffHd4Tq/7E/AmsNX0fe7KMVPrVpM2KJ6V3QlvrFyorT8lEb6dWe9D2c+hHngs0g88W0ULpoMNfopKs4oP4vc0/FkC7lv5W7Cnc3w49KdhztMNe/xYNrTQ85Ggx4g1iZTUdIjzhqlF1fD34XSNpauV6vUrLJkmWkA0jDPktHnrwu2KjjHmtYruS2SnZTewZDhW7DlRq806H63jrQ7TYOjFocUKahPcvVU6nztyifphuXUBVItnWM8TU5igw1y5zqdyef+6Jf/KN70hoWZ1U2HDhqaKjG80X0eE5ZwZq93S9dwe9jTDDon0RvWlK+Xp/7vHAV5+5xt1Mqjitt8v/ybgGESvaoB149oMCX1cH/TSflUNe8sP4ojl9OM0pisir6CLaeVp7dR8yNQ8nkWawX+tNLhQmfd0atGN/ubTyWj/bTt+PTPU5K5UkuDpUPl5BdhlsH5zlqxvom8iMe/2m3WXlejLSsnX2YouBpsGVGh5ldmaxqBM9cy8NvBKw6nRzHsJSz1lS1QENhWWnSpdoPcwzpr9cH09/W2xNwDSvKzShFCICW7AOuPJyHDSspAqabZjnPX0ezTzXYdW8eRTlJ7Lt5C5+mRSMkuUPX+X+k5ylLv0fbTzCfGnr3tAtJyCjHj77u9oMesPo59Mbfx5bYL+Df2Nk4nZaDLZ5F4ZcEBfL/TtkB1jEF84hcLPf4NybVxv91xCSnZBfjKSg/4sorBVjejG2brCKnhYtYsMhmaWqQVWKaw8pQmSxdOuZOoK5PqO+qIwURDV1NzMHXDGaP8LO/b0RNy6b/Wv1fnTJBl/G+tViBTIr+iPSz1hpu64QzOXVOWEF0NSW7+ZE7tew9X35QeMAiuWvrsVx++guT0PKP15aw9ouyG78Uf9ytar7wavPgQ/jyW6NBDwnwb0mcoGbYtm09PoowXb2Tiy20XrDb81XLBIKWPYb0Zvfr4nZEM8lzZgzM5Pc/o87PlOzt/PRPr7zy82B9722jouuF+pCYFM+UOE2QZHdupO1d37/b0Ol2ioK3gDr1ZDeUVavHB78ry6lrqRWPYk9SVdU+n0Er6J9PJXICSYfyGxq05geG/HDZ7cGOJIz9Dw9m8P1AwR8Kv+xMwcsUxs9elUiy9t/I4pm44o89d+fmdHlm6yeh0SvP0YOk6OHnDGXxnEMywds2UG3bvzJ6tvx28LDuJn6MPJoUALlyXbocr3bc99yVK9n0uORPpuSVtALnco7by8777EMq02RJ3Kxv/26ReGghb5gew9hGa9qCXmrTakflKXlmwH+/+dtTowYMUqe967Grjc4iuLtiSzufJ7yynI1Cyp3PXMuy+R35n+RGEj99o9JrSXZlWfcN/W/rdSu0/Lfdumzf7TttrpY2pFS3FY+Q+RyU/HSEERvx6GG8tO+ySyd+cgcFWN+OMPDyGcmWeSJnm4wCASQ4MP3IWS4FTuWVKcpSWBYaJ6J19AiqN4bXOSFwv5aeoOKNGrr2UlrbH5zsdPlZZYuvXOG6N+qkSdEOOLJ0fpm08i/98s0fR01t3uLnW0WhKhv58YGFiFFcaueKY3TfktqarURIjkPsNyE1m9nXkRUz+q3SudYb5nU3Pf098GyV7khFCWBz652x/HlOWX1COrgfrSxYeTrxnkNtNju66V6wViLp4q9SC5HKc2WtO7T07a0LP7IJitwu4yjmXnKF4Qhi1J1dylKWf2rEraUaBTR25a52lYOt5w4CY0P/PYUompdl5Xjr3Y5NP/pYNjMfeyjabhNDw+q3kvFIaNprk4rW3OjqzaZJbWIwBPx+UXNZP4vdlq/n/mE/YLITAkyYTQsp9NKYPD5TQipJgvemwdEMnrqbj/slbJedkUIPuOzuckIJJf51WvUOG0l6FgPVr1reR5ukT1h65iie/i3K43EJAP5nu71bOB1LFVONBw2kLHdqmbzqLG5nWJzQVsL/+SuXktrezgtLPQ2rvajQHLN1vSZVNA2VtpltZBfj7VDI2nUxG5NkbOJ2UXubzuXq5ugBUuuzN1+YuLFU4uYqfkatO70lXMDxfpeXeHZaSW1iMAJ+yVX1NLyilObxLDVLFtTZc0ZYhe6Vl53npngv2srWhYCkXob3Dk56duw/xM/5r9rrpk9qU7AIE+dpWb+buisEXz7e2q1zq0GDMavfO5WtvTbZ11IFjOVvlj3XUwo2Ys2RJ5DyTq0vnZXoElaaNFoYsW6PknmCzgnOl7oHjz1FxqvYMspczL2Fq79uZAdHb2QWoVsHPaft31LxdMagY4INn58rnQLTELZoqMmW4dCMTT3+/V3KZVG9XAPA0qJAH41Iw3WCYqW5SlpJDCqfdM2g05u35XQrbJobncl06lTqh/pLryk045gy2/E7sbf86s2erjlS79XCC49fItUcS0bJmCJ5tVwsZuUWoFOiNYq0wm2RW7lot1SnImmIhMH7tCUWjlT5adxK/j+hi9vr55ExUr+CHkABvxceVeg/PzSsJWKsxMaEt5bickoO6oQHQaDT6HrxyjhjUfx3dCMppG8/gs2daYZ2dE/AZtm+yrYxkkaofZj077SqFtPPJmYpzFr+/6hj+GRuh2rEN3+rm08k4EHsbnRqEma2ngQYp2QXw9NAgxN/brD1sSycsqc/X1jOSpfb7jnPSD84UzZlhUJIhS0tSYr4d0RD3166IbWeuY+rTLY17jpcBZStaQ/c8eybIGrzYebntJq+3PjOqIwx7o3y/8+5T4cX74vHWw42cemx7yQUMDsQaN/zXHy+9RrCz3D95q8Xlht+ZUnI3SGpZr3AWXKXUuhFNzS6QzDuk1BWJCXWkhuIpCQ6bNuJcm0vXHe70LXvtp9KZuEvJfabc79FSXDfWBb3Y2k41/63LlV1q1tnSdtxKXmxL7L0pupKSY9QbVpezeI2CXnKlwZkT5Nnz4MlSQNWZI5WUPARxpUPxqXYHWgFXn/9LrDh0Gd0aVTZ7vdec3Tbvy7Bnq1SPWJ1D8ak4FB8tu9wRWiHw4Vpl6R5MSQ6LNQgKl0ZAUoot80NYajfN2yXfbjTtxesMw2XyGKthyoYzmHKnd3mgjyd2jzMPWql5VtVqheK0UKYKirQ4FJ+CVxceQICPJ85M6a14W8P34GHSk1yNwLUSP0XFYcOJJBy9nIZxvZvitc710H2m+ci7iFm78Mdb3awGkzedTLaYJ9oaW+4VpFY17ZD/w+5YDOoWbraeLfUwp6AIAT5eyMpX/lDpekY+Np1U7/7V9HN58cf9kp1H8oqK0e5OuzH2s/+YtclkU2hJvS7xmq2dH6Lt+B0rCQhLpYUwvJcOrxyItyPcM/4hh2kEqEyxlJdESfd/tS3aG+/U/ctNFmVp0hFX6/3VHsnX3/3tqP5vAbh9bz1z5hcJZwRBdLNJlhVqNYy3nbnu0NCWVxYqy5+qpMFnerMmNcN1aXGLXlVWHJXoDaGEtfe2aG8cHpm1S/8AQskttFwPybLwORrmXstWKZ+1O9BoNNhuIa+2nB6f75QcPuzMIKdNnFiMhVZm75by/HzHh/raQwiB1OwCTPzTuQ+f72XT/z5XkmpEBc6cZEmpRXvjsTLathyBOlIdK2zJI+0sFid3NGEYdMgpKMKY1cfx++Gr6PP1HrMZwV3JmWfa7IJiyTafqhNk2dio1PU+BYBec/7Bq3cmYNW19e2aKBeumYxz6oYz+rbZ55vPyw6hj7uVjYGLDtqVpsEWtnwCUpMeSz3U6zJ9h1k6qnUGwXXDwJ1UKolv70zKbevXczVVvVQQUu0ZqYet1wwm5N16JhmXTTqYyPeWN39dql5Y622/L8Y4ldVNmZiLpWdd1j5nAet1rDQmJlYbg61EZVBZTxqtu8CVJWX8I3d7pk//bXUlJVfZEBU7vsfSHPblCmo+pbeFta9i8voziL2Vre99paRBLNVrRAg3CtApNGX9Gf1NRFk/92hwdziYGtzl83CTYricVpSM8rF3eCmVLkevta4mVe+kJtEqbfG3lY+QMAx8zd8VgzWHr2L06uN2TXBclkkFeJYdUG9y5hHLjti9rWkwKyktV3FQzvAeTaPR4Ln59vesLw3HrqThDYlRoGpe4y4ZzF5vjVT6ErkRFFNM8nDrRsCY+ksi97yuvWhrMHzuLvXuYaUe0kjN/WEY3JRKj2JLahKpda19Bq8sUDaCzdJIFyXtcGv55ctaWx5gGgEi96OgHbzk3wQM6dHA+WVxkhSZi6E7K3un99KhVuBDjfu/MSpNvlW2b0Vt95YDNySOsKVxmJlXiIo25EwzdPFGFr7abj7xgztbGX0FK6OvYMITzfGPzAzRZYaKFepGRp7NvZWcJV+lmavLOq0Q+olPiJzN2nXDVZ0RbBnhsS/mtv7vJDfuqeXstpDUuVzNEYNqBq+7ztiBFzrUVrSu4U9QA/tH/6jJWrU4IJHCzJ3meZGbJO/EVfnOEIaTenlJ3GQE+pTk/rS1SWFplK0a5HqN6hRJBGjlvl+p16V6Mas1abWl9pm1z3nN4atoV7eSxXWEKOn5Gx3v3JR7amKwlcjdKDzf3Wuz3ruau/SmcjdqTXT2qw0zqspRku9WrsFmKFImuTupy5ZZ7m9k5iM00MeJpXFPpr02yiI1b9hTcwrdYsgwAOy5VMaD4CrRCqHovEruQepGvSyx1uRwk2cxVr2z/AievL+m2+SgliIAvO3Eh7HOnvdCbZYmeDVk2PvOXVJaj1p51PpK5YBhGg7DXq5enuaDuXUTBpe1CZulen/KvQepdQ0f9ui3d3Kz6q/jSfhLwZwhH62znMtbABi14hg2HVE2oZk7YLCVyM24YrIWInupNXxPahZUKt8+/eOUTeuXlZtoMqb2hDW2TIDhTO7Sw9bV5u6KQWZe+ckxXN45Ozejs1kLjJSVwMmGE9ckhwO7kyspOdjoxDRDW07bnsu7LHDHn+B1N57rw5k2n0pGzM0syZ6tABA+fmMpl8hxRRKRUbnmiNKRUQVu0q6yRgj5uRncFYOtREQKlMU8MaWBnwuVBg3Kzk00GVMz1Ho4IdUlk2FKmVAOJoTKVWGSx+UHLqtQEiJlBv580OJyPgNRD/Mw24dtFffx5q+HAQA9Gld2cUnUIzU64V75zblDfm5bcYIsIiIF3ll+bwzBsRV7d1Gp4U/tnmdtiBnZ5v4pW11dBCKbHIo3nwTRUFmfQJbKPsOfIH+O7qEszhUiR6oXqq2TfJVVvx9x37QrchhsJSIiuzHY6lz8dO/6N9Y8zxS5v4s2zEJMpass9hIhsmTerhhXF4HucYbN4tvlKMhXlnmqMQuvmyiSCLZauhfjAyjXYrCViIjstv0sJ5NyJjaSSpxMTOeQRiIqkyY80dzVRbhnSE0IQ1SarqbmuLoIZKI8BVul5riwFGx96vu9TiwNWcNgKxEREbm1kSuOuboIRER28fXm7RbRveJccqari0AmyvszGEsPmU5cTceh+JRSLA0Z4tWfiIiIiIjICTSqThNHRES2OH4lzdVFcKpuM3dYXN5v/r+lVBIyxWArERGRm8orZE5FIqKyTMNYKxEROQkzjrkvBluJiIjcVG5hsauLQEREDihH6QKJiIhIIQZbTWRmZmLSpElo1aoVgoKCEBISgo4dO2L27NkoKOCMgkRERAQE+3m5ugh0D/tv6xquLoJTff5ca1cXQTUadm0lIiK65zDYaiAhIQGtW7fG5MmTcerUKQghkJ+fj+joaIwZMwadO3dGamqqq4tJRERELtagcqCri0D3sBY1K7i6CE5Vnh5meDDYSkREdM9hsPWOoqIi9O3bF/Hx8ahRowa2bduG7Oxs5OTkYMWKFQgODsbRo0fRv39/VxeViIiIXMyTY4OpFLSsJR1ULe+TLpWn3qAaAD6evOWisq1KsK+ri0BEVKbwyn/HkiVLcPLkSQDA77//jl69egEAPDw88OKLL+KHH34AAGzatAmRkZEuKycRqadzg1CM79PM1cUgBSqUo15O94JjEx51dRGcrrz1VhvQpZ6ri6CaPi2ru7oIqvD00KBB5SBXF8Ml/H08XV0E1Xh4AN6e5et8YWrx4I7Y9H898N9W5Tu9xb3Ms5xd84iInI3B1juWLFkCAIiIiECXLl3Mlr/00kuoX78+AGDp0qWlWjZbjX28qauLQFQmrBjWBQ2r3Js3smVNeerlZEl4WACeblPT1cVwSMfwSqgY4OPqYjjdK53quroIdul7v/Tv6/Vu9Uu5JM7zv2da4YfX2ru6GA7z1GjgJdODuryfEmuG+Lm6CKrRQFPuHs6YerhpVTSvWQHfvdLW1UVxWP/OZfPcbovt7z9o8zYczEFEZBsGWwHk5ORg7969AIA+ffpIrqPRaNC7d28AwNatW0utbPaoUY4aqGqoVdHf1UVwOXe4aQnx93Z1ESS1rVvR1UUgGV+92Aa1K/ljylMtyn1gAQCiP+mFre89hDZ1Krq6KHbbMy4Cv7zRCQDwXLvaLi6NcsuHdrJ5m2fa1nJCSZyvWKuVfD2gHPUk9Pf2xGPNq7m6GA4rKNaW+3QV/2lVHeem9jZ7vXKQb5n+DlvVCtH/rdEA5Tzrg54tD0attb9c0Xb19fLAtKdblfpxne2RZlWN2lGNqgZLrhcWKP+g1Np36+oHXPakOfjrnW5OKIlyg7uFu/T4RORcDLYCOHv2LLR3bj5atmwpu55uWXJyMlJSUuw61ruPNLJrO1vkFUrfSDmTs3JRqREo3TyqB34a2EGF0tz1ToT09+jn7SE7hKpeWIBDx/xtaGe7t21Y1f7em//Xs7Hd2xo6PvEx2YCGadAi2Lf0hoxXDvLFwY974vTkx0vtmM2qSzdynemljnWsruNOQ/U9NMDTbWsh6oNHMKBLuFPuUxtVDULd0ACMeLihw/ua9rT8tUOJTvVDUTnIFz5eHvDxkg56GU4Y071RZYeO5yx1QgPg511S/keaVXVxaZTr2rAyLkzrgwlPNFe8jUajQbiD5/X29SrZva2lm2JL+nWwfi6Q8k5EI6MAkjWz+91v13Hs9dWLbTCpb3NMfaoF/H08rQYGdox+CA+Eh+KVTnXRqGoQ/tu6htl1tqrCm/eRPRtj+RDlAft2dSviz7eV3eS/92gThEp817aeEzUaIG76f2zcyvkm9m2hP2cYCvD1xI8DOuD4xMfg5229jTn8wQY2HVeqh/f/mbTRB9qZWiMs0Ad/GHy/HhrzDLsNq7h+gr3PnmmF9vUq4Y3urunV/l6vJrLL/Lw9sPj1B1Q7VrdGYRaX//Bae9xXowJ+G1ZyDtCdC0Y/Kl9GZ6oa7IthNvym64ZavhZNeaoFhLC+nz0fROCXN6Q/9ydam9/ffP1SG0zs2xyrhnfB4y1cl77l1zc6Ycfoh2zuQNGkWjAimlZxTqGs6NmsKib2bYFRvczvs55t55qHuWo/2yuLD8w+VJhezhX3cra6FzqquDv3ubN2oaSkJP3ftWrJn9wMlyUlJSE0NNRsnfz8fOTn5+v/nZ6eDgDQ5ucgLNAHQzvXgDY/B+tPJOHFDnUwc/N52eM1qBKI2JvZ+n+/2qku2tWthNGrj+tfu792CI5fTTfarnaQgDY/R3a/ADC+T1O81LEu2kzZZnG9BlUC8UB4KFYcuoLalfyRlJaLnwd1xO6Lt7A6+jIy84oBAJOfa4UPfj+p3+7j/zbDjnM3Mb5PM4z49TCS0vKwclhnvPjjfgAljYj+netizraLFo///kONMX93LM4nZxq9/soDddCtcWV8vPYk0nKLzLaLaFYFO8/dxOLBHSEKctG+hh8mPh6ORlWDMOHP04i587l+2Kcppv8t/x0AJQ0+wwB2n5bVMaBDVURfSkSwjxeOX03Du480xhN3Gu6eHhq0reGLBpUDcfZahv47/vKZNhi18hjib5V8NyMeaoh5/8QAAHrdVxVTnm6JYF8vjFxxFDvO3TQqw8znWqFFFW/Z73XPuAjsPHcDAgI5BcVoWSsElYN88friQ/D38cTEx+vjkVOXLb5PoOQJ66K98UavDexQFdX8ivHh2lMAgA96N8UXW85Da9Jo2/h/3fHfb6Jk952RkYGWVXywbkgbPPX9PqNlEx6vjzERddBl+g4AwIaREXj3t6PIzi/CxRtZVsttqFKAN1JzCs1e/3tkD2w6eQ3f7rgEoOQ3lJGRAQDwA1CcX4DKPsW4kZlvtN0zbWtiXO9m8Pf2tFhfnm9fC+892hT5hcU4l5yJt5YdMVunor8X2tSthIZVAnAm4brRsrqh/rickmvTe139ZmcMWRKN9NwirBjWCS1rVYQQAm/+ehh7L902WndA+6pYHlXyW+zTsjpC/L2x4tCVkmVd6mFc72bIKyxG/K1shAb64JHZ/8ged2iP+rivRgW8v+ruuej4xMdwOikdRxJS8XCzqhi6JBrFWi0WDX4A0zacxYCu9fBN5EWcvXa3Lg/uHo5FUfEASoJOhxNS0bVRGLo1DMPDTavqvx8AaBjiiYOpxuc6AFj3Vld8uO4kvDQanErKwJo3u2DW1vPYH2v5gdgLHWrj0yeaQ6PRIL+oGKv3nceNzAL9cg9NSbDo0o1sDOoWjrzCYjw/fx+uZ9z9fQzoUg9L/01AnVB/PNm8Ej5aafm8CwC7x0Ug4XYW3l1+FJ/3ux9aAVy6nokBXcL177dJqIdZXd/23oPYffEmpm44ixohfpj7wn3ov/AAjl1JQ4CPB3IK7p6jVgzrhJd+PKD/95Du9bEwKs5q2YCS813nBmG4kpKjP0/WCfVH7YoBCPLzxLYzN2S3rV85wOg7614vAPNeuA+FxVo0r1HB6Dfl7emBwmLzB4NDutfH0ctpOHw5VbJs9j5MfK1LPfzyb4Lksvtrh+jL/XzrMDzWpDNWR19Bh3qV8MPuWOy+cAsAUL2CL955pDGOXUnF021rISMjA4tfbYkHv9gFAGheowIKiovRo0kVvNi+Dnp/vcdquX5+pTuWHUjADInrUKtaIQgPC8D6E9ckt13zfw9g48lr+FyiHfHRf5rh7LUMrDuaZLasXXVffN+vGUb8WnKO8vHyQPUKvvDR5uPIB93w+5GrmLrhrNE2TaoFYViXGuhSxx8vLzhgtk9TLWtWwKONK+h/x7Uq+SExNc9sPW9PD/zxdlcE+HjiamoOXvvpkNV9v969Po5dTsWRy2lGrz/S8O6Nj+77rBsExN8uKcPS1ztiwM9391/ZV4uFr7S486+7QTU/kYecAi3mv9YebepUxKaT1zBl/RkAQN/7a+BA7G39ueKUyUO6Ex/1wNQNp7H2aBICvD2QLtE+AYCIBrVRP+RuPe91X1VsP2tetzrUq4Qgj0Ls/L8H0GrS3VFVQX6e6N0kBN9tKcBDTariv61rYPgvh2U/M6Dk3JOZmYltb3fE6ugrmL87Vr/sje71USXYBzP+Po//69kI9SsH4r2Vxy3szZi1utm5QSh6Na+GaQa/q94tquGlB+rCHwXIyCjAnGeaYFX0Fey7c+3Kz8lGPkqCygfHlgQur6bk4NsdF7HxZLLR/j994j6EBnjqP891b3XF3F2XjM5Xm0f2QIFWi80nk/Fw0ypoXjMEfx68ZLSf9jX9jM6991f3Rec6/voy6VSr4IvI0Q9j8KKDOBSfiu6Nw1AvLBDPtauFJtVKJjTLzsrU7ysABXi4fhD+OFZSH1/tVBcf/uc+DF0ajX9jTK7VXerhVGI6+nWorW93OUuvRsF44r6WOJKQggWRd9/32re6wNNDgw9/P4Uz1+6e08c81gT1wgLh6anBd5GX8GSbGpjx93k8376W0bn/l9da4kBcCmqE+OHstUw0rBKIiX+d0S/38/bA8IcaoG4QZNu1r3VqgGp+Wqv3M3KC/DyRdece5f1HG2Nwt/o4nZSOWhUD0OPznWbrd6kTgC6DSx4QZWRkoEUVbxwa1xW+Xp64ePUG/rlwE3VCA3Diqnk7BAD8fTyQe+c6XDXYF/XCAvDOI43w9vIj+nL071wPv+6Xvg6Z2vZhd2g0GgzrXANv/nIYB+Pvtmm6NAyDRgO89XBDvLPsCP73bElP3LeXHZXc108DO6CCZxHCK2gQezMb3l4eyMjIwKFxXfHh7yew9U49WfJ6RxTl5SA82LwNAgAvt6mMedvu/ibDAr0R0cD4vGvv9wUA4ZUD9PdIUk5NfhzpOYXoNnOH/rWWtULQrHoQWlfzgSjIxZynm+CnqDgslbjej+rVGPXCAvHeymP61wpys/FQeBAiT0h/L10ahqGCrxdCArwxsmcTFBQXI2LWP0bLR/VsrL+/NdW1URi+eaktOkzbbrYsokEQMjIy0KdJCOZsvPu+v3qpDSKaVkWzUE8E+Hjho3WWzwNhgT6oWsHXqH1tzejHGmP2VvN78Hd6NkKDKkEYteKY1X2sGNYJdcMCkV9o/JnodGsUhplPNsbmo/Ltz6bVg7F6eBdM3XgGq6Ov6l9vXDXIpvu/5jUqYPHrHfHA/8zn1YloVgUnrqThdrb5vSFQUl8N7/36Nq+E//1h/Xe8bGA37Dp3A++tMr5WtqwVglnPt8Z7q46ZfSfvP9pYH/v4oHdTfLPjIj7scx++23ERNzIL8Fy7WmhaPRh1QgOw+8JN1Krojy+3X0SxyU33wK710Kx6sNE1YvqzLY3+Xa2CL5YMfgALo2Kx5nCi0fbBfp76GA4APNy0Ct6OaIhBiw4hO78Y1hjW1b7310DczWycSsqQXb917RD9uVOuPWjN8+1rYXC3+pi64Qz2nS1plwslT5BcTCPKQimdbPny5Xj11VcBABcvXkSjRtK9Frdt24bHHnsMALBv3z7J3K6TJk3C5MmTnVdYIiIiIiIiIiKie1BMTAwaNLBtVEtpY89WlX344Yd4//339f/WarVISUlBWFjYPTPBS3mTkZGBOnXq4MqVK6hQoYKri0N0T2I9JHIt1kEi12IdJHI91kMi10pPT0fdunUlR5m7GwZbAQQH3x0CkZMj323ccJnhNoZ8fX3h62uc46tixYqOFZDcQoUKFXhRJXIx1kMi12IdJHIt1kEi12M9JHItDw/3n37K/UtYCmrWvJskPzExUXY9w2WG2xAREREREREREREx2Argvvvu00fGT52ST0KtW1a9evUy0W2ZiIiIiIiIiIiISg+DrQACAgLQrVvJTKebN2+WXEcIgS1btgCAfpIsujf4+vpi4sSJZukhiKj0sB4SuRbrIJFrsQ4SuR7rIZFrlaU6qBFCCFcXwh389NNPGDJkCDQaDf7991906tTJaPmqVavw4osvAgC2b9+Onj17uqKYRERERERERERE5KbYs/WOgQMHolWrVhBC4LnnnkNkZCQAQKvVYvXq1Rg6dCgAoE+fPgy0EhERERERERERkRn2bDUQHx+PiIgIxMfHAyhJL6DVapGXlwcAaNu2LSIjI1GpUiUXlpKIiIiIiIiIiIjcEYOtJjIzMzFr1iysXbsWcXFx8PDwQJMmTfDyyy/j3XffhY+Pj6uLSERERERERERERG6IwVYiIiIiIiIiIiIiFTBnK5U5OTk5+PvvvzFt2jQ8++yzqFevHjQaDTQaDSZNmqR4P8nJyfj000/Rvn17hIaGwt/fH/Xq1UPv3r0xY8YMFBYWym57/fp1jB49Gk2bNoW/vz9CQ0PRo0cPLFy4EEqeX8TExGD48OGoX78+/Pz8UKVKFTz++OP4/fffFZefyFXUqINr1qxB3759UbNmTfj4+CAwMBBNmzbF0KFDcezYMavbsw7Sve727dtYtGgR+vfvj+bNmyMwMBC+vr6oXbs2nn76aaxbt87qPjIzMzFp0iS0atUKQUFBCAkJQceOHTF79mwUFBRY3Z71kO5ljtRBNeov4Po6TORqatUlQzNmzNC3azUajdX1WQ/pXqZWHczIyMDMmTPRtWtXVKlSRb+PiIgITJo0CWlpabLbum0dFERlzM6dOwUAyf8mTpyoaB8rVqwQFSpU0G/n5+dn9G8AIjU1VXLb6OhoERYWpl8vKChIeHl56f/9+OOPi/z8fNljb9y4UQQEBOjXr1ChgvDw8ND/e/DgwUKr1drxyRCVDkfqYF5enujbt6/RNkFBQcLHx0f/bw8PDzFnzhzZfbAOEgmj37zuOhYYGGj0Wp8+fUR2drbk9vHx8SI8PFy/bkBAgPD19dX/u23btiIlJUX2+KyHdK9zpA46Wn+FcH0dJnIHatQlQ+fOnRN+fn5G21vCekj3OjXq4I4dO0S1atX06/v4+IiKFSsa7ePo0aOS27pzHWSwlcqcnTt3ikqVKomePXuKsWPHit9++01Ur15dcbB11apV+hu6YcOGidOnT+uXZWRkiN27d4v33ntPZGVlmW2blpamP1azZs3EoUOHhBBC5Ofni++++054e3sLAGLEiBGSx46NjdWffLp16ybOnz8vhBAiMzNTTJgwQV+pZ86caccnQ1Q6HKmDhr/zt956S1y9elUIIURxcbGIjo4W3bt3FwCERqMR0dHRZtuzDhKVACAeeOABMXfuXBETE6N/PS4uTrzxxhv633L//v3Nti0sLBStWrUSAESNGjXEtm3bhBAl9XDFihUiODhYABD/+c9/JI/NekjkWB10ZFshXF+HidyFo3XJUHFxsejatasAILp06WI12Mp6SOR4HYyKihL+/v4CgHj22WfFoUOH9A/bs7OzxcGDB8XHH38sYmNjzbZ19zrIYCuVOUVFRWav1atXT1GgJykpSVSqVEkAELNnz7b52J988okAIPz9/SUr/GeffSYACE9PT/3No6H+/fsLAKJ69eqSPWeHDRum7+Fj6QkMkSs5Ugd1Tx4feughyeVpaWkiKChIABDjx483W846SFRix44dFpcPHz5c38C9fPmy0bKFCxfql+3bt89s2+XLl+uXb9++3Ww56yGRY3XQkW2FcH0dJnIXjtYlQ1999ZUAIF599VUxceJEq8FW1kMix+pgdna2aNCggQAg3n33XZuP7e51kMFWKheUBnrGjx+v705uz/DEunXrCqBkeKOUzMxMfaBowoQJRsuysrL0T20mT54suX1cXJz+hPDzzz/bXD4iV1FaB3XDOkaPHi27Trt27QQA8c4775gtYx0kUubgwYP63/LatWuNlvXo0UMAEBEREZLbarVaUb9+fQFADBgwwGw56yGRdZbqoKPburIOE5UlSuuhbsRFWFiYuHHjhqJgK+shkXWW6uD8+fP1D99zc3Nt3re710FOkEX3lKVLlwIA+vfvryjhuaHz58/j8uXLAIA+ffpIrhMUFIQePXoAALZu3Wq0LCoqCrm5uRa3Dw8Px3333Se5PVF50KBBAwDA4cOHJZenp6fjwoULAIAOHToYLWMdJFLOz89P/3dxcbH+75ycHOzduxeAfD3QaDTo3bs3APN6wHpIpIxcHXR0W1fXYaKyRGk9HDp0KLKzszFnzhxUqVLF6n5ZD4mUsVQHdbGZfv36Ga2nRFmogwy20j0jLi4OSUlJAID27dvj5MmTeOWVV1CjRg39bHcvvviivtKaOnXqlP7vli1byh5Ht+zMmTMObX/69Gkr74io7BkxYgQAYNeuXXj77beRmJgIABBC4MiRI3jiiSeQlZWFLl26oH///kbbsg4SKbdr1y79361atdL/ffbsWWi1WgDK6kFycjJSUlL0r7MeEikjVwcd3dbVdZioLFFSDxcsWIDIyEj06tULAwYMULRf1kMiZeTqYH5+PqKjowGUxGYuX76MYcOGoU6dOvDx8UG1atXQt29fbNy4UXK/ZaEOMthK9wxdbzkA2Lt3Lzp06IDffvsN6enp8PPzQ2JiIlatWoUePXpg6tSpZtvrArUAUKtWLdnj6JZlZGQgKyvLbPtKlSrB39/f6vaGxyMqL95++22MGzcOHh4emDt3LmrXro3g4GD4+fmhffv2uHTpEsaPH4/IyEh4enoabcs6SKRMWloapk+fDgDo0aMHmjZtql9maz0y3Yb1kMg6S3XQ0W1dXYeJygol9TAxMRFjx46Fv78/fvjhB8X7Zj0kss5SHYyPj0dBQQEAIDY2Fi1btsSCBQtw48YNBAYG4saNG9iwYQOeeOIJDB06FEIIo32XhTrIYCvdM1JTU/V/f/rpp6hZsya2bduGrKwspKen4/Tp03j44YchhMCECROwdu1ao+0zMzP1fwcEBMgex3CZ4Ta6vy1ta7jccFui8sLDwwPTp0/Hzz//jKCgIABAVlaW/mKbl5eH9PR0ZGdnm23LOkhknVarxWuvvYZr167Bz88P3333ndFyteqRo9uzHlJ5Za0OOrqtq+swUVmgtB4OHz4c6enpmDRpkj7VlRKsh0SWWauDhrGZadOmwdvbG6tXr0ZWVhZSU1ORkJCAfv36AQAWLlyIL7/80mj7slAHGWyle4aumzlQMmT5999/R69eveDhUVINmjdvjvXr16N69eoAgMmTJ7uknETl2a1bt9CzZ08MGjQIXbp0QVRUFNLS0nDt2jWsXbsWVapUwbx589CpUyd9igEiUm7kyJHYsGEDAOD7779H69atXVwionuLI3WQ9ZdIHUrq0q+//oqNGzeiTZs2eP/990u7iETlmrU6aBib0Wq1+Omnn/D888/D29sbAFC3bl2sWLEC999/PwDgs88+Q1FRUSmVXh0MttI9Izg4WP93z5490a5dO7N1goKC8PbbbwMATpw4gevXr0tun5OTI3scw2WG2+j+trSt4XLDbYnKi4EDB2LXrl146KGHsGXLFnTr1g0hISGoXr06nnnmGURFRaFy5cqIjY3F+PHjjbZlHSSybMyYMfqeA19++SVef/11s3XUqkeObs96SOWRkjro6LaursNE7k5JXbp+/TpGjRoFT09PLFiwAF5eXjYdg/WQSJ6t7dHGjRvj6aefNlvHw8MDY8aMAQDcvn3baILlslAHGWyle4ZhLg7dLMdSmjdvrv87ISFB/3fNmjX1f1vqcadbVqFCBf0wacPtU1NT9TMxW9re8HhE5cHZs2exadMmAMDo0aOh0WjM1qlatap+coK1a9ca5edhHSSSN27cOMyePRsAMGvWLIwaNUpyPVvrkek2rIdE0pTWQUe3dXUdJnJnSuvS+PHjcfv2bQwbNgzNmjVDVlaW0X+69FYAJF9jPSSSprQOGsZmmjVrJrs/tWIzptuURh1ksJXuGc2bNzebcEeKYXDHMBhkOEud4ex1pnTLDE8M9mzfokULq2UlKksMZ3Fs2LCh7HqNGzcGUPIk8caNG/rXWQeJpI0dOxZffPEFAODzzz/H6NGjZde977779OlzlNSD6tWrIzQ0VP866yGROVvqoKPburoOE7krW+pSXFwcAGDevHkIDg42+083qQ8A/Wvjxo3Tv8Z6SGTOljoYGhpqcWIqHbnYTFmogwy20j3Dz88PDz74IICSHnZydAEhjUaD8PBw/etNmjRB3bp1AQCbN2+W3DY7Oxt79uwBADz22GNGy7p3766feVlu+4SEBH3ZTLcnKut0F0TA+MmkKcP0HYZPEFkHicyNGTMGs2bNAlDSsB07dqzF9QMCAtCtWzcA8vVACIEtW7YAMK8HrIdExmytg45u6+o6TOSOHKmH9mA9JDJmTx3U/a6VxGYAoH79+vq/y0QdFETlQL169QQAMXHiRIvrLV26VAAQGo1GHD582Gx5ZmamqF69ugAgOnfubLb8k08+EQBEQECAiIuLM1s+c+ZMAUB4enqK8+fPmy3v37+/ACBq1Kgh0tLSzJaPGDFCABDBwcEiJSXF4nshcidK6mB8fLwAIACIvn37Sq6TlZUlGjRoIACI1q1bmy1nHSS6a/To0fo6NWvWLMXbLVy4UH8t3L9/v9nylStX6ve7fft2s+Wsh0Ql7K2Djm7r6jpM5E4cqUtyJk6cqN+nHNZDohL21sHdu3frt1u3bp3Z8uLiYtG6dWsBQNSqVUsUFxcbLXf3OshgK5VJKSkp4ubNm/r/6tSpIwCIsWPHGr2emZlptF1xcbF44IEHBAARHh4utm/frq+0Z86cEREREQKA8PDwEJGRkWbHTUtL0wdjmzdvLqKjo4UQQuTn54u5c+cKHx8fAUCMGDFCstyxsbEiMDBQABA9evQQFy5cEEKUBJgmT54sNBqNACBmzpyp5sdFpDp762Dfvn31F73+/fuLS5cuCa1WKwoKCsTevXtFhw4d9MuXLFlidlzWQaISY8eO1deVOXPm2LRtYWGhaNWqlb7xqmuAFhcXi1WrVokKFSoIAKJPnz6S27MeEjlWBx3ZVgjX12Eid+FoXZKjJNjKekjkeB18/vnnBQARFhYm1qxZIwoLC4UQQiQkJIgXXnhBv+/FixebbevudZDBViqTdL3orP03cOBAs22vXbsmmjdvrl8nICBAhISE6P/t7e0tfvzxR9ljR0dHi7CwMP36wcHBwtvbW//vxx57TOTl5cluv3HjRhEQEKBfPyQkRHh6eur/PXjwYKHVatX4mIicxt46ePPmTdG+fXujdQICAoSXl5fRa2PHjpU9Nusg3esSEhL0v1cPDw9RrVo1i/998cUXZvuIi4sT4eHhRvXQz89P/++2bdta7FXKekj3MkfqoBr1VwjX12EiV1OrLklREmwVgvWQ7m1q1MGsrCzx4IMP6vfj6+srKlWqZHRfaGnkpDvXQQZbqUxyJNgqhBB5eXli1qxZomPHjiIkJET4+PiI8PBw8frrr4uTJ09aPX5ycrJ47733ROPGjYWfn5+oWLGi6N69u1iwYIFZ93Yply5dEkOHDhXh4eHC19dXVK5cWTz66KNizZo1tn4URC7hSB0sLCwUCxcuFI8//rioVq2a8Pb2Fn5+fqJBgwaif//+Ys+ePVaPzzpI97K4uDhF9c9aIzUjI0NMmDBBtGzZUgQGBorg4GDRvn17MWvWLJGfn2+1HKyHdK9ypA6qVX+FcH0dJnIlNeuSKaXBViFYD+nepVYdLC4uFgsWLBAPPvigCA0NFd7e3qJWrVripZdeEnv37rVaDnetgxohDKb3IiIiIiIiIiIiIiK7eFhfhYiIiIiIiIiIiIisYbCViIiIiIiIiIiISAUMthIRERERERERERGpgMFWIiIiIiIiIiIiIhUw2EpERERERERERESkAgZbiYiIiIiIiIiIiFTAYCsRERERERERERGRChhsJSIiIiIiIiIiIlIBg61EREREREREREREKmCwlYiIiIiIiIiIiEgFDLYSERERERERERERqYDBViIiIiIiIiIiIiIVMNhKREREROTGwsPDodFojP77448/XF0s1Y0aNcrsfQ4aNMjVxSIiIiKyCYOtRERERC4wadIks8CS3H9EAFChQgVUq1YN1apVg5+fn+r7N/xNWhMfH69fd/Hixaoc39nvj4iIiKg0eLm6AERERET3umrVqrm6CFQGfP311+W6p+eUKVMwZcoUAMCgQYOwZMkSF5eIiIiIyHYMthIRERG5WHJysquLQEREREREKmAaASIiIiIiIiIiIiIVMNhKREREVIY8/PDD0Gg0mDRpEgoLCzF79mx06NABFStWhEajwa5du4zWj4+Px6hRo9CiRQsEBQUhICAAzZo1w8iRI3H58mWrx1u2bBm6deuG4OBghISEoFOnTvjxxx8hhMCgQYMkJzHSTehkKZen3LZqlN3wMxJCYMGCBejUqRMqVKiA4OBgdOnSBb/++qvF933lyhWMGzcObdq0QUhICPz9/dGwYUM89dRTWLp0KfLy8gAAqampCAgIgEajwapVqyzu89NPP4VGo0GDBg0ghLC4rr0eeughaDQaTJgwQXad0aNHQ6PRoG/fvqof35ZcxPHx8aofn4iIiMjVmEaAiIiIqAzKy8vDww8/jH379sHLywvBwcFmExstW7YMb7zxBvLz8wEAvr6+8PDwwPnz53H+/HksWrQIa9aswWOPPWa2fyEE3njjDSxatAgAoNFoULFiRURHR+PgwYPYuXMnfH19nfb+HCm7TnFxMZ555hn8+eef8PLyQkBAADIzM7F//37s378fFy9exOTJk822++WXXzBs2DB9QNXHxwfBwcG4fPkyYmNj8ddff6F169Zo06YNKlWqhBdeeAFLlizBjz/+iBdeeEG2LLrPcsiQIU6Z+EwIgWPHjgEA2rVrJ7vekSNHrK5jr6CgIIs5iHNycpCZman6cYmIiIjcBXu2EhEREZVB33//PU6cOIFFixYhIyMDKSkpuHnzJlq3bg0A2LZtGwYMGIDi4mKMGzcOcXFxyM3NRXZ2Ns6dO4d+/fohMzMT/fr1k+wl+u233+qDg++88w5u3LiBlJQUpKSkYNKkSVi5ciX+/PNPp7w3R8tu+Bnt2rULixcvRkZGBtLT03HlyhV9j85p06bh4sWLRtts3LgRAwcORF5eHrp164Y9e/YgNzcXt27dQnZ2Nvbs2YOhQ4fCx8dHv82IESMAADt27EBsbKxkWTZt2oTExER4eXnh9ddfd/QjkhQTE4OMjAwAQPv27WXX0wVk27Ztq3oZxowZg+TkZMn/zp07h9q1awMoCfTWqFFD9eMTERERuRp7thIRERG5WPXq1WWXRUZGokWLFmavZ2Vl4a+//jIaCh4WFgYA0Gq1ePvtt6HVajFv3jwMGzbMaNumTZti1apVeOqpp/DXX39hzpw5+Oqrr/TL8/Ly9D0+X3vtNXz77bf6ZSEhIZg4cSLy8vIwY8YMu96vJY6W3VBqaip27NiBiIgI/Wu1a9fG6tWr0aBBAyQlJWHVqlX4+OOPAQBFRUV49913IYRA9+7dERkZaRRU9fHxQffu3dG9e3ej43Tq1Alt2rTBsWPHsGDBAkyfPt2sLD/++CMA4Mknn7T4fTvi6NGjAIAqVaqgTp06kuvExMQgLS0NgOWerdbKWFxcbFPZioqK8Pzzz+Ps2bOoVasW/vrrL6f2jCYiIiJyFfZsJSIiInKx69evy/5XWFgouU2LFi1kc27u3r0bFy9eROXKlTFkyBDZ4w4YMAAAsGXLFqPXt27dipSUFACQzf05fvx4+Pn5WX1vtnK07Ia6detmFGjV8fX1xeOPPw4AOHHihP71nTt3Ii4uDgDw5ZdfGgVardH1bl28eLHZd5aYmIi///4bADB8+HDF+7SVkvQAuoBsWFgY6tatK7uepd/k9evXcevWLZvKNmLECERGRiIwMBDr169HrVq1bNqeiIiIqKxgz1YiIiIiF7NnsqRu3brJLtu7dy8AID09HTVr1pRdr6CgAACQkJBg9Hp0dDQAoE6dOmjUqJHktiEhIWjfvr3+WGpxtOyGOnXqJLtMt29dUBkA9u3bB6CkV2eHDh2UFxrAK6+8oh9Cv379ejz77LP6ZT///DOKi4tRv359PProozbt1xa6QKqSfK3WUghY+03Gx8ejfv36iso1c+ZMLFy4EB4eHli2bJlT0hcQERERuQsGW4mIiIjKoKpVq8ouS0pKAgAUFhbi+vXrVveVm5tr9O8bN24AgNXeh7r8m2pytOyGgoODZZd5eXnpj6OTnJwMAKhXr56ishoKCgpC//79MW/ePPz444/6YKtWq8VPP/0EABg6dKhTJsbS0QVbLeVrVRpsVcuaNWvw4YcfAigJuj711FOlclwiIiIiV2EaASIiIqIyyNPTU3aZLp9mp06dIIRQ9J+7cGXZHQ2E6lIJbNu2DfHx8QBKUjIkJCTAy8sLgwcPdrSIshITE/VBciVpBCyto5aDBw9iwIABEEJgyJAhGDNmjNOPSURERORqDLYSERERlTO6yY0sDbG3RNdrNjEx0eJ6cst1vUbz8vJkt01PT5d83dGyO8LRY7dq1Qpdu3Y16s26YMECAMBTTz3ltImxgLs9VkNCQmSH98fExOgDss7u2ZqQkIAnn3wSubm5eOSRRzB37lynHo+IiIjIXTDYSkRERFTO6PK5Jicn6/Ov2kKXr/TKlSuIiYmRXCcjIwOHDx+WXFapUiX99lK0Wq1suRwtuyO6du3q8LF1vVt//vlnJCYmYv369QCAYcOGqVNIGboeq9WqVZNdZ+XKlQBKUh40btzYaWXJyMjAE088gevXr6NJkyZYs2YNvL29nXY8IiIiInfCYCsRERFRORMREaGf2Oq9997TTyYlx3CSKAB49NFH9QHTqVOnSm7z+eefy+ZLvf/++wEA69atkxzmv2TJEly9etUpZXdEREQEGjRooPjYUvr164ewsDAkJSXhlVdeQWFhodMnxgLu9my9evWqZI/i+Ph4fPHFFwCA1q1bw8PDObcBRUVF6NevH06dOoWwsDBs3LhR/1siIiIiuhcw2EpERERUznh5eWH+/Pnw8vJCVFQUHnzwQURGRhpNBhUbG4v58+ejY8eOZkO8/f398emnnwIoCYyOGjUKt2/fBlDSa3Hq1Kn47LPPULFiRcnjv/zyywCAs2fPYtiwYUbbfvnll3jzzTcRGhrqlLI7wtPTE9999x00Gg2ioqLQs2dPREVFQavVAgAKCgqwa9cu9O/fH2fOnJHch6+vLwYNGgQA2L17NwDnT4wF3O3ZmpOTg9dee02fCiErKwsrVqxA9+7dkZaWBqDkM05PT5dN5eCIkSNHYuvWrfDx8cHatWv1gXMiIiKiewWDrURERETlUM+ePbF69WoEBwfjwIED6NWrFwIDA1G5cmX4+fmhYcOGGDFiBKKjoyUDgSNHjsRrr70GAPj6669RtWpVhIaGIjQ0FBMmTMCLL74oO7N8z5499dsuXLgQlStXRqVKlVCpUiW8//77GD58OPr27eu0sjuiT58+WLx4MXx9fREVFYUePXogICAAlStXRmBgICIiIrBs2TKLvV7ffPNNfbmcPTEWANy+fRuXL18GAPTt2xdr1qxBeHg4fHx8EBwcjJdffhm1a9dGly5dAJQEgStVqoT8/HzVy7Jx40YAgBACL7zwAqpXry77n1yaCSIiIqKyjMFWIiIionLq6aefxqVLlzBx4kQ88MADCAoKQlpaGnx9fXH//fdjyJAhWLduHcaOHWu2rYeHB5YuXYqlS5eic+fO8Pf3R1FREdq1a4f58+dj+fLlFo+9ePFifP3112jTpg38/f2h1WrRrVs3rFq1Ct98841Ty+6oAQMG4Ny5cxg1ahSaN28OLy8v5Obmol69enj66afxyy+/4L777pPdvlGjRmjTpg0A50+MBdzt1ern54dVq1bhk08+QZ06daDRaNCwYUN89NFH2LFjByZMmIAqVaqgSpUqGD58uH4iNGcoLCzE9evXLf5XXFzstOMTERERuYpGSCXSIiIiIiKyYtCgQViyZAkGDhyIxYsXu7o4biM5ORl16tRBUVERtmzZgscee8yh/YWHhyMhIQGLFi3Spygw9Pnnn+ODDz5Ax44dcfDgQYeO5S742yIiIqKyij1biYiIiIhUNH/+fBQVFaFRo0ZOnxgLuNuztW3btk4/FhERERFZxmArEREREZFKoqOjMXv2bADA+++/r2pO2cGDB0Oj0UCj0eCPP/7Qv37kyBEA0KcuKKtGjRqlf39LlixxdXGIiIiI7OLl6gIQEREREZV14eHhyM/PR3JyMoCSXqZDhgxRZd9VqlRBXl6e0Wt+fn4AgKysLFy6dEl/zLKsQoUKqFatmtFrISEhLioNERERkX0YbCUiIiIiclBCQgIAoHr16ujduzdmzJgBb29vVfZ96NAh2WXHjx+HVquFh4cHWrdurcrxXGXKlCmYMmWKq4tBRERE5BBOkEVERERERERERESkAuZsJSIiIiIiIiIiIlIBg61EREREREREREREKmCwlYiIiIiIiIiIiEgFDLYSERERERERERERqYDBViIiIiIiIiIiIiIVMNhKREREREREREREpAIGW4mIiIiIiIiIiIhUwGArERERERERERERkQoYbCUiIiIiIiIiIiJSAYOtRERERERERERERCr4f2i4oNn8qq+vAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1600x900 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "downsampling = 10\n",
    "\n",
    "f = f[::downsampling]\n",
    "\n",
    "s = s[::downsampling]\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['font.size'] = 18 \n",
    "\n",
    "fig, ax = plt.subplots(figsize=(16,9))\n",
    "ax.plot(f,s)\n",
    "ax.set_xlim(1500, 2600)\n",
    "ax.set_ylim(0, 60)\n",
    "ax.set_ylabel(r'PSD [ppm$^2/\\mu$Hz]')\n",
    "ax.set_xlabel(r'Frequency [$\\mu$Hz]')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a5c04b2",
   "metadata": {},
   "source": [
    "## Observational parameters\n",
    "As with the previous version of PBjam we supply some basic initial parameters to guide the mode ID process. These are $T_{\\mathrm{eff}}$, Gaia $G_{\\mathrm{bp}}-G_{\\mathrm{rp}}$ color index, an estimate of $\\nu_{\\mathrm{max}}$ and $\\Delta \\nu$. You need to look up these parameters for the target you've chosen. \n",
    "\n",
    "[Simbad](https://simbad.u-strasbg.fr/simbad/sim-fid) is very helpful for looking up the $G_{\\mathrm{bp}}-G_{\\mathrm{rp}}$ color information, you can use either the GDR2 or GDR3 values here, since they don't need to be super precise or accurate.\n",
    "\n",
    "Each bit of information is added to a dictionary, where each entry consists of a list-like entry (tuples or arrays are also fine) which contain the mean and error of the relevant parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "26146e1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "obs = {'teff' : (5846, 80), # K\n",
    "       'bp_rp': (0.819174, 0.05),\n",
    "       'numax': (2089.3, 20.0), # muHz\n",
    "       'dnu'  : (95.545, 0.5)} # muHz"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "526e07e6",
   "metadata": {},
   "source": [
    "## Priors\n",
    "With the update to PBjam we now apply a PCA-based dimensionality reduction method to the sample of prior observations for each of the model parameters. This means that the priors are constructed in a lower-dimensional latent space, instead of the full model space. We then compute 1D KDEs of the latent parameter distributions, which are then used in the sampling.\n",
    "\n",
    "This relies on the sample of prior observations being fairly dense in the parameter space around the target star. Due to the increased model complexity ($\\ell=1$ modes for example) PBjam doesn't currently have complete coverage of all parameters in all parts of parameter space. We therefore added an option to include manually specified priors for some of the parameters, in which case they are no longer included in constructing the latent space priors. An example of how to do this is shown at the end of this tutorial. For now let's get going with the mode ID.\n",
    "\n",
    "## Initializing the mode ID sampler\n",
    "The mode identification in PBjam 2.0 is done in two main steps:\n",
    "\n",
    "1. We first identify the $\\ell=2,0$ modes by sampling a simple model of asymptotic p-modes. We also include a background model at this stage. This is almost identical to what was used in the first release of PBjam.\n",
    "\n",
    "2. We then divide the PSD by this model to obtain a residual SNR spectrum which mostly just contains $\\ell=1$ modes. We then apply a choice of three different models to this residual to find the $\\ell=1$ modes. \n",
    "\n",
    "Let's start by initializing the mode ID and running the $\\ell=2,0$ model. The mode ID class takes the frequency and spectrum arrays we defined above, as well as the dictionary of observational parameters. \n",
    "\n",
    "As in the first version of PBjam the choice of the number of radial p-mode orders, $N_p$, is left to the user. Including an excessive amount of orders will probably not improve the result since we rely on the 2nd order asymptotic relation for p-modes. This is only reliable within a few orders of $\\nu_{\\mathrm{max}}$ anyway, and so a suitable number for most stars is between 5-10."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dac16603",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pbjam.modeID import modeID\n",
    "\n",
    "N_p = 7\n",
    "\n",
    "\n",
    "M = modeID(f, s, obs, N_p=N_p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9196ab4",
   "metadata": {},
   "source": [
    "### The $\\ell=2,0$ mode ID\n",
    "With the class initialized we can run the $\\ell=2,0$ model. For those who have used PBjam 1.0 this stage should be fairly familiar, however here we are now using `Dynesty` to do the sampling instead of `emcee`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c32077ac",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "18863it [04:31, 69.38it/s, +600 | bound: 197 | nc: 1 | ncall: 557270 | eff(%):  3.496 | loglstar:   -inf < -7096.259 <    inf | logz: -7125.379 +/-  0.196 | dlogz:  0.000 >  0.100]\n"
     ]
    }
   ],
   "source": [
    "M.runl20model();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7c60558-0271-41b0-8016-259d9e66f33e",
   "metadata": {},
   "source": [
    "After running the l20model the model class instance will be updated to contain a `samples` attribute, which are the likelihood-weighted samples from the nested sampling process. Note that these are samples drawn from the prior that is defined in the latent space, this means that they don't necessarily correspond to any physical parameters, but are instead linear combinations of these parameters. \n",
    "\n",
    "Don't worry though, The mode ID class instance will also be updated to contain the results in a form that is a bit more human readable. These are contained in the `result` attribute, which is just a big dictionary of summary statistics and samples of the model parameters. The parameters have all been converted to linear scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "239da780-3378-42c3-a700-094a3506d9a4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'ell': array([0., 0., 0., 0., 0., 0., 0., 2., 2., 2., 2., 2., 2., 2.]),\n",
       " 'enn': array([17., 18., 19., 20., 21., 22., 23., 16., 17., 18., 19., 20., 21.,\n",
       "        22.]),\n",
       " 'emm': array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]),\n",
       " 'zeta': array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]),\n",
       " 'summary': {'freq': array([[1.75621656e+03, 1.85150156e+03, 1.94694313e+03, 2.04261580e+03,\n",
       "          2.13844582e+03, 2.23440717e+03, 2.33053276e+03, 1.74930398e+03,\n",
       "          1.84459969e+03, 1.94005485e+03, 2.03571457e+03, 2.13155937e+03,\n",
       "          2.22753129e+03, 2.32365879e+03],\n",
       "         [4.55413513e-01, 2.90680532e-01, 1.48138285e-01, 9.99539760e-02,\n",
       "          1.72728390e-01, 2.78185343e-01, 3.69628177e-01, 5.40605848e-01,\n",
       "          3.71327937e-01, 2.20871121e-01, 1.20323316e-01, 1.31073835e-01,\n",
       "          2.07954646e-01, 3.00153762e-01]]),\n",
       "  'height': array([[3.40373183, 5.12217864, 6.36293712, 6.10159078, 4.74017102,\n",
       "          2.99006372, 1.45942334, 2.33236341, 3.55656203, 4.48092988,\n",
       "          4.37493199, 3.44996245, 2.21239206, 1.09984957],\n",
       "         [0.43395854, 0.37763105, 0.41302235, 0.44004288, 0.33133705,\n",
       "          0.1898756 , 0.12195641, 0.31105897, 0.26859343, 0.29480716,\n",
       "          0.31448696, 0.2414944 , 0.13964041, 0.09007912]]),\n",
       "  'width': array([[1.47584377, 1.47584377, 1.47584377, 1.47584377, 1.47584377,\n",
       "          1.47584377, 1.47584377, 1.47584377, 1.47584377, 1.47584377,\n",
       "          1.47584377, 1.47584377, 1.47584377, 1.47584377],\n",
       "         [0.02288559, 0.02288559, 0.02288559, 0.02288559, 0.02288559,\n",
       "          0.02288559, 0.02288559, 0.02288559, 0.02288559, 0.02288559,\n",
       "          0.02288559, 0.02288559, 0.02288559, 0.02288559]]),\n",
       "  'rotAsym': array([[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "         [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]]),\n",
       "  'dnu': array([95.60011937,  0.16415247]),\n",
       "  'numax': array([1979.76227256,   12.76013212]),\n",
       "  'eps_p': array([1.36608988, 0.0361177 ]),\n",
       "  'd02': array([6.90729371, 0.1032822 ]),\n",
       "  'alpha_p': array([0.00171483, 0.00018255]),\n",
       "  'env_width': array([209.33328481,   9.82147204]),\n",
       "  'env_height': array([3.22559739, 0.21508694]),\n",
       "  'mode_width': array([1.47584377, 0.02288559]),\n",
       "  'teff': array([5928.20304439,  110.14789371]),\n",
       "  'bp_rp': array([0.78161292, 0.01161255]),\n",
       "  'H1_nu': array([1197.57921981,  768.32510693]),\n",
       "  'H1_exp': array([2.04652181, 1.83179639]),\n",
       "  'H_power': array([1397.63266683,  178.4591665 ]),\n",
       "  'H2_nu': array([399.62551993,  52.21881081]),\n",
       "  'H2_exp': array([2.94603882, 0.59286184]),\n",
       "  'H3_power': array([8617.41310727, 4108.05309297]),\n",
       "  'H3_nu': array([0.77586561, 0.45106515]),\n",
       "  'H3_exp': array([4.53051833, 1.3349547 ]),\n",
       "  'shot': array([0.83646182, 0.04979915]),\n",
       "  'nurot_e': array([0.02388869, 0.02063556]),\n",
       "  'inc': array([0.74782077, 0.34498241])},\n",
       " 'samples': {'freq': array([[1755.4437615 , 1851.0341552 , 1946.76619107, ..., 2131.66269706,\n",
       "          2227.81965943, 2324.11826396],\n",
       "         [1756.46359128, 1851.7215938 , 1947.15264206, ..., 2131.72876008,\n",
       "          2227.67894554, 2323.80217672],\n",
       "         [1755.41776398, 1850.99829636, 1946.72201926, ..., 2131.59713308,\n",
       "          2227.75042751, 2324.04691246],\n",
       "         ...,\n",
       "         [1756.1572393 , 1851.49169038, 1946.99175917, ..., 2131.50026008,\n",
       "          2227.49718199, 2323.6597216 ],\n",
       "         [1755.3714635 , 1850.95304791, 1946.67543848, ..., 2131.54627802,\n",
       "          2227.69108702, 2323.97670218],\n",
       "         [1756.32998244, 1851.64358947, 1947.12458485, ..., 2131.63882579,\n",
       "          2227.62198624, 2323.77253503]], shape=(5000, 14)),\n",
       "  'height': array([[2.99974512, 4.93146762, 6.32829456, ..., 3.59797248, 2.22161056,\n",
       "          1.06607295],\n",
       "         [3.75432431, 5.35161398, 6.2096966 , ..., 3.27070745, 2.06360308,\n",
       "          1.05510839],\n",
       "         [3.17602407, 5.12527842, 6.50324604, ..., 3.69338995, 2.30383356,\n",
       "          1.12507617],\n",
       "         ...,\n",
       "         [3.32420016, 4.9030869 , 5.89878786, ..., 3.3532533 , 2.20826494,\n",
       "          1.18114715],\n",
       "         [3.0440651 , 4.97299756, 6.35429781, ..., 3.6035255 , 2.22872387,\n",
       "          1.0734619 ],\n",
       "         [3.46532403, 5.08369443, 6.09063803, ..., 3.4461285 , 2.26793395,\n",
       "          1.21373934]], shape=(5000, 14)),\n",
       "  'width': array([[1.45672781, 1.45672781, 1.45672781, ..., 1.45672781, 1.45672781,\n",
       "          1.45672781],\n",
       "         [1.48608894, 1.48608894, 1.48608894, ..., 1.48608894, 1.48608894,\n",
       "          1.48608894],\n",
       "         [1.466041  , 1.466041  , 1.466041  , ..., 1.466041  , 1.466041  ,\n",
       "          1.466041  ],\n",
       "         ...,\n",
       "         [1.49561911, 1.49561911, 1.49561911, ..., 1.49561911, 1.49561911,\n",
       "          1.49561911],\n",
       "         [1.45491654, 1.45491654, 1.45491654, ..., 1.45491654, 1.45491654,\n",
       "          1.45491654],\n",
       "         [1.49643402, 1.49643402, 1.49643402, ..., 1.49643402, 1.49643402,\n",
       "          1.49643402]], shape=(5000, 14)),\n",
       "  'rotAsym': array([[0., 0., 0., ..., 0., 0., 0.],\n",
       "         [0., 0., 0., ..., 0., 0., 0.],\n",
       "         [0., 0., 0., ..., 0., 0., 0.],\n",
       "         ...,\n",
       "         [0., 0., 0., ..., 0., 0., 0.],\n",
       "         [0., 0., 0., ..., 0., 0., 0.],\n",
       "         [0., 0., 0., ..., 0., 0., 0.]], shape=(5000, 14)),\n",
       "  'dnu': array([95.87403099, 95.55556281, 95.86492614, ..., 95.64963924,\n",
       "         95.86225478, 95.62965851], shape=(5000,)),\n",
       "  'numax': array([1994.92423044, 1968.12741788, 1993.30710251, ..., 1985.53530824,\n",
       "         1993.94815689, 1984.22900696], shape=(5000,)),\n",
       "  'eps_p': array([1.30527233, 1.37713362, 1.30675012, ..., 1.35531256, 1.30682785,\n",
       "         1.3609634 ], shape=(5000,)),\n",
       "  'd02': array([6.99249224, 6.80511569, 7.00190351, ..., 6.98848979, 6.99636003,\n",
       "         6.94991488], shape=(5000,)),\n",
       "  'alpha_p': array([0.00147738, 0.00181094, 0.00149367, ..., 0.0017315 , 0.00146884,\n",
       "         0.00175038], shape=(5000,)),\n",
       "  'env_width': array([191.98864748, 209.95047023, 194.85381639, ..., 211.12706095,\n",
       "         192.75034356, 211.72628417], shape=(5000,)),\n",
       "  'env_height': array([3.26527321, 3.12038131, 3.34589155, ..., 2.99895515, 3.27415196,\n",
       "         3.09244314], shape=(5000,)),\n",
       "  'mode_width': array([1.45672781, 1.48608894, 1.466041  , ..., 1.49561911, 1.45491654,\n",
       "         1.49643402], shape=(5000,)),\n",
       "  'teff': array([6128.5365841 , 5913.28097489, 6105.14429563, ..., 5958.23734281,\n",
       "         6114.58172715, 5938.6832663 ], shape=(5000,)),\n",
       "  'bp_rp': array([0.76452778, 0.78528912, 0.76774448, ..., 0.77342595, 0.76722184,\n",
       "         0.77885832], shape=(5000,)),\n",
       "  'H1_nu': array([2667.46562177, 1129.54342317, 2627.14469037, ..., 1202.81322221,\n",
       "         2661.52851126, 1208.45154109], shape=(5000,)),\n",
       "  'H1_exp': array([5.60478316, 2.06504413, 5.51922301, ..., 1.96032355, 5.56950369,\n",
       "         1.92396647], shape=(5000,)),\n",
       "  'H_power': array([1102.22232902, 1385.44272929, 1107.69649263, ..., 1415.76550255,\n",
       "         1110.4882708 , 1447.41890278], shape=(5000,)),\n",
       "  'H2_nu': array([497.01580405, 393.09609184, 498.93554013, ..., 398.51846692,\n",
       "         498.77058953, 399.63319227], shape=(5000,)),\n",
       "  'H2_exp': array([1.82594179, 3.0162425 , 1.9919362 , ..., 2.85604736, 1.91246849,\n",
       "         2.8832552 ], shape=(5000,)),\n",
       "  'H3_power': array([ 8627.41302942, 11807.75245977,  7372.96681695, ...,\n",
       "          4767.81585139,  9198.00871568,  5535.88960394], shape=(5000,)),\n",
       "  'H3_nu': array([1.64206018, 0.57638377, 1.38722135, ..., 0.756377  , 1.38039724,\n",
       "         0.72432482], shape=(5000,)),\n",
       "  'H3_exp': array([2.19633475, 4.92352727, 2.06965272, ..., 4.72773144, 2.19553537,\n",
       "         4.61065602], shape=(5000,)),\n",
       "  'shot': array([0.78081294, 0.84982394, 0.76182058, ..., 0.84635884, 0.76220298,\n",
       "         0.85480019], shape=(5000,)),\n",
       "  'nurot_e': array([0.04680386, 0.01347233, 0.05509917, ..., 0.03266987, 0.05369844,\n",
       "         0.02498345], shape=(5000,)),\n",
       "  'inc': array([1.19734749, 0.48226059, 1.23138016, ..., 0.88380619, 1.20072887,\n",
       "         0.65404719], shape=(5000,))}}"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M.result"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d18b5509",
   "metadata": {},
   "source": [
    "### The $\\ell=1$ mode ID\n",
    "At this stage, if you only really care about the l=2,0 mode pairs you can skip to the detailed peakbagging. Otherwise we can run the $\\ell=1$ model. \n",
    "\n",
    "Here the user must make a choice for which model to apply. The choices of models are:\n",
    "- `MS`: when there are no g-modes near the p-modes. Typically applicable to main-sequence stars. This simply determines the frequencies based on the asymptotic relation for p-modes. \n",
    "- `SG`: when you expect a few g-modes to couple to many p-modes. Typically applicable to sub-giants. Here the mode frequencies are determined by computing the coupling between all combinations of p- and g-modes (see e.g., Ong & Basu 2020).\n",
    "- `RGB`: when you expect lots of g-modes to couple to few p-modes. Typically applicable to low-luminosty red giants. In this model the coupling is computed using the single-value coupling model (e.g. Shibahashi 1979)\n",
    "\n",
    "Exceptions to the above suggestions of course exist, so if in doubt try different models. All the models should run, but results may vary, and also the time it takes to evaluate the model!\n",
    "\n",
    "PBjam does have some limited ability to suggest a model, but take care with this as it is just based on drawing some lines in the ($\\Delta\\nu$, $T_{\\mathrm{eff}}$) plane. So when stars start to transition between main-sequence and sub-giants it is bound to be inaccurate.\n",
    "\n",
    "**Note:** unless you're running the `MS` model, the sampling will probably take a while longer than for the $\\ell=2,0$ since it's a more complicated model. So consider going for a cup of coffee at this point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c3182129",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Input Teff=5846K and dnu=95.545muHz suggests the appropriate l=1 model is: ms\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "1540it [00:06, 234.58it/s, +150 | bound: 44 | nc: 1 | ncall: 28711 | eff(%):  5.917 | loglstar:   -inf < -5958.266 <    inf | logz: -5966.166 +/-  0.207 | dlogz:  0.001 >  0.100]\n",
      "4309it [00:16, 269.31it/s, +150 | bound: 165 | nc: 1 | ncall: 92398 | eff(%):  4.834 | loglstar:   -inf < -5958.258 <    inf | logz: -5966.166 +/-  0.194 | dlogz:  0.000 >  0.000]\n"
     ]
    }
   ],
   "source": [
    "M.runl1model();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96de5c4a",
   "metadata": {},
   "source": [
    "After running the $\\ell=1$ model the result dictionary will be updated to contain the most recent values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "1c72d720",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'ell': array([1., 1., 1., 1., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 2., 2., 2.,\n",
       "        2., 2., 2., 2.]),\n",
       " 'enn': array([-1., -1., -1., -1., -1., -1., -1., 17., 18., 19., 20., 21., 22.,\n",
       "        23., 16., 17., 18., 19., 20., 21., 22.]),\n",
       " 'emm': array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "        0., 0., 0., 0.]),\n",
       " 'zeta': array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]),\n",
       " 'summary': {'freq': array([[1.80086678e+03, 1.89615178e+03, 1.99159334e+03, 2.08726602e+03,\n",
       "          2.18309604e+03, 2.27905739e+03, 2.37518298e+03, 1.75621656e+03,\n",
       "          1.85150156e+03, 1.94694313e+03, 2.04261580e+03, 2.13844582e+03,\n",
       "          2.23440717e+03, 2.33053276e+03, 1.74930398e+03, 1.84459969e+03,\n",
       "          1.94005485e+03, 2.03571457e+03, 2.13155937e+03, 2.22753129e+03,\n",
       "          2.32365879e+03],\n",
       "         [9.85093899e-02, 9.85093899e-02, 9.85093899e-02, 9.85093899e-02,\n",
       "          9.85093899e-02, 9.85093899e-02, 9.85093899e-02, 4.55413513e-01,\n",
       "          2.90680532e-01, 1.48138285e-01, 9.99539760e-02, 1.72728390e-01,\n",
       "          2.78185343e-01, 3.69628177e-01, 5.40605848e-01, 3.71327937e-01,\n",
       "          2.20871121e-01, 1.20323316e-01, 1.31073835e-01, 2.07954646e-01,\n",
       "          3.00153762e-01]]),\n",
       "  'height': array([[5.46272976e+00, 7.26705210e+00, 7.85789748e+00, 6.89811267e+00,\n",
       "          4.91045026e+00, 2.83209815e+00, 1.32186648e+00, 3.40373183e+00,\n",
       "          5.12217864e+00, 6.36293712e+00, 6.10159078e+00, 4.74017102e+00,\n",
       "          2.99006372e+00, 1.45942334e+00, 2.33236341e+00, 3.55656203e+00,\n",
       "          4.48092988e+00, 4.37493199e+00, 3.44996245e+00, 2.21239206e+00,\n",
       "          1.09984957e+00],\n",
       "         [2.19692119e-03, 1.36599489e-03, 2.08875927e-04, 1.66699819e-03,\n",
       "          2.24455906e-03, 1.90554424e-03, 1.17507827e-03, 4.33958544e-01,\n",
       "          3.77631051e-01, 4.13022351e-01, 4.40042881e-01, 3.31337051e-01,\n",
       "          1.89875600e-01, 1.21956414e-01, 3.11058973e-01, 2.68593426e-01,\n",
       "          2.94807164e-01, 3.14486955e-01, 2.41494405e-01, 1.39640406e-01,\n",
       "          9.00791207e-02]]),\n",
       "  'width': array([[1.47584377, 1.47584377, 1.47584377, 1.47584377, 1.47584377,\n",
       "          1.47584377, 1.47584377, 1.47584377, 1.47584377, 1.47584377,\n",
       "          1.47584377, 1.47584377, 1.47584377, 1.47584377, 1.47584377,\n",
       "          1.47584377, 1.47584377, 1.47584377, 1.47584377, 1.47584377,\n",
       "          1.47584377],\n",
       "         [0.02288559, 0.02288559, 0.02288559, 0.02288559, 0.02288559,\n",
       "          0.02288559, 0.02288559, 0.02288559, 0.02288559, 0.02288559,\n",
       "          0.02288559, 0.02288559, 0.02288559, 0.02288559, 0.02288559,\n",
       "          0.02288559, 0.02288559, 0.02288559, 0.02288559, 0.02288559,\n",
       "          0.02288559]]),\n",
       "  'rotAsym': array([[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "          0., 0., 0., 0., 0.],\n",
       "         [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "          0., 0., 0., 0., 0.]]),\n",
       "  'd01': array([44.65021781,  0.09850939]),\n",
       "  'nurot_e': array([0.02388869, 0.02063556]),\n",
       "  'inc': array([0.74782077, 0.34498241]),\n",
       "  'dnu': array([95.60011937,  0.16415247]),\n",
       "  'numax': array([1979.76227256,   12.76013212]),\n",
       "  'eps_p': array([1.36608988, 0.0361177 ]),\n",
       "  'd02': array([6.90729371, 0.1032822 ]),\n",
       "  'alpha_p': array([0.00171483, 0.00018255]),\n",
       "  'env_width': array([209.33328481,   9.82147204]),\n",
       "  'env_height': array([3.22559739, 0.21508694]),\n",
       "  'mode_width': array([1.47584377, 0.02288559]),\n",
       "  'teff': array([5928.20304439,  110.14789371]),\n",
       "  'bp_rp': array([0.78161292, 0.01161255]),\n",
       "  'H1_nu': array([1197.57921981,  768.32510693]),\n",
       "  'H1_exp': array([2.04652181, 1.83179639]),\n",
       "  'H_power': array([1397.63266683,  178.4591665 ]),\n",
       "  'H2_nu': array([399.62551993,  52.21881081]),\n",
       "  'H2_exp': array([2.94603882, 0.59286184]),\n",
       "  'H3_power': array([8617.41310727, 4108.05309297]),\n",
       "  'H3_nu': array([0.77586561, 0.45106515]),\n",
       "  'H3_exp': array([4.53051833, 1.3349547 ]),\n",
       "  'shot': array([0.83646182, 0.04979915])},\n",
       " 'samples': {'freq': array([[1800.91800855, 1896.203005  , 1991.64457353, ..., 2131.66269706,\n",
       "          2227.81965943, 2324.11826396],\n",
       "         [1801.05115263, 1896.33614909, 1991.77771762, ..., 2131.72876008,\n",
       "          2227.67894554, 2323.80217672],\n",
       "         [1800.90252919, 1896.18752564, 1991.62909417, ..., 2131.59713308,\n",
       "          2227.75042751, 2324.04691246],\n",
       "         ...,\n",
       "         [1800.90312155, 1896.188118  , 1991.62968653, ..., 2131.46228486,\n",
       "          2227.40020322, 2323.50202836],\n",
       "         [1800.86431731, 1896.14931377, 1991.5908823 , ..., 2131.50285171,\n",
       "          2227.44096778, 2323.55290985],\n",
       "         [1800.74920869, 1896.03420515, 1991.47577368, ..., 2131.46922559,\n",
       "          2227.43713675, 2323.56232498]], shape=(4459, 21)),\n",
       "  'height': array([[5.46387219, 7.26776224, 7.85778856, ..., 3.59797248, 2.22161056,\n",
       "          1.06607295],\n",
       "         [5.46684096, 7.26960619, 7.85750329, ..., 3.27070745, 2.06360308,\n",
       "          1.05510839],\n",
       "         [5.463527  , 7.2675477 , 7.85782153, ..., 3.69338995, 2.30383356,\n",
       "          1.12507617],\n",
       "         ...,\n",
       "         [5.46354021, 7.26755591, 7.85782026, ..., 3.12164498, 1.98563627,\n",
       "          1.02281077],\n",
       "         [5.46267484, 7.26701796, 7.85790271, ..., 3.74821676, 2.46381313,\n",
       "          1.32016318],\n",
       "         [5.46010752, 7.26542093, 7.85814568, ..., 3.10471094, 2.01350759,\n",
       "          1.05473943]], shape=(4459, 21)),\n",
       "  'width': array([[1.47584377, 1.47584377, 1.47584377, ..., 1.45672781, 1.45672781,\n",
       "          1.45672781],\n",
       "         [1.47584377, 1.47584377, 1.47584377, ..., 1.48608894, 1.48608894,\n",
       "          1.48608894],\n",
       "         [1.47584377, 1.47584377, 1.47584377, ..., 1.466041  , 1.466041  ,\n",
       "          1.466041  ],\n",
       "         ...,\n",
       "         [1.47584377, 1.47584377, 1.47584377, ..., 1.4696376 , 1.4696376 ,\n",
       "          1.4696376 ],\n",
       "         [1.47584377, 1.47584377, 1.47584377, ..., 1.50430534, 1.50430534,\n",
       "          1.50430534],\n",
       "         [1.47584377, 1.47584377, 1.47584377, ..., 1.46569286, 1.46569286,\n",
       "          1.46569286]], shape=(4459, 21)),\n",
       "  'rotAsym': array([[0., 0., 0., ..., 0., 0., 0.],\n",
       "         [0., 0., 0., ..., 0., 0., 0.],\n",
       "         [0., 0., 0., ..., 0., 0., 0.],\n",
       "         ...,\n",
       "         [0., 0., 0., ..., 0., 0., 0.],\n",
       "         [0., 0., 0., ..., 0., 0., 0.],\n",
       "         [0., 0., 0., ..., 0., 0., 0.]], shape=(4459, 21)),\n",
       "  'd01': array([44.7014464 , 44.83459049, 44.68596704, ..., 44.6865594 ,\n",
       "         44.64775516, 44.53264655], shape=(4459,)),\n",
       "  'nurot_e': array([0.04680386, 0.01347233, 0.05509917, ..., 0.03266987, 0.05369844,\n",
       "         0.02498345], shape=(5000,)),\n",
       "  'inc': array([1.19734749, 0.48226059, 1.23138016, ..., 0.88380619, 1.20072887,\n",
       "         0.65404719], shape=(5000,)),\n",
       "  'dnu': array([95.87403099, 95.55556281, 95.86492614, ..., 95.64963924,\n",
       "         95.86225478, 95.62965851], shape=(5000,)),\n",
       "  'numax': array([1994.92423044, 1968.12741788, 1993.30710251, ..., 1985.53530824,\n",
       "         1993.94815689, 1984.22900696], shape=(5000,)),\n",
       "  'eps_p': array([1.30527233, 1.37713362, 1.30675012, ..., 1.35531256, 1.30682785,\n",
       "         1.3609634 ], shape=(5000,)),\n",
       "  'd02': array([6.99249224, 6.80511569, 7.00190351, ..., 6.98848979, 6.99636003,\n",
       "         6.94991488], shape=(5000,)),\n",
       "  'alpha_p': array([0.00147738, 0.00181094, 0.00149367, ..., 0.0017315 , 0.00146884,\n",
       "         0.00175038], shape=(5000,)),\n",
       "  'env_width': array([191.98864748, 209.95047023, 194.85381639, ..., 211.12706095,\n",
       "         192.75034356, 211.72628417], shape=(5000,)),\n",
       "  'env_height': array([3.26527321, 3.12038131, 3.34589155, ..., 2.99895515, 3.27415196,\n",
       "         3.09244314], shape=(5000,)),\n",
       "  'mode_width': array([1.45672781, 1.48608894, 1.466041  , ..., 1.49561911, 1.45491654,\n",
       "         1.49643402], shape=(5000,)),\n",
       "  'teff': array([6128.5365841 , 5913.28097489, 6105.14429563, ..., 5958.23734281,\n",
       "         6114.58172715, 5938.6832663 ], shape=(5000,)),\n",
       "  'bp_rp': array([0.76452778, 0.78528912, 0.76774448, ..., 0.77342595, 0.76722184,\n",
       "         0.77885832], shape=(5000,)),\n",
       "  'H1_nu': array([2667.46562177, 1129.54342317, 2627.14469037, ..., 1202.81322221,\n",
       "         2661.52851126, 1208.45154109], shape=(5000,)),\n",
       "  'H1_exp': array([5.60478316, 2.06504413, 5.51922301, ..., 1.96032355, 5.56950369,\n",
       "         1.92396647], shape=(5000,)),\n",
       "  'H_power': array([1102.22232902, 1385.44272929, 1107.69649263, ..., 1415.76550255,\n",
       "         1110.4882708 , 1447.41890278], shape=(5000,)),\n",
       "  'H2_nu': array([497.01580405, 393.09609184, 498.93554013, ..., 398.51846692,\n",
       "         498.77058953, 399.63319227], shape=(5000,)),\n",
       "  'H2_exp': array([1.82594179, 3.0162425 , 1.9919362 , ..., 2.85604736, 1.91246849,\n",
       "         2.8832552 ], shape=(5000,)),\n",
       "  'H3_power': array([ 8627.41302942, 11807.75245977,  7372.96681695, ...,\n",
       "          4767.81585139,  9198.00871568,  5535.88960394], shape=(5000,)),\n",
       "  'H3_nu': array([1.64206018, 0.57638377, 1.38722135, ..., 0.756377  , 1.38039724,\n",
       "         0.72432482], shape=(5000,)),\n",
       "  'H3_exp': array([2.19633475, 4.92352727, 2.06965272, ..., 4.72773144, 2.19553537,\n",
       "         4.61065602], shape=(5000,)),\n",
       "  'shot': array([0.78081294, 0.84982394, 0.76182058, ..., 0.84635884, 0.76220298,\n",
       "         0.85480019], shape=(5000,))}}"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M.result"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "429f69e8",
   "metadata": {},
   "source": [
    "# 3.5 Specifying manual prior probability densities in the mode ID stage\n",
    "In the event that your target is in a part of parameter space that is poorly sampled it may sometimes be necessary to provide a prior manually for the mode ID to work.\n",
    "\n",
    "In PBjam 2.0 this is done by supplying a distribution class instance from the `pbjam.distributions` module. These function more or less the same as the corresponding `scipy.stats` versions, with the important distinction that they can be incorporated into a function that is jit-compiled with `Jax`. `Jax` is an accelerated linear algebra library that can compile individual functions. `Jax` is continually being developed, so these functions may become redundant in the future. \n",
    "\n",
    "An example could be if we are reasonably certain about the p-mode phase offset $\\epsilon_p$, for a star, and so don't want PBjam to define a prior on its own."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19137483",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pbjam.distributions as dist\n",
    "\n",
    "addPriors = {'eps_p': dist.normal(loc=1.45, scale=0.1)}\n",
    "\n",
    "M_manual = modeID(f, s, obs, N_p=N_p, addPriors=addPriors)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3cdca366-756b-4dc6-bd3e-ef82559915fc",
   "metadata": {},
   "source": [
    "Before running the sampling it is often useful to do what is known as a prior predictive check. Here samples are drawn from the prior distribution and the model is evaluated. You can then check that the prior you have chosen makes sense. Most of the plotting functions in PBjam will have an option to do this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57d1723b-4cef-40b9-bdd0-ccb5b76154bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "M_manual.spectrum(stage='prior')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0066429c-aa86-475d-8238-e03147af73c6",
   "metadata": {},
   "source": [
    "If this looks vaguely sensible you can try running the sampling for that model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa9eab26",
   "metadata": {},
   "outputs": [],
   "source": [
    "M_manual.runl20model();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe3f80f4-0fd5-448d-b7ff-aa9bf4115bee",
   "metadata": {},
   "outputs": [],
   "source": [
    "M.echelle()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acd46b4a-4871-4a50-906d-30c604e59015",
   "metadata": {},
   "source": [
    "**Warning:** You should be careful though since many parameters, like $\\nu_{\\mathrm{max}}$ or $\\Delta \\nu$, need to be specified in terms of their logarithms. You can check the list of model parameters for the l20 model class by inspecting the `variables` attribute. Unfortunately this is currently only accessible once the model class has been initialized, this will probably be changed in future updates."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.12.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}