{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "code is copied in hastings-2step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numba import jit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import cos, sin, pi\n",
    "from scipy.optimize import minimize\n",
    "from random import random\n",
    "\n",
    "@jit(nopython=True, nogil=True)\n",
    "def qaoa2(beta_2, gamma_2, beta_1, gamma_1, D):\n",
    "    \"\"\"Cut fraction with p=2 QAOA on D-regular girth > 5 graphs\"\"\"\n",
    "    c = cos(2*beta_2)\n",
    "    s = sin(2*beta_2)\n",
    "    m = cos(gamma_2)\n",
    "    n = sin(gamma_2)\n",
    "    r = cos(2*beta_1)\n",
    "    t = sin(2*beta_1)\n",
    "    G = cos(gamma_1)\n",
    "    H = sin(gamma_1)\n",
    "\n",
    "    A = -2*c*c*r*t*H*G**(D-1)\n",
    "    bpt1 = 0.5*s*c*(\n",
    "        (1 + r)*(-m*r*H - n*G)*((m*G - n*r*H)**(D-1))\n",
    "        + (1 - r)*(m*r*H - n*G)*((m*G + n*r*H)**(D-1))\n",
    "    )\n",
    "    bpt2 = 0.5*s*c*t*(\n",
    "        (m*t*(G**(D-1))*H + (1j)*n)*((m + (1j)*n*t*(G**(D-1))*H)**(D-1))\n",
    "        + (m*t*(G**(D-1))*H - (1j)*n)*((m - (1j)*n*t*(G**(D-1))*H)**(D-1))\n",
    "    )\n",
    "    Ea = s*s*t*H * 0.5 * \\\n",
    "        ((1+r)*(m*G - n*H*r)**(D-1) - (1-r)*(m*G + n*H*r)**(D-1))\n",
    "    Eb = (m + (1j)*n * G**(D-1) * H*t)**(D-1) + \\\n",
    "        (m - (1j)*n * G**(D-1) * H*t)**(D-1)\n",
    "    return 0.5 - 0.5 * (A + 2*(bpt1 + bpt2) + Ea*Eb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "qaoa2_d_deg = lambda i, d: -qaoa2(*i, d).real"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_best(d):\n",
    "    best = None\n",
    "    iter_times = 100 if d < 500 else 200\n",
    "    for _ in range(iter_times):\n",
    "        init_val = [random()*(i % 2 + 1)*pi for i in range(4)]\n",
    "        result = minimize(qaoa2_d_deg, init_val, args=(d), options={'fatol': 1e-20}, method='Nelder-Mead')\n",
    "        if not best or result.fun < best.fun:\n",
    "            best = result\n",
    "    return best"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = range(2, 1000 + 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import multiprocessing\n",
    "pool = multiprocessing.Pool(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 162 ms, sys: 44.8 ms, total: 206 ms\n",
      "Wall time: 5min 24s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "results = pool.map(get_best, ds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "out = {d: -r.fun - 0.5 for d, r in zip(ds, results)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "pool.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.4, 0.4745933547367115)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEGCAYAAABy53LJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dfZRcdZ3n8fen+iEJSRMeAg0mkQ5DQBCVSBtwcJ3WXRDdneDsMDvB8YFd5uTsjCw6O5zdeHYFBWfmyNkjrkd2xuyOux5Hja7OuG0mLuswlk9HMEEiJoFICJE04Sk8BDohJN393T/urerbXdXVVdUPtx8+r3PqVP1+93dvfX91Id/+3d99UERgZmZWr0LeAZiZ2ezixGFmZg1x4jAzs4Y4cZiZWUOcOMzMrCGteQcwHZYtWxZdXV1NrXvkyBEWL148uQHNcO7z/OA+zw8T6fP9999/KCLOGF0/LxJHV1cX27dvb2rdYrFIT0/P5AY0w7nP84P7PD9MpM+Sfl2t3oeqzMysIU4cZmbWECcOMzNriBOHmZk1xInDzMwa4sRhZmYNceKo4eVjJ3jx1aG8wzAzm1GcOGq483uPsPGHr+QdhpnZjOLEUYOUdwRmZjNPLolD0tWS9kjaK2ljjXbXSgpJ3Wn5DyTtyLyGJF0yfZGbmdm0Jw5JLcBdwLuBi4DrJF1UpV0HcBNwX6kuIr4SEZdExCXAB4D9EbFjymIF/HxEM7OR8hhxrAX2RsS+iDgObAauqdLuduAO4NgY27kO+NrUhJiQnDjMzEbL4yaHy4EDmXIfcFm2gaQ1wMqI2CLp5jG28/tUTzilbWwANgB0dnZSLBYbDvTAgeNERFPrzmb9/f3u8zzgPs8PU9HnPBJHtSnn8h/2kgrAncD1Y25Augw4GhE7x2oTEZuATQDd3d3RzN0hf3r0IXh8n++mOQ+4z/OD+zw58jhU1QeszJRXAAcz5Q7gYqAoaT9wOdBbmiBPrWeKD1MB4ENVZmYV8hhxbANWS1oFPEGSBN5XWhgRh4FlpbKkInBzRGxPywXg94C3T3WgcuYwM6sw7SOOiBgAbgTuBh4CvhERuyTdJmldHZt4O9AXEfumMk7w5LiZWTW5PAEwIrYCW0fV3TJG255R5SLJ4asp5+v/zMwq+crxGjziMDOr5MRRgxDhzGFmNoITRw2+V5WZWSUnjhp8yxEzs0pOHLV4yGFmVsGJo4ZS2ghPdJiZlTlx1FAacDhvmJkNc+KoQemYw3nDzGyYE0cNwyMOpw4zsxInjhrKcxy5RmFmNrM4cdTgOQ4zs0pOHDVIpTkOZw4zsxInDjMza4gTRx18qMrMbJgTRw2+cNzMrJITRw3l6zg84jAzK3PiqKF8VpUnx83Mypw4ahi+V1WuYZiZzShOHDUMjzjMzKzEiaOG4TkOpw4zs5JcEoekqyXtkbRX0sYa7a6VFJK6M3VvlPRTSbsk/VLSwqmLM3l32jAzG9Y63V8oqQW4C7gS6AO2SeqNiN2j2nUANwH3Zepagb8BPhARv5B0OnBiqmP2gMPMbFgeI461wN6I2BcRx4HNwDVV2t0O3AEcy9RdBTwYEb8AiIjnImJwqgKVhxxmZhWmfcQBLAcOZMp9wGXZBpLWACsjYoukmzOLzgdC0t3AGcDmiLij2pdI2gBsAOjs7KRYLDYc6KP7k8HMj378Y5a0z5+rAfv7+5v6vWYz93l+cJ8nRx6Jo9q/wOW/6SUVgDuB66u0awXeBrwFOArcI+n+iLinYoMRm4BNAN3d3dHT09NwoI/95DF4eDdXXHEFpy5ub3j92apYLNLM7zWbuc/zg/s8OfI4VNUHrMyUVwAHM+UO4GKgKGk/cDnQm06Q9wE/iIhDEXEU2Aq8eaoC9fM4zMwq5ZE4tgGrJa2S1A6sB3pLCyPicEQsi4iuiOgC7gXWRcR24G7gjZJOSifKfwvYXfkVk6N8W3XPjpuZlU174oiIAeBGkiTwEPCNiNgl6TZJ68ZZ9wXgMyTJZwfw84j4+6mK1XPjZmaV8pjjICK2khxmytbdMkbbnlHlvyE5JXfK+ZYjZmaVfOV4LX4CoJlZBSeOGsqnfzlvmJmVOXHU4DkOM7NKThw1+EFOZmaVnDhq8IOczMwqOXGYmVlDnDhq8Om4ZmaVnDhq8OS4mVklJ44a/ARAM7NKThy1lEYczhtmZmVOHDXMnydwmJnVz4mjhuG74+YciJnZDOLEUcPw8zicOczMSpw4apDnOMzMKjhx1ODTcc3MKjlx1ODTcc3MKjlx1OARh5lZJSeOOnjAYWY2zImjhtLpuB5zmJkNc+KowTc5NDOrlEvikHS1pD2S9kraWKPdtZJCUnda7pL0iqQd6euvpjbO5N15w8xsWOt0f6GkFuAu4EqgD9gmqTcido9q1wHcBNw3ahOPRsQl0xKrnwBoZlYhjxHHWmBvROyLiOPAZuCaKu1uB+4Ajk1ncFl+AqCZWaVpH3EAy4EDmXIfcFm2gaQ1wMqI2CLp5lHrr5L0APAS8J8j4kfVvkTSBmADQGdnJ8ViseFAdz01AMC2bdt5qmP+TAf19/c39XvNZu7z/OA+T448Eke1m86W/6SXVADuBK6v0u5J4LUR8ZykS4FvS3p9RLxUscGITcAmgO7u7ujp6Wk40GM7n4QdP+fSS7u56DUnN7z+bFUsFmnm95rN3Of5wX2eHHn8Gd0HrMyUVwAHM+UO4GKgKGk/cDnQK6k7Il6NiOcAIuJ+4FHg/KkLNZ3j8KEqM7OyPBLHNmC1pFWS2oH1QG9pYUQcjohlEdEVEV3AvcC6iNgu6Yx0ch1J5wKrgX1TFaj8QA4zswrTfqgqIgYk3QjcDbQAX4yIXZJuA7ZHRG+N1d8O3CZpABgE/m1EPD/1MU/1N5iZzR55zHEQEVuBraPqbhmjbU/m87eAb01pcBkecJiZVZo/pwo1wU8ANDOr5MRRg58AaGZWyYmjBj8B0MyskhNHDb5XlZlZpXETh6TB6QhkJvITAM3MKtUz4pi/Jxd5xGFmVqGexDFv/9308zjMzCo1PcchaYWkN0haPJkBzSR+AqCZWaWGE0f6MKWfkzwn49vAM5K+I2kK7xmVD484zMwqNTPi+DTwhYhYHhG/ASwFvgN8V9LqSY0uZz6rysysUjOJ4/yI+EKpEBED6S3M/wioetuQ2cpPADQzq9RM4qj6z2hE/D/gwomFM7MMXwDozGFmVtJM4jhL0g2SLpO0ZNSyOfUvrKfGzcwqNXN33E8Aa4APAhdLegnYBewEzpq80GYA33LEzKxCw4kjnc8ok7QCeCPwBuAHkxTXjFCQrxw3MxutnsQhSW+IiF9WWxgRfSSPg91abfls1lJIEseQ84aZWdm4iSMiCpLOlfQv06rngJ9ExMDUhpa/NG8w5BGHmVlZXYeqImIf6bO9JZ0CvFtSK3CcJIm8OHUh5qd05figE4eZWVkzcxwvklzwR5o8rpB0GrAAuDci9k9qhDlq8RyHmVmFiT5zvAt4C8n1G88CvwL2T3CbM0ZpcnxwKOdAzMxmkInc5PCTwO8CWyLihojYGBE/b2D9qyXtkbRX0sYa7a6VFJK6R9W/VlK/pJub7cN4Cumv4zkOM7NhExlx3AG8EhEN/z0uqQW4C7iS5IysbZJ6I2L3qHYdwE0kN1Qc7U7guw1H3QCfjmtmVqnpEUdEHGkmaaTWAnsjYl9EHAc2A9dUaXc7SYI6lq2U9F6SyfpdTX5/XXyoysys0kTnOJq1HDiQKfcBl2UbSFoDrIyILdnDUenzP/4jyWhlzMNUkjYAGwA6OzspFosNB3mwP8kYO3ftYvHzexpef7bq7+9v6veazdzn+cF9nhx5JY5qj6MtHw+SVCA5FHV9lXafBO6MiP7hBy1V2VhyhfsmgO7u7ujp6Wk4yEef7Ycf/4DXXXghPZcsb3j92apYLNLM7zWbuc/zg/s8OfJKHH3Aykx5BXAwU+4ALgaKaXI4C+iVtI5kZHKtpDuAU4AhScci4vOTHWTpdFxPjpuZDWtqjkPS+9L39U1+7zZgtaRVktqB9UBvaWFEHI6IZRHRFRFdwL3AuojYHhH/JFP/WeDPpyJpgOc4zMyqaXZyfLmkf0UyUmhYeruSG4G7gYeAb0TELkm3paOKGcGn45qZVWr4UJWkW4GFwF8Ad0i6JSJua3Q7EbGVUTdGjIiqTxCMiJ4x6j/R6Pc2wqfjmplVanjEERGfBJ4H3g8830zSmC18qMrMrFKzh6oORsRm4InJDGam8aEqM7NKTSWOiPhK+v61yQ1nZin4rCozswpNXzk+H5RPx/WTnMzMypw4aijPcThvmJmVTegCQEmnAqtJzrICICJ+ONGgZgqladVnVZmZDWs6cUj6Q+AjJNdy7AAuB34KvHNyQsufrxw3M6s0kUNVHyF5iNOvI+IdwBqShznNGT4d18ys0kQSx7GIOAYgaUFEPAxcMDlhzQw+HdfMrNJE5jj6JJ0CfBv4nqQXGHmjwlmv4LOqzMwqNJ04IuJ30o+fkPR9YClT/ES+6TY8x5FzIGZmM8hEnjn+6dLniPhBRPQCn5qUqGaI0uM+fKjKzGzYROY4rqxS9+4JbG/GkYRw4jAzy2rm7rh/BPwxcK6kBzOLOoCfTFZgM0VBThxmZlnNzHF8lWQu4y+AjZn6lyPi+UmJagYRPh3XzCyr4cQREYeBw8B1kx/OzFOQrxw3M8tq5lDVy0DpX9J0+phIP0dEnDxJsc0IEgz6tCozs7JmRhwdUxHITFUQDHrEYWZWNpHTcSXp/ZI+npZXSlo7eaHNDC2CAd8e18ysbCKn4/434K3A+9JyP3DXhCOaYVoKYmDIs+NmZiUTSRyXRcSHgWMAEfEC0F7vypKulrRH0l5JG2u0u1ZSSOpOy2sl7Uhfv5D0O2OtOxk84jAzG2ki96o6IamFdKJc0hlAXX+ap+vdRXIRYR+wTVJvROwe1a4DuAm4L1O9E+iOiAFJZwO/kPSdiBiYQF/G1CIY8OS4mVnZREYcnwP+DuiU9GfAj4E/r3PdtcDeiNgXEceBzcA1VdrdDtxBOqoBiIijmSSxkOEzvKZEi+CEL+QwMyubyE0OvyLpANADPAe8NyIeqnP15cCBTLkPuCzbQNIaYGVEbJF086hllwFfBM4BPlBttCFpA7ABoLOzk2KxWGdoow3x1NPPTGD92ae/v39e9Rfc5/nCfZ4czVzHIeBW4EaSazcKwABwOnBbvZupUlceOUgqAHcC11dbOSLuA14v6ULgS5K+W3o2SKbNJmATQHd3d/T09NQZ2khtP/kup5y2jJ6e7qbWn42KxSLN/l6zlfs8P7jPk6OZQ1UfBa4A3hIRp0fEqSSjhSsk/Umd2+gDVmbKKxj5LI8O4GKgKGk/yWNpe0sT5CXpCOdI2nZKtAqfVWVmltFM4vggcF1EPFaqiIh9wPvTZfXYBqyWtEpSO7Ae6M1s73BELIuIrojoAu4F1kXE9nSdVgBJ55A8dXB/E/2oS8FXjpuZjdDMHEdbRBwaXRkRz0pqq2cD6RlRNwJ3Ay3AFyNil6TbgO3psz3G8jZgo6QTJGdx/XG1eCZLS8GT42ZmWc0kjuNNLhshIrYCW0fV3TJG257M5y8DX673eybK13GYmY3UTOJ4k6SXqtSL5PTYOaVF8nUcZmYZzdzksGUqApmpWgqeHDczy5rIBYDzgg9VmZmN5MQxjoJvOWJmNoITxzhaCzDgs6rMzMqcOMbRVhCvDjhxmJmVOHGMo60Fjp0YzDsMM7MZw4ljHO0FOHbCIw4zsxInjnG0tYhjA4OEnztuZgY4cYyrvQARcNwT5GZmgBPHuNpbkjvA+3CVmVnCiWMcbekv9KonyM3MACeOcZUSh0ccZmYJJ45xlA9VDXjEYWYGThzjak9v6ehrOczMEk4c42greHLczCzLiWMcHnGYmY3kxDGO9vLkuBOHmRk4cYyrrTw57kNVZmbgxDEujzjMzEbKJXFIulrSHkl7JW2s0e5aSSGpOy1fKel+Sb9M39851bG2pXMcvgDQzCzR8DPHJ0pSC3AXcCXQB2yT1BsRu0e16wBuAu7LVB8CfjsiDkq6GLgbWD6V8bb7rCozsxHyGHGsBfZGxL6IOA5sBq6p0u524A7gWKkiIh6IiINpcRewUNKCqQx2QTriOHJ8YCq/xsxs1pj2EQfJCOFAptwHXJZtIGkNsDIitki6eYzt/C7wQES8Wm2hpA3ABoDOzk6KxWJTwb5y9AgLW8SuRx6j2Hpw/BXmgP7+/qZ/r9nKfZ4f3OfJkUfiUJW68sMuJBWAO4Hrx9yA9Hrg08BVY7WJiE3AJoDu7u7o6elpKthiscjpHUN0nHYaPT2XNLWN2aZYLNLs7zVbuc/zg/s8OfI4VNUHrMyUVwDZP+U7gIuBoqT9wOVAb2aCfAXwd8AHI+LR6Qh46aI2Dr9yYjq+ysxsxssjcWwDVktaJakdWA/0lhZGxOGIWBYRXRHRBdwLrIuI7ZJOAf4e+FhE/GS6Aj7lJCcOM7OSaU8cETEA3EhyRtRDwDciYpek2yStG2f1G4HzgI9L2pG+zpzikD3iMDPLyGOOg4jYCmwdVXfLGG17Mp8/BXxqSoOrYumiNl504jAzA3zleF2W+lCVmVmZE0cdli5q4/jAkG87YmaGE0ddTl/cDsCzL1e9ZMTMbF5x4qjDWUsXAfD0S8fGaWlmNvc5cdTh7KULAXjysBOHmZkTRx2GE8crOUdiZpY/J446dCxsY8mCVo84zMxw4qjbWUsXcvBFjzjMzJw46nTussXse/ZI3mGYmeXOiaNO5525hMcOHeHEoB/oZGbzmxNHnVZ3LmFgKPj1cx51mNn85sRRp/M7OwDY/eTLOUdiZpYvJ446XdDZwaK2Fn7+6xfyDsXMLFdOHHVqbSnwppVLud+Jw8zmOSeOBlx+7unsPHjY96wys3nNiaMBV110FhHwDw89nXcoZma5ceJowIVnd/Da007iuzufyjsUM7PcOHE0QBLvveQ1/OiRZ3nskE/LNbP5yYmjQe9/6zm0FQps+uG+vEMxM8uFE0eDzuxYyPq1K/n6tsfZffClvMMxM5t2uSQOSVdL2iNpr6SNNdpdKykkdafl0yV9X1K/pM9PX8Qj/emVF3DqSe38+2/s4MirA3mFYWaWi2lPHJJagLuAdwMXAddJuqhKuw7gJuC+TPUx4OPAzdMQ6piWntTGZ37/En719Mt89Os7OD7g+1eZ2fyRx4hjLbA3IvZFxHFgM3BNlXa3A3eQJAsAIuJIRPw4W5eX3zr/DG797dfzvd1Pc8OXtvH8keN5h2RmNi3ySBzLgQOZcl9aVyZpDbAyIrZMZ2CN+tBvdnHHtW/k3n3PcdWdP6T3FwcZGoq8wzIzm1KKmN5/6CT9HvCuiPjDtPwBYG1E/Lu0XAD+Ebg+IvZLKgI3R8T2zDauB7oj4sYa37MB2ADQ2dl56ebNm5uKt7+/nyVLltRsc+DlIf77g6/y+MtDrOwocNU5rXSf1cqiVjX1nXmrp89zjfs8P7jPjXnHO95xf0R0j65vnXBUjesDVmbKK4CDmXIHcDFQlARwFtAraV02eYwnIjYBmwC6u7ujp6enqWCLxSL1rPu+fx5sefAgn7vnEf565xG+umeQngvO4O3nn8HbzlvGilMXkfZnxqu3z3OJ+zw/uM+TI4/EsQ1YLWkV8ASwHnhfaWFEHAaWlcrVRhwzUUtBXHPJcta96TX8/PEX+Ob9T/D9h58pX2W+bEk7b1i+lIuXL+XcMxbTdfpiVi1bzCknteccuZlZY6Y9cUTEgKQbgbuBFuCLEbFL0m3A9ojorbW+pP3AyUC7pPcCV0XE7qmOu16SuPSc07j0nNOICB55pp979z3Hg32H+WXfYX7wq2fJToMsXdRG58kL6Dx5IWd2LKTz5AWc2bGAU05qZ+miNk5e1MYpJ7WxdFHyamvxpTdmlq88RhxExFZg66i6W8Zo2zOq3DVlgU0ySZzf2VF+CBTAqwODHHj+KI8dOsr+Q0d4/PmjPPPyMZ5+6VUefeYQz7z8KgM1JthPam/hpPbW9L2FRaX3tpF1i9paaGsp0N5aYEFr8t7WUqC9ZfjziPrW0jLx1JEhDjx/lJaCyq/W8nuBQoHkXcyaw29mNnlySRzz2YLWFs47s4PzzuyounxoKHjh6HFefOUEh9PXS+n7i0eTz0eOD/LK8QGOHh/klRODHD0+yPNHXhmuS+trJaBx/ej7dTVrzSSX4QRTGFHfWhASFCQKynwukJZFobycUeXsulXaF0rtM8thjDYgku2V2gBI8ETfq/yof3danywTQLpOqV15eWY7lL9z9PZHJtaqy9MyDH9ntfjG/Z4k0JHrlWNIt8/IWHYfHODwjidG7M8R8WbaDteN3MbotiPr62lb+X3Z9mN99xgfR/S92nfvOjRI6yOHxo6p4b5qnOVQbcnwPqncVq1tZP9brCW7/OiJyT8ByoljhikUxOlLFnD6kgUT3tbQUHB8cCh5DQxxIn0/PjBcl9QHxwcH0/pg567dnH/B6xgcGmJgKBgaCgaGgsHM+/DnpM3gYDAYmfrB5H0oghODQwQQEQwNwVAEQ5GW089DEUSUliV1g0ND5WUj2g9Vto9Mu1I5MuXs+kksSTkAAk4MDNDy5OPDy0i3UfoxR9VFZjuz2oM78o5g+m2/b/w2c8hH3ryA90zyNp045rBCQSwstLCwraWh9U5+4Vf0XLpiiqKamSZ65knE2EmllHCgMimNaFslYWWXR9KganIrff9wPKO+d1SsAD/72c94y9q1I9bJtKqoG7mNbH1U1I+3fLRq7cf+7hijvrJ29DYeeOAB1qxZM+q3ihrbGhl/tY/19LVa/FV+7trbqLIvs9ur9l0ArxyY/ClgJw6zSSCp6iGamezxxQV+44z5dU3Dkf0tvKXrtLzDmFbFZx+e9G36FB0zM2uIE4eZmTXEicPMzBrixGFmZg1x4jAzs4Y4cZiZWUOcOMzMrCFOHGZm1hAnDjMza4gTh5mZNcSJw8zMGuLEYWZmDXHiMDOzhjhxmJlZQ5w4zMysIU4cZmbWECcOMzNrSG6JQ9LVkvZI2itpY41210oKSd2Zuo+l6+2R9K7pidjMzCCnR8dKagHuAq4E+oBtknojYveodh3ATcB9mbqLgPXA64HXAP8g6fyIGJyu+M3M5rO8Rhxrgb0RsS8ijgObgWuqtLsduAM4lqm7BtgcEa9GxGPA3nR7ZmY2DXIZcQDLgQOZch9wWbaBpDXAyojYIunmUeveO2rd5aO/QNIGYENa7Je0p8lYlwGHmlx3tnKf5wf3eX6YSJ/PqVaZV+JQlbooL5QKwJ3A9Y2uW66I2ARsajK+4S+TtkdE9/gt5w73eX5wn+eHqehzXomjD1iZKa8ADmbKHcDFQFESwFlAr6R1daxrZmZTKK85jm3AakmrJLWTTHb3lhZGxOGIWBYRXRHRRXJoal1EbE/brZe0QNIqYDXws+nvgpnZ/JTLiCMiBiTdCNwNtABfjIhdkm4DtkdEb411d0n6BrAbGAA+PMVnVE34cNcs5D7PD+7z/DDpfVZExfSAmZnZmHzluJmZNcSJw8zMGuLEMYZ6b4ky20haKen7kh6StEvSR9L60yR9T9Ij6fupab0kfS79HR6U9OZ8e9A8SS2SHpC0JS2vknRf2uevpydqkJ548fW0z/dJ6soz7mZJOkXSNyU9nO7vt871/SzpT9L/rndK+pqkhXNtP0v6oqRnJO3M1DW8XyV9KG3/iKQPNRKDE0cVmVuivBu4CLguvdXJXDAA/GlEXAhcDnw47dtG4J6IWA3ck5Yh+Q1Wp68NwF9Of8iT5iPAQ5nyp4E70z6/ANyQ1t8AvBAR55FcT/TpaY1y8vxX4P9GxOuAN5H0fc7uZ0nLSW5R1B0RF5OceLOeubef/xdw9ai6hvarpNOAW0kuvF4L3FpKNnWJCL9GvYC3Andnyh8DPpZ3XFPU1/9Dcs+wPcDZad3ZwJ708xeA6zLty+1m04vkep97gHcCW0guJD0EtI7e5yRn+701/dyatlPefWiwvycDj42Oey7vZ4bvSHFaut+2AO+ai/sZ6AJ2NrtfgeuAL2TqR7Qb7+URR3XVbolScVuT2S4dmq8huYlkZ0Q8CZC+n5k2myu/xWeB/wAMpeXTgRcjYiAtZ/tV7nO6/HDafjY5F3gW+J/p4bn/IWkxc3g/R8QTwH8BHgeeJNlv9zO393NJo/t1QvvbiaO6um5rMptJWgJ8C/hoRLxUq2mVuln1W0j6F8AzEXF/trpK06hj2WzRCrwZ+MuIWAMcYfjwRTWzvs/poZZrgFUkd85eTHKoZrS5tJ/HM1YfJ9R3J47q5vRtTSS1kSSNr0TE36bVT0s6O11+NvBMWj8XfosrgHWS9pPcifmdJCOQUySVLoLN9qvc53T5UuD56Qx4EvQBfRFReiTBN0kSyVzez/8MeCwino2IE8DfAr/J3N7PJY3u1wntbyeO6mreEmU2kyTgr4GHIuIzmUW9QOnMig+RzH2U6j+Ynp1xOXC4NCSeLSLiYxGxIpLb16wH/jEi/gD4PnBt2mx0n0u/xbVp+1n1l2hEPAUckHRBWvVPSe62MGf3M8khqsslnZT+d17q85zdzxmN7te7gasknZqO1K5K6+qT9yTPTH0B7wF+BTwK/Ke845nEfr2NZEj6ILAjfb2H5NjuPcAj6ftpaXuRnGH2KPBLkjNWcu/HBPrfA2xJP59Lcp+zvcD/Bhak9QvT8t50+bl5x91kXy8Btqf7+tvAqXN9PwOfBB4GdgJfBhbMtf0MfI1kDucEycjhhmb2K/Bv0r7vBf51IzH4liNmZtYQH6oyM7OGOHGYmVlDnDjMzKwhThxmZtYQJw4zM2tIXs8cN5vXJA2SnB7ZRnLjyS8Bn42IoZorms0AThxm+XglIi4BkHQm8FWSK5dvzTUqszr4Og6zHEjqj4glmfK5JHcsWBb+n9JmOM9xmM0AEbGP5P/HM8dra5Y3Jw6zmaPaHUvNZhwnDrMZID1UNcjwXU3NZiwnDrOcSToD+Cvg857fsNnAk+NmOahyOu6Xgc/4dFybDZw4zMysIT5UZWZmDXHiMDOzhjhxmJlZQ5w4zMysIU4cZmbWEK/y6+wAAAASSURBVCcOMzNriBOHmZk15P8DgAUL+r4fGrYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.plot([k for k in sorted(out)], [out[k]*k**0.5 for k in sorted(out)])\n",
    "plt.grid()\n",
    "plt.xlabel(\"D\")\n",
    "plt.ylabel(\"Delta * $\\sqrt {D} $ \")\n",
    "plt.ylim(0.4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np \n",
    "sum (  np.array([out[k]*k**0.5 for k in sorted(out)]) < 0.405 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.40762859, 0.40762851, 0.40762843, 0.40762834, 0.40762826,\n",
       "       0.40762817, 0.40762809, 0.40762801, 0.40762792, 0.40762784])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array([out[k]*k**0.5 for k in sorted(out)])[-10:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}