{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# Computing Bounds on the Local Density of States using General Constraint Descent\n",
    "\n",
    "This notebook demonstrates the use of general constraint descent (GCD) for evaluating dual bounds, using local density of states (LDOS) maximization as the example problem. For more background on LDOS optimization and bounds, see LDOS.ipynb. \n",
    "\n",
    "The tightest dual bound for a photonic inverse design problem is found by imposing all possible local power conservation constraints in the corresponding field optimization QCQP. Unfortunately, due to the large number of constraints, evaluating this tightest dual bound can be extremely expensive.\n",
    "\n",
    "The basic idea of GCD is to approximate the tightest dual bound with a smaller, more manageable number of constraints. This is pursued by iteratively adding new constraints to the QCQP to tighten the bounds, and merging old constraints to keep the total number of constraints fixed. \n",
    "\n",
    "GCD is implemented for all shared projection QCQPs in dolphindes, and this notebook demonstrates the high-level use of the available GCD functionality. For more mathematical details on GCD, see Appendix B of https://arxiv.org/abs/2504.10469."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.sparse as sp\n",
    "import matplotlib.pyplot as plt\n",
    "import time\n",
    "from dolphindes import photonics, geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x79e75a97e9b0>"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAGkCAYAAAABnUPEAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAGQpJREFUeJzt3Q2QVXX9+PHPKrCCLaSQPKQgzDDjA5oPKBMyso2Ck2Q6TD6hhmPTQPiElDyEJTIJQUVMkRpMYxYxOE2a1tgE+YAyVCKGD9hgjoSbypDFAAqBwPnP98x/d9gVNfst3+Xufb1mDrv33LPL2TN373u/537vvTVFURQBAJkclus/AgDhASA7Ix4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALKq6PDcdddd0b9//zjiiCPizDPPjKeeeiqqzezZs+Oss86Kurq6OOaYY+KSSy6J9evXN9smvSrSjBkzok+fPtG5c+eor6+PdevWRTUeq5qampg4cWLTumo/Nq+//npcffXV0b179+jSpUucdtppsWbNmqj247Nnz5647bbbyvuX9HMPGDAgZs6cGfv27YtqPzatoqhQS5cuLTp27FgsWrSoeOmll4qbb765OPLII4uNGzcW1eSCCy4o7r333uLFF18s1q5dW4waNaro27dv8fbbbzdt8+1vf7uoq6srfvWrXxUvvPBCcfnllxe9e/cutm3bVlSLp59+ujj++OOLU089tbytNKrmY/Pvf/+76NevX3HttdcWf/7zn4sNGzYUf/jDH4pXXnmlqPbj861vfavo3r178dvf/rY8Lr/85S+Lj33sY8X8+fOLaj82raFiw3P22WcX48ePb7buhBNOKKZOnVpUs82bN6cXfS1WrFhRXt63b1/Rq1ev8pek0X/+85+iW7duxT333FNUg+3btxcDBw4sli9fXgwfPrwpPNV+bKZMmVIMGzbsfa+v5uOT/oC77rrrmq0bPXp0cfXVVxfVfmxaQ0Weatu9e3d5OmDkyJHN1qfLq1atimq2devW8uPRRx9dftywYUNs2rSp2bGqra2N4cOHV82xuv7662PUqFFx/vnnN1tf7cfm4YcfjsGDB8ell15anqY9/fTTY9GiRU3XV/PxGTZsWDz66KPx8ssvl5efe+65WLlyZVx44YVR7cemNXSICvTWW2/F3r17o2fPns3Wp8vpxlCt0gh20qRJ5S/NoEGDynWNx+NAx2rjxo3R3i1dujSeffbZWL169Xuuq/Zj8+qrr8bdd99d3ma+/vWvx9NPPx033XRTeQf6xS9+saqPz5QpU8o/4k444YQ4/PDDy/ubO++8M6688sry+mo+NlUbnkbpgeKWd7wt11WTG264IZ5//vnyL7OWqvFYNTQ0xM033xzLli0rJ6C8n2o8Nkl6oDyNeGbNmlVeTiOe9OB4ilEKTzUfn/vvvz8WL14cS5YsiZNPPjnWrl1bTkpJEwnGjh1b1cemNVTkqbYePXqUf4W0HN1s3rz5PX+BVIsbb7yxPHXy+OOPx7HHHtu0vlevXuXHajxW6XRs+jnTjMcOHTqUy4oVK+IHP/hB+Xnjz1+Nxybp3bt3nHTSSc3WnXjiifHaa69Ftd92br311pg6dWpcccUVccopp8Q111wTt9xySzkzstqPTdWGp1OnTuWdyfLly5utT5eHDh0a1ST9hZVGOg888EA89thj5fTP/aXL6Zdk/2OVHiNLd8Dt/Vidd9558cILL5R/rTYu6S/8q666qvw8TZGt1mOTnHPOOe+Zep8e0+jXr19U+21nx44dcdhhze8e0x+7jdOpq/nYtIqiwqdT/+QnPymnU0+cOLGcTv33v/+9qCZf+cpXypk0TzzxRPHmm282LTt27GjaJs28Sds88MAD5bTPK6+8smqnfe4/q63aj02aYt6hQ4fizjvvLP72t78Vv/jFL4ouXboUixcvLqr9+IwdO7b45Cc/2TSdOv38PXr0KCZPnlxU+7FpDRUbnuRHP/pR+TyETp06FWeccUbTFOJqkv52ONCSntvTKE39vP3228vpn7W1tcW5555b/qJUo5bhqfZj85vf/KYYNGhQ+bOnpyMsXLiw2fXVenxSPNLtJD0n7ogjjigGDBhQTJ8+vdi1a1dR7cemNdSkf1pn7AQA7fQxHgAql/AAkJXwAJCV8ACQlfAAkJXwAJBVRYdn165d5RsxpY84Pm47frfc71SGin4ez7Zt26Jbt27lq8h27dq1rXfnkOP4ODZuO36vDkUVPeIBoPK0aXjuuuuu8sX20kvWpxf9fOqpp9pydwBoz+/Hk97vIr2/RYpPepXcH//4x/HZz342Xnrppejbt+8Hfm16hdg33nijfGXmxlNKvFfjcXF8HJuPym3Hsfmo0v3x9u3by/csavnK3gfauE2cffbZxfjx45utSy9SOHXq1A/92oaGhvd9cUyLY+A24DbgNhBtdgzS/fOHaZMRT3rfivQmXemNlvaX3r/8QO9Xnmat7T9zrXGkMywujA7RMcMeA/BB9sS7sTIeibq6uvgwbRKet956q3wP8wO9X3nLd/RL0rv+3XHHHe9Zn6LToUZ4ANrc/58f/d+89XebTi74b9+vfNq0aeWU6caloaEh414C0JraZMTTo0eP8m1k/9v3K6+trS0XACpfm4x4OnXqVE6f3v/9ypN02fuVA7RvbTadetKkSXHNNdfE4MGD49Of/nQsXLgwXnvttRg/fnxb7RIA7Tk8l19+efzrX/+KmTNnxptvvhmDBg2KRx55JPr169dWuwRABhX5Wm2Nr0FWHxeb1QZwCNhTvBtPxEP/1Wtneq02ALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeACo7PLNnz46zzjor6urq4phjjolLLrkk1q9f32yboihixowZ0adPn+jcuXPU19fHunXrWntXAKiG8KxYsSKuv/76+NOf/hTLly+PPXv2xMiRI+Odd95p2mbu3Lkxb968WLBgQaxevTp69eoVI0aMiO3bt7f27gBwiKkp0vDjIPrnP/9ZjnxSkM4999xytJNGOhMnTowpU6aU2+zatSt69uwZc+bMiXHjxr3ne6Tr09Jo27Ztcdxxx0V9XBwdajoezN0H4L+wp3g3noiHYuvWrdG1a9e2fYwn7URy9NFHlx83bNgQmzZtKkdBjWpra2P48OGxatWq9z19161bt6YlRQeAynRQw5NGN5MmTYphw4bFoEGDynUpOkka4ewvXW68rqVp06aVAWtcGhoaDuZuA3AQdTiY3/yGG26I559/PlauXPme62pqat4TqZbr9h8RpQWAynfQRjw33nhjPPzww/H444/Hscce27Q+TSRIWo5uNm/e/J5READtT6uHJ41c0kjngQceiMceeyz69+/f7Pp0OcUnzXhrtHv37nLywdChQ1t7dwBo76fa0lTqJUuWxEMPPVQ+l6dxZJMmBaTn7KTTaWlG26xZs2LgwIHlkj7v0qVLjBkzprV3B4D2Hp677767/JieFLq/e++9N6699try88mTJ8fOnTtjwoQJsWXLlhgyZEgsW7asDBUA7dtBfx7PwZCex5NGUJ7HA3BoOKSexwMA+xMeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgCyEh4AshIeALISHgDaV3hmz54dNTU1MXHixKZ1RVHEjBkzok+fPtG5c+eor6+PdevWHexdAaC9h2f16tWxcOHCOPXUU5utnzt3bsybNy8WLFhQbtOrV68YMWJEbN++/WDuDgDtOTxvv/12XHXVVbFo0aI46qijmo125s+fH9OnT4/Ro0fHoEGD4r777osdO3bEkiVLDtbuwCHr92+sbbZAe3fQwnP99dfHqFGj4vzzz2+2fsOGDbFp06YYOXJk07ra2toYPnx4rFq16oDfa9euXbFt27ZmCwCVqcPB+KZLly6NZ599tjyN1lKKTtKzZ89m69PljRs3vu/jRHfcccfB2FUAKn3E09DQEDfffHMsXrw4jjjiiPfdLk042F86BddyXaNp06bF1q1bm5b0fwBQmVp9xLNmzZrYvHlznHnmmU3r9u7dG08++WQ5mWD9+vVNI5/evXs3bZO+puUoaP9TcWmB9uiCPqe19S5AZY94zjvvvHjhhRdi7dq1TcvgwYPLiQbp8wEDBpSz2JYvX970Nbt3744VK1bE0KFDW3t3AGjvI566urpyptr+jjzyyOjevXvT+vScnlmzZsXAgQPLJX3epUuXGDNmTGvvDgDVMLngw0yePDl27twZEyZMiC1btsSQIUNi2bJlZbQAaN9qivSofoVJ06m7desW9XFxdKjp2Na7A1D19hTvxhPxUDkBrGvXrh94PLxWGwBZCQ8AWQkPAFkJDwBZCQ8A7X86NbQlrwD9wbySAgebEQ8AWQkPAFkJDwBZCQ8AWQkPAFkJDwBZCQ8AWQkPAFkJDwBZCQ8AWQkPAFkJDwBZCQ8AWQkPAFkJDwBZCQ8AWQkPAFkJDwBZCQ8AWQkPAFkJDwBZCQ8AWQkPAFkJDwBZCQ8AWQkPAFkJDwBZCQ8AWQkPAFkJDwBZCQ8AWQkPAFkJDwBZCQ8AWQkPAFkJDwBZCQ8AWQkPAFkJDwBZCQ8AWQkPAFkJDwBZCQ8AWQkPAFkJDwBZCQ8AWQkPAFkJDwBZCQ8AWQkPAFkJDwBZCQ8AWQkPAFkJDwBZCQ8AlR+e119/Pa6++uro3r17dOnSJU477bRYs2ZN0/VFUcSMGTOiT58+0blz56ivr49169YdjF0BoL2HZ8uWLXHOOedEx44d43e/+1289NJL8b3vfS8+/vGPN20zd+7cmDdvXixYsCBWr14dvXr1ihEjRsT27dtbe3cAOMR0aO1vOGfOnDjuuOPi3nvvbVp3/PHHNxvtzJ8/P6ZPnx6jR48u1913333Rs2fPWLJkSYwbN661dwmA9jziefjhh2Pw4MFx6aWXxjHHHBOnn356LFq0qOn6DRs2xKZNm2LkyJFN62pra2P48OGxatWqA37PXbt2xbZt25otAFSmVg/Pq6++GnfffXcMHDgwfv/738f48ePjpptuip/97Gfl9Sk6SRrh7C9dbryupdmzZ0e3bt2aljSiAqAytXp49u3bF2eccUbMmjWrHO2kU2df/vKXyxjtr6amptnldAqu5bpG06ZNi61btzYtDQ0Nrb3bAFRqeHr37h0nnXRSs3UnnnhivPbaa+XnaSJB0nJ0s3nz5veMgvY/Fde1a9dmCwCVqdXDk2a0rV+/vtm6l19+Ofr161d+3r9//zI+y5cvb7p+9+7dsWLFihg6dGhr7w4A7X1W2y233FIGJJ1qu+yyy+Lpp5+OhQsXlkuSTqdNnDixvD49DpSW9Hl6vs+YMWNae3cAaO/hOeuss+LBBx8sH5eZOXNmOcJJ06evuuqqpm0mT54cO3fujAkTJpTP+xkyZEgsW7Ys6urqWnt3ADjE1BTpUf0Kk6ZTp9lt9XFxdKjp2Na7Q4X5/Rtr23oXDmkX9DmtrXeBCrSneDeeiIfKCWAf9ji812oDICvhASAr4QEgK+EBICvhAaCyp1PDoc6sLWhbRjwAZCU8AGQlPABkJTwAZCU8AGQlPABkJTwAZCU8AGQlPABkJTwAZCU8AGQlPABkJTwAZCU8AGQlPABkJTwAZCU8AGQlPABkJTwAZCU8AGQlPABkJTwAZCU8AGQlPABkJTwAZCU8AGQlPABkJTwAZCU8AGQlPABkJTwAZCU8AGQlPABkJTwAZCU8AGQlPABkJTwAZCU8AGQlPABkJTwAZCU8AGQlPABkJTwAZCU8AGQlPABkJTwAZCU8AGQlPABkJTwAZCU8AGQlPABkJTwAZCU8AGQlPABUdnj27NkTt912W/Tv3z86d+4cAwYMiJkzZ8a+ffuatimKImbMmBF9+vQpt6mvr49169a19q4AUA3hmTNnTtxzzz2xYMGC+Otf/xpz586N73znO/HDH/6waZu0bt68eeU2q1evjl69esWIESNi+/btrb07ALT38Pzxj3+Miy++OEaNGhXHH398fOELX4iRI0fGM8880zTamT9/fkyfPj1Gjx4dgwYNivvuuy927NgRS5Ysae3dAaC9h2fYsGHx6KOPxssvv1xefu6552LlypVx4YUXlpc3bNgQmzZtKmPUqLa2NoYPHx6rVq064PfctWtXbNu2rdkCQGXq0NrfcMqUKbF169Y44YQT4vDDD4+9e/fGnXfeGVdeeWV5fYpO0rNnz2Zfly5v3LjxgN9z9uzZcccdd7T2rgLQHkY8999/fyxevLg8bfbss8+Wp9G++93vlh/3V1NT0+xyOgXXcl2jadOmlTFrXBoaGlp7twGo1BHPrbfeGlOnTo0rrriivHzKKaeUI5k0ahk7dmw5kaBx5NO7d++mr9u8efN7RkH7n4pLCwCVr9VHPGmSwGGHNf+26ZRb43TqNM06xWf58uVN1+/evTtWrFgRQ4cObe3dAaC9j3guuuii8jGdvn37xsknnxx/+ctfyqnT1113XXl9Op02ceLEmDVrVgwcOLBc0uddunSJMWPGtPbuANDew5Oer/ONb3wjJkyYUJ4+S08SHTduXHzzm99s2mby5Mmxc+fOcpstW7bEkCFDYtmyZVFXV9fauwPAIaamSI/qV5g0nbpbt25RHxdHh5qObb07AFVvT/FuPBEPlRPAunbt+oHHw2u1AZCV8ACQlfAAkJXwAJCV8ACQlfAAkJXwAJCV8ACQlfAAkJXwAJCV8ACQlfAAkJXwAJCV8ACQlfAAkJXwAJCV8ACQlfAAkJXwAJCV8ACQlfAAkJXwAJCV8ACQlfAAkJXwAJCV8ACQlfAAkJXwACA8ALRfRjwAZCU8AGQlPABkJTwAZCU8AGQlPABkJTwAZCU8AGQlPABkJTwAZCU8AGQlPABkJTwAZCU8AGQlPABkJTwAZCU8AGQlPABkJTwAZCU8AGQlPABkJTwAZCU8AGQlPABkJTwAZCU8AGQlPABkJTwAZCU8AGQlPABkJTwAZCU8ABza4XnyySfjoosuij59+kRNTU38+te/bnZ9URQxY8aM8vrOnTtHfX19rFu3rtk2u3btihtvvDF69OgRRx55ZHz+85+Pf/zjH//3nwaA9heed955Jz71qU/FggULDnj93LlzY968eeX1q1evjl69esWIESNi+/btTdtMnDgxHnzwwVi6dGmsXLky3n777fjc5z4Xe/fu/b/9NAAc8mqKNET5X7+4pqYMyCWXXFJeTt8qjXRSWKZMmdI0uunZs2fMmTMnxo0bF1u3bo1PfOIT8fOf/zwuv/zycps33ngjjjvuuHjkkUfiggsu+ND/d9u2bdGtW7eoj4ujQ03H/3X3AWgle4p344l4qLyP79q1a77HeDZs2BCbNm2KkSNHNq2rra2N4cOHx6pVq8rLa9asiXfffbfZNilWgwYNatqmpRSvFJv9FwAqU6uGJ0UnSSOc/aXLjdelj506dYqjjjrqfbdpafbs2eUIp3FJoyMAKtNBmdWWTsHtL52Ca7mupQ/aZtq0aeXwrXFpaGho1f0FoELDkyYSJC1HLps3b24aBaVtdu/eHVu2bHnfbVpKp+vSOcP9FwAqU6uGp3///mVYli9f3rQuRWbFihUxdOjQ8vKZZ54ZHTt2bLbNm2++GS+++GLTNgC0Xx0+6hekqc+vvPJKswkFa9eujaOPPjr69u1bzmibNWtWDBw4sFzS5126dIkxY8aU26fHaL70pS/FV7/61ejevXv5dV/72tfilFNOifPPP791fzoAKj88zzzzTHzmM59pujxp0qTy49ixY+OnP/1pTJ48OXbu3BkTJkwoT6cNGTIkli1bFnV1dU1f8/3vfz86dOgQl112WbnteeedV37t4Ycf3lo/FwDt8Xk8bcXzeAAOLW32PB4A+DDCA0BWwgNAVsIDwKE9q+1Q0DgfYk+8G1FxUyMA2p/y/ni/++d2F57Gt1hYGY+09a4A0OL+OT1fs91Np963b1/5Vgpp19OTVtNrt3kZnQNPO08vqOr4ODYflduOY/NRpfvjFJ30bgOHHXZY+xvxpB/q2GOPbXp7BK/f9sEcH8fmf+W249h8FB820mlkcgEAWQkPAFlVdHjS2yXcfvvt5UccH7cdv1vudypDRU4uAKByVfSIB4DKIzwAZCU8AGQlPABkJTwAZCU8AGQlPABkJTwARE7/D7JwyRQbFFHUAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 461.538x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# First, let's define the relevant parameters for the simulation. \n",
    "\n",
    "wavelength = 1.0 # Dolphindes uses dimensionless units. \n",
    "omega = 2 * np.pi / wavelength\n",
    "chi = 4+1e-4j # Design material \n",
    "px_per_length = 40 # pixels per length unit. If wavelength = 1.0, then this is pixels per wavelength.\n",
    "dl = 1/px_per_length \n",
    "Npmlsep = int(0.5 / dl) # gap between design region and PML. Not required to be defined, it is just convenient.\n",
    "Npmlx, Npmly = int(0.5 / dl), int(0.5 / dl) # PML size.\n",
    "Mx, My = int(0.5 / dl), int(0.5 / dl) # design mask size \n",
    "Dx = int(0.1 / dl) # distance from the design region to the source region.\n",
    "Nx, Ny = int(Npmlx*2 + Npmlsep*2 + Dx + Mx), int(Npmly*2 + Npmlsep*2 + My) # grid size. This includes the pml layer!\n",
    "\n",
    "cx, cy = Npmlx + Npmlsep, Ny//2\n",
    "\n",
    "ji = np.zeros((Nx, Ny), dtype=complex) # current density\n",
    "ji[cx, cy] = 1.0/dl/dl # a delta function source in 2D is approximated by amplitude 1/dl/dl so that integration int(ji)dxdy = 1.0. \n",
    "design_mask = np.zeros((Nx, Ny), dtype=bool) # design mask\n",
    "design_mask[Npmlx + Npmlsep + Dx: Npmlx + Npmlsep + Dx + Mx, Npmly + Npmlsep: Npmly + Npmlsep + My] = True # design mask\n",
    "ndof = np.sum(design_mask) # number of degrees of freedom in the design region\n",
    "\n",
    "chi_background = np.zeros((Nx, Ny), dtype=complex) # background material\n",
    "chi_background[:, :] = 0\n",
    "\n",
    "plt.matshow(design_mask + np.real(ji)*dl*dl) # visualize where the mask and the source are"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Photonics_TM_FDFD(omega=6.283185307179586, geometry=CartesianFDFDGeometry(Nx=104, Ny=100, Npmlx=20, Npmly=20, dx=0.025, dy=0.025, bloch_x=0.0, bloch_y=0.0), chi=(4+0.0001j), des_mask=True, ji=True, ei=False, chi_background=True, sparseQCQP=True)\n",
      "Vacuum LDOS:  0.7878298937193577\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAGhCAYAAAC3eRkzAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAQFNJREFUeJztnX2wVdV5xhdfXu5FQIHKR0XFGTpqMI1Ry5TYQEclk5q0DtN8+JGYpn9oMUZCG5SSVnQqVNtaJrUxo5OxtpbqdGJa2+kHpGlIHNpKMCaGdLSZUGKNlFEJoMBF4XT2nt7re17uefZe++x17z33/H4zZ7jrrn32Xnvtve9iPc+73zWu0Wg0AgAAQELGp9w5AAAAgw0AAAwLzGwAACA5DDYAAJAcBhsAAEgOgw0AACSHwQYAAJLDYAMAAMlhsAEAgOQw2AAAwNgebL7whS+EBQsWhMmTJ4eLL744fPOb3xzJ5gAAwFgbbB5//PGwatWqsG7duvDtb387/MIv/EJ4//vfH370ox+NVJMAACAR40YqEefixYvDu9/97vDAAw8M/u78888PV199ddi4caP87okTJ8KPf/zjMHXq1DBu3LhhaC0AAAxFNoQcOnQozJs3L4wf33r+MjGMAMeOHQs7d+4Mt99+e9Pvly9fHrZv337S9v39/flngJdeeilccMEFw9JWAAAo5sUXXwxnnnnm6BpsXnnllXD8+PEwe/bspt9n5b179560fTbTufPOO0/6/Yt33hmmTZ6ctK0AANCag0ePhvl33JErTYoRGWwG8BJYNh0bShZbu3ZtWL169WD54MGDYf78+flAM623d1jaCgAArSmyNEZksJk1a1aYMGHCSbOYffv2nTTbyejp6ck/AADQmYxINNopp5yShzpv3bq16fdZecmSJSPRJAAASMiIyWiZLPaxj30sXHLJJeHnf/7nw4MPPpiHPd90000j1SQAABhrg81HPvKR8Oqrr4a77rorvPzyy2HRokXhH/7hH8LZZ589Uk0CAIBEjGiAwMqVK/MPAACMbciNBgAAyWGwAQCA5DDYAABAchhsAAAgOQw2AACQHAYbAABIDoMNAAAkh8EGAACSw2ADAABjO4PAmGWkVw8dmcVXAbqHkX7GO/CZZ2YDAADJYbABAIDkIKMNxzR6woQwohw/nv4YHTCNhzHMcMhaI/0cq2danf8oeTaZ2QAAQHIYbAAAIDkMNgAAkBw8m7o0YavnjndjuC37uhScONFcnjixdV3MfurwhdrRj2PaA6OHqvd8Kl+0bHvUdu08x+PHt/88+uP7/QyHTxsJMxsAAEgOgw0AACQHGa0qftpupSo/xS1b1w5q2mylKz/dVtN4ta2vU9JA2Sl9kaRQVnKoS24bi7JdXTLuSOxHSdVqn0XlVnV+OyXrKRlPHe+EeFa85Gy3feut5jpfbrXPofY7TDCzAQCA5DDYAABAchhsAAAgOXg2CqvRek1WeS/ehznllKF/9vtpJx1GWV/Ga7vKh1HbqjqP8mzU8WM8ozra0o6ePdL+Torw4qo+RF2eiapT+/TPX4yHWnY/MX6OohHhyxw71vr4HvtdldpmGP0bZjYAAJAcBhsAAEgOgw0AACQHz6YsRe/H2PLkyc11tuzrlNcTkza8rPei4vOLYvdVLL+ti/FT1PHVtqou5r2CmPeO6kgBEqOR15U2P8Z7UV6A8jBj3iVTvogtx7Tbfs/7ouoYVetUv8Vct0aj9X1sPZqMo0eH/rnoPvb9OEKpbJjZAABAchhsAAAgOchoFj/9VekxVHizl8r6+ob+2X+vKGTTEhPe3N/fegptp+p2u3ayzr75ZnOdLXtpwNap7/m2x9TFpOtR6ULU98rW1UXV8GF1j8eE906a1Ho/qs7X+21t2cthrdpZ1G7/XNlj9PQ019lj+nbbbdsJi26I+1E9j1X/Hqh7XLWtZpjZAABAchhsAAAgOQw2AACQHDybqihdOMazsduqVDYx4c3eF7H78XXK6/HeR9kU5/57KmTTln3bvGZt9xvjS8Wk2bFlFUKdasXTFCtMxqRdKpuuJcbP8L6M3dbf8/Z5UP3mj2+va9GyHfa7MamlVJ0KmVacEJ7N4cOt7z+VyqbIQ7LnP4xh0MxsAAAgOQw2AACQHGS0ssSEPitpIEZG8ygZzU6ji7I3t8JPqdVU3cthdsrvp/+2rL7n65Ss5qU6u23VkOkiGa1sBoeRoOwb7EUyWtkQZl9n710lm/ltleTsr3/ZLBVFz2pZ6SgmS0FdrzBMFH+SlVTuv1d2VdNhZHS0AgAAxjQMNgAAkBwGGwAASA6ejaLsyoBF4aUqBYfyepS+6zXbsigfxocaK3/l9deb62zZezaqzh5D1fm2xvg5SutvJ7XNKMikG50huSiVir1XlWcRE84c42Haa+7rYjKEW9Sz688jJmy51T5VtmiParc/x7Kh1+2u+psIZjYAAJAcBhsAAEgOMlpV/NRcZbJVbyyrBZliQp/VW/JKflIylpfKDh5sXWfLMRKbLb/xRuu2FWWPtuWYBeJGW5aAkcwuEHOvxoT+++fBymj+GpcNfVYyWsy1UOeowotjjqFkvBMVZbQiqa6sBUAGAQAAGEsgowEAQHIYbAAAIDl4NimIWUUwJpWF0qyr+jKHDpXzaHy9ryvr58SksonJEG39HK9D27p20szU4cvEpA4Z7mzRflsVFh0T3u9T0qjVWdV9XDX0WeHPsWwIc0ydb9t48feg1XYx3xuldF6LAQCg42CwAQCA5CCjlaVo2hozHbbEhD6qrK9q0SWVkdmGGyv5y0tlXkY7cGDofXqp7siR6jJa2fDmmAXR6qLsG9upQqTVeVU9Z39Odj/+ni4blh4jlanF+9q5jkqOsnKgCieOkdFipPOqDMcx2mT0tQgAAMYcDDYAAJAcBhsAAEgOns1wUDbNTNEKm1VXyiybSqYo9NmWrUdTFEJd1jNSYbC+rNLO+PDmGMquchmjkat9xqDOS6VAiVlxtOwKr/6c7LZFob+2rHwZ/72qqYRiPBuVPkp5NmXT/LQTwkzoMwAAgAYZDQAAksNgAwAAycGzGW6KfBmLSjuj3klR79koz8Z7NNaH8T6Nr1Pej/VlVLuLPBuVdqYsKqV+XXp+0TGqonypst5f0eqjZb0vX2eP4T2iqu89Fe2njpQ06hor70WtlBmz+ufEiFU8R3oZi9E2s9m4cWO49NJLw9SpU8MZZ5wRrr766vD88883bdNoNML69evDvHnzQm9vb1i2bFnYtWtX3U0BAICxOths27Yt3HzzzeHf//3fw9atW8Nbb70Vli9fHt4w/7u99957w3333Rfuv//+sGPHjjBnzpxw5ZVXhkP+f8sAADAmqF1G+6d/+qem8sMPP5zPcHbu3Bne+9735rOaTZs2hXXr1oUVK1bk2zzyyCNh9uzZYfPmzeHGG28MHUnVVTRbbRcro1mpSqWLqSp/FR3DftfX2RQ1vt39/eWyNcdIZ0rG8rKJyl4cU6dCX1OEPqvQb9WPRX1cNoQ6RmLzx2i1T//dGKnU9rFKQVN0rWydusZKKvMSW9VUMieEbNaBklryAIED/6/zz5gxI/939+7dYe/evflsZ4Cenp6wdOnSsH379tTNAQCAsRYgkM1iVq9eHS677LKwaNGi/HfZQJORzWQsWXnPnj1D7qe/vz//DHDQG9kAANC9M5tPfepT4bvf/W74q7/6q5PqxjlpIRuY/O9s0MH06dMHP/Pnz0/WZgAA6KCZzS233BKefPLJ8I1vfCOceeaZg7/PggEGZjhz584d/P2+fftOmu0MsHbt2nyGZGc2wz7gtLOio9WlY1bfVGHCvs76Ij6Nvy0rr8V7Nt7fsSHUfj+2Pf74ts56NEV+QlXNXmn0vs7r67ZebRuTYt62LSYMWvWHWmJBpTny/okKN1d+Toyfpjycqt6DCllXXovfVl3Hnp7mOltWdf6e8sefKP7sKs+sA32apDObbIaSzWieeOKJ8LWvfS0sWLCgqT4rZwNOFqk2wLFjx/IotiVLlgy5z8zTmTZtWtMHAAC6eGaThT1nUWV/+7d/m79rM+DRZPJX9k5NJpWtWrUqbNiwISxcuDD/ZD/39fWFa6+9tu7mAADAWBxsHnjggfzf7EVNHwL9iU98Iv95zZo14ciRI2HlypVh//79YfHixWHLli354DQmUfKHnTZ7iUmFPsesxmnlsaqrePqyCr3252HLXn6pmq055u1uW/bSmJdDJk9uvR9VV/bN85gwaPWWfswqrvZaqbqia2W/68N3VdvUecWsuGn7TmUFULKZvwfUvVJXnb8fJk0qFxYdk727GwebTEYrIpvdZBkEsg8AAIx9SMQJAADJYbABAIDkkPW5KjGZbK2GrTR7r5GrEOKYlTqr1qnwarWt8mXaCW9WmXVV6Kn1Wvr6Wtdl9Pa2rrNl7/WU1ezr8mzUKqbeMysbMu/r/TVW2ZL9MVu1u+gc7f3hj6HS3qiURP5esecV472U9fN8H/ttjwnvy94rpKsBAACIAxkNAACSg4yWArWwlapTi4X5eiVjqRDmmEzO6hiqrXUteqZCRpXE5eumTHn751NPbS2b+W295GbL6vjqDfKisFwlx9pyTEbwmPB2K6vFZESOkQfLympKNvPY/vDnr87DS2wqg4S6x2y/+jol1U2MWDxNhYx3QLYBZjYAAJAcBhsAAEgOgw0AACQHz6YuVGZntVKnrfMadUzWZ+W9qIzMahVNVfaegT2vqitsKo+mSDO3qY6812J9Gp8SyZetZ6Pq/DGs96M0+qJVG9WKr7bPff+rsHR7P/hM3v48bL1Ku+OvTVViwqLtfeX7rWw4sW+7fx5sSLu/x2w/qjRH/rnx+zl2rJxn4xmFPkwMzGwAACA5DDYAAJAcZLQUqAXSPEp+iglvtW9wq7fEY6Q5f3x7DP/GuD2PurICqPBiJXFNn95cZ9c/Ou201nV+v34/SmJTEl/ZUFePWgRNXUcvlVn558CB5jq/rW27vx5VMyGkkNX8dvZ+9G3zsqbtKy+HKcnZ9k1MlgYlq01094OXAMs+Yx2QEZqZDQAAJIfBBgAAksNgAwAAycGzqYpKJRGzwl7ZTL4xunDVuiLPSGVvLqsZez3datYqnNT7JNY/8V6M92VU3YwZrbdV3o9Pe2PLI+HZ2NVY7c8ZBw+27ref/KS5bH0an5LGltvxbOyzorwH5f35+816Hb7O38f2eqR4joq8z2PHWvtJ9prHhMV3AMxsAAAgOQw2AACQHAYbAABIDp6NIkYjVduq9wNUChIfy6+2VXq+ej/GlouObz0ddY4eq+/7dzfKrozoU6t4P8V6LTNntvZlvEcza1br/fhtrWfj38+xno1PAaM8G6XLq3vFp6SxPo31aHxZpdkpSrHvPZyyPoxauVWldlJpj9RSHf57/j6210M9D+pdGvX8FT1H/f3lVnVV90YHwswGAACSw2ADAADJQUYbDsqmq4lZqTMmhFmFZaoVNlV4c0x6DLXCo5UR1AqbXqryMpaSv6xUdsYZzXVecrNlJbG5EOr+E5NaL3D6k3JJvzPGj3+7fyZObJaxJk9+W/Lqm9EsI/bMerN1OLMt+z724eZKxil7b/iT9GUloympVj1Hdp++3SoUWqWyUc9KTEZ0/zy++Wa5vvKSa1mpfqj6snUJYWYDAADJYbABAIDkMNgAAEBy8Gyq4nVgpacqjVTptyq8MiYs2h7Dt9se0x/fa832PMqGOnsN3YfT2rJaYdOnq1EpabxnY30a79n4svVpnGdzZPzbHtJPXimfLcbK+cWeTWvJ3totJ2fLebsfTzvtp5rqem2/qlVEi8KbVXixva98nfIl/Lb2fvD3qr2v/P2n0kWp9ig/SXmmKry5KO3UW+L4ZX2pGN9llIRMM7MBAIDkMNgAAEBykNFSoGQ0lRG6KOuyLXuJQU3jVXi1mtJXXQ1QZQv2so0NvfVhuf7tdiuj+dBnK535kGUbzuxlszlzmsumfv/rzZLfK6+0ji62L+l7Gc2GQrcjo1k1zMtotjv88WfNelv+O32eeGPdoyRXFc5btFJl2fBiJeOWlZuGOkd1z9vnQ9V5Gdsesyj0ub+/nMTn5WglXY8SqUzBzAYAAJLDYAMAAMlhsAEAgOTg2dSF1UzVKoLKz1EauS+r1DL+e1YjVvssMhSUZmzx6UJUJmHr03iPxqersWXr3xStxqlS0DgP59WDb7fv5ZebN3311bd/fu215jrr4SjPRiUyLkqQrTwbe8o+XY69/CfmNu90pvew1D1md3TkSHOdPegbb7SuK/Js7DF8B5T1DItSt5QNPfa+jK1TKXCKXiE4LtI+qdcLWKkTAABAg4wGAADJQUYbDspKZUVvE1cN2SybkVdN6Yuw+o+X0Wzos1o8TS2W5qUzv3iaLausz05G8+HNVjr73/9t3s2+fUNLal5GO3SofEJgFfqs1pLzKqJVrrwypaKEx49vPv/Tbf/4xtode2nMNsDLaKrsT9LeHz6E3nZOO2HBVZ+Hsq8QqLqYZ15l3S4iJiv7MMHMBgAAksNgAwAAyWGwAQCA5ODZVEVpxCrNi/JF2vFs1H7Kpsupmp7G6+s+PYjV4ZURUeTZ2NBnHxZt87WIVTxt5mafgsZ7Mdaj8WX/PevZ2NQ13t5ox7Ox3aHCm9Vl85fGL9Q5+cy3+6fXh5DbE/N9fOBA62vjr6O65vakvb9XdoXLmPQ1Mdnb7TGq1sVkiFc3R0xm+VECMxsAAEgOgw0AACSHwQYAAJKDZ1MXVbVWq0MrPyVmVc2YlRHVCocx2HNU79l408Bq9CcZCMLD8S+a2PwtwrPxK2z6pQJsGhr/Lo31abxnY783HJ5NzLs0tvuLbDHbjb2zTmvdWWoVVeXR+OvsT9LeH+o9m3beQVEpYao+R2VX/yzyd5Th1gG+jIKZDQAAJIfBBgAAkoOMNtLYabyfQqusr2pbFRZZlBFXHd/i04VYYmQ0W1Zh0T4rtK+zso6TePpPTGqZkdlLXlYp8hKbLfusz0pGs8dsR0azp6VkMxXerCLGfVv7ZzSHHveIPpbZu1V4s7oflIym7kd/38aksqn6HNnnsUgaOy5W51XHV3J8B8DMBgAAksNgAwAAyWGwAQCA5ODZlKUdjbTsCntFfkrVVfzUMeta/U+lf7c6fIxno8reB7BlF3prQ4+9Z6PKfqkA68V4X8aWvddTV+hzWZ9GhTfHnL9PidNjd6T6P+Y6qvvB16llLGIoe8/HpJIp6/XEeEgd6MsomNkAAEByGGwAACA5yGgjTYz8pcIrVVi02mcKfMiqlT9iJDa/Hyu/KGnG1R39Ses3771UZMtqW/89JT/VFfps8V1lVSx/fLXAptrWn384rVr/n3Qdy0pl/l7x+0mBej7qeP6K5PAxDDMbAABIDoMNAAAkh8EGAACSg2eTgrrSTMRsq3TimBU3q6LS11T1c1S6koiQWZWQVyXP9v6KLddVp24VTx3HV8mKffmkUOuyIewqXVFdvkzZ+60dqj5TMSmhuojkM5uNGzeGcePGhVWrVg3+rtFohPXr14d58+aF3t7esGzZsrBr167UTQEAgLE42OzYsSM8+OCD4Z3vfGfT7++9995w3333hfvvvz/fZs6cOeHKK68Mh/xbdAAAMCZINti8/vrr4brrrgsPPfRQOP3005tmNZs2bQrr1q0LK1asCIsWLQqPPPJIOHz4cNi8eXOq5owNsun5wKcbyOSYVp+I72UqRtmP7eKY73XKx55f7DlW7f+oa9ftnLAdXvF7o1S2S3blb7755nDVVVeFK664oun3u3fvDnv37g3Lly8f/F1PT09YunRp2L59e6rmAADAWAsQeOyxx8IzzzyTS2SebKDJmD17dtPvs/KePXuG3F9/f3/+GeCgT0wFAACjmtpnNi+++GK49dZbw6OPPhom+2R9hixowJLJa/53Nshg+vTpg5/58+fX3WwAAOikwWbnzp1h37594eKLLw4TJ07MP9u2bQuf//zn858HZjQDM5wBsu/42c4Aa9euDQcOHBj8ZANaV5INxgOfbkAaCOW/V2Qh2I/t4pjvdcrHnl/sOVbt/07wE0YN4yt6Wx3gi9Uuo11++eXhueeea/rdr/3ar4Xzzjsv3HbbbeHcc8/No8+2bt0aLrroorz+2LFj+YB0zz33DLnPzNPJPgAA0JnUPthMnTo1jzCzTJkyJcycOXPw99k7Nxs2bAgLFy7MP9nPfX194dprr627OQAA0K0ZBNasWROOHDkSVq5cGfbv3x8WL14ctmzZkg9UYwI1jY2Z4sZsq968Ho5suWXDsf2b1/Z7fh9+WyvBRLz6PnHiJPNz89d8edKkaomNy9Z5YrI+13F8e36xa5eVTsWgspXHXvNWDEf4f9Vnysvco1TWGpODzde//vWmchYIkGUQyD4AADD2YcgFAIDkMNgAAEByyPo80lg912u7Knuu2lbpycOhH8do9MoH8PupmPZ48uRe83Pz1/r6WpfVtv57p57a+jRaNTPWs7HHtMcraps9D3W+ftuTXpOrmnbaX0d1zevwc9pBPR91PH9+2255jYGZDQAADAfIaAAAkBwGGwAASA6eTVna8TqULxMTn6+03hjvp0rbirDmQ4wvY8vKB/Dlo0eb62z58OGmqr4Z01t6HarsX/uyu/WHVz6NfV+lHc/Gtm3atOY6W/Z19jxizt/7OeE10QG2HHMdY5YKtfdVO2lvyt7z6jnydfZ5LHr+xgmfpupz3AF0dusBAKAjYLABAIDkIKONNHZK7UMmlVSmtq0qv6m2FYWlWlS6khgZzUs1R460rnv99aF/zhK5znpz8OdTT23O1+Ilp9NOe/vnN95orjNLKpWWzXwIcV2hz77dM2YMfQ6+7Ov8fqyM1jP+7X4r6uOm62Gvk6+LkdFUuiKPuh/LylbtPEcqLNo/qxMmtM4fpI6v2t0BdF6LAQCg42CwAQCA5DDYAABAcvBs6qJqWKTVbIu0XpX/3e5H5Y33GrHV0+sKfY7xbKx+b00R3zYfe3zoUHOd9RAOHmyu+8lPBn887bSfavk1fwh/eJUtxXaxT/Ni95nKs7FezKxZzXUzZw7t7fjvnVQ2/XZSv/qOs9fDhZ6f1JH2OvsOKevZ1BX67L2Xqs+RejbVczy+4Jlv1e4OpLNbDwAAHQGDDQAAJAcZrSoqnFKFRaqpcdEykmqqrvajJD7btqJ2q3hfJZW9+Wa58GYlm/lYZB+XbCUeIaP1utfiZ82aUim82XdNT0/rN+/Vy/UxMpqV53x2Ayt/Wdks44wzWtd5ya33xBvlZDTfx+raKFlNhUXb+6ZIYlOo1wR8nXpWlPxVti4mS0DV7AajlNHfQgAA6HgYbAAAIDkMNgAAkBw8m+FA6anKh4nxbKzAr0Iv1T69fh2ji6uMvHY/XodXmZy91m/Daw8caK6bMqXc8pPOCDl9XnP5xFwX0loyvNkeXoVT+9P3WVasneCja9VKndaz8eHN1qeZO7e57vRTXYN+/MrbP79ifs547bXW/W/LPixdeTbexLIdpEKfVXoaj/IpY54HW+cNtbJ1Mc+8r4sJ9y5aLXQEYGYDAADJYbABAIDkIKPVhcrIXDaEUU3bi6bqShqwcbk+I6+atvvYX6vxxGR9ttKI15FUtmAfQmt1JC/V2DDd3t7W518gccw0ccLjx08qHd5s3+hXMprvUhX67C+HktHs8X1WABvefJJstm9fc9lKZ6++2rqPfVi0vR7+uvmyyt6t7pWyMlJRZvOyMpa94L7OP6tKxvZ66IQJ5f5WxCyQ2AF0XosBAKDjYLABAIDkMNgAAEBy8GxSEJOCQq3a5/0FW1b5UqqGZXqtWaWvKZu6xp+zD3W1+WHUypzeF/Cxx7bvVEbeopBQc16nu1wuk8+c0tIzKbuIZTuejT1lf3xb9p5NUwoaG9o8VHiz9XC8Z2O3tWHQPn2N99PUyp0+07e9P/x9VDYUv2jF26qvENg67+fYY/pnzN+PPT3V7lV1/h3g4Yz+FgIAQMfDYAMAAMlBRqtKTLZYNcVVizUpGc1P422dmsZXXZDNT+O9NGFDoX1YtNWKfDirLfs3zWMWr1JypOp/L01YGcdJfr1Gn+qd1axV9c+YVOqF+bpkNB963TP+zdZhybbsZTMvldmyD4u20pk6hlrYriilgi37zlHh9mUXKIxZ6EzJ2v4ZUzJ2jHQ9UdRVlc1GicQ2OloBAABjGgYbAABIDoMNAAAkB89mOFDaq9Wli0Imbb0P/bVl/z1b579nDQV/fB+WavVsr6eXDYX2+rH1RVRYeIzWrsKbi9qtQrFteK/zLHpM7HGPN1ROE+G0/pxt+3zbbF+9dri1L+JX0bRlH7LsvRcV3qw8G3sM79H4dDUq67O9V6quxumvv0otE/Os2G3V81f0HE8SPq1a8VP5wkUpelrVxfRxmzCzAQCA5DDYAABAchhsAAAgOXg2ipiU3mX1VOUneI0+5l0a5efY/fh92rJ/58Fvq/R0tfyALftjlE3l47Vur1H7cqt2q/c6fGoV/8KMzeOv8sX4/lcpUBTKs/F+ksqXY/0Uv8Km917Kvkvj6+x+1YtGBe8yNV2Psu/V+HtHpYcpel9NPSvKlyn7Dk5Mainl53UgzGwAACA5DDYAAJAcZLSqKKlMyUFKRilKc6HCMuuo86HOKhTaS1wqZNfipRG7rT++ksYU/hhK/vMympWAfMiulYqmTi3fx8Mho9myTxdjz8PLaH5bEd4t5Th7DL9P31Z7nf05KunMosKbVahzTAhzXXXquZ4oQp89rNQJAACgQUYDAIDkMNgAAEBy8GxS4D0blWZFpY6I8WxsOKWv6+1trZ8rz0aFQqvwZn8eSodX/k5VlC/kz8n3h/UevGczZUrrHP+2j5VGX7T8hPW+fKirChlWIdvKT/HnaOvVUgH+GHZb5Sf5axBz/a2H558jFbJc1fu019TXqbBotc+iMOmyqa3UKrqjFGY2AACQHAYbAABIDjJaXdhprJJK/DTZTo1VBlovnfg6K/F4OczKGF4asNsWyWhW8lDZk31d2bfC65LUVLv9MbzEo0KfrXQWk/VXheUqGc3LkbbsZTR7HuqclMTm5TiVvVlJZUo2i7nOXioqu6ptkYylpDJbVs+Y/bmd0OdThIymlnHtQDq79QAA0BEw2AAAQHIYbAAAIDl4NilQGWlVyKLXb73Wbet96K3VyVWdynrs67xnoLwP/91W+O3q8nBU6LUtF3k2tuz9DBX6WlaHj0nBo0K4vWdjy8qz8+drPRpfr7I3++PbY1b1aHz/qJVqY1ax9c+D8l5UnfLsyvp5MelqlPdZ5P2NQn9n9LUIAADGHAw2AACQHGS0qqhpq5/i2qmxn/6r6a4Py7TTaC9j2Gm8CotV8ouq88f3U3wvXbWqU98ryvhbNpu0P4YtF2UQsGX15rmSypQ0UpeM5vtCSWzq+quMzL6v1L2i2qZQsrJaPE/JWOrNf7+tl9iq1tUV+jx+fLl+HIUyWRGd12IAAOg4GGwAACA5DDYAAJAcPJvhQGWrVdqr1+ytFqxWcazqyyiPpiiViq1TGaEV7WSLVr6QPQ/vmflztJ6FCq9Vvoyvs9c/Jjuv6g/l2ah7Q4W+F21r2xMTwu5R/WH73IeX27JKF6PCmWN8mVNPrfY9lZG96D6yfw/8fdwBmZ2HfWbz0ksvheuvvz7MnDkz9PX1hXe9611h586dg/WNRiOsX78+zJs3L/T29oZly5aFXbt2pWgKAACMxcFm//794T3veU+YNGlS+Md//Mfw/e9/P/zRH/1ROO200wa3uffee8N9990X7r///rBjx44wZ86ccOWVV4ZDfv0MAAAYE9Quo91zzz1h/vz54eGHHx783TnnnNM0q9m0aVNYt25dWLFiRf67Rx55JMyePTts3rw53HjjjWFUEhNq6LdVGaHtNLroGEoqseGe/g1yW/aSgtqnClNWMpoKg1b484+RamxdjPznt7UyjpJ4fJ0tKxmtHZSMpq6NPcei868rTL1s9mbfj2UXAVQZmb3ENXVqc9nKY14qs/tRGaFVeHU7MlqHZQWIofbWP/nkk+GSSy4JH/rQh8IZZ5wRLrroovDQQw8N1u/evTvs3bs3LF++fPB3PT09YenSpWH79u11NwcAAMbiYPPDH/4wPPDAA2HhwoXhn//5n8NNN90UPv3pT4c///M/z+uzgSYjm8lYsvJAnae/vz8cPHiw6QMAAF0so504cSKf2WzYsCEvZzObzPzPBqCPf/zjg9uNc1PrTF7zvxtg48aN4c4776y7qQAA0KmDzdy5c8MFF1zQ9Lvzzz8/fPnLX85/zoIBMrJZTLbtAPv27TtptjPA2rVrw+rVqwfL2cwm84VGFVZPjVmpU2m2HhX6bMtKh4/JiKz0fa/nx2j4FvsfjKLM0WVX/FRpXvx/aJQvoTwkdY1VCpZ2UJ6JPY+YVVTr8mXKejRF2Zutv+F9EeunKK+lKPRZeTbW3/FeT9WVOn15ogiTt5CuRpNFoj3//PNNv3vhhRfC2Wefnf+8YMGCfMDZunXrYP2xY8fCtm3bwpIlS4bcZ+bpTJs2rekDAABdPLP5zGc+kw8amYz24Q9/ODz99NPhwQcfzD8ZmVS2atWqvD7zdbJP9nP2Ps61115bd3MAAGAsDjaXXnpp+MpXvpJLX3fddVc+k8lCna+77rrBbdasWROOHDkSVq5cmb+Xs3jx4rBly5Yw1U9bAQBgTJAkXc0HPvCB/NOKbHaTZRDIPl2Hes+mCKuvK8/Gv2dQNjV9kZ5fVcNX/kXMOygq/XpZr0H5OR7fNruteudB1bWzxIC6Nql9GI/ypYp8SbVUgPVCfF1Zr2X69NZ1GVaG9/ux5ZiUNNZfUksK+LK/V+y14z0bAACAODr7lVQAAOgIyPo8msKii6bNdortp/hlMzK3Exar9lNWcolZ4VRlqPYpeVS23KoyUjuS23BTlzxWVipTKZlUehZfViteqhBmL43ZspfGfPSqkuPKprLx7VapdNSqruNdP8asctphMLMBAAAGGwAA6HyY2QAAQHLwbEYalcrGew9W+/V1KpVN2aUBPMr78JQNBVap+f0+jh5tvR/vJ6hVJMv6Oe14Hyk8k+Ggqi/j08woz0alb1EpabwvYut8eLP1V9SSAt7D8X5OHb6M8mhiXndQz1SMZzpKYGYDAADJYbABAIDkIKMNN+2EPitUWLKqs/JPjMTm5TAlcamVGe33vDSj5AcvVaiwaCWxKekyRkYcaRlNhV4racyX7X5ipDIlI3nJSb2Jr8KbVZ3KCqCkspjwZtVuWy6S0caXlJyLXkVQdaNQZmNmAwAAyWGwAQCA5DDYAABAcvBsylKkgVbVSGM8m6qhkMprGI60O8ozsL7A4cOtv+e39WHRyjOwfo7Keu3Lvt+s9h6jn6cgJuu09WF8nypfzNepPlaZnJUvo3wRlUpGpaBRHo0Pm1bH8GHZyrOx/RET6nyi4r0yCj2ZIpjZAABAchhsAAAgOchoKagahphKRqtryu0lmLIyWh3STMaRI61Db62s5kOf7X5jZLSYBdqqZsSuiro32lnYTF0rez1UlgC1sFhdMlpVia1oP7atXn5T4c1KRot5bk+UfIXBS7wxry2MEMxsAAAgOQw2AACQHAYbAABIDp5NVapqrTFaqtJ6lffhqUu/VStuKq1feQYxno31CXzos/JsbJ3Xtv1qoMqzsf3o096odDV19b+6H1SaGZVZXF0P79FZz2I4PBvvmajVOG1ZrfDpt1Vtq5rZOWbFXU/ZVXX9/RZz/42Qn8PMBgAAksNgAwAAyUFGq4rKCKzCEtuZ0ioZRclodaEWQVNSTdVswV4OsrKGDYOOkdG8/OWlMvtdFU6qJLa6sjmUDTX3ZRV66/fpr4f9rup/JaN5aSomg4CSuKyMpr6npLmY8OYYGU2F/itOnKjnHovJUD5CMLMBAIDkMNgAAEByGGwAACA5eDZlKfJWbOihCkusqvWn8m/KZmseDs+mKFvwG2+0rrO+jA+Ltj6ND3WOCX0uGxadSj8vu+JmTAqaFJmdlUfjPRTlr8SETNs6la3Zl2N8GVsW/d8IzemCxgWR2kjdYyr0Wd1/eSMaoy51DTMbAABIDoMNAAAkh8EGAACSg2ejiHlfoqwOG6PR+rLShYfjHZyyvozS+mPes/F6uv2uepfG+zCqzu/HXle/bdl0NTGp4RXKF1Orcap0Nb6P1Xs36nqoOn+N1bssVT0bf28oH0aVlS8T8S6NfR6lRxPz96BqKqWhyqMAZjYAAJAcBhsAAEgOMlpZVGhhOyGzqs5LVaYN49Q0PkZSU9JMTHht2dUgY0KfDx9uLtt6lfVZyWhKNvPfTbGKZzvUsRpnUbZwFaZeVkYrkrGU5FVWDvMSmy0XybE1ZG/2MrYkRkY7dqxcXVHapeFeObYEzGwAACA5DDYAAJAcBhsAAEgOnk1VYlKDW59ApUspSiVSMo35SXry+Anl/Jy60tj70NuyKUC8DxOTkkZ5NrYck65GhT5XTVejVlT0KF+m6rVRHoWvjwkLjvFsym4bE7JcNpw5In1PjC/TFO6svJaie6zsveq9Ru/h1LXMRY0wswEAgOQw2AAAQHKQ0cpS9IauCmcsm3W4SEYTjCs5/T8p80DJUM+TylWzDvs31mMyCFQNb1b97+UHW1ahzzGrsVrqyuytMnT7PlaykZc8bTkmE0RZicvX+zp7DB/eXFXiU/ejeIVAcVKWgLKvN8RIZX41WvU9ldFE3XMxsm6bMLMBAIDkMNgAAEByGGwAACA5eDYKFSIYk61Vaa3Wa6gp9HlcRJoNmeYm5vhVV4q0/eH9A++nWH2/asioCkP15brS1aTwbOpKV6PKMb5I2bp2soDbsgrZVvetKyuPRmZvVp6tesZjQvi991g2lVI7mcYTwswGAACSw2ADAADJQUazePnDSxWWmGm0Lftps5I4Yt72t+0Ri6WpbNGFC7KVDcVV2atj3livGjIak0lXhal7+aHsQncqg0A7lM0goORQHxbtpcuy16pq5gH/3ap16lkpeGWgtHSmrnHVrACpMmGo+3gEw50tzGwAACA5DDYAAJAcBhsAAEgOno3CaptFunvZ0MeYsNQYz0b4NOoY40qGRXua0tx4PTsmLNr2m6rz5xiziqG9du14NsqXKVvXDmVDz1VdO56N/a7K7N1OeLVKn5RgFU0Z3lzX6w1+xdnDpvzGG63r/PeUn+Pv+RHyZRTMbAAAIDkMNgAAkBxktLL4aal6g9xTVuJQodbtoGQcc/yTpDFH6TBpdY6+LWUlNl9WEltMRm51HeuS0VptV2fW56oyWkyG6LISV4yMVlEq80RlAij75n1MeHPZcGYvnfnMzlY683Uqu4C6V2Puo4QwswEAgOQw2AAAQHIYbAAAIDl4NmWJ0dq9nmrxvkzJTM6F7VGegfIzVMi0a0/lMOmqfo7ycJQvE1NXdduY8OYUoc+qrmoqm5htq9bFtNVReRVNda1iUhuptFPKs/HhzUeOtA5vVn6O3a/ymooYRp+m6bB17/Ctt94Kn/vc58KCBQtCb29vOPfcc8Ndd90VTpjOaDQaYf369WHevHn5NsuWLQu7du2quykAADBKqH2wueeee8IXv/jFcP/994f//M//DPfee2/4gz/4g/Anf/Ing9tkv7vvvvvybXbs2BHmzJkTrrzyynDo0KG6mwMAAGNxsPm3f/u38Cu/8ivhqquuCuecc0741V/91bB8+fLwrW99a3BWs2nTprBu3bqwYsWKsGjRovDII4+Ew4cPh82bN9fdHAAAGIuezWWXXZbPbF544YXwMz/zM+E73/lOeOqpp/IBJmP37t1h7969+QA0QE9PT1i6dGnYvn17uPHGG8OooS6NPkXqCJXGXrXVezSqTmjtVVcDjfJz1Ds5yt9px2sZ7ndp1DkO9zs4MdvWdYwUvkyML1f1XRrvp1if1tepd2kOHWq9rfKFis5Rvb82Qp5N7YPNbbfdFg4cOBDOO++8MGHChHD8+PFw9913h2uuuSavzwaajNmzZzd9Lyvv2bNnyH329/fnnwEOHjxYd7MBACAhtQ9xjz/+eHj00UdzSeyZZ57JJbI//MM/zP+1jHNRWZm85n83wMaNG8P06dMHP/Pnz6+72QAA0Ekzm89+9rPh9ttvDx/96Efz8oUXXpjPWLIB44YbbsiDAQZmOHPnzh383r59+06a7Qywdu3asHr16qaZzbAPOGoqPlS9pSANTKU2eGmubF1MmpeKGarrktjC+OZUKifttw6Jq6o8OhyhzilktJhtazpGTBbmWqSyopREVbM3l01Jo7I8x6Sr8cewmZ2LpPmysrInocRW+54zo3+8a3Ampw2EPmch0dmAs3Xr1sH6Y8eOhW3btoUlS5YMuc/M05k2bVrTBwAAunhm88EPfjD3aM4666zwjne8I3z729/Ow5w/+clP5vWZVLZq1aqwYcOGsHDhwvyT/dzX1xeuvfbaupsDAABjcbDJ3qf5nd/5nbBy5cpcGste3MwizH73d393cJs1a9aEI0eO5Nvs378/LF68OGzZsiVMnTq17uYAAMAoYFwjc+Y7jMyzyQIFDtxzT5jW21vfjuvSgVUYog8vnjx56J8z+vpal32d2k/ZuqLQ57KrMbYTFmtoa8XFqn5KHV7MSHs2MdvWlH6+tmtVNfRcLSORYqmAdjybwxVX47Tt9ksKqBVX/XNtyzH+XgsOHjkSpv9/FLKyOEjECQAAyWGwAQCA5JD1ua7QZ1u2IYoe/z27bVE4dVlZT60+qc7DS2x++m23VZJbTRKbCnWOCqFW4dRe0qka+jkc0pmippDVqmHKUZJmCqmsbFYAX+/lMFunMgH476mQZSWjHRGZCPx5lF1909erPmalTgAAGEsgowEAQHIYbAAAIDl4NoqYTLJqhT/rmfiQRZWCQmV2Vv5OjJ9kfRr/PR9O2dPTetuyYdEqA21EdlqVEkdR2euJ9ILaDtGus20ViWlr6VQ+vr7qfax8GO91+JVz7TNQNbxZ+TIxWZ+PHm3dVv+3Qt3vKpv7KIGZDQAAJIfBBgAAkoOMVhY/pY0JvbRTWi9jWeko5hhq2xj5QdX5qbmS8WxZyWh+nzEyWtkQao/NSC3q2qGs5JVKGouSvFJkUCgbsuzrVZ0KYY6R0VR4c11ZAsrWFR3Dtk319wQn2/rlWWy9l8NHSGJjZgMAAMlhsAEAgOQw2AAAQHLwbCwqvNmHJSvPxPsyKl2M9QxUKhtfjjm+1YF9GKiq8xm17bYqtY1KZeP1dBUyXdWzqbr6ZFG9CsUOo5g6fJi60sz4euW9VF1hU/kgPhTZ3/N11MV4Nm++WU/2eP/3wP69Uon91TFqhpkNAAAkh8EGAACSw2ADAADJwbOpGnOuUsvE6LAW/z2bHsYf02u06h2YqnUqtY3yXtTKgD7m3+5HvZ/jyyn8nKHKZeuqbDdc9+1w+zIxy3Goe055j8qzGY73bGKWGFDteatgBeBW97h/r0alqylKtdWqrub7mJkNAAAkh8EGAACSg4ymiJmKqrDksplcfQoKFd7sQ49te1R4c4z84Kfmtt7XWXnMy39lw6J9ne8Pe4yq2aNjZLQ6JLXYbeuQzqrKZr4cs1KsqlOvAqiUNOpeVZmcU8loNrzZn1PZFDQxGeIV/h4v+70RhJkNAAAkh8EGAACSw2ADAADJwbOpitLFfVi0LXutt6y27Y8Z4+eo8Garfbfj2diy19OrejYxK36qsFDr/RSlxKnDsxnp0Oe6PBvlJ/h7vOwyFr5c1bOpK/TZr6KpVvFUXo+951XIdpG/q+4jFV5etMrvKICZDQAAJIfBBgAAkoOMNtLZo6tOd1WW15i3u204sd9nVRlNhUX7DAIq63NVGS1V6LMte6lOhWynQMkv6h5LFfpcVUZTWZ9V1vOqEltMCLWS0fw+1esNSio8fryee0fdj6MEZjYAAJAcBhsAAEgOgw0AACQHz6YqMaGvZbX+Iu1dhU2rbLEKpR/HhKzWEd6s0tP4betKV1MUCj3WPBtPjL9XNV2NSl+j7rmqYdHKB/Lluo5hz8M/p6r/J4h7JSazud/PcITfRzL6WgQAAGMOBhsAAEgOMlpZikJmreTjp9F2+uun1Hb6W5S51W4bM01WEkcMtn1e4lLhrPb8vTRhz8nvU8kIXv6yspaXFMourDZUuQ45tC5ZrWxm33ZktLILpKnQ+yIZTS00qLJdlJXYijJxKDms7EKH6t4out7jxT1m7yslIxdl21B/K0ZIYmNmAwAAyWGwAQCA5DDYAABAcvBsFEqj9xqp1VCV1uv1UpXJWaF8Cd/WqhqtCrVWKF9K+TL+eCq8U4V6Vk1BM1TZkiLcWd0PqUKh60hXo+p8W1T28qp+TsxKuVXvY39t7PX3z3+r7Yq2jfFsyr5e4L87SlLZMLMBAIDkMNgAAEBykNHKSip+SuulibKLIPn9KGlAEZOloGwYpP9eUSiyxbbdb6cWj7Nt8+dfVQ5T8ocn5s3rqnJkXaGmVSW2mNDnsnJcVYmtaFu1QJvaZ9nXEvw9qDJIKDk4JrN2DONFCLU9Dy+j+QUT7Xmp53gYw6CZ2QAAQHIYbAAAIDkMNgAAkBw8G6XRxmQL7u0tF5asssV6lGZdNQy3nYzUZb2fkfAolGdk8XW+/8u2fZSEk0bfN+34OVXr2vEwqmT2LvKh6vBeqj63MfeR91rUKwTqVYyqmc1rhpkNAAAkh8EGAACSw2ADAADJwbOpqhGrbb0Oq1LZpIjPr7pd0T7K+hQxnlHM8evQk+vyD8YKKfqjneukvJeyPoyn6vs6VT3D4Vjxd2KBDxOzrMYwwcwGAACSw2ADAADJQUYrm4G3KPS5bOhtqrBQRV3T5uFO15Jiuj9CEsKoJeZeLfu9om3LPiuqPTGZlOt6xkb6WR0f8XpB2W1JVwMAAGMJ/psHAADJYbABAIDk4NnUkcqmaNvRFIo72jyL0dYeSH9thsMnqoORDpMf34YvVnU/CRkdrQAAgDENgw0AACQHGW00SwrdDn0zNhlNUpliNLVlDLQ1uoXf+MY3wgc/+MEwb968MG7cuPA3f/M3TfWNRiOsX78+r+/t7Q3Lli0Lu3btatqmv78/3HLLLWHWrFlhypQp4Zd/+ZfD//zP/7R/NgAAMDYGmzfeeCP87M/+bLj//vuHrL/33nvDfffdl9fv2LEjzJkzJ1x55ZXh0KFDg9usWrUqfOUrXwmPPfZYeOqpp8Lrr78ePvCBD4Tjao0IAADoHhnt/e9/f/4ZimxWs2nTprBu3bqwYsWK/HePPPJImD17dti8eXO48cYbw4EDB8KXvvSl8Bd/8RfhiiuuyLd59NFHw/z588NXv/rV8L73va/dcwIAgLHs2ezevTvs3bs3LF++fPB3PT09YenSpWH79u35YLNz587w5ptvNm2TSW6LFi3KtxkTg00H6KcAIwbPR1dS62CTDTQZ2UzGkpX37NkzuM0pp5wSTj/99JO2Gfi+J/N4ss8ABw8erLPZAACQmCT/Bc8CB7y85n/nUdts3LgxTJ8+ffCTSW4AANClg00WDJDhZyj79u0bnO1k2xw7dizs37+/5TaetWvX5l7PwOfFF1+ss9kAANBJg82CBQvywWTr1q2Dv8sGlm3btoUlS5bk5YsvvjhMmjSpaZuXX345fO973xvcxpP5PtOmTWv6AADAGPZssjDlH/zgB01BAc8++2yYMWNGOOuss/Kw5g0bNoSFCxfmn+znvr6+cO211+bbZzLYr//6r4ff/M3fDDNnzsy/91u/9VvhwgsvHIxOAwCALh9svvWtb4Vf/MVfHCyvXr06//eGG24If/ZnfxbWrFkTjhw5ElauXJlLZYsXLw5btmwJU6dOHfzOH//xH4eJEyeGD3/4w/m2l19+ef7dCWodcQAA6FjGNTJnvsPIotGyGdKBe+4J03p7R7o5AABdy8EjR8L0227L/XRlcfBCCAAAJIfBBgAAksNgAwAAyWGwAQCA5DDYAABAcjpy8bSBALqDR4+OdFMAALqag///d7gosLkjQ5+zhdbIjwYAMHrI0oideeaZY2uwOXHiRPjxj3+cj6RZ1oLsJElhc/K7SNmATN+cDH3TGvqGvokl+zucLY6ZLRUzXiwf0ZEyWnZC2Qg6sNQA+dJaQ9/QN1XgvqFvYshesi+CAAEAAEgOgw0AACSnowebbOmBO+64I/8X6BvuG54p/t6MXjoyQAAAADqLjp7ZAABAZ8BgAwAAyWGwAQCA5DDYAABAcjp6sPnCF74QFixYECZPnhwuvvji8M1vfjN0Exs3bgyXXnppvuT2GWecEa6++urw/PPPN22TxX+sX78+f7u3t7c3LFu2LOzatSt0G1lfjRs3LqxatWrwd93cNy+99FK4/vrrw8yZM0NfX19417veFXbu3DlY381989Zbb4XPfe5z+d+W7NzPPffccNddd+WZSwbo5v6pTKNDeeyxxxqTJk1qPPTQQ43vf//7jVtvvbUxZcqUxp49exrdwvve977Gww8/3Pje977XePbZZxtXXXVV46yzzmq8/vrrg9v8/u//fmPq1KmNL3/5y43nnnuu8ZGPfKQxd+7cxsGDBxvdwtNPP90455xzGu985zvz+6Tb++a1115rnH322Y1PfOITjf/4j/9o7N69u/HVr3618YMf/KDR7X2T8Xu/93uNmTNnNv7+7/8+75u//uu/bpx66qmNTZs2DW7Tzf1TlY4dbH7u536ucdNNNzX97rzzzmvcfvvtjW5l3759WRh7Y9u2bXn5xIkTjTlz5uQPxgBHjx5tTJ8+vfHFL36x0Q0cOnSosXDhwsbWrVsbS5cuHRxsurlvbrvttsZll13Wsr6b+yYj+0/bJz/5yabfrVixonH99dfnP3d7/1SlI2W0Y8eO5VP+5cuXN/0+K2/fvj10KwcOHMj/nTFjRv7v7t27w969e5v6KXsBdunSpV3TTzfffHO46qqrwhVXXNH0+27umyeffDJccskl4UMf+lAuv1500UXhoYceGqzv5r7JuOyyy8K//Mu/hBdeeCEvf+c73wlPPfVU+KVf+qW83O39U5WOTMT5yiuvhOPHj4fZs2c3/T4rZzdBN5LNUlevXp0/KIsWLcp/N9AXQ/XTnj17wljnscceC88880zYsWPHSXXd3Dc//OEPwwMPPJDfL7/9278dnn766fDpT386/4P58Y9/vKv7JuO2227L/+N23nnnhQkTJuR/a+6+++5wzTXX5PXd3j9dNdgMkBm+/g+u/1238KlPfSp897vfzf8H5unGfsqWVrj11lvDli1b8gCSVnRj32RGdzaz2bBhQ17OZjaZuZ0NQNlg0819k/H444+HRx99NGzevDm84x3vCM8++2weWJIFA9xwww2h2/unKh0po82aNSv/H4efxezbt++k/210A7fccksujfzrv/5r0+JFc+bMyf/txn7KZNbsPLMoxYkTJ+afbdu2hc9//vP5zwPn3419M3fu3HDBBRc0/e78888PP/rRj0K33zcZn/3sZ8Ptt98ePvrRj4YLL7wwfOxjHwuf+cxn8ojGjG7vn64abE455ZT8j8jWrVubfp+VlyxZErqF7H9S2YzmiSeeCF/72tfyUE1LVs4eDNtPmd+V/dEd6/10+eWXh+eeey7/X+nAJ/vf/HXXXZf/nIWzdmvfvOc97zkpRD7zJ84+++zQ7fdNxuHDh09aBCz7z+1A6HO3909lGh0e+vylL30pD31etWpVHvr83//9341u4Td+4zfyCJivf/3rjZdffnnwc/jw4cFtsoiZbJsnnngiD9G85pprujZE00ajdXPfZKHgEydObNx9992N//qv/2r85V/+ZaOvr6/x6KOPNrq9bzJuuOGGxk//9E8Phj5nfTBr1qzGmjVrBrfp5v6pSscONhl/+qd/mr8vcMoppzTe/e53D4b8dgvZ/xWG+mTv3gyQhWnecccdeahmT09P473vfW/+cHQjfrDp5r75u7/7u8aiRYvy885eGXjwwQeb6ru5b7IBI7tPsnfWJk+e3Dj33HMb69ata/T39w9u0839UxWWGAAAgOR0pGcDAACdBYMNAAAkh8EGAACSw2ADAADJYbABAIDkMNgAAEByGGwAACA5DDYAAJAcBhsAAEgOgw0AACSHwQYAAJLDYAMAACE1/weznc2BmPOkTwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Let's now initiate the photonics TM FDFD class. Leave the objective empty for now, let's use the class to compute the source field first. \n",
    "# s0 and A0 do not have to be passed now, and in general don't need to be passed to do some EM calculations. \n",
    "\n",
    "# Setup geometry\n",
    "geometry = geometry.CartesianFDFDGeometry(\n",
    "    Nx=Nx, Ny=Ny, Npmlx=Npmlx, Npmly=Npmly, dx=dl, dy=dl\n",
    ")\n",
    "\n",
    "ldos_problem = photonics.Photonics_TM_FDFD(\n",
    "    omega=omega, geometry=geometry, chi=chi,\n",
    "    des_mask=design_mask, ji=ji, chi_background=chi_background, \n",
    "    sparseQCQP=True\n",
    ")\n",
    "\n",
    "# You can print the ldos problem to see the attributes.\n",
    "print(ldos_problem)\n",
    "\n",
    "ei = ldos_problem.get_ei(ji, update=True) # update = true sets the ei to the source field. Not required if you just need to do a Maxwell solve. \n",
    "plt.imshow(np.real(ei), cmap='bwr')\n",
    "\n",
    "vac_ldos = -np.sum(1/2 * np.real(ji.conj() * ei) * dl * dl)\n",
    "print(\"Vacuum LDOS: \", vac_ldos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now let's set s0. We need to restrict ei to the design region. \n",
    "ei_design = ei[ldos_problem.des_mask] # restrict the field to the design region\n",
    "c0 = vac_ldos\n",
    "s0_p = - (1/4) * 1j * omega * ei_design.conj() #* dl * dl  # the dl*dl factor is not needed here, as set_objective will take care of it.\n",
    "A0_p = sp.csc_array(np.zeros((ndof, ndof), dtype=complex))\n",
    "\n",
    "# We set the objective with set_objective(). \n",
    "ldos_problem.set_objective(s0=s0_p, A0=A0_p, c0=vac_ldos, denseToSparse=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Precomputed 2 A matrices and Fs vectors.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/alessio/code/dolphindes/dolphindes/photonics/_base_photonics.py:287: UserWarning: If both ji and ei are specified then ji is ignored.\n",
      "  warnings.warn(\"If both ji and ei are specified then ji is ignored.\")\n"
     ]
    }
   ],
   "source": [
    "# We are ready to set up the QCQP for calculating limits. We will use Pdiags = 'global': this represents two constraints (extinction and real power global conservation). We will show how to refine these constraints below, or you may pass Pdiags = 'local' to directly do the local problem (often slower).\n",
    "ldos_problem.setup_QCQP(Pdiags = 'global', verbose=1) # verbose has a few levels. 0 is silent, 1 is basic output, 2 is more verbose, 3 is very verbose.\n",
    "\n",
    "# get a copy of the ldos_problem QCQP for testing and comparing with GCD\n",
    "import copy\n",
    "gcd_QCQP = copy.deepcopy(ldos_problem.QCQP)\n",
    "gcd_QCQP_hs = copy.deepcopy(ldos_problem.QCQP)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found feasible point for dual problem: [5.41788396e-08 1.00000000e-01] with dualvalue 816.2080132586462\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Splitting projectors: 2 → Precomputed 4 A matrices and Fs vectors.\n",
      "previous dual: 1.990600873418748, new dual: 1.9906008734187208 (should match)\n",
      "4\n",
      "at step 1, number of constraints is 4, bound is 1.9551398986638036\n",
      "Splitting projectors: 4 → Precomputed 8 A matrices and Fs vectors.\n",
      "previous dual: 1.9551398986638036, new dual: 1.9551398986637052 (should match)\n",
      "8\n",
      "at step 2, number of constraints is 8, bound is 1.8715794692206325\n",
      "Splitting projectors: 8 → Precomputed 16 A matrices and Fs vectors.\n",
      "previous dual: 1.8715794692206325, new dual: 1.8715794692206873 (should match)\n",
      "16\n",
      "at step 3, number of constraints is 16, bound is 1.86061948342903\n",
      "Splitting projectors: 16 → Precomputed 32 A matrices and Fs vectors.\n",
      "previous dual: 1.8564669252943178, new dual: 1.856466925294244 (should match)\n",
      "32\n",
      "at step 4, number of constraints is 32, bound is 1.8346626259272552\n",
      "Splitting projectors: 32 → Precomputed 64 A matrices and Fs vectors.\n",
      "previous dual: 1.830639882980767, new dual: 1.8306398829804047 (should match)\n",
      "64\n",
      "at step 5, number of constraints is 64, bound is 1.7564576543713888\n",
      "Splitting projectors: 64 → Precomputed 128 A matrices and Fs vectors.\n",
      "previous dual: 1.753141079194479, new dual: 1.7531410791945412 (should match)\n",
      "128\n",
      "at step 6, number of constraints is 128, bound is 1.6056830226992571\n",
      "Splitting projectors: 128 → Precomputed 256 A matrices and Fs vectors.\n",
      "previous dual: 1.6047359180016323, new dual: 1.6047359180016803 (should match)\n",
      "256\n",
      "at step 7, number of constraints is 256, bound is 1.5765102131704438\n",
      "Splitting projectors: 256 → Precomputed 512 A matrices and Fs vectors.\n",
      "previous dual: 1.572773343069001, new dual: 1.5727733430692474 (should match)\n",
      "512\n",
      "at step 8, number of constraints is 512, bound is 1.5557944042003922\n",
      "Splitting projectors: 512 → Precomputed 800 A matrices and Fs vectors.\n",
      "previous dual: 1.5547870055683592, new dual: 1.5547870055685937 (should match)\n",
      "800\n",
      "at step 9, number of constraints is 800, bound is 1.5454809807765972\n",
      "Reached maximum projectors or pixel-level constraints.\n",
      "Newton iterative splitting took 43.612792015075684s to reach pixel level.\n"
     ]
    }
   ],
   "source": [
    "## iterative splitting with Newton\n",
    "ldos_problem.QCQP.solve_current_dual_problem(method='newton')\n",
    "results = []\n",
    "result_counter = 0\n",
    "t = time.time()\n",
    "for result in ldos_problem.QCQP.iterative_splitting_step(method='newton'): # When we reach pixel level constraints, the generator will return and stop this loop.\n",
    "    result_counter += 1\n",
    "    num_constr = ldos_problem.QCQP.get_number_constraints()\n",
    "    print(f'at step {result_counter}, number of constraints is {num_constr}, bound is {result[0]}')\n",
    "\n",
    "    results.append((num_constr, result[0]))\n",
    "print(f'Newton iterative splitting took {time.time()-t}s to reach pixel level.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found feasible point for dual problem: [2.055518e-07 1.000000e-01] with dualvalue 789.381050098289\n",
      "Precomputed 2 A matrices and Fs vectors.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "At GCD iteration #1, best dual bound found is             1.9906008734198435.\n",
      "At GCD iteration #2, best dual bound found is             1.7469384268205286.\n",
      "At GCD iteration #3, best dual bound found is             1.6904015371494634.\n",
      "At GCD iteration #4, best dual bound found is             1.669117289771638.\n",
      "At GCD iteration #5, best dual bound found is             1.6611703027651563.\n",
      "At GCD iteration #6, best dual bound found is             1.6452481749189358.\n",
      "At GCD iteration #7, best dual bound found is             1.6361471371787073.\n",
      "At GCD iteration #8, best dual bound found is             1.6302547947455026.\n",
      "At GCD iteration #9, best dual bound found is             1.6182099497837665.\n",
      "At GCD iteration #10, best dual bound found is             1.6086168530098797.\n",
      "At GCD iteration #11, best dual bound found is             1.603718146915681.\n",
      "At GCD iteration #12, best dual bound found is             1.589971296147006.\n",
      "At GCD iteration #13, best dual bound found is             1.5853419862753348.\n",
      "At GCD iteration #14, best dual bound found is             1.584176382488307.\n",
      "At GCD iteration #15, best dual bound found is             1.5808257653080457.\n",
      "At GCD iteration #16, best dual bound found is             1.5805080948875025.\n",
      "At GCD iteration #17, best dual bound found is             1.5782609239532965.\n",
      "At GCD iteration #18, best dual bound found is             1.5771810785746125.\n",
      "At GCD iteration #19, best dual bound found is             1.576804060346324.\n",
      "At GCD iteration #20, best dual bound found is             1.5767441307951382.\n",
      "gcd took time 3.398860216140747 to reach 1.0202287510538184 of pixel dual.\n"
     ]
    }
   ],
   "source": [
    "from dolphindes.cvxopt import gcd\n",
    "\n",
    "### now compare with tightening the bounds using GCD\n",
    "\n",
    "## gcd parameters, play around and see how the result changes\n",
    "\n",
    "# maximum number of QCQP constraints before merging, larger values may lead to tighter final bounds but makes GCD slower\n",
    "max_cstrt_num = 10\n",
    "\n",
    "# maximum number of GCD iterations\n",
    "max_gcd_iter_num = 50\n",
    "\n",
    "# check to see how much the bound improved after gcd_iter_period number of GCD iterations\n",
    "gcd_iter_period = 5\n",
    "\n",
    "# relative tolerance for required minimum improvement of bounds or GCD terminates\n",
    "gcd_tol = 1e-2\n",
    "\n",
    "t = time.time()\n",
    "\n",
    "gcd_params = gcd.GCDHyperparameters(\n",
    "    max_proj_cstrt_num=max_cstrt_num,\n",
    "    orthonormalize=True,\n",
    "    opt_params=None,\n",
    "    max_gcd_iter_num=max_gcd_iter_num,\n",
    "    gcd_iter_period=gcd_iter_period,\n",
    "    gcd_tol=gcd_tol,\n",
    "    ortho_metric=\"euclidean\"\n",
    ")\n",
    "gcd_QCQP.run_gcd(gcd_params=gcd_params)\n",
    "print(f'gcd took time {time.time()-t} to reach {gcd_QCQP.current_dual/ldos_problem.QCQP.current_dual} of pixel dual.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found feasible point for dual problem: [6.21838174e-07 1.00000000e-01] with dualvalue 726.018736123337\n",
      "Precomputed 2 A matrices and Fs vectors.\n",
      "At GCD iteration #1, best dual bound found is             1.9906008743566286.\n",
      "At GCD iteration #2, best dual bound found is             1.7469424006141143.\n",
      "At GCD iteration #3, best dual bound found is             1.6904022844291133.\n",
      "At GCD iteration #4, best dual bound found is             1.669118646168242.\n",
      "At GCD iteration #5, best dual bound found is             1.6612231598651965.\n",
      "At GCD iteration #6, best dual bound found is             1.6450098782106113.\n",
      "At GCD iteration #7, best dual bound found is             1.633875885312799.\n",
      "At GCD iteration #8, best dual bound found is             1.6269193755699332.\n",
      "At GCD iteration #9, best dual bound found is             1.6169850384810138.\n",
      "At GCD iteration #10, best dual bound found is             1.610662001684383.\n",
      "At GCD iteration #11, best dual bound found is             1.603090291203674.\n",
      "At GCD iteration #12, best dual bound found is             1.5913546683547983.\n",
      "At GCD iteration #13, best dual bound found is             1.5862866331491277.\n",
      "At GCD iteration #14, best dual bound found is             1.5792711288131023.\n",
      "At GCD iteration #15, best dual bound found is             1.5773331231907368.\n",
      "At GCD iteration #16, best dual bound found is             1.5773213500490737.\n",
      "At GCD iteration #17, best dual bound found is             1.5693330528214062.\n",
      "At GCD iteration #18, best dual bound found is             1.568301511814236.\n",
      "At GCD iteration #19, best dual bound found is             1.564442968095412.\n",
      "At GCD iteration #20, best dual bound found is             1.5655280951436892.\n",
      "gcd took time 3.877166271209717 to reach 1.0129714403583396 of pixel dual.\n"
     ]
    }
   ],
   "source": [
    "# We may also use the Hilbert-Schmidt metric, which is slower per iteration but may lead to faster or better convergence\n",
    "gcd_QCQP = copy.deepcopy(ldos_problem.QCQP)\n",
    "t = time.time()\n",
    "gcd_params = gcd.GCDHyperparameters(\n",
    "    max_proj_cstrt_num=max_cstrt_num,\n",
    "    orthonormalize=True,\n",
    "    opt_params=None,\n",
    "    max_gcd_iter_num=max_gcd_iter_num,\n",
    "    gcd_iter_period=gcd_iter_period,\n",
    "    gcd_tol=gcd_tol,\n",
    "    ortho_metric=\"hilbert_schmidt\"\n",
    ")\n",
    "gcd_QCQP_hs.run_gcd(gcd_params=gcd_params)\n",
    "print(f'gcd took time {time.time()-t} to reach {gcd_QCQP_hs.current_dual/ldos_problem.QCQP.current_dual} of pixel dual.')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "dolphindes",
   "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.10.19"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}