2024-10-30 One bunch optimizer.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "4a911c31-caf9-4d1d-8688-69200b63de83",
   "metadata": {},
   "outputs": [],
   "source": [
    "from UTILITY_quickstart import *\n",
    "\n",
    "with open('setLattice_defaults.yml', 'r') as file:\n",
    "    importedDefaultSettings = yaml.safe_load(file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f4da901a-3e70-4bec-9469-a9d342d9b74d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Environment set to:  /Users/nmajik/Documents/SLAC/FACET2-Bmad-PyTao\n",
      "Tracking to end\n",
      "CSR off\n",
      "Overwriting lattice with setLattice() defaults\n",
      "No defaults file provided to setLattice(). Using setLattice_defaults.yml\n",
      "Number of macro particles = 1000.0\n",
      "Loaded activeBeamFile.h5\n",
      "Set track_start = L0AFEND, track_end = L0BFEND\n",
      "Tracking!\n",
      "trackBeam() exiting\n"
     ]
    }
   ],
   "source": [
    "csrTF = False\n",
    "evalElement = \"PENT\"\n",
    "\n",
    "\n",
    "inputBeamFilePathSuffix = importedDefaultSettings[\"inputBeamFilePathSuffix\"]\n",
    "tao = initializeTao(\n",
    "    #inputBeamFilePathSuffix = '/beams/nmmToL0AFEND_2bunch_2024-02-16Clean/2024-02-16_2bunch_1e5Downsample_nudgeWeights.h5',\n",
    "    inputBeamFilePathSuffix = inputBeamFilePathSuffix,\n",
    "    \n",
    "    csrTF = csrTF,\n",
    "    numMacroParticles=1e3,\n",
    ")\n",
    "\n",
    "#Set aside the initial beam for later reference\n",
    "trackBeam(tao, trackEnd = \"L0BFEND\", verbose = True)\n",
    "PInit = ParticleGroup(data=tao.bunch_data(\"L0AFEND\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "44e61c34-d532-48b8-8c96-6b1ca6b106fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# PInit = ParticleGroup(\"./beams/2024-10-14_Impact_TwoBunch/2024-10-14_TwoBunch.h5\")\n",
    "# makeBeamActiveBeamFile(PInit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d4958339-61a6-442c-9213-a99470a64549",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using checkpoint: CB1LE\n",
      "centerBC14ForCheckpoint = True, centerBC20ForCheckpoint = True\n",
      "centerBC14ForRemainder = False, centerBC20ForRemainder = False\n"
     ]
    }
   ],
   "source": [
    "targetBunchSpacing = 150e-6\n",
    "\n",
    "#masterToleranceScaling = 1.0 #Higher is looser; generally tighten for early sims and loosen for refinement\n",
    "masterToleranceScalingStart = 1.0\n",
    "masterToleranceScalingEnd = 1.0\n",
    "masterToleranceScalingEvolutionSteps = 1000\n",
    "masterToleranceScaling = masterToleranceScalingStart\n",
    "\n",
    "stageOneOptimization    = 0  # Disabling transverse terms and length limits, per https://docs.google.com/presentation/d/1b6WoEwmDz5cA0fm9FbbGiZoMwBcbCNNNCsD7JiFDxRc/edit#slide=id.g2f91233284b_0_0\n",
    "enableAlignmentTerms    = 0  # Enable centroid offset and angle penalties. Useful for removing final focus kickers from free parameters\n",
    "enableLongitudinalTerms = 0  # Enable bunch spacing penalty and bunch length objectives \n",
    "\n",
    "usingCheckpoint = True\n",
    "#checkpointElement = \"MFFF\"\n",
    "#checkpointElement = \"IM1988\" #Shortly upstream of BEGBC20\n",
    "checkpointElement = \"CB1LE\" #Shortly downstream of BEGBC20; there's an intervening dipole though\n",
    "checkpointElementS = tao.ele_param(checkpointElement,\"ele.s\")['ele_s']\n",
    "centerBC14 = True \n",
    "centerBC20 = True\n",
    "\n",
    "#If the checkpoint is DOWNSTREAM of the centering location AND you want centering, set True\n",
    "centerBC14ForCheckpoint = True and (checkpointElementS > tao.ele_param(\"BEGBC14_1\",\"ele.s\")['ele_s'])\n",
    "centerBC20ForCheckpoint = True and (checkpointElementS > tao.ele_param(\"BEGBC20\",\"ele.s\")['ele_s'])\n",
    "\n",
    "#If the checkpoint is UPSTREAM of the centering location AND you want centering, set True\n",
    "centerBC14ForRemainder = centerBC14 and not centerBC14ForCheckpoint\n",
    "centerBC20ForRemainder = centerBC20 and not centerBC20ForCheckpoint\n",
    "\n",
    "if usingCheckpoint:\n",
    "    print(f\"Using checkpoint: {checkpointElement}\")\n",
    "    print(f\"centerBC14ForCheckpoint = {centerBC14ForCheckpoint}, centerBC20ForCheckpoint = {centerBC20ForCheckpoint}\")\n",
    "    print(f\"centerBC14ForRemainder = {centerBC14ForRemainder}, centerBC20ForRemainder = {centerBC20ForRemainder}\")\n",
    "else:\n",
    "    print(\"No checkpoint. Start to end tracking\")\n",
    "    print(f\"centerBC14 = {centerBC14}, centerBC20 = {centerBC20}\")\n",
    "\n",
    "\n",
    "pbounds = {\n",
    "    # \"QA10361kG\": eval(importedDefaultSettings[\"QA10361kGBounds\"]),\n",
    "    # \"QA10371kG\": eval(importedDefaultSettings[\"QA10371kGBounds\"]),\n",
    "    # \"QE10425kG\": eval(importedDefaultSettings[\"QE10425kGBounds\"]),\n",
    "    # \"QE10441kG\": eval(importedDefaultSettings[\"QE10441kGBounds\"]),\n",
    "    # \"QE10511kG\": eval(importedDefaultSettings[\"QE10511kGBounds\"]),\n",
    "    # \"QE10525kG\": eval(importedDefaultSettings[\"QE10525kGBounds\"]),\n",
    "\n",
    "    # 'L0BPhaseSet': (-30, 30),\n",
    "    # 'L1PhaseSet': (-60, 0),\n",
    "    # 'L2PhaseSet': (-60, 0),\n",
    "    # 'L3PhaseSet': (-15, 15),\n",
    "    \n",
    "    # 'L0BEnergyOffset': (-5e6, 5e6),\n",
    "    # 'L1EnergyOffset': (-10e6, 10e6),\n",
    "    # 'L2EnergyOffset': (-500e6, 500e6),\n",
    "    # 'L3EnergyOffset': (-500e6, 500e6),\n",
    "\n",
    "    # \"Q1EkG\":  eval(importedDefaultSettings[\"Q1EkGBounds\"]),\n",
    "    # \"Q2EkG\":  eval(importedDefaultSettings[\"Q2EkGBounds\"]),\n",
    "    # \"Q3EkG\":  eval(importedDefaultSettings[\"Q3EkGBounds\"]),\n",
    "    # \"Q4EkG\":  eval(importedDefaultSettings[\"Q4EkGBounds\"]),\n",
    "    # \"Q5EkG\":  eval(importedDefaultSettings[\"Q5EkGBounds\"]),\n",
    "    # \"Q6EkG\":  eval(importedDefaultSettings[\"Q6EkGBounds\"]),\n",
    "    \n",
    "    # \"S1ELkG\": eval(importedDefaultSettings[\"S1ELkGBounds\"]),\n",
    "    # \"S2ELkG\": eval(importedDefaultSettings[\"S2ELkGBounds\"]),\n",
    "    # \"S3ELkG\": eval(importedDefaultSettings[\"S3ELkGBounds\"]),\n",
    "    # \"S3ERkG\": eval(importedDefaultSettings[\"S3ERkGBounds\"]),\n",
    "    # \"S2ERkG\": eval(importedDefaultSettings[\"S2ERkGBounds\"]),\n",
    "    # \"S1ERkG\": eval(importedDefaultSettings[\"S1ERkGBounds\"]),\n",
    "\n",
    "    \"S1EL_xOffset\" : ( -0.004, 0.004 ),  \n",
    "    \"S1EL_yOffset\" : ( -0.004, 0.004 ),  \n",
    "    \"S2EL_xOffset\" : ( -0.004, 0.004 ),  \n",
    "    \"S2EL_yOffset\" : ( -0.004, 0.004 ),  \n",
    "    \"S2ER_xOffset\" : ( -0.004, 0.004 ),  \n",
    "    \"S2ER_yOffset\" : ( -0.004, 0.004 ),  \n",
    "    \"S1ER_xOffset\" : ( -0.004, 0.004 ),  \n",
    "    \"S1ER_yOffset\" : ( -0.004, 0.004 ),\n",
    "\n",
    "    # # 'Q5FFkG': eval(importedDefaultSettings[\"Q5FFkGBounds\"]),\n",
    "    # # 'Q4FFkG': eval(importedDefaultSettings[\"Q4FFkGBounds\"]),\n",
    "    # # 'Q3FFkG': eval(importedDefaultSettings[\"Q3FFkGBounds\"]),\n",
    "    # # 'Q2FFkG': eval(importedDefaultSettings[\"Q2FFkGBounds\"]),\n",
    "    # # 'Q1FFkG': eval(importedDefaultSettings[\"Q1FFkGBounds\"]),\n",
    "    # # 'Q0FFkG': eval(importedDefaultSettings[\"Q0FFkGBounds\"]),\n",
    "    # # # 'Q0DkG':  eval(importedDefaultSettings[\"Q0DkGBounds\"]),\n",
    "    # # # 'Q1DkG':  eval(importedDefaultSettings[\"Q1DkGBounds\"]),\n",
    "    # # # 'Q2DkG':  eval(importedDefaultSettings[\"Q2DkGBounds\"]),\n",
    "\n",
    "    \"XC1FFkG\" : tuple(2 * x for x in eval(importedDefaultSettings[\"XC1FFkGBounds\"])), #2024-10-11: Extending bounds as proxy for tuning final chicane dipole strength\n",
    "    \"XC3FFkG\" : eval(importedDefaultSettings[\"XC3FFkGBounds\"]),\n",
    "    \"YC1FFkG\" : eval(importedDefaultSettings[\"YC1FFkGBounds\"]),\n",
    "    \"YC2FFkG\" : eval(importedDefaultSettings[\"YC2FFkGBounds\"]),\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1c8c357-a2a7-4e99-bdb1-4e53ed6e5393",
   "metadata": {},
   "source": [
    "# Optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "704ad00f-d1d2-42f6-a734-50be8ad97a45",
   "metadata": {},
   "source": [
    "## Utility functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "63ed0772-8ef8-4cdb-a890-83a3d7c5fd1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# #2024-08-23: Now to IM1988 to enable centering at BC20\n",
    "# def propagateToBEGBC20(): \n",
    "\n",
    "#     ##################################\n",
    "#     #Propagate to BEGBC20 and save result\n",
    "#     ##################################\n",
    "#     makeBeamActiveBeamFile(PInit)\n",
    "    \n",
    "\n",
    "#     trackBeam(tao, trackStart = \"L0AFEND\", trackEnd = \"IM1988\", centerBC14 = True, verbose = True)\n",
    "\n",
    "#     P = ParticleGroup(data=tao.bunch_data(\"IM1988\"))\n",
    "    \n",
    "#     makeBeamActiveBeamFile(P)\n",
    "\n",
    "#     return\n",
    "\n",
    "\n",
    "#Generalizing\n",
    "def propagateFromStartToCheckpoint(\n",
    "    checkpointElement,\n",
    "    centerBC14,\n",
    "    centerBC20\n",
    "): \n",
    "    makeBeamActiveBeamFile(PInit)\n",
    "\n",
    "    trackBeam(tao, trackStart = \"L0AFEND\", trackEnd = checkpointElement, centerBC14 = centerBC14, centerBC20 = centerBC20, verbose = False)\n",
    "\n",
    "    P = ParticleGroup(data=tao.bunch_data(checkpointElement))\n",
    "    \n",
    "    makeBeamActiveBeamFile(P)\n",
    "\n",
    "    return"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac1e0d14-7eb3-4bcc-ba7d-0c71a5a83895",
   "metadata": {},
   "source": [
    "## Optimizer functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "92e8cd7a-9adb-4abd-a60d-010d2ca631a4",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if usingCheckpoint:\n",
    "    setLattice(tao) #Set lattice to current default config\n",
    "    propagateFromStartToCheckpoint(checkpointElement, centerBC14ForCheckpoint, centerBC20ForCheckpoint)\n",
    "\n",
    "def rampToZero(val, thresh, scale = 1):\n",
    "    return (max(val, thresh) - thresh) / scale\n",
    "\n",
    "def rampToZeroFlip(val, thresh, scale = 1):\n",
    "    #This ensures some minimum val is reached and linearly penalizes values which are lower\n",
    "    return (max( -1 * val, -1 * thresh) + thresh) / scale\n",
    "\n",
    "\n",
    "def updateMasterToleranceScaling(totalNumEvals):\n",
    "    \"\"\"\n",
    "    I'm not sure if this is actually a good idea or not. \n",
    "    The general idea I'm aiming for is that early optimization iterations will get closer to the thresholded constraints than strictly necessary, \n",
    "    allowing later optimization to proceed with more wiggle room and less likelihood of falling off a \"cliff\".\n",
    "    Maybe it'd make more sense to make this conditional on something other than evaluation count?\n",
    "    Alternatively, what if it's useful to start with looser constraints and gradually tighten them? For something like COBYQA, it's less likely to \"overlearn\" when trying to satisfy, e.g. the first condition\n",
    "    An insane person might try using an exponentially decaying sine to try both...\n",
    "    \"\"\"\n",
    "    \n",
    "    global masterToleranceScaling\n",
    "\n",
    "    if totalNumEvals < masterToleranceScalingEvolutionSteps: \n",
    "        masterToleranceScaling = masterToleranceScalingStart + (masterToleranceScalingEnd - masterToleranceScalingStart) * totalNumEvals / masterToleranceScalingEvolutionSteps\n",
    "    else:\n",
    "        masterToleranceScaling = masterToleranceScalingEnd\n",
    "    \n",
    "    return\n",
    "\n",
    "\n",
    "\n",
    "def specificOptimizer(\n",
    "    self,\n",
    "    **kwargs\n",
    "):\n",
    "\n",
    "    self.totalNumEvals += 1\n",
    "    self.displayEvals()\n",
    "\n",
    "    updateMasterToleranceScaling(self.totalNumEvals)\n",
    "\n",
    "    \n",
    "    savedData = kwargs\n",
    "    \n",
    "    badValue = -1e30  #The value returned for illegal config. Should be colossal. Double limit ~= 1e308\n",
    "    bigCost  = 1e20   #Should be large enough to dominate any \"normal\" return value but be dominated by badValue\n",
    "    \n",
    "    try: #This try block deals with bad configurations. Instead of causing the optimizer to halt we now 'except' a low value\n",
    "        setLattice(tao, **kwargs)\n",
    "\n",
    "    except:\n",
    "        print(f\"specificOptimizer() excepted'd on setLattice()\")\n",
    "        return badValue * 5\n",
    "\n",
    "    try:\n",
    "        if usingCheckpoint: \n",
    "            trackBeam(tao, trackStart = checkpointElement, trackEnd = evalElement, centerBC14 = centerBC14ForRemainder, centerBC20 = centerBC20ForRemainder, verbose = False)\n",
    "        else:\n",
    "            trackBeam(tao, trackStart = \"L0AFEND\", trackEnd = evalElement, centerBC14 = centerBC14, centerBC20 = centerBC20, verbose = False)\n",
    "\n",
    "    except:\n",
    "            print(f\"specificOptimizer() excepted'd on trackBeam()\")\n",
    "            return badValue * 4\n",
    "    \n",
    "\n",
    "    #BEGBC20NumLiveParticles = tao.bunch_params(\"BEGBC20\")['n_particle_live']\n",
    "    #PENTNumLiveParticles = tao.bunch_params(\"PENT\")['n_particle_live']\n",
    "\n",
    "    if tao.bunch_params(evalElement)['n_particle_live'] < 10:\n",
    "        print(f\"specificOptimizer() got ~no particles after tracking\")\n",
    "        return badValue * 2 \n",
    "\n",
    "    \n",
    "\n",
    "    \n",
    "    try: \n",
    "        P = getBeamAtElement(tao, evalElement)\n",
    "        #PDrive, PWitness = getDriverAndWitness(P)\n",
    "        PDrive = P.copy()\n",
    "    \n",
    "\n",
    "        for PActiveStr in [\"PDrive\"]: #[\"PDrive\", \"PWitness\"]: #, \"PSystem\"]:\n",
    "            PActive = locals()[PActiveStr]\n",
    "            for val in [\"mean_x\", \"mean_y\", \"sigma_x\", \"sigma_y\", \"mean_xp\", \"mean_yp\"]:\n",
    "                savedData[f\"{PActiveStr}_{val}\"] = PActive[val]\n",
    "                #print(f\"\"\"{PActiveStr}_{val} = {PActive[val]}\"\"\")\n",
    "\n",
    "\n",
    "            savedData[f\"{PActiveStr}_median_x\"] = np.median(PActive.x)\n",
    "            savedData[f\"{PActiveStr}_median_y\"] = np.median(PActive.y)\n",
    "\n",
    "            savedData[f\"{PActiveStr}_median_xp\"] = np.median(PActive.xp)\n",
    "            savedData[f\"{PActiveStr}_median_yp\"] = np.median(PActive.yp)\n",
    "            \n",
    "            savedData[f\"{PActiveStr}_sigmaSI90_x\"] = smallestIntervalImpliedSigma(PActive.x, percentage = 0.90)\n",
    "            savedData[f\"{PActiveStr}_sigmaSI90_y\"] = smallestIntervalImpliedSigma(PActive.y, percentage = 0.90)\n",
    "            savedData[f\"{PActiveStr}_sigmaSI90_z\"] = smallestIntervalImpliedSigma(PActive.t * 3e8, percentage=0.9)\n",
    "\n",
    "            savedData[f\"{PActiveStr}_sigmaSI90_xp\"] = smallestIntervalImpliedSigma(PActive.xp, percentage = 0.90)\n",
    "            savedData[f\"{PActiveStr}_sigmaSI90_yp\"] = smallestIntervalImpliedSigma(PActive.yp, percentage = 0.90)\n",
    "\n",
    "            savedData[f\"{PActiveStr}_emitSI90_x\"] = smallestIntervalImpliedEmittance(PActive, plane = \"x\", percentage = 0.90)\n",
    "            savedData[f\"{PActiveStr}_emitSI90_y\"] = smallestIntervalImpliedEmittance(PActive, plane = \"y\", percentage = 0.90)\n",
    "\n",
    "            savedData[f\"{PActiveStr}_norm_emit_x\"] = PActive[\"norm_emit_x\"]\n",
    "            savedData[f\"{PActiveStr}_norm_emit_y\"] = PActive[\"norm_emit_y\"]\n",
    "\n",
    "            #savedData[f\"{PActiveStr}_zLen\"] = smallestIntervalImpliedSigma(PActive.t * 3e8, percentage=0.9)\n",
    "\n",
    "            savedData[f\"{PActiveStr}_zCentroid\"] = np.median(PActive.t * 3e8)\n",
    "\n",
    "        # savedData[\"bunchSpacing\"] = savedData[\"PWitness_zCentroid\"] - savedData[\"PDrive_zCentroid\"]\n",
    "\n",
    "        # savedData[\"transverseCentroidOffset\"] = np.sqrt(\n",
    "        #         (savedData[\"PDrive_mean_x\"] - savedData[\"PWitness_mean_x\"])**2 + \n",
    "        #         (savedData[\"PDrive_mean_y\"] - savedData[\"PWitness_mean_y\"])**2\n",
    "        #     )\n",
    "\n",
    "        \n",
    "        #savedData[\"lostChargeFraction\"] = 1 - (tao.bunch_params(\"DTOTR\")['n_particle_live'] / tao.bunch_params(\"BEGBC20\")['n_particle_live'])\n",
    "        savedData[\"lostChargeFraction\"] = 1 - (P.charge / PInit.charge)\n",
    "            \n",
    "\n",
    "    except:\n",
    "        print(f\"specificOptimizer() excepted'd while getting beam and compiling savedData\")\n",
    "        return badValue\n",
    "\n",
    "\n",
    "    if stageOneOptimization: \n",
    "        enableTransverse = 0\n",
    "        lengthLimitMultiplier = 0\n",
    "        stageOneSpacingToleranceMultiplier = 0.5\n",
    "\n",
    "    else:\n",
    "        enableTransverse = 1\n",
    "        lengthLimitMultiplier = 1\n",
    "\n",
    "    \n",
    "\n",
    "    \n",
    "    tolerableBeamLossFraction  = 0.02  * masterToleranceScaling\n",
    "    tolerableBunchSpacingError = 100e-6 * masterToleranceScaling * (stageOneSpacingToleranceMultiplier if stageOneOptimization else 1.0)\n",
    "    \n",
    "    tolerableBeamOffset  = 5e-6 * masterToleranceScaling\n",
    "    tolerableAngleOffset = 5e-3 * masterToleranceScaling\n",
    "    \n",
    "\n",
    "    driveEmittanceThreshold   = 50e-6 * masterToleranceScaling\n",
    "    witnessEmittanceThreshold = 20e-6 * masterToleranceScaling\n",
    "\n",
    "    driveSpotThreshold     = 10e-6 #* masterToleranceScaling\n",
    "    witnessSpotThreshold   = 30e-6 #* masterToleranceScaling\n",
    "    \n",
    "    driveLengthThreshold   = lengthLimitMultiplier * 10e-6 #* masterToleranceScaling\n",
    "    witnessLengthThreshold = lengthLimitMultiplier * 20e-6 #* masterToleranceScaling\n",
    "\n",
    "\n",
    "    savedData[\"errorTerm_lostChargeFraction\"] = 1e3 * rampToZero( savedData[\"lostChargeFraction\"], tolerableBeamLossFraction, scale = 0.01)**2\n",
    "    \n",
    "    #savedData[\"errorTerm_bunchSpacing\"] = enableLongitudinalTerms * 1e3 * rampToZero( abs(savedData[\"bunchSpacing\"] - targetBunchSpacing), tolerableBunchSpacingError, scale = 10e-6)**2\n",
    "    \n",
    "    savedData[\"errorTerm_transverseOffset\"] = enableAlignmentTerms * enableTransverse * 1e3 * np.mean([\n",
    "                #rampToZero(abs(savedData[\"PWitness_median_x\"]), tolerableBeamOffset, scale = 1e-6) ** 2,\n",
    "                #rampToZero(abs(savedData[\"PWitness_median_y\"]), tolerableBeamOffset, scale = 1e-6) ** 2,\n",
    "                rampToZero(abs(savedData[\"PDrive_median_x\"]  ), tolerableBeamOffset, scale = 1e-6) ** 2,\n",
    "                rampToZero(abs(savedData[\"PDrive_median_y\"]  ), tolerableBeamOffset, scale = 1e-6) ** 2,\n",
    "    ])\n",
    "    \n",
    "    savedData[\"errorTerm_angleOffset\"] = enableAlignmentTerms * enableTransverse * 1e3 * np.mean([\n",
    "                #rampToZero(abs(savedData[\"PWitness_median_xp\"]), tolerableAngleOffset, scale = 100e-6) ** 2,\n",
    "                #rampToZero(abs(savedData[\"PWitness_median_yp\"]), tolerableAngleOffset, scale = 100e-6) ** 2,\n",
    "                rampToZero(abs(savedData[\"PDrive_median_xp\"]  ), tolerableAngleOffset, scale = 100e-6) ** 2,\n",
    "                rampToZero(abs(savedData[\"PDrive_median_yp\"]  ), tolerableAngleOffset, scale = 100e-6) ** 2,\n",
    "    ])\n",
    "    \n",
    "    # savedData[\"errorTerm_mainObjective\"] = max([\n",
    "    #             savedData[f\"PDrive_sigmaSI90_x\"],savedData[f\"PDrive_sigmaSI90_y\"],savedData[f\"PDrive_zLen\"],\n",
    "    #             savedData[f\"PWitness_sigmaSI90_x\"],savedData[f\"PWitness_sigmaSI90_y\"],savedData[f\"PWitness_zLen\"]\n",
    "    # ]) / 10e-6\n",
    "\n",
    "    #2024-11-25-15-58-15: A desperate, ad hoc move... Trying to prevent optimizer from \"rolling up\" phase space at PENT\n",
    "    #savedData[\"errorTerm_sigma_xp_rule\"] = 1e3 * ( 1 / savedData[f\"PDrive_sigmaSI90_xp\"] if savedData[f\"PDrive_sigmaSI90_xp\"] < 400e-6 else 0) ** 2\n",
    "    savedData[\"errorTerm_sigma_xp_rule\"] = 1e3 * ( rampToZeroFlip(savedData[f\"PDrive_sigmaSI90_xp\"], 400e-6, 10e-6) ) ** 2\n",
    "    \n",
    "    savedData[\"errorTerm_mainObjective\"] = np.mean([\n",
    "                enableTransverse *        rampToZero(savedData[f\"PDrive_sigmaSI90_x\"],   driveSpotThreshold,     10e-6) ** 2,\n",
    "                enableTransverse *        rampToZero(savedData[f\"PDrive_sigmaSI90_y\"],   driveSpotThreshold,     10e-6) ** 2,\n",
    "                enableLongitudinalTerms * rampToZero(savedData[f\"PDrive_sigmaSI90_z\"],   driveLengthThreshold,   10e-6) ** 2,\n",
    "\n",
    "                #enableTransverse *                          rampToZero(savedData[f\"PWitness_sigmaSI90_x\"], witnessSpotThreshold,   10e-6) ** 2,\n",
    "                #enableTransverse *                          rampToZero(savedData[f\"PWitness_sigmaSI90_y\"], witnessSpotThreshold,   10e-6) ** 2,\n",
    "                #enableLongitudinalTerms *                   rampToZero(savedData[f\"PWitness_sigmaSI90_z\"], witnessLengthThreshold, 10e-6) ** 2,\n",
    "\n",
    "                # rampToZero(savedData[f\"PDrive_emitSI90_x\"], 0, 10e-6) ** 2,\n",
    "                # rampToZero(savedData[f\"PDrive_emitSI90_y\"], 0, 10e-6) ** 2,\n",
    "\n",
    "                # rampToZero(savedData[f\"PDrive_norm_emit_x\"], 0, 10e-6) ** 2,\n",
    "                # rampToZero(savedData[f\"PDrive_norm_emit_y\"], 0, 10e-6) ** 2,\n",
    "\n",
    "    ]) \n",
    "\n",
    "    #Secondary objective includes all \"ramp\" terms with thresholds disabled. Intended to gently nudge all specs to better values if nothing else is going on; mostly expect this to do anything once thresholds are hit\n",
    "    #2024-10-15 comment: Advise not setting this weight above 1e-6. For quite-good mainObjective settings, even 1e-4 is too much\n",
    "    savedData[\"errorTerm_secondaryObjective\"] = 1e-6 * np.mean([\n",
    "                #rampToZero( abs(savedData[\"bunchSpacing\"] - targetBunchSpacing), 0 * tolerableBunchSpacingError, scale = 10e-6)**2,\n",
    "        \n",
    "                #rampToZero(abs(savedData[\"PWitness_median_x\"]),                  0 * tolerableBeamOffset, scale = 1e-6) ** 2,\n",
    "                #rampToZero(abs(savedData[\"PWitness_median_y\"]),                  0 * tolerableBeamOffset, scale = 1e-6) ** 2,\n",
    "                rampToZero(abs(savedData[\"PDrive_median_x\"]  ),                  0 * tolerableBeamOffset, scale = 1e-6) ** 2,\n",
    "                rampToZero(abs(savedData[\"PDrive_median_y\"]  ),                  0 * tolerableBeamOffset, scale = 1e-6) ** 2,\n",
    "\n",
    "                #rampToZero(abs(savedData[\"PWitness_median_xp\"]),                 0 * tolerableAngleOffset, scale = 100e-6) ** 2,\n",
    "                #rampToZero(abs(savedData[\"PWitness_median_yp\"]),                 0 * tolerableAngleOffset, scale = 100e-6) ** 2,\n",
    "                rampToZero(abs(savedData[\"PDrive_median_xp\"]  ),                 0 * tolerableAngleOffset, scale = 100e-6) ** 2,\n",
    "                rampToZero(abs(savedData[\"PDrive_median_yp\"]  ),                 0 * tolerableAngleOffset, scale = 100e-6) ** 2,\n",
    "\n",
    "                rampToZero(savedData[f\"PDrive_sigmaSI90_x\"],                     0 * driveSpotThreshold,     10e-6) ** 2,\n",
    "                rampToZero(savedData[f\"PDrive_sigmaSI90_y\"],                     0 * driveSpotThreshold,     10e-6) ** 2,\n",
    "                rampToZero(savedData[f\"PDrive_sigmaSI90_z\"],                     0 * driveLengthThreshold,   10e-6) ** 2,\n",
    "                #rampToZero(savedData[f\"PWitness_sigmaSI90_x\"],                   0 * witnessSpotThreshold,   10e-6) ** 2,\n",
    "                #rampToZero(savedData[f\"PWitness_sigmaSI90_y\"],                   0 * witnessSpotThreshold,   10e-6) ** 2,\n",
    "                #rampToZero(savedData[f\"PWitness_sigmaSI90_z\"],                   0 * witnessLengthThreshold, 10e-6) ** 2,\n",
    "    ])\n",
    "    \n",
    "    savedData[\"maximizeMe\"] = 1 / np.mean([\n",
    "        savedData[\"errorTerm_lostChargeFraction\"], \n",
    "        #savedData[\"errorTerm_bunchSpacing\"],\n",
    "        savedData[\"errorTerm_transverseOffset\"],\n",
    "        savedData[\"errorTerm_angleOffset\"],\n",
    "        savedData[\"errorTerm_sigma_xp_rule\"],\n",
    "        savedData[\"errorTerm_mainObjective\"],\n",
    "        savedData[\"errorTerm_secondaryObjective\"],\n",
    "        1e-20 #Avoid infinities \n",
    "    ])\n",
    "\n",
    "    savedData[\"inputBeamFilePathSuffix\"] = inputBeamFilePathSuffix\n",
    "    savedData[\"csrTF\"] = csrTF\n",
    "    \n",
    "\n",
    "    \n",
    "    #Collect desired data as a pandas Series\n",
    "    tmpData = pd.Series( savedData ) \n",
    "    self.history = pd.concat([self.history, tmpData.to_frame().T])\n",
    "\n",
    "    #Optional: Write to file\n",
    "    self.history.to_json('optimizerHistory.json', orient='records')\n",
    "    \n",
    "    self.updatePlot()\n",
    "\n",
    "    return savedData[\"maximizeMe\"]\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# Function to create optimizerWrapper based on pbounds\n",
    "def create_optimizer_wrapper(pbounds):\n",
    "    param_names = list(pbounds.keys())\n",
    "    \n",
    "    def optimizerWrapper(self, **kwargs):\n",
    "        params = {name: kwargs.get(name, None) for name in param_names}\n",
    "        if None in params.values():\n",
    "            raise ValueError(\"All parameters must be provided\")\n",
    "        return specificOptimizer(self, **params)\n",
    "    \n",
    "    return optimizerWrapper\n",
    "\n",
    "# Create the optimizerWrapper function\n",
    "optimizerWrapper = create_optimizer_wrapper(pbounds)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "125b726c-ad45-4c40-870a-5d34b3d0d194",
   "metadata": {
    "editable": true,
    "jp-MarkdownHeadingCollapsed": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "### General optimizer setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b7955b78-5c54-4205-80b8-cf73062fa36b",
   "metadata": {
    "editable": true,
    "jupyter": {
     "source_hidden": true
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "class OptimizationProblem:\n",
    "    def __init__(self):\n",
    "        self.history = pd.DataFrame()\n",
    "        self.totalNumEvals = 0\n",
    "        self.plot_display_handle = None\n",
    "        self.evals_display_handle = None\n",
    "\n",
    "\n",
    "    def updatePlot(self):\n",
    "        plt.figure()\n",
    "        plotKey = \"maximizeMe\"\n",
    "        plt.plot(np.arange(len(self.history[plotKey])), self.history[plotKey], '-')\n",
    "        \n",
    "        plt.title('Optimization History')\n",
    "        plt.xlabel('Evaluation #')\n",
    "        plt.ylabel(plotKey)\n",
    "        plt.yscale('log')\n",
    "        plt.grid(True)\n",
    "        \n",
    "        \n",
    "        if self.plot_display_handle is None:\n",
    "            self.plot_display_handle = display(plt.gcf(), display_id=True)\n",
    "        else:\n",
    "            update_display(plt.gcf(), display_id=self.plot_display_handle.display_id)\n",
    "        plt.close()\n",
    "\n",
    "    def displayEvals(self):\n",
    "        if self.evals_display_handle is None:\n",
    "            self.evals_display_handle = display(f\"Total Num Evals: {self.totalNumEvals}\", display_id=True)\n",
    "        else:\n",
    "            #Can't use '\\n' in this context. <br> doesn't work either\n",
    "            update_display(f\"Total Num Evals: {self.totalNumEvals}, masterToleranceScaling: {masterToleranceScaling}\", display_id=self.evals_display_handle.display_id)\n",
    "\n",
    "# Attach the function to the class as a method\n",
    "OptimizationProblem.optimizerWrapper = optimizerWrapper\n",
    "\n",
    "# Instantiate the optimization problem\n",
    "problem = OptimizationProblem()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "636c9a69-ac6f-406c-b7c7-1a0801da07bc",
   "metadata": {
    "editable": true,
    "jupyter": {
     "source_hidden": true
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from scipy.optimize import differential_evolution\n",
    "\n",
    "# Convert pbounds to the format required by differential_evolution\n",
    "bounds = [(low, high) for (low, high) in pbounds.values()]\n",
    "param_names = list(pbounds.keys())\n",
    "\n",
    "#scipy.optimize wants to optimize a function which is passed a vector of all the parameters\n",
    "#This programmatically wraps the existing wrapper (ugh...) to handle this format\n",
    "def create_DE_wrapper(obj, param_names):\n",
    "    def wrapperDE(params):\n",
    "        param_dict = dict(zip(param_names, params))\n",
    "\n",
    "        #bayes_opt is a maximizer but differential_evolution is a minimizer... hence the inversion\n",
    "        return -1*obj.optimizerWrapper(**param_dict)\n",
    "    \n",
    "    return wrapperDE\n",
    "\n",
    "wrapperDE = create_DE_wrapper(problem, param_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e86c1a72-1d2e-4dd7-9056-c34bfa63a40f",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "defaultSettingsVector = [importedDefaultSettings[key] for key in list(pbounds.keys())]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "477c9433-724f-429d-9ea4-54c010f39ffe",
   "metadata": {},
   "source": [
    "## Run optimizer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c43e77d2-b6ce-4963-aff9-239f373446f8",
   "metadata": {
    "editable": true,
    "jupyter": {
     "source_hidden": true
    },
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<span style=\"color: red; font-size: 18px;\">Checkpoint enabled: CB1LE</span>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<span style=\"color: green; font-size: 18px;\">No masterToleranceScaling</span>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<span style=\"color: red; font-size: 18px;\">CSR disabled</span>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<span style=\"color: green; font-size: 18px;\">Evaluation at PENT</span>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<span style=\"color: red; font-size: 18px;\">Offset and angle penalties disabled</span>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<span style=\"color: red; font-size: 18px;\">Bunch spacing and length terms disabled</span>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.display import HTML, display\n",
    "\n",
    "if usingCheckpoint:\n",
    "    display(HTML(f'<span style=\"color: red; font-size: 18px;\">Checkpoint enabled: {checkpointElement}</span>'))\n",
    "else:\n",
    "    display(HTML('<span style=\"color: green; font-size: 18px;\">No checkpoint. Start-to-end simulation</span>'))\n",
    "\n",
    "if masterToleranceScalingStart == masterToleranceScalingEnd == 1.0:\n",
    "    display(HTML('<span style=\"color: green; font-size: 18px;\">No masterToleranceScaling</span>'))\n",
    "else:\n",
    "    display(HTML('<span style=\"color: red; font-size: 18px;\">masterToleranceScaling is changing</span>'))\n",
    "\n",
    "if csrTF:\n",
    "    display(HTML('<span style=\"color: green; font-size: 18px;\">CSR enabled</span>'))\n",
    "else:\n",
    "    display(HTML('<span style=\"color: red; font-size: 18px;\">CSR disabled</span>'))\n",
    "\n",
    "if evalElement == \"PENT\":\n",
    "    display(HTML('<span style=\"color: green; font-size: 18px;\">Evaluation at PENT</span>'))\n",
    "else:\n",
    "    display(HTML(f'<span style=\"color: red; font-size: 18px;\">Evaluation at {evalElement}, not PENT</span>'))\n",
    "\n",
    "if stageOneOptimization:\n",
    "    display(HTML('<span style=\"color: red; font-size: 18px;\">Stage 1 optimization disables spot size, angle, and offset penalties </span>'))\n",
    "\n",
    "if not enableAlignmentTerms:\n",
    "    display(HTML('<span style=\"color: red; font-size: 18px;\">Offset and angle penalties disabled</span>'))\n",
    "\n",
    "if not enableLongitudinalTerms:\n",
    "    display(HTML('<span style=\"color: red; font-size: 18px;\">Bunch spacing and length terms disabled</span>'))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "294e9bdb-4ae7-46c9-b0fd-9f729563b7da",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Total Num Evals: 296, masterToleranceScaling: 1.0'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAHFCAYAAAAe+pb9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAACL8klEQVR4nO2dd5wTdf7/X5OyvQG7sLv0Ll2aCIKCHigop6hgP6xfC3pnO89y/sTzvqdfvbPcudh7b2C5UxAVBKS3owpKW8rSlrK9JJnfH9mZzEwmySSZSdl9PR8P2GQy+Xw++Ww288q7CqIoiiCEEEIIaYHY4r0AQgghhJB4QSFECCGEkBYLhRAhhBBCWiwUQoQQQghpsVAIEUIIIaTFQiFECCGEkBYLhRAhhBBCWiwUQoQQQghpsVAIEUIIIaTFQiFESBKyfPlyTJ06FUVFRUhJSUFhYSEuvfRSLFu2LKpxZ82ahTfffNPv+O7duyEIgu5jkWLFmEq2bNmCmTNnYvfu3X6PXXvttejSpYsl84ZCEATcfvvtuo99+umnEAQBCxculI/NnDkTgiCENUdNTQ1mzpypGocQog+FECFJxr/+9S+cccYZ2LdvH5588kl89913+Pvf/479+/dj9OjReP755yMeO5AQKioqwrJly3D++edHsXLrx1SyZcsWPProo7pC6OGHH8acOXMsmddsbrzxxrAFbk1NDR599FEKIUIM4Ij3Agghxvnpp59w5513YtKkSZgzZw4cDt+f8OWXX44pU6bgD3/4AwYPHowzzjjDtHlTU1Nx+umnmzaeVWMapXv37nGZNxI6dOiADh06xHsZAIDa2lqkpaWFbaEiJJGhRYiQJOLxxx+HIAh44YUXVCIIABwOB2bNmgVBEPDEE0/IxyXXyrp163DxxRcjJycHubm5uPrqq3HkyBH5vC5dumDz5s348ccfIQgCBEGQ3Ud6bixp3A0bNmDq1KnIzc1F69atcffdd8PlcmHbtm0477zzkJ2djS5duuDJJ59UrVdvTGlevX+SZWf16tW4/PLL0aVLF6Snp6NLly644oorsGfPHnmcN998E1OnTgUAjBs3Th5DmkvPNVZXV4cHHngAXbt2RUpKCtq3b48ZM2bgxIkTqvO6dOmCCy64AHPnzsWQIUOQnp6OU045Ba+//nqoX19E6LnGfvjhB4wdOxZt2rRBeno6OnXqhEsuuQQ1NTXYvXs3CgoKAACPPvqo/NqvvfZa+flLlizBOeecg+zsbGRkZGDUqFH4z3/+o5rjzTffhCAI+Pbbb3H99dejoKAAGRkZWLJkCQRBwAcffOC31rfffhuCIGDVqlXmbwQhFkGLECFJgtvtxoIFCzBs2LCAFoKOHTti6NCh+OGHH+B2u2G32+XHpkyZgmnTpuGWW27B5s2b8fDDD2PLli1YsWIFnE4n5syZg0svvRS5ubmYNWsWAK/VJhTTpk3D1VdfjZtvvhnz58/Hk08+icbGRnz33Xe47bbbcO+99+L999/Hn/70J/To0QMXX3xxwLG0LqDa2lpcc801cLvdaN26NQCvgOrduzcuv/xytG7dGmVlZXjhhRcwfPhwbNmyBfn5+Tj//PPxt7/9DQ8++CBKSkowZMgQAIEtQaIo4qKLLsL333+PBx54AGPGjMGGDRvwyCOPYNmyZVi2bJlqL/773//innvuwf3334927drh1VdfxQ033IAePXrgzDPPDLlnoijC5XL5Hfd4PCGfu3v3bpx//vkYM2YMXn/9deTl5WH//v2YO3cuGhoaUFRUhLlz5+K8887DDTfcgBtvvBEAZHH0448/Yvz48Rg4cCBee+01pKamYtasWZg8eTI++OADXHbZZar5rr/+epx//vl45513UF1djVGjRmHw4MEoKSnBFVdcoTr3+eefx/DhwzF8+PCQr4OQhEEkhCQFBw8eFAGIl19+edDzLrvsMhGAeOjQIVEURfGRRx4RAYh33XWX6rz33ntPBCC+++678rF+/fqJZ511lt+Yu3btEgGIb7zxhnxMGvcf//iH6txTTz1VBCDOnj1bPtbY2CgWFBSIF198cdAxlbhcLvHCCy8Us7KyxDVr1gR8vS6XS6yqqhIzMzPF5557Tj7+ySefiADEBQsW+D1n+vTpYufOneX7c+fOFQGITz75pOq8jz76SAQgvvzyy/Kxzp07i2lpaeKePXvkY7W1tWLr1q3Fm2++OeA6JQCE/Kdcs7TPEp9++qkIQFy/fn3AOY4cOSICEB955BG/x04//XSxbdu2YmVlpXzM5XKJ/fv3Fzt06CB6PB5RFEXxjTfeEAGIv/vd7/zGkB5bt26dfGzlypUiAPGtt94KuQeEJBJ0jRHSzBBFEQD83ClXXXWV6v60adPgcDiwYMGCqOa74IILVPf79OkDQRAwceJE+ZjD4UCPHj1U7qtQ3H777fjPf/6DTz75RLboAEBVVZVsXXI4HHA4HMjKykJ1dTW2bt0a0Wv44YcfAEDlPgKAqVOnIjMzE99//73q+KmnnopOnTrJ99PS0tCrVy/Dr2/atGlYtWqV37//+7//C/ncU089FSkpKfif//kfvPXWW9i5c6ehOQGguroaK1aswKWXXoqsrCz5uN1uxzXXXIN9+/Zh27ZtqudccsklfuNcccUVaNu2LUpKSuRj//rXv1BQUOBnUSIk0aEQIiRJyM/PR0ZGBnbt2hX0vN27dyMjI0N2JUkUFhaq7jscDrRp0wbl5eVRrUs7T0pKCjIyMpCWluZ3vK6uztCYf/3rX/Hiiy/ipZdewnnnnad67Morr8Tzzz+PG2+8EfPmzcPKlSuxatUqFBQUoLa2NqLXUF5eDofDIbuPJARBQGFhod8etWnTxm+M1NRUw/MXFBRg2LBhfv+6desW8rndu3fHd999h7Zt22LGjBno3r07unfvjueeey7kc48fPw5RFFFUVOT3WHFxMQD4vVa9c1NTU3HzzTfj/fffx4kTJ3DkyBF8/PHHuPHGGw25UwlJJCiECEkS7HY7xo0bh9WrV2Pfvn265+zbtw9r1qzB2WefrYoPAoCDBw+q7rtcLpSXl+te1OPJm2++iYcffhgzZ87E9ddfr3rs5MmT+Pe//4377rsP999/P8455xwMHz4cAwYMwLFjxyKes02bNnC5XKrgccBrXTt48CDy8/MjHtsKxowZg6+++gonT57E8uXLMXLkSNx555348MMPgz6vVatWsNlsKCsr83vswIEDAOD3WgNliN16661obGzE66+/jldeeQUulwu33HJLhK+IkPhBIURIEvHAAw9AFEXcdtttcLvdqsfcbjduvfVWiKKIBx54wO+57733nur+xx9/DJfLhbFjx8rHwrFqWMHcuXNx00034frrr8cjjzzi97ggCBBF0c/q8Oqrr/rth3SOkddzzjnnAADeffdd1fHPPvsM1dXV8uOJht1ux4gRI2QX1dq1awEEfu2ZmZkYMWIEZs+erXrM4/Hg3XffRYcOHdCrVy9DcxcVFWHq1KmYNWsWXnzxRUyePFnlLiQkWWDWGCFJxBlnnIFnn30Wd955J0aPHo3bb78dnTp1QmlpKUpKSrBixQo8++yzGDVqlN9zZ8+eDYfDgfHjx8tZY4MGDcK0adPkcwYMGIAPP/wQH330Ebp164a0tDQMGDAgJq9t165dmDp1Krp164brrrsOy5cvVz0+ePBg5OTk4Mwzz8RTTz2F/Px8dOnSBT/++CNee+015OXlqc7v378/AODll19GdnY20tLS0LVrV10L2Pjx43HuuefiT3/6EyoqKnDGGWfIWWODBw/GNddcY9nrDpcXX3wRP/zwA84//3x06tQJdXV1cur+b37zGwBAdnY2OnfujC+++ALnnHMOWrduLe/X448/jvHjx2PcuHG49957kZKSglmzZmHTpk344IMPwqoR9Ic//AEjRowAALzxxhvmv1hCYkE8I7UJIZGxbNky8dJLLxXbtWsnOhwOsW3btuLFF18sLl261O9cKetozZo14uTJk8WsrCwxOztbvOKKK+TMMondu3eLEyZMELOzs0UAcmZVsKyxI0eOqMaYPn26mJmZ6beOs846S+zXr598XzvmggULgmZS7dq1SxRFUdy3b594ySWXiK1atRKzs7PF8847T9y0aZPYuXNncfr06ao5n332WbFr166i3W5XzaXNGhNFb+bXn/70J7Fz586i0+kUi4qKxFtvvVU8fvy46rzOnTuL559/vu7r08u40wJAnDFjhu5jeplu2qyxZcuWiVOmTBE7d+4spqamim3atBHPOuss8csvv1SN9d1334mDBw8WU1NTRQCqvVm8eLF49tlni5mZmWJ6erp4+umni1999ZXq+VJm2KpVq4K+ni5duoh9+vQJ+boJSVQEUWxKMSGENEtmzpyJRx99FEeOHEm4WBeS3GzYsAGDBg1CSUkJbrvttngvh5CIoGuMEEJIWOzYsQN79uzBgw8+iKKiIr+yA4QkEwyWJoQQEhaPPfYYxo8fj6qqKnzyySfIyMiI95IIiRi6xgghhBDSYqFFiBBCCCEtFgohQgghhLRYKIQIIYQQ0mJh1lgQPB4PDhw4gOzs7LCKjBFCCCEkfoiiiMrKShQXF8NmC27zoRAKwoEDB9CxY8d4L4MQQgghEbB371506NAh6DkUQkHIzs4G4N3InJwcU8dubGzEt99+iwkTJsDpdJo6dnOFexYe3K/w4Z6FD/csPLhf4RPJnlVUVKBjx47ydTwYFEJBkNxhOTk5lgihjIwM5OTk8I/BINyz8OB+hQ/3LHy4Z+HB/QqfaPbMSFgLg6UJIYQQ0mKhECKEEEJIi4VCiBBCCCEtFgohQgghhLRYKIR0KCkpQd++fTF8+PB4L4UQQgghFkIhpMOMGTOwZcsWrFq1Kt5LIYQQQoiFUAgRQgghpMVCIUQIIYSQFguFECGEEEJaLBRChBBCCGmxUAgRQgghpMVCIUQIIYSQFguFECGEENJCqW1wQxTFeC8jrrD7PCGEkBbP/hO1WPrrUVwwsBjpKfag5z773XasLT2Bkd3a4LLhHdE6M8Xy9ZWdrMWBE7Uoyk1HUW6aX1d1t0fEc9//gjd+2oUGlwdnn9IWL1w91G+cgyfr8Nz3v2B0j3zUNbrxwJyNGNqpFWZdNQStMlNQ0+DCseoGtM9LN9S5PRx+OVSJzQcqcP7AIjjtiWOHoRAihJBmQL3LjU/X7MPZp7RFUW66oedsOVCBAydqUZyXjj5F2SEvfJv2n8SWAxUY0CEXNkFAbroTbTK8okEURXy9sQxVdS50b5uFIZ3yIIrA3uM1aJ+XjmPVDVhbegJjexcgzWmHKIp4dfEu7D9Ri/snnoI0Z3DxEYzaBjd++vUoTu/eBlmpDjS6Pdh3vBZtslKQk+YEADS6PZizdj827D8BALh2VBf0aJsNANiw7wSue2MVyqsb8P7KUjwwsQ/qGt0Y0a01Uh3qda3ZcwzPfvcLAGDR9iP4YGUp3rtxBDq2zlCd53J7UOfyICtV/zJb3+iG3e6AzSbI+7fraDU+X38AFbWNGNu7ALnpTizfeQyfrNmLnUeq5eeO612AF64eilSHDbvLa7Bo+xF8+d8DWLPnuHzON5sOoq7RrdrX6noXrntzFbaWVeCDlaXy8WU7yzHxucXo0TYLa0uPo6bBjcGd8jChbyHat0rHuf3a4WRNIz5evRebD1SgXU4abhvXHW2z0/xe1/ZDlThSWY8TNY3Ye7wG43q3RZf8DDz27y14f0UpPCKwtvQ4Tu2Yh49X70VhThpOKcrB9JFdQgpQq6AQ0qGkpAQlJSVwu93xXgohxABLdxzFutITyE5z4JrTO/td0KvrXcgMcEFqLtz27lp8//NhXHhqMZ67fHDQc0XRaz2QLugAMHNyX1x7RteAz5m76SB+/8E6NLg98jGbAHx752gAwDsr9uKx//wsP/Z/lwzA9kNVeG3JLmSnOlDT6IbbI2Ji/0LMumoIPl69F//79VYAQGWdC3+fOjCkEDtR04DVu49jbO8C1Lk8WLGzHB1aZeBPn23A+r0n0KFVOoZ1boWvNx1Eg8uDVIcNL1w9BH2LcvH7D9Zh5e5j8ljvryhFUW46qupdOFnbKB9fV3oC015aBgDo1S4L/++Cfjita2ukOGwQRRFPfON9jaO6t8He4zUoPVaDy19eju/vOQs2QcC8zQfxxfr9WL7zGKrqXTijRxu4PSKOVjXg1d8NQ0NjI17aasO9K39AbroTd5zdE+tKj2PJr+U4WlUvr+PNpbtVr90mAIU5aThUWY8F245gyqylqKxrxL7jtfI5aU4bZk7uh/tnbwQANLg9+NNnG1Be1YC3rz8ND83ZiK1lFchJc6CmwQ2XR8SlQztg2Y5y7D9Ri4MVdap9WFd6AgDQo20WyqvqcbzGt0+frN6LoV1aQxRFlJ2sw4OTTkGqw46rXl2hWvfLi3bijB75+Oq/B+Rjby/bg7eX7ZHvp2w6iBtHB37vWU3z/mSIkBkzZmDGjBmoqKhAbm5uvJdDSIvmRE0D8jJS4PaI2HW0Ct3ys+Rv0QDw2pJdeOzfW+T7Qzq1Qv/23r/bepcbM95bh4XbDuOzW0dhUMc8Q3O6PSKcYa5z2Y5yvLZkJ2aM64HBnVoBAH49XAmPCPRql+13/oZ9J3DT26sBeB9/cFIfHDhRi//uPYHehTnoXZiF4rx0ZKSE/pg+eLIO3/98GADwxfoD+N8pA/D/vtiEgqxUXDq0A3q2y8ai7Ufw0eq9aHR5sKWsQr6AdivIxM4j1Xhq3jaM7J6PBdsO470Ve3CiuhF9inPQrzgHB0/WYd7mg/CIQNf8TByprEd1gwseEdh5pBq7KoHnV2xrei1Z2H6oCo9/8zMqmgRGZb0LACAIXkvFXR+txzebDsrr/2ztPqQ6bbj+jK44WduIgyfrkJFix7hT2uL9FaX4aFUpXr92OP7+7XZ8sLIUs64agi0HKvD8gl9V+7DveK38uhw2AfUuD256ew1EUYRHBLJTHbhyRCfsPFqN+VsOYf8Jn4j4TZ+2mDGuBx7+YhMOnqxHo9uD7YeqcPVrK5DisKFtdqo8R6rDhqennQoAOPfZRdh/ohabD5zEW0v34EvFBR8Afvq1XL69YNth/HqoEltO2AB4xdEjX26WH7fbBIzpmY+i3DQs2n4UggAU56Xj0qEdcF7/QuSkObFq9zFMf30ltpZVAACcdgHDOrfGmF75uGBAMTq0SpeFUH2jB1+s967nUGUdvt/qfY/Mumoo2uWkYt+JWoztVYCqehdW7DyGYzUN6F6QhY6t0/Hxqr3YebQai7Yfwa+HqwAAfYpycOGpxfhm00H8d+8JLNp+RF775+sOYGhn7/s+O82BHm2zcKy6AXvKa2QR9OLVQ7Bpv+/3dv0ZXdEmKwWVdS444ugqoxAihETFsh3l+GTNXtxxdk90zc/UPafR7cH7K0pxVq8CdAlwjpbqehf+9NkG/HtDGQZ1zENVXSN2HKnG/RNPwS1ndQcA7DhShSfn/qx6XlXTRXfvsRo8/MUmLNzm/bDefKBCFkKiKOLlRTvxwo874LAJ6NImE0M7t8K2gxVYt8uOe1Z8h9+eWozHLuwvW5IaXB6s3HUM+dkp6N1O7UaqbXDjro/W42BFHZbvPIY3rxuOrDQHznt2MbLTHFjz5/FIcdhU59/54XocqvBaAA5V1GPic4t196F1ZgrG9i7AJUM64LF/b4FNEHDr2O6Y2L8QDrsNbo+Ip+Ztk8/vV5yDuZsOYvba/QCAVxbvxEPn98Uz87fLewMAqQ4bHv1tP0wb1hGXvLgU60pP4NxnF6nmXrnrGFbu8llRrjitEx67sB8cdhumvrgUq3YfR12jG9/us8HlEXH+wCI8M+1UTHjmR+wurwEATOjbDr8/pycyUuz44efD+Ot/tuLzpovzuN4FGN2zQHabvL+iVDX/69cOw4NzvBf1T9bsw5FKr8Xi4Mk67DteI5+Xn5WKl64Zgq/+W4bKOheuPr0T+hXn4u6P1+PfG8oAAIM75eHvUwehe0EWAGDX0WpU1DYiPcWO4rx02YX17zvGAACOVTfgqXk/49vNh1Be3aCyvNw6tjsKc71uoYLsVJysbUS9y4NtByub9qkjrhrRGdlpDnyz6aBsRapr9Mi/g9vO6ganw44PV5ViTM8CTBvWEQPa54Z0Dw3v0hof3HQ65m0+iKGdW+H0bm38rJ1Ou4BGt4hqxe+7tsGNmkavl6NH2ywU5qahZ5NAz05z4jd926nGuP3snvI+PPfddrTKTMGtY7sj1WHH/4zphvX7TuDnskp8tKoU/913EtX1Lvm1TexfiCcvHYQDJ2rx2+d/wtGqevz+nJ44r38RxvctRF6GEx1bZ+DcfoVBX2usoBAipIVxoqYBm/ZXYHjXVn7xD0o27T+JrvmZQV1K/9lQhjs/WodGtze+Yfato3TdG/d9ugFz1u3HoI55+GLGGX6PL9p+xOsqmTIAuelOVNW7cOkLS/Fz04Xlv3tPyOeu3n0cOMtrtbn3k/+i3uXBmJ75OHiyDr8croJHFPHKop14/Jut8CiSYdwer0vH4xHxt6+34tUlu+THjlY1YLUcXyEAEDF77X4s/uUo+hblwGkXsH7vSdl10b0gE09eOhBDO7cG4HVjSG6FqnoXrn5tBfKzvBaEyjoXahvdeG/FHpSdrMMDE0/Bk/N+xs6j1SjMScM/rxiMlxftxHdbDyHNacP4voXYdbQKe47WoLLeG7g6e+1+WdgAwB0frEObzBQM7tQKe8qr8UvTN3YAOFHTqHKxeETIFrNBHfNwyZD26NQ6A8O6tJYv/n+9qD8ufP4nuDwi+hbl4NozuqB/cS62llVg04GTcNptuHhIe5xSmCOPK8We1Ls8qGr0/s4vHtweKQ6bV6y+uxYpdhseOr8POrfxit+u+Zk4XFmP0vIanNW7ABcPaY9Uhx2922XjqW+3YcuBk3Lcyf4TtXj4c5+1pHVmClxNv9A6lxt1jd7f5z3je+G60V2RleqQfx8Sz10+GFOHdUTXNpno1EYdwxNItCvne/zigfjfi0TsPV4j72nb7DR0aOWLwUppsmQ0uDyy23DK4A6yVfKWs7pj77EavLeiFPUuN+pd3nMKslNw3ejuuGt8r6Dr0GNQx7yg1k2n3YZGt1slfCvrXHA37V96GPFYrTNT8OiF/VXHbDYBQzq1wpBOrZCV5sDvP1iH6gYXahq880lWzOK8dHx1xxnYvL8CZ5/SFoDX6nXjmG6G548FFEKEtBBEUcTDX2zChyv3wuURcWrHPLw6fZh8wVYyd1MZbnl3Lc45pS1eu3a47njvrdiDP3++CVLm7brSE/jb11txoqYRZSfr0Lc4Bw9O6gNRFDFnnfcirhQ0EtX1Lvzu9ZUAgNO6tsbvRnbB//5nK34+WIn8rFQ8eekA7DhcjRW7yvHd1sPYedR70X918U5vXFCqA/93yUBc/+aqptcJfL2pDB4RGN6lFQ5V1KP0WA0a3CJcbg/+9NlGfLZ2HwDgwUmnYHSPAqzfewIb9p1A1zbpqN23BYOHjcD9czaj7GQdfqz0mf/zs1JQUefCjiPVuOyl5Xj84gGY0K8QLyz0mvr/d0p/fLflEBZsO6KyIHg8Ip6cuw21jW5cO6oLvt7otVI8eqE3/mR4l1ZYW3ocHVtnqAJQT9Y2YsuBCjzy5SZsP1SFc05piwEdcvHOsj0or27Ad1sPAQBy0hy44rROeGnRTpyoacDx6gYA3oDgPeXVWLDtCDJS7Pjn5afKokRJv+JcfHH7GWhweXBqxzxZzPYtzsElQzvo/v4lIVTb6EZDU9iQdIE9t18h/u+SASjMTVfNJwgCHpzUx2+s0T3zMbpnPkRRhCAI2HKgApP+uVjluoIIuNzeN1t9owf1Lq91o11uWsCAZLtNwFm9CnQfM4rNJqBzm0zdfQMgW/oaXB40NIkcpfUPgPyFo97lQV2TVSbYl5BoSXHYUNPgVlmEjtU0yLfTUsxzQ2Wlel9Hdb0b1fXe15aZ6ntt3iw3Y8H78YJCiJBmhnQx0bJ85zG8u9zrevBaOE5g4nOL8bvTO+PMXgXoXZgtX9ykIFYp7kTLx6v34qE5mwAAV47ohKKcNPxj/na8sthnZVny61HcNKYbDlf6AjBP1fkWqwwKtdsEORMHAP51xWCM7N4GZ58CnD+wCN9t/QGl5TXYfqgS/5i/HQDw8AV9UaxI9fU0xYMA3m/j/9lYhtJjNWh0e/Dy4p34bO0+2G0CnrxkoHyR71ucgytHdEJjYyO+/noLRnVvgx/uGYv1e09gT7k3W6dtTirG9CxAbaMbD87eiH9vKMP/+2Kz18VS50LPtlm4fHgnXD68E56c+zM+WFmKijrvhcgtirK1oKbBhZqmC4YUOyQIgp81AwBy050Y2b0N/n3HGOw8WiW75GaM64FlO8qx93gNUh12nH1KW9gE4KVFO1Hd4MahJutUQXYq/nhub7yyeCdGdc8PeDEHvGIoHCTRU9foQZNxBmlNbh1BEHDZ8E5hjSc9D/D+PvoV52DzgQr5MZdHhKvJqlfv8siWlVRHfNOwZSHk9q0pRRPvkub03q9rdMtiycp1S6npSouQJI5tgv/6okGy/lQ3uGThlWyJCcm1WkLigCiKOFHTiFZBaoWIoohfD1fheE0jOrXOkOMHYs0z87ejZMGvOL1bG9w8povqsQ9XecXFFad1xE1juuHGt1Zj59Fq/GP+dvxj/na0y0nFgnvHIiPFgb3HanVG9+LxiHj+B68F5KYxXfHgpD6od3nw4/YjOFRZh0n9i/DSop0AvO4rpUsnO80BURSxbEc5ehdmI9Vpx8tN50p83mQ9unJEJ4zs3kY+XpiThnSnHbWNbjz73XY0uDwY3SMfU4d5xYwUP+0RIReIswmC/KHf6PJg2yGvq+33Z/cMaOmQSE+xY2T3Nqo1AN6LzL+uGIzd5dXYtL8CsxbuAADcMLor7E2LeGBSH9w/8RT0eOgbuD2i/A/wCoe6JmuGdIEMRYrDpnJLOe02nKmxdHg8IgTBaxHbddQr3lpnpiAz1YE7fxO++yUUqou7JIRMtHJcOrQDNh/wBcG7PB7ZIlTX6JYtK9Gk3ZtBqsoi1GTtcQa2CEliyejvPhJSdITQsSYhlO60m1ofSLLGVde7UNPQZBEyEOCfSCRORSNCEoiqehc27jsJt0fEHR+sw+DH5uOJb36Gx6NfgfXDVXsx/plFmPbSMoz9+wJVOu7q3cewevcxlFfVW17B9dM1++DyiFjy61E89IUvvuJ4dQO+2ejN0rnytM7oVpCFb+4cg2cvOxVn9iqATfAG7O47Xit/cwSAjq29Ju3FvxyR3RQrdh1D6bEaZKU6cNf4XhAEAWlOOz69dRQW33c2HpjUR/4g9ogivt3iyw4Sm2qIXPnqCjwweyM27T+p2iu3R5QvFKcUqjOtbDZBDrSev8XrEpo8qEj+UBdkISTCI+2z4Pt23Oj2oLHJKpOfHV0BPEEQcOtZPeT7rTKcuGhwe79z7E2LalSknFfXu9DYdEE30z1is3nr+gCQa860yrCu0F+ajkXIzDowlw/vhKtP7yQLDZdblGOEEsoipBMjpLW4SMKoXiGCrVy3ZKWSXFUAcLzJNWZ2rR7J+lNT74tJokWIkCRn2Y5y3PPxehw4WYeOrdNl68iLP+5AXaMbM3/bT3V+XaMbzzS5abz3PThSWY/cdCeW/HIUV7/mq6uRk+bA1ad3xh/P7W161dayk7WqmArp2xngTU9ucHvQv30OBnTwukBSHXZcNLg9LhrcHkMfm4/y6gaIIrB8py/dt01mKrYdrMQ1r63EiK6t8dHNI/HJmr0AgAsGFgVO7W56aSKAqjrft1KPKOJIZVOmVGW97CaQ8F7svMfsNv/96VaQia1lFbKQULqTbIKvMF3TELAJgk8IeUT5eWZUtT2vfyG65mdi19FqXDmik65lwmYD4IY8LwCV8DPbKtAqIwUnahrldPU2WdYJIZ9rzD9GyJTxU+z460UDUF3vxpx1+70WIck11uhWWFbiaxFSusYCub3Smu7Xu9yob5TOsW7dTrv3b0EVI1Ttfd+ZvV+ZTcJK5RqLU2HESKFFiLRoPl2zD//z9mrc+u4aLNp+BOtKj+OqV5fjwElvjIUkgi5u+rb/6Zp9fmO8u3wPDlfWo31euuIDwHvh+6EpxsZrjgYq6lyYtXAHHvlys+wqkRBFETuPVAW0OoVi9W5v1pP0IexRWJ8k91SguA1JlIkQsXSHTwiJoojypmyZo1X1qK53yZYlySWlh6RhRFGE8tUo43dcbt+FTULpQnLoCKHuikyfVhlOdC9QB+J65/S9dpsAOB1NVhmXzyJkRoyE3SbgX1cMxi1ndcetY3vonuOw+awFEiohZPLFMC9DXf3ISotQatMFtareBY/o3WMzhZCE9D5weURfsLQq6DgxYoRqG9zye9svWFphPauTxJKVrrGm+St1YoTM/h1J1h+PCJQ3zZFBixAhicO2g5V46ccduG1cD/Ro660f8uvhSqwtPYHN+0/iLUV1059+PYpuBVnwiN7iaveddwreXb4H/YtzcUbPfMxet1/l4gC8F7gXf/TGiPz+nB544pufUa34QJSsK09eOhDj+7bD7LX78dDnG/H2sj1Ys+c4BnbIw7Hqevzx3N54YaE3kPfUjnm4e3wvDOvSylAxPYnVTVVzT+vaGot/OSpbRX45XIUtZRVw2ARcMKBI97myW8kDRRq5V85JQkYUgd3l1ahtdKNVhhNDmooG6o4HhShRCDuPQqS43KLKUgJIAbHeY3ab/4WiW1MNGMBrDVJa1ZQxQr5jihghRTCrWX2O+rfPldOk9ZDWpCeEUuw2VWFIM8hLVwshK3tgSRdUZbVhM7ORJBxN1g21a8ytcI3F2SKkE4/jnzXmswjFMli6RidrzGzXmFJYHW4K0s9KTS6LEIUQabacrGnE9W+uwv4TtThW04A3rzsN/95wAHd/9F9Vm4CbxnTFkl/LsbWsAuv3nkCK3YZHL+yP9nnp+EtT/YyDTRYijybG54efD+NoVQPaZqfi4iEd8ORcb2E7UfTOv/WgN+tlRLfWSHPaceWITshIsePhLzZh84EKOStm0fajqG36hrt+7wn87vWVSHHY8MTFA3DxkOBBvRKSgBnRJITcTWv96r/eVO2xvdsGDPiWLsciRLkWiPQ6pJfsUbic0kIEXAqyRQgqi5CosAg1KoJfJdweT1CLUDeFBWh4F7UQs6myxkR5HdJFocEt+ixCMbIiSO495fvtRJNwsMIioLQA2QTIMUNWILn1TjQJO7OzkSQkq5pLG3DeGF7AuVVI7yWVENLGCKmCpWOQPt80f3WDv0XIbNeYzSYgM8WO6gY3qpvc8eF8gUsEkmu1MYK9xpIfURRx76f/lWNmFm47gn99/wue/m47RNFbfTe7qfbKhae2x6b9J/Hb55fAI3qzldrnqeteSMYJrTtLcpVNGdIeTrvNJwAgYuXuYxBF78VbWRvmosHtMaZnPt5etgeNbg8W/XIEm/Z7BdHvz+6Bo9UN+GHrYRysqMMDszeib3GOKmNIwuX24KPVezGiaxsU5qbJJfeHd/HGzXg8XtHxZVNl3YsGFwfcL5vCraTUeiJ8osKtEBi2EPFNNoWrTamElBldyngg+TUp3B96MULKInjD/IRQ05oVYksVI6QIlpZiKKzGHsQ1ZkVsS67CNZaXkaK7h2Yhrf9Eje8Ca3bcG+B7H7gUvz+VRShBYoSkWDibAL92EWlysLSvEKSlWWOSa0wRnycHS1uwXxmpDlkEAcmXNZZcq40R7DWWHKza7S3/f90ZXfy+gXyyeh/mbzmEFLsN3dtmYWtZhVx35soRnfDYhf1VF4n+7XPx/y7oiwXbjuCOs/3jPeyytcF7oX1q3jYcq27Awm3eGKBLm6w2ci0bj88tdnq3Nn7jtclKlSvK3jK2Ox77aguK8tJx1296QhAEeC4Ucf1bq7Bw2xFc9tJy9CnKxv0T+yAr1YHr31yF28f1QIPbgz9/7q3l8/HNI+ERgeLcNBQ3iTiPKOJkA7D/RB0cNgG/6dPObx0SSguO0uqltOh4PIrYmxCf4dLOejTjKa01LreeRSh4jFB2mhO/G9kZhyvqMbBDnmZO3+9IKdgk0dPo9qDR5T1uheVCD2kalUWoSQhZ4RpRWoRaZVhnDQL8XWNWXdil35+2BIEkLtMSJEZIsgjpWRsl609dY2yy3WSLkDJGqOn3ZIUQykp1yEkQgLqgYjJAIUSSkg9XluKhzzfB7RGxp7waT146SH6s7GSt3FLg7gm90KcoB9ObKheP613gJ4Ikrj2ja8Du28rzT9Y2yrVjAGBQh1y5Z4/SxbSqKWZHTwgpyUlz4qmpg1THbDYB/5g6CFNmLUXpsRos33kM76/Yg8GdWqH0WA2+3XIQBdm+itBSjZDC3DQ57sQtev8B3g/dYBYI5brFABYcMQyLEBTWGXWwNGT3WqMngEUoSNYYANld6TelIn1ezp7XpM9LgsQZK9eYlD4fI4uQMljayvggQOEas/ACC/isao2KmLKKOl9cUrwtQqmaGCE9kS25QdW93qzMGvNPn5fba1iQ0ZWhGZPp84RYzLKd5XJ3ZQD4ePU+tMtJQ/u8dJzVuwC3vLsWlfUunNoxDzeN6QYB3uDnEzWNePbywRG5C5RBrZJpGwDa5aTiD7/pKd9XWlYkcaLsSxQObbJS8e1dZ+KJb37Gm0t3w+X2CRFRVJufK5suDFlpTp/1yqNwEYV4zYLG4iWhFDJujcspGD7XmNrVBo1FSBssrbIIhem+sgXKGrP7LqQNASr/WoVNJ0bopOxKMn8NeQqLkPVCyHvxkwKYrUpj91mEPHK/OGXmXaJkjUluqBQdgSNlB8Zq3XpxS/JaLPg9KYWP3SbE/XcSLhRCJKkQReCZ77xVjacN64DcdCdeWbwL/2qqdCxV1m2V4cQzl50qi55Xp+v3yzKKXXHhly6mggCsePA3qvP0Ym1CWk+CkOa0y/FKSlHhEUXVh4/0gZed5vDFM4UhXIQAFhx5YkjWHJ/AMDqe2jXmew2upt5fSpSZQXpZY8GQTlcHS2tcY25zs8ZCIbn3lDFCkmvM7NR5QO0Oi5UQkrDOItT0+1PEj1U07aHdJsTsdxkIrejQEwFai5AA0dI4NadOsLSEFb8nZd2gjBRrYsWshEKIJDQutwevLtmFhdsOo7bBDaHWhvXlJ5HqsOHeCb2Rk+5EvcuDwxX12Hm0CtsPVSHNacNr1w4P2V06HJRWJOnbvV3nj13pYpKINl5V7fIRm26r/fAVtd4PvJw0h9oyorN+PWwqi5DvuNK6IoqinIkWUlgpnu9fR8h7pNHjkUWPRKisMSOvQStCpQuVUgilOGLzQa1nETpRY6FrLF0ZIxRbIWRdjJB3XLdCJEtvm0SwPMjp83XBYoTUx5w2WCoWfJWldYSQBSUOlF/Kki1QGqAQIgnE1rIKrNhZjoo6Fy4d2gFOuw03vb0a61Udy71/xFeN6Iy2Od5MLClmxOPxFgNsm5MqN7M0C5uORUhPXChdTLJVAtF94KkLBaJpfLVF6FBTY9NshWtMWod2/fpzSLfUFhxlzJBHNG7lUrrGtFlj0prcnlB1hMLbN2XTVd86Fa4xl881FisrgvS7UFqEpPgWa1xjsbMIaS0LVrnGHLJFyL8AZ7yrSgM+V1jQGCGN9c9q/ZbSZG1SZo1JWGMRUgihJAuUBiiESIJQXlWPS15YKreF+HpjGdrlpGH93hPITnPg7vG90Crdgbe/X4+MVgWYMa673xg2m4DRPfMtWZ/yoixZFfSFkPenqAnYjQabjkVIFNUfuAeaygRkpzpU8UDS9TfUdT+wBcd33+1RBEuHjDlSjqeNOWqyCOm4xkJljRl5Df5ZY74WCJLwinkdIYUQEmVrRvMIlpaw2jWmV4AzESxCkosrWNaYdq+cFhskU+QCjh6/x6yOEUq2QGmAQojEiePVDZi/9RAm9G2HvIwUvLV0N2oa3OjQKh01DW78fLASPx+sRIrDhk9vGYXehdlobGyEsNeDSZOGwum0NjVYi/Ka3BjMNSbXEVJYhKJ1jUE5pve2RxPLIwuhNIdKoEnXDb21KlG7xjTCRXbH+URKaI2iLG7oO6q936CNEfKIQYVm8NegnMO399KFStkUM2YWIcma4fa/IFlRUDEr1QGHTYDLIwYsnmkW/hYha11j2r50QGIIIa3wCZY+L2F1Dchg728rssaUViBtBlkyQCFEYka9y40TNY1olZGC6W+sxIZ9J/FkVgruOLun3OrioUl9kOa047o3VwEA/jihN3oXmuvmigRBEGATvEJBupjqWUWUjT+l6320rjFpHuWYSusQAOw/Lgkhp0qkuJt6QIVs5aBjyQL8W2wYdY2pLEJ+wdK++7UN6qKlaotQmMHSChUqiS0BvhYbUuVuIPZCSO8ibsU3c0EQ0CYrBYcq6lGQlRr6CVGgTVu32iIkdW1XkgiuMa0Y03eNqY9Zrd+Cvb+tECpKK1AWLUKE+ONye/D8gl/xzrI9KK9uQPu8dLni89GqBjzy5WYAQLf8TEzoVwi7TcD/XTIAZSfrcP1o/bo+8cBuE+BRmOd1XWNNP5UCIMzruaExlbE2AOSqrtmKYGlAYREyGCytFD7SAWnOcFxjcpVn6DRdVSxcEidOu+B1lZkWI+Tbe4dOcbnYFVT0D5aWsCJrDAAemNgHG/efRN8i/2rkZqK1AFmdPl/X6C+EksUiZLN5BblskbQ6RijIvlgSI6QQP8nWXgOgECIms/nASczfcgh7j9Wifat0DO6Uhzlr9+PL/x6Qz5FE0KyrhmD/8Vp8u+Ugdh6pxoOT+sgXjkBd0uOJVyyIcnG8UMHSosIqEQ16QcDatHSJ7DSn2jXmUa49yBxNPz2iNlhanbJv1DWmbLrqn4Xmu1/b6Ou71Oh2qbPGwq4jJK1TbbmSLqRK61PMYoSEIELIoqvhRYPb46LB7S0ZW0mK3SZbSQHrCyoq63dJxLuYIgCk2NVrCPTeSnXEUAgFEfqWxAgprEyMESItil1Hq7FiZznqGt04r38RVuwqx50frYfO9RkOm4DHLx6AEV3b4O1lu9GtIAuTmjqh33RmtxivPDK03+6Dxggp3FjRps/rFQrUprlLaC1CLkX2VDDkp/gJF1HlGgu/LpH6uFZoKRtnVtWrawuFHyOkjEvyCTZtTRWbEP7YkWKLsWsslgiCgDSnXU5wsCLmCUg+i1CgNaU6bZC6UDhsOn+8Fq5JidUWoUzGCJGWQmVdIy55YalcPXnWwh2obXRDFIExPfMxvEtr7D1Wg8W/HMXxmgY8e9mpmNgkfP58Qd94Lj1itIGvQV1jMDFYWmnpkMbXxAhJ5GhihCThEk4doUAFEN2KOY0GX7s1a1RatQCflUYKJo0ma0xyQSqbrgqKOkLSsVgW4NNLn5dIhIt4tKQrhJDVMUJ6GVBWtqkwihHXGKBeqwWlfFTEPFha4Q7LoEWINEfWlR7HL4erUNvgRnWDC73bZWNrWQWOVTegbXYqMlMd2HW0GgAwqGMe3rh2uByXIYoiGtyehPjAihZt4Kte7I9Nx42FaIOlFVYmpUXIE8AipAzsdhu04Ehou8WLiigfj6KgolFx59YsUtQIrVqFRQjwCqeIY4QCNF21aVxssYoPAnzuPb2ssWS3CAHq15BukUXIIbvG9IKl4y8mte+nYBYhCcvrCMXcIuQbM4t1hJoHJSUlKCkpgdvt/4eXqEjfgqMx+YuiqKp2urWsAo98uRkrdx3zO1f6tv7Q+X1wZs8C3PzOGuw5Vo2npw2SRRDg/UbeHEQQ4B/vEcw1BtXFOLp55XgbIGSMUFaa909aCuz21REyGCyt1kEqd5goQhEjZNQ1pieEfPeluA99i1B4Vwu9Gk56brBYxQd55w/mGov/RTxalBd3qwsq6rvG4v/Z4mcRCiC0lWu1uo5QsPYdVtcRYrB0M2HGjBmYMWMGKioqkJubG+/lhKTR7cFvn/8JO45UoXe7bDx+8QD0b2983R+uLMW/fvgVRyrrMbpnPv4xdRBW7j6Guz5aj5oGN1LsNozo1ho5ad7aPd9sKoPLI6JbQSYuGFgMu03Ax7eMhMvtUYmg5oYU7xE0WFphlZAURbSl9AO12NBahNKddtkkLjQFdhvNGlPOoQqW9kt/N2atkV1jmkUGixECvO/lSC1Ceu49Zfq8RExdY7I71V+0NgeLULrKImSREJJjhBJTTPqlzwcJlpawetnB3K7W1BFiZWkSZ3YcqcLWsgoAwMb9J/HK4p147vLBhp67Zs9xPDhno3xR/eHnwxjx+PfyN9gzerTB36cOQlGur4P6T78exZtLd+PWsd1VF6vmLIIAHYtQsMrSMC9YWtliQ5nBpbW2SNYg5VrdBsWYshCkXhC2hHRBD6XtpMf1Y4T0XGN21fhAJL3GfHOoCypqhFCM+owBiR/fEi1KMZdmUZCs7BrTqSOUCHtoNEYoLYausaAxQlZnjdEiROLBziPe+By7TYDbI2LR9iNwe8SQ36ir6l2466P18IjA5EHFmD6yM/7w4XrsP1GLNKcNV57WGQ9MOsXvj+qMHvk4o4c1rSwSGe23ez33kLovmM8qEQ3qC7x0G36usWylELKphVCoTHRlIUiPxgKknMYV5LUrUdY+UqK1ZEnB0tIFtV5xsbOHnT7vLxhtNgFObcPLeARLxzB9PpakxyJGyO77vWqxKlMtHLQWR206vYTKNRbHGCErLJEZKWyxQeKMFKh8wcAi/PDzYRyvacSGfScwuFMrv3N3H63G8wt+xdjeBVi8/ShKj9WgfV46/npRf+SmO/Gf34/Gmj3HMaJbm6SsEGolUshKsKar6r5g3ttmZY1J4wLqzCiJ7DRf2xFpHUZjhJTCRRsjpCyJ6PYYS20P5BrTCi1tsHS9wv0Rdq+xAE1XtePEMljarnGnKkkEa0a0KMWcVQUig70PrJozHIxnjcXONRZI7Kc4bJaUjkhx2OSCkRRCJC7sOFIFAOjZNguNbg++3ngQC7cd8RNC32wswx8/3YCqehc+XbMPgPci+49pg5Cb7r2I5mWk4Jw+7WL7ApIE6du9oaarMC99XpWJ1nRMG2sDADkKi5BNYxEKHdystKYoYoQgQtnw22cNC7HogK4xbbB0kxBySBYh32ThV5ZWzuGzXGkvCrEMlg5aWToBrBnRoixoaFmwdBDhmhAWIcOuMd/+xKuOkFVxXABQkJ2K/Sdq0cbiHndWEP93EZGpqGtUtR8wimQR6laQhbG92gLwxvpI38Yr6xrx6Febcet7a1FV70K3gkz5uf9zZjec3q2NCatv/sjB0kHrCPlcTNJvMtpgaQltiw2tqyBbJ0bI09RrLJxgaW0laOU0riZVFDLmSH5+8BghSVil6rjGwu811jQn1DFCdpugev3xCJZujgUVAY1rzKLiOMEtQvG/hDlsgurLTiJYhPzddd77Vgqhf14xGP+8YjA6ts6wbA6roEUoQVj8yxH87vWVuPOcXvjDb3oafp4oinKMUNf8TLRuUuMb95/EuL8vRPu8dGzYd0LuRXXzWd1w74TeWLHzGH4+WIHfjexi+mtprvgsQoGLCvpq/ih6jZlpEQoSLJ2dqnCNNU0qVZY27BqDxjUGtZiRM7oMWpi0hhCPRz/tX3aNuSShFXnWmDJFXxKmTrsgfzEIllpsNjaNFVFJcxBCabFInw/y+0qEFhuC4M1MlN67qYHS550xFEIaMZaT7sTRqnpLMsYkhnZuhaGd/cMxkoH4y2kCAJi9dj9EEfh49V6/C1x5Vb3uN0oAOF7TiJO1jQC8QqhdTpoc71N6rAbLdpajusGNnm2z8Nr0YXhgYh847TaM7pmPG8d0i6mbINnRZgDpNh7VCdiNvtcY5DF9MUL+6fNKi5C0NKO9xgIFSystW4AiWDrE28Y3f/A6QhJSvIwUIxRufBCgFF++CaRhlFaglBjGlUjTNtc6QjFJnw/yZkuU6tzKz9GArrGY1hFSryE33fvZ0BzEtxXQIpQAeDwiFv9yBIC3Ien2Q1XoXZgNAFi9+xiueGU5Tu/WBm9dd5rfxXfXUW98UPu8dPlNfvXpnTFlcHss/uUI6l0eFOelY1jnVqa5aFoq2m/3wSxCyngeM3uNSWPqFVRUBktLawu315hSwEn3lUpI6gMW2jXmE1ZKtK43CZ9FyGu5jCSg0xcgrhRC3oNKV0FKDC1CUsNQvfT5RAj0jZa0WMQIBXONJciFPdVhQ2XT7YCusRhahLRWz5ymGFCrMvuSHQqhBGBLWQWOVjXI9xdsO4zehdnweEQ89u8taHSLWPzLUXywqhRXjeiseu4OhVtMSWaqA+f1L7J+8S2ISHuNRWkQ0q2Po1dQMVsnWNporzFVIUgFqtcBoNGwa8z7UxssrXW1SUgXNGn+cOODAIULUUcIKb8hxzZGyPtTzzWWCIG+0SL93gSIlgnMoK6xRLEIqYR26MrSsW6xISXDWOkaS2YohOJAbYMb3205hHn7BEwCsKjJGuS0C2h0i/jh58O45azu+GrDAfx330n5eY9//TPOOaUdqupdeGb+dizfWY6KOq9bTBkATazBULC0wsXkS+GO9gLRNKb8X4AYoSB1hIy2xNBamfxdY01uwZAWJv30eb1sN8A/6DUSi5CgYxESmoZVFlGMRx0h3RihZmQRctrMSwrQEkwUJ4pFyGnANRbPYGlZCCXIfiUaiSGnWxgnaxtxx4f/xTd7bSivbsDi7UcBANeO6gLAW+35ZG0jXlm8EwBw9/heGNwpD1X1Ljw572fc+NYq/GdjGcqrG+TA3WFdWsfltbQkHJoMoGB1hJTX+mgvD3oWIW3zUkBbR6jJNWa0jlAgIQT1a5FEhlkFFSW0F7RIYoR8tYs8fsecttAXKiuQXGPaGKEUu00/xizJkFyaVnZTTzqLkIH0eWeM0+el9kiJIhwTDVqE4kBhbhpOKczGzwcr8e8NZVi9x9vU9KoRnTF/yyHsLq/B2j3Hse2g1+s8ZXB7nNGjDS55YRlmr90PAMjPSsULVw9BqwwnGt0iTmmKKSLWIX27l4Ol9SpLQ7LEKKwSUfca8w/AVlaZdtoF9G+fi9O7+cSwHKysqLAcjIC9wTyiqqCiL1jaoGssREFFCe0HdDQWIaXxRS9YOh6uMW0doebgFgN8FgYrPS7BRHGiFKU0EiytFG2xbrEheQzat0rXO73FQyEUJ8b1ysfPByvx7Pc70OgW0acoB13yM9GvOBe7y2vw9cYyNLpFZKbY0T4vHR1bZ+CcU9ri+58PAwDuO7c3htMKFFMko4LPNaZzko4AiD5Y2vtTWU1aaR2aPrIL/nxBX9Vz/FtsRChcEMgiFGrN/oLQu279VgnaDCqzLEJy+rzCNRbLYGmfO1X9opvLN3Ola8wqghVUTJTMO6X4CWSlimWwtFTbSPpbmzqsI3q3y8apnfKsnThJSYx3UQvkrF7eXl1V9S4AwORB3sBmybIzd9NBAEDPdtnyh+l9552CdKcdQzu3wiVDO8R6yS0eba8x464xk9LnAVWwtBjE2uMTIup1BZxDDpb2j5ZW1RFyB7aG6a5ZJ2tMzyKk/WYfbp8x75qa1qiMEYqzRSiQoEuUC3i05DXFnmRY+JU6KSxC9tBCKJbp84Kmonq6045RPfJVPcGID+5KnDi1Yx4yHCJqXN6/iMkDiwEAfYpyAACVTQKpdzufy6t3YTaWPXA20px2S/rFkOD4pc/rBHHKrjGdgN1IUfbQkhAVwdJ6mkQbI2TUlaWXNaauLG3QNdb0UxsjLIr+VifA/+IRVdaYYp/0ssZiGiMUQDAmygU8WkZ0a4N7x/eEp2yrZXPo9YqTXI2JIihVrrFATVdjaBECvIUdG1wev8rqxJ/EeBe1QOw2Aafkej+wT+2YJ5cl71Ocozqvlyb2Jy8jpdmY1ZMNbbsEPaOFdP1WxQhFOa8y8FgvfV7P4iStVW6xYbgStE7WWASuMSGAa8y7Jh0hZEKMkLQNLrdSCHl/psTJIhRIMCbKBTxa7DYBN5/ZFZ0tDFHUvhcyUn3vlUQRlKmGYoRi130e8GWyxbLJcLLCHdKhpKQEffv2xfDhwy2d5+xiD3q3y8Jd43vJx4pz01TNM5UWIRJftKnQehc52cWkcs9EJ4V8lg5tsHRgUaJtuhq6W7xvXCUek11j3jH0gqXNjBHSswgpYoQSwCLUHFLnY4UgCKr3Q6bCvZMoQeeGKksr1hqLt6AkgGLZUiZZSYx3UYIxY8YMbNmyBatWrbJ0no5ZwL9vH4WzehXIxwRBwClFPqtQr8IsS9dAjCOJi4YglaX1CglGa5VWigrZIuQJXqdIjpcxmDXmc2Vpg6XVdYR83eeNjafXRFjPNWZG1phNZ+/1Y4Rid2EIFOtEq254KFPos1IVQiiJ0udjbxESmtbD91ooEuNdRFT0bRJCrTKcKMhKjfNqiIRsEXIFDpYWFNYb+ViUzjH9FhsIahGyy9YR9f1ABHaNqV+LlJFltHeZjvFHFcwsYWYdIWl8QfC9LlXBuzgUVNTSXFxjsUIZM5bZ5BpLcdgSpm1QioH3lyTanHYh6i9HhtbUtI5EEYuJDHcoARnQPhcA0K84N2H+0El4LTY8OplLkaIsdii5mpQZZHrvkUhdYzpJY7rd5432LtO1CFlWWbrpNTe9aOUITsV4zpgWVPQP9AWAdGbvhIXSItQuJw2CgIT6kqgUQoEsjvlZqXDaBRTmpMVkTU66xgzDv8YE5MJTi3GsugHjTikIfTKJGX6usSDp824zhZCixYZeQUU964xdY5EJvYYAwc2au8Hio/TW7BNrvrW7dWOEtBahSLLGvD/1ql/Hr6Ciep8mDihEbYMbV43oFLM1NAeUFsKi3HS8Pn04CrITSAjZQ1upWmWm4LNbRyHDKWDz8oXWr0kKlqZFKCQUQgmIw27DTWd2i/cyiAbpi1VjkIBhvVT3aHuN6Tdd9VmH9IOlvT+NFlTUE3C+eXz33ToiQxdNvI5dEOBquq3nGtOa7yOyCGnEl0oIxcs1pnkdHVtl4N5ze8ds/uaCUhg77ALGndI2jqvxRxIbqSHeWwM75KGxsRGbY7EmO4WQUbhDhBhEtggF6TWmF3Qcdfq8oqKizyIEeILUCJLjZQy6xgK5srx1hHzHfMHSwdds04ynnF9Z+VlCW+skWH+pUHMqY4QkjASzWoFWMDaH/mLxQPXeSMA9TETri881ljhrSlS4Q4QYxC5be5ruGw2WNtUi5L2tzCDTG16uI2S0oGIA15iyrQcAuAwGS/usM+r1eOdQn+u0CxAEtRCKKGtMfs16rjF9N5nVaC/aiXgRTwaUv79E3MPUBBRCKawjZBjuECEG0Vop9IWQ96faNRbdvMoWG5J1xqPIIAsWIyRbhEK5xpo+CXQtQkohFGHTVeX8WouQ5PZwRPmtX5rCJ9Z8j8UrfV67T6zwGxlqa2HiXbYS0Q3lTMA1JSrcIUIMohUcunWEmn6qg6Wju/gp4470Cirqje6r7Oy9H7lFSO0acxvMGtO2u1DOry2oKAnMqC1CmhIAgYKl45k+H228WEtF+ftLRDGZiNYXyUrF9PnQcIcIMYj2Azhos1OdOJVIkYZQ1g4SRX0XkG+t3p/uphYbhtPdder+6LfYMBhzpBPQrQ3Ili5yaotQ9FljSoVopPKvFWjfM4no1kkGlPuYiOngiRkjJDT9TJw1JSrcIUIMYsgipBEAZnxkqwoqKjO4gqTGay/ARgsq6tX9UR4LVkNJD6mgoyAI/kKlCUkcKJvYRpM1pm8Rik+MkBHxTELjUFmEEu+ylYhuqEQUZ4kKd4gQg/iJC52/Hq1LyIyCmHotNoDgVZ6180baYgPQVpY2Zunyc40J+r3AgEAWochjhPTcd0oLUzyFUAIaM5ICR4JbhHLTnQCAnDRnnFfiQxZntAiFhHWECDGIvxDS+YCRL8ben2YYAFRNVxXHXUFS2bUWIKOVpUN1izfca0yyjMmiR2h6juhvEYpFjJDKNRbDYGnt74EXpYhwRPnesJoze+XjnvG9Eqq+kWQJimUl9WSFQogQg/hd1ALrIIVrzLwPbWURRUBxwde5MJjpGlO32AjccFY1nrRGRYp/oFpF0lqVWXnR1BHyWa18Y6Qou8/bY9eEUmvZCrVvRB9V1l8CusZSHXbccU7PeC9DhdTklcHSoaEQIsQgWuGj3/Xd/GBpOUYI2hihIL3GInSN6QVLK61Ech8vo64xj79rzKVJn5cubGbVEdJzjakupDG0CDFY2hyifW+0RC4YWIT1e4/jt4OK472UhIdCiBCDGHE3aeNUTBFCTddw/xihIC02NMdCGVh86fY6TVIVuqXRaEFFTRaaAMHPYiMhWX+izRoTZKGVOL3GtAKUwdKR4YzSWtgS6d8+Fx/+z8h4LyMpoM2MEIMYKY7ny/Ay2JPLAHLTVTFAcUPd9PkwLUIB3FbeeRWuMbd/y4xAqwbUrjGbRqhISBlBytcRWdaYFz0RGi8h5Oca4yduRNijFMmEBIPvKEIMYsgi1PTTbWr6vPenssWGcg79pqsaIRRm01Xla/OoYoTC6zWmEoQBLEJOvRghk4OllQHSsYyZ8I8r40duJCjT52kRImbDv0pCDOJnEdITF5qsMVPT5xX/e+cIHCMUbtaYtjeY8nS9uKFQr0vrIvSOGcgiJGWNRVdHSJv5lggWoXCD1ok+0ZZWICQYFEKEGCScytIeE2OElBld6gaoga1O/jFCRmN6/K0perWFjAornyD0tzpJmFVHKJhFSHKn2ITYBttqNRddY5GhdIclYq8xktzwHUWIQfysLMF6jZnoGpPG8GaNGSuoGG6QrhBERHh0AqiNtuxQusZki5Bb23RVsggpY4QiCZaGanzlGiXXWKzbDWhfB11jkUGLELES/lUSYhD/goqhLUJmZAkpW2yoLEJyJ3idtRqoeaREdmWJwWOEfOeHijnyD5bWii3f2nQsQhHVEQos5uLVAiHc3wPRJ9r4MUKCwT9LQgxiRAhpBYUpFiGFdUU/fd6ARchgsLRHJxhaL0YoZKyLxg0mwD+GR8JpUmVp2SKk45ZMd3qLy6U5Y1dMEfAXqew+HxnRimRCgsE6QjqUlJSgpKQEbrc73kshCYSR9Hm/+jkmXPiULTaUBAuW9iuoGLIStDSHv0VI1HONhbIwyc+FvEaftUx9rsPkGCG9Pm892mbh2lFd0KcoO+xxo0Gb6s3U78hQZY1xD4nJUAjpMGPGDMyYMQMVFRXIzc2N93JIgqC1ggRrdqpnWYkWEaLhVHb/IF2jWV5oGjN4sHRoC5PaTaUMlvavLO2fNWZ201VBEDDzt/3CHjNa/KqR8xoeEYnea4wkN/yzJMQgRsSFHCytaDYaLZIlyqMpqBjMNWZEtOk9rmdN0cQ2GxpPz0UoBLBs6VWWjqQ5qdZylghuKL86QgmwpmRE6Q6LdcA7af7wHUWIQQw1XdUIAHO6z3t/htNiQ+suM/ot2ldQ0XdM1zVmsOmqqBBWgawhktvDblJBRXkNCSA6/FxjjG+JiGhrTBESDAohQgziHyzt/+djSR2hAC02gjVd9V9r8Dm0FiGl5SKS9Hn/DC7f69AiV5aOtumqEPx+PGCwtDk4Fb9MJ8UkMRkKIUIMYqRKsHRErzBhpChbbBh2jYWZNaYtqCiohJDOmkIWEmpao07TVS2yRSjKYGntS0wE0WEk05CERhkszT0kZkMhRIhB/DKxdF1jTZYQWQCYgOQagzdgWkIqHKg3R/hZY16UgdFaN19440lWLP+mq1ocJlmEtJaxRLhehtvqhOjjUFmEeNki5sJ3FCEGMWQRkiwrQVLbwyVQQUXZIqTzV+zXYiPEBVgOyJayxmwK65Zu1lioNTc9V2FhCrQVer3GIkmR9hNatAg1G5SxVdxDYjYUQoQYxD9YWs815p82Hi3KIfTS543ECBm2CMlZXoIiy0uvjpDBrDFVQcVAFiH/OkLNJUaITVfNQWURYg0CYjJ8RxFiEGMtNrw/9TqgR4qq75fCOhO0snS43ed1gpt9Vp3ga9IdTyrQqGq6GiBYWqeytBlZY4kQI2SkCCcJjSpGiMHSxGQohAgxiKE6QhrXmBkXY+UQboMFFf1dY8bmUKa7a+N8go3v93jTfMqg8cCuMR2LUAQXO/9g6bCHMB2toKMQioxoRTIhwaAQIsQgRiwOfoUJTZhX6fqSGq2GWkf0rjHfQb3K0qFdPFLQuPFg6WjrCGldhIlQRygRrVTJiJNNV4mFUAgRYhBDga+ya6zprskWIf1O8P7PCb/pquQa843pE0d6azIWI+RRxggF+LRx6lmEmmmMEAsqRgYLKhIrYa8xQgxiqPu8HBtjTYyQ0b5f4aZt62W7ad1l6jmDDqeopySNLwT81mVV1ligAo6xxO/3QItQRCiLbiaCpY80L2gRIsQgRsSFHCztMdE1prhtNHA54jpCojJYOnDWWMh0fJ2mq4GeIbvGFGuMrI6QZg0J8Olms6ljo0IWoiS6SO8HWtSIFSTARwUhyUFYdYRMrSwdyiKk85wwg3S1qfKCwp7i0Wm6atg1ptN0VYucPm9yr7FEicdRvk8Y3xIZPvcpL1nEfPiuIsQgRsSF7BozMX1elTWmI4T06wgFvx9oDqVLL2gdIcOuMZ8gDNxiQ3KNRRsjlHjB0oD6fUOLUGTQIkSshEKIEINoLUB6FzWta8zs9HnDwdJhWkdkV5YqfT7wnOHWJdLLGpMygfSCpSO54CVi+jygCQJPEHGWbEjvB1rUiBVQCBFiEENVgmUrinnzRhIsHXZBxaafcgFExUGjc6rGk11j0niC33PSnXYAQKrDnKar/lljiXHRjDb2idA1RqyFWWOEGMSIuNDrsRUt6mBpY26qSLvPK/uX2YKIulAvy4iL8LZxPbD9YCVGdm8DQJs+H/4FLxGbrgJqyyGFUGSkNAmhFAeFEDEfCiFCDBJO+rw7SNXncDHDIhS6N5hOsHTTU/TS50NnjUE9nuCfzj6hbzvcclZ3xZjK9PnoY4TMydmLHrrGoqdfcQ4uGdIBp3VtFe+lkGYIhRAhBvELQA6WNWZm09WQwdL+z/FrsRFmlpdN8K8FpB7f4Hiyq03wS2fXjtEcCyoCPhEqCAyWjhSH3YZ/TBsU72WQZgrtjIQYxN/KoneO96fbxPR5pcvHqCjxc42F+EvX1v2BIPgFPKvPDz6eXxNXm87+BXE1Nsf0eVqDCElMKIQIMYgh15i2VYVJc2tjeJTousaMBHYrx4d6fAEImjVmPFha6WoTdM+RUGaKmRFLkyhxtdJroTWIkMQkQT4qCEl8jARL+8XVmGQF8LPYqB7zPz/iFhuib0ytmNFbT8DxdIKltc8IVpfJlBYbCWKBkWvgUAgRkpBQCBFiEEOVpTVd18269skWG906QqGDpUOJAm07DUHhGosmRkgp3ELF8KhihCKoIxQqBileSO8busYISUwohAgxiDHXmPenmb3GgPAtQlpRYNTVpMx2C+oaC1WpWn6utEb/OkL+FjZzs8YSxQAjrYOuMUISEwohQgyiFBPKFhRK5LRxEytLeycM8pBesLRWZIRrEVKkz3siKKgoj6fTssO3bvVzzM4aSxTZIfdSoxAiJCGhECLEIKoKwQGEgNY1ZpYOCnYNNdJ0NaQFR+vKEvxfi3pOg64xRdPVUJWfVb3GItg4/4KKiSE8pN8FLUKEJCYUQoQYxEiFYL/WEiYHSxt9zE9kGMwa0w2W1uk+b9w1pmy6GlyoSBYTmxCZaNA+I3GCpZt+Jsh6CCFqKIQIMYiR7uiC1iVk0tzBxtG7vhqJZ1KiFR7K3mB6laVDW4SkffCt0T+YWX/NkfaTStQYISn2ie01CElMKIQIMYgx15gX811j4VmElOsLFM+kRPuoUotElTUmiyidOkKaWaU6QpEKhsQtqNj0k0KIkISEQogQgyjFQSDXjXRYVGRLmULQGKHA2WtG16AnUvzFjI/QrjZ18LVN0Enp98tsiy6oWLukRCuoSCFESGKSIB8ViUVJSQn69u2L4cOHx3spJIFQWYRCuMZ8982ZO7hFyP+Ycn1Grr/a4QVVjJBO7aKQLTu8P5V1JUMFS0sCKJIaQoCOey9BLELS66QQIiQxoRDSYcaMGdiyZQtWrVoV76WQBMJYjJD6vlkWoWDD6KbPh5mKrnVTCYLgZ9VRYtQ1phzfz/2mOZDi8H4cpWi72xokYdPnJZdfgggzQogadp8nxCDeasteK0fAGCGLLnbhWoSU5xu5AOuJCLkmkk6MUOhu9pp4HQNNV3sUZOHy4R3Rv31uyPXqzong48cLaR1MnyckMaEQIiQM7IIAlygGtghp7ptmEQrymH76vOK2EYuQnyVLkfml29YjvPGMNF212QQ8ccnAkGsNRKgWHvGCvcYISWzoGiMkDHzF8QI8blGMUDBLk27WWJjFCfVdY14iqSzt5wgzECMULf5CKzGEh1wfiUKIkISEQoiQMJBERWDXmOa+SfMGvaaHcI0Z0QO66w6WNWawm71yPVant1sttCJFDpZOjOUQQjRQCBESBqFSoa1yjUXTYsNQsLSONUW2COnWEQo+nm7MUYiCitGSuAUVJdcYP24JSUT4l0lIGEgX10DiwirXWDQFFSMKllbECOlVlg5doNF/H6x2XfkJoQRRQvYQ7lRCSHzhnyYhYeBoSu0OKEz8BEWcgqUNFH9Uja8TaBwoa8yYhcl/jaoAbis0SoKmz7OgIiGJDYUQIWEQqjien0XIpHmDCSrdXmMqi5CB8bUWHEXlH7dGCUVyPReg3hsr4nf8rVqJITzkuDKahAhJSPiXSUgYyJ3EDcYImZc1FvixUN3nI7II2QLPaURg+K1JiIUQSswYIWn/GSxNSGJCIURIGNhDWISsqiwddkHFcNPnQzREVWJsPM16LIqdCjZHomSNOegaIyShoRAiJAx83+5jGywdrkVIedE1ZsHRThj4eYZ6l+nct9oipBfnlAjYKIQISWgohAgJA18GkDGLkFlxKsGEQ8gYIQN/5XrtKQJdt42IGP8GqLA8WNpfhCaG8AhlRSSExBcKIULCIFRBRS2mBUsHOh7AcqM8FIkrS9A5JmEo5kj7HEFQPc8KkWKVWzJaZPGcIOshhKihECIkDKSLuSNA5KtVcSoGs/VllNYHI8JFr09XoDghQ4YNHReh8lAsYoQSRXew1xghiQ2FECFhIFlXAgkcf9eYOfMGjtfRP652jRlSLn7zBbQIRRJzBMHyGKFEb7qaKAUeCSFqKIQICYNQga9W1REKN15HeTgS4SIg8NqNucZiHyOUqE1X5YKKCbIeQogaCiFCwkAKPA4oQDT3TXONBZAlwWr9SGLDUJaXXq+xaLLGdKwzVluEpHmsniNcJAEUyJ1KCIkvFEKEhIEvA0j/cb9rr2muMf3jwS724bR28Fu2ENmcvnO046uFlVXWGrXYsmSKsJEEEJuuEpKYOOK9AEKSCTlYOsBFTXuBt7qgYvCu9AIA0WC6u/+4AV1jRrLQ4uAak+bxzZEYSmjSgCKs2n0MFw0ujvdSCCE6UAgREgZysLTRFhsmzRuJdUZaojGLkDa2KYhrzIhhQyeVXTmHVSJFaBJ/3tuWTBE2vdpl470bT4/3MgghAaCtlpAw8AW+6j9uVQq30Sw11XPCqF+jl+0WjUVIv/+Z/m0zUY6bKMHShJDEhkKIkDAIt7K01XWEgmVwhYpnUo9vPH3eUIFGv/FheUFFIDFjhAghiU1UQujXX3/FvHnzUFtbCwAQRdGURRGSqIRKhU6UOkJAeBWN9YOlw7dCBTrHphFWVsUNxyIzjRDSvIjo46i8vBy/+c1v0KtXL0yaNAllZWUAgBtvvBH33HOPqQskJJGQLq6BUqGtqmUT2E0V5Dly+nz4rqyog6X9ssZiI1LUwdKWTEEIaWZEJITuuusuOBwOlJaWIiMjQz5+2WWXYe7cuaYtjpBEI5SVxapg6cAXdSOusUiESxDXmKGWHaGyxiwSQqo5qYQIIaGJKGvs22+/xbx589ChQwfV8Z49e2LPnj2mLIyQRMQWQlxY12ssgvR5WbQZGF9n3MBFHMN/Td5u9soYobCHMDaPsscahRAhxAARWYSqq6tVliCJo0ePIjU1NepFEZKohKwsbVGMULgtNpSPGbMIGQ+WjqRSNTQxR9ZVlrZebBFCmhcRCaEzzzwTb7/9tnxfEAR4PB489dRTGDdunGmLIyTRCFWt2b8ejzlE0gleCmOKJKYn2PMirSwd6/R5xggRQowQkWvsqaeewtixY7F69Wo0NDTgvvvuw+bNm3Hs2DH89NNPZq+RkIRBDpYO6BpT3zctWDrAMMHGD9UgVnWujksvkpR9eV3QjherYGnrU/QJIc2LiCxCffv2xYYNG3Daaadh/PjxqK6uxsUXX4x169ahe/fuZq+RkIQhVB0hreHG6oKKwdLQ5SrYEaS7By3UGOF4sdAlidh0lRCS2ETcYqOwsBCPPvqomWshJOGRM7ECZo1ZFSytfzzY+JJFxJgFRztuZLWLlM9Xjy/EJEZI3cbDkikIIc2MsIRQaWmpofM6deoU0WIISXRsISxC/gLApHkjKqjY9NNQjJBOsHSgcQ2JGD3XmOK+ZQUVlbephAghoQlLCHXt2lW+LVWRVn6AiqIIQRDgdrtNWh4hiUVIi5BFvcYCxwgFfk6oBrHBxhEQ2KISSWVpaNLnYxMjZMkUhJBmRlhCSBAEdOjQAddeey0mT54Mh4PN60nLYnjX1vjiv/sxuFOe7uNa8WB9HSEDwdIRVJb2ps9H4xrzz56LRUNUpaWJwdKEECOEpWT27duHt956C2+++SZefPFFXH311bjhhhvQp08fq9ZHSEJx6dAOuOjUYjgCdDL1t4SYM2+gYYINbwvHIqS9LwQe21BdIp21qGOEQg4REWy6SggJl7A89YWFhfjTn/6ErVu34tNPP8Xx48cxYsQInH766XjllVfg8XisWichCUMgEeTFmmDpiAoqyjFCocf3b5IamTsu0DlCjNLn2XSVEBIuEYcsjh49Gq+99hp++eUXZGRk4JZbbsGJEydMXBohyYdVwdKRdIL3pc9H4soSEGj1ZrjGrLLWKIelRYgQYoSIhdDSpUtx4403olevXqiqqkJJSQny8vJMXBohyYdWsMTVIhSGa0xLsLo/RlxjWg1ls2l7jVkVLK28TSVECAlNWDFCZWVlePvtt/HGG2/g+PHjuOqqq7B06VL069fPqvURklT4V5Y2a+QA1pkgX2WiDZYOLL5CDqe7WiEGFiG6xggh4RKWEOrcuTOKi4sxffp0/Pa3v4XT6YTb7caGDRtU5w0cONDURRKSLFjVaywSi5Dca8yA3VcvfT6a7vN6ZQRi3XSVrjFCiBHCEkIulwulpaV47LHH8Ne//hWAr56QBOsIkZaMf5CwWenzgY6bkz7vHywduNeYMQuTZny/pquxcI1ZMgUhpJkRlhDatWuXVesgpFkQTs+ucAhcWTr0c4zECPm7xoK09TBiYQrRdNUqkRKLOCRCSPMibNcYISQwVvUai6jFRhhZY3rlj6Jzjfnfj4VFSCnSGCNECDFCxFljixcvxtVXX42RI0di//79AIB33nkHS5YsMW1xhCQberE25gysfzioRSisOkIaAWeLzjXmvw+xKajIpquEkHCJSAh99tlnOPfcc5Geno5169ahvr4eAFBZWYm//e1vpi6QkGRCz8Vkxbi+8Q1YhCLoNRZsbGNZY/77EJuCisrbVEKEkNBEJIT++te/4sUXX8Qrr7wCp9MpHx81ahTWrl1r2uIISTYsC5YO87hy7ohbYgQ415CrTWcfYtFrjE1XCSHhEpEQ2rZtG84880y/4zk5OawuTVo0VtURiih93mY8Rii8YOkIxtMcsy5YOvAaCCFEj4iEUFFREX799Ve/40uWLEG3bt2iXhQhyYs1wdIB3VRB/oKLc9MAAEVNP4OPrxlXCGxtMuQa0xtPUN+3ApXYsmYKQkgzI6ysMYmbb74Zf/jDH/D6669DEAQcOHAAy5Ytw7333ov/9//+n9lrJCRpsCpYOqB1JojQuvs3PZBXsQMT+rQNPb5fIUgh4NgRZaEJQoxihBRzMFqaEGKAiITQfffdh5MnT2LcuHGoq6vDmWeeidTUVNx77724/fbbzV5j1EyZMgULFy7EOeecg08//TTeyyHNGO0F3jSLUASp7KlOOzpnRxYsLXjz53UxNl5wV5tlbqsYWJ0IIc2LiNPn//d//xdHjx7FypUrsXz5chw5cgSPPfaYmWszjd///vd4++23470M0gLwt4SYM240fb+MoBfcHEh8ReIa01qEYhEjxIKKhBAjRCSE3n77bWzduhUZGRkYNmwYTjvtNGRlZaGuri4hBce4ceOQnZ0d72WQFkAsWkdYMZ9ecHM0c+oWaIxFQUU2XSWEhElEQujaa6/Faaedhs8++0x1/OTJk7juuuvCGmvRokWYPHkyiouLIQgCPv/8c79zZs2aha5duyItLQ1Dhw7F4sWLI1k2IZaj17PLDCJpsREO4QVLR5aFFouGqGy6SggJl4hdY48++iiuueYazJw5M6oFVFdXY9CgQXj++ed1H//oo49w55134qGHHsK6deswZswYTJw4EaWlpfI5Q4cORf/+/f3+HThwIKq1ERItZhklImm6Gtb4fgUQA1eWjqSOkC1GwdIxiUMihDQrIgqWBoCrr74ao0aNwpQpU7Bp0ya88847EY0zceJETJw4MeDjTz/9NG644QbceOONAIBnn30W8+bNwwsvvIDHH38cALBmzZqI5tZSX18vV8kGgIqKCgBAY2MjGhsbTZlDQhrP7HGbM8mwZx63S3Vf9HhMWa8oioEeCDh+OPvl1qzb43EDAecM/ZrcLrd6PLdbtTciAq87GgT41uxyucKeIxneY4kG9yw8uF/hE8mehXNuREJI+hZ6+umnY8WKFfjtb3+LUaNG4cUXX4xkuIA0NDRgzZo1uP/++1XHJ0yYgKVLl5o6FwA8/vjjePTRR/2Of/vtt8jIyDB9PgCYP3++JeM2ZxJ5z47XA8o/qy1bNuPrY5uiHndvqQ16BtxDhw7i66+/DvpcI/tV0QCo1r15M/ZVC7pz7t69C19/vSP4eqvU423cuBFHMkT52IF9+/D116W6z42GI0d8+7T0pyXYkxnZOIn8HktUuGfhwf0Kn3D2rKamxvC5EQkh5bfTTp06YenSpbjqqqswfvz4SIYLyNGjR+F2u9GuXTvV8Xbt2uHgwYOGxzn33HOxdu1aVFdXo0OHDpgzZw6GDx/ud94DDzyAu+++W75fUVGBjh07YsKECcjJyYn8hejQ2NiI+fPnY/z48ao2JSQwybBnByvqMHPtIvl+/379MGlEp6jHXfnVVvx0aK/f8eKiIkyaNEj3OeHsV3lVPR5e86N8v3///rCVVWLZ4X1+5/bo3g2TJvQKOt6Wsgr8feNy+f6gQQPRoyATz25aCQDo1KkjJk3qF3SMSJhTvhZbTxwFAJw5Zgx6F4aXJJEM77FEg3sWHtyv8IlkzySPjhEiEkKPPPIIsrKy5PsZGRmYM2cOHnnkESxatCjIMyNDGwchimJYsRHz5s0zdF5qaipSU1P9jjudTsvesFaO3VxJ5D1LcapdQg6Hw5S1Ouz64Xx2uy3k+Eb2y+n0qO87HLBr5rQJgEcEHHZ76PEc6scddjtSFM9xGFh3JCj3KSUl8vdJIr/HEhXuWXhwv8InnD0LZ28jFkJ66LmVoiE/Px92u93P+nP48GE/KxEhiYBVdYQCd4K3poWHoJM15rDZ0OD2GGviGrKOkPVNV5k1RggxgmEh9OWXX2LixIlwOp348ssvA54nCAImT55syuJSUlIwdOhQzJ8/H1OmTJGPz58/HxdeeKEpcxBiJtoLvHlCRf+4WRd7v2ax8F+73SYAbmMixr/Cdqx6jflus6AiIcQIhoXQRRddhIMHD6Jt27a46KKLAp4nCALcbnfAx7VUVVWpGrju2rUL69evR+vWrdGpUyfcfffduOaaazBs2DCMHDkSL7/8MkpLS3HLLbcYnoOQWGFZr7GAVZ6tSZ+36aTP+7rZGxhPp2UHCyoSQhIRw0LI4/Ho3o6W1atXY9y4cfJ9KVh5+vTpePPNN3HZZZehvLwcf/nLX1BWVob+/fvj66+/RufOnU1bAyFmob30mldQMcB8ZrnetCFIOq4xaQ32iCpLx76OEGUQIcQIEdcRMouxY8cGrpHSxG233YbbbrstRisiJHL8LvAWF1Q0zyLkP67WtSQFIkfadDUWFhqBFiFCSJhELIRWrlyJhQsX4vDhw34WoqeffjrqhcWTkpISlJSUhOXiIwSwzjVmfYsN/15jWiTXmBF9oR8s7bsfC9cYdRAhxAgRCaG//e1v+POf/4zevXujXbt2qg/R5hCgOGPGDMyYMQMVFRXIzc2N93JIEmFVsHQgRWWV681m0wmWbrofmWssNhldKrHFtDFCiAEiEkLPPfccXn/9dVx77bUmL4eQ5EYvSNgMtJYOyZtsWa8xBAuWDj9rzOsaUzxukUhh01VCSLhE1HTVZrPhjDPOMHsthCQ9sQiWdtjUosgMdLO8NOdE4xrTNl21ynDMpquEkHCJSAjdddddKCkpMXsthCQ9epYQM1BabKyweujF9ASyCBkqqOhnYYpNartyXuogQogRInKN3XvvvTj//PPRvXt39O3b16+U9ezZs01ZHCHJhp6gMAOtRahePm6Va8x/7eG4xkLXEYpomSFRFVRkAj0hxAARCaE77rgDCxYswLhx49CmTZtmESBNiBn4WYTMGlgxrtIiY1mwtOAvIxxRFVQUYuK2YowQISRcIhJCb7/9Nj777DOcf/75Zq+HkGaFecHSvtt2S2KE/F16gTLgIqojBGhihCwSQgpnP2OECCFGiChGqHXr1ujevbvZayEk6dELEjZlXFhrEdJPd1cfC8s1pr0vaCtLh79GI7CgIiEkXCISQjNnzsQjjzyCmpoas9dDSFJjlWssUDFCS4OlNeeE02tMr+lqbAoq+m77tQ0hhBAdInKN/fOf/8SOHTvQrl07dOnSxS9Yeu3ataYsLl6wsjSJFD1LiCnjBkifN6+7vZ5rTH2OI+pg6RhYhAJk1xFCSCAiEkLBus83B1hZmkSKnqAwe1ybKkbIvIu9TQA8onRb8Mu6Sk+xAwBSnfaQY/m72tQtNiyLEVJljRFCSGgiEkKPPPKI2esgpFmgtXSYdTEObBEyaQI0iZOmktWCzti3ntUd3fIzMa53gbGxVPdjVEeIMUKEkDCJe/d5QpoTVvUaswW0CJkyvHcs5W2daOmhXVphVI98Y2PpxBzFIliaTVcJIeFiWAi1bt0a27dvR35+Plq1ahXUtH3s2DFTFkdIMqLuBWbSmIrbVsQI+cYS5dt67i2j6GahxSC1PRYB2YSQ5oVhIfTMM88gOztbvs0iioTo45MTFlmErKrHI6hv+5cCMD6Uf9ZYbHqN2SxyGxJCmi+GhdD06dPl2+w6T0hgBJVJyKwxfbcddmsu9lpril/bjTDUS6gmrlZ9kVKOSosQIcQIEVXaeO2113SPu1wuPPDAA1EtiJBkxwr3jFI42C0KCFY1LEV0gd+hm66Gvz5D8zJGiBASJhEJoXvuuQeXXHKJKhbo559/xmmnnYaPP/7YtMURkoxoBYUZBGqxYW7WmPq2nlXH8FiaT5bY9RpTz0kIIaGISAitW7cOhw4dwoABAzB//nyUlJRgyJAh6N+/P9avX2/yEglJLrSCwpQxFbftltURUlua/FPgowiW9kufj2iJIZH7oVEDEUIMElH6fNeuXbFo0SLcddddOO+882C32/H222/j8ssvN3t9hCQdVlg+bAEyxcx1jenfBsIXdPpNVwM/bhbSHIwPIoQYJeJuPP/+97/xwQcfYNSoUcjLy8Mrr7yCAwcOmLm2uFFSUoK+ffti+PDh8V4KSUKscI2p0uctCpYWNEJFdT/csTT3bTZtHSFrCypSCBFCjBKRELr55psxbdo03HfffVi0aBE2bNiA1NRUDBgwoFnECM2YMQNbtmzBqlWr4r0UkoRYYfkIVDHZVIuQJtBYJejCnEev+azaUhbREkMizUEdRAgxSkSusZ9++gkrVqzAoEGDAACFhYX4+uuvUVJSguuvvx7Tpk0zdZGEJBNWZC4FarFhahkhjUsvmr5dgZquSpUFrAuWpkWIEBIeEQmhNWvWIDU11e/4jBkz8Jvf/CbqRRGSzETjUgqE8sIei2BpbdZY9NP4BIpbFK0rqCiofxJCSCgico3piSCJ3r17R7wYQpoDqqJ+Jl2RA2WNmRojpLltpmtMK1CsjhFi6jwhxCgRN1399NNP8fHHH6O0tBQNDQ2qx9auXRv1wghJVlSuMZPGDGQRsi5GKMpgaT/XWJNAaWpAYos4TSM4NlkIWTM+IaT5EdHH0T//+U9cd911aNu2LdatW4fTTjsNbdq0wc6dOzFx4kSz10hIUmFJmriqoKLvz9bKgoqBHjM0VoD7gsUWIabPE0LCJSIhNGvWLLz88st4/vnnkZKSgvvuuw/z58/H73//e5w8edLsNRKSVFgRLK2yCFlUj0fl0tM2SQ3TJqTXdFX507JeY4wRIoSESURCqLS0FKNGjQIApKeno7KyEgBwzTXX4IMPPjBvdYQkIVpBYQa2gBYhi4KlEV2wdCCLktXBzMwaI4SES0RCqLCwEOXl5QCAzp07Y/ny5QCAXbt2QZS6bhPSQrEiRkhQCSHfbUt7jSkeC1dYBLL4WC1UGCxNCAmXiITQ2Wefja+++goAcMMNN+Cuu+7C+PHjcdlll2HKlCmmLpCQZMOKXmMxCZbWzBetoNNrNWK164rp84SQcIkoa+zll1+Gx+MBANxyyy1o06YNFi9ejMmTJ+PWW281dYGEJBs2HQFgJnbLCioKmvvKOxGMB0CyD8sCyGatxYZZY4SQcIlICNlsNjQ0NGDt2rU4fPgwUlNT5UKKc+fOxeTJk01dJCHJRLiBxUZQB0tblT6vni96i5DgLSMN+GWLWaVTmDVGCAmXiITQ3Llzcc0118hxQkoEQYDb7Y56YfGkpKQEJSUlSf86SHywovu8akyLLEJ+laVV84c/kU0A3PJtQT6mnctUGCxNCAmTiGKEbr/9dkybNg1lZWXweDyqf81BPLDpKokGraAwfUwIlggKv2BplaCLYDz4W5Tk7vCWFVSU5rFmfEJI8yOij6PDhw/j7rvvRrt27cxeDyHNCgvqKcImWBML4xcsHUWLDe2A0tO1gshsmD5PCAmXiITQpZdeioULF5q8FEKaBzZVertZrjG1lcmKC37QOkIRjKfnWrNaqDBrjBASLhHFCD3//POYOnUqFi9ejAEDBsDpdKoe//3vf2/K4ghJRvRcQlGP6RfI7LttGioLjqARMuEPpxVW3mPqn2ajFVyEEBKKiITQ+++/j3nz5iE9PR0LFy70a9ZIIURaMoJGUJiBTWOe0QYfmz2H0urkvR/+RHr7YLVQkcelDiKEGCQiIfTnP/8Zf/nLX3D//ffDZlXUIyFJijXB0urxfUHBJgZLa2+b6BqTLUFNHxeWJY3J81EJEUKMEZGKaWhowGWXXUYRRIgO2qBjU8bUiBIrLEJ+7rcAjxlFr2mr5TFCNmluS4YnhDRDIlIy06dPx0cffWT2WghpHkRpSdEdUlNE0YoYIb86QjpCJix0ssakqth2i5QKs8YIIeESkWvM7XbjySefxLx58zBw4EC/YOmnn37alMURkoxY4RrTWmekoorW1eOJ3iKk9/wbRnfFj9uOYGCH3GiWF3hONl0lhIRJREJo48aNGDx4MABg06ZNqsf4AURaOla4xrSBy1bU49GuNdoyAOoK2N7bV43ojKtGdI5sgUbmZPo8ISRMIhJCCxYsMHsdhDQbrHDLBIoRMnMmbRuPaHum+QVfxwC6xggh4cJoZ0JMJlBfsGhQXtiVDVEta7GhvR+Ja8yi5rBB55Tnjsl0hJBmAIUQIRZiWrC04ra3xYZ027pgadX8EWWNRff8SGCMECEkXCiECDEZS4KlA7bYMGd8QKfXWLRZYxZU2A5FqtP7kZZq50cbIcQYEcUINXdKSkpQUlICt9sd76WQJERbj8fsMQWrCipqe40pHouo+3wcgoRGdmuDa07vjAn92BCaEGIMfm3SYcaMGdiyZQtWrVoV76WQJESvx5apYypq/FhVUNErtvyzvsJBWw07FqQ57Xjsov4Y07MgJvMRQpIfCiFCTMaKXmNaK5NcQdlEJaSt+xN99/nYu8YIISRcKIQIMZloCxHqYdOIErvVTVehES9RusYYvEwISVQohAgxGSvSxrVjThveEad1bY1+xeZVaPbrNRa1RUg5XqSrIoQQa2GwNCEmE62A0B1TM/5tY3vgtrE9TBpdmkOb7RadoIs+64wQQqyHFiFCTMaK9PloA5eNECgzTftYJONRBxFCEhUKIUJMRm29sSB93pQRQ8whRG/RUbvaolkZIYRYB4UQISYTbWsKPdQtNswZM9gc2mDpyCpLW2/FIoSQaKEQIsRkrO6xFQvXmDZYOqLxAtwmhJBEgkKIEJOxQgDEwiKkDZaOtkJ2PJquEkJIuFAIEWIyVrjGYlGTx98iFF3QtxX7QAghZkMhRIjJ2CywhFiRiaZFK7CijRGi9iGEJAMUQoSYjBVCxYpGrlq0vcGizxqja4wQkvhQCBFiMtZYhHy3LUufV94WzMgai+75hBASCyiECLEQ82KErLeuqC1A0cclsekqISQZoBAixGQs6TUW8I55+LnGohQysXDnEUJItFAIEWIyVrixrHC3+aMOyDbTtUUdRAhJVCiEdCgpKUHfvn0xfPjweC+FJCHRxtbojhmDGCGb1hUW5ZysLE0ISQYohHSYMWMGtmzZglWrVsV7KSQJsUIAqCxCFv3VSlPIPxHd69CORwghiQiFECEmY7UAiCSVPZxxpdGjbZrqE1SEEJK4UAgRYjreS7+ZsTw2m9I6Y9qwmjmafgpqQeS9Hf6k8jg0CRFCEhgKIUJMxmaBJUQdd2SxRahpeFuUUd/SU6zqjUYIIWZAIUSIychCwkyLUCyarsouPT2LUOQDWuXKI4QQM6AQIsRk5Au/idf/WNTk0QqgaJumyoKNOogQksBQCBFiMlKsjamusRi22PDNFV3tIuogQkgyQCFEiMkIFgRLR5vKbgTJguMLclbMH1HWmPn7QAghZkMhRIjJWJE+H4sGplrXmKoeUkRZY9K40a6MEEKsg0KIEJOxwhISixYbWgEUbYVsbV0iQghJRCiECDEZK2JjYhIjpKmAGLXesiB7jhBCzIZCiBCTsSJbStXR3uoWG9L9KIOlmTVGCEkGKIQIMRkrXGNqi5DFwdI2k4Kl6RojhCQBFEKEmIw1wdLWt9gIJlwimVK3QjUhhCQYFEKEmIw16fOK25YVVPT+tOlYtKLqPh/1ygghxDoohAgxGSsEQGxabKhdYtEGaLPpKiEkGaAQIsRkfPVzkitGyK/XmCpGKPI56RkjhCQyFEKEmIy2i7spY6p6jZk3rhKbxpKlrmYd/ng+8UQlRAhJXCiECDEZS3qNIUo/VRhzmOca8x+HEEISDQohQkzHisrSytuxCpb2fyys8SCNF926CCHESiiECDEZK9LnVQUVLWu6qk2fj67XmK93GZUQISRxoRAixGS0XdzNHBOw3tWkFywdSTVrusYIIckAhZAOJSUl6Nu3L4YPHx7vpZAkxAoLiNoiZPrwTeNqYoSU80dWUlE1LiGEJCIUQjrMmDEDW7ZswapVq+K9FJKEWOEaU49ndfq89DO6aGnqH0JIMkAhRIjJ6FVmNndcU4eV8QU36wRLRzAeXWOEkGSAQogQizDdIiSPa1GwtNRsVZ4vugBtK1qNEEKI2VAIEWIyzcUiFHX3eVqECCFJAIUQISZjWbNReVyrYoQCLzyqXmORL4kQQiyHQogQk/G5sMwd1+qYG60OirrXmAU91wghxGwohAgxGTnWxmQBYEUPM/X4XvS6xkeePE/XGCEksaEQIsRkrGo1akWhRvX4aqFli9IiJNA1RghJAiiECDEZwaJgab2Kz+aO7/3pi+0R/B4LBxtdY4SQJIBCiBCTsbqgomVNVwPMp/dYOOOx6SohJJGhECLEZKyq92N1FpagiQ1StdiIKH2eTVcJIYkPhRAhJmOVYDn7lLboXpCJjq0zTB7Zi19skMoiFEmMkPonIYQkIo54L4CQ5oZVAuCZy06FKIrWVZb2C5ZWVJaO4CuTL8uNSogQkrjQIkSIyWjT0E0d20JR4Zc+r/toGONZVViSEEJMhEKIEJOxOrvLKvwLKpqVNRbdugghxEoohAgxGZ9rLLkUQNBg6UjGY9NVQkgSQCFEiMkka48tbWwTm64SQloCFEKEmEyy1s/RCjilRSsSqw4rSxNCkgEKIUJMJmldY00/tdljysfCGi9J94EQ0rKgECLEZJLVEuLnGlM9FoFFSDMuIYQkIhRChJiM1a0wrELbdT5aS06yxkoRQloWFEKEmIxNz6SSROgUlo4qWDrZBCEhpGVBIUSIySR9sLROtldEwdLSzyTbB0JIy4JCiBCT8RmEkksBaC04SvETWbA0m64SQhIfCiFCTEYv6yoZCGbBYR0hQkhzhUKIEItIttgYm01tCVIXVIzENZacgpAQ0rKgECLEZJLdIuS7H61rzH8cQghJNCiECDGZZC0k6NdrLEqLkBQsbuOnDCEkgeFHFCEmo5d+ngz4gqWln9F1n2ewNCEkGaAQIsRkfLE2cV5ImMgxPfJ95WORjNf0M8n2gRDSsqAQIsRk9JqWJgP2pk8D/WDp8MdLttdPCGmZUAgRYjZJ2lpiVPd8DOqYh4sGtwegFjKRuLdYWZoQkgw44r0AQpobkkss2SwiHVtn4IsZZ+g+Fombj64xQkgyQIuQDiUlJejbty+GDx8e76WQJCQnzen9mZ783zNsUagZNl0lhCQDFEI6zJgxA1u2bMGqVavivRSShIzv2w5/mzIA90zoHe+lRI0QhZiha4wQkgwk/1dWQhKMNKcdV47oFO9lmEI07i26xgghyQAtQoSQgERVHVpWQFRChJDEhUKIEBIQSQBFEiztK8xo4oIIIcRkKIQIIYGJooM8m64SQpIBCiFCSECiKQXApquEkGSAQogQEpBoRAybrhJCkgF+RBFCAiJE4xpj01VCSBJAIUQICYgkYaKqBUQdRAhJYCiECCEBiaagoiSeWFCREJLIUAgRQgISnWus6ad5yyGEENOhECKEBMRXEjGCrDHpJ5UQISSBoRAihAREdo1FVFCRTVcJIYkPhRAhJCBCFHWERvVog675mZjQr9DkVRFCiHmw6SohJCDRdAvrV5yLBfeONXE1hBBiPrQIEUICYovCNUYIIckAhRAhJCDM/CKENHcohAghQZAsQpRChJDmCYUQISQgkv6xUQcRQpopFEKEkIDI+ocWIUJIM4VCiBASEMYIEUKaOxRChJCAMGuMENLcoRAihATElO7zhBCSwFAIEUICEk33eUIISQYohAghIaFBiBDSXKEQIoQExBcsTSVECGmeUAgRQgJiY9oYIaSZQyFECAmIr6AilRAhpHlCIUQICUg03ecJISQZoBAihAREYB0hQkgzh0KIEBIQ2SJEIUQIaaZQCBFCAsKsMUJIc4dCiBASELrGCCHNHQohQkhAfK4xKiFCSPOEQogQEhCWESKENHcohAghAZFig2gQIoQ0VyiECCEBYbA0IaS5QyFECAmIFBtkow4ihDRTKIQIIQFhHSFCSHOHQogQEhCfAKISIoQ0TyiECCEBkWOEqIMIIc0UCiFCSEDkrLE4r4MQQqyCQogQEhCbbBGiFCKENE8c8V4AISRxmTSgCFX1LpzaMS/eSyGEEEugECKEBOTms7rj5rO6x3sZhBBiGXSNEUIIIaTFQiFECCGEkBYLhRAhhBBCWiwUQoQQQghpsTR7IbR3716MHTsWffv2xcCBA/HJJ5/Ee0mEEEIISRCafdaYw+HAs88+i1NPPRWHDx/GkCFDMGnSJGRmZsZ7aYQQQgiJM81eCBUVFaGoqAgA0LZtW7Ru3RrHjh2jECKEEEJI/F1jixYtwuTJk1FcXAxBEPD555/7nTNr1ix07doVaWlpGDp0KBYvXhzRXKtXr4bH40HHjh2jXDUhhBBCmgNxF0LV1dUYNGgQnn/+ed3HP/roI9x555146KGHsG7dOowZMwYTJ05EaWmpfM7QoUPRv39/v38HDhyQzykvL8fvfvc7vPzyy5a/JkIIIYQkB3F3jU2cOBETJ04M+PjTTz+NG264ATfeeCMA4Nlnn8W8efPwwgsv4PHHHwcArFmzJugc9fX1mDJlCh544AGMGjUq6Hn19fXy/YqKCgBAY2MjGhsbDb8mI0jjmT1uc4Z7Fh7cr/DhnoUP9yw8uF/hE8mehXOuIIqiGPaqLEIQBMyZMwcXXXQRAKChoQEZGRn45JNPMGXKFPm8P/zhD1i/fj1+/PHHkGOKoogrr7wSvXv3xsyZM4OeO3PmTDz66KN+x99//31kZGSE9VoIIYQQEh9qampw5ZVX4uTJk8jJyQl6btwtQsE4evQo3G432rVrpzrerl07HDx40NAYP/30Ez766CMMHDhQjj965513MGDAAL9zH3jgAdx9993y/YqKCnTs2BETJkwIuZHh0tjYiPnz52P8+PFwOp2mjt1c4Z6FB/crfLhn4cM9Cw/uV/hEsmeSR8cICS2EJARBUN0XRdHvWCBGjx4Nj8dj6NzU1FSkpqb6HXc6nZa9Ya0cu7nCPQsP7lf4cM/Ch3sWHtyv8Alnz8LZ27gHSwcjPz8fdrvdz/pz+PBhPysRIYQQQki4JLRFKCUlBUOHDsX8+fNVMULz58/HhRdeaPn8UvhUOCY2ozQ2NqKmpgYVFRX8VmAQ7ll4cL/Ch3sWPtyz8OB+hU8keyZdt42EQcddCFVVVeHXX3+V7+/atQvr169H69at0alTJ9x999245pprMGzYMIwcORIvv/wySktLccstt1i+tsrKSgBg3SFCCCEkCamsrERubm7Qc+KeNbZw4UKMGzfO7/j06dPx5ptvAvAWVHzyySdRVlaG/v3745lnnsGZZ55p+do8Hg8OHDiA7OxswzFJRpECsffu3Wt6IHZzhXsWHtyv8OGehQ/3LDy4X+ETyZ6JoojKykoUFxfDZgseBRR3IdRSqaioQG5urqHUPuKFexYe3K/w4Z6FD/csPLhf4WP1niV0sDQhhBBCiJVQCBFCCCGkxUIhFCdSU1PxyCOP6NYtIvpwz8KD+xU+3LPw4Z6FB/crfKzeM8YIEUIIIaTFQosQIYQQQlosFEKEEEIIabFQCBFCCCGkxUIhRAghhJAWC4VQHJg1axa6du2KtLQ0DB06FIsXL473khKGmTNnQhAE1b/CwkL5cVEUMXPmTBQXFyM9PR1jx47F5s2b47ji2LNo0SJMnjwZxcXFEAQBn3/+uepxI3tUX1+PO+64A/n5+cjMzMRvf/tb7Nu3L4avInaE2q9rr73W7z13+umnq85pSfv1+OOPY/jw4cjOzkbbtm1x0UUXYdu2bapz+B5TY2TP+D5T88ILL2DgwIHIyclBTk4ORo4ciW+++UZ+PJbvMQqhGPPRRx/hzjvvxEMPPYR169ZhzJgxmDhxIkpLS+O9tIShX79+KCsrk/9t3LhRfuzJJ5/E008/jeeffx6rVq1CYWEhxo8fL/eFawlUV1dj0KBBeP7553UfN7JHd955J+bMmYMPP/wQS5YsQVVVFS644AK43e5YvYyYEWq/AOC8885Tvee+/vpr1eMtab9+/PFHzJgxA8uXL8f8+fPhcrkwYcIEVFdXy+fwPabGyJ4BfJ8p6dChA5544gmsXr0aq1evxtlnn40LL7xQFjsxfY+JJKacdtpp4i233KI6dsopp4j3339/nFaUWDzyyCPioEGDdB/zeDxiYWGh+MQTT8jH6urqxNzcXPHFF1+M0QoTCwDinDlz5PtG9ujEiROi0+kUP/zwQ/mc/fv3izabTZw7d27M1h4PtPsliqI4ffp08cILLwz4nJa8X6IoiocPHxYBiD/++KMoinyPGUG7Z6LI95kRWrVqJb766qsxf4/RIhRDGhoasGbNGkyYMEF1fMKECVi6dGmcVpV4/PLLLyguLkbXrl1x+eWXY+fOnQCAXbt24eDBg6r9S01NxVlnncX9a8LIHq1ZswaNjY2qc4qLi9G/f/8Wu48LFy5E27Zt0atXL9x00004fPiw/FhL36+TJ08CAFq3bg2A7zEjaPdMgu8zfdxuNz788ENUV1dj5MiRMX+PUQjFkKNHj8LtdqNdu3aq4+3atcPBgwfjtKrEYsSIEXj77bcxb948vPLKKzh48CBGjRqF8vJyeY+4f4ExskcHDx5ESkoKWrVqFfCclsTEiRPx3nvv4YcffsA//vEPrFq1CmeffTbq6+sBtOz9EkURd999N0aPHo3+/fsD4HssFHp7BvB9psfGjRuRlZWF1NRU3HLLLZgzZw769u0b8/eYI4rXQCJEEATVfVEU/Y61VCZOnCjfHjBgAEaOHInu3bvjrbfekgMLuX+hiWSPWuo+XnbZZfLt/v37Y9iwYejcuTP+85//4OKLLw74vJawX7fffjs2bNiAJUuW+D3G95g+gfaM7zN/evfujfXr1+PEiRP47LPPMH36dPz444/y47F6j9EiFEPy8/Nht9v91Orhw4f9lC/xkpmZiQEDBuCXX36Rs8e4f4ExskeFhYVoaGjA8ePHA57TkikqKkLnzp3xyy+/AGi5+3XHHXfgyy+/xIIFC9ChQwf5ON9jgQm0Z3rwfQakpKSgR48eGDZsGB5//HEMGjQIzz33XMzfYxRCMSQlJQVDhw7F/PnzVcfnz5+PUaNGxWlViU19fT22bt2KoqIidO3aFYWFhar9a2howI8//sj9a8LIHg0dOhROp1N1TllZGTZt2sR9BFBeXo69e/eiqKgIQMvbL1EUcfvtt2P27Nn44Ycf0LVrV9XjfI/5E2rP9Gjp7zM9RFFEfX197N9jEQZ3kwj58MMPRafTKb722mvili1bxDvvvFPMzMwUd+/eHe+lJQT33HOPuHDhQnHnzp3i8uXLxQsuuEDMzs6W9+eJJ54Qc3NzxdmzZ4sbN24Ur7jiCrGoqEisqKiI88pjR2Vlpbhu3Tpx3bp1IgDx6aefFtetWyfu2bNHFEVje3TLLbeIHTp0EL/77jtx7dq14tlnny0OGjRIdLlc8XpZlhFsvyorK8V77rlHXLp0qbhr1y5xwYIF4siRI8X27du32P269dZbxdzcXHHhwoViWVmZ/K+mpkY+h+8xNaH2jO8zfx544AFx0aJF4q5du8QNGzaIDz74oGiz2cRvv/1WFMXYvscohOJASUmJ2LlzZzElJUUcMmSIKsWypXPZZZeJRUVFotPpFIuLi8WLL75Y3Lx5s/y4x+MRH3nkEbGwsFBMTU0VzzzzTHHjxo1xXHHsWbBggQjA79/06dNFUTS2R7W1teLtt98utm7dWkxPTxcvuOACsbS0NA6vxnqC7VdNTY04YcIEsaCgQHQ6nWKnTp3E6dOn++1FS9ovvb0CIL7xxhvyOXyPqQm1Z3yf+XP99dfL18GCggLxnHPOkUWQKMb2PSaIoiiGZ0MihBBCCGkeMEaIEEIIIS0WCiFCCCGEtFgohAghhBDSYqEQIoQQQkiLhUKIEEIIIS0WCiFCCCGEtFgohAghhBDSYqEQIoQkHbt374YgCFi/fr3lc7355pvIy8uzfB5CSHygECKEmMq1114LQRD8/p133nnxXlpIunTpgmeffVZ17LLLLsP27dtjuo5Vq1ahuLgYAHDgwAGkp6ejoaEhpmsgpKXgiPcCCCHNj/POOw9vvPGG6lhqamqcVhMd6enpSE9Pj+mcy5YtwxlnnAEAWLx4MYYNG4aUlJSYroGQlgItQoQQ00lNTUVhYaHqX6tWrQAAV1xxBS6//HLV+Y2NjcjPz5fF09y5czF69Gjk5eWhTZs2uOCCC7Bjx46A8+m5rz7//HMIgiDf37FjBy688EK0a9cOWVlZGD58OL777jv58bFjx2LPnj246667ZCtWoLFfeOEFdO/eHSkpKejduzfeeecd1eOCIODVV1/FlClTkJGRgZ49e+LLL780tnkAli5dKguhJUuWyLcJIeZDIUQIiSlXXXUVvvzyS1RVVcnH5s2bh+rqalxyySUAgOrqatx9991YtWoVvv/+e9hsNkyZMgUejyfieauqqjBp0iR89913WLduHc4991xMnjwZpaWlAIDZs2ejQ4cO+Mtf/oKysjKUlZXpjjNnzhz84Q9/wD333INNmzbh5ptvxnXXXYcFCxaoznv00Ucxbdo0bNiwAZMmTcJVV12FY8eOBVzfkiVLkJeXh7y8PHz66ad46KGHkJeXhxdffBH//Oc/kZeXhyeeeCLi108ICUCUDWQJIUTF9OnTRbvdLmZmZqr+/eUvfxFFURQbGhrE/Px88e2335afc8UVV4hTp04NOObhw4dFAHL36V27dokAxHXr1omiKIpvvPGGmJubq3rOnDlzxFAfcX379hX/9a9/yfc7d+4sPvPMM6pztGOPGjVKvOmmm1TnTJ06VZw0aZJ8H4D45z//Wb5fVVUlCoIgfvPNNwHXUltbK+7atUv85ptvxFatWok7d+4UV69eLaakpIhbt24Vd+3aJR4/fjzo6yGEhA8tQoQQ0xk3bhzWr1+v+jdjxgwAgNPpxNSpU/Hee+8B8Fp/vvjiC1x11VXy83fs2IErr7wS3bp1Q05ODrp27QoAsvUmEqqrq3Hfffehb9++yMvLQ1ZWFn7++eewx9y6daufq+qMM87A1q1bVccGDhwo387MzER2djYOHz4ccNy0tDR06dIFGzZswMSJE9G1a1f8/PPPGDNmDE455RR06dKF2WuEWACDpQkhppOZmYkePXoEfPyqq67CWWedhcOHD2P+/PlIS0vDxIkT5ccnT56Mjh074pVXXkFxcTE8Hg/69+8fMHPKZrNBFEXVscbGRtX9P/7xj5g3bx7+/ve/o0ePHkhPT8ell14aUTaWMvYIAERR9DvmdDr9nhPMtZeVlQUAqK+vh81mwxdffIGGhgaIooisrCyMGTMG33zzTdhrJYQEh0KIEBJzRo0ahY4dO+Kjjz7CN998g6lTp8pZUeXl5di6dSteeukljBkzBoA3fiYYBQUFqKysRHV1NTIzMwHAr8bQ4sWLce2112LKlCkAvDFDu3fvVp2TkpICt9sddK4+ffpgyZIl+N3vficfW7p0Kfr06RPydQdj/fr1cLlcOPXUU/Hdd9+hsLAQY8aMwaxZszBgwICYZ64R0lKgECKEmE59fT0OHjyoOuZwOJCfnw/Aax258sor8eKLL2L79u2qQONWrVqhTZs2ePnll1FUVITS0lLcf//9QecbMWIEMjIy8OCDD+KOO+7AypUr8eabb6rO6dGjB2bPno3JkydDEAQ8/PDDfhaaLl26YNGiRbj88suRmpoqr1fJH//4R0ybNg1DhgzBOeecg6+++gqzZ89WZaBFQo8ePbB8+XK0a9cOo0ePRmlpKSorK3HBBRf4WZcIIebBGCFCiOnMnTsXRUVFqn+jR49WnXPVVVdhy5YtaN++vSrmxmaz4cMPP8SaNWvQv39/3HXXXXjqqaeCzte6dWu8++67+PrrrzFgwAB88MEHmDlzpuqcZ555Bq1atcKoUaMwefJknHvuuRgyZIjqnL/85S/YvXs3unfvjoKCAt25LrroIjz33HN46qmn0K9fP7z00kt44403MHbsWOMbFICFCxfizDPPBAD8+OOPGDlyJEUQIRYjiFrHOiGEEEJIC4EWIUIIIYS0WCiECCGEENJioRAihBBCSIuFQogQQgghLRYKIUIIIYS0WCiECCGEENJioRAihBBCSIuFQogQQgghLRYKIUIIIYS0WCiECCGEENJioRAihBBCSIuFQogQQgghLZb/D4CdJtwwq/V0AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/nmajik/Documents/SLAC/FACET2-Bmad-PyTao/UTILITY_quickstart.py:407: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  popt, pcov = curve_fit(smallestIntervalImpliedEmittanceModelFunction, zValues, sigmaXResults, p0=[1, 1, 0])\n"
     ]
    }
   ],
   "source": [
    "from scipy.optimize import minimize\n",
    "\n",
    "#optimizationType = \"differentialEvolution\" #Classic, standard for exploration\n",
    "optimizationType = \"nelderMead\"            #Classic, standard for refinement\n",
    "#optimizationType = \"L-BFGS-B\"              #This never seems to work; I think the problem is too noisy to get a reasonable Hessian\n",
    "#optimizationType = \"powell\"                #This rarely seems to help; I think the orthogonal approach is not appropriate for this problem\n",
    "#optimizationType = \"COBYQA\"                #Seems promising. Not timid about making big, coordinated moves. I like it more for refinement though. Seems (much) more likely than nelderMead to give up too early though...\n",
    "#optimizationType = \"bayesian\"              #This uses a different library, different signature, different objective, different bounds, etc... beware!\n",
    "\n",
    "match optimizationType:\n",
    "\n",
    "    case \"differentialEvolution\":\n",
    "        totalPopSize = 10            # Total population size\n",
    "        numDimensions = len(bounds)  # Number of parameters\n",
    "        \n",
    "        #initializationType = \"uniform\"\n",
    "        initializationType = \"normal\"\n",
    "        \n",
    "        match initializationType:\n",
    "        \n",
    "            case \"uniform\":\n",
    "                populationDE = np.random.rand(totalPopSize, numDimensions)\n",
    "                for i in range(numDimensions):\n",
    "                    low, high = bounds[i]\n",
    "                    populationDE[:, i] = low + populationDE[:, i] * (high - low)  \n",
    "                    \n",
    "            case \"normal\":\n",
    "                populationDE = np.zeros((totalPopSize, numDimensions))\n",
    "                for i in range(numDimensions):\n",
    "                    low, high = bounds[i]\n",
    "                \n",
    "                    #Define mean based on bounds\n",
    "                    #mean = (high + low) / 2\n",
    "                \n",
    "                    #Optional: Define mean based on defaultSettingsVector\n",
    "                    mean = defaultSettingsVector[i]\n",
    "                    \n",
    "                    std_dev = (high - low) / 50\n",
    "                    populationDE[:, i] = np.random.normal(mean, std_dev, totalPopSize)\n",
    "\n",
    "        #Optional: Add specific points to initial evaluation list\n",
    "        populationDE = np.vstack([[ defaultSettingsVector ], populationDE])\n",
    "        \n",
    "        result = differential_evolution(\n",
    "            wrapperDE, \n",
    "            bounds,\n",
    "            maxiter=500, \n",
    "            disp=True,\n",
    "            polish = False, \n",
    "            init = populationDE\n",
    "        )\n",
    "\n",
    "    case \"nelderMead\":\n",
    "        #Optional: Manually define starting simplex. Need to enable in minimize() options\n",
    "        numDimensions = len(defaultSettingsVector)\n",
    "        initial_simplex = np.tile(defaultSettingsVector, (numDimensions + 1, 1))\n",
    "        for i in range(1, numDimensions + 1):\n",
    "            initial_simplex[i][i - 1] += 1e-2 * (bounds[i-1][1] - bounds[i-1][0])\n",
    "        \n",
    "        result = minimize(\n",
    "            wrapperDE, \n",
    "            defaultSettingsVector,\n",
    "            method = \"Nelder-Mead\",\n",
    "            bounds = bounds,\n",
    "            options={'initial_simplex': initial_simplex, 'adaptive': True},\n",
    "            #options={'adaptive': True}\n",
    "        )\n",
    "\n",
    "    case \"L-BFGS-B\":\n",
    "        result = minimize(\n",
    "            wrapperDE, \n",
    "            defaultSettingsVector,\n",
    "            method = \"L-BFGS-B\",\n",
    "            bounds = bounds,\n",
    "            #options = {\"eps\": [1e-4*x for x in defaultSettingsVector]}\n",
    "            options = {\"eps\": [1e-3*(x[1]-x[0]) for x in pbounds.values()]}\n",
    "        )\n",
    "\n",
    "\n",
    "    case \"powell\":\n",
    "        result = minimize(\n",
    "            wrapperDE, \n",
    "            defaultSettingsVector,\n",
    "            method = \"powell\",\n",
    "            bounds = bounds,\n",
    "        )\n",
    "\n",
    "    case \"COBYQA\":\n",
    "        result = minimize(\n",
    "            wrapperDE, \n",
    "            defaultSettingsVector,\n",
    "            method = \"COBYQA\",\n",
    "            bounds = bounds,\n",
    "            options = {\n",
    "                \"scale\": True, \n",
    "                \"initial_tr_radius\": 1e-1, #1e-1\n",
    "                \"final_tr_radius\": 1e-8, \n",
    "                \"disp\": True } #\"*_tr_radius\" appears to constrain the scaled, not raw, variables (if scaling is enabled)\n",
    "        )\n",
    "\n",
    "    case \"bayesian\":\n",
    "        optimizer = bayes_opt.BayesianOptimization(\n",
    "                f=problem.optimizerWrapper,\n",
    "                pbounds=pbounds,\n",
    "                random_state=1,\n",
    "                allow_duplicate_points=True, #2024-04-26 it was whining about this\n",
    "                verbose = 0\n",
    "        )\n",
    "        \n",
    "        #Initial point(s) to check\n",
    "        optimizer.probe(\n",
    "            params={key: importedDefaultSettings[key] for key in pbounds.keys()},\n",
    "            lazy=True,\n",
    "        )\n",
    "        \n",
    "        #Refer to https://bayesian-optimization.github.io/BayesianOptimization/exploitation_vs_exploration.html\n",
    "        #and https://github.com/bayesian-optimization/BayesianOptimization/blob/master/bayes_opt/bayesian_optimization.py\n",
    "        boInitialKappa = 10 + 0*10.\n",
    "        boFinalKappa = 1e-6 + 0*0.1\n",
    "        boNumIter = 1000\n",
    "        boKappaDecay = (boFinalKappa / boInitialKappa)**(1/boNumIter)\n",
    "        \n",
    "        acquisition_function = bayes_opt.util.UtilityFunction(\n",
    "                                           kind='ucb',\n",
    "                                           kappa=boInitialKappa,         #Default 2.576\n",
    "                                           xi=0.0,              #Default 0\n",
    "                                           kappa_decay=boKappaDecay,       #Default 0\n",
    "                                           kappa_decay_delay=0  #Default 0                     \n",
    "                                        )\n",
    "        \n",
    "        optimizer.maximize(\n",
    "            init_points=100, #Initial, random points. Unwise to omit. \"Often, 2*dim to 5*dim random points are recommended before switching to model-driven exploration.\"\n",
    "            n_iter=boNumIter,\n",
    "            acquisition_function=acquisition_function\n",
    "        )\n",
    "\n",
    "print(problem.history) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "537d5d94-ee16-4a6d-af06-53139a015f2b",
   "metadata": {},
   "source": [
    "## Check out results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b50f42b1-5f6b-4a17-99ac-3902f40189bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "#problem.history.to_csv('output_data.csv', index=False)\n",
    "\n",
    "problem.history = problem.history.sort_values(by='maximizeMe', ascending=False)\n",
    "\n",
    "\n",
    "bestConfigData = problem.history.iloc[0]\n",
    "bestConfigDict = bestConfigData.to_dict()\n",
    "\n",
    "print( bestConfigData ) \n",
    "\n",
    "\n",
    "\n",
    "setLattice(tao, **bestConfigDict)\n",
    "\n",
    "\n",
    "trackBeam(tao, centerBC14 = True, centerBC20 = True, verbose = True)\n",
    "\n",
    "P = getBeamAtElement(tao, \"PENT\")\n",
    "PDrive, PWitness = getDriverAndWitness(P)\n",
    "\n",
    "print(f\"\"\"P, sigma x: {P[\"sigma_x\"]}\"\"\")\n",
    "print(f\"\"\"PDrive, sigma x: {PDrive[\"sigma_x\"]}\"\"\")\n",
    "print(f\"\"\"PWitness, sigma x: {PWitness[\"sigma_x\"]}\"\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64dec3b8-84cd-4210-955d-3ff84b7e60b6",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "bestConfigDict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "041feb20-cbd4-4875-a53a-27e1673b6381",
   "metadata": {},
   "outputs": [],
   "source": [
    "display(plotMod(P, 'x', 'y',  bins=300))\n",
    "display(plotMod(P, 'x', 'pz', bins=300))\n",
    "display(plotMod(P, 'x', 'xp', bins=300))\n",
    "display(plotMod(P, 'y', 'yp', bins=300))\n",
    "display(plotMod(P, 'delta_t', 'pz', bins=300))\n",
    "display(slicePlotMod(P, 'norm_emit_x',n_slice=100))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24e145e5-5de1-4039-8a0e-7615383bd0cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "display(plotMod(PWitness, 'x', 'y',  bins=300))\n",
    "display(plotMod(PWitness, 'x', 'pz', bins=300))\n",
    "display(plotMod(PWitness, 'x', 'xp', bins=300))\n",
    "display(plotMod(PWitness, 'y', 'yp', bins=300))\n",
    "display(plotMod(PWitness, 'delta_t', 'pz', bins=300))\n",
    "display(slicePlotMod(PWitness, 'norm_emit_x',n_slice=100))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e873d9e9-c62f-40e3-859e-798648f88cc2",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"\"\"P, sigma x: {P[\"sigma_x\"]}\"\"\")\n",
    "print(f\"\"\"PDrive, sigma x: {PDrive[\"sigma_x\"]}\"\"\")\n",
    "print(f\"\"\"PWitness, sigma x: {PWitness[\"sigma_x\"]}\"\"\")\n",
    "\n",
    "display(plotMod(P, 'x', 'y',  bins=300))\n",
    "display(plotMod(PDrive, 'x', 'y',  bins=300))\n",
    "display(plotMod(PWitness, 'x', 'y',  bins=300))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce9d5824-6cf6-4e11-8aae-00d366d69c77",
   "metadata": {},
   "outputs": [],
   "source": [
    "P = getBeamAtElement(tao, \"MFFF\")\n",
    "\n",
    "display(plotMod(P, 'x', 'xp', bins=300))\n",
    "display(slicePlotMod(P, 'norm_emit_x',n_slice=100))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e7c2188-44ab-4d6b-8d52-bbb3c47f758d",
   "metadata": {},
   "source": [
    "#### "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}