binary-search-on-alpha.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "import networkx as nx\n",
    "import numpy as np\n",
    "import scipy\n",
    "import pandas as pd\n",
    "import random\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "from data_helpers import make_polarized_graphs\n",
    "from core import query_graph_using_dense_matrix, sweep_on_x\n",
    "from helpers import flatten\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "random.seed(12345)\n",
    "np.random.seed(12345)\n",
    "\n",
    "size = 10\n",
    "k = 2 \n",
    "g, true_comms, true_groupings = make_polarized_graphs(k, [(size, size) for i in range(k)])\n",
    "A = nx.adj_matrix(g, weight='sign')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# check if the correlation constraint is tight (from SDP)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARN: m less than n, problem likely degenerate\n",
      "alpha= 0.198893609536022\n",
      "beta= 1.0038596859017648\n",
      "[13 18 17 10 11 14 12 15 19 16 30 39 21 35 33 36 37 31 38 34 25 22 20 26\n",
      " 24 29 32 23 27 28  8  9  0  5  2  1  7  6  3  4]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.sparse.linalg import eigs\n",
    "\n",
    "import cvxpy as cp\n",
    "\n",
    "from helpers import (\n",
    "    signed_laplacian,\n",
    "    prepare_seed_vector,\n",
    "    degree_diag,\n",
    "    is_rank_one,\n",
    "    sbr_by_threshold\n",
    ")\n",
    "\n",
    "kappa = 0.7\n",
    "\n",
    "seeds = [\n",
    "    [4], [13],\n",
    "]\n",
    "\n",
    "\n",
    "n = g.number_of_nodes()\n",
    "L = signed_laplacian(g).A\n",
    "D = degree_diag(g).A\n",
    "\n",
    "s = prepare_seed_vector(seeds, D)\n",
    "\n",
    "# requirement check\n",
    "_, v1 = eigs(np.real(L), k=1, which='SM')\n",
    "assert not np.isclose((s.T @ np.real(v1))[0, 0], 0)\n",
    "\n",
    "DsDsT = (D @ s) @ (D @ s).T\n",
    "\n",
    "X = cp.Variable((n, n))\n",
    "constraints = [\n",
    "    X >> 0,\n",
    "    cp.trace(D @ X) == 1,\n",
    "    cp.trace(DsDsT @ X) >= kappa\n",
    "]\n",
    "prob = cp.Problem(cp.Minimize(cp.trace(L @ X)), constraints)\n",
    "\n",
    "opt_val = prob.solve(eps=1e-10)\n",
    "\n",
    "print('alpha=', constraints[1].dual_value)\n",
    "print('beta=', constraints[2].dual_value)\n",
    "assert np.isclose(constraints[2].residual, 0), \"correlation constraint should be tight!\"\n",
    "\n",
    "x_opt = X.value[0, :]\n",
    "D_sq = np.diag(np.sqrt(np.diag(D)))\n",
    "# normalize to such that x^T D x = 1\n",
    "x_opt = x_opt / np.sqrt(x_opt @ D @ x_opt[:, None])\n",
    "\n",
    "# x^T D x = 1\n",
    "assert np.isclose(x_opt @ D @ x_opt[:, None], 1)\n",
    "\n",
    "# s^T D x = 1\n",
    "assert np.isclose((x_opt @ D @ s)**2, kappa, atol=1e-3)\n",
    "\n",
    "print(np.argsort(x_opt))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## conclusion\n",
    "\n",
    "- correlation constraint is tight"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# solve linear system usnig $\\alpha^{*}$ (numerical stability issue could arise)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = g.number_of_nodes()\n",
    "L = signed_laplacian(g)\n",
    "D = degree_diag(g)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "y/x_opt [[ -1.96139505  -1.4524639   -2.1238421   -1.47171285  -0.16989237\n",
      "   -2.06202635  -1.55462491  -2.13929172  -2.79051998  -2.07333373\n",
      "   -1.29733388  -1.47127753  -1.89263965   0.06284787  -1.94613782\n",
      "   -1.8347658   -5.35781326  -1.49651444  -1.52940226  -1.1050006\n",
      "    8.30868916 -10.11231598   8.78899135   1.1801487    6.45196518\n",
      "   11.58674796   9.17816119   1.14380505   0.0609099    8.67645437\n",
      "    0.96192606  16.85314891  -9.46436153   3.00677005  15.52706097\n",
      "    7.91917134  12.56938482  13.18641059  14.73612463   1.56945348]]\n",
      "y [[-0.0574527  -0.0580015  -0.07641299 -0.06339318 -0.0275918  -0.06883934\n",
      "  -0.06314604 -0.08681728 -0.06516751 -0.05643026  0.0498793   0.05500156\n",
      "   0.06466768 -0.01044185  0.07144404  0.05782207  0.06849034  0.05941406\n",
      "   0.06218344  0.02328673  0.00994615  0.01359103  0.01028192  0.0045221\n",
      "   0.00884245  0.01195837  0.01203348  0.00447351  0.00024927  0.0123085\n",
      "  -0.00375429 -0.01382279 -0.01417022 -0.00304701 -0.01104062 -0.01059764\n",
      "  -0.01267915 -0.01238104 -0.01176287 -0.00457377]]\n",
      "node order [[ 7  2  5  8  3  6  1  0  9  4 32 31 36 37 38 34 35 13 39 30 33 28 27 23\n",
      "  24 20 22 25 26 29 21 19 10 11 15 17 18 12 16 14]]\n"
     ]
    }
   ],
   "source": [
    "from numpy import linalg as LA\n",
    "from scipy.sparse import diags\n",
    "from scipy.sparse.linalg import spsolve\n",
    "b = D @ s\n",
    "\n",
    "alpha_opt = constraints[1].dual_value\n",
    "A = L - alpha_opt * D  \n",
    "\n",
    "# is it numerically stable? as we assume L - \\alpha^{*} D has rank n-1\n",
    "y = spsolve(A, b)\n",
    "\n",
    "# normalize so that y^T D y = 1\n",
    "y /= LA.norm(y, 2)\n",
    "y = diags(1 / np.sqrt(D.diagonal())) @ y[:, None]\n",
    "\n",
    "assert np.isclose((y.T @ D @ y)[0, 0], 1), 'y not normalized w.r.t D'\n",
    "\n",
    "print('y/x_opt', y.T / x_opt)\n",
    "print('y', y.T)\n",
    "print('node order', np.argsort(y.T))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lambda_1= (-0.8519099893484773+0j)\n",
      "lambda_n= (15.92531716412138+0j)\n",
      "condition number = 18.693661728630186\n"
     ]
    }
   ],
   "source": [
    "lam_1 = eigs(A, k=1, which='SM')[0][0]\n",
    "lam_n = eigs(A, k=1, which='LM')[0][0]\n",
    "print('lambda_1=', lam_1) # <--almost 0 because of rank reduction\n",
    "print('lambda_n=', lam_n)\n",
    "print('condition number =', abs(lam_n / lam_1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# relation  between $\\alpha$ and correlation constraint residual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(40, 40) (40, 1)\n"
     ]
    }
   ],
   "source": [
    "from scipy.sparse.linalg import eigs, spsolve, cg\n",
    "\n",
    "s = prepare_seed_vector(seeds, D)\n",
    "\n",
    "n = g.number_of_nodes()\n",
    "L = signed_laplacian(g)\n",
    "D = degree_diag(g)\n",
    "\n",
    "# lb = - g.number_of_edges() * 2\n",
    "lb = -2\n",
    "lambda1 = eigs(L, k=1, which='SM')[0][0]  # THIS IS WRONG!!!! Use the **normalized** Laplacian\n",
    "ub = np.real(lambda1)\n",
    "\n",
    "b = D @ s\n",
    "gamma_list = np.linspace(lb, ub, 1000)\n",
    "gap_list = []\n",
    "y_list = []\n",
    "print(A.shape, b.shape)\n",
    "for gamma in gamma_list:\n",
    "    A = L - gamma * D\n",
    "    if True:\n",
    "        # conjugate gradient solver\n",
    "        y, info = cg(A, b)\n",
    "        if info != 0:\n",
    "            print('error, info=', info)\n",
    "    else:\n",
    "        # linear system solver\n",
    "        y = spsolve(A, b)\n",
    "        \n",
    "    y /= LA.norm(y, 2)\n",
    "    y = diags(1 / np.sqrt(D.diagonal())) @ y[:, None]\n",
    "\n",
    "    assert np.isclose((y.T @ D @ y)[0, 0], 1), 'y not normalized w.r.t D'\n",
    "\n",
    "    gap = (np.sqrt(kappa) - y.T @ D @ s)[0, 0]\n",
    "    \n",
    "    y_list.append(y)\n",
    "    gap_list.append(gap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7f02bb11b5f8>"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEPCAYAAACneLThAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXmYVNWZ/79vb6y3WWWHBtmxEgkqoD/UnjETQVHMaBJ1\nYhxcY1wzGkmiEyGaRNEZNRo1bokmOpiJkwQXNCo2S2RTVCh2EZpV2eGyNd3V5/fHW4dz6tat6tq6\ntn4/z9NP1b11l3Or4HvPfVdSSkEQBEEoTkpyPQBBEASh+RCRFwRBKGJE5AVBEIoYEXlBEIQiRkRe\nEAShiBGRFwRBKGIyIvJE9BwRfUlEy5rY7jQiaiCif83EeQVBEIT4ZGom/zsA58bbgIhKANwP4K0M\nnVMQBEFogoyIvFJqPoC9TWx2M4A/A9iRiXMKgiAITZMVmzwR9QJwEYCnAFA2zikIgiBkz/H6CIAp\nytRQEKEXBEHIAmVZOs+pAGYQEQHoCmACEdUrpWbaGxGRFNIRBEFIAaWU7+Q5kzN5QowZulLqxPDf\nALBd/gdegbe2Lfi/e+65J+djkGsonusohmsoluvI5TUcOKDwzjv8d+BA5GfxyMhMnoheBlANoAsR\nbQJwD4AK1mz1tFfHM3FOQRCEloLrAqefDqxYwcvDhwOLFgGO0/S+GRF5pdTlSWx7VSbOKQiC0FII\nBoGVK83yqlXA4sXAOec0va9kvDYD1dXVuR5C2hTDNQDFcR3FcA1AcVxHrq4hEAA6dEhtX2rKnpNN\niEjl03gEQRDyAdcF+vUD9u3j5aFDgSVLjLmGiKBiOF6zFV2TFv3790dtbW2uh9FiqKqqwsaNG3M9\nDEEQwgSDRuABoL4+8X0LYiYfvkvlYEQtE/m+BSG/cF2gRw/g8GFeLi8H5s4Fxo7l5XgzebHJC4Ig\n5DmOA4wcCZxwAi+PGAGcdFJi+4rIC4IgFACNjcBTTwGtWgEPPZT4fiLygiAIBcDhw0C3bkBdHTBh\nAnDmmWzGaQoR+Txh8+bNqKysjGkLnzZtGq644oqMnKukpASff/55Ro4lCELz47rA7t3AZ5/xckMD\nx83r5Kh4FER0TUugb9++OHDgQNxtuPRP+mTqOIIgND+uC5xxBrB1K/Czn7HTVanE7fIyk88woVAo\n10NoEomcEYTC4dNPOYQSADZv5tfHHgPmzUusrEFRiLzrAgsWJGafao5jDBgwANOnT8fJJ5+M9u3b\nY/Pmzbj44ovRrVs3DBw4EI899tjxbZcsWYLTTjsNHTp0QM+ePXHHHXcAAGpra1FSUoLGxkYAwMaN\nG1FdXY0OHTrg3HPPxa5du44fY86cOejbt2/UGGbPnn38HGeccQY6deqE3r174+abb0ZDQ0PyFyYI\nQs5p0yZyub6enbCJCDxQBCLvuuyAOOusxB0RzXGMGTNmYNasWdizZw+++c1vYtSoUdi+fTvee+89\nPProo3jnnXcAALfeeituu+027N+/H+vXr8e3v/3t48ewzSiXX345TjvtNOzatQt33303XnjhhYjz\nxTO5lJaW4pFHHsGePXuwYMECzJ49G0888UTyFyUIQk5xXeA73zHLZWVA+/Zsn09Upwpe5INBdj40\nNPBjTWUlQJTcX2Ul75uMM8PLrbfeil69emHZsmXYtWsX7rrrLpSWlqJ///645pprMGPGDABAeXk5\nPvvsM+zevRtt27bF6NGjo461adMmfPjhh/j5z3+O8vJynHnmmbjgggsSHsuoUaMwevRoEBH69euH\n6667DnPmzEn+ogRByCnBIGAnn//850BJCTBtWguKrgkE2PlQXg6cfDJw4AA7JZL5O3CA9y0vTy7J\nwKZPnz4AWKC3bt2Kzp07o3PnzujUqRN+9atfYccObm37/PPPY82aNRg2bBjGjBmDN954I+pY27dv\nR6dOndDGek6rqqpKeCzr1q3DBRdcgJ49e6Jjx4646667Isw9giAUBoEAZ7pqnnmGhT0USnxCWvAi\n7zjsgJg7N3FHRHMcQ5tP+vbtixNPPBF79uzBnj17sHfvXuzfvx+vvfYaAGDgwIF4+eWXsXPnTtx5\n55245JJLcOTIkYhj9ezZE3v37o1Yv2nTpuPv27Vrh8M6vxns7N25c+fx5RtuuAHDhw/H+vXrsW/f\nPvziF78QZ6sgFCCOA0yaZJY3b2bLQ2lpC4uucRyu4ZCKOGfyGAAwevRoVFZWYvr06Th69ChCoRBW\nrFiBDz/8EADw0ksvHZ9Vd+jQAUSE0tJSACbqpV+/fjj11FNxzz33oL6+HvPnzz9+kwCAIUOG4OjR\no5g1axYaGhpw33334dixY8c/d10XlZWVaNu2LVavXo0nn3wyvYsSBCFntG9v3o8YAVx+OXD99WZC\n2pTJpihEPtfYTtCSkhK89tpr+OSTTzBgwAB069YN11577fEY+LfeegsnnXQSKisr8cMf/hCvvPIK\nKioqoo7z8ssvY+HChejSpQvuvfdeXHnllcc/q6ysxBNPPIGrr74affr0geM4x81FAPDQQw/hpZde\nQmVlJa6//npceumlMccrCEL+4rrA735nlt94A+jaFeje3Qj8aafFP4ZUoRSikO9bEPKDBQuAceM4\nZBIAPvgAeP994OBB4Je/5M/PPBMIhaQKpSAIQsERCLAN3l5u3Ro4etQse1JmosiIyBPRc0T0JREt\ni/H55UT0KRF9QkTziegrmTivIAhCMeM4XNKgY0ez3KaNEXnHAaZOjX+MTM3kfwfg3Diffw7gLKXU\nSAD3AXgmQ+cVBEEoag4fBrp0McutWwN2QF5TLraMiLxSaj6AvXE+X6iU2h9eXAigdybOKwiCkAyZ\nKIGSyjnffZf/Us3It002trkGiHzvRy6qUF4DYFYOzisIQgtGly9ZsYLjy1PNiUn2nKefbpKWAgF2\nniZz3gMHgJ49zbJSwLZtfGzHyTORJ6J/AjAZwLhY20y1DEzV1dWorq5u9nEJglD8BIPAsmUskjpb\nVPdIbc5zrl5tlletSv68+/dzsxCAhf3uu4ENG2oweHANJk8Gwik4McmayBPRVwE8DWC8UiqmaWeq\njxehqqpKYruzSDIlFAShUAgEWOABYOjQ1MqXpHLOgQOBtWt5uawM6Ncv8f1dF9i5k/8AYNEigJPf\nq7FnTzUmTeJaNu++Oy3mMTIp8hT+i/6AqB+AVwFcoZRan+yBN9oVegRBEFKkvJxL9WYLxwEefNCU\nJmhsZJHu1Sux/YNBc2MC2Mnavz/w+eec/dqvH/Dss/GPkakQypcBfABgCBFtIqLJRHQ9EV0X3uQ/\nAXQG8AQRfUxEizNxXkEQhEQJBo3Ar1mTWrXZVBg1il+Jki+AGAjwa0lYqUeP5mbeI0awT6G2Fmiq\n9mBGZvJKqcub+PxaANdm4lyCIAipEAhwvPm+falXm00FXXuma1fgzTeTc7o6DhcjO/dcs2/Hjhxh\n4zh8TZ06cX35WEjGqyAILQLHAXQZp/vvz955devmnTuB885LLYyye3fzvqzMPJE4DjBhQvx9ReQF\nQWgxaLPHBRek3gUuWVatMu9TaUqkVGTCU3k5NzjSNNVWWkReEIQWg7Zfp9MFLlmGDTM3l8GDkzcT\n+Ym87TzOqzh5QRCEXNK5M7+WlWXPLt++PWestm8P/Md/JL+/UuYmASQv8jKTFwShxRBu3YDLLkve\nCZoqeia+Zw/w/e+nbyYSkRcEQYhBXR0L7ksvpe4ETRal2G5+5EjqZiLbXGM7XgEReUEQhOPs2sWi\n29iYPZu8UizMyfZm1ft6kZm8IAg5qbZYCHTuzDHmRNkrbaAUi/ukScAdd0T2ZvX7jez1sUTejq45\ndCj++UXkBaHI0JUPzzyTG06I0BuUYpFMpbtlOjdOIr7BdOtmBP6UU6J/owMHONnprLP4swMHeN9Y\n0TWuC2zYEP/cIvKCUGQsWsRmiFCIU/kXSxGR4+zebWa+yZQ2cF3OLtXim4zQ6xuKFneAf6N166J/\no48/ZlNSQwOPLRiMbgpii7xdqiEWIvKCUOQcPpzrEeQPnTsDJ5zA76uqEq8IGQxyYbFUHKc6uqa8\nnKtRxrtB2B2gysq4f6tX5EtL+eZw4ADfsEqaUHEReUEoMsaMATp0MMtTpojJRlNeDtxyCwvnxo2J\nR9jYhcKSja9Xiv9++1tgxgx+EhgxgscCAMOHc+ExwJQkBljIN22KFnkiFvpx47ikQWNj/POLyAtC\nkeE4LPSaVavEZKNRCti7l1+TmZXrePrRo5OPr9chlF98YaJ6Vq0yMfulpWbbPn34Vd9Mhg/37+Fa\nWgosXx7pgI2FiLwgFCG2cAgGpYDevYG2bZPLerVt6cnG12tnb58+ptywUsaMZvsG9Kz83/+do3Da\nt48WeddNria+iLwgFCEdOxpbbXk523YFFtc2bTii5frrE5+VB4Nmf7/Zf7zIGx1C+eCD/IQ1bx6/\n6hLEgwaZG82OHfyqSwl769bosWhnbmkpX088ROQFoQixY6zr64Hzzxe7PMDfSV0d8NFHwJNPJm+T\nB6Jn/zryZty42JE3ROxU1eJtm9Sef95E3tx6K6/7wx+4WbefyNtj6duXn0riISIvCEVIp07cpEKz\ncWP2OiHlO1u3cnPsZLJe9Wy/b1+TzKTRkTexjqdvtm3acGkDjX7SWrqUBV4fB+Dls882cfJ+YwGA\nH/1IRF4QWiRlZTwr1AJRWppcA+liRSm2jXfrllqkTGNjtHknEADateP3fsfTs3GvyGu7+i238BNA\nly6RkTK1tXzT8HO8alauNA7cWIjIC0IRYkeRAGyisJtXFCrplmvQNvmf/ASYODF6Vt4UfiUEHIdt\n/ID/8WKJvL6GUIjFes6cyP2qqrgWvTfj1ebJJ/nJJB6ZauT9HBF9SUTL4mzzayJaR0SfENHITJxX\nEAR/GhujhaHQk6JcF/jqV1PLOtVowe3cmXu9JsvBg/7n7diRX/1uGLFEvl07dorromXnnWc+a9WK\nRT9WdI2msTF7Bcp+B+DcWB8S0QQAA5VSgwFcD+CpDJ1XEAQflAIGDgSGDDHrfvrTwna+BoPsW0in\nq5NSLIpTpwJz5yZ+s9DbNDT479OqVeR23nMSmSxVvU1pKXDiicBNN/ETQM+evL53by6e1qtX7Oia\nZMiIyCul5gPYG2eTSQBeDG+7CEAHIuoeZ3tBENJAKXbI/eY3Zt3q1YXtfA0EWPDS7eq0bRvbu4HE\nbxa2sHr3cV1g1ix+73cD0BmvF17IIq+3UYozk3XRMm2P79jR2Oubiq5JhGzZ5HsD2Gwtbw2vEwSh\nGdDiMGKEieIoKyts56vj8DU8/3zytnSNTobSTziJ3iwCAfM9evcJBtn/AcS+adTX803W3qaxkc0x\n2s7vF4XjJ/L2dbdtG98xC2Svx6vfMHyLfU6dOvX4++rqalRXVzfPiAShiNHiUFtrxKOhgUP0evXK\n7djSQSl2SKbatk8/4cydC3Tvzq+JHMtxOPpl587oG0wgwEXPvvwydnRNRQV/tmwZlyo46SRe3769\n8ZXomXzr1sbO7ifyTA2AGlRW8sx/27bYY8+WyG8BYOfc9QHgOyxb5AVBSA3d/LmqyrSLK/SZPMBC\nuGgR8LWvpT6TtxOTFi0Cxo5N7Fg6s9S7reNw79Zp00yEzIIFLP46a7W0FJg/n/MXZs8269u1i57J\nt27Nwr9gAdvp/UW+GkA1unbl2P1t26bFHHcmzTUE/xk7AMwE8D0AIKKxAPYppb7M4LkFQbDQ0TX2\nTL6+vrDDKF2Xr+snP0k/usZ1gWPHOKIl0WPpekB790aHcWrH69GjnMlqRwDpczoOb9e6tRlLRQX/\nRnaGclkZL591Fod5xuPIkSyVGiailwF8AGAIEW0ioslEdD0RXQcASqk3AWwgos8A/BbADzJxXkEQ\n/NHCEghwpAbA0R0331y4ETba+anjylONriEyzTaSidTRInz22dFhnPqzpUv5Rmof1za5lJfzzQXg\n87/xBvDOO3ys/ft5/dGjpkrm2rV8vbFIROQzYq5RSl2ewDY3ZeJcgiA0jT17/OUvub8oYMoOn3NO\nbseXCtrWnUqmqo2++bVvz4Ka6LG0kOsCYVrEx4419vT+/dnmf/iwOe7WrZEiryNnDh8G9uwxx1q5\nktdXVpooosGDuURxLLI2kxcEIb+wZ49N1TYpFHTVxn/7t/SiawDe9xvfAO6+O/Fj6X39omzsqpCj\nRvF7+7j6t6ioMCLfujWHT+ooKP3EVV7OQv/668DMmfHLRh850nRZaRF5QShCbJEfMYJnhQDbhIcP\nz9240kHPlg8eTP0Y9vfSuTPHrSezb0mJiU6yyxRrkT9yxNSx8X4GRJprSkqAO+7gUsO6drymTRtT\nEydeiOTRo/HNOYCIvCAUDXZdFy1IQHQYZaE6X7Ug//WvmXG8zpwJPPJI4sdqbOSb5JfhkBG7TLG+\nAS1eHG0+8drk7USnLl2MWc2+GbRpY2zz8WrXAEBNTfxxi8gLQhHgusBXvmIcgnV1Rhi8ztfbbitM\n5+vy5fwaq3FHItiO1127kis3HArx96pn4vZ+Oq79hhuADz7wPyfA5hq9v46T1yGUdgXK1q35qSB2\nnLyhqd9SRF4QioBgENi82UR12HXIHQd49FGzbaGWNxgxgl+1DTsVx6sdddS7d3JO3Pr6SCGuqoru\n6BQK+Zc10L9FSQnw4YeRT1u6nk28mXw8Kivjfy4iLwhFQCDAdVC0aDlOpDgUQ3kDbeseODB1xytg\nzCP33sv+CW1bb6qMcWlpZO329983Y+jSxWzToUPkfraJaP164Npr+Wnr6FHgzjt5vddklMxMfvz4\n+J+LyAtCEeA4bCO+9loWwNLSSHHwK29QaNjjT6esAcCCet99/ERz3nlcFmDgwPhljIl49q4LhHXq\nZG4M+rh33cWtFvU59Dm1iejIEfO0tX9/ZKE0HUKpa94nOpPXyVWxEJEXhCKhtBTo0cO/AbQubwAU\n7kxem0r27Uu9cYgtuBs38rqVKzkpaefO+MlRjY0cjlpezstffAGcdhrfGP7wB17XpQvw1lv83r5Z\n6BuE3amrbVtgwABeHjTI+E0A/o0+/pj3b0rkJU5eEFoI2jEIRIt8MZQ3UIqjW/btS71xiF8m8IgR\nPPuuqIhfxliXIdBO1o8/Btas4RvDrl287vPPOcEJiMx4BaKfpkIhrqg5YAC3atQ3sfp6YOFCLt9w\nySWRtnobXUpBRF4QWgiNjZEib//nL4YIm8ZG8zSSauMQOxN43jz+jmpqOPa9Z0/gwQdj2/uV4lm8\n/o61mJeUGJu863L8PWBuFvaNRRc5GzKERZqInbY33ghcdBF/tn8/O2NDIS5rcPiwicixEZEXhBaG\nLfLe9n+OwzHhmkKMsFGKRb6sjMU2lQgb+wmnUyf+XhYtMtEtdhljryNWz+R1rfcbb+TXyy/nGTcR\n8Pvfm3Ppm4V9Y9GJaP/7vyzOn3/OIh4K8VMBwNEynTrx50RcEO0vfzFj0mhbvGS8CkILwTuT99py\nR4wwglCIdvnGRh5///7Aiy+mHmGjvxfXZcGeOJFNP/X1pra76wLjxkVXk6yoMNs0NPBrWRmbkJTi\nMeoZvi4n4U2GAvg6Ghu5No0OgbQzVy+5BPjWt8w59DFtE5XM5AWhhdGUyNfWGrtvIUbY6Gvq1Int\n2OnUrgHY+XrsmDH91NWZxKRgkJ90bLOQzni1xR3gph1durCAl5QAXbvyejvpSf8W+vXwYZMM9f/+\nX+QYDx7ka/N7SrFNVCLygtDCaErkq6rMTLJQZ/IlJTxDXrgwvegawFSi1M5WgEXUdfkz/f3YtvWK\nCvMd//Sn/KrNRwMGABMmAFdcwev9RF6zZAnfLHTUjU379pFhke3b83i8Jio98xdzjSC0EBob4wuL\nPZMvxAgb3RD744+B229PL7oG4NnyhAkc2/7mm+zwfOwxc9xrruHtbNt6ebmZyffowa+hEH+vlZW8\nnb6R2s5SfU79/f/gB/x7HD4cnTxVXm66QwF8E/ngA25VOG+e2W7DBn7V44mFiLwgFAnHjnHstrdA\nmSYQAIYN4/eF2ECksZH/Dh5MvXGI9+bXvTtHw9TWGrFesYJt8XffbbZzXTPzBliYdTXMUMi08tNZ\nqkDkDVejnbY63HX9ehNxY2OLfGMj3zx0m8JgMPK31Q7bWIjIC0IR4LrcQ3TxYp6JHjsWPZN3HOAX\nvzDLuoFIoaBn0p07p944xCvy7dsDn3zCJpOSEv7r358TpbQ46+/0yBHg7bd5XWWlsd9rkW/blrfR\ns3W/pyrtjC0tZZv64MHRIq9UtMjbBAJcjE77BN57L/41i8gLQhEQDPLMUldoPHjQP1OykBuI6Oia\nq68GJk9Or34NwDfGF18EXniBSxuUl3OkzZw5kf4KpYBPP+X3e/fya+vWHMMOmCeARERem3J++lMu\nkOY4/r9J69b+1SkBE+P/xBP+n3sRkReEIkA3mAB4htu2rb/IF3IDEW2C6taNhS7V6Br9vQSDXBte\nlxsOhfiYvXoB999v9hkzxrzX+372GfDKK/xem2jatePZdzyR16+9e8euE69bAtpPCl4cB7j0Uv4t\ns9XIezwRrSaitUQ0xefzvkQ0m4iWEtEnRDQhE+cVBIFxHODUU4ETT/QvUKbRBbGAwmsgohO8KipY\nZDMRXdOvnyld3NhoTCT27Nq+mWgB1zHxALB7t0nU2r3blD2IJ/L6xnDoEDB9euQY58/nCpmbN0ee\n04vjcJbulVfGv+a0RZ6ISgA8DuBcACcBuIyIhnk2uxvAK0qpUQAuA/BEuucVBCGSkhIWd78CZZpA\ngFPqgcIrb6Bt5A89xAXFMhFd8+yz7IyeN4/FVM+eY6GfgnQ2qj5OXR3w2mtsy//Tn3j9Rx9FFiiz\nX7XIr1vHjb5tGhtZ4D/+mJfjtfdr0yY6BNNLJmbyowGsU0rVKqXqAcwAMMmzTSMAXdq+IwDPZQmC\nkC7xatdoCrmBSGMjh35u2ZJ6dyjvza9zZw6d1LNlPZOPhRbcykpg5Eh+v3s3RzXt3s3LOjv1hhuA\nM86I7EnrFflhwyLt/61aRce9x7O5V1Rkx1zTG8Bma3lLeJ3NNABXENFmAK8DuDkD5xUEwcKOk/fW\nrrEp1AYiOupEl+dNNLrGW4PGLmvw3e9yLfmzzuJ1O3bw8urV/sfSAnzggAmHfOcdPv4JJ/Cybsjd\n2Mh2/zfeMEJti7wOjXzCsmv0788RPPoam6JVq6ZFviyxQ8XF75+StzjmZQB+p5R6mIjGAvgj2LQT\nxdSpU4+/r66uRnV1dQaGKAjFT7xSwzZ+DUR69crOGNNBR9f8z/8AF1+cWHSN6wJf/SrP/k86yZQc\nAFiAdYx5MMiva9dy8xDdbNvLiBHAsmUs6FrkleKngRtuAH77W56961BLAPj1r9msYteG37CBfy8i\nU8ES4Oibc87hEsSXXGKeDrzU1NSgpqYGX37JHarikQmR3wLAngv0AbDNs83VYJs9lFILiag1EXVV\nSu3yHswWeUEQEqepsgYa3UCkvr7wZvIlJXxDqq9PLLrG2xzka1+LdLwOG8bb2CYR7Tj1Y84crp1z\n443Ac8/xOiJOjho2jIVbR//Y1SuPHGHTkj72//0f37AOHYp08uqxdekS+0YDmAnwnDnAN74BvPfe\ntJjbZsJcswTAICKqIqIKAJcCmOnZphbA1wGAiIYDaOUn8IIgpE6iIl+oDUS0CaqyMvHuULpVH8Cz\n8HbtIh2mdnJYIuiKkR07Gvv9uHHcdKRNG/6rr+fa8N27A3368DZt2vDNdPlyXlaKn6LWrYtMhtKm\nFzsZKh6JmGvSFnmlVAjATQD+DmAFgBlKqVVENI2IJoY3uwPAtUT0CYCXADQR9CMIQrI0NvJMUmdg\nxhL5Qm0gomfyjY08M06kO5Se7ffo4R9a6k1EKi2NL5r2vvv382uHDiZevU0bXv/22xyDv2ULbxMI\nmFo19rmGDo0UeX381q2brkkDZM/xCqXUW0qpoUqpwUqp+8Pr7lFKvR5+v0opNU4pNVIpNUop1UQi\nriAIyaJNDnV18UW+UBuI6Jm8Hmsy3aHatPEPLR0zhmPNdWJSvHBFwOz74IPGyX30qBlbmzbsvN25\nM3K/Y8dY6AMBFvc2bTiyp7IyUsxtkU8Ev2gcL5LxKghFghb5efNYOOLN8AqxgYieyeuyu379WL2R\nNBr7u/B2zPrZz1jsO3Ximb03ZNGvx6od275/v7l5tGnDr15Hdtu2phzB448DJ5/MYzp0iE093nMl\nKvIlJZEJbr7bJHYoQRDyHe2omziRZ+fxbLqF2EBEz5Ydh52czz4bGWHjusCoURzd4jXjaGH3E+wu\nXUxNeL86Mn5x6raIt25tbkCtWnFEzJ13AmefDVx7LW+ja9bozORjx3ifNWsifSLaMZuIyLsuh2/a\neQ9+iMgLQpFgt6U7csTUG/ejEBuI2AleXbrwmO0Im2CQe6YC0WYcW+S9Zqy2bVkwS0tN/R8bP5G/\n4w5TT16XODh6lM+/YQPwX//FNwLtePWeT3eGGj488klE32TKEoh7DAYjE7liISIvCEWCLdqtW3O8\ndywKsYGIneDVrp1pwK0JBExCkteMk6jI65Z63vN62bTJ3GD0rHzrVhM7v3Urm2L0jdd+gmjTxmS8\nahOO5ujRyHj6eOiSw/p3j4WIvCAUCVoYfv97YNAgk3nph7eBSCFE2OiZvOtyp6S77oo0yzhOdDcn\njX4C8BN5pUzlx82bEYXtjNXnevRRs60WefvJondvnun71YS3Z/La/KRZuZKvaZsn08jvt9E3iLlz\noz+zEZEXhCJB1zQfNiyygJYfjgM8/LBZLoQIGz2Tt80UtlnGdU1UizdRKpbIuy5nqu7bx/Vn9Ezc\ne15NMMhPSqGQsZ/rEgVt2wJf/zqv+/nPOXLGbyYfCnE9G92z1svKlVwKwSbWb6M7RsVDRF4QioTG\nRjYF6BDKpuKn7Qib0tL8t8vb0TUnnBDZHcp1OWLl6ad5W+/M1xZ2+71tx29o8G/FZ4c4BgJ8vrKy\nyBIFuhOX7tfarh2bfuzSB3pcEybw+v37Wey9Yx06NDLiBki+A5aNiLwgFAmNjWyLr6uLX6BMoxtl\nALxPvtvHxJu2AAAgAElEQVTl7eiaW28FvvUtY5YJBiNDCb0zX3smb2MnhgH+Med2eQFvVyaABX73\nbj6HdtzqSBuvuSYY5O9dr1uzxjwdaB55JDoEM50OWCLyglAkeGfyiTjvbBJJo88l9tNJjx4mwQlg\nse7WzWzrnfkePmwanHvj5H/5S7NslwXWeGvI6K5Mts/j00/ZfNOUyOsnAc3w4ZHrKiqA0aOjx5CO\nv0REXhCKBG0X1o7ApkR+zJjI9n8//Wl+O1/tp5OKCjaz2E7XH/zAbGvHzgPA+vXs0KyvT/7mpzNb\nvdhPBa7LJQz0jPzo0cj6M3pb/STQuTMvV1ZGRtjo4mZeUmmQohGRF4QiIdmZvONwGVxNvjtf7eia\nn/2MhdEWP78Yd11CGGAzyeHD/iGU8fCrBqkbp9vH6NaNSwQDfMNUyt/x6jhc4AwwTyZa2GOVVUil\nQYpGRF4QioRQiGOmg0FTq7wpCqm8gR1ds2lTdHcovwQibxVKXXbAZswYU1OmKZu8pqoq0rF96qkc\n2fPll7y8eTPHyut68F5fgL6xeMdin8seS6INUvwQkReEIqGhgW3D993HIpiIjb2QyhvY0TWDBvE6\nW/z8RF7PkHUfV78G547DcfdPPulf9sDPXFNbGx2K2b8/0LcvL/fuzc1NtDPYe6PQUTzesdjnsq8n\nkQYpsRCRF4Qiob7exF/X1bEduikKqbyBHV0zaxbbs23xi1cKwHFMFcpYn196KVek9GILtDYN6YqS\n5eXmxtGuHfCb3/Dn3/8+sH272c8bf19Rwa9eR68+l+vG9gUki4i8IBQJRBzxUVLCIjJ4cNP7FFJ5\nAzu6pndvFsjly43wxhP5eBmvGsfhJCYv+/aZ99oHYGeb9ujBN1Ui41A9diyybo0WdYD3X7aM3//L\nv0Q6VLXIL1oUeUNavDj2tTWFiLwgFAGua5JrOnRg4Unk8b6QyhvY0TVHj7IInn22Ed54Ih+vdo2N\nLjpms2SJeW/7AHS2abt2bBo7cgSYPJk/u+eeSJG27fe209br7E6kUUiyiMgLQhGg+5TqptKJhgoW\nUnkDW6CDQdNCTwtvIh2dmhJ5v2JfdtvpoUOjHaAVFcDevdzK77PPzHq75rxdGiEQMNE1w4aZ42lH\nq+saZ7DGL3Y+UUTkBaEI0LPx0lIWkPLyxOPBR4wwApnP5Q3sWi+BACcb2Y1D4pXc1eUD4om86wI3\n3xy93rbJP/JI5BOS63K8/t69XFJhyBDzWf/+5r03hPLSS/l9TY1p+q3Hf+aZ/PrBB1yqQe+TKglU\nLRYEId8pKWF7/O23s1190aLERX7lysjWgatWRafV5wO2QDsOMHIkm0cuv5yX7Rhzva22py9fzuIZ\nr3BbMOjvrC4rM2YU74w6GIwsL/zWW+azdu2A008347Hp2pVfdWNw/WQCmCeTsWO5W1W6ZGQmT0Tj\niWg1Ea0loikxtvk2Ea0gouVE9MdMnFcQBObQIRb5oUM5qiMUarpAWSzytbyBt+ha587s3NSzXHsm\nr5tsf/KJWaeToWIRCPg7q+2yB94ZdSBghLpfP74JnHMO/2knrHdsQHR0TSBg1jUVFposaYs8EZUA\neBzAuQBOAnAZEQ3zbDMIwBQApyulvgLgtnTPKwiC4eBBTrC5+25g/nwum3voUGL7Fkp5g4MHeWZu\nZ7h++KFZtq/3rLN4vS3aI0ZwqYF40TXf/370+j17Yo/JcYBvfpPfP/xw5E3AbkBii7zrAr/9Lb8/\n+2wTraN/g0TDQhMlEzP50QDWKaVqlVL1AGYAmOTZ5loAv1FKHQAApdSuDJxXEIQwhw7xLFdngjY0\nsCMwEQqhvIHrct/UxYtNU4333+eQRx1dYzf80Nege6X278/i2VSd/Usuie6vOn16/LHp8sLe49qR\nMra5Jhg0MfSrVpnv2q/JS76IfG8Adj+VLeF1NkMADCWi+UT0ARGdm4HzCoIQZscOftzXLf/KyiJL\n6DZFvpc38JYyeOMNdnbajUPs8EcdBaOdpmVlJhkqnsj36hUZIQPEd+i6LvDKK/z+hhtMRyfXBS66\nyGxn+wsCAU5CAyLr4S9fzuvsejx+ZRaSJROOV7+vzJtXVgZgEICzAPQDMI+ITtIze5upVrxSdXU1\nqqurMzBEQShedHejjRv5kV8LWjLVFv3KG+ST8zUQYBHfvp2F8fzzuSDYjh1GKO2Eoddf5+9AO151\n9mgihdt6h6eoRLEzZDXBIJvGAHa8nn02sHQpr1+92mxXV2feOw7w1FPAN74BvPkmLy9YYMxNtuM1\n1ky+pqYGNTU18QcXJhMivwUs3Jo+ALb5bLNAKdUIYCMRrQEwGMBH3oPZIi8IQtPYs9y1a3nWuHcv\n8K//CixcmFj4XVUVPwnU1bFJw3Ya5gOOw5E0mzcDjz3Gy1OmcAjiH/4QHV0TDPI11NezbdyvyFgs\n9CxaC3xlJXDggPnM/j6rqiKjb2prWaB1jfhPP+X1tli7Ljc9AYDzzmMzki6TsHJlYo5X7wR42rRp\nMa8nE+aaJQAGEVEVEVUAuBTATM82fwXwzwBARF3BAv95Bs4tCC2eQIBn3SUlbHvW4rR2beK29dpa\nI1THjvFMOd+cr/v2RdqtO3SIrP1im1UuuojNHnv3soM2mZl8MBhpJjntNBP54q3rXlsbOduvqmKB\ntmvEA9EtB7W/RM/a7TIJeed4VUqFANwE4O8AVgCYoZRaRUTTiGhieJu3AewmohUA3gNwh1Jqb7rn\nFgSBBeG224BJk9hWbQtKojPyQAAYMMAsb9yYX85X1wVeeokrRWrH6733ArNnG+G1Z/I6E3bVKo46\n0jP5RETeLj5WWck3CL2/t6673rasjCtjzpljBNqe8ds3Aj3LLy+PnLXrMgn2fvrGlc4NNyPJUEqp\ntwAM9ay7x7N8O4DbM3E+QRAiaWhgp6s22+h1EyZwrHhTJhvH4RvE8OEsLPnmfA0GeSZvO143b46s\nKW+LvM6E1Q2/tU080WYq8+bxMR94gAW2Y0d+avDWdbe31TN4P+yKkonu47omuerMM1MvNyxlDQSh\nCDh0yL8z0qZNic/Id++OvEHkU235QIBn5Fq8zz8/uqa8ba55+WV2at50E19HXV3TZQ1s9Ky6fXsu\nhjZpUrQpxbutd/02yzP5xReRy7H2sQkGjS9AOkMJQgtnzx5g1y4WPDuxyS6A1RTa+Qrkn/PVcXh8\nzz/PQturF/Duu9x846GHeBtvmGJtLZudNLq8bzJRRxUVJgEr3qzbj9dfj1x+883E9wUizUbpdIaS\n2jWCUOC4Lnch2ruXs13ffdfUhR89OnFh8nO+Ll2aXnGsTOJ1vDoO142ZMIEFcPx481lDAwtkv37A\nhg08gx8ypOmQSC9E7MBetYrj55MxmUycyIlVR4/y8nnnJXfuRM06TSEzeUEocIJBnsnrxKBNm0z9\nlGSEIZ+dr67LMfLf/rZxtOqxaSer7q8KsJnDcYD77+eyDe3aseAmaq6xz1tfH5l0lSi9epmCZ126\npJZ3kIhZpylE5AWhwPHaq1N9rNfOVy2C+VR2+KOPouvHe8sN27P8yZNNI5Hu3dns8o9/8P7JiLxu\nB1hSktp3q4U9FMpdSKqIvCAUOI7D8fHPPptew2cgMjpHlx3OB3r2ZMG27dOOA5x6Ktd4nzePzVW6\nSuWGDXwj0MldrgtccEHiDc412kdx662pfbda2Pfti46xzxYi8oJQBOzZY0reZpJ8KTu8ezeXFZ41\nK1Jsu3QxpQu6dOFZO5FJSqqrY7t9QwP/1dVF16aJhesCL77I7199NbVxB4PmfToRMukgIi8IBc6B\nAxyiZ9urUyUfyw67LnDFFewYvv32yPX/+Afws5/xdYdCwI03cjclXfqgro5r3GizTqtWiTU4B1ig\nd+7k91u3ptZM227hl44pLR1E5AWhwPnww2h7darkY9nhYNCUD7CvLxiMrES5cyc7WDt3NmaWujq+\npmHDuD1fnz7+JX390OUigNSbnDsO+0uA9E1pqSIiLwgFTo8e0fbqdLDLDueD81VXoPQ6P73r27Y1\nNwM7FLSigk1ZJ57IppxkHK92gtWaNand8HS9+VyFokqcvCAUOHv2sBA//XRycfGxWLnSJBblQ89X\nxwG+9z02ST36aGRtmClTOLv1xReBM84wBcn0jLuujk00bdtyVnAyIZTBYGRYprbzJ4vORPZWsMwW\nMpMXhALGdVkAN2yItFdnknxwvu7ZY2bEGtflyJrFizlxS+cKHDzI3wdgRL5dO76OZEQ+XvGxRHFd\ntucDEl0jCEIKxLJXp4PX+TplSm6dr7r70mOPRQplMGiKlG3cyEJeWspCrJ88Dhzg2TgR13YPhRIX\neZ1xOm8eZ/6m8jQTDMauYJktxFwjCAVMIAB07cp1a3TLu3TRmaKTwp2aV63i2fI556R/7FQIBoH9\n+yNvZGPH8rUPGcLr+vdnm/tZZwFLlrB/wnWBP/6Rwy9LSoC//IVfP/iAI24SmZXrjNNU0WWFvc1A\nsonM5AWhwKmvT74mS1PoiJB8wJvZatdf16V4N27kFnqDB3PRslDIlHtQipcbG9khO3ly9kwnsZqB\nZBMReUEoYOw666lGf/gxYoQJQywvB/r2zcxxU8FxgE6duEmI7omq2bKFXxsa2Bb/xRemHV8gwFE1\nJSWRHZZCoeyaTjJRfyYdROQFoYAJBFiMS0szZ64B2M6vI2zq67nSY67s8jrZ6z//kys52uMIBLjw\nWGkpvw4bZkTecYBx44D77uOEKX3TAjL7XeU7IvKCUMAoFdl1KFMEAhwyqPn889QyPjPB4sXxk73K\ny/lzpdisU1oaGSc/dCjw1FPmeyot5aicfCmh3NyIyAtCAfP228bmnElzjeMADz+cmWOli7bH+yV7\nBYMc/97YyHXbt2zh9+vW8Yz/0CGOrtm+3ewTCuXW/JRtMiLyRDSeiFYT0VoimhJnu0uIqJGIRmXi\nvILQknFdDm/UZNoEceqpubfLuy47Shsa+MnCa5MPBLiPKxHfCFwX+NvfgP/6L+NctcsTaObOze51\n5JK0RZ6ISgA8DuBcACcBuIyIhvls1x7AzQAWpntOQRBMjDzQPCYIr13+/POzb5cPBk0Lv9pa/76z\nrVoZc85NN0U2UNm9m5ORZs40N6zWrZPv0lTIZGImPxrAOqVUrVKqHsAMAJN8trsXwAMA6jJwTkFo\n8VRVmcSe8vLIBKZMkA+dogIBLh8cq2lHMBjZIFvflEpK+Mnmyy+5guVVVwHLlwPPPMPdmnJZpiHb\nZELkewPYbC1vCa87DhGNBNBHKZVkK1tBEGJhz7QbGvxnuemQD52iHAe46CLgwgujTTWASYiyOeEE\n4OKL2aeg68ivXMkz/GuuaVkCD2Qm49UvSfh4agYREYCHAVzZxD4AgKlTpx5/X11djerq6rQHKAjF\niB39UlbWPALs7RT10UfZFUndpHz/fq5H400ochwOrbzsMrPu/PM5metrX+PlTFXnzCdqampQU1OT\n0LaZEPktAOx/Xn0AWA9QcMC2+pqw4PcA8DciulAptdR7MFvkBUGIzULLu6Vn8s0twLfcAlRXZy/8\n0E72sksa2OzYEbl87Jix0XfpArz+umkXWCx4J8DTpk2LuW0mzDVLAAwioioiqgBwKYCZ+kOl1AGl\nVDel1IlKqQFgx+sFfgIvCEJiuG5k1clhw5pnpjpmDNeF0WzcmN14+USSvS65hJ2vAL+OG8cVJzdv\n5qJlxSbwyZK2yCulQgBuAvB3ACsAzFBKrSKiaUQ00W8XxDHXCILQNMGgscE3Z3KP43AN91yiqzjG\nolcvTtZ65hl+7dKF/RXf+Q5/R7kq8ZsvkMp0ZaM0ICKVT+MRhHxl2zaOW29s5Nnr5583n6lm2zaO\nsjl2jO3bwWC0s7O5mDsXOPtsfl9ezsvxqkK6LjByJH8fmkT2K3SICEop38mzZLwKQgGycqVpTae7\nNzUXuaxjU1aWXGtD+wkH4MigllSnxg8ReUEoQOzY8OYmV3VsXJdj3GNlu/phh1Qm08u1mBGRF4QC\nw3WBu+4yy8OHc2/X5sJxgF/8InLdrl3Ndz5NItmuXhyH6/kApmhZJmv6FCIi8oJQYCxaZOqoA8AD\nDzR/9EjXrpHLP/5x85tsdNcrnb2aqMll8+bIZaLsJ3HlEyLyglBg5KKxdq5CKXXMezJ06RK5HApl\nPhu4kBCRF4QCo85T/SkbrfpyEUpp93ZNxuTiTQTt2VMcr4IgFAiuC9xxh1lubnu8zamnmqSjbJQe\nDgT4PMl2vZo4kStNApxI9d57kgwlCEKBsGhRpOkhG/Z4TW2t6biUjVDKUCi1rle9enGlyWee4Xo3\n2Yrpz1dE5AWhgMhGVEsssh1K+frr/JpK16tevVpmxUk/ROQFoUBwXY5qscmGPV6TzVBK1wXuvtss\nt/SEpnQQkReEAsGbzXniidmzx2uyFUoZDJpQyJbWeDvTiMgLQoHg7QQ1a1b2hS9boZTN3fWqJSEi\nLwgFwqJFpl5NfX100k82cBzgV7+KXNcc42jurlctCRF5QSgAXBe44YZcj4Lxmmyuuy7ztXSy0fWq\npSAiLwgFwKJF3JRa079/9u3xmjFjgB49zHJ9PfDyy5k9h1/XKyE1ROQFoQDwRrHcf3/uHJGOw45Q\nm8ceS90B67rAggVm/2x1vWopiMgLQp7jusB//EfkOq/JJNucd16kCWXTptQcsK4LfPWr3LJPd3DK\nVterloKIvCDkOYsWAdu3m+Vcmmo0jsPZtjapOGCDQXayNjaaRt1ee7xE1qSHiLwg5Dn5ZKqxyYQD\nNhAw16I7P2Wz61VLICMiT0TjiWg1Ea0loik+n/+QiFYQ0SdE9A4RNXNpI0EoDvLRVKPJhAPWcYAz\nzuD38+bxsu1gFtInbZEnohIAjwM4F8BJAC4jomGezZYCOEUpNRLAqwAeTPe8gtASyEdTjcbPAfvo\no8k7YMvLzfG8pRuyWWWzWMnETH40gHVKqVqlVD2AGQAm2RsopeYopY6GFxcC6J2B8wpC0eOd1eaL\nqUZz3nlAnz5mecsW4NVXUz9eLrpeFTuZEPneAGyXyxbEF/GrAczKwHkFoahxXeDWWyPX5YupRuM4\nwM03R6675prUk6NyWWWzWCnLwDH8eqL7Nuwiou8COAXA2bEONnXq1OPvq6urUV1dnd7oBKFAefNN\nYPdus9yvX36aLr77XeAnPzHO0lAIeP75yCqSieBtiAJkt8pmIVFTU4MabwusGJBKtoGi9wBEYwFM\nVUqNDy//GIBSSj3g2e7rAB4FcJZSanf0kQAiUumORxCKAdcFBg4Edu406x58MFoE84Xf/x6YPNks\nl5ZylEwiDTsuvBB47TXgb38DJlmG3v79gWXLxFyTCEQEpZTfhDsj5polAAYRURURVQC4FMBMzwC+\nBuApABfGEnhBEAyLFkUKfHk5cPnluRtPU1x8MdCtm1kOhUxyU6J4TVO//rUIfCZIW+SVUiEANwH4\nO4AVAGYopVYR0TQimhjebDqAdgD+l4g+JqK/pnteQShmNmyIXH766fzucuQ4LMo2O3Yk5oTVJYU3\nbjTr+vcHxFKbGdI212QSMdcIArB2Lddrsf8rvPsucM45uRtTIrguMHhwZERQSQknM8Uz20yaBMyc\nGbluxgzgO99pnnEWI81trhEEIUO4LvBP/xQp8L165afD1YvjAHPnsrBrGhuBsWPjR9vs2xe9Lt+i\niAoZEXlByCPefz9SEMvKeF2h2KaHDAGeey5y3d69wMiR/kK/di3fGGx69CiMm1qhICIvCHmC6wI3\n3RS57plnEotQyScuvhjo3j1y3c6d3Ix76VKzbts24PTTo/eXqpOZJRNx8oIgZID334+u5Ni3AKs8\nabPN8OEmdh4ADh4ETjkFOPdcoHVr4J13gMOHI/cdPJizaIXMIY5XQcgDXJerMNop/f36cSneQp3V\nLl0KnHpqpH8hHtddBzz0UOFeby4Rx6sg5Dlvvhkp8ADw+OOFLXijRgEffmhCJONRWgrcc09hX2++\nIiIvCDnGdYFbbolc16NHccSJjxrFNy9vopNNaSl3lcrnPIBCRsw1gpBj0ikJUEisXcvmmP37ebm8\nnPMBrrpKBD5d4plrROQFIYf4JT7de2/yxb2Elo2IvCDkIa7Ls/UvvjDrSku5ibXMbIVkEMerIOQh\n778fKfAA8OyzIvBCZpGZvCDkgG3bgJNPjmyS0acP2+IlwkRIFpnJC0Ie4brcBNvbBemJJ0Tghcwj\nIi8IWebVV6Nj4oslZFLIP8RcIwhZxC+aJpFyvIIQDzHXCEIesG0bl921BZ4IWLJEBF5oPkTkBSEL\naDv83r2R659/nrNCBaG5EJEXhCzwxz9G2+G7deOyvILQnIhNXhCambffBsaPj1wndnghkzS7TZ6I\nxhPRaiJaS0RTfD6vIKIZRLSOiBYQUb9MnFcQ8h0/gQfEDi9kj7RFnohKADwO4FwAJwG4jIiGeTa7\nGsAepdRgAI8AmJ7ueQUh34kl8HPnih1eyB6ZmMmPBrBOKVWrlKoHMAPAJM82kwC8EH7/ZwB53nde\nEFLHdYH//m9/gX/rLeDMM7M/JqHlkon2f70B2E3LtoCF33cbpVSIiPYRUWel1J4MnF8Q8oa1a4Ez\nzgB2747+7K23uPWdIGSTTIi8n7Hf6z31bkM+2wAApk6devx9dXU1qiUNUCgQ5s0DzjrL/zMReCGT\n1NTUoKamJqFt046uIaKxAKYqpcaHl38MQCmlHrC2mRXeZhERlQLYrpTq5nMsia4RCg7XBZ55Brj9\ndv/PReCF5iZedE0mZvJLAAwioioA2wFcCuAyzzavAbgSwCIA3wIwOwPnFYScE888QwTMmSM2eCG3\npO14VUqFANwE4O8AVgCYoZRaRUTTiGhieLPnAHQlonUAbgPw43TPKwi5xHW5auTQof4Cf8IJwOrV\nIvBC7pFkKEFIAtcF/vxn4Ic/NL1KvTz8MHD11VI2WMgezW2uEYSiJxFxF/OMkI+IyAtCHBIRd4DN\nM/PnSxarkH9IgTJB8GHtWuB73wN69wauuir+7P3JJ4H160XghfxEZvKCEGbbNuC554CFC4E332x6\n+x/9CLjtNmm8LeQ3IvJCi0YL+0cfAX/7W2L7iLgLhYSIvNDiWLsWmD4d2L49sRm7RsRdKERE5IWi\nxnWBN94AZs0CDh8GDhwA/v73xPefOJE7Ol11lYi7UJiIyAtFhZ6l798P1NcDs2ez0CeD43DHpp/8\nRJypQuEjIi8UJLaYa5Kdpds4DnDhhfw3YYIkMgnFg4i8kJdoM8t77wFt2wI7d/LMHEhPzG3Gjwd6\n9GBRF2EXihUReSHr6IiWVauMcNukamZpiokT+YbRsSNXjBRTjNASEJEX0sbPdBKLw4eTi2hJlfHj\ngcpKoLwcGDZMHKdCy0VEvoXijTpJlUyZTlJFi7lGZumCEImIfI5JdBZcXs71UWzbdKo0lzmkubj6\nauDIEaChwawTMReExMg7kb/yyvgzy0yKXa7J9Sw41zgO8PWv82/qRcwsgpAZ8q6efIzWr0IB4DWd\nxKJtW4loEYRMIvXkhbjoqJNUEdOJIOQvIvJ5QlOz4PJyoGtXNlPZtulUEXOIILQM0hJ5IuoE4BUA\nVQA2Avi2Umq/Z5uTATwJwAEQAvBLpdSfYh3zyivZyRaLTItdrpFZsCAIzUlaNnkiegDAbqXUdCKa\nAqCTUurHnm0GAVBKqfVE1BPARwCGKaUO+ByvKHq81tTUoLq6OtfDSItiuAagOK6jGK4BKI7ryNdr\niGeTT7cz1CQAL4TfvwDgIu8GSqnPlFLrw++3A9gB4IQ0z5vX1NTU5HoIaVMM1wAUx3UUwzUAxXEd\nhXgN6Yp8N6XUlwCglPoCTYg3EY0GUK5FXxAEQWhemrTJE9E7ALrbq8Bxjncnc6KwqeZFAFcks58g\nCIKQOuna5FcBqFZKfUlEPQC8r5Qa7rOdA6AGwC+UUv8X53iFb5AXBEHIAc0VJz8TwL8DeADAlQCi\numQSUTmAvwJ4IZ7AxxukIAiCkBrpzuQ7A/gTgL4ANgH4llJqHxGdAuB6pdR1RPRvAJ4HsALG1PPv\nSqllaY9eEARBiEtelTUQBEEQMku60TVpQUTTiWgVEX1CRK8SkW/OJxGNJ6LVRLQ2HI+fVxDRJUQU\nJKIQEY2Ks91GIvqUiD4mosXZHGNTJHEN+f5bdCKivxPRGiJ6m4g6xNguRERLw7/FX7M9Tj+a+m6J\nqIKIZhDROiJaQET9cjHOeCRwDVcS0Y7wd7+UiK7KxTjjQUTPEdGXRBTT2kBEvw7/Dp8Q0chsji9p\nlFI5+wPwdQAl4ff3A/iVzzYlAD4DZ9WWA/gEnEyV07F7xjgUwGAAswGMirPd5+CEsZyPOZVrKJDf\n4gEAd4bfTwFwf4ztDuR6rMl+twBuAPBE+P13AMzI9bhTuIYrAfw612Nt4jrGARgJYFmMzycAeCP8\nfgyAhbkec7y/nM7klVLvKqUaw4sLAfTx2Ww0gHVKqVqlVD2AGeAkrLxBKbVGKbUO7HOIByHHT0+x\nSPAa8v63QAIJemHyzcmfyHdrX9ufAZyTxfElQqL/PvLtu49AKTUfwN44m0wCh4NDKbUIQAci6h5n\n+5yST4JzFYBZPut7A9hsLW8JrytEFIC3iWgJEV2b68GkQCH8Fokm6LUiosVE9AER5cONKpHv9vg2\nSqkQgH3h4Id8IdF/H/8aNnP8iYj8Jnb5jvc6tyL//h8cp9mrUMZJprpLKfVaeJu7ANQrpV72O4TP\nuqx7ixO5jgQ4Qyn1BRGdAOAdIloVnjVkhQxcQ77/Fskk6PUL/xYDAMwmomVKqQ2ZHGeSJPLderch\nn21ySSLXMBPAy0qpeiK6Hvxkkm9PJE2RF/8PEqXZRV4p9S/xPieiKwGcB+CfY2yyBYDtYOoDYFtm\nRpc4TV1Hgsf4Ivy6k4j+An68zZrIZ+Aa8v63CDvMuiuToLcjxjH0b7GBiGoAfA1ALkU+ke92Mzhc\neeG077IAAAMGSURBVBsRlQKoVErFMytkmyavwTPeZ8A+lEJjC/h30OTk/0Gi5Dq6ZjyAOwFcqJSq\ni7HZEgCDiKiKiCoAXAqeDeQrvvZGImpLRO3D79sB+AaAYDYHlgSxbKaF8FvoBD0gdoJex/D4QURd\nAZwBYGW2BhiDRL7b18DXBADfAjvJ84kmryF849VMQu6/91gQYv8/mAngewBARGMB7NMmwrwkx17s\ndQBqASwN/+nIgZ4AXre2Gw9gTXj7H+faW+1zHReBZ1lHAGwHMMt7HQAGgKMNPgawPN+uI5FrKJDf\nojOAd8NjfAdAx/D6UwA8HX5/OoBl4d/iU3ByXj6MPeq7BTANwMTw+1bg5MN14ECF/rkecwrX8Evw\n5OZjAO8BGJLrMftcw8vgmXkdOMlzMoDrAVxnbfM4OJLoU8SJqMuHP0mGEgRBKGLyKbpGEARByDAi\n8oIgCEWMiLwgCEIRIyIvCIJQxIjIC4IgFDEi8oIgCEWMiLwgWBDRhqbqwSSyjSDkCyLyghBJIokj\nklwiFAwi8kKLhYj+Eq4IupyIrtGrw59VhRva/JGIVoYrJra2trmFiD4KN4EZEt7nNCL6R3j9fCIa\nnIPLEoQIROSFlsxkpdRpAE4DcKuPCWYogMeVUiMAuAB+YH22Qyl1CoCnAPwovG4VgDPD6+8B8Ktm\nHb0gJICIvNCSuY2IPoFpWDMYkaaYTUqpheH3fwR3DNL8Jfz6EbgTEgB0BPBnIloO4GEAI5pr4IKQ\nKCLyQouEiM4Gl7ceo5QaCS4e1zr+XhE3AF01NQRTsvteALOVUl8BcEECxxOEZkdEXmipdACwVylV\nR0TDAIwNr7fLy/YjojHh95cBmJfAMbeG30/O2EgFIQ1E5IWWylsAyoloBbj87Qfh9fZsfQ2AG4lo\nJYBOYPu7dxub6QDuJ6KPIP+3hDxBSg0Lgg9EVAWuo/+VXI9FENJBZhuCEBuZAQkFj8zkBUEQihiZ\nyQuCIBQxIvKCIAhFjIi8IAhCESMiLwiCUMSIyAuCIBQxIvKCIAhFzP8HgHYO7pUanMQAAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f02bb0d9dd8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df = pd.DataFrame.from_dict({'y': y_list, 'alpha': gamma_list, 'residual': gap_list})\n",
    "df['abs_residual'] = df['residual'].abs()\n",
    "df.plot(x='alpha', y='residual', style=['.-'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### questions\n",
    "\n",
    "- why the fluctuations from 0.0 upwards?\n",
    "- why the best $\\alpha$ here is different from the $\\alpha$ in SDP?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>alpha</th>\n",
       "      <th>residual</th>\n",
       "      <th>y</th>\n",
       "      <th>abs_residual</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>536</th>\n",
       "      <td>-0.183790</td>\n",
       "      <td>0.000621</td>\n",
       "      <td>[[0.03156807432367406], [0.04112115231024713],...</td>\n",
       "      <td>0.000621</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>535</th>\n",
       "      <td>-0.187179</td>\n",
       "      <td>-0.001451</td>\n",
       "      <td>[[0.03138253376730166], [0.04093080458090798],...</td>\n",
       "      <td>0.001451</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>537</th>\n",
       "      <td>-0.180402</td>\n",
       "      <td>0.002731</td>\n",
       "      <td>[[0.03175543312517406], [0.04131290404695293],...</td>\n",
       "      <td>0.002731</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>534</th>\n",
       "      <td>-0.190567</td>\n",
       "      <td>-0.003486</td>\n",
       "      <td>[[0.03119879457733433], [0.040741857883688865]...</td>\n",
       "      <td>0.003486</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>538</th>\n",
       "      <td>-0.177013</td>\n",
       "      <td>0.004880</td>\n",
       "      <td>[[0.03194462688327459], [0.0415060621549289], ...</td>\n",
       "      <td>0.004880</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>533</th>\n",
       "      <td>-0.193955</td>\n",
       "      <td>-0.005484</td>\n",
       "      <td>[[0.03101683973298032], [0.04055430867099873],...</td>\n",
       "      <td>0.005484</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>539</th>\n",
       "      <td>-0.173625</td>\n",
       "      <td>0.007069</td>\n",
       "      <td>[[0.032135672115284425], [0.0417006283455375],...</td>\n",
       "      <td>0.007069</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>532</th>\n",
       "      <td>-0.197344</td>\n",
       "      <td>-0.007446</td>\n",
       "      <td>[[0.030836652094924023], [0.0403681528598232],...</td>\n",
       "      <td>0.007446</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>540</th>\n",
       "      <td>-0.170236</td>\n",
       "      <td>0.009299</td>\n",
       "      <td>[[0.03232858511559456], [0.041896603634072056]...</td>\n",
       "      <td>0.009299</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>531</th>\n",
       "      <td>-0.200732</td>\n",
       "      <td>-0.009374</td>\n",
       "      <td>[[0.030658214427329314], [0.040183385866682644...</td>\n",
       "      <td>0.009374</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        alpha  residual                                                  y  \\\n",
       "536 -0.183790  0.000621  [[0.03156807432367406], [0.04112115231024713],...   \n",
       "535 -0.187179 -0.001451  [[0.03138253376730166], [0.04093080458090798],...   \n",
       "537 -0.180402  0.002731  [[0.03175543312517406], [0.04131290404695293],...   \n",
       "534 -0.190567 -0.003486  [[0.03119879457733433], [0.040741857883688865]...   \n",
       "538 -0.177013  0.004880  [[0.03194462688327459], [0.0415060621549289], ...   \n",
       "533 -0.193955 -0.005484  [[0.03101683973298032], [0.04055430867099873],...   \n",
       "539 -0.173625  0.007069  [[0.032135672115284425], [0.0417006283455375],...   \n",
       "532 -0.197344 -0.007446  [[0.030836652094924023], [0.0403681528598232],...   \n",
       "540 -0.170236  0.009299  [[0.03232858511559456], [0.041896603634072056]...   \n",
       "531 -0.200732 -0.009374  [[0.030658214427329314], [0.040183385866682644...   \n",
       "\n",
       "     abs_residual  \n",
       "536      0.000621  \n",
       "535      0.001451  \n",
       "537      0.002731  \n",
       "534      0.003486  \n",
       "538      0.004880  \n",
       "533      0.005484  \n",
       "539      0.007069  \n",
       "532      0.007446  \n",
       "540      0.009299  \n",
       "531      0.009374  "
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.sort_values(by='abs_residual', ascending=True).head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "min_pos 536\n",
      "alpha                                                    -0.18379\n",
      "residual                                              0.000620934\n",
      "y               [[0.03156807432367406], [0.04112115231024713],...\n",
      "abs_residual                                          0.000620934\n",
      "Name: 536, dtype: object\n"
     ]
    }
   ],
   "source": [
    "min_pos = df['abs_residual'].idxmin()\n",
    "print('min_pos', min_pos)\n",
    "print(df.iloc[min_pos])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x_opt(SDP) / y: [[1.07771194 1.02974908 1.00218376 0.96095702 0.96349426 1.00176999\n",
      "  0.96157006 1.24317896 0.97846892 1.03732259 0.99158415 0.96054483\n",
      "  0.96688204 1.03701798 1.0786745  0.99908685 1.0834897  0.93138703\n",
      "  0.93142046 1.06775708 1.14102158 0.9437951  1.09758494 1.10534923\n",
      "  1.09095222 1.03157305 1.14280023 1.10258652 1.10265224 1.18984706\n",
      "  1.10529734 1.21617174 0.93849277 1.04306492 1.16284959 1.05691426\n",
      "  1.15427166 1.2014746  1.07693653 0.93726841]]\n"
     ]
    }
   ],
   "source": [
    "print('x_opt(SDP) / y:', df.iloc[min_pos].y.T / x_opt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### what is strange?\n",
    "\n",
    "- `x_opt` (by SDP) and `y` (by linear solver) really close, even though the $\\alpha$ they used is different"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha from SDP:  0.198893609536022\n",
      "alpha from linear solver -0.18379012753640223\n",
      "**their signs are **different!**\n"
     ]
    }
   ],
   "source": [
    "print('alpha from SDP: ', alpha_opt)\n",
    "print('alpha from linear solver', df.iloc[min_pos].alpha)\n",
    "print('**their signs are', alpha_opt * df.iloc[min_pos].alpha > 0 and \"same\" or \"**different!**\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "y: [[-0.17229521 -0.03959887 -0.03812398 -0.03787031 -0.03697758 -0.03590856\n",
      "  -0.03303641 -0.03148591 -0.02250187 -0.01385053 -0.00431385 -0.00273143\n",
      "  -0.00141439 -0.00126847 -0.00116435 -0.00112809 -0.00105702 -0.00099749\n",
      "  -0.00085965 -0.00082685  0.00106466  0.00128402  0.00136589  0.00140513\n",
      "   0.00149516  0.00149832  0.00168793  0.00423548  0.0043123   0.00451254\n",
      "   0.02285036  0.02823298  0.03156807  0.0334434   0.03605723  0.03905723\n",
      "   0.04112115  0.04139267  0.050451    0.15647871]]\n",
      "node order:  [[13 14 10 18 17 11 12 15 19 16 30 39 35 21 36 37 33 31 38 34 25 22 20 32\n",
      "  24 26 29 23 27 28  8  9  0  5  2  6  1  3  7  4]]\n"
     ]
    }
   ],
   "source": [
    "print('y:', np.sort(df.iloc[min_pos].y.T))\n",
    "print('node order: ', np.argsort(df.iloc[min_pos].y.T))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# binary search on best $\\alpha$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "matrices loading done\n",
      "found lambda_1= (0.10625088207286182+0j)\n",
      "at iteration 5 (alpha=-16.52207)\n",
      "residual: sqrt(kappa) - y' D s=-0.16243300775250202\n",
      "at iteration 10 (alpha=-0.41338)\n",
      "residual: sqrt(kappa) - y' D s=-0.08410848378875846\n",
      "terminates after 14 iterations:\n",
      "  - alpha=-0.18604\n",
      "  - residual=-0.00076\n"
     ]
    }
   ],
   "source": [
    "from core import query_graph_using_sparse_linear_solver\n",
    "\n",
    "x_opt_sp, val_opt = query_graph_using_sparse_linear_solver(g, seeds, kappa, verbose=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "count    40.000000\n",
       "mean      1.045820\n",
       "std       0.081768\n",
       "min       0.928034\n",
       "25%       0.970553\n",
       "50%       1.041302\n",
       "75%       1.095096\n",
       "max       1.238488\n",
       "dtype: float64"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# difference w.r.t to the x_opt by SDP\n",
    "ratio = x_opt_sp / x_opt\n",
    "pd.Series(ratio).describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# running time comparison with SDP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "from core import query_graph_using_dense_matrix\n",
    "from helpers import sample_seeds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "random.seed(12345)\n",
    "np.random.seed(12345)\n",
    "\n",
    "size = 10\n",
    "k = 6\n",
    "g, true_comms, true_groupings = make_polarized_graphs(k, [(size, size) for i in range(k)])\n",
    "A = nx.adj_matrix(g, weight='sign')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def wrapper_sp():\n",
    "    seeds, target_comm = sample_seeds(true_comms, true_groupings)\n",
    "    query_graph_using_sparse_linear_solver(g, seeds, kappa, verbose=0)\n",
    "    \n",
    "def wrapper_sdp():\n",
    "    seeds, target_comm = sample_seeds(true_comms, true_groupings)\n",
    "    query_graph_using_dense_matrix(g, seeds, kappa, verbose=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "43.9 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "% timeit -n10 -r1 wrapper_sp()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARN: m less than n, problem likely degenerate\n",
      "WARN: m less than n, problem likely degenerate\n",
      "WARN: m less than n, problem likely degenerate\n",
      "WARN: m less than n, problem likely degenerate\n",
      "WARN: m less than n, problem likely degenerate\n",
      "WARN: m less than n, problem likely degenerate\n",
      "WARN: m less than n, problem likely degenerate\n",
      "WARN: m less than n, problem likely degenerate\n",
      "WARN: m less than n, problem likely degenerate\n",
      "WARN: m less than n, problem likely degenerate\n",
      "553 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "% timeit -n10 -r1 wrapper_sdp()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## conclusion\n",
    "\n",
    "- time difference is 10 fold on a small graph"
   ]
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}