{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-12-19T12:08:39.496262Z",
     "start_time": "2018-12-19T11:10:49.113830Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==============================\n",
      "Epoch #:1\n",
      "==============================\n",
      "----------------- I-1---------------------\n",
      "19 0.2329746928277702\n",
      "39 0.23296960625078644\n",
      "59 0.2329683031558343\n",
      "79 0.23296830231118462\n",
      "99 0.23296830229031726\n",
      "119 0.23296830228982082\n",
      "139 0.23296830228979928\n",
      "159 0.23296830228979892\n",
      "179 0.23296830228979884\n",
      "199 0.23296830228979884\n",
      "----------------- I-2---------------------\n",
      "19 0.797467250763613\n",
      "39 0.797467250763613\n",
      "59 0.797467250763613\n",
      "79 0.797467250763613\n",
      "99 0.797467250763613\n",
      "119 0.797467250763613\n",
      "139 0.797467250763613\n",
      "159 0.7974672507636129\n",
      "179 0.7974672507636129\n",
      "199 0.7974672507636129\n",
      "----------------- I-3---------------------\n",
      "19 0.23296063519219054\n",
      "39 0.23295898638401366\n",
      "59 0.2329588988009353\n",
      "79 0.23295889805188805\n",
      "99 0.23295889789682397\n",
      "119 0.2329588978941721\n",
      "139 0.232958897894096\n",
      "159 0.23295889789409135\n",
      "179 0.23295889789409077\n",
      "199 0.23295889789409024\n",
      "----------------- I-4---------------------\n",
      "19 0.9524239664387182\n",
      "39 0.9524239664387182\n",
      "59 0.9524239664387182\n",
      "79 0.9524239664387182\n",
      "99 0.9524239664387182\n",
      "119 0.9524239664387182\n",
      "139 0.9524239664387182\n",
      "159 0.9524239664387182\n",
      "179 0.9524239664387182\n",
      "199 0.9524239664387182\n",
      "----------------- I-5---------------------\n",
      "19 0.7389694176862874\n",
      "39 0.7389693923368936\n",
      "59 0.738969375147863\n",
      "79 0.7389693463220711\n",
      "99 0.7389693463220711\n",
      "119 0.7389693463220711\n",
      "139 0.7389693433340967\n",
      "159 0.7389693433340967\n",
      "179 0.7389693423334336\n",
      "199 0.7389693423334336\n",
      "Round robin best migration\n",
      "==============================\n",
      "Epoch #:2\n",
      "==============================\n",
      "----------------- I-1---------------------\n",
      "19 0.23295889789409024\n",
      "39 0.23295889789409024\n",
      "59 0.23295889789409024\n",
      "79 0.23295889789409024\n",
      "99 0.23295889789409024\n",
      "119 0.23295889789409024\n",
      "139 0.2329588978940902\n",
      "159 0.2329588978940902\n",
      "179 0.2329588978940902\n",
      "199 0.2329588978940902\n",
      "----------------- I-2---------------------\n",
      "19 0.23296817479060336\n",
      "39 0.23296245325841197\n",
      "59 0.23295953733840824\n",
      "79 0.232958899616963\n",
      "99 0.23295889793705846\n",
      "119 0.2329588978952983\n",
      "139 0.23295889789421256\n",
      "159 0.23295889789409466\n",
      "179 0.23295889789409124\n",
      "199 0.23295889789409044\n",
      "----------------- I-3---------------------\n",
      "19 0.23295889789409033\n",
      "39 0.23295889789409033\n",
      "59 0.23295889789409027\n",
      "79 0.23295889789409027\n",
      "99 0.23295889789409027\n",
      "119 0.23295889789409027\n",
      "139 0.23295889789409027\n",
      "159 0.23295889789409027\n",
      "179 0.23295889789409022\n",
      "199 0.23295889789409022\n",
      "----------------- I-4---------------------\n",
      "19 0.23296830228979878\n",
      "39 0.23296830228979878\n",
      "59 0.23296830228979878\n",
      "79 0.23296830228979878\n",
      "99 0.23296830228979878\n",
      "119 0.23296830228979878\n",
      "139 0.23296830228979878\n",
      "159 0.23296830228979878\n",
      "179 0.23296830228979878\n",
      "199 0.23296830228979876\n",
      "----------------- I-5---------------------\n",
      "19 0.23296830228979884\n",
      "39 0.23296830228979884\n",
      "59 0.23296830228979884\n",
      "79 0.23288175962275934\n",
      "99 0.2327043026152132\n",
      "119 0.23266538612106336\n",
      "139 0.23266382214800546\n",
      "159 0.23265739742314323\n",
      "179 0.23265739742314323\n",
      "199 0.23265694156412065\n",
      "[0.2329589  0.2329589  0.2329589  0.2329683  0.23265694]\n"
     ]
    }
   ],
   "source": [
    "%reload_ext autoreload\n",
    "%autoreload 2\n",
    "import numpy as np\n",
    "import de2 as de\n",
    "import multiprocessing as mp\n",
    "import makeqx as mkq\n",
    "import matplotlib.pyplot as plt\n",
    "from jupyterthemes import jtplot\n",
    "# currently installed theme will be used to\n",
    "# set plot style if no arguments provided\n",
    "jtplot.style()\n",
    "\n",
    "def rndtop5(x):\n",
    "    return np.round(x*2.0)/2\n",
    "\n",
    "def tmm_wrapper2(arg):\n",
    "    args, kwargs = arg\n",
    "    return mkq.vgdr_eval_wsweep(*args, **kwargs)\n",
    "\n",
    "def arc_par(pop, **kwargs):\n",
    "    jobs = []\n",
    "    pool=mp.Pool(90)\n",
    "    #num_layers = int(pop[0].size/2)\n",
    "    for indiv in pop:\n",
    "        #jobs.append((indiv.reshape(num_layers,2), 0))\n",
    "        jobs.append((indiv, 0))\n",
    "    arg = [(j, kwargs) for j in jobs]\n",
    "    answ = np.array(pool.map(tmm_wrapper2, arg))\n",
    "    pool.close()\n",
    "    return answ\n",
    "\n",
    "num_layers = 8\n",
    "norlam_min = 0.25   # this is lam/d\n",
    "norlam_max = 1.5 # this is lam/d\n",
    "lam_pts = 256\n",
    "\n",
    "\n",
    "bnds = [(0, 1)]*num_layers\n",
    "its_first = 0\n",
    "psnew = 160\n",
    "islands = 5\n",
    "its_second = 200\n",
    "\n",
    "numgens = 2\n",
    "\n",
    "# isl = np.ones((psnew*islands, num_layers))\n",
    "# for indiv in range(int(psnew*islands)):\n",
    "#     isl[indiv] = np.random.uniform(0,1, num_layers)\n",
    "# # #     nxhold = isl[indiv][::2]\n",
    "# # #     nxhold[::2] = 1\n",
    "# # #     nxhold[1::2] = 0\n",
    "# # #  isl[indiv] = qn.make_qxn(num_layers, dlevels=2)\n",
    "# poplist = np.split(isl, islands)\n",
    "\n",
    "bids = np.zeros(islands, dtype=int)\n",
    "bfits = np.zeros(islands)\n",
    "bests = isl[0:islands]\n",
    "\n",
    "for gen in range(numgens):\n",
    "    print(\"==============================\")\n",
    "    print(\"Epoch #:\" + str(gen + 1))\n",
    "    print(\"==============================\")\n",
    "    \n",
    "    for isln in range(islands):\n",
    "        print(\"----------------- I-\" + str(isln + 1)+ \"---------------------\")\n",
    "        poplist[isln], bids[isln], bfits[isln], bests[isln] = de.de_cp(\n",
    "            fobj=arc_par, \n",
    "            bounds=bnds, \n",
    "            pop=poplist[isln], \n",
    "            history=[], \n",
    "            it_start=its_first, \n",
    "            popsize=psnew, \n",
    "            its=its_second, \n",
    "            lam_low=norlam_min, \n",
    "            lam_high=norlam_max, \n",
    "            lam_pts=lam_pts)\n",
    "        \n",
    "#     if its_second > 64:\n",
    "#         its_second = int(its_second/2)\n",
    "    \n",
    "    if gen < (numgens - 1):\n",
    "        print(\"Round robin best migration\")\n",
    "        stmp = np.copy(poplist[islands-1][bids[islands-1]])\n",
    "        for isln in range(1, islands):\n",
    "            poplist[isln][bids[isln]] = poplist[isln-1][bids[isln-1]]\n",
    "        poplist[0][bids[0]] = stmp \n",
    "    \n",
    "print(bfits)\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-12-19T16:54:28.866885Z",
     "start_time": "2018-12-19T16:54:28.402050Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fd20c5d90f0>]"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAFWCAYAAACVXsE9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XmYVPWdLvD3d5bai6YXlgYaodlkEUXQBiQJMeqgQS6TYExAx5mbZPQ6TpyYu8wkhtzHG3N1cjUm4T5e75Mn5k4iEqIRRyLEDaMoiDQuiCI7zdKs3V1de50659w/qrtt6K26u6rO0u/neeaZ2FV1zsuh6W//djFp0iQTREREeZCsDkBERM7BokFERHlj0SAioryxaBARUd5YNIiIKG8sGkRElDcWDSIiypuSz5u+8IUvYMGCBRg7diwOHz6MRx99tMf3er1erFy5ErNnz0Y2m8Vbb72FDRs2FCwwERFZJ6+iEYlEsHnzZkyYMAG1tbW9vvdrX/saQqEQfvCDH8Dv9+Pee+9FJBLBli1bChKYiIisk1f31Pvvv4/3338fra2tvb5PVVVcffXVeP7555FIJHD+/Hm8/PLLuOaaawoSloiIrJVXSyNfo0aNgizLOH78eMfXGhoaUF1dDSEETLP7HUsCgQA0TStkFCIiGiBVVZFIJLp9raBFw+fzIZPJwDCMjq8lk0nIsgyPx4N0Ot3lM4FAAKtWrSpkDCIiGqSnnnqq28JR0KKRSqXg8XggSVJH4fD7/dB1HZlMptvPtLcwnnrqqQG3NsLhMMTl30G64U2kG94YWPgiC4fDiEajVsfokd3zAfbPaPd8ADMWQqnzearnwT91OSJ/uT/vzwwmo6qqWLVqVY8/jwtaNE6fPg1d1zFu3Dg0NDQAAGpqatDY2Nhj11Q7TdMGXDQ0TQO0LDRTtW0312D+fKVg93yA/TPaPR/AjIVQ6nwCCuRstl/3LGbGvAbCJUmCoiiQZRlCiI7/fTFN07Bjxw4sW7YMfr8flZWVuP766/HWW28VPHgXegZCVot/HyKiEhKyF9C776mxQl4tjZtuuglLly7t+O81a9Zg3759ePTRR3HPPffgwIED2Lx5MwBg/fr1WLlyJX7yk59A13Vs3bq1JNNtTT2Te7hERC4iZA9MpxWNjRs3YuPGjd2+tmbNmgv+O51O48knnxx8sn4y9QzAlgYRuY3shal3nURkFfdsI8KWBhG5kN1aGq4pGrnuKY/VMYiICkrIHluNabiqaIBFg4hcRrB7qkjY0iAiN2L3VHGYeppFg4hcJzflli2NgjMNjd1TROQ6QmFLoziybGkQkQvJXhaNYjANjUWDiFxHyB6YWXZPFZypp9k9RUSuk5tyy6JReHoGQmLRICJ3EeyeKg4u7iMiV+KU2+Iw9QyEwm1EiMhd2D1VLO2VWOKmhUTkFoJ7TxVL+0NlFxURuUbbzzNuI1IEHZWYRYOIXKJj5+4sWxqFx5YGEblM+zgtu6eK4LPuKQ6GE5E7CHZPFZGpwzSybGkQkXu0/zwzNGtzdOKeooG2asxpt0TkEkL22moLEcBlRYNHvhKRm+Sm27JoFI3JnW6JyE1YNIordxATWxpE5A65A5jsM3MKcF3RyHBMg4hcg91TxZZlS4OIXET2wrTRwj7AZUWD54QTkZvY7ahXwJVFgy0NInKH3JgGu6eKhmMaROQmdtvhFnBZ0QCn3BKRm9js1D7AZUWD3VNE5CZ2O4AJcF3RyAAsGkTkEkLxcsptMZnZNI98JSLXELKPe08VFafcEpGbsKVRXBzTICI3ye1ym7I6xgVcVjQ45ZaI3EMoXoDdU0XEKbdE5CJC9rF7qpjYPUVEbiIUdk8V1Wd7TwmroxARDZLgeRrF1vFwOa5BRE4neyCExCm3RdW2hTDHNYjI6Tq62tnSKJ72lgbHNYjI6doXKnNMo4jaN/Zi0SAip2v/OcYxjWIyNJimwTENInI8ofhgGlnAyFod5QLuKhoAj3wlInew4RYigAuLBo98JSI3ELLPdqvBAdcWDbY0iMjZ7LgtOuDKopEBFLY0iMjZcpsVsmgUH8c0iMgF7LiFCODComHqGRYNInI+G25WCLiyaLClQUTOZ8dt0QGXFg2u0yAipxMyB8JLg2dqEJELCMVnyzENJZ83SZKEFStWoK6uDkII7Nq1C+vWrUM223WlYllZGW699VZMnToVAHDw4EE8/fTTaGlpKWzyHph6GkLxl+ReRETFImQvjHTE6hhd5NXSWLJkCaZMmYIHHngAq1evxpgxY7B8+fJu3/uNb3wDsizjBz/4Af7lX/4F6XQat912W0FD94YD4UTkCk6ePbVo0SJs2rQJkUgEsVgMGzduxIIFCyBE18OOqqqqUF9fj3Q6DU3TsGPHDowdO7bgwXtiZjmmQUTOJ2Sv7bZFB/IoGn6/HxUVFWhoaOj4WkNDA4LBIMrLy7u8/9VXX8XcuXPh9/vh9XpRV1eH3bt3FzZ1b/QMxzSIyPFyYxr2Kxp9jmn4fD4AQDKZ7PhaIpG44LXODh48iIULF+KRRx4BABw/fhyPPfZYn0HC4TA0TcsvdQ+fBwBJBYQn0PHfdmLHTJ3ZPR9g/4x2zwcwYyGUIp+s+uFVBTwDvNdAM6qq2uvrfRaNVCrXp+b3+xGPxwEAgUDggtfaCSFw7733or6+Hr/4xS8AADfffDPuuece/Ou//muv94lGowMuGuFwGNFoFADgibfCL5SO/7aLzhntyO75APtntHs+gBkLoVT5hksqkrEWaAO412Ay9lU0+uyeSiaTaGpqQk1NTcfXampqEI/H0dzcfMF7A4EAKisr8dprr0HTNGiahi1btqC2thbBYHBAf4D+MrNpgN1TROR0Tt6wcOvWrbjxxhtRVlaGUCiEpUuXYtu2bTBN84L3xeNxnDlzBosXL4aiKFAUBV/84hfR1NTU0UopOq4IJyIXsOvW6Hmt09i8eTNCoRBWr14NSZJQX1+PDRs2AABWrlwJAFi7di0A4PHHH8eKFSvw0EMPQQiBY8eO4fHHHy9S/K445ZaIHE9SICTZllNu8yoahmFg/fr1WL9+fZfX2otFu8bGRvzyl78sTLoByC3u8wIQAMy+3k5EZDtCzk0ycmz3lJN0PGS598EcIiK7Em1rzew45dZ1RaO9D5BdVETkWG0/v9jSKAFTzwBg0SAi5xKKF6ZpAG0/z+zEhUWjrTLzyFcicqjcFiIZ2HFc1nVFA2xpEJHD2XVbdMCNRQO5wSMWDSJyKrsewAS4tWhwgR8ROZlNNysEXFw0uJUIETmVsOkWIoBbi0Y23THPmYjIaYTsteUWIoBLiwa4lQgROViupcGB8JLhmAYROZmQOaZRUqae4ZGvRORcNj0fHHBp0UA2zSNficixOOW2xNg9RUROJhR7nqUBuLposKVBRM7ElkaJmdk0xzSIyLEExzRKjFNuicjJ2NIoLVNPsWgQkWNxw8ISy60I91kdg4hoQHJbo7OlUTJmNsVtRIjIsXJjGiwaJWNmUwBbGkTkREKCkD0c0yglU0+xe4qIHKl9PJYtjRIysykImUWDiByorWudGxaWUjYFIcmApFqdhIioXzp+4WVLo3Tap6qxi4qInEZ0tDRYNEqmvS+QRYOInEa0H/VqGlZH6ZY7i0ZbXyAX+BGR0wjFb9uFfYBLiwaMLEwjC6H4rU5CRNQvuZZG0uoYPXJn0UD7Wg22NIjIWdjSsEhuVTjHNIjIWVg0LMKtRIjIidg9ZRU9zTENInIcFg2L5FaFs6VBRA6j+mFqLBolxzENInIiO5+lAbBoEBHZCgfCLcLt0YnIiTimYRHudEtETsSWhkVMnUe+EpHzCMUPsKVhAY5pEJEDsXvKIlzcR0ROxNlTFuHsKSJyHEmBkFWu07ACzwknIqdp/5nFloYFOHuKiJymfesjFg0LcGt0InKaz4oGu6dKzsymIYQEcP8pInIIofhg6hnA1K2O0iMXF41cpeYMKiJyCrsv7ANcXDSQTQMAxzWIyDHsvkYDcHHRMPW2osEZVETkEHZfowG4uGjA1GHqGRYNInIMYfOzNAA3Fw1wVTgROQzHNKxlZtMAxzSIyCGc0D2l5PMmSZKwYsUK1NXVQQiBXbt2Yd26dchms92+/7LLLsPNN9+MkSNHIp1O45VXXsHLL79c0OD5MPUUhMqiQUTOkJs9Ze/uqbyKxpIlSzBlyhQ88MAD0HUdd999N5YvX45nnnmmy3tnzJiBVatW4cknn8T+/fvh8XhQUVFR8OD54KpwInISofhgpJqtjtGrvLqnFi1ahE2bNiESiSAWi2Hjxo1YsGABhBBd3rts2TK8+OKL+PTTT2EYBlKpFE6ePFnw4HnhpoVE5CBOWKfRZ0vD7/ejoqICDQ0NHV9raGhAMBhEeXk5mpqaOr7u8Xgwfvx47Ny5Ez/60Y8QDAZx6NAh/P73v0dzc+mrJ3e6JSInEYrP1gcwAXkUDZ8v90M3mfzsD5JIJC54rV0gEIAkSZg/fz7WrFmD1tZW3HLLLbjzzjvx0EMP9XqfcDgMTdP6/Qfo/PmLyUKH7A9D6eY1K3SX0U7sng+wf0a75wOYsRCKlU/xBiEUQC3A9QeaUVXVXl/vs2ikUrmmkt/vRzweB5ArDp1fu/i9W7Zswfnz5wEAzz//PB555BGUlZUhEon0eJ9oNDrgohEOhxGNRrt8PZCKQUgq4t28Vmo9ZbQLu+cD7J/R7vkAZiyEYuYrkzxIRZuRGeT1B5Oxr6LR55hGMplEU1MTampqOr5WU1ODeDzepcsplUrh/PnzME1zQGELLdc95bc6BhFRXpwwppHXQPjWrVtx4403oqysDKFQCEuXLsW2bdu6LQ5vvPEGrr32WpSXl0NRFCxbtgxHjhzptZVRLNwenYicxAl7T+U15Xbz5s0IhUJYvXo1JElCfX09NmzYAABYuXIlAGDt2rUAgJdeegmBQADf//73IYTAwYMH8cQTTxQpfu84EE5EjiFkCNlj+5ZGXkXDMAysX78e69ev7/Jae7FoZ5omnnvuOTz33HOFSTgY2RQEz9MgIgf47KhXe7c0XL6NCMc0iMgZWDRswNS5YSEROYMTzgcH3F40smmOaRCRIwjFD9PIAkb3e/rZhcuLRrKtaHTd7oSIyE6E6rP9WRqA64tG7vQ+TrslIttzwBoNwO1FQ8/9BXAGFRHZnRPWaABuLxptVZszqIjI7pxwlgbg8qKBtu4pzqAiIrtzwql9gNuLBkzOoCIiR2BLwyZyM6jYPUVE9saWhk18Nu2WiMi+hOIHOOXWemxpEJETsKVhE6aWglBZNIjI3oTKdRq2wJYGETmBUAMwtbjVMfrEokFEZANCCcDg7CnrmVqS3VNEZHtC9cPUElbH6JP7iwZbGkTkAFynYRMsGkRke7IXQpI5pmEH7J4iIrsTagAAuDW6HXBxHxHZndRRNDimYTl2TxGR3Qkl0HZqn2Z1lD65v2iwe4qIbM4pM6eAoVA0OloaPPKViOwpt7CPRcMW2qewcVyDiOxKqAGYWRYNWzA1nt5HRPYmFL8jZk4BQ6FotLc0OK5BRDbF7ik7MTSYRpYtDSKyLRYNmzE1TrslIvsSSsARW4gAQ6VoZJMAu6eIyKY45dZmuMCPiOxMqAEYHAi3D3ZPEZGd5Xa4ZUvDNswsV4UTkX1xINxmzGwSElsaRGRTLBo2k9t/KmB1DCKibggIxcfZU3Ziagl2TxGRLQnFDyEkRxzABAyVopFNsKVBRLYkPEEAgJlh95Rt5GZPsWgQkf10nNrH2VP2YWpxdk8RkS0JNZjbrNA0rI6SlyFSNJIQatDqGEREXeQW9jljPAMYKkUjm+DiPiKyJUkNOmYQHBgqRUNLQMgqIKlWRyEiukCue8oZ4xnAECkaRttfCGdQEZHdCE+ALQ274UFMRGRXQg3CzLBo2Es2BdM0OO2WiGyHA+E2lVsVzqJBRPbCMQ2b4v5TRGRHnD1lU6YW5063RGQ7wsOWhi3lztTgAj8ishehcPaULXGnWyKyI+EJOmazQmDIFQ2OaRCRjcheCEnh7Ck7ynVPsWgQkX107HDroKKh5PMmSZKwYsUK1NXVQQiBXbt2Yd26dchmsz1+RlVV/PCHP0QoFMJ9991XsMADZWoJSP5Kq2MQEXWQ2sZZXTcQvmTJEkyZMgUPPPAAVq9ejTFjxmD58uW9fubmm29GS0tLQUIWAruniMhuhBqAqWcAQ7M6St7yKhqLFi3Cpk2bEIlEEIvFsHHjRixYsABCiG7fX1NTg1mzZmHz5s0FDTsYLBpEZDdOm24L5NE95ff7UVFRgYaGho6vNTQ0IBgMory8HE1NTRe8X5Ik3HbbbVi3bl2/goTDYWjawKttOBzu9XWhmJC9oT7fV0xW3jsfds8H2D+j3fMBzFgIhconwlUQ2URR/rwDvaaq9r4beJ9Fw+fzAQCSyWTH1xKJxAWvdXbdddfhxIkT2LdvH6ZOnZp30Gg0OuCiEQ6HEY1Ge32P4j2HsBzo833Fkk9GK9k9H2D/jHbPBzBjIRQyn2+kDDXVWvA/72Ay9lU0+uyeSqVSAHItjnaBQOCC19pVVVVh8eLFePbZZ/sdtNg6jnwVQ2bCGBHZnPCEYGoxq2P0S58tjWQyiaamJtTU1ODcuXMAcmMW8Xgczc3NF7x3ypQpCIVCWL16de7iigKfz4eHH34YTzzxBA4dOlSEP0J+2vsNc6svnfWXRETuJNSQo7ZFB/Kccrt161bceOONOHToEHRdx9KlS7Ft2zaYpnnB+3bu3Ik9e/Z0/HdtbS3uuOMOPPjgg4jHrX0w7X8xuQNPWDSIyHrCE4KROGd1jH7Jq2hs3ry5owUhSRLq6+uxYcMGAMDKlSsBAGvXroWmaReMSyQSCZimidbW1iJE7x8zm8ydqcH9p4jIJiRPEHrkiNUx+iWvomEYBtavX4/169d3eW3t2rU9fm7fvn22WNiXY7ZNu2XRICJ7cGL31JAaFTa1RMcKTCIiqwlPCEbGvjPFujPEikYcwsOiQUT2IHmcdQATMBSLBlsaRGQHQoZQ/DAzzpqYM7SKRoZjGkRkD8ITAgAYLBr2ZWhxd+8/JXshl10CyV9ldRIi6kN7V7nTuqfymj3lFqYWhxwcaXWMgpMCI+CfvgKeMXUQkgwA0GONSO5/AZmGNy1OR0TdkdRQbodbPWN1lH4ZckXDbd1TStV0hK66F3rLQUTf+gmyLYcgPCF4xy5EcNbt8Iyag1j9447aeploKBCekOPGM4Ah1j3ltu3RlfLJCM//z0gf+jOi236KbNM+wMjCTLUgdfBFRF7/AeRh4xCu+y4g9b4JGRGVVm66LYuGrZkZ97Q0hK8cobrvInXoJSQ/fa7b9xiJs2jd+iCk4EgEZ99R4oRE1BtJDTluPAMYYkXDcFH3VPCKb0GPHEXy466r9Dsz0xHEdvwcnrHz4Rn/+RKlI6K+CE+Q3VN217E9Oro/cdApPDWLoFRMRvy9XwEw+3y/3noMiY9+h8CsVRC+8uIHJKI+CU+Y3VN2Z2pxCCG1FQ6Hkr0IzPgaknv/CCPV1Pf726SPvg695RC7qYhsQvKEYDpsCxFgqBWN9u3R1ZDFSQbON/kmmNkM0odf6fdn4x/8Buqoy6FUzShCMiLqD+EdBiPNomFrppbIbY/ucWbREIofvklLkNz7DGDq/f68ET+N9OFXEJi1Ek7voiNyOrY0HMGEqcUhObRoeCd+CWY6gsyJdwZ8jeSnGyAFRsAz5qoCJiOi/hKeYSwaTmBmYs5saUgKfLV/heT+PyGfwe+emFoc6UN/hm/aX4OtDSKrCAhP0HHbogMsGo7hGXM1ICRkjr816GulDv4Zsr8SavWVBUhGRP0lPEEIIcHkmIb9GZkYJAcOhHsnXof00b8ARnbQ1zK1ONINf4Gv9q8KkIyI+kt4hgEAWxpO4MSWhjxsPJTySUgfea1g10wdehlK5TTIw8YX7JpElB/J075ZYdrqKP3GouEA3vGfh3bmIxjJcwW7ppE4A+3U+/DV3lCwaxJRfoR3mCNbGcAQLBqGFofkpCNfhQzPuAXIHCv8FuepQ3+GZ9wCCE+44Ncmop5JnrAjxzOAIVg0nNbSUEfNBiQFmcb6gl87e+5j6PEz8F6yuODXJqKeCW/YkdNtgSFZNKIQqnN+s/aMXQjt5LtFOw8jfegleCd+CZx+S1Q6kifsyNXgwBAsGkYm5pzFfZIKz+jLkTm5o2i3SJ/YBkkNQZRPK9o9iOhCwqGrwYEhWDTMTCy3YaGQrY7SJ3XkLJimAe3snuLdJJtCpvFdSNXzi3cPIrqA8HAg3DHa968XDhgM91RfBe3UewPaZ6o/0g1vQoy4AkJx8O6/RA4iednScAyj7aQs2w+GCxnq6DnInNxZ9Ftlz30CaDF4xtYV/V5E1NbSSLdaHWNAhlzRgJ6GqWdsvypcqZoOIanQzu4uwd1MGI3b4an5XAnuRUSStwwmi4Zz5Kbd2nsGlad6HrTTHwB6piT3M069A7VyKqTg6JLcj2ioEkoAQlZhpCNWRxmQIVk0DNuv1RDwVM9DpvHd0t0yeRba+U/hrVlYunsSDUHCl9t3ii0NBzEzUUhe+7Y05PJaCE8QmdPvl/S+mRPb4RnLWVRExSR5hsHUMzCzSaujDMiQLBpGOmrr7inPyNnInt8HZFMlvW/m5LuQgqO4iSFREQlvmWMHwYEhWjTMTCsk7zCrY/RIHTm7RAPgFzLTEWTP7WVrg6iIJF8ZTIeOZwBDtGjYuaUh1BDk8lpopz+05P65LipOvSUqFsnr3Om2wBAtGnYe01BHzISZboXe2mDJ/TON70LyV0IeXmvJ/YncTjh4ui0AKFYHsIKRbu04Octu1FGzoZ0pfddUOzMTg3Z2Dzxj65BsOWRZDhrahK8cnlGXQw6NASQFevwUtFPvwUictTraoEneMujRE1bHGLCh29KwafeUOnI2tDPWdE21+6yLijvfUmkJbxnk6bdj+PWPwjf5yxC+cgg1AG/N5zD8+kcRuuofIfkrrY45KMI7zLFrNICh2tLItEIoXkD22uq4RXnYeAjvMGhnPrI0h3ZqF4JXfBPy8InQ2dqgElFHzUHwyr8HEo2Ibv0Jss37L3hdDo9DYNZKDFv8IGI71yB71tp/JwMleYc5untqaLY02vaxt1trQx15GfTmQzC1mKU5TC2O7LlP4KmeZ2kOGjq8E69H6OrvILX/Bei7HutSMABAjx5HdNu/InVwE8Lzvwd19JUWJB08yVvGloaVhAB8IYFsv+pfAl6lBYGKEPTWpqJl6yyfjKFxE5GN1MMfLn0tvzifHNmGwITFwHH7/F7R/7/n0rJ7PsCeGb2XfAH+KV9G/MOfQZz5CL5QqPeMjS/AVCOoWvA3iO1sRbbEreFBPUPZA58/DU2NQC/iv3NfSCAWA0yz8Nd2fNHwhSR8+Z4AgEA/P/kgcCUAjCx8qB71lfGZtv9fykyddc63L/d/i6zK0pP+/j2Xmt3zAfbL+Enu/+YBn33v95VxL4CHLvpMKQ3mGf4YmAMUO/fvH4ohGTUKfl3HF41UzMCf1iQQi/WvS2fYNd9H6vArRT0Vr7NQKNRrRnnYeITr7kPLa/+1ZJsUdtZdvnDdd5E5uwfpQy+VPE93+nqGVrN7PsBeGeVQNcJ19yG5/wWkG97o+Hr+GQVC8+4GAMR2/u8ipexqMM9QGV6L0Lx70PLKfQVOdaFQKIRUrPAFA3BB0TBNIBUz+11RlagJLVuBVBEqcbf3Q+8ZfSOnI36uCcmW0m4d0q67fOaRvfCMXYhkdLMlmS7W1zO0mt3zAfbJKNQQhtXdh8iRj5HY8/oFr/UnY2rb71B27f+EUbkY6SOvFSFpV4N5hp5wBaSYUfS/AwVmUbqmgCE6EA7kthIRNtpKRKmcBu3cXqtjXCDTuBPK8ImOn+JI9hOc8y2Y6QgSH/6/QV3HTDUj8dFT8M+41Vb/nnsifMNhpFqsjjEoQ7ZoGOlWG60KF1AqpyLb9KnVQS5gxE8j23oMavVcq6OQi3jGXQN1xEzE6v8PYGQHfb1Mw1YY0RPwT7+lAOmKS/KVw0w1Wx1jUIZ00bDLqnA5PBaSJ5zb2dZmtMad8FRfZXUMcgnhK0dg9u1I7HkaRuJMga5qIr77t/CO/7ztd2iW2NJwLjMdgeQrszoGgFzXVLb1GMy288vtJNNYD6Vyqm03eCRnCV7xTejNhwo+/qC3HEbm5A74L/1qQa9baCwaDmakIpC8Nioa5+3VNdVOjxyFkTjv2IVUZB+esfOhVk5F/L1fFeX6yb1/hDrqcsjlk4py/UKQfOUw2D3lTEa6pW3gzPr9ldTKacjabBC8s1wXFVeH0yAoPgRmrURy73MwUsVZUGvEGpE58Q78k79clOsXgvCypeFYZroVQlIgPEFLc0iBEZD8FdBsNgjeWaZxJ9QRMwHFZ3UUcqjApV+FkYkjVeQ1P6mDL0KtngspaLdFqQBkDyRPECaLhjO17/0ieYdbmkOpvBR67LStv5GyTfthZpNQR862Ogo5kDysBt6J1+em15p6Ue+lR44ie+4T+CbdWNT7DET7zxqnd0/ltbhPkiSsWLECdXV1EEJg165dWLduHbLZC6fLKYqCr3/965g2bRrC4TBaW1vx+uuv47XXSrPopl/0DAwtkeuiiloXQ62chux5+3ZN5ZjQTr0Hz+i50Eq0gp7cIzDrNmROvFOy7/PUgRcRuuo7SH7yrOWbf3Ym+YbD1DO2nPDSH3kVjSVLlmDKlCl44IEHoOs67r77bixfvhzPPPPMBe+TJAmRSAS/+MUvcO7cOYwdOxbf+c53EIlEUF9fX5Q/wGDkZlBZ3dKYhuT+f7c0Qz4yjTsRnHsXIOSi/7boSEKCGD4VvuqxkAMjAUmBmYki23oM2bN7YCTPW53QEuroK6FUTEbklf9SsntqZz6EnjgL78QvIbXv+ZLdty+5mVPO3d22XV7dU4sWLcKmTZsQiUQQi8WwceNGLFiwAEJcOIicyWTwwgsv4OzZszBNE8ePH8fu3bsxaZI9ZzNYPYNK+IZDDo1G9px9xzPaaWf3QAgZStV0q6PYi+yBb8pSDL/h55Avvxtq1UyYpg4zE4PwlsE/5WYMv+ExhBf+M5QRs6xboQi3AAAc4UlEQVROW1pCRmDm15E6sKlog989SR3cBN/E63K/5NiE5K8o+XMohj5bGn6/HxUVFWho+OzM6oaGBgSDQZSXl6OpqeeHIEkSJk+ejJde6nvwKxwOQ9O0PGN3//n+kvU4pPAIqAP47EBcnFGMvAJmugVBOQmUKENv+nqGZtPHCI5fACN1tESJuhrI33OxiPJLIU9fBZgG9CN/gnl6J6Cn0fnHlAHACFZDGvt5hOffB7NpL/R9vwcs/OFRqmcojVsMSQ1AOfV6v+856IyRjyDkVRg2YQHMcx8M7lrdGEg+adhoCK21ZM9/oPdRVbXX1/ssGj5fbsZMMpns+FoikbjgtZ7ceuutSCaT2L59e59Bo9HogItGOBxGNNr/gYlA/DyEGkB8AJ/tr+4yBiaMh352b0nu35d8nqHn2DsIzPw6ovW/AlCk3dB6MdC/52LwT/sKfFNvRmr/C0juewEwtJ7zRaPAqX2QPnkegdl/C2XePyNe/7glx/qW6hkKJYCyS25E4pP1SLec69dnC5UxcOwtSKPmI3Z466Cv1dlA84WUYdCjp5G06OdNvvoqGn12T6VSuV1X/X5/x9cCgcAFr3Xnq1/9KiZNmoRf/vKX0HV79oEbqRZLu6cURwyCf0Y7/T6EJwy5vNbqKNYREoJX3glv7fWIvv0wknv/CBj5/bJjJM4itv1/IXXgTwjVfRf+aX9d5LDW8U1dBiPdgvTRv1iWIX30dagjL4Pkq7AsQ2eSv8IVY1t9Fo1kMommpibU1NR0fK2mpgbxeBzNzd1PHbvlllswY8YMPPbYY4jH7TtTwEi3QlhUNIQahDxsnG1XgnfH1BLIntsLz+ghuoGhkBC88j9BrZqB1jcfGGDBN5HavxHRbT+Ft/YGBC7/O9hhgWkhSYER8NXegORHT8OKFmk7vbUBeuQIPJd8wbIMnUn+ShhJ549p5DUQvnXrVtx4440oKytDKBTC0qVLsW3bNpjdbNj+ta99DZdeeil+9rOf2eawl56Y6RbLZk8pFVNhagno0ROW3H+gMqfq4akemluKBGauhFo1Da1vPQgj1jioa2XPfYzoWz+Bp3oegnO+DQj3LJnyT/8atHOfQDu72+ooSB95Hd7xn4flhVlScrOnXNDSyGvK7ebNmxEKhbB69WpIkoT6+nps2LABALBy5UoAwNq1a1FRUYFrr70Wmqbhxz/+ccfnDxw4gDVr1hQh/uAYyWZI3mGApBRki+b++Gy/Ket+ExuITGM9grPvgBSqHvQPTifxTrgW3ksW5wpGvDC7s+qtx9C69UEMu+afEZzzbcR3/V847fvhYnL5JHjGXIXW1++3OgoAIH1iGwKzVkEdOQvaGeuKmOQrB4ChUzQMw8D69euxfv36Lq+tXbu24383NTXhrrvuKly6ImtfmSl5h8NI9m+wbrDUymnInHy3pPcsBDPVjGzzQXiq5yK1f6PVcUpCqZqBwGW3I1b/OPSWwwW9thE7ida3foJhi1bDP+NWJD9eV9Drl1pg5jeQPvYm9Ohxq6PkZFPInNwBT80ia4uGvzK3sC9j796XfLinTTwAZiYKU9cg+ctLe2PZC3n4BGgOGgTvLNNYP2TGNYS3DKF5/4DUvn8v2mp4I3YK0XcegW/idfDW3lCUe5SCWj0PStkEJPc+a3WUC6SPv537fpW9lmVwyyA4MMSLBtA2g8pX2qKhVEwGDB16xLr1DoOROVUPpWIyRImfmxWCc74NPXoCyU83FPU+evNBxHauQWDmN6COnlPUexWFkBGYcStSB/5ku33Usuc+hplNwmPhCZRuGQQHWDRgpJogSjwlT62YhmzTfsdux2FET0KPNcLj8jM2vBO+BKViMuK7nkApxhq00+8jsef3CF35nyCFxhT9foXknfglCMWH5IEXrY7SlWkgc2I7POMWWhYhVzTY0nAFI9Vc8u4ppXKaY7um2mUa3T2LSgqNRmDWSiQ++E1J/7GnD21G5lQ9wnX/BKEESnbfwRBKAP5py5H85BlAT1sdp1vp429DHTGr7Qyd0pMCVTASLBquYCabS9s9JWQoFZMdtT6jO1pjPZSqGY75wdZfwcv/I7RT7yNzou/dDAot/v6vYWpJBOfdDcuniubBN3UZjFQL0g1vWB2lR3rLYRiJM/CMnW/J/eXACOgFOxPdWkO+aORaGqXrnlKG1wJCQrb5YMnuWQzZ5oMwMzGooy63OkrBeWo+B7nsEsQ/+q01AQwNsR0/h1I2Af7pK6zJkCfJXwVf7fVI7rF2IV8+0sfehteSLiqRa2kUaKq21Vg0UqVtaShV05BtPpT31hP2ZSLTWA/VwsHFYhBqCIGZ30Dykz9YOqBrpJoQe/cX8E2+CaqNj9r1z/w6suc/tXQ6a74yJ96GUj4JUqC0p/oJ33AI2QM9cbak9y0WFo1kU0lXhSsV0xzfNdVOO1UPz6jZgNT7BmdOEpj5dRiJs0gfftXqKMg27UPio7UIzfl7SKFqq+N0oVTNgKd6LuIf/c7qKHkx4meQbTkEz9i6kt5XDoyAqWu2m1U2UCwaqWYI2QPhCZXgbgJK5VRHbVLYG+3cJ4AJqCNmWB2lIOThE+GpWYT4B0/CLl0t6cMvI3PqPYSvvtfSdQZdCBmBy25H6tBLMKInrU6Tt8yJd0peNKTgCBiJs7DL99RgsWi0Vf9S7IQpl42HUHy56bZuYGSROf0+VJcs9AvMWoXMsa3QI0esjnKB+Ae/hmnoCM75ltVROngnXg/JE0Ly0+esjtIvmRPvQCm7pKQtNzkw0jVdUwCLBmBouQV+gcqi30qpnAY90gAzm+z7zQ6hNda3rdew/yyf3qhjroZSdgkSn/zB6ihd6WnE3v051JGXwTdpidVpIHzlCFz6FST2rAOyPR+PYEdG8jy0pv0lnUUlBUe2tTTcgUUDgJE4Bykwouj3UR12fkY+Mmc+gFADUCqmWB1l4CQVgRlfR3L/Rphpe57hbMTPIL7rCfhn3Aql8lJLswQv/ztkWw4ic/wtS3MMVObEdnhL2EUlBUbAiLNouIqePAfJX4KWRsU0aC4ZBO+QTUE7uweeMVdZnWTAfLXXA0IgZcfVzJ1op95Dav+fEJp3j2VbuHjGLoBaNQPx939tyf0LIXNyB6RQNeTwuJLcz01rNAAWDQCAkTgPudgtjcBISL4yZM/vK+59LJA5uR2eMVfDiV1UQgnAN+Xmfp3AZ6Xk3mehtx5FaN49gJD7/kABCU8YgctuR+KT9Y7ubjFTLcie/xSecSXoolJ8uc0KY6eKf68SYdFA7hjOYo9piOFToEdPwsy0FvU+VtAad0F4wlAqp1odpd98U26CkWpB5lhhz5EuHhOxnY9D8pcjcPnflvTOwcv/Dnr8FNKHXi7pfYshc+IdeMYUv2jIwdEwTQN6/HTR71UqLBoAjOQ5SP6qot5DKpvsvq6pNmY2Ce3Mh5Zt0TBQwlsGX+0SJPc+AydNhzS1GGLbH4VnzNXwTf5ySe7pnXg9lBEzEa//P3DSs+pJ5uS7kIIjIJdNKOp95FA1jMQ5R7Ri88WigVz3lOQdVtR58GL4ZNcNgneWObEdnuqrHHVsqX/qMujR49Aa662O0m969Dhi7/4S/ukroFYXdzxJHj4RgZnfQPy9X8FwSd+8mWlF9uzHRf9FRw5VQ3fZCZfO+RdeRHrbqX3F6qKS/JUQ/krXrATvTubUexCKz/KZPfmSAiPgnXAtEh93PY3SKbJnP0Liw39DaO5dUEbMKso9hBpEaN4/IN3wF2iNzjtpsjfpE9uLvtDPjccis2gAQDYFIxOF7C/OYLhSOQ1m6rxr9tPvlp6GdvoDx3RR+af9NbLn9yJ77mOrowxK+ugWJPf+EeGr/6nwBVv2IDT/ezDTrUh8tLbv9zuM1rgTkm84lPLJRbuHHBrNloZbGYnzkALFGddQKqfBbDlQlGvbSfrEttzUWymvo+ctI4fHwlNzDRIf23Ah3wCkDvwJqQN/Qnj+fVAqpxXmokJCaN49kNQAotsfdVWffDtTi0M781ERWxuC3VNuZiTOQAoWZ/dLtXIajCFQNLTTHwAQUEddYXWUXvmnr4DWWA+95ZDVUQom+elzSB3YjPCC/zr41p6QEJzzbShllyC67acwtVhhQtpQ5mR7F1Xhp4tL/goIxcui4VZ67DTk4KiCX1d4yyCHx8Jsdsl+U70xtNxq25pFVifpkVw+CeqoOUjsfcbqKAWX/PSPiH/4/xC88k74pi7HgH4QKj6Erv4nKFXT0fr2w+7uUkXbdHE1WLgWWidyeCwMLe6a3W3bsWi00eOnIBWhaKiVl8JINQNJd8w66Uvm2Faooy6H8IStjtKtwPRbkDn+lqN2Zu2PTMMbiG77KXyTbkB44X/r1/Y48vCJKPvCA5D8FYi++T9gxNz5jDr7bLp44buo5GE10CPHCn5dq7FotDFipyAHR6LQzVSlajq0s58U9Jp2lm0+ACNxFp5xC6yO0oUyYiaUiqm51d8ulj33MSJbvg9Tz6Ds2ofgn/mNXrcdkfxVCMz+Wwz73I+gnf0YrW8+4PoWRme5hX5XF3y6uFw2HnprQ0GvaQf2HrEsIT1+GkL25Jb8F/AfjFo1HamDm+Ap2BXtL31sK7w1n0P60EtWR7lAYPrXkD7y2pD4gWimWhB751GoIy+Hf9py+CYtQfb8PkixI/BGTgEwIfkroVRMbTvjZR+ib/0E2Sb3bXPTl8yp9xC84ltQqmYge/ajgl1XHjYe6YObC3Y9u2DRaGOmIzCzSUih0QX7oZIbzxgD7dwn8AyhNl3m2FvwX/pVyOFx0KPHrY4DILf1uRyuRnT7I1ZHKSntzAfQznwAedh4qKOvhK9qCrzDc9u9GKkItLN7kPjwN9CjJyxOaiE9jczp9+AZW1e4oiGpkEPVyLa6r3uKRaOT3GD4aGTP7inI9dSq6TCSTTDip4GwPfv4i8FInkf23MfwXvIFJD56yuo4uVPmpq9A8sAmV+79lQ+9tQF6awPUxjCi0ajVcWwnc+IdBK/4JhIf/AYw9UFfTw6PAYSwzS9NhTSEfv/tW24wvHDTbpWq6bkjUYeg1JEt8NR8DpCt75jzjv88hBpE6uAmq6OQTWmnP4CQZKgjC7OyXh42PvfLop4pyPXshEWjEyN2GnIBj4FUq6YjO0SLhtZYDxgZeMZaPCAue+G/9CtIfrrBcafMUQkZGjKN9QXb0UAZPhHZliMFuZbdsGh0osdO5pqVBSB8wyGHqqGdH5pFA6aO9JEt8E28ztIYvtobYOoZpI+8ZmkOsr/MiXfgGT0XkNRBX0upmIRsszsX9LJodKK3Hs8t8CvAbrdq1Yy28YyhsT6jO6mjWyAPq4FcPsmS+ws1BN+UpUh+8kxB+qnJ3bQzu2GaBtSRswd3IdkDedh4ZJtYNFxPj52EaegFaW2oIy+DdmZ3AVI5l5lqgdZYb1lrwz91GYz4GWRObLfk/uQwpg6tceegu6iUsgmAaUCPHC1MLpth0ejMyMKIn4I8bLBnBwuoI2dDO/NhQWI5Werwy/CMqYPwlpX0vlJoNLy11yOx52m44dAgKo30ie3wjJ4DoQQGfA2lYnJuPMOlrVsWjYtkW49DCdcM6hpy2SUQnhC0Ai4Ucqrs+b3QWxvgm3xTSe8bmLkK2qn3HL/1OZVW9uweGOkIPDXXDPgaSvlk145nACwaXeitxwfd0lBHXY5s0wGYWqJAqZwt+enz8E34EoQnVJL7qSMvhzpiRlsrg6g/TKSPvAbvhGsH+HmRW1nu4lM6WTQuokcLUDRGXsauqU600+9Bj5+Cr3ZJ8W8mZAQuW4XUwU0wEmeLfz9ynXTDm5BDo6FUTO33Z+XyWgjF6+r1WSwaF9EjDZB85QPugxdqAErFlLazJahdct+/w1t7/aD6ivPhm7IUkD1I7nuhqPch9zIzrcicfBfeCV/q92fVkbNzs6ZcvCaIReMiRuIMjHQUSnntgD6vjJgFMx117cyJgdJOvgsz1QzvpL8q2j2k8Bj4p/4HJD54EtDTRbsPuV/6yGvwjLmq31v8qyNmub6XgUWjG9mWQ1CGD2xtgWfkbGhnd4Mzdi5mIrH3j/BPvqlIM6kEgld8C5nGd9nKo0HLnt8LI3EG3ksW5/0Z4QlDKZ/k+qn2LBrd0JsPQhnQgjQBddRsaKfd/ZvGQGkndyAbaYB/+oqCX9tbez3k4Ggkdv+u4NemoSl5YBN8k5bkvX+aZ2wdjMRZ6JEjxQ1mMRaNbmSbD0Iun4j+HsiklE+G8IShneFvuj1JfPQ7eGs+B6V8SsGuKZddgsCMW5HY/W8wM9zBlQojc2wrTD2Td2vDO24hMsffLm4oG2DR6Ea25RAkNQgpNLpfn1PHzIN2dg+n2vZCbzmM9JFXEbjiPwLS4HfmF4ofoXn/iHTDm1z5TYVl6kjt3wj/lJsBxdfrW6XgSCgVU5Bm0RiazEwMevQk1MpL+/U5z5irkDn5bpFSuUfikz9AKF74p98y6GsF53wbZjZpj3M7yHXSR1+HqcVzhaMXvkk3QTv7cW47dJdj0eiBdvYjqCPy31tfLp8EyVcB7VR9EVO5RDaF+M7H4au9AeqoKwZ8Gd/U5VBGzEDs3V8ChlbAgERtTB2Jj56Gb9ISyOHu129Jvgp4L/kCkvs2lDicNVg0eqCd+QjKiJnId1zDW7MI2pkPYWZixQ3mEtnm/Uh+/AeE5t3d4z/G3ngnfAn+qcsQ2/FzGImhu5MwFZ925gNkTmxHcN7d3W6b7p/5dWSbDw6Zs3NYNHqgnf8EQvFCzme9hqTCM3YB0g1vFD+Yi6QOvoj08e0IL/xvkPqxs7Bv0k0IzFqFWP3/HjL/UMla8d2/hRASwnXfvaBweCdeD8/oOYi//2sL05UWi0ZPsilkmw7ktbe+Z/QcACa00+8XP5fLJD54EtqZDzFs0f1QR17W+5tlH4Jz/h7+S5cjtuPnudMBiUohm0LrW/8TwlcOpe6H8F/6VYSu+kcEZq1EbNf/hRE7aXXCkmHR6EWm8V14x/V9XKm39gZkGt4EjGwJUrmNifh7v0LqwCaE6u5DcM7fd521JnvhnXAtlPmrIQ+fiMgb/53TmqnkzHQErW/8dxjH34A8vBamrqH19fuhNQ6tyS+Dn/PoYpnj2xGY+Q3I5ZOgNx/s9j1y+SQo5ZMRr3+8xOncxERq/wvQzuxGYMatGP6ln0KPnYKRPA+h+CGXjYeZboXR8CpaP3nBtecUkAPoaRjHXkEs+pzVSSzDotELM9MK7fSH8NYsQqKHouGf/GVkTu6AkTxf4nTuo0eOILrtYUj+KihV0yD5ymFmU9B3H0W25RDCoQALBpHFWDT6kD76OkLz/gHJTzfATEcueE2pvBTq6DmIvP4Di9K5k5E8h8yxc1bHIKJu5FU0JEnCihUrUFdXByEEdu3ahXXr1iGb7dqH35/3OoF2+j1kI0fhv/SrSHzQaYaEpCJw2e1IH34FRnToDIIR0dCW10D4kiVLMGXKFDzwwANYvXo1xowZg+XLlw/6vU6R2LMW3vGfgzrm6ravCITm3gUhq0juHbp9m0Q09OTV0li0aBGeeeYZRCK57pmNGzfim9/8Jp599lmYpjng93amql0XzeRLVdVBfb5PsQZoe/4Nw+d9G9kz8yGHRgNKANFtD0MRGpDHvYuecZDsng+wf0a75wOYsRDsng8YXMa+Ptdn0fD7/aioqEBDQ0PH1xoaGhAMBlFeXo6mpqYBvffigKtWrer7T2O5HcAoADic+8+apVaGISIqGlVVoWldt+fps2j4fLndHZPJZMfXEonEBa8N5L2dX3/qqae6DUdERKWnqmrHz+6L9Vk0UqncWbd+vx/xeBwAEAgELnhtIO/trKdwRERUer39Et/nQHgymURTUxNqamo6vlZTU4N4PI7m5uYBv5eIiJwnr9lTW7duxY033oiysjKEQiEsXboU27Zt63Zguz/vJSIiZxGTJk3q86d557UXkiShvr6+Y+3FypUrAQBr167t871ERORseRUNIiIiwCHbiKiqih/+8IcIhUK47777APS98ryUK9O7y3fHHXfgqquuuuB+P/vZz3D06NGS5RtsBjtktMNzBIDLLrsMN998M0aOHIl0Oo1XXnkFL7/8si2eYV8Z7fAMH3vssQv+W1VVNDY24sc//nFeGYqdsa98dniGAFBWVoZbb70VU6dOBQAcPHgQTz/9NFpaWkr2DB3R0vjKV76CCRMmYNy4cR0/lG+66SbMmTMHa9asga7ruPvuu3Ho0CE888wzeb1e7Hx33HEHEokE/vCHP3T7mVLkG2wGO2S0w3OcMWMG/uZv/gZPPvkk9u/fD4/Hg4qKCpw8edIWz7CvjHZ4hhe7//77sXPnTmzevDmvDKXOeHE+uzzDu+66C0II/PrXv4ZhGLj99tsRCASwZs2akj1D25+nUVNTg1mzZnX85bVbtGgRNm3ahEgkglgsho0bN2LBggUQQuT1erHz9aVU+QaTwQ4Z+1KKjMuWLcOLL76ITz/9FIZhIJVK4eTJk3ndv1TPsLeMfSn13/OECRNQXV2Nbdu25Z2hlBm7y9eXUuWrqqpCfX090uk0NE3Djh07MHbs2LwyFCqjrbunJEnCbbfdhnXr1l3w9b5WnieTyX6vTC9kvnbz58/H/PnzEYlE8Pbbb+PVV1+FaZoDWjk/UAPNUKpn2FvGwf4ZCpHR4/Fg/Pjx2LlzJ370ox8hGAzi0KFD+P3vf49UKmWLZ9hbxvap7nb4Xmy3cOFC7Nmzp2OrITt9L3aXr50dnuGrr76KuXPnYvfu3TAMA3V1ddi9e3dJn6Gti8Z1112HEydOYN++fR19eEDfK8/bf+D0Z2V6IfMBwGuvvYZnn30W8XgcEyZMwLe//W2YpolXX311QCvnB2IwGUr1DHvLONg/QyEEAgFIkoT58+djzZo1aG1txS233II777wTTzzxRK/3L9Uz7C3jQw89ZPkz7ExVVVx11VX4zW9+0/E1u3wv9pQPsP77sN3BgwexcOFCPPLIIwCA48eP47HHHivpM7Rt91RVVRUWL16MZ599tstrnVeet+u88ryv14udDwCOHTuGWCwG0zRx+PBh/PnPf8bcuXPzyl8og8lgh4yD/TMUQvt1tmzZgvPnz0PTNDz//POYMGECDMPo9f6leoa9ZSwrK7P8GXY2b948ZDIZ7N69u0t+q59jT/kA678PAUAIgXvvvReHDx/Gvffei3vvvRd79+7FPffcU9JnaNuiMWXKFIRCIaxevRoPP/ww7rzzTvh8Pjz88MOorq7udeV5KVam95avtra2y/vbf8AA1q2c708GO2Ts6/VSZEylUjh//nyPi1Pt8Az7yngxK78Xr7nmGmzbts2234vd5euOFc8wEAigsrISr732GjRNg6Zp2LJlC2prayFJUsmeoW2Lxs6dO3H//ffjwQcfxIMPPojf/va3SKfTePDBB3H06NE+V54Xe2V6X/nmzp3b0ewbP348lixZgvfee6/j86VYOT/YDHbIaIfn+MYbb+Daa69FeXk5FEXBsmXLcOTIEUQiEVs8w74y2uEZAsCoUaNQW1uLt99+u8trdniOveWzwzOMx+M4c+YMFi9eDEVRoCgKvvjFL6KpqQnxeLxkz9ARU24BYOrUqbjrrru6XafR3crzUq9Mvzjf9773PYwdOxaSJKGlpQVvv/02Xn755Y6/oFLkG2wGO2S0w3MUQmD58uVYuHAhhBA4ePAg1q1bh+bmZls8w74y2uEZAp9NTX/00Ue7vGaH59hbPrs8w+rqaqxYsQKXXHIJhBA4duwYnn32WRw7dqxkz9AxRYOIiKxn2+4pIiKyHxYNIiLKG4sGERHljUWDiIjyxqJBRER5Y9EgIqK8sWgQEVHeWDSIiChvLBpERJS3/w/CTYztUydq4AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 460.8x403.2 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%reload_ext autoreload\n",
    "%autoreload 2\n",
    "import numpy as np\n",
    "import makeqx as mkq\n",
    "import matplotlib.pyplot as plt\n",
    "from jupyterthemes import jtplot\n",
    "# currently installed theme will be used to\n",
    "# set plot style if no arguments provided\n",
    "jtplot.style()\n",
    "\n",
    "qx = bests[np.argmin(bfits)]\n",
    "\n",
    "\n",
    "norlam_min = 0.25   # this is lam/d\n",
    "norlam_max = 1.5       # this is lam/d\n",
    "lam_pts = 256\n",
    "#lams = np.linspace(lam_low, lam_high, endpoint=True, num=lam_pts)\n",
    "lam_inv = np.linspace(1/norlam_min, 1/norlam_max, num=lam_pts, endpoint=True)\n",
    "lams = 1.0/lam_inv\n",
    "\n",
    "\n",
    "\n",
    "Rs = np.zeros(lams.size)\n",
    "Rs = mkq.tmm_eval_wsweep(qx, 0, lam_low=norlam_min, lam_high=norlam_max, lam_pts=lam_pts)\n",
    "\n",
    "vgdr = 100*np.ones(lam_pts)\n",
    "for idx, lam in enumerate(lams):\n",
    "    if lams[idx] <= 1:\n",
    "        cuts = Rs[np.logical_and(lams >= lam, lams <= 2*lam)]\n",
    "        #print(cuts)\n",
    "        vgdr[idx] = np.mean(cuts)\n",
    "#         print(vgdr)\n",
    "\n",
    "\n",
    "\n",
    "#400.0/np.amin(vgdr)\n",
    "\n",
    "\n",
    "# lams\n",
    "c, cthick = mkq.vgdr2_eval_wsweep(qx = qx, inc_ang=0, lam_low=norlam_min, lam_high=norlam_max, lam_pts=lam_pts)\n",
    "\n",
    "\n",
    "# num_layers = 10\n",
    "# # qx = np.random.uniform(0,1, (num_layers,2))\n",
    "# # qx[:,0] = qx[:,0]/np.sum(qx[:,0])\n",
    "\n",
    "norlam_min = 400.0/cthick   # this is lam/d\n",
    "norlam_max = 800.0/cthick       # this is lam/d\n",
    "lam_pts = 256\n",
    "\n",
    "lam_inv = np.linspace(1/norlam_min, 1/norlam_max, num=lam_pts, endpoint=True)\n",
    "lams = (1.0/lam_inv)\n",
    "lam_ac = np.linspace(400, 800,  num=lam_pts, endpoint=True)\n",
    "lams2 = np.linspace(400, 800, endpoint=True, num = lam_pts)\n",
    "#lams = np.flip(lams, axis=0)\n",
    "\n",
    "# # # for thick in [200, 400, 800]:\n",
    "# # qx[:,1] = mkq.digitize_qx(qx[:,1], dlevels=2)\n",
    "# d_x, n_x = mkq.make_nxdx(qx=qx)\n",
    "\n",
    "Rs2 = mkq.tmm_eval_wsweep(qx = qx, inc_ang=0, lam_low=norlam_min, lam_high=norlam_max, lam_pts=lam_pts)\n",
    "# lams = cthick*lams\n",
    "meansl = np.mean(Rs2)*np.ones_like(lams)\n",
    "# #     #plt.plot(1.0/lams, Rs)\n",
    "# plt.plot( Rs)\n",
    "# #plt.plot(lams, meansl)\n",
    "plt.ylim([0,1])\n",
    "#plt.xlim([0.75,1.1])\n",
    "    \n",
    "\n",
    "plt.plot(lams2, Rs2)\n",
    "plt.plot(lams2, meansl)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-12-19T16:54:36.746337Z",
     "start_time": "2018-12-19T16:54:36.580467Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.23265694156412065,\n",
       " 637.9084967320262,\n",
       " array([  0.        ,  93.94781328,  29.2463678 ,   6.49963796,\n",
       "         72.73994762,  34.09921199,  12.7687302 , 388.60678789]))"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c, cthick, cthick*qx/(np.sum(qx))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-12-10T13:42:38.721768Z",
     "start_time": "2018-12-10T13:42:37.792591Z"
    }
   },
   "outputs": [],
   "source": [
    "%reload_ext autoreload\n",
    "%autoreload 2\n",
    "import numpy as np\n",
    "import makeqx as mkq\n",
    "import matplotlib.pyplot as plt\n",
    "from jupyterthemes import jtplot\n",
    "# currently installed theme will be used to\n",
    "# set plot style if no arguments provided\n",
    "jtplot.style()\n",
    "\n",
    "\n",
    "num_layers = 50\n",
    "\n",
    "layers = np.arange(num_layers)\n",
    "nlr = np.random.uniform(0,1, num_layers)\n",
    "#nlr = mkq.make_qx(num_layers)\n",
    "#nlr = mkq.digitize_qx(nlr, dlevels=2)\n",
    "#plt.bar(layers, nlr)\n",
    "#plt.ylim([0,1])\n",
    "\n",
    "lam_inv = np.linspace(0.2, 2, num=50, endpoint=True)\n",
    "lams = (1.0/lam_inv)\n",
    "lams = np.flip(lams, axis=0)\n",
    "\n",
    "for thick in [200, 400, 800]:\n",
    "    d_x, n_x = mkq.make_nxdx(qx=nlr, cthick=thick, wavelen=550)\n",
    "    Rs = mkq.tmm_eval_wsweep(qx = nlr, cthick=thick, inc_ang=0, lam_low = 0.5*thick, lam_high=5*thick, lam_pts=50)\n",
    "    #plt.plot(1.0/lams, Rs)\n",
    "    plt.plot(lams, Rs)\n",
    "    plt.ylim([0,100])\n",
    "\n",
    "\n",
    "vgdr = np.zeros(25)\n",
    "for idx, lam in enumerate(lams):\n",
    "    if idx < 25:\n",
    "        cuts = Rs[np.logical_and(lams >= lam, lams <= 2*lam)]\n",
    "        vgdr[idx] = np.mean(cuts)\n",
    "# plt.plot(vgdr)\n",
    "# plt.ylim([0,100])\n",
    "\n",
    "size =  400/lams[np.argmin(vgdr)]\n",
    "val = np.amin(vgdr)\n",
    "size, val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-12-11T14:38:47.524293Z",
     "start_time": "2018-12-11T14:38:46.239283Z"
    }
   },
   "outputs": [],
   "source": [
    "#\\][;/++plt.plot(Rs)\n",
    "# for fd in zip(lams, Rs):\n",
    "#     print(fd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-12-11T14:57:18.592232Z",
     "start_time": "2018-12-11T14:57:17.305396Z"
    }
   },
   "outputs": [],
   "source": [
    "c "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-12-11T14:57:29.843041Z",
     "start_time": "2018-12-11T14:57:28.560532Z"
    }
   },
   "outputs": [],
   "source": [
    "np.mean(Rs2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-12-11T14:40:31.683888Z",
     "start_time": "2018-12-11T14:40:30.373859Z"
    }
   },
   "outputs": [],
   "source": [
    "vgdr[178]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-12-11T14:41:59.272354Z",
     "start_time": "2018-12-11T14:41:57.982383Z"
    }
   },
   "outputs": [],
   "source": [
    "400.0/lams[178]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}