aistats_2022_bilinear.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<torch._C.Generator at 0x7f5c88604eb0>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from itertools import product\n",
    "import matplotlib\n",
    "import matplotlib.patches as mpatches\n",
    "import matplotlib.lines as mlines\n",
    "from matplotlib import pyplot as plt\n",
    "from collections import defaultdict\n",
    "import numpy as np\n",
    "import torch\n",
    "from torch import nn, autograd\n",
    "import optimizers\n",
    "import math\n",
    "from tqdm import tqdm\n",
    "np.random.seed(0)\n",
    "torch.manual_seed(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Configurations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A conditioning: 2.5051719389033593, smax: 1.4629746556141772, smin: 0.583981735103817\n",
      "gd_alt_ro: 1.0\n",
      "gd_sim_ro: 1.010644818138249\n",
      "pp_sim_ro: 0.998299175731378\n",
      "eg_sim_ro: 0.9434110166401143\n"
     ]
    }
   ],
   "source": [
    "COLORS = ['#7DB1F1', '#FC82AB', '#7FD97B', '#FFC000', '#000000']\n",
    "\n",
    "# game settings\n",
    "LR = 0.1\n",
    "N = 1000\n",
    "EPS = 0.01\n",
    "\n",
    "# sample a matrix\n",
    "A = (torch.eye(N, dtype=torch.float64) + EPS * torch.randn(N, N, dtype=torch.float64)).cuda()\n",
    "smax = A.svd()[1][0].cpu().numpy()\n",
    "smin = A.svd()[1][-1].cpu().numpy()\n",
    "sis = A.svd()[1].cpu().numpy()\n",
    "print('A conditioning: {}, smax: {}, smin: {}'.format(smax/smin, smax, smin))\n",
    "\n",
    "\n",
    "def spow(z, k):\n",
    "    return z ** k\n",
    "\n",
    "def gd_alt_radius(k, alpha):\n",
    "    return abs(1- alpha + alpha * spow(\n",
    "        1 - LR**2 * sis**2 / 2 + LR * sis * np.sqrt(1- LR**2 * sis**2 / 4) * 1j\n",
    "    , k)).max()\n",
    "\n",
    "\n",
    "def gd_sim_radius(k, alpha):\n",
    "    return abs(1 - alpha + alpha * spow(1 + LR * sis * 1j, k)).max()\n",
    "\n",
    "\n",
    "def pp_sim_radius(k, alpha):\n",
    "    return abs(1 - alpha + alpha * spow((1+LR*sis*1j)/(1+LR**2*sis**2), k)).max()\n",
    "\n",
    "\n",
    "def eg_sim_radius(k, alpha):\n",
    "    return abs(1 - alpha + alpha * spow(1- LR*sis + LR*sis*1j, k)).max()\n",
    "\n",
    "\n",
    "gd_alt_ro = gd_alt_radius(1, 1)\n",
    "gd_sim_ro = gd_sim_radius(1, 1)\n",
    "pp_sim_ro = pp_sim_radius(1, 1)\n",
    "eg_sim_ro = eg_sim_radius(1, 1)\n",
    "print(f'gd_alt_ro: {gd_alt_ro}')\n",
    "print(f'gd_sim_ro: {gd_sim_ro}')\n",
    "print(f'pp_sim_ro: {pp_sim_ro}')\n",
    "print(f'eg_sim_ro: {eg_sim_ro}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GD_SIM BETAS: (6.09, 37.514)\n",
      "PP_SIM BETAS: (1, 97.95)\n",
      "EG_SIM BETAS: (20.31, 76.16)\n"
     ]
    }
   ],
   "source": [
    "GD_SIM_BETA1 = 6.09\n",
    "GD_SIM_BETA2 = 37.514\n",
    "print(f'GD_SIM BETAS: {GD_SIM_BETA1, GD_SIM_BETA2}')\n",
    "\n",
    "PP_SIM_BETA1 = 1\n",
    "PP_SIM_BETA2 = 97.95\n",
    "print(f'PP_SIM BETAS: {PP_SIM_BETA1, PP_SIM_BETA2}')\n",
    "\n",
    "EG_SIM_BETA1 = 20.31\n",
    "EG_SIM_BETA2 = 76.16\n",
    "print(f'EG_SIM BETAS: {EG_SIM_BETA1, EG_SIM_BETA2}')\n",
    "\n",
    "C = {\n",
    "    'GD_SIM': lambda k: min(\n",
    "        (2 * np.cos(np.angle(1-spow(1+LR*sis*1j, k))) / abs(1-spow(1+LR*sis*1j, k)))\n",
    "    ),\n",
    "    'PP_SIM': lambda k: max(\n",
    "        (\n",
    "            np.cos(np.angle(1-((1+LR*sis*1j)/(1+LR**2*sis**2))**k)) - \n",
    "            np.sqrt(\n",
    "                1/(1+LR**2*smin**2)**k - \n",
    "                np.sin(np.angle(\n",
    "                    1-spow((1+LR*sis*1j)/(1+LR**2*sis**2), k)\n",
    "                )) ** 2\n",
    "            )\n",
    "        ) /  abs(1-((1+LR*sis*1j)/(1+LR**2*sis**2))**k)\n",
    "    ),\n",
    "    'EG_SIM': lambda k:  max(\n",
    "        (\n",
    "            np.cos(np.angle(1-spow(1-LR*sis+LR*sis*1j, k))) -\n",
    "            np.sqrt(\n",
    "                (1-2*LR*smin+2*LR**2*smin**2)**k -\n",
    "                np.sin(np.angle(\n",
    "                    1-spow(1-LR*sis+LR*sis*1j, k)\n",
    "                )) ** 2\n",
    "            )\n",
    "        ) / abs(1-spow(1-LR*sis+LR*sis*1j, k))\n",
    "    )\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GD_ALT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-4-f6fee4ad5edd>:14: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3.  Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading'].  This will become an error two minor releases later.\n",
      "  plt.pcolormesh(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rmin: 0.6418117802935815, rmax: 0.9999999914846447\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASUAAAErCAYAAACYbZUHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAbMUlEQVR4nO3dfXhU5ZkG8PtJJglJCCEgEj4kAYQEFhMgikCBWtqK2BVQLF+7ULXCLopVqyu6uqIrrkutu3IJgogfVbblUgqKWiyIWlfQFQwkBEiAliAYAiQgCZAAyTz7x5zAGCZhgMx73pncv+vKxXyceedJznDnnJNz3kdUFUREtohyuwAiIn8MJSKyCkOJiKzCUCIiqzCUiMgqDCUisorH7QIuVG5u7giPxzNLVVPBUCUKN14RKa2pqXmyf//+fw60gITTeUq5ubkj4uLi5qWnp5+Kj4+vjoqKCp/iiQher1eqqqpaFBcXx548eXJGoGAKqy0Nj8czKz09/VRiYmIVA4ko/ERFRWliYmJVenr6KY/HMyvgMqaLuhSqmhofH1/tdh1EdGni4+OrnUMw5wirUAIQxS0kovDn/D8OmD/hFkpEFOEYSkRkFYYSEVmFoRSmvF4vMjIyer/wwgttL+R1kydP7jJu3Li0UNUFAAUFBXGTJk1Ky8jI6B0dHZ0zYMCAjFC+X6Rrbus67E6eDOSZVRXZVafNfy/xMah5ZGSrPNPvCwCvvPJKSkVFRfS0adMOX8jrHnvssdKsrKw+BQUFpX369DkZito2b94c//HHHyf369fv2OnTp6Upx6797cpsVJ0y/7mNj62JfnAU13U9oVjXEbGl5EYgufm+ALBgwYL2t9566+G4uLgL+mtkRkbGqZycnGNz585tdynvv3Xr1riGnps4ceJ3paWl+atWrfpbz549qy7lfc7hRiC5+b5wf103JhTrOiJCKVwdOHAgety4cWmtW7fuGxsb279///6ZhYWFsed7XUFBQdymTZsSJ0yY8L3fnAcOHIieNGlSWkpKSnZycnLf6dOndwKAnJycjKeffvryuuVGjx59ZMWKFW1ra2svqN4TJ07IggUL2gwcOLDnVVdd1aeh5aKjoy9o3ObA1nW9Zs2axOeff75tcXFxzMV8X6FY1wwll+zfv98zcODAzH379sXNmzev+Lnnntuza9euFlOmTOl6vteuWrUqKT4+3jtw4MAzv5nKysqiBw4cmJmfn58wb9684unTp5cuXLgwdc6cOe2++eabuPvvv/9Q3bLDhg07Vl5e7vnqq6/ig6l1/fr18VOmTOmSmpqafe+996anpKTULF26dOfFfefNj83ruqKiInr27Nmdu3btmtWrV6/eM2bM6PThhx+2PH36dNN88xchIo4phaNf/epXnVu2bOn97LPPdng8vtVQUlISM2fOnE5Hjx6NSk5O9jb02tzc3MRu3bpV+/+WevTRRzuUlZXFrF+/vqBDhw41AI6+9NJL7WfPnt3p17/+9f6EhIQzm/5XX311VXR0NNavX584aNCggJvc5eXl0S+//HKbN99887Jt27YlZGZmVj300EMlU6dOLW/fvv2FbWI1czav67Fjx1bcfPPNeZ9//nnCypUrkz/66KPkBQsWpCYmJtYOGTKkYsSIEUfHjBlTkZaWZiylGEouOHToUPSKFSvavvHGG7vqPqQA0L1795OqimPHjkUVFRXF3XXXXV02btxYVP/1Bw4ciGnTpk1N3X2v14u333677aRJk8qcDykAICkpqbaysjL6gQceOOT/+piYGCQlJdWUlpYG3GRftmxZq8mTJ3ePi4vTMWPGlL/88svFgwcPbtpjQ82E7esaAKKiojBs2LATw4YNOwFgf2lpafTSpUtbP/vssx1XrVqVcv/998Pr9X59qT+LYDGUXPDee++1UlWMGjWqwv/x/fv3xyQkJHhTU1Nrli9fnpyVlXUi0OtPnjwp8fHxZ3675uXltThy5IhnxIgR3xuvtrZWpk6derBVq1bn/CaOiYnR6urqgLvvcXFx2qJFC291dXVURUVF9OHDhz1erxdRUdzbv1C2r2t/hYWFse+8807y6tWrk7/88suk2tpaGTx4cMX1119/NPjv+NLxU+aCDRs2JCQlJdW2bNnye39NWbNmTfKwYcOORkdHY9OmTQk5OTnHAWDWrFntr7vuuisPHDgQDQApKSm1FRUVZ36h7Ny5Mxbw/fb1GyuxpKQktl+/fgE/7JWVldH+v4H93XTTTZUlJSX5c+fOLS4tLY0ZNWpUzy5dulz14IMPdtixY8d5D87SWbav64KCgrg77rjjivT09D69evW66oUXXkjt2LHjqcWLF/+trKxs87p163bOmjXr4KX/JILHUHLBli1bEiorK6PLysrOHChYs2ZN4rp161pNnTq1DADy8/MTMjMzT44cObLboUOHPGvXrt1VdyynR48e1fv27TsTDnW7BYcOHfIAvk38hx9+uDMABJovq6SkxFNdXR2VkZHR4IwL8fHxOm3atCPr16/fuX379i1jx44tX7JkSbvMzMyrBg8e3HPBggVtmuanEdlsX9efffZZYmFhYYvbbrvtUG5u7tZvv/12y5IlS76ZNGnS0aSkpAaPdYUSQ8kFhYWFCW3btj3985//vOuyZctaPf3005ePHTu2x7hx48puueWWCq/Xi6Kiovjp06enjRs37siLL774rf+BzqFDhx7bv39/bElJiQcABg8efDwuLk5nzpzZedmyZa3GjRuX/t1333nS0tJOLl26tE1RUdH3tm4+//zzBBHB8OHDjwVTb0ZGxqm5c+eW7Nu3L//3v//9rsTExNp77rknvaHlKysro1577bWU1157LaW0tDT28OHDnrr7lZWVzeozZ/u6Hj9+/NH58+fvHT169FEA2LRpU4tAXw19f6FY1zymZNiePXtiysvLPUuWLNn1xhtvXDZ58uTuSUlJtXfeeefB3/zmNyWA78TEtLS0kwkJCd6amnO3um+88cbK5OTk2hUrVrS6++67D6emptbOmzdv96xZszpPnDjxyp49e1Z98MEHO9esWZM0c+bMLm+++eaJ2bNnH6h7/apVq5KvueaaytTU1Av6K5rH48GECROOTpgw4ejevXsb/OyUlJR47rjjjm7+j9XdHzx48JaMjIxTF/K+4Soc1vXvfve7lAceeOC8l6KoasAD3aFY12E1HW5eXl5xdnZ2Wf3Hw+kyk7feeqvV+PHjexw8eHBzu3btAn5QFi9enPLJJ58kPffcc98OHTo086WXXioePnz4cf9lbr/99it2794d9+mnn+66kPevqalBp06dsmbNmrXvrrvuuqDLFmwQTpeZcF03Li8v77Ls7Oz0+o9HxJaSW9efXYxNmzYldOzY8VRDH9K6Zfr27XsiNTW1dunSpX8dO3bslatXr97Rs2fPM791HnvssdI+ffr0yc/Pj8vKygr6uqZXX301pUWLFt6pU6da9yENhlvXn10MruuLExGhFE62bNmS0KtXr4B/Jakzf/78b+tu5+TkVBcXFxfUX6Z79+6n586dW7x3796YC/mgqqosXLiwOCbmoq4qoAvAdX1xImL3jYjCT0O7b83qLyFEZD+GEhFZhaFERFZhKBGRVcItlLxer7dJp1clIvOc/8cBL2MJq1ASkdKqqqoGT3knovBQVVXVQkRKAz0XVqFUU1PzZHFxcezx48fjucVEFH68Xq8cP348vri4OLampubJQMuE1XlKAJCbmzvC4/HMcvqQh1WoEhG8IlJaU1PzZP/+/f8caIGwCyUiimzc0iAiqzCUiMgqEXFBrohwH5QofJSpaoMNMiMilADg3969+LnNnxqdzDGacAwbauAYTT9GE9awp7FluPtGRFZhKBGRVRhKRGQVhhIRWYWhRERWYSgRkVWMh5KI3CAiRSKyS0QeDvB8ioisEJF8EflKRPqYrpGI3GM0lEQkGsB8ACMB9AYwUUR611vsXwFsVtUsAFMAzDVZIxG5y/SW0gAAu1T1b6p6CsBSAKPrLdMbwFoAUNVCAOki0t5smUTkFtOh1AnAXr/7+5zH/OUBuAUARGQAgDQAnY1UR0SuMx1KgSZmq3/d2n8CSBGRzQDuAbAJwDlN1kVkmohsFJGNTV4lEbnG6HxKIjIIwBOqOsK5/wgAqOozDSwvAHYDyFLVikbG5QW5ROHja1W9uqEnTW8pbQDQQ0S6ikgsgAkAVvovICKtnecA4E4AnzUWSEQUWYzOEqCqNSIyA8CfAUQDeFVVt4rIPzvPLwTQC8AbIlILYBuAXwYzdrhfgR1JY9hQA8do+jGasobGGJ+6RFX/BOBP9R5b6Hf7CwA9TNdFRHbgGd1EZBWGEhFZhaFERFZhKBGRVRhKRGQVhhIRWYWhRERWYSgRkVUYSkRkFYYSEVnF6CwBocJZAojCilWzBBARNcr4BbmhEu5XYEfSGDbUwDGafgxTswRwS4mIrMJQIiKr2Nj3LVlE3hORPBHZKiK3m66RiNxjY9+3uwFsU9VsANcBeM5velwiinA29n1TAElO04CWAA4jQDcTIopMNvZ9mwffPN0lALYAuFdVvfUHYosloshkY9+3EQA2A+gIoC+AeSLS6pwXqS5S1asbOwmLiMKP6VDaB+AKv/ud4dsi8nc7gOXqswu+vm+ZhuojIpdZ1/cNwDcAfgwAItIeQAaAvxmtkohcY2Pft6cAvC4iW+Db3ZupqmUm6yQi99jY960EwPWm6yIiO3CWACIyjbMEEFH44CwBsOMK7Egaw4YaOEbTj8FZAoioWWIoEZFVGEpEZBWGEhFZhaFERFZhKBGRVRhKRGQVhhIRWYWhRERWYSgRkVUYSkRkFRtbLP2LiGx2vgpEpFZE2piuk4jcYXTqEqfF0g4AP4VvatwNACaq6rYGlr8JwP2qOvw843LqEqLwYdXUJcG0WPI3EcAfjFRGRFYwPXVJoBZL1wZaUEQSANwAYEYwA4f7tBCRNIYNNXCMph/D1NQlpkMpmBZLdW4CsE5VDwccSGQagGlNVRgR2cF0KAXTYqnOBDSy66aqiwAsAnhMiSiS2NhiCSKSDOCHAN41XB8RuczGFksAcDOA1ap63GR9ROQ+61osOfdfB/C6uaqIyBY8o5uIrMJQIiKrMJSIyCoMJSKyCkOJiKzCUCIiqxidJSBUeEY3UVixapYAIqJGGT95MlTC/QrsSBrDhho4RtOPYWqWAG4pEZFVGEpEZBWGEhFZhaFERFZhKBGRVRhKRGQV6/q+Octc5/R92yoifzFdIxG5x+h5Sk7ft/nw6/smIiv9+76JSGsALwK4QVW/EZHLTdZIRO6yse/bJADLVfUbAFDVg4ZrJCIXmQ6lQH3fOtVbpieAFBH5VES+FpEpxqojItfZ2PfNAyAHwI8BxAP4QkS+VNUd3xuIfd+IIpLRWQJEZBCAJ1R1hHP/EQBQ1Wf8lnkYQAtVfcK5/wqAD1X17UbG5SwBROHDqlkCgun79i6AoSLicVp3Xwtgu+E6icgl1vV9U9XtIvIhgHwAXgCLVbXgfGOH+xXYkTSGDTVwjKYfw9QsAbb2fXsWwLMm6yIiO/CMbiKyCkOJiKzCUCIiqzCUiMgqDCUisgpDiYiswlAiIqswlIjIKgwlIrIKQ4mIrGJ0loBQ4SwBRGHFqlkCiIgaZfyC3FAJ9yuwI2kMG2rgGE0/hqlZArilRERWYSgRkVWCCiURuduZKfKSna/vm9Pz7ajT922ziDzeFO9LROEh2C2liQB2i8h9IhLv/4SIdBGRB4MZxK/v20gAvQFMFJHeARb9X1Xt63z9e5A1ElEECCqUVHUIgKcAPA6gWEQeEZEbReRnAP7DeTwYwfR9I6JmLKi/vonIMwAeAlAOoALAVADp8LVH2glgRpDvF6jv27UBlhskInkASgA8qKpbA9TEFktEESjY3bdpAB5V1ctVtYeqdgMwBMAXAGoA/CXIcYLp+5YLIE1VswG8AOCdQAOp6iJVvbqxk7CIKPwEG0peAN/rKKKq653dug0AGuzJVs8+AFf43e8M39aQ/7gVqnrMuf0nADEiclmQ4xNRmAs2lJYDeFxEUgI8txRAVpDjnLfvm4ikiog4twc4NZYHOT4Rhblgz+h+GL5dtb+KyCIAHwHYA6ANgH8F8G0wgwTT9w3ArQCmi0gNgCoAEzQSLtAjoqAEFUqqekRErgXwGIDb4TvorfAdI6oGMCXYNzxf3zdVnQdgXrDjEVFkuahZAkQkA0BXALXwXfF7uKkLu8B6uCVFFD4anSXgoi7IVdUiAEUXXRIRUQM4SwDsuAI7ksawoQaO0fRjcJYAImqWGEpEZBWGEhFZhaFERFZhKBGRVRhKRGQVhhIRWYWhRERWYSgRkVUYSkRkFYYSEVnFeCidr8WS33LXiEitiNxqsj4ictdFTV1y0W/ma7G0A8BP4ZsadwOAiaq6LcBya+Cbq+lVVV12nnE5dQlR+Gh06hLTW0rBtli6B8AfARw0WRwRuc/01CXnbbEkIp0A3AxgOIBrgh043KeFiKQxbKiBYzT9GKamLjEdSsG0WHoewExVrXX6BwQeiH3fiCKS6VA6b4slAFcDWOoE0mUAbhSRGlV9x38hVV0EYBHAY0pEkcR0KJ1psQRfB5QJACb5L6CqXetui8jrAN6vH0hEFLmMhlKQLZaIqBkzPkf3+Vos1Xv8NhM1EZE9eEY3EVmFoUREVmEoEZFVGEpEZBWGEhFZhaFERFYxOktAqPCMbqKwYtUsAUREjTJ+8mSohPsV2JE0hg01cIymH8PULAHcUiIiqzCUiMgqDCUisgpDiYiswlAiIqswlIjIKtb1fROR0SKSLyKbRWSjiAwxXSMRucfoeUpOP7f58Ov7JiIr6/V9WwtgpaqqiGQBeAtApsk6icg91vV9U9Vjevbal0Sc2+2EiCKY6VAK1PetU/2FRORmESkE8AGAOwzVRkQWMB1KwfR9g6quUNVMAGMAPBVwIJFpzjGnjU1bIhG5yegsASIyCMATqjrCuf8IAKjqM428ZjeAa1S1rJFluItHFD6smiXgTN83EYmFr+/bSv8FRORKcTpRikh/ALEAyg3XSUQusbHv21gAU0TkNIAqAOM1iM25cL8CO5LGsKEGjtH0Y5iaJcC6vm+qOgfAHNN1EZEdeEY3EVmFoUREVmEoEZFVGEpEZBWGEhFZhaFERFZhKBGRVRhKRGQVhhIRWYWhRERWMTpLQKhwlgCisGLVLAFERI0yfkFuqIT7FdiRNIYNNXCMph/D1CwB3FIiIqswlIjIKjb2ffsHp+9bvoisF5Fs0zUSkXuMhpJf37eRAHoDmCgivestthvAD1U1C76mAYtM1khE7rKx79t6VT3i3P0SQGfDNRKRi6zs++bnlwBWBXqCLZaIIpPpUwKC6vsGACLyI/hCaUig51V1EZxdO548SRQ5TIfSPgBX+N3vDKCk/kIikgVgMYCRqsr2SkTNiI1937oAWA5gsqruMFwfEbnMxr5vjwNoC+BFpydlTWPXyRBRZLGx79udAO40XRcR2YGzBBCRaZwlgIjCB2cJgB1XYEfSGDbUwDGafgzOEkBEzRJDiYiswlAiIqswlIjIKgwlIrIKQ4mIrMJQIiKrMJSIyCoMJSKyCkOJiKzCUCIiq9jYYilTRL4QkZMi8qDp+ojIXUanLnFaLO0A8FP4psbdAGCiqm7zW+ZyAGkAxgA4oqq/DWJcTl1CFD6smrokmBZLB1V1A4DThmsjIguYnrokUIula5ti4HCfFiKSxrChBo7R9GOYmrrE2hZL5x1IZBqAaZdWDhHZxsoWS8Fg3zeiyGRdiyUiat6sa7EkIqkANgJoBcArIvcB6K2qFSZrJSJ32NhiqRS+3ToiaoZ4RjcRWYWhRERWYSgRkVUYSkRkFYYSEVmFoUREVjE6S0Co8IxuorDS6CwBxs9TCpEyAHsucYzLnHHcxjrsqgFgHU1dQ1pjT0bEllJTEJGNjaU362ieNbAO8zXwmBIRWYWhRERWYSidtcjtAhys4ywbagBYh7+Q18BjSkRkFW4pEZFVmmUoicirInJQRAr8HmsjImtEZKfzb0qIa7hCRD4Rke0islVE7nWpjhYi8pWI5Dl1POlGHc57RovIJhF538UaikVki4hsFpGNLtbRWkSWiUih8xkZ5MJnI8P5OdR9VYjIfaGuo1mGEoDXAdxQ77GHAaxV1R4A1jr3Q6kGwAOq2gvAQAB3i0hvF+o4CWC4qmYD6AvgBhEZ6EIdAHAvgO1+992oAQB+pKp9/f707UYdcwF8qKqZALLh+7kYrUNVi5yfQ18AOQBOAFgR8jpUtVl+AUgHUOB3vwhAB+d2BwBFhut5F75+eK7VASABQC58HWaM1gHfxH5rAQwH8L5b6wRAMYDL6j1m+mfRCsBuOMd83aqj3ntfD2CdiTqa65ZSIO1VdT8AOP9ebuqNRSQdQD8A/+dGHc5u02YABwGsUVU36ngewEMAvH6PubFOFMBqEfna6ZjjRh3dABwC8JqzO7tYRBJdqMPfBAB/cG6HtA6GkstEpCWAPwK4T12ah1xVa9W3id4ZwAAR6WPy/UXk7wEcVNWvTb5vA36gqv0BjIRvl3qYCzV4APQHsEBV+wE4DnO7rudwmnyMAvC2ifdjKJ11QEQ6AIDz78FQv6GIxMAXSP+jqsvdqqOOqn4H4FP4jreZrOMHAEaJSDF8XZOHi8gSwzUAAFS1xPn3IHzHTwa4UMc+APucLVYAWAZfSLn12RgJIFdVDzj3Q1oHQ+mslQB+4dz+BXzHeEJGRATAKwC2q+p/uVhHOxFp7dyOB/ATAIUm61DVR1S1s6qmw7eb8LGq/qPJGgBARBJFJKnuNnzHUQpM16G+5hl7RSTDeejHALaZrsPPRJzddUPI6zB1oMymL+cHvB/Aafh+K/0SQFv4DrTudP5tE+IahsB3/CIfwGbn60YX6sgCsMmpowDA487jRuvwq+c6nD3Qbfpn0Q1AnvO1FcCjbv0s4PtL6EZnvbwDIMWlOhIAlANI9nsspHXwjG4isgp334jIKgwlIrIKQ4mIrMJQIiKrMJSIyCoMJSKyCkOJrCMiE0TkpHPGOzUzDCWyUTaAbap62u1CyDyGEtmoL3xnmVMzxFAiG2XDd9kNAEBEhotIuYg8LyLR7pVFJjCUyCoi0g6+icM2O/dnAHgfwCOqep+q1rpYHhkQKW27KXL0he9C5W0isgjAWAA/U9VPXK2KjGEokW2yARwGsBy+GQ2vVdVd7pZEJjGUyDbZAATAUAA/YSA1P5y6hKwiIlvgm2mxO3xzK+Wo6iFXiyKjeKCbrOHMBZ0J38Rm/wTf5Plviwi36JsRhhLZ5O/gO6SwRVWrANziPPbfrlZFRjGUyCbZ8HXu+CsAqOoeAOMBTBeR21ysiwziMSUisgq3lIjIKgwlIrIKQ4mIrMJQIiKrMJSIyCoMJSKyCkOJiKzCUCIiqzCUiMgq/w8X7/qcvJR1awAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "GD_ALT_KS = list(range(5, 80, 5))\n",
    "GD_ALT_ALPHAS = list(reversed([x/100 for x in range(5, 95 + 1, 90 // 15)]))\n",
    "GD_ALT_HPS = np.meshgrid(GD_ALT_KS, GD_ALT_ALPHAS)\n",
    "\n",
    "# compute\n",
    "gd_alt_radii = np.array([\n",
    "    gd_alt_radius(k, alpha) for\n",
    "    (k, alpha) in zip(GD_ALT_HPS[0].flatten(), GD_ALT_HPS[1].flatten())\n",
    "]).reshape(GD_ALT_HPS[0].shape)\n",
    "\n",
    "# and plot (matplotlib choses the first color when there is only one value)\n",
    "cmap = matplotlib.colors.ListedColormap([COLORS[0], COLORS[1]])\n",
    "plt.figure(figsize=(4, 4))\n",
    "plt.pcolormesh(\n",
    "    GD_ALT_HPS[0], GD_ALT_HPS[1], gd_alt_radii < 1, \n",
    "    edgecolors='k',\n",
    "    cmap=cmap,\n",
    ")\n",
    "plt.xlabel(r'$k$', fontsize=15)\n",
    "plt.ylabel(r'$\\alpha$', fontsize=15)\n",
    "\n",
    "#\n",
    "handles, labels = plt.gca().get_legend_handles_labels()\n",
    "patch_pos = mpatches.Patch(\n",
    "    color=COLORS[0], \n",
    "    label=r'$\\rho_k(\\alpha)<1$'\n",
    ")\n",
    "patch_neg = mpatches.Patch(\n",
    "    color=COLORS[1],\n",
    "    label=r'$\\rho_k(\\alpha)\\geq1$'\n",
    ")\n",
    "line_suff = mlines.Line2D([0], [0], color=COLORS[3], label='Thm. 6', linewidth=3)\n",
    "\n",
    "handles.append(patch_pos) \n",
    "handles.append(patch_neg) \n",
    "#handles.append(line_suff)\n",
    "\n",
    "# plot the legend\n",
    "plt.legend(\n",
    "    handles=handles, loc='upper center', \n",
    "    bbox_to_anchor=(-0.01, 1.20, 1.0, 0), \n",
    "    ncol=2, fontsize=15, columnspacing=1,\n",
    ")\n",
    "plt.tight_layout(pad=0.1)\n",
    "plt.savefig('./bilinear_sweep_gd_alt.png', format='png', dpi=200)\n",
    "print(f'rmin: {gd_alt_radii.min()}, rmax: {gd_alt_radii.max()}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GD_SIM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-5-79eac1bd85ae>:14: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3.  Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading'].  This will become an error two minor releases later.\n",
      "  plt.pcolormesh(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "beta1: 6.09, beta2: 37.514\n",
      "=> c: [0.7912922703296353, 0.8805040334283156, 0.8993097944210748, 0.8930392790987423, 0.8603330786306569, 0.7689185514601906, 0.4629918667201008, -1.1749286262903889, -3.638956650736324, -4.137990641887426, -4.638074932663175, -5.139347824807506, -5.641375267247767, -6.143891516596959, -6.646947012114167]\n",
      "=> a: [0.95, 0.89, 0.83, 0.77, 0.71, 0.65, 0.59, 0.53, 0.47, 0.41, 0.35, 0.29, 0.23, 0.17, 0.11, 0.05]\n",
      "=> k: [5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASUAAAErCAYAAACYbZUHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAbOUlEQVR4nO3dfXhU5ZkG8PtJJglJCCEgEj4kAYQEFhMgikCRWtqK2BVQlK9dqFphF8Wq1RWtrtGK61LrrlyCIOJHlW25lIKiFgui1hV0BQOE8BGgJQiGAAlIAoSPZJ79Y05gDJMwCZn3vDO5f9eVy8zk5J2HzHjnzMk5zyOqCiIiW0S5XQARkT+GEhFZhaFERFZhKBGRVRhKRGQVhhIRWcXjdgENlZeXN9zj8eSqaioYqkThxisiJVVVVU/279//L4E2kHA6TykvL294XFzcnPT09NPx8fEno6Kiwqd4IoLX65XKysoWRUVFsadOnZoeKJjCak/D4/Hkpqenn05MTKxkIBGFn6ioKE1MTKxMT08/7fF4cgNuY7qoi6GqqfHx8SfdroOILk58fPxJ5xDMecIqlABEcQ+JKPw5/x8HzJ9wCyUiinAMJSKyCkOJiKzCUApTXq8XGRkZvV944YW2Dfm+SZMmdRk7dmxaqOoCgIKCgriJEyemZWRk9I6Ojs4ZMGBARigfL9I1t+c67E6eDOSZFeXZlWfM/1viY1D1yIhWm0w/LgC88sorKeXl5dFTp0493JDve+yxx0qysrL6FBQUlPTp0+dUKGrbuHFj/Mcff5zcr1+/Y2fOnJGmXLv6d8uzUXna/Os2PrYq+sGRfK5rCcVzHRF7Sm4EkpuPCwDz5s1rf8sttxyOi4tr0F8jMzIyTufk5BybPXt2u4t5/C1btsTV9bUJEyZ8V1JSkr9ixYq/9+zZs/JiHuc8bgSSm48L95/r+oTiuY6IUApXBw4ciB47dmxa69at+8bGxvbv379/5vbt22Mv9H0FBQVxGzZsSBw/fvz3fnMeOHAgeuLEiWkpKSnZycnJfadNm9YJAHJycjKefvrpS2u2GzVq1JFly5a1ra6ublC9J06ckHnz5rUZOHBgzyuuuKJPXdtFR0c3aN3mwNbnetWqVYnPP/9826KiopjG/LtC8VwzlFyyf/9+z8CBAzP37dsXN2fOnKLnnntuz65du1pMnjy564W+d8WKFUnx8fHegQMHnv3NVFpaGj1w4MDM/Pz8hDlz5hRNmzatZP78+amzZs1q980338Tdf//9h2q2HTp06LGysjLPV199FR9MrWvXro2fPHlyl9TU1Ox77703PSUlpWrx4sU7G/cvb35sfq7Ly8ujZ86c2blr165ZvXr16j19+vROH374YcszZ840zT++ESLimFI4+uUvf9m5ZcuW3s8++2yHx+N7GoqLi2NmzZrV6ejRo1HJycneur43Ly8vsVu3bif9f0s9+uijHUpLS2PWrl1b0KFDhyoAR1966aX2M2fO7PSrX/1qf0JCwtld/yuvvLIyOjoaa9euTRw0aFDAXe6ysrLol19+uc2bb755ydatWxMyMzMrH3rooeIpU6aUtW/fvmG7WM2czc/1mDFjym+66aZNn3/+ecLy5cuTP/roo+R58+alJiYmVg8ZMqR8+PDhR0ePHl2elpZmLKUYSi44dOhQ9LJly9q+8cYbu2pepADQvXv3U6qKY8eORRUWFsbdddddXdavX19Y+/sPHDgQ06ZNm6qa216vF2+//XbbiRMnljovUgBAUlJSdUVFRfQDDzxwyP/7Y2JikJSUVFVSUhJwl33JkiWtJk2a1D0uLk5Hjx5d9vLLLxcNHjy4aY8NNRO2P9cAEBUVhaFDh54YOnToCQD7S0pKohcvXtz62Wef7bhixYqU+++/H16v9+uL/VkEi6Hkgvfee6+VqmLkyJHl/vfv378/JiEhwZuamlq1dOnS5KysrBOBvv/UqVMSHx9/9rfrpk2bWhw5csQzfPjw761XXV0tU6ZMOdiqVavzfhPHxMToyZMnA759j4uL0xYtWnhPnjwZVV5eHn348GGP1+tFVBTf7TeU7c+1v+3bt8e+8847yStXrkz+8ssvk6qrq2Xw4MHl11133dHg/8UXj68yF6xbty4hKSmpumXLlt/7a8qqVauShw4dejQ6OhobNmxIyMnJOQ4Aubm57a+99trLDxw4EA0AKSkp1eXl5Wd/oezcuTMW8P329Vsrsbi4OLZfv34BX+wVFRXR/r+B/d14440VxcXF+bNnzy4qKSmJGTlyZM8uXbpc8eCDD3bYsWPHBQ/O0jm2P9cFBQVxd9xxx2Xp6el9evXqdcULL7yQ2rFjx9MLFy78e2lp6cY1a9bszM3NPXjxP4ngMZRcsHnz5oSKioro0tLSswcKVq1albhmzZpWU6ZMKQWA/Pz8hMzMzFMjRozodujQIc/q1at31RzL6dGjx8l9+/adDYeatwWHDh3yAL5d/IcffrgzAATql1VcXOw5efJkVEZGRp0dF+Lj43Xq1KlH1q5du3Pbtm2bx4wZU7Zo0aJ2mZmZVwwePLjnvHnz2jTNTyOy2f5cf/bZZ4nbt29vcdtttx3Ky8vb8u23325etGjRNxMnTjyalJRU57GuUGIouWD79u0Jbdu2PXPrrbd2XbJkSaunn3760jFjxvQYO3Zs6c0331zu9XpRWFgYP23atLSxY8ceefHFF7/1P9B5zTXXHNu/f39scXGxBwAGDx58PC4uTmfMmNF5yZIlrcaOHZv+3XffedLS0k4tXry4TWFh4ff2bj7//PMEEcGwYcOOBVNvRkbG6dmzZxfv27cv/w9/+MOuxMTE6nvuuSe9ru0rKiqiXnvttZTXXnstpaSkJPbw4cOemtsVFRXN6jVn+3M9bty4o3Pnzt07atSoowCwYcOGFoE+6vr3heK55jElw/bs2RNTVlbmWbRo0a433njjkkmTJnVPSkqqvvPOOw/+9re/LQZ8JyampaWdSkhI8FZVnb/XfcMNN1QkJydXL1u2rNXdd999ODU1tXrOnDm7c3NzO0+YMOHynj17Vn7wwQc7V61alTRjxowub7755omZM2ceqPn+FStWJF911VUVqampDformsfjwfjx44+OHz/+6N69e+t87RQXF3vuuOOObv731dwePHjw5oyMjNMNedxwFQ7P9e9///uUBx544IKXoqhqwAPdoXiuw6od7qZNm4qys7NLa98fTpeZvPXWW63GjRvX4+DBgxvbtWsX8IWycOHClE8++STpueee+/aaa67JfOmll4qGDRt23H+b22+//bLdu3fHffrpp7sa8vhVVVXo1KlTVm5u7r677rqrQZct2CCcLjPhc12/TZs2XZKdnZ1e+/6I2FNy6/qzxtiwYUNCx44dT9f1Iq3Zpm/fvidSU1OrFy9e/LcxY8ZcvnLlyh09e/Y8+1vnscceK+nTp0+f/Pz8uKysrKCva3r11VdTWrRo4Z0yZYp1L9JguHX9WWPwuW6ciAilcLJ58+aEXr16BfwrSY25c+d+W/N5Tk7OyaKiooLa23Tv3v3M7Nmzi/bu3RvTkBeqqsr8+fOLYmIadVUBNQCf68aJiLdvRBR+6nr71qz+EkJE9mMoEZFVGEpEZBWGEhFZJdxCyev1epu0vSoRmef8fxzwMpawCiURKamsrKzzlHciCg+VlZUtRKQk0NfCKpSqqqqeLCoqij1+/Hg895iIwo/X65Xjx4/HFxUVxVZVVT0ZaJuwOk8JAPLy8oZ7PJ5cZw55WIUqEcErIiVVVVVP9u/f/y+BNgi7UCKiyMY9DSKyCkOJiKwSERfkigjfgxKFj1JVrXNAZkSEEgBU/ebtRn+v5/FbuUYTrmFDDVyj6ddowhr21LcN374RkVUYSkRkFYYSEVmFoUREVmEoEZFVGEpEZBXjoSQi14tIoYjsEpGHA3w9RUSWiUi+iHwlIn1M10hE7jEaSiISDWAugBEAegOYICK9a232awAbVTULwGQAs03WSETuMr2nNADALlX9u6qeBrAYwKha2/QGsBoAVHU7gHQRaW+2TCJyi+lQ6gRgr9/tfc59/jYBuBkARGQAgDQAnY1UR0SuMx1KgRqz1b5u7T8BpIjIRgD3ANgA4Lwh6yIyVUTWi8j6Jq+SiFxjtJ+SiAwC8ISqDnduPwIAqvpMHdsLgN0AslS1vJ51eUEuUfj4WlWvrOuLpveU1gHoISJdRSQWwHgAy/03EJHWztcA4E4An9UXSEQUWYx2CVDVKhGZDuAvAKIBvKqqW0TkX52vzwfQC8AbIlINYCuAXwSzti1XYP/7u0cbvcZTo5KbbI0IuZqca1i0RlPWUO82jV69kVT1zwD+XOu++X6ffwGgh+m6iMgOPKObiKzCUCIiqzCUiMgqDCUisgpDiYiswlAiIqswlIjIKgwlIrIKQ4mIrMJQIiKrGO0SECrsEkAUVqzqEkBEVC/jF+SGSiRd4W/LGjZcTc417FnDVJcA7ikRkVUYSkRkFRvnviWLyHsisklEtojI7aZrJCL32Dj37W4AW1U1G8C1AJ7za49LRBHOxrlvCiDJGRrQEsBhBJhmQkSRyca5b3Pg69NdDGAzgHtV1Vt7IY5YIopMNs59Gw5gI4COAPoCmCMirc77JtUFqnplfSdhEVH4MR1K+wBc5ne7M3x7RP5uB7BUfXbBN/ct01B9ROQy6+a+AfgGwI8BQETaA8gA8HejVRKRa2yc+/YUgNdFZDN8b/dmqGqpyTqJyD02zn0rBnCd6bqIyA7sEkBEprFLABGFj4jpEmDLlfWRtIYNV5NzDXvWYJcAImqWGEpEZBWGEhFZhaFERFZhKBGRVRhKRGQVhhIRWYWhRERWYSgRkVUYSkRkFYYSEVnFxhFL/yYiG52PAhGpFpE2puskIncYbV3ijFjaAeCn8LXGXQdggqpurWP7GwHcr6rDLrAuW5cQhQ+rWpcEM2LJ3wQAfzRSGRFZwXTrkkAjlq4OtKGIJAC4HsD0YBa2pd0H17j41ieAHa06uEboaqh3m0av3jjBjFiqcSOANap6OOBCIlMBTG2qwojIDqZDKZgRSzXGo563bqq6AMACgMeUiCKJjSOWICLJAH4I4F3D9RGRy2wcsQQANwFYqarHTdZHRO6zbsSSc/t1AK+bq4qIbMEzuonIKgwlIrIKQ4mIrMJQIiKrMJSIyCoMJSKyitEuAaHCM7qJwopVXQKIiOpl/OTJUAn3K+sjaQ12CYjMNUx1CeCeEhFZhaFERFZhKBGRVRhKRGQVhhIRWYWhRERWsW7um7PNtc7cty0i8lfTNRKRe4yep+TMfZsLv7lvIrLcf+6biLQG8CKA61X1GxG51GSNROQuG+e+TQSwVFW/AQBVPWi4RiJykelQCjT3rVOtbXoCSBGRT0XkaxGZbKw6InKdjXPfPAByAPwYQDyAL0TkS1Xd8b2FOPeNKCIZ7RIgIoMAPKGqw53bjwCAqj7jt83DAFqo6hPO7VcAfKiqdV5wwy4BRGHFqi4Bwcx9exfANSLicUZ3Xw1gm+E6icgl1s19U9VtIvIhgHwAXgALVbXgQmuH+5X1kbQGuwRE5hqmugTYOvftWQDPmqyLiOzAM7qJyCoMJSKyCkOJiKzCUCIiqzCUiMgqDCUisgpDiYiswlAiIqswlIjIKgwlIrKK0S4BocIuAURhxaouAURE9TJ+QW6ohPuV9ZG0BrsEROYaproEcE+JiKzCUCIiqwQVSiJyt9Mp8qJdaO6bM/PtqDP3baOIPN4Uj0tE4SHYPaUJAHaLyH0iEu//BRHpIiIPBrOI39y3EQB6A5ggIr0DbPq/qtrX+fhNkDUSUQQIKpRUdQiApwA8DqBIRB4RkRtE5GcA/sO5PxjBzH0jomYsqL++icgzAB4CUAagHMAUAOnwjUfaCWB6kI8XaO7b1QG2GyQimwAUA3hQVbcEqIkjlogiULBv36YCeFRVL1XVHqraDcAQAF8AqALw1yDXCWbuWx6ANFXNBvACgHcCLaSqC1T1yvpOwiKi8BNsKHkBfG+iiKqudd7WrQMQ7IkL+wBc5ne7M3x7Q/7rlqvqMefzPwOIEZFLglyfiMJcsKG0FMDjIpIS4GuLAWQFuc4F576JSKqIiPP5AKfGsiDXJ6IwF+wZ3Q/D91btbyKyAMBHAPYAaAPg1wC+DWaRYOa+AbgFwDQRqQJQCWC8RsIFekQUlKBCSVWPiMjVAB4DcDt8B70VvmNEJwFMDvYBLzT3TVXnAJgT7HpEFFka1SVARDIAdAVQDd8Vv4eburAG1sM9KaLwUW+XgEZdkKuqhQAKG10SEVEd2CUAdlxZH0lrsEtAZK7BLgFE1CwxlIjIKgwlIrIKQ4mIrMJQIiKrMJSIyCoMJSKyCkOJiKzCUCIiqzCUiMgqDCUisorxULrQiCW/7a4SkWoRucVkfUTkrka1Lmn0g/lGLO0A8FP4WuOuAzBBVbcG2G4VfL2aXlXVJRdYl61LiMJHva1LTO8pBTti6R4AfwJw0GRxROQ+061LLjhiSUQ6AbgJwDAAVwW7cLi3+4ikNZqyhnBv9xFJa5hqXWI6lIIZsfQ8gBmqWu3MDwi8EOe+EUUk06F0wRFLAK4EsNgJpEsA3CAiVar6jv9GqroAwAKAx5SIIonpUDo7Ygm+CSjjAUz030BVu9Z8LiKvA3i/diARUeQyGkpBjlgiombMeI/uC41YqnX/bSZqIiJ78IxuIrIKQ4mIrMJQIiKrMJSIyCoMJSKyCkOJiKxitEtAqPCMbqKwYlWXACKiehk/eTJUwv3K+khag10CInMNU10CuKdERFZhKBGRVRhKRGQVhhIRWYWhRERWYSgRkVWsm/smIqNEJF9ENorIehEZYrpGInKP0fOUnHluc+E3901Eltea+7YawHJVVRHJAvAWgEyTdRKRe6yb+6aqx/TctS+JOH/aCRFFMNOhFGjuW6faG4nITSKyHcAHAO4wVBsRWcB0KAUz9w2qukxVMwGMBvBUwIVEpjrHnNY3bYlE5CajXQJEZBCAJ1R1uHP7EQBQ1Wfq+Z7dAK5S1dJ6tuFbPKLwYVWXgLNz30QkFr65b8v9NxCRy8WZRCki/QHEAigzXCcRucTGuW9jAEwWkTMAKgGM0yB258L9yvpIWoNdAiJzDVNdAqyb+6aqswDMMl0XEdmBZ3QTkVUYSkRkFYYSEVmFoUREVmEoEZFVGEpEZBWGEhFZhaFERFZhKBGRVRhKRGQVo10CQoVdAojCilVdAoiI6mX8gtxQCfcr6yNpDXYJiMw1THUJ4J4SEVmFoUREVrFx7ts/OXPf8kVkrYhkm66RiNxjNJT85r6NANAbwAQR6V1rs90AfqiqWfANDVhgskYicpeNc9/WquoR5+aXADobrpGIXGTl3Dc/vwCwItAXOGKJKDKZPiUgqLlvACAiP4IvlIYE+rqqLoDz1o4nTxJFDtOhtA/AZX63OwMorr2RiGQBWAhghKpyvBJRM2Lj3LcuAJYCmKSqOwzXR0Qus3Hu2+MA2gJ40ZlJWVXfdTJEFFlsnPt2J4A7TddFRHZglwAiMo1dAogofLBLAOy4sj6S1mCXgMhcg10CiKhZYigRkVUYSkRkFYYSEVmFoUREVmEoEZFVGEpEZBWGEhFZhaFERFZhKBGRVRhKRGQVG0csZYrIFyJySkQeNF0fEbnLaOsSZ8TSDgA/ha817joAE1R1q982lwJIAzAawBFV/V0Q67J1CVH4sKp1STAjlg6q6joAZwzXRkQWMN26JNCIpaubYuFwb/cRSWuwdUlkrmGqdYm1I5YuuJDIVABTL64cIrKNlSOWgsG5b0SRyboRS0TUvFk3YklEUgGsB9AKgFdE7gPQW1XLTdZKRO6wccRSCXxv64ioGeIZ3URkFYYSEVmFoUREVmEoEZFVGEpEZBWGEhFZxWiXgFDhGd1EYaXeLgHGz1MKkVIAey5yjUucddzGOuyqAWAdTV1DWn1fjIg9paYgIuvrS2/W0TxrYB3ma+AxJSKyCkOJiKzCUDpngdsFOFjHOTbUALAOfyGvgceUiMgq3FMiIqs0y1ASkVdF5KCIFPjd10ZEVonITue/KSGu4TIR+UREtonIFhG516U6WojIVyKyyanjSTfqcB4zWkQ2iMj7LtZQJCKbRWSjiKx3sY7WIrJERLY7r5FBLrw2MpyfQ81HuYjcF+o6mmUoAXgdwPW17nsYwGpV7QFgtXM7lKoAPKCqvQAMBHC3iPR2oY5TAIapajaAvgCuF5GBLtQBAPcC2OZ3240aAOBHqtrX70/fbtQxG8CHqpoJIBu+n4vROlS10Pk59AWQA+AEgGUhr0NVm+UHgHQABX63CwF0cD7vAKDQcD3vwjcPz7U6ACQAyINvwozROuBr7LcawDAA77v1nAAoAnBJrftM/yxaAdgN55ivW3XUeuzrAKwxUUdz3VMKpL2q7gcA57+XmnpgEUkH0A/A/7lRh/O2aSOAgwBWqaobdTwP4CEAXr/73HhOFMBKEfnamZjjRh3dABwC8JrzdnahiCS6UIe/8QD+6Hwe0joYSi4TkZYA/gTgPnWpD7mqVqtvF70zgAEi0sfk44vIPwI4qKpfm3zcOvxAVfsDGAHfW+qhLtTgAdAfwDxV7QfgOMy9dT2PM+RjJIDGD3xrAIbSOQdEpAMAOP89GOoHFJEY+ALpf1R1qVt11FDV7wB8Ct/xNpN1/ADASBEpgm9q8jARWWS4BgCAqhY7/z0I3/GTAS7UsQ/APmePFQCWwBdSbr02RgDIU9UDzu2Q1sFQOmc5gJ87n/8cvmM8ISMiAuAVANtU9b9crKOdiLR2Po8H8BMA203WoaqPqGpnVU2H723Cx6r6zyZrAAARSRSRpJrP4TuOUmC6DvUNz9grIhnOXT8GsNV0HX4m4NxbN4S8DlMHymz6cH7A+wGcge+30i8AtIXvQOtO579tQlzDEPiOX+QD2Oh83OBCHVkANjh1FAB43LnfaB1+9VyLcwe6Tf8sugHY5HxsAfCoWz8L+P4Sut55Xt4BkOJSHQkAygAk+90X0jp4RjcRWYVv34jIKgwlIrIKQ4mIrMJQIiKrMJSIyCoMJSKyCkOJrCMi40XklHPGOzUzDCWyUTaArap6xu1CyDyGEtmoL3xnmVMzxFAiG2XDd9kNAEBEholImYg8LyLR7pVFJjCUyCoi0g6+xmEbndvTAbwP4BFVvU9Vq10sjwyIlLHdFDn6wneh8lYRWQBgDICfqeonrlZFxjCUyDbZAA4DWApfR8OrVXWXuyWRSQwlsk02AAFwDYCfMJCaH7YuIauIyGb4Oi12h6+3Uo6qHnK1KDKKB7rJGk4v6Ez4Gpv9C3zN898WEe7RNyMMJbLJP8B3SGGzqlYCuNm5779drYqMYiiRTbLhm9zxNwBQ1T0AxgGYJiK3uVgXGcRjSkRkFe4pEZFVGEpEZBWGEhFZhaFERFZhKBGRVRhKRGQVhhIRWYWhRERWYSgRkVX+Hz6u+pzRVMLbAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "GD_SIM_KS = list(range(5, 80, 5))\n",
    "GD_SIM_ALPHAS = list(reversed([x/100 for x in range(5, 95 + 1, 90 // 15)]))\n",
    "GD_SIM_HPS = np.meshgrid(GD_SIM_KS, GD_SIM_ALPHAS)\n",
    "\n",
    "gd_sim_improvements = np.array([\n",
    "    1 > gd_sim_radius(k, alpha) for\n",
    "    (k, alpha) in zip(GD_SIM_HPS[0].flatten(), GD_SIM_HPS[1].flatten())\n",
    "]).reshape(GD_SIM_HPS[0].shape)\n",
    "gd_sim_cs = [C['GD_SIM'](k) for k in GD_SIM_KS]\n",
    "\n",
    "edge_colors = [ ]\n",
    "cmap = matplotlib.colors.ListedColormap([COLORS[1], COLORS[0]])\n",
    "plt.figure(figsize=(4, 4))\n",
    "plt.pcolormesh(\n",
    "    GD_SIM_HPS[0], GD_SIM_HPS[1], gd_sim_improvements, \n",
    "    edgecolors='k',\n",
    "    cmap=cmap,\n",
    ")\n",
    "plt.xlabel(r'$k$', fontsize=15)\n",
    "plt.ylabel(r'$\\alpha$', fontsize=15)\n",
    "\n",
    "# where some data has already been plotted to ax\n",
    "handles, labels = plt.gca().get_legend_handles_labels()\n",
    "\n",
    "# manually define a new patch \n",
    "patch_pos = mpatches.Patch(\n",
    "    color=COLORS[0], \n",
    "    label=r'$\\rho_k(\\alpha)<1$'\n",
    ")\n",
    "patch_neg = mpatches.Patch(\n",
    "    color=COLORS[1],\n",
    "    label=r'$\\rho_k(\\alpha)\\geq1$'\n",
    ")\n",
    "line_suff = mlines.Line2D([0], [0], color=COLORS[3], label='Thm. 7', linewidth=3)\n",
    "\n",
    "# handles is a list, so append manual patch\n",
    "handles.append(patch_pos) \n",
    "handles.append(patch_neg) \n",
    "#handles.append(line_suff)\n",
    "\n",
    "# plot the legend\n",
    "plt.legend(\n",
    "    handles=handles, loc='upper center', \n",
    "    bbox_to_anchor=(-0.01, 1.20, 1.0, 0), \n",
    "    ncol=2, fontsize=15, columnspacing=1,\n",
    ")\n",
    "plt.tight_layout(pad=0.1)\n",
    "plt.savefig('./bilinear_sweep_gd_sim.png', format='png', dpi=200)\n",
    "print(f'beta1: {GD_SIM_BETA1}, beta2: {GD_SIM_BETA2}')\n",
    "print(f'=> c: {gd_sim_cs}')\n",
    "print(f'=> a: {GD_SIM_ALPHAS}')\n",
    "print(f'=> k: {GD_SIM_KS}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PP_SIM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-6-c87b9cf73d40>:16: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3.  Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading'].  This will become an error two minor releases later.\n",
      "  plt.pcolormesh(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=> beta1: 1, beta2: 97.95\n",
      "=> c: [0.20136786039250618, 0.08069798598520445, 0.054038407564570276, 0.04402563597496977, 0.1074618793654902, 0.2114666490134871, 0.28242476102585523, 0.3632591686497847, 0.45400520323347043, 0.5549581220320206, 0.6672550989793833, 0.7962543581332878, 0.9999999999999906, 0.9999999999999956, 0.9999999999999807]\n",
      "=> a: [0.95, 0.89, 0.83, 0.77, 0.71, 0.65, 0.59, 0.53, 0.47, 0.41, 0.35, 0.29, 0.23, 0.17, 0.11, 0.05]\n",
      "=> k: [5, 13, 21, 29, 37, 45, 53, 61, 69, 77, 85, 93, 101, 109, 117]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASoAAAErCAYAAABpZs6KAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAby0lEQVR4nO3de3xUhZUH8N8hkzcPAyoR0AQQEqhLFLACq1SxFbGrgCgvC7VbYOWj1IJdFXVFq121Vi0fQCg+VsS1fFaqFR+oqN26ig0okRBegiUIJIEkUMIjPJI5+8dccIyTZJKZ3Htm8vt+PvnITO69cyZcf7lzufccUVUQEVnWxusCiIgaw6AiIvMYVERkHoOKiMxjUBGReQwqIjLP5+WLr1u3brjP55ujqplgaBK1Rn4RKaupqXmwf//+79a3kHh1HdW6deuGJycnz8/Ozj6Rmpp6rE2bNrygi6iV8fv9Ul1dnVJcXJx0/Pjx2+oLK8+OYnw+35zs7OwT6enp1QwpotapTZs2mp6eXp2dnX3C5/PNqXc5N4sKpqqZqampx7x6fSKyIzU19ZhzCigkL88LteGRFBEBgSMrNJBHPIFNROYxqIjIPAYVRcWaNWtSBw4cmON1HWRDtPcHBhVFRX5+flq/fv2Oel0H2RDt/YFBZYDf70dOTk7fefPmdWrKepMmTTpv7NixWS1VV1MUFBSkDRgw4AgAzJkzp/Pll19+/t69exO8risWcX/4LgaVAc8991xGVVVVwrRp0/Y3Zb377ruv7PXXX+9UVFSU3FK1hauwsDAtNzf3+IgRI3qUl5f7Pvjgg+2dO3eu9bquWMT94bsYVAYsXLiw8w033LA/OTm5SZdr5OTknBgwYMDhuXPnnhXJ62/cuDGiHdvv92Pr1q2p06dPzxo7duyBp59+ek9CAg+mmsvr/SFSLbE/eHqvXyiPrKzKqz7pfl2piaiZPaL9+uauv3fv3oQZM2Z0e++99zKOHj3a5oILLjj68ssv/z03N/dEQ+sVFRUlFxQUpM+bN29n3e3NnDmz28qVK8/w+/0yceLE8oULF+4ZMGBAzvXXX3/g3nvv3QcAI0eOPPD44493WbRo0e6m7AxHjx6VJUuWZCxZsuTMNWvWtPP7/Z839/1t3LgxOSsr63haWpq/pqYm7BrCUfu7FXmoPuH+fpqaVJPwq+vibn9YtWpV+saNG1NGjRpVlZ2dfbIl3l9L7A/mjqi8CKlIX7e0tNQ3aNCg3N27dyfPnz+/+Iknnti5ffv2lMmTJ3dvbN2VK1e2S01N9Q8aNKj61HMVFRUJgwYNyi0sLEybP39+8fTp08sWLVqU+dhjj5319ddfJ8+cObP81LJDhw49XFlZ6VuzZk1qOLWuXr06dfLkyedlZmbm3X777dkZGRk1y5Yt2xbJ+8vPz0/Ly8s7smLFiu2//e1vu3z44Yfp4dQSFi9CKsLXtbw/VFVVJTz88MPdunfv3q9Pnz59b7vttq7vvPNO25Mnw88sL/YHc0dUsegXv/hFt7Zt2/o/+uijL32+wI+0pKQk8bHHHut68ODBNh06dPDXt+66devSe/TocSz4t9+99957TkVFReLq1auLzjnnnBoAB//whz90fvjhh7vOmjWrNC0t7fRHgoEDB1YnJCRg9erV6YMHD64O9RqVlZUJzzzzTMelS5eeuWnTprTc3NzqO++8s2Tq1KmV4Zw3aOz9FRQUpF144YVHMzMza5ctW/bVmDFjzn/vvfe+7N27d4NHD/HK8v4wZsyYqtGjR6//+OOP01asWNHh/fff77Bw4cLM9PT02ksvvbRq+PDhB0eNGlWVlZVVb3J5sT8wqCJUXl6e8Nprr3V68cUXt5/6SwOAnj17HldVHD58uMEdc+/evYkdO3Y8fXzs9/vxyiuvdJo4cWKFs1MCANq1a1d76NChhDvuuKM8eP3ExES0a9eupqysLDHU9pcvX95+0qRJPZOTk3XUqFGVzzzzTPGQIUNCBlpz39+CBQv2nHp+wIABx4qLi4vC3X68sb4/AECbNm0wdOjQo0OHDj0KoLSsrCxh2bJlZzz++ONdVq5cmTFz5kzUdyrAq/3B3Ee/WPPGG2+0V1Vcd911VcHPl5aWJqalpfkzMzNrGrr47fjx45KUlHR6x12/fn3KgQMHfMOHD//W9mpra2Xq1Kn72rdv/52dPDExUY8dOxby7zI5OVlTUlL8x48fl6qqqoT9+/f7/P56/z9p1vsLe2OtgPX9IdiWLVuSHn300bMmTpzY/Z577jmvoqIicciQIVVz5szZFcn7a+x1m4NHVBFau3ZtWrt27Wrbtm37rX+hWbVqVYehQ4ceTEhIaPDit4yMjNp9+/ad/u23bdu2JCDwGypoW+klJSVJF110UchtHDp0KCH4t3Cwa6+99lBJSUnh0qVLz3jhhRfOvO6663p36dLlxPjx4yumTZtW2djheDjvj75hfX8oKipKfvLJJ8/+8MMPO+zcuTO5S5cuJ6644oqDzz777N+vvfbaQ+3atWvwt5hX+wOPqCK0YcOGtEOHDiVUVFSc/htatWpV+ieffNJ+6tSpFUDDF7/16tXr2O7du5NOrXvqcLq8vNwHBA7977777m4AEKrJYUlJie/YsWNtcnJy6m2Zk5qaqtOmTTuwevXqbZs3b94wZsyYypdeeums3NzcfxoyZEjvhQsXdozk/a1duzalb9++fbKysi4YNmzY+QcOHGi1+5X1/eGjjz5K37JlS8rNN99cvm7duo179uzZ8NJLL309ceLEg42FVLjvryX2h1a7Q0XLli1b0jp16nTyxhtv7L58+fL2v/nNb84eM2ZMr7Fjx1Zcf/31VUDDF79ddtllh0tLS5NKSkp8ADBkyJAjycnJetddd3Vbvnx5+7Fjx2b/4x//8GVlZR1ftmxZx61btyYFv/7HH3+cJiIYNmzY4XDqzcnJOTF37tyS3bt3F7788svb09PTa2fMmJEdyfu75ZZbsubMmVOyc+fOot69ex974IEH6u0rFO+s7w/jxo07uGDBgl0jR448CAAFBQUpob4ieX8tsT8wqCKwc+fOxMrKSt9TTz21MykpSSdNmtRz7ty5mVOmTNm3dOnSnUDjF79dc801hzp06FD72muvtQeAzMzM2vnz5+8oKytLmjBhwvmbN29Ofeutt7bdcccdpe+///4ZS5cuzQiuYeXKlR0uvvjiQ5mZmU266tfn82H8+PEHV61a9dWOHTsKm/v+du3a5du1a1fyuHHjDgLA9OnTy998882MUNuLd7GwPyxZsiSjf//+32vsq7nvr6X2B56jikB+fn4qAFx11VWHb7rppoOhlmns4reUlBQdPXp05SuvvNLx1ltv3Q8AU6ZMOTBlypQDwcvl5uZWzpgxozL4uZqaGrz99tsZc+bM2R3J+zj33HNDns8I5/0VFxcnZWZmnj7P1bNnzxNlZWVJoZaNd7GwP8yaNati1qxZFc14e57uD+aOqFIT4cm/IjXndQsKCtK6dOly4qyzzqr3aCaci9/uu+++svz8/HaFhYVNupXl+eefz0hJSfFPnTq1SfeEhSuc96eqEJGWePmA1CRv/lWxGa/L/aHl9gdzR1SR3Mbitg0bNqT16dOnwVYW4Vz81rNnz5Nz584t3rVrV2K/fv2ON7S9YKoqixYtKk5MrPeSmYiE8/6ys7NPlJaWnv6N+dVXXyV17tw5ahd6RnIbi9u4P7Tc/uDZuKz169cX5+XlNesQlGzp379/7l133VU6bty4g7fccku3xMREnTdv3p7G16R41Nz9Yf369Wfm5eVlh/oeg4oilp+fn3rzzTd3P3r0aJsePXocW758+Y5OnTqxxUsr1dz9gUFFROY1FFTmTqYTEdXFoCIi8xhURGQeg4qIzGNQEZF5XgaV3+/3t+AlzUQUK5wsqLd7g2dBJSJl1dXV9d6lTUStR3V1dYqIlNX3fc+Cqqam5sHi4uKkI0eOpPLIiqh18vv9cuTIkdTi4uKkmpqaB+tbzrMLPgFg3bp1w30+3xxVzUQEoVlaWmpiOiwRASLiz8zMrLedcR1+ESmrqal5sH///u/Wu00vgypaRET/4/WQXSca9dDIDgCA5q4fjW3E+voWauB7+Gb9ml+/0qz1AcB3/40RbcNZ/3NVHdjsIkLgv/oRkXkMKiIyj0FFROYxqIjIPAYVEZnHoCIi81wPKhG5WkS2ish2Ebk7xPczROQ1ESkUkTUicoHbNRKRLa4GlYgkAFgAYASAvgAmiEjfOovdA+ALVe0HYDKAuW7WSET2uH1E9X0A21X176p6AsAyACPrLNMXwAcAoKpbAGSLSGd3yyQiS9wOqq4Agi+t3+08F2w9gOsBQES+DyALQDdXqiMik9wOqlA3H9e9h+dRABki8gWAGQAKgO8OBxWRaSLymYh8FvUqicgUV+/1E5HBAB5Q1eHO49kAoKqP1LO8ANgBoJ+qVjWw3di/YZEofsT8vX5rAfQSke4ikgRgPIAVwQuIyBnO9wBgCoCPGgopIop/ro50V9UaEbkNwLsAEgA8r6obReQW5/uLAPQB8KKI1ALYBODn4Wybd93zZxDJ+hZqiKPuCVHnalABgKq+DeDtOs8tCvrzpwB6uV0XEdnFK9OJyDwGFRGZx6AiIvMYVERkHoOKiMxjUBGReQwqIjKPQUVE5jGoiMg8BhURmRc3k5K9roGITov57glERE3m+k3JLYV33fNnEMn6Fmpg94T68YiKiMxjUBGReRbn+nUQkTdEZL2IbBSRn7ldIxHZYnGu360ANqlqHoDLATwR1JqYiFohi3P9FEA7Z7BDWwD7EWIKDRG1Hhbn+s1HoG96CYANAG5XVX/dDXFcFlHrYXGu33AAXwDoAuBCAPNFpP13VlJdrKoDo31hGRHZ43ZQ7QZwbtDjbggcOQX7GYBXNWA7AnP9cl2qj4gMMjfXD8DXAK4EABHpDCAHwN9drZKITLE41+8hAC+IyAYEPirepaoVbtZJRLZYnOtXAuAqt+siIrvYPYGIoo3dE4io9WH3BN51z5+BkRoi7X4QaeeDaGyD3ROIqNViUBGReQwqIjKPQUVE5jGoiMg8BhURmcegIiLzGFREZB6DiojMY1ARkXkMKiIyz+K4rH8XkS+cryIRqRWRjm7XSUR2uNrmxRmX9SWAHyHQlngtgAmquqme5a8FMFNVhzWyXbZ5IbIj5tu8hDMuK9gEAH90pTIiMsvtNi+hxmVdEmpBEUkDcDWA28LZMNuD8GcQyfoWamCblwa22yJbrV8447JOuRbAJ6q6P+SGRKYBmBatwojILreDKpxxWaeMRwMf+1R1MYDFAM9REcU7i+OyICIdAPwAwOsu10dEBlkclwUAowG8p6pH3KyPiGwyNy7LefwCgBfcq4qILOOV6URkHoOKiMxjUBGReQwqIjKPQUVE5jGoiMg8V7sntBRemU5kSsx3TyAiajLXL/hsKbzrnj+DSNaPZg1edT+I5+4JPKIiIvMYVERkHoOKiMxjUBGReQwqIjKPQUVE5pmb6+csc7kz12+jiPzV7RqJyBZXr6Ny5votQNBcPxFZETzXT0TOAPA0gKtV9WsROdvNGonIHotz/SYCeFVVvwYAVd3nco1EZIzbQRVqrl/XOsv0BpAhIv8rIp+LyGTXqiMikyzO9fMBGADgSgCpAD4Vkb+p6pff2hDn+hG1Gq52TxCRwQAeUNXhzuPZAKCqjwQtczeAFFV9wHn8HIB3VLXem4/YPYHIlJjvnhDOXL/XAVwmIj5nrPslADa7XCcRGWJurp+qbhaRdwAUAvADeFZVixrbdjzcdR+r61uowdJ7YPeE6LM61+9xAI+7WRcR2cUr04nIPAYVEZnHoCIi8xhURGQeg4qIzGNQEZF5DCoiMo9BRUTmMaiIyDwGFRGZ52r3hJbC7glEpsR89wQioiZz/abklhIPd93H6voWaojme4jV7gfx3D2BR1REZB6DiojMCyuoRORWpyNnxBqb6+fM9DvozPX7QkTuj8brElHsCveIagKAHSLySxFJDf6GiJwnIr8KZyNBc/1GAOgLYIKI9A2x6P+p6oXO16/DrJGI4lRYQaWqlwJ4CMD9AIpFZLaIXCMiPwbwn87z4Qhnrh8R0beE9a9+IvIIgDsBVAKoAjAVQDYCo662AbgtzNcLNdfvkhDLDRaR9QBKAPxKVTeGqInjsohaiXA/+k0DcK+qnq2qvVS1B4BLAXwKoAbAX8PcTjhz/dYByFLVPADzAPw51IZUdbGqDoz2hWVEZE+4QeUH8K1JMKq62vlIuBZAuBdd7AZwbtDjbggcNQVvt0pVDzt/fhtAooicGeb2iSgOhRtUrwK4X0QyQnxvGYB+YW6n0bl+IpIpIuL8+ftOjZVhbp+I4lC4V6bfjcDHvK9EZDGA9wHsBNARwD0A9oSzkXDm+gG4AcB0EakBUA1gvMbDDYlE1GxhBZWqHhCRSwDcB+BnCJxYVwTOOR0DMDncF2xsrp+qzgcwP9ztEVH8a1b3BBHJAdAdQC0Cd0rvj3ZhTayHR1xEdkS9e0KzbkpW1a0AtkazECKi+rB7AjsH8GdQZxux2v2A3ROIiDzEoCIi8xhURGQeg4qIzGNQEZF5DCoiMo9BRUTmMaiIyDwGFRGZx6AiIvMYVERknutB1di4rKDlLhaRWhG5wc36iMieZrV5afaLBcZlfQngRwi0JV4LYIKqbgqx3CoEel09r6rLG9ku27wQ2RH1Ni9uH1GFOy5rBoA/AdjnZnFEZJPbbV4aHZclIl0BjAYwDMDF4W7YQnuQWG1xEk8/AwstTmJ1/WjWEG1uB1U447J+D+AuVa11ZjyE3hDn+hG1Gm4HVaPjsgAMBLDMCakzAVwjIjWq+ufghVR1MYDFAM9REcU7t4Pq9LgsBCbXjAcwMXgBVe1+6s8i8gKAN+uGFBG1Lq4GVZjjsoiIvsX1numNjcuq8/zNbtRERLbxynQiMo9BRUTmMaiIyDwGFRGZx6AiIvMYVERknqvdE1oKr0wnMiXmuycQETWZ6xd8thR2DoiPn0E8dA6I1fWjWUO08YiKiMxjUBGReQwqIjKPQUVE5jGoiMg8BhURmWdurp+IjBSRQhH5QkQ+E5FL3a6RiGxx9ToqZ17fAgTN9RORFXXm+n0AYIWqqoj0A/A/AHLdrJOIbDE3109VD+s39/Wk47tTaoiolXE7qELN9etadyERGS0iWwC8BeBfXaqNiIxyO6jCmesHVX1NVXMBjALwUMgNiUxzzmF9Ft0SicgaV7sniMhgAA+o6nDn8WwAUNVHGlhnB4CLVbWigWX48ZDIjpjvnnB6rp+IJCEw129F8AIicr4400dFpD+AJACVLtdJRIZYnOs3BsBkETkJoBrAOA3jsC8eOgfE6vrRrCEeOgfE6vrRrCHazM31U9XHADzmdl1EZBevTCci8xhURGQeg4qIzGNQEZF5DCoiMo9BRUTmMaiIyDwGFRGZx6AiIvMYVERknqvdE1oKuycQmRLz3ROIiJrM9ZuSW0o8dA6I1fWDt8HOAbG7fjRriDYeURGReQwqIjLP4ly/m5y5foUislpE8tyukYhscTWogub6jQDQF8AEEelbZ7EdAH6gqv0QGOyw2M0aicgei3P9VqvqAefh3wB0c7lGIjLG5Fy/ID8HsDLUNzgui6j1cPvyhLDm+gGAiFyBQFBdGur7qroYzsdCXvBJFN/cDqrdAM4NetwNQEndhUSkH4BnAYxQVY7KImrlLM71Ow/AqwAmqeqXLtdHRAZZnOt3P4BOAJ525pDWRPu+ISKKLRbn+k0BMMXtuojILnZPIKJoY/cEImp94qZ7QqR3e7N7go277uOhc0Csrh/NGqKNR1REZB6DiojMY1ARkXkMKiIyj0FFROYxqIjIPAYVEZnHoCIi8xhURGQeg4qIzGNQEZF5Fsdl5YrIpyJyXER+5XZ9RGSPq21enHFZXwL4EQJtidcCmKCqm4KWORtAFoBRAA6o6u/C2C7bvBDZEfNtXsIZl7VPVdcCOOlybURklNttXkKNy7okGhtmmxe2eYlkfQs1eL1+NGuINrPjshrdkMg0ANMiK4eIYoHJcVnh4Fw/otbD3LgsIqK6zI3LEpFMAJ8BaA/ALyK/BNBXVavcrJWI7LA4LqsMgY+EREQAeGU6EcUABhURmcegIiLzGFREZB6DiojMY1ARkXmudk9oKbwynciUqHdPcP06qhZSAWCnR699pvP6VlmvD2CN0WCpvqxobzAujqi8JCKfRfu3RzRZrw9gjdFgvb5I8RwVEZnHoCIi8xhUkVvsdQGNsF4fwBqjwXp9EeE5KiIyj0dURGQegypMInKuiPxFRDaLyEYRud15vqOIrBKRbc5/MzyuM0FECkTkTaP1nSEiy0Vki/OzHGywxpnO33GRiPxRRFK8rlFEnheRfSJSFPRcvTWJyGxnJN1WERnuZq0tgUEVvhoAd6hqHwCDANwqIn0B3A3gA1XtBeAD57GXbgewOeixtfrmAnhHVXMB5CFQq5kaRaQrgF8AGKiqFyDQ4HG8gRpfAHB1nedC1uTsl+MBfM9Z52lnVF3sUlV+NeMLwOsIzCfcCuAc57lzAGz1sKZuCOywwwC86Txnqb72AHbAOTca9LylGk9NSuqIwAXRbwK4ykKNALIBFDX2cwMwG8DsoOXeBTDYq59pNL54RNUMIpIN4CIA+QA6q2opADj/PdvD0n4P4E4A/qDnLNXXA0A5gP9yPp4+KyLplmpU1T0AfgfgawClAA6q6nuWagxSX02hxtJ1dbm2qGJQNZGItAXwJwC/VEN93EXkXwDsU9XPva6lAT4A/QEsVNWLAByB9x9Fv8U5zzMSQHcAXQCki8hPvK2qyaI2ls4KBlUTiEgiAiH136r6qvP0XhE5x/n+OQD2eVTePwO4TkSKEZhAPUxEXjJUHxD4zb5bVfOdx8sRCC5LNf4QwA5VLVfVkwBeBTDEWI2n1FdT1MbSWcGgCpOICIDnAGxW1SeDvrUCwE+dP/8UgXNXrlPV2araTVWzETiR+qGq/sRKfcDpwR27RCTHeepKAJtgqEYEPvINEpE05+/8SgRO+Fuq8ZT6aloBYLyIJItIdwC9AKzxoL7o8fokWax8AbgUgcPnQgBfOF/XAOiEwAnsbc5/Oxqo9XJ8czLdVH0ALkRgHFohgD8DyDBY44MAtgAoArAUQLLXNQL4IwLnzE4icMT084ZqAnAvgK8QOOE+wut9MtIvXplORObxox8RmcegIiLzGFREZB6DiojMY1ARkXkMKiIyj0FFnhKR8SJy3LnqnygkBhV5LQ/AJg3crkIUEoOKvHYhgAKviyDbGFTktTwEbkcCAIjIMBGpFJHfx3yzN4oaBhV5RkTOQqDh2xfO49sQaFQ3W1V/qaq1HpZHhsTLSHeKTRcicKP3JhFZDGAMgB+r6l88rYrMYVCRl/IA7Eeg59PZAC5R1e3elkQWMajIS3kIdKO8DMAPGVJUH7Z5Ic+IyAYEunz2RKCH1gBVLfe0KDKJJ9PJEyKSBCAXgQZ6/4bA0IdXRIRH+fQdDCryyvcQOPWwQVWrAVzvPPeUp1WRSQwq8koeAlNovgIAVd0JYByA6SJys4d1kUE8R0VE5vGIiojMY1ARkXkMKiIyj0FFROYxqIjIPAYVEZnHoCIi8xhURGQeg4qIzPt/JvDu6AWMS8kAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "PP_SIM_KS = list(range(5, 120 + 1, 120 // 15))\n",
    "PP_SIM_ALPHAS = list(reversed([x/1000 for x in range(50, 950 + 1, 900 // 15)]))\n",
    "PP_SIM_HPS = np.meshgrid(PP_SIM_KS, PP_SIM_ALPHAS)\n",
    "\n",
    "# compute\n",
    "pp_sim_improvements = np.array([\n",
    "    pp_sim_radius(k, 1) > pp_sim_radius(k, alpha) for\n",
    "    (k, alpha) in zip(PP_SIM_HPS[0].flatten(), PP_SIM_HPS[1].flatten())\n",
    "]).reshape(PP_SIM_HPS[0].shape)\n",
    "pp_sim_cs = [C['PP_SIM'](k) for k in PP_SIM_KS]\n",
    "\n",
    "# and plot\n",
    "edge_colors = [ ]\n",
    "cmap = matplotlib.colors.ListedColormap([COLORS[1], COLORS[0]])\n",
    "plt.figure(figsize=(4, 4))\n",
    "plt.pcolormesh(\n",
    "    PP_SIM_HPS[0], PP_SIM_HPS[1], pp_sim_improvements, \n",
    "    edgecolors='k',\n",
    "    cmap=cmap,\n",
    ")\n",
    "plt.xlabel(r'$k$', fontsize=15)\n",
    "plt.ylabel(r'$\\alpha$', fontsize=15)\n",
    "\n",
    "# where some data has already been plotted to ax\n",
    "handles, labels = plt.gca().get_legend_handles_labels()\n",
    "\n",
    "# manually define a new patch \n",
    "patch_pos = mpatches.Patch(\n",
    "    color=COLORS[0], \n",
    "    label=r'$\\rho_k(\\alpha)<\\rho_0^k$'\n",
    ")\n",
    "patch_neg = mpatches.Patch(\n",
    "    color=COLORS[1],\n",
    "    label=r'$\\rho_k(\\alpha)\\geq\\rho_0^k$'\n",
    ")\n",
    "line_suff = mlines.Line2D([0], [0], color=COLORS[3], label='Thm. 8', linewidth=3)\n",
    "\n",
    "# handles is a list, so append manual patch\n",
    "handles.append(patch_pos) \n",
    "handles.append(patch_neg) \n",
    "#handles.append(line_suff)\n",
    "\n",
    "# plot the legend\n",
    "plt.legend(\n",
    "    handles=handles, loc='upper center', \n",
    "    bbox_to_anchor=(-0.01, 1.20, 1.0, 0), \n",
    "    ncol=2, fontsize=15, columnspacing=1,\n",
    ")\n",
    "plt.tight_layout(pad=0.1)\n",
    "plt.savefig('./bilinear_sweep_pp.png', format='png', dpi=200)\n",
    "print(f'=> beta1: {PP_SIM_BETA1}, beta2: {PP_SIM_BETA2}')\n",
    "print(f'=> c: {pp_sim_cs}')\n",
    "print(f'=> a: {PP_SIM_ALPHAS}')\n",
    "print(f'=> k: {PP_SIM_KS}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# EG_SIM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-7-b5ec719c936b>:16: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3.  Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading'].  This will become an error two minor releases later.\n",
      "  plt.pcolormesh(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "beta1: 20.31, beta2: 76.16\n",
      "=> c: [0.9999999999999999, 1.0000000000000004, 0.9386372877575877, 0.88189405962177, 0.8538707136828575, 0.8755727534510303, 0.90238232625363, 0.923621131781686, 0.9404464925478963, 0.95377865415444, 0.9643498736793074, 0.9727451081119008, 0.9794374463024417, 0.9848302882829575, 0.9896994732761079, 0.9979913799116246, 1.0000000000000002]\n",
      "=> a: [0.95, 0.94, 0.93, 0.92, 0.91, 0.9, 0.89, 0.88, 0.87, 0.86, 0.85, 0.84, 0.83, 0.82, 0.81, 0.8]\n",
      "=> k: [15, 19, 23, 27, 31, 35, 39, 43, 47, 51, 55, 59, 63, 67, 71, 75, 79]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS0AAAErCAYAAACLutXzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAguUlEQVR4nO3dfZAV9Z3v8feHGQYGRERjGBHDIEEe1isKrkFiKEN2RbIbNZoouMHVWuRigSHqLRX1LslNdtesaxJKCIREb4wmS62sXtFgEDW5xpBCkOfHFQVkHFBAAyhPDvO9f5zGezIOMMCcPqeHz6uqy9N9+nf6e2baD7/u6f61IgIzs6xoVewCzMyOhkPLzDLFoWVmmeLQMrNMcWiZWaY4tMwsU8qLXUC+RYsWDS0vL58YEVU4UM1ONPWSttTV1X2nf//+cw61kkrlOq1FixYNbdOmzeTq6ur9lZWVe1u1alUahZlZKurr67Vnz562GzZsqNi3b9+4QwVXyfRmysvLJ1ZXV+9v3779HgeW2YmnVatW0b59+z3V1dX7y8vLJx5yvTSLOpyIqKqsrNxb7DrMrLgqKyv3JqeIGlUyoQW0cg/LzJIcOGQ2lVJomZkdkUPLzDLFoWUF8eqrr1ZeeOGFvYpdhxVfc+8LDi0riPnz57c777zzdhe7Diu+5t4XHFolqL6+nl69evV96KGHTjuadiNHjvzMtdde261QdR2NxYsXtxswYMCHABMnTux86aWXfvadd94pK3ZdWeN94ZMcWiXo4Ycf7rRz586y0aNHv3c07e67774tTz/99GkrVqxoU6jammrZsmXtevfuvW/YsGFnb926tfzFF19c17lz5wPFritrvC98kkOrBE2dOrXz1772tffatGlzVJeA9OrVa/+AAQM+mDRp0unHs/2VK1ce145eX1/P2rVrK2+55ZZu11577fs//vGP3y4rcyfrWBR7XzhehdgXSurew8b8y3M7++35KP06K1tTN2HYyUuPtf0777xTduutt3Z9/vnnO+3evbvVueeeu/tXv/rVm717995/uHYrVqxos3jx4vYPPfTQxoafd9ttt3V97rnnTqmvr9f111+/derUqW8PGDCg19VXX/3+vffe+y7AlVde+f4DDzzQZdq0aTVHs3Ps3r1bjz76aKdHH330U6+++mqH+vr61471+61cubJNt27d9rVr166+rq6uyTUcyYF/m9WPPfvT32crK+rK/scVLW5fmDt3bvuVK1e2veqqq3ZWV1d/VIjvV4h9oeR7WsUIrOPd7ubNm8sHDhzYu6amps3kyZM3PPjggxvXrVvX9oYbbuh+pLbPPfdch8rKyvqBAwfuObhs27ZtZQMHDuy9bNmydpMnT95wyy23bJk2bVrV97///dPfeuutNrfddtvWg+sOHjz4g+3bt5e/+uqrlU2pdd68eZU33HDDZ6qqqvqNHz++ulOnTnUzZsx4/Xi+3/z589v169fvw1mzZq3713/91y4vvfRS+6bUckTFCKzj3G4p7ws7d+4s+973vte1e/fu5/Xp06fvuHHjzvzNb35z0kcfNT2/irEvlHxPK4u++c1vdj3ppJPqX3755f8qL8/9iGtra1t///vfP3PHjh2tOnbsWH+otosWLWp/9tln783/l/Hee+89Y9u2ba3nzZu34owzzqgDdvzkJz/p/L3vfe/M22+/fXO7du0+PnS48MIL95SVlTFv3rz2F1988Z7GtrF9+/ayn/70p6c+9thjn1q1alW73r1777nzzjtrb7755u1NOddwpO+3ePHidueff/7uqqqqAzNmzHjjmmuu+ezzzz//X+ecc85hexYtUSnvC9dcc83Or371q0tfeeWVdrNmzer4wgsvdJw6dWpV+/btD1xyySU7hw4duuOqq67a2a1bt0OmWDH2BYdWM9u6dWvZU089ddovfvGLdQd/iQA9evTYFxF88MEHh91R33nnndannnrqx/3o+vp6nnjiidOuv/76bclOCkCHDh0O7Nq1q+yOO+7Ymt++devWdOjQoW7Lli2tG/v8mTNnnjxy5Mgebdq0iauuumr7T3/60w2DBg1qNNyO9ftNmTLl7YPLBwwYsHfDhg0rmvr5LUmp7wsArVq1YvDgwbsHDx68G9i8ZcuWshkzZpzywAMPdHnuuec63XbbbRzqVEGx9oWSPzzMmmeeeebkiOCKK67Ymb988+bNrdu1a1dfVVVVd7iL7fbt26eKioqPd+SlS5e2ff/998uHDh36Z5934MAB3Xzzze+efPLJn9jpW7duHXv37m30d9umTZto27Zt/b59+7Rz586y9957r7y+/pD/3xzT92vyh7Vwpb4v5FuzZk3F/ffff/r111/f/Z577vnMtm3bWg8aNGjnxIkTNx3P9zvSdo+Fe1rNbMGCBe06dOhw4KSTTvqzv/bMnTu34+DBg3eUlZUd9mK7Tp06HXj33Xc//pfx9ddfr4Dcv155n9W+tra24oILLmj0M3bt2lWW/y90vq985Su7amtrlz322GOn/PznP//UFVdccU6XLl32Dx8+fNvo0aO3H6nb3pTvZzmlvi+sWLGizQ9+8INPv/TSSx03btzYpkuXLvu/+MUv7vjZz3725le+8pVdHTp0OOy/ZsXaF9zTambLly9vt2vXrrJt27Z9/BubO3du+z/84Q8n33zzzdvg8Bfb9ezZc29NTU3FwbYHu91bt24th9whwt13390VoLEBHGtra8v37t3bqlevXocc5qeysjJGjx79/rx5815fvXr18muuuWb7448/fnrv3r3/26BBg86ZOnXqqcfz/RYsWNC2b9++fbp163bukCFDPvv++++fkPtZqe8LL7/8cvs1a9a0vfHGG7cuWrRo5dtvv7388ccff+v666/fcaTAaur3K8S+cELuTIW0Zs2adqeddtpHX//617vPnDnz5H/6p3/69DXXXNPz2muv3Xb11VfvhMNfbPeFL3zhg82bN1fU1taWAwwaNOjDNm3axF133dV15syZJ1977bXVf/rTn8q7deu2b8aMGaeuXbu2In/7r7zySjtJDBky5IOm1NurV6/9kyZNqq2pqVn2q1/9al379u0P3HrrrdXH8/3GjBnTbeLEibUbN25ccc455+z99re/fcixkVqyUt8Xrrvuuh1TpkzZdOWVV+4AWLx4cdvGpuP5foXYFxxazWjjxo2tt2/fXv7DH/5wY0VFRYwcObLHpEmTqkaNGvXuY489thGOfLHdl7/85V0dO3Y88NRTT50MUFVVdWDy5Mnrt2zZUjFixIjPrl69uvLXv/7163fcccfmF1544ZTHHnusU34Nzz33XMe//Mu/3FVVVXVUVxyXl5czfPjwHXPnzn1j/fr1y471+23atKl806ZNba677rodALfccsvWZ599tlNjn9eSZWFfePTRRzv179//L440Hev3K9S+4HNazWj+/PmVAJdddtkHf/d3f7ejsXWOdLFd27Zt46tf/er2J5544tSxY8e+BzBq1Kj3R40a9X7+er17995+6623bs9fVldXx+zZsztNnDix5ni+x1lnndXoOZCmfL8NGzZUVFVVfXxerEePHvu3bNlS0di6LVkW9oXbb7992+23377tGL5eUfeFku9pVbamKH+NOpbtLl68uF2XLl32n3766Yfs5TTlYrv77rtvy/z58zssW7bsqG6neeSRRzq1bdu2/uabbz6q+9SaqinfLyKQVIjNQ2VFcf4yeQzb9b5QuH2h5Htax3MrTdqWL1/erk+fPocdgqMpF9v16NHjo0mTJm3YtGlT6/POO2/f4T4vX0Ro2rRpG1q3PuRlOcelKd+vurp6/+bNmz/+1/SNN96o6Ny5c7NcVHo8t9KkzftC4faFknmE2NKlSzf069fvmLqqVlr69+/f+6677tp83XXX7RgzZkzX1q1bx0MPPfT2kVtaS3Os+8LSpUs/1a9fv+rG3nNoWbObP39+5Y033th99+7drc4+++y9M2fOXH/aaad5WJoT0LHuCw4tM8uUw4VWyZ+INzPL59Ays0xxaJlZpji0zCxTHFpmlimlFFr19fX1BbqU2syyIsmBQ44yUTKhJWnLnj17DnlHuZmdGPbs2dNW0pZDvV8yoVVXV/edDRs2VHz44YeV7nGZnXjq6+v14YcfVm7YsKGirq7uO4dar2QuLgVYtGjR0PLy8okRUUUzBermzZtL4im7ZicKSfVVVVWHHKb5MOolbamrq/tO//795xzy80sptApBUtT9ryeavH75P34dgKNpk9/ufz7d6Cgdh/TdKzu22HbFqvFYf3dptEtzW0Vs91pEXHhUDY9CyRwempk1hUPLzDLFoWVmmeLQMrNMcWiZWaY4tMwsUxxaZpYpDi0zyxSHlpllSuqhJelySWslrZN0dyPvd5L0lKRlkl6VdG6D98skLZb0bHpVm1mpSDW0JJUBU4BhQF9ghKS+DVa7B1gSEecBNwCTGrw/Hlhd6FrNrERFRGoTcDEwJ29+AjChwTq/Bi7Jm38D6Jy87gq8CAwBnm3iNsOTJ0+pTgsLmSNpHx6eCeTf/V2TLMu3FLgaQNJFQDdyYQXwI+BODjNAWNJutKSFkhY2Q81mVkLKU95eY+NkRYP5+4FJkpYAy4HFQJ2kvwXejYjXJF16uI1ExHRgOniUh2K2a8mjNRxruyzU2BztCint0KoBzsqb7wrU5q8QETuBmwAkCVifTMOBKyR9GWgLnCzp8Yj4RhqFm1lpSPvwcAHQU1J3SRXkgmhW/gqSTkneAxgFvBwROyNiQkR0jYjqpN1LDiyzE0+qPa2IqJM0DpgDlAGPRMRKSWOS96cBfYBfSDoArAL+Ic0azay0pX14SETMBmY3WDYt7/UfgZ5H+IzfAb8rQHlmVuJ8RbyZZYpDy8wyxaFlZpni0DKzTHFomVmmOLTMLFMcWmaWKSfEE6aLXYPZCcZPmDYzOyj1K+KLIQt34ZfyaA3H2s6jPBR3W8VsV0juaZlZpji0zCxTHFpmlikOLTPLFIeWmWWKQ8vMMsWhZWaZ4tAys0xxaJlZpqQeWpIul7RW0jpJdzfyfidJT0laJulVSecmy8+S9FtJqyWtlDQ+7drNrPhSDS1JZcAUYBjQFxghqW+D1e4BlkTEecANwKRkeR1wR0T0AQYCYxtpa2YtXUSkNgEXA3Py5icAExqs82vgkrz5N4DOjXzW08BfN2Gb4cmTp1SnhYXMkbQPD88ENuXN1yTL8i0FrgaQdBHQjdyTqD8mqRq4AJjf2EYkjZa0UNLC5inbzEpF2qM8qJFl0WD+fmCSpCXAcmAxuUPD3AdIJwH/CXwrInY2tpGImA5MT9aPLNyFX8qjNRxru5Y8WsOxtstCjc3RrpDSDq0a4Ky8+a5Abf4KSRDdBCBJwPpkQlJrcoH1y4h4Mo2Czay0pH14uADoKam7pApgODArfwVJpyTvAYwCXo6InUmAPQysjogfpFq1mZWMVHtaEVEnaRwwBygDHomIlZLGJO9PA/oAv5B0AFgF/EPS/PPASGB5cugIcE9EzE7zO5hZcaU+cmkSMrMbLJuW9/qPQM9G2r1C4+fEzOwE4ivizSxTHFpmlikOLTPLFIeWmWWKQ8vMMsWhZWaZ4tAys0xxaJlZpigZvqXFktSyv6BZ6XktIi4s1Ie7p2VmmZL6bTzFkIWhQ0p5iJljbeehaYq7rWK2KyT3tMwsUxxaZpYpDi0zyxSHlpllikPLzDLFoWVmmeLQMrNMcWiZWaakHlqSLpe0VtI6SXc38n4nSU9JWibpVUnnNrWtmbV8qYaWpDJgCjAM6AuMkNS3wWr3AEsi4jzgBmDSUbQ1sxYu7Z7WRcC6iHgzIvYDM4ArG6zTF3gRICLWANWSOjexrZm1cKmO8iDpa8DlETEqmR8JfC4ixuWt889A24i4XdJFwDzgc0D3I7U9xDY9yoNZulrUKA+NPbewYajcD3RKHsh6K7AYqGti29xGpNGSFkpaeBy1mlkJSnuUhxrgrLz5rkBt/goRsRO4CUCSgPXJ1O5IbfM+YzowPfmMyMJd+KU8WkN+uyz8LEu5XRZqbI52hZR2T2sB0FNSd0kVwHBgVv4Kkk5J3gMYBbycBNkR25pZy5dqTysi6iSNA+YAZcAjEbFS0pjk/WlAH+AXkg4Aq4B/OFzbNOs3s+JLfRDAiJgNzG6wbFre6z8CPZva1sxOLL4i3swyxaFlZpni0DKzTHFomVmmOLTMLFMcWmaWKQ4tM8sUh5aZZUqqozwUg0d5MEtdixrlwczsuKR+G08xZOEufI/ycGK0y0KNzdGukNzTMrNMcWiZWaY4tMwsUxxaZpYpDi0zyxSHlpllikPLzDLFoWVmmeLQMrNMST20JF0uaa2kdZLubuT9jpKekbRU0kpJN+W9d1uybIWkf5fUNt3qzazYUg0tSWXAFGAY0BcYIalvg9XGAqsioh9wKfCgpApJZwLfBC6MiHPJPUZseGrFm1lpiIjUJuBiYE7e/ARgQoN1JgA/BgR0B9aRC9czgU3AqeTumXwWuKwJ2wxPnjylOi0sZI6kfXh4MHgOqkmW5ZtM7oGttcByYHxE1EfE28C/AW8Bm4EdEfF8YxuRNFrSQkkLm/sLmFlxpT3KgxpZFg3mhwJLgCFAD2CupN+TOxy8klzv60/AE5K+ERGPf+IDI6YD0yE3nlZLvgvfozxkq10WamyOdoWUdk+rBjgrb74ruR5VvpuAJyNnHbAe6A38FbA+IrZGxEfAk8CgFGo2sxKSdmgtAHpK6i6pgtyJ9FkN1nkL+BKApM5AL+DNZPlASe0kKVlndWqVm1lJaFJoSRqbhMxxiYg6YBwwh1zg/EdErJQ0RtKYZLXvAoMkLQdeBO6KiG0RMR+YCSwid66rFckhoJmdOJp6TmsEcI+kB4CfRMSeg29I+gxwbUT8W1M+KCJmA7MbLJuW97oWuOwQbScCE5tYs5m1QE0KrYi4JOkJ/TMwQdKPgKXkTqyPAK4g95c9M7OCalJoSfoX4E5gO7ATuBmoJveXv9fJHfKZmRVcU0/EjwbujYhPR0TPiDgbuAT4I1AH/N9CFWhmlq+poVUPrMhfEBHzIuIScn8RPLqLOczMjlFTQ+tJ4B8ldWrkvRnAec1XkpnZoTX1r4d3kzsUfEPSdOAFYCO5+wDvAd4uTHlmZn+uqX89fF/S54D7yF2xfie5k/AC9gI3FKxCM7M8SkZCOLpGUi9y9wAeAF6LiPeau7DmIunov6CZHY/XIuLCQn34Md0wHRFrgbXNXIuZ2RGlPcpDUbTku/DTGK0hf3st+Wfpn0nztSskjxFvZpni0DKzTHFomVmmOLTMLFMcWmaWKQ4tM8sUh5aZZYpDy8wyxaFlZpmSemhJulzSWknrJN3dyPsdJT0jaamklZJuynvvFEkzJa2RtFrSxelWb2bFlmpoSSoDpgDDgL7ACEl9G6w2FlgVEf2AS4EH854ENAn4TUT0BvrhR4iZnXDS7mldBKyLiDcjYj+5AQSvbLBOAB2SZxueBLwH1Ek6GRgMPAwQEfsj4k+pVW5mJeGYhqY55o1JXwMuj4hRyfxI4HMRMS5vnQ7kHuDaG+gAXBcRv5Z0PrnnHK4i18t6DRgfER8eYZsemsYsXQUdmibtnpYaWdYwVIYCS4AuwPnA5KSXVQ70B6ZGxAXAh+RGVP3kRqTRkhZKWthMdZtZiUh7aJoa4Ky8+a5AbYN1bgLuj1wXcJ2k9eR6XW8BNcmTpiH3tOlGQysippM8fVpStOShQzw0TbbaZaHG5mhXSGn3tBYAPSV1T06uDyd3KJjvLeBLAJI6A72ANyNiC7ApGTWVZJ1V6ZRtZqUi1Z5WRNRJGgfMAcqARyJiZfL0aiJiGvBd4OeSlpM7nLwrIrYlH3Er8Msk8N4k1yszsxNI6iOXRsRsYHaDZdPyXtcClx2i7RKgYCf4zKz0+Yp4M8sUh5aZZYpDy8wyxaFlZpni0DKzTHFomVmmOLTMLFMcWmaWKamO8lAMHuXBLHUtapQHM7PjkvptPMXgu/CL0y4LNabdLgs1Nke7QnJPy8wyxaFlZpni0DKzTHFomVmmOLTMLFMcWmaWKQ4tM8sUh5aZZYpDy8wyJfXQknS5pLWS1kn6xHMLJXWU9IykpZJWSrqpwftlkhZLeja9qs2sVKQaWpLKgCnAMKAvMEJS3warjQVWRUQ/4FLgweSRYQeNB1anUK6ZlaKISG0CLgbm5M1PACY0WGcC8GNyzzzsDqwDWiXvdQVeBIYAzzZxm+HJk6dUp4WFzJG0Dw/PBDblzdcky/JNBvoAtcByYHxE1Cfv/Qi4E6jnMCSNlrRQ0sLmKNrMSkfaozyokWXRYH4osIRcb6oHMFfS74HBwLsR8ZqkSw+3kYiYDkyH3Hhavgu/OO2yUGPa7bJQY3O0K6S0e1o1wFl5813J9ajy3QQ8GTnrgPVAb+DzwBWSNgAzgCGSHi98yWZWStIOrQVAT0ndk5Prw4FZDdZ5C/gSgKTOQC/gzYiYEBFdI6I6afdSRHwjvdLNrBSkengYEXWSxgFzgDLgkYhYKWlM8v404LvAzyUtJ3c4eVdEbEuzTjMrXamPXBoRs4HZDZZNy3tdC1x2hM/4HfC7ApRnZiXOV8SbWaY4tMwsUxxaZpYpDi0zyxSHlpllikPLzDLFoWVmmaJkJIQWS1LL/oJmpee1iLiwUB/unpaZZUrqV8QXg+/CL067LNSYdrss1Ngc7QrJPS0zyxSHlpllikPLzDLFoWVmmeLQMrNMcWiZWaY4tMwsUxxaZpYpDi0zy5TUQ0vS5ZLWSlon6e5G3u8o6RlJSyWtlHRTsvwsSb+VtDpZPj7t2s2s+FINLUllwBRgGNAXGCGpb4PVxgKrIqIfcCnwYPK4sTrgjojoAwwExjbS1sxauFRHeZB0MfDtiBiazE8AiIh/yVtnArkHuo4FqoG5wDkRUd/gs54GJkfE3CNs06M8mKWrRY3ycCawKW++JlmWbzLQh9yTp5cD4xsJrGrgAmB+YxuRNFrSQkkLm6luMysRaY/yoEaWNewJDQWWAEOAHsBcSb+PiJ0Akk4C/hP41sFln/jAiOnA9GT98F34xWmXhRrTbpeFGpujXSGl3dOqIXfod1BXcj2qfDcBT0bOOmA90BtAUmtygfXLiHgyhXrNrMSkHVoLgJ6Suicn14cDsxqs8xbwJQBJnYFewJuSBDwMrI6IH6RYs5mVkFRDKyLqgHHAHGA18B8RsVLSGEljktW+CwyStBx4EbgrIrYBnwdGAkMkLUmmL6dZv5kVX+ojl0bEbGB2g2XT8l7XApc10u4VGj8nZmYnEF8Rb2aZ4tAys0xxaJlZpji0zCxTHFpmlikOLTPLFIeWmWWKQ8vMMiXVoWmKwUPTmKWuRQ1NY2Z2XFK/jacYPHRIcdploca022WhxuZoV0juaZlZpji0zCxTHFpmlikOLTPLFIeWmWWKQ8vMMsWhZWaZ4tAys0xJPbQkXS5praR1ku5u5P2Okp6RtFTSSkk3NbWtmbV8qYaWpDJgCjAM6AuMkNS3wWpjgVUR0Q+4FHhQUkUT25pZC5d2T+siYF1EvBkR+4EZwJUN1gmgQ/Kcw5OA94C6JrY1s5YuIlKbgK8BP8ubHwlMbrBOB+C3wGbgA+Bvmtr2ENsMT548pTotLGSOpH3DdGPPLYwG80OBJcAQoAcwV9Lvm9g2txFpNDA6mf0AWHssxeb5FLDtOD+jOZRKHeBaGlMqdUBxa+lWyA9PO7RqgLPy5rsCtQ3WuQm4P3LdpHWS1gO9m9gWgIiYDkxvrqIlLSzk+EBZqwNcSynXAaVVS3NL+5zWAqCnpO6SKoDhwKwG67wFfAlAUmegF/BmE9uaWQuXak8rIuokjQPmAGXAIxGxUtKY5P1pwHeBn0taTu6Q8K6I2AbQWNs06zez4kt9EMCImA3MbrBsWt7rWuCyprZNSbMdah6nUqkDXEtjSqUOKK1amlWLHyPezFoW38ZjZpni0Moj6SxJv5W0OrmFaHyy/FRJcyW9nvy3Uwq1tJX0at7tTN8pVi3JdsskLZb0bJHr2CBpuaQlkhYWqxZJp0iaKWlNsr9cXKQ6eiU/i4PTTknfKtbvJw0OrT9XB9wREX2AgcDY5Fahu4EXI6In8GIyX2j7gCHJ7UznA5dLGlikWgDGA6vz5otVB8AXI+L8vD/pF6OWScBvIqI30I/czyb1OiJibfKzOB8YAOwGnipGLalJ84r4rE3A08Bfk7s49Yxk2RnA2pTraAcsAj5XjFrIXRP3IrkLfp9NlhXlZwJsAD7VYFmqtQAnA+tJzgkXq45G6roM+EMp1FLIyT2tQ5BUDVwAzAc6R8RmgOS/n06phjJJS4B3gbkRUaxafgTcCdTnLSvKz4TcXRDPS3otufOhGLWcDWwF/ndyyPwzSe2LUEdDw4F/T14Xu5aCcWg1QtJJwH8C34qIncWqIyIORK7b3xW4SNK5adcg6W+BdyPitbS3fQifj4j+5Eb7GCtpcBFqKAf6A1Mj4gLgQ4p8+JVccH0FcHQPKswgh1YDklqTC6xfRsSTyeJ3JJ2RvH8GuZ5PaiLiT8DvgMuLUMvngSskbSA3ssYQSY8XoQ7g4+v4iIh3yZ27uagItdQANUnPF2AmuRAr5n4yDFgUEe8k80XdZwvJoZUnGQ7nYWB1RPwg761ZwN8nr/+e3LmuQtdyuqRTkteVwF8Ba9KuJSImRETXiKgmd/jxUkR8I+06ACS1l9Th4Gty53BWpF1LRGwBNknqlSz6ErAq7ToaGMH/PzSkyLUUVrFPqpXSBFxC7pzJMnIjTSwBvgycRu5E9OvJf09NoZbzgMVJLSuAf0yWp15LXk2X8v9PxBfjZ3I2sDSZVgL3FrGW84GFye/n/wCdivW7IfeHmu1Ax7xlRdtPCj35ingzyxQfHppZpji0zCxTHFpmlikOLTPLFIeWmWWKQ8vMMsWhZSVF0nBJ+5I7E8w+waFlpaYfuSeMf1TsQqw0ObSs1JxP7k4As0Y5tKzU9CN3+xQAkoZI2i7pR5LKileWlQqHlpUMSaeTG7BuSTI/DngWmBAR34qIA0Usz0pE6o8QMzuM88ndsL5K0nTgGuBvIuK3Ra3KSopDy0pJP+A94ElyI21+LiLWFbckKzUOLSsl/cg9VfwLwF85sKwxHprGSoak5eRGAe1BbuyuARGxtahFWcnxiXgrCckY573JDar338k9OOIJST4asD/j0LJS8RfkTlcsj4g9wNXJsh8WtSorOQ4tKxX9yD3V5g2AiNgIXAfcIunGItZlJcbntMwsU9zTMrNMcWiZWaY4tMwsUxxaZpYpDi0zyxSHlpllikPLzDLFoWVmmeLQMrNM+X/tqz8lKS8AygAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "EG_SIM_KS = list(range(15, 80 + 1, 65 // 15))\n",
    "EG_SIM_ALPHAS = list(reversed([x/1000 for x in range(800, 950 + 1, 150 // 15)]))\n",
    "EG_SIM_HPS = np.meshgrid(EG_SIM_KS, EG_SIM_ALPHAS)\n",
    "\n",
    "# compute\n",
    "eg_sim_improvements = np.array([\n",
    "    eg_sim_radius(k, 1) > eg_sim_radius(k, alpha) for\n",
    "    (k, alpha) in zip(EG_SIM_HPS[0].flatten(), EG_SIM_HPS[1].flatten())\n",
    "]).reshape(EG_SIM_HPS[0].shape)\n",
    "eg_sim_cs = [C['EG_SIM'](k) for k in EG_SIM_KS]\n",
    "\n",
    "# plot\n",
    "edge_colors = [ ]\n",
    "cmap = matplotlib.colors.ListedColormap([COLORS[1], COLORS[0]])\n",
    "plt.figure(figsize=(4, 4))\n",
    "plt.pcolormesh(\n",
    "    EG_SIM_HPS[0], EG_SIM_HPS[1], eg_sim_improvements, \n",
    "    edgecolors='k',\n",
    "    cmap=cmap,\n",
    ")\n",
    "plt.xlabel(r'$k$', fontsize=15)\n",
    "plt.ylabel(r'$\\alpha$', fontsize=15)\n",
    "\n",
    "# where some data has already been plotted to ax\n",
    "handles, labels = plt.gca().get_legend_handles_labels()\n",
    "\n",
    "# manually define a new patch \n",
    "patch_pos = mpatches.Patch(\n",
    "    color=COLORS[0], \n",
    "    label=r'$\\rho_k(\\alpha)<\\rho_0^k$'\n",
    ")\n",
    "patch_neg = mpatches.Patch(\n",
    "    color=COLORS[1],\n",
    "    label=r'$\\rho_k(\\alpha)\\geq\\rho_0^k$'\n",
    ")\n",
    "line_suff = mlines.Line2D([0], [0], color=COLORS[3], label='Thm. 9', linewidth=3)\n",
    "\n",
    "# handles is a list, so append manual patch\n",
    "handles.append(patch_pos) \n",
    "handles.append(patch_neg) \n",
    "#handles.append(line_suff)\n",
    "\n",
    "# plot the legend\n",
    "plt.legend(\n",
    "    handles=handles, loc='upper center', \n",
    "    bbox_to_anchor=(-0.01, 1.20, 1.0, 0), \n",
    "    ncol=2, fontsize=15, columnspacing=1\n",
    ")\n",
    "plt.tight_layout(pad=0.1)\n",
    "plt.savefig('bilinear_sweep_eg.png', format='png', dpi=200)\n",
    "print(f'beta1: {EG_SIM_BETA1}, beta2: {EG_SIM_BETA2}')\n",
    "print(f'=> c: {eg_sim_cs}')\n",
    "print(f'=> a: {EG_SIM_ALPHAS}')\n",
    "print(f'=> k: {EG_SIM_KS}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-8-f6fee4ad5edd>:14: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3.  Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading'].  This will become an error two minor releases later.\n",
      "  plt.pcolormesh(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rmin: 0.6418117802935815, rmax: 0.9999999914846447\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASUAAAErCAYAAACYbZUHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAbMUlEQVR4nO3dfXhU5ZkG8PtJJglJCCEgEj4kAYQEFhMgikCBWtqK2BVQLF+7ULXCLopVqyu6uqIrrkutu3IJgogfVbblUgqKWiyIWlfQFQwkBEiAliAYAiQgCZAAyTz7x5zAGCZhgMx73pncv+vKxXyceedJznDnnJNz3kdUFUREtohyuwAiIn8MJSKyCkOJiKzCUCIiqzCUiMgqDCUisorH7QIuVG5u7giPxzNLVVPBUCUKN14RKa2pqXmyf//+fw60gITTeUq5ubkj4uLi5qWnp5+Kj4+vjoqKCp/iiQher1eqqqpaFBcXx548eXJGoGAKqy0Nj8czKz09/VRiYmIVA4ko/ERFRWliYmJVenr6KY/HMyvgMqaLuhSqmhofH1/tdh1EdGni4+OrnUMw5wirUAIQxS0kovDn/D8OmD/hFkpEFOEYSkRkFYYSEVmFoRSmvF4vMjIyer/wwgttL+R1kydP7jJu3Li0UNUFAAUFBXGTJk1Ky8jI6B0dHZ0zYMCAjFC+X6Rrbus67E6eDOSZVRXZVafNfy/xMah5ZGSrPNPvCwCvvPJKSkVFRfS0adMOX8jrHnvssdKsrKw+BQUFpX369DkZito2b94c//HHHyf369fv2OnTp6Upx6797cpsVJ0y/7mNj62JfnAU13U9oVjXEbGl5EYgufm+ALBgwYL2t9566+G4uLgL+mtkRkbGqZycnGNz585tdynvv3Xr1riGnps4ceJ3paWl+atWrfpbz549qy7lfc7hRiC5+b5wf103JhTrOiJCKVwdOHAgety4cWmtW7fuGxsb279///6ZhYWFsed7XUFBQdymTZsSJ0yY8L3fnAcOHIieNGlSWkpKSnZycnLf6dOndwKAnJycjKeffvryuuVGjx59ZMWKFW1ra2svqN4TJ07IggUL2gwcOLDnVVdd1aeh5aKjoy9o3ObA1nW9Zs2axOeff75tcXFxzMV8X6FY1wwll+zfv98zcODAzH379sXNmzev+Lnnntuza9euFlOmTOl6vteuWrUqKT4+3jtw4MAzv5nKysqiBw4cmJmfn58wb9684unTp5cuXLgwdc6cOe2++eabuPvvv/9Q3bLDhg07Vl5e7vnqq6/ig6l1/fr18VOmTOmSmpqafe+996anpKTULF26dOfFfefNj83ruqKiInr27Nmdu3btmtWrV6/eM2bM6PThhx+2PH36dNN88xchIo4phaNf/epXnVu2bOn97LPPdng8vtVQUlISM2fOnE5Hjx6NSk5O9jb02tzc3MRu3bpV+/+WevTRRzuUlZXFrF+/vqBDhw41AI6+9NJL7WfPnt3p17/+9f6EhIQzm/5XX311VXR0NNavX584aNCggJvc5eXl0S+//HKbN99887Jt27YlZGZmVj300EMlU6dOLW/fvv2FbWI1czav67Fjx1bcfPPNeZ9//nnCypUrkz/66KPkBQsWpCYmJtYOGTKkYsSIEUfHjBlTkZaWZiylGEouOHToUPSKFSvavvHGG7vqPqQA0L1795OqimPHjkUVFRXF3XXXXV02btxYVP/1Bw4ciGnTpk1N3X2v14u333677aRJk8qcDykAICkpqbaysjL6gQceOOT/+piYGCQlJdWUlpYG3GRftmxZq8mTJ3ePi4vTMWPGlL/88svFgwcPbtpjQ82E7esaAKKiojBs2LATw4YNOwFgf2lpafTSpUtbP/vssx1XrVqVcv/998Pr9X59qT+LYDGUXPDee++1UlWMGjWqwv/x/fv3xyQkJHhTU1Nrli9fnpyVlXUi0OtPnjwp8fHxZ3675uXltThy5IhnxIgR3xuvtrZWpk6derBVq1bn/CaOiYnR6urqgLvvcXFx2qJFC291dXVURUVF9OHDhz1erxdRUdzbv1C2r2t/hYWFse+8807y6tWrk7/88suk2tpaGTx4cMX1119/NPjv+NLxU+aCDRs2JCQlJdW2bNnye39NWbNmTfKwYcOORkdHY9OmTQk5OTnHAWDWrFntr7vuuisPHDgQDQApKSm1FRUVZ36h7Ny5Mxbw/fb1GyuxpKQktl+/fgE/7JWVldH+v4H93XTTTZUlJSX5c+fOLS4tLY0ZNWpUzy5dulz14IMPdtixY8d5D87SWbav64KCgrg77rjjivT09D69evW66oUXXkjt2LHjqcWLF/+trKxs87p163bOmjXr4KX/JILHUHLBli1bEiorK6PLysrOHChYs2ZN4rp161pNnTq1DADy8/MTMjMzT44cObLboUOHPGvXrt1VdyynR48e1fv27TsTDnW7BYcOHfIAvk38hx9+uDMABJovq6SkxFNdXR2VkZHR4IwL8fHxOm3atCPr16/fuX379i1jx44tX7JkSbvMzMyrBg8e3HPBggVtmuanEdlsX9efffZZYmFhYYvbbrvtUG5u7tZvv/12y5IlS76ZNGnS0aSkpAaPdYUSQ8kFhYWFCW3btj3985//vOuyZctaPf3005ePHTu2x7hx48puueWWCq/Xi6Kiovjp06enjRs37siLL774rf+BzqFDhx7bv39/bElJiQcABg8efDwuLk5nzpzZedmyZa3GjRuX/t1333nS0tJOLl26tE1RUdH3tm4+//zzBBHB8OHDjwVTb0ZGxqm5c+eW7Nu3L//3v//9rsTExNp77rknvaHlKysro1577bWU1157LaW0tDT28OHDnrr7lZWVzeozZ/u6Hj9+/NH58+fvHT169FEA2LRpU4tAXw19f6FY1zymZNiePXtiysvLPUuWLNn1xhtvXDZ58uTuSUlJtXfeeefB3/zmNyWA78TEtLS0kwkJCd6amnO3um+88cbK5OTk2hUrVrS6++67D6emptbOmzdv96xZszpPnDjxyp49e1Z98MEHO9esWZM0c+bMLm+++eaJ2bNnH6h7/apVq5KvueaaytTU1Av6K5rH48GECROOTpgw4ejevXsb/OyUlJR47rjjjm7+j9XdHzx48JaMjIxTF/K+4Soc1vXvfve7lAceeOC8l6KoasAD3aFY12E1HW5eXl5xdnZ2Wf3Hw+kyk7feeqvV+PHjexw8eHBzu3btAn5QFi9enPLJJ58kPffcc98OHTo086WXXioePnz4cf9lbr/99it2794d9+mnn+66kPevqalBp06dsmbNmrXvrrvuuqDLFmwQTpeZcF03Li8v77Ls7Oz0+o9HxJaSW9efXYxNmzYldOzY8VRDH9K6Zfr27XsiNTW1dunSpX8dO3bslatXr97Rs2fPM791HnvssdI+ffr0yc/Pj8vKygr6uqZXX301pUWLFt6pU6da9yENhlvXn10MruuLExGhFE62bNmS0KtXr4B/Jakzf/78b+tu5+TkVBcXFxfUX6Z79+6n586dW7x3796YC/mgqqosXLiwOCbmoq4qoAvAdX1xImL3jYjCT0O7b83qLyFEZD+GEhFZhaFERFZhKBGRVcItlLxer7dJp1clIvOc/8cBL2MJq1ASkdKqqqoGT3knovBQVVXVQkRKAz0XVqFUU1PzZHFxcezx48fjucVEFH68Xq8cP348vri4OLampubJQMuE1XlKAJCbmzvC4/HMcvqQh1WoEhG8IlJaU1PzZP/+/f8caIGwCyUiimzc0iAiqzCUiMgqEXFBrohwH5QofJSpaoMNMiMilADg3969+LnNnxqdzDGacAwbauAYTT9GE9awp7FluPtGRFZhKBGRVRhKRGQVhhIRWYWhRERWYSgRkVWMh5KI3CAiRSKyS0QeDvB8ioisEJF8EflKRPqYrpGI3GM0lEQkGsB8ACMB9AYwUUR611vsXwFsVtUsAFMAzDVZIxG5y/SW0gAAu1T1b6p6CsBSAKPrLdMbwFoAUNVCAOki0t5smUTkFtOh1AnAXr/7+5zH/OUBuAUARGQAgDQAnY1UR0SuMx1KgSZmq3/d2n8CSBGRzQDuAbAJwDlN1kVkmohsFJGNTV4lEbnG6HxKIjIIwBOqOsK5/wgAqOozDSwvAHYDyFLVikbG5QW5ROHja1W9uqEnTW8pbQDQQ0S6ikgsgAkAVvovICKtnecA4E4AnzUWSEQUWYzOEqCqNSIyA8CfAUQDeFVVt4rIPzvPLwTQC8AbIlILYBuAXwYzdrhfgR1JY9hQA8do+jGasobGGJ+6RFX/BOBP9R5b6Hf7CwA9TNdFRHbgGd1EZBWGEhFZhaFERFZhKBGRVRhKRGQVhhIRWYWhRERWYSgRkVUYSkRkFYYSEVnF6CwBocJZAojCilWzBBARNcr4BbmhEu5XYEfSGDbUwDGafgxTswRwS4mIrMJQIiKr2Nj3LVlE3hORPBHZKiK3m66RiNxjY9+3uwFsU9VsANcBeM5velwiinA29n1TAElO04CWAA4jQDcTIopMNvZ9mwffPN0lALYAuFdVvfUHYosloshkY9+3EQA2A+gIoC+AeSLS6pwXqS5S1asbOwmLiMKP6VDaB+AKv/ud4dsi8nc7gOXqswu+vm+ZhuojIpdZ1/cNwDcAfgwAItIeQAaAvxmtkohcY2Pft6cAvC4iW+Db3ZupqmUm6yQi99jY960EwPWm6yIiO3CWACIyjbMEEFH44CwBsOMK7Egaw4YaOEbTj8FZAoioWWIoEZFVGEpEZBWGEhFZhaFERFZhKBGRVRhKRGQVhhIRWYWhRERWYSgRkVUYSkRkFRtbLP2LiGx2vgpEpFZE2piuk4jcYXTqEqfF0g4AP4VvatwNACaq6rYGlr8JwP2qOvw843LqEqLwYdXUJcG0WPI3EcAfjFRGRFYwPXVJoBZL1wZaUEQSANwAYEYwA4f7tBCRNIYNNXCMph/D1NQlpkMpmBZLdW4CsE5VDwccSGQagGlNVRgR2cF0KAXTYqnOBDSy66aqiwAsAnhMiSiS2NhiCSKSDOCHAN41XB8RuczGFksAcDOA1ap63GR9ROQ+61osOfdfB/C6uaqIyBY8o5uIrMJQIiKrMJSIyCoMJSKyCkOJiKzCUCIiqxidJSBUeEY3UVixapYAIqJGGT95MlTC/QrsSBrDhho4RtOPYWqWAG4pEZFVGEpEZBWGEhFZhaFERFZhKBGRVRhKRGQV6/q+Octc5/R92yoifzFdIxG5x+h5Sk7ft/nw6/smIiv9+76JSGsALwK4QVW/EZHLTdZIRO6yse/bJADLVfUbAFDVg4ZrJCIXmQ6lQH3fOtVbpieAFBH5VES+FpEpxqojItfZ2PfNAyAHwI8BxAP4QkS+VNUd3xuIfd+IIpLRWQJEZBCAJ1R1hHP/EQBQ1Wf8lnkYQAtVfcK5/wqAD1X17UbG5SwBROHDqlkCgun79i6AoSLicVp3Xwtgu+E6icgl1vV9U9XtIvIhgHwAXgCLVbXgfGOH+xXYkTSGDTVwjKYfw9QsAbb2fXsWwLMm6yIiO/CMbiKyCkOJiKzCUCIiqzCUiMgqDCUisgpDiYiswlAiIqswlIjIKgwlIrIKQ4mIrGJ0loBQ4SwBRGHFqlkCiIgaZfyC3FAJ9yuwI2kMG2rgGE0/hqlZArilRERWYSgRkVWCCiURuduZKfKSna/vm9Pz7ajT922ziDzeFO9LROEh2C2liQB2i8h9IhLv/4SIdBGRB4MZxK/v20gAvQFMFJHeARb9X1Xt63z9e5A1ElEECCqUVHUIgKcAPA6gWEQeEZEbReRnAP7DeTwYwfR9I6JmLKi/vonIMwAeAlAOoALAVADp8LVH2glgRpDvF6jv27UBlhskInkASgA8qKpbA9TEFktEESjY3bdpAB5V1ctVtYeqdgMwBMAXAGoA/CXIcYLp+5YLIE1VswG8AOCdQAOp6iJVvbqxk7CIKPwEG0peAN/rKKKq653dug0AGuzJVs8+AFf43e8M39aQ/7gVqnrMuf0nADEiclmQ4xNRmAs2lJYDeFxEUgI8txRAVpDjnLfvm4ikiog4twc4NZYHOT4Rhblgz+h+GL5dtb+KyCIAHwHYA6ANgH8F8G0wgwTT9w3ArQCmi0gNgCoAEzQSLtAjoqAEFUqqekRErgXwGIDb4TvorfAdI6oGMCXYNzxf3zdVnQdgXrDjEVFkuahZAkQkA0BXALXwXfF7uKkLu8B6uCVFFD4anSXgoi7IVdUiAEUXXRIRUQM4SwDsuAI7ksawoQaO0fRjcJYAImqWGEpEZBWGEhFZhaFERFZhKBGRVRhKRGQVhhIRWYWhRERWYSgRkVUYSkRkFYYSEVnFeCidr8WS33LXiEitiNxqsj4ictdFTV1y0W/ma7G0A8BP4ZsadwOAiaq6LcBya+Cbq+lVVV12nnE5dQlR+Gh06hLTW0rBtli6B8AfARw0WRwRuc/01CXnbbEkIp0A3AxgOIBrgh043KeFiKQxbKiBYzT9GKamLjEdSsG0WHoewExVrXX6BwQeiH3fiCKS6VA6b4slAFcDWOoE0mUAbhSRGlV9x38hVV0EYBHAY0pEkcR0KJ1psQRfB5QJACb5L6CqXetui8jrAN6vH0hEFLmMhlKQLZaIqBkzPkf3+Vos1Xv8NhM1EZE9eEY3EVmFoUREVmEoEZFVGEpEZBWGEhFZhaFERFYxOktAqPCMbqKwYtUsAUREjTJ+8mSohPsV2JE0hg01cIymH8PULAHcUiIiqzCUiMgqDCUisgpDiYiswlAiIqswlIjIKtb1fROR0SKSLyKbRWSjiAwxXSMRucfoeUpOP7f58Ov7JiIr6/V9WwtgpaqqiGQBeAtApsk6icg91vV9U9Vjevbal0Sc2+2EiCKY6VAK1PetU/2FRORmESkE8AGAOwzVRkQWMB1KwfR9g6quUNVMAGMAPBVwIJFpzjGnjU1bIhG5yegsASIyCMATqjrCuf8IAKjqM428ZjeAa1S1rJFluItHFD6smiXgTN83EYmFr+/bSv8FRORKcTpRikh/ALEAyg3XSUQusbHv21gAU0TkNIAqAOM1iM25cL8CO5LGsKEGjtH0Y5iaJcC6vm+qOgfAHNN1EZEdeEY3EVmFoUREVmEoEZFVGEpEZBWGEhFZhaFERFZhKBGRVRhKRGQVhhIRWYWhRERWMTpLQKhwlgCisGLVLAFERI0yfkFuqIT7FdiRNIYNNXCMph/D1CwB3FIiIqswlIjIKjb2ffsHp+9bvoisF5Fs0zUSkXuMhpJf37eRAHoDmCgivestthvAD1U1C76mAYtM1khE7rKx79t6VT3i3P0SQGfDNRKRi6zs++bnlwBWBXqCLZaIIpPpUwKC6vsGACLyI/hCaUig51V1EZxdO548SRQ5TIfSPgBX+N3vDKCk/kIikgVgMYCRqsr2SkTNiI1937oAWA5gsqruMFwfEbnMxr5vjwNoC+BFpydlTWPXyRBRZLGx79udAO40XRcR2YGzBBCRaZwlgIjCB2cJgB1XYEfSGDbUwDGafgzOEkBEzRJDiYiswlAiIqswlIjIKgwlIrIKQ4mIrMJQIiKrMJSIyCoMJSKyCkOJiKzCUCIiq9jYYilTRL4QkZMi8qDp+ojIXUanLnFaLO0A8FP4psbdAGCiqm7zW+ZyAGkAxgA4oqq/DWJcTl1CFD6smrokmBZLB1V1A4DThmsjIguYnrokUIula5ti4HCfFiKSxrChBo7R9GOYmrrE2hZL5x1IZBqAaZdWDhHZxsoWS8Fg3zeiyGRdiyUiat6sa7EkIqkANgJoBcArIvcB6K2qFSZrJSJ32NhiqRS+3ToiaoZ4RjcRWYWhRERWYSgRkVUYSkRkFYYSEVmFoUREVjE6S0Co8IxuorDS6CwBxs9TCpEyAHsucYzLnHHcxjrsqgFgHU1dQ1pjT0bEllJTEJGNjaU362ieNbAO8zXwmBIRWYWhRERWYSidtcjtAhys4ywbagBYh7+Q18BjSkRkFW4pEZFVmmUoicirInJQRAr8HmsjImtEZKfzb0qIa7hCRD4Rke0islVE7nWpjhYi8pWI5Dl1POlGHc57RovIJhF538UaikVki4hsFpGNLtbRWkSWiUih8xkZ5MJnI8P5OdR9VYjIfaGuo1mGEoDXAdxQ77GHAaxV1R4A1jr3Q6kGwAOq2gvAQAB3i0hvF+o4CWC4qmYD6AvgBhEZ6EIdAHAvgO1+992oAQB+pKp9/f707UYdcwF8qKqZALLh+7kYrUNVi5yfQ18AOQBOAFgR8jpUtVl+AUgHUOB3vwhAB+d2BwBFhut5F75+eK7VASABQC58HWaM1gHfxH5rAQwH8L5b6wRAMYDL6j1m+mfRCsBuOMd83aqj3ntfD2CdiTqa65ZSIO1VdT8AOP9ebuqNRSQdQD8A/+dGHc5u02YABwGsUVU36ngewEMAvH6PubFOFMBqEfna6ZjjRh3dABwC8JqzO7tYRBJdqMPfBAB/cG6HtA6GkstEpCWAPwK4T12ah1xVa9W3id4ZwAAR6WPy/UXk7wEcVNWvTb5vA36gqv0BjIRvl3qYCzV4APQHsEBV+wE4DnO7rudwmnyMAvC2ifdjKJ11QEQ6AIDz78FQv6GIxMAXSP+jqsvdqqOOqn4H4FP4jreZrOMHAEaJSDF8XZOHi8gSwzUAAFS1xPn3IHzHTwa4UMc+APucLVYAWAZfSLn12RgJIFdVDzj3Q1oHQ+mslQB+4dz+BXzHeEJGRATAKwC2q+p/uVhHOxFp7dyOB/ATAIUm61DVR1S1s6qmw7eb8LGq/qPJGgBARBJFJKnuNnzHUQpM16G+5hl7RSTDeejHALaZrsPPRJzddUPI6zB1oMymL+cHvB/Aafh+K/0SQFv4DrTudP5tE+IahsB3/CIfwGbn60YX6sgCsMmpowDA487jRuvwq+c6nD3Qbfpn0Q1AnvO1FcCjbv0s4PtL6EZnvbwDIMWlOhIAlANI9nsspHXwjG4isgp334jIKgwlIrIKQ4mIrMJQIiKrMJSIyCoMJSKyCkOJrCMiE0TkpHPGOzUzDCWyUTaAbap62u1CyDyGEtmoL3xnmVMzxFAiG2XDd9kNAEBEhotIuYg8LyLR7pVFJjCUyCoi0g6+icM2O/dnAHgfwCOqep+q1rpYHhkQKW27KXL0he9C5W0isgjAWAA/U9VPXK2KjGEokW2yARwGsBy+GQ2vVdVd7pZEJjGUyDbZAATAUAA/YSA1P5y6hKwiIlvgm2mxO3xzK+Wo6iFXiyKjeKCbrOHMBZ0J38Rm/wTf5Plviwi36JsRhhLZ5O/gO6SwRVWrANziPPbfrlZFRjGUyCbZ8HXu+CsAqOoeAOMBTBeR21ysiwziMSUisgq3lIjIKgwlIrIKQ4mIrMJQIiKrMJSIyCoMJSKyCkOJiKzCUCIiqzCUiMgq/w8X7/qcvJR1awAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "GD_ALT_KS = list(range(5, 80, 5))\n",
    "GD_ALT_ALPHAS = list(reversed([x/100 for x in range(5, 95 + 1, 90 // 15)]))\n",
    "GD_ALT_HPS = np.meshgrid(GD_ALT_KS, GD_ALT_ALPHAS)\n",
    "\n",
    "# compute\n",
    "gd_alt_radii = np.array([\n",
    "    gd_alt_radius(k, alpha) for\n",
    "    (k, alpha) in zip(GD_ALT_HPS[0].flatten(), GD_ALT_HPS[1].flatten())\n",
    "]).reshape(GD_ALT_HPS[0].shape)\n",
    "\n",
    "# and plot (matplotlib choses the first color when there is only one value)\n",
    "cmap = matplotlib.colors.ListedColormap([COLORS[0], COLORS[1]])\n",
    "plt.figure(figsize=(4, 4))\n",
    "plt.pcolormesh(\n",
    "    GD_ALT_HPS[0], GD_ALT_HPS[1], gd_alt_radii < 1, \n",
    "    edgecolors='k',\n",
    "    cmap=cmap,\n",
    ")\n",
    "plt.xlabel(r'$k$', fontsize=15)\n",
    "plt.ylabel(r'$\\alpha$', fontsize=15)\n",
    "\n",
    "#\n",
    "handles, labels = plt.gca().get_legend_handles_labels()\n",
    "patch_pos = mpatches.Patch(\n",
    "    color=COLORS[0], \n",
    "    label=r'$\\rho_k(\\alpha)<1$'\n",
    ")\n",
    "patch_neg = mpatches.Patch(\n",
    "    color=COLORS[1],\n",
    "    label=r'$\\rho_k(\\alpha)\\geq1$'\n",
    ")\n",
    "line_suff = mlines.Line2D([0], [0], color=COLORS[3], label='Thm. 6', linewidth=3)\n",
    "\n",
    "handles.append(patch_pos) \n",
    "handles.append(patch_neg) \n",
    "#handles.append(line_suff)\n",
    "\n",
    "# plot the legend\n",
    "plt.legend(\n",
    "    handles=handles, loc='upper center', \n",
    "    bbox_to_anchor=(-0.01, 1.20, 1.0, 0), \n",
    "    ncol=2, fontsize=15, columnspacing=1,\n",
    ")\n",
    "plt.tight_layout(pad=0.1)\n",
    "plt.savefig('./bilinear_sweep_gd_alt.png', format='png', dpi=200)\n",
    "print(f'rmin: {gd_alt_radii.min()}, rmax: {gd_alt_radii.max()}')"
   ]
  }
 ],
 "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.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}