Survivors from strain pools.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1928b4d3",
   "metadata": {},
   "source": [
    "### Fine-scale diversity due to satellite niches in ecosystems with boom-and-bust cycles\n",
    "\n",
    "It has long been known that more than one strains with similar functions can coexist in both natural and lab environments, yet the mechanism behind remains as a puzzle. Here, we propose a model to study this coexitence, where strains that obey Monod's growth law undergo boom and bust cycles. \n",
    "\n",
    "Here, we take strain pools of different species (defined by having a distinct function, i.e. portion of protein it use to consume different resources) and investigate the survivors from such pools. We see that there can be multiple strains from the same species coexisting in the final state. We say these strains are in satellite niches of their species. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "52da93ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.integrate import solve_ivp\n",
    "import copy\n",
    "import pickle\n",
    "from math import *\n",
    "from matplotlib.colors import ListedColormap, LinearSegmentedColormap\n",
    "import seaborn as sns\n",
    "import scipy.stats as stats\n",
    "import time\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "cdec2ede",
   "metadata": {},
   "outputs": [],
   "source": [
    "# g: maximum growth rate\n",
    "# K: Monod constant\n",
    "def OneStep(t, u, p):\n",
    "    N, R, gmax, K, x = p\n",
    "    Ni, Rj = u[:N], u[N:]\n",
    "    xmat = np.hstack([x.T, 1-x.T])\n",
    "    cmat = np.vstack([Rj for i in range(N)])\n",
    "    gmat = xmat*(gmax*(cmat/(K+cmat)))\n",
    "#     print(gmat)\n",
    "    return np.hstack([np.sum(gmat, axis=1)*Ni, -Ni@gmat])\n",
    "\n",
    "\n",
    "def Dilute(p, p_dilute, u_end):\n",
    "    N, R, gmax, K, x = p\n",
    "    c0, d, t, n, fluc_thres, elim_thres = p_dilute\n",
    "    output = u_end\n",
    "    u_add = np.zeros(N+R)\n",
    "    u_add[-R:] = np.array(c0)\n",
    "    stop, count = 1, 1\n",
    "    utransf = copy.deepcopy(u_end)\n",
    "    while stop>0 and count<n: \n",
    "        u = copy.deepcopy(utransf)\n",
    "        u = u/d + u_add\n",
    "        N_alive = u[:N]>0\n",
    "        N_temp = sum(N_alive)\n",
    "        gmax_temp = gmax[N_alive]\n",
    "        K_temp = K[N_alive]\n",
    "        u_temp = np.hstack([u[:N][N_alive], u[N:]])\n",
    "        x_temp = x[0:,N_alive]\n",
    "        p_temp = N_temp, R, gmax_temp, K_temp, x_temp\n",
    "        sol = solve_ivp(OneStep, t_span=[0, t], y0=u_temp, method=\"LSODA\", args=[p_temp], rtol=1e-6)\n",
    "        utransf_temp = sol.y[:, -1] * (sol.y[:, -1]>elim_thres*d)\n",
    "        u[:N][N_alive] = utransf_temp[:N_temp]\n",
    "        u[-R:] = utransf_temp[-R:]\n",
    "        utransf = u\n",
    "        output = np.vstack([output, utransf])\n",
    "        stop = sum(abs((output[-1, :-R]-output[-2, :-R])) > fluc_thres*output[-1, :-R])\n",
    "        count+=1\n",
    "#     print(\"c\", count)\n",
    "    return output"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59d1d46b",
   "metadata": {},
   "source": [
    "### We start from a single resource and single species\n",
    "\n",
    "There could be 2 or even 3 satellite niches for a single species. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "e9f2fa71",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of coexistence community 2\n",
      "coexistence community [[1.0129039322408333, 0.16894171079681455], [1.0921105875522974, 0.48642814731585643], [1.1799589005680498, 0.9124767440293353]]\n",
      "index of that community 695\n",
      "final profile [3.62131667 0.01137843 7.47841602 0.        ]\n",
      "coexistence community [[0.9674590707974304, 0.45056709952411456], [1.0151638305652246, 0.693363919124809], [1.158425502180989, 1.4575872437223094]]\n",
      "index of that community 1319\n",
      "final profile [ 2.48084328  0.          8.63026764 -0.        ]\n"
     ]
    }
   ],
   "source": [
    "# randomly generate 2000 pools, each one has 100 strains\n",
    "N=100\n",
    "Nsamples = 2000\n",
    "all_paretobugs = []\n",
    "thres=1e-2\n",
    "all_params = []\n",
    "c0, D, t, n, fluc_thres, elim_thres = 10.0, 10.0, 24.0, 10000, 1e-5, 2.0e-7\n",
    "for i in range(Nsamples):\n",
    "    gs_all = np.exp(np.random.normal(log(1), 0.05, size=N))\n",
    "    gs_all *= (gs_all>0)\n",
    "    Ks_all = 10**np.random.normal(log10(4), 0.5, size=N)\n",
    "    x = 1\n",
    "    paretobugs = findCoexistPareto(findPareto(gs_all, Ks_all), c0, D)\n",
    "    all_paretobugs.append(paretobugs)\n",
    "    all_params.append([gs_all, Ks_all])\n",
    "    \n",
    "print(\"number of coexistence community\", len([i for i in all_paretobugs if len(i)>2]))\n",
    "    \n",
    "for idx, i in enumerate(all_paretobugs):\n",
    "    if(len(i)>2):\n",
    "        print(\"coexistence community\", i)\n",
    "        paretobugs = i\n",
    "        gs_pareto = np.array([i[0] for i in paretobugs])\n",
    "        Ks_pareto = np.array([i[1] for i in paretobugs])\n",
    "        p = len(gs_pareto), gs_pareto, Ks_pareto, 0\n",
    "        p_dilute = c0, D, t, n, fluc_thres, elim_thres\n",
    "        u_end = np.array([1.0 for i in gs_pareto]+[0.0])\n",
    "        profile = Dilute(p, p_dilute, u_end, x)\n",
    "        all_paretobugs[idx] = [j for jdx, j in enumerate(i) if profile[-1][jdx]>0]\n",
    "        print(\"index of that community\", idx)\n",
    "        print(\"final profile\", profile[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "7bd96a74",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEVCAYAAADdFfNTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAqVUlEQVR4nO3de9icdX3n8fcngQRCSUJIIFkbBEIwgu16iFW2bEFrTcEqlILpqtvqYrPqZbN1abtXPYGHbldalVIr8nisrdUqKoqFDWo5lDYgYfGAEAlIBBQwIQcwJyD57h+/GZ/JZOZ55n7mPs3M53Vdc80zv/swv5k7ub/zOysiMDMzy9u0qjNgZmbDyQHGzMwK4QBjZmaFcIAxM7NCOMCYmVkhHGDMzKwQpQcYSSdIulzSdyTtlXR9D8c8X9InJd0jaaekH0i6UNIhHfY9S9L3JO2WdKeklYV8EDMzm9BBFbznycCZwM3AjB6PWQksAd4HbAB+GXhP4/l3mjtJOhX4IvBhYHXjfT4raWtEXDvRG8yfPz+OPfbYTB/EzGzU3XbbbZsjYkGnbSp7oKWkaRGxr/H3FcD8iDh9kmMWRMSmtrRVwOXAsRHxo0baGuDgiHhxy35XA7Mj4tSJ3mP58uWxbt26qXwkM7ORJem2iFjeaVvpVWTN4JLxmE0dkm9vPB8FIGkm8CLg8237fQ44RdKcrO9rZmZTN8iN/P8J2Af8oPF6CXAwsL5tv7tIn/PE8rJmZmYDGWAkLQTeBvx9RDzWSD6i8bytbfetbdtbz7NK0jpJ6zZt6lRIMjOzqRq4ACNpBqka7GfAWzrs0t6opC7pRMRYRCyPiOULFnRsozIzsymqohfZlEkS8GlST7RfjYitLZubf89tO6z5eluReTMzs/0NVIABPgicBfxGRLS3tdwLPAksA25oSV9Gaqu5u5QcWrW2bIENG2DbNpg7F5YuhXnzqs6V2UgamCoySX8G/CHwmoi4qX17ROwBrgPOa9u0ElgbEduLz6VVassWuOUW2LMHjjwyPd9yS0o3s9KVXoKRNIs0ABLgacBsSec2Xl8dETsl3QPcEBHnN455FfC/gU8BP5b0wpZT3tvSjfk9wPWSLgGubLzPmcBvFveJrDY2bIDDDksPGH/esAFe8ILq8mU2oqqoIjsK+EJbWvP1ccBGUr6mt2x/aeP5tY1Hq9eRAg8RcVMjWL0XeCNwH/CqyUbx25DYti2VXFrNmgWPPlpJdsxGXekBJiI2Mt6zq9s+x7a9fi0HBpZux15JKr3YqJk7F3buHC+5QHo9d25VOTIbaQPTBmM2qaVLYceO9IgY/3vp0qpzZjaSHGBseMybl9paZs5M1WIzZ6bX7kVmVgkHGCvX2BgsWgTTpqXnsbF8z98MMitWOLiYVWzQxsHYIBsbg9WrU/dhgIcfTq8BVq2qLl9mVgiXYKw8F144Hlya9uxJ6WY2dBxgrDyPPJIt3cwGmgOMlefoo7Olm9lAc4Cx8rzrXalnV6uZM1O6mQ0dBxgrz6pVcOmlsHAhSOn50kvdwG82pNyLzMq1apUDitmIcAnGzMwK4QBjZmaFcIAxM7NCOMCYmVkhHGDMzKwQDjBmZlYIBxgzMyuEA4yZmRXCAcbMzArhAGNmZoUoPcBIOkHS5ZK+I2mvpOt7OGaGpL+U9K+SdkmKLvt9SlJ0eCzL/YOYmdmEqpiL7GTgTOBmYEaPx8wCXg98C/h34MUT7LseeF1b2sZsWTQzs35VEWCuioivAEi6Apg/2QERsU3SvIgISW9m4gCzIyJuzimvZmY2RaVXkUXEvike17FazMzM6mkYG/lPkvSYpD2SbpJ0WtUZMjMbRcMWYG4HLgBeDrwamA58XdKvdNpZ0ipJ6ySt27RpU4nZNDMbfkO14FhE/HXra0n/DNwJvBU4u8P+Y8AYwPLly10FZ2aWo2ErwewnInYBVwPPrTovZmajZqgDTAuXTszMSjbUAUbSocAZwG1V58XMbNSU3gYjaRZpoCXA04DZks5tvL46InZKuge4ISLObznuDOAw4NmN181jbo2IH0maA3wN+AfgHtL4mrc03uOVxX4qMzNrV0Uj/1HAF9rSmq+PI426P4jUA6zVZcDTOxzzOuBTwB5gE/D2xnvsBtYCp0XEunyybmZmvSo9wETERkCT7HNsL2lt23cD5/SRNTMzy9FQt8GYmVl1HGDMzKwQDjBmZlYIBxgzMyuEA4yZmRXCAcbMzArhAGNmZoVwgDEzs0I4wJiZWSEcYMzMrBAOMGZmVggHGDMzK4QDjJmZFcIBxszMCuEAY2ZmhahiwTGzamzZAhs2wLZtMHcuLF0K8+ZVnSuzoeUAY6Nhyxa45RY47DA48kjYuTO9fsELBivIOEjaAHEVmeWveTNfsyY9b9lSdY7STfmww9JDGv97w4aqc9a75ve6Z08Kknv21Of7NevAAcbyVdeb4LZtMGvW/mmzZqX0QTEMQdJGSukBRtIJki6X9B1JeyVd38MxMyT9paR/lbRLUkyw71mSvidpt6Q7Ja3M9QPYxOp6E5w7N1WLtdq5M6UPimEIkjZSqijBnAycCdzdePRiFvB6YCfw7912knQq8EXgOuAM4J+Bz0p6aT8ZtgzqehNcuhR27EiPiPG/ly6tNl9ZDEOQtJFSRYC5KiIWR8R5wPd7OSAitgHzImIF8OUJdn0HcGNErI6I6yLiT4D/C7yz30xbj+p6E5w3LzXoz5wJjz6angetgX8YgqSNlNJ7kUXEvike17VaDEDSTOBFwOq2TZ8DPilpTkRsn8p7WwZLl6Y2F0gll507003wWc+qNl8wHmSKVlRPr2b+N2xIQXLu3PS9DlKQtJEyTN2UlwAHA+vb0u8ildROBG4tO1MjZ9RvgkV3hy4rSJrlYJgCzBGN521t6Vvbtv+cpFXAKoBjjjmmsIyNnFG+CbZ2coDx5w0bqv1OPH7GKjBMAaapvSpNXdKJiDFgDGD58uUTVsHZECjjJrttWyq5tJo1K5Xmssgzr8MyyNQGzjCNg2mWVOa2pTdfbysrIwNlbAwWLYJp09Lz2FjVOSpGWeNz8ujkMJW8TjS4ta5dx23oDVOAuRd4EljWlr4M2EfvXaJHx9gYrF4NDz+ceiU9/HB6PYxBpqybbB49vbLmdbKAVNeu4zb0hibARMQe0viX89o2rQTWugdZBxdemG5GrfbsSenDpqybbB7dobPmdbKAVNeu4zb0Sm+DkTSLNNAS4GnAbEnnNl5fHRE7Jd0D3BAR57ccdwZwGPDsxuvmMbdGxI8af78HuF7SJcCVjfc5E/jNwj7QIHvkkWzpg6x5k202ukNxN9l+Ozlkzetk7T517jpuQ62KEsxRwBcajxcCJ7W8Pqqxz0HA9LbjLmvs0ww6zWNe1NwhIm4CzgVeAqwBXgG8KiKuLeKDDLyjj86WPsgGaZBi1rxOVkIZhkGmNpA0yfjFkbF8+fJYt25d1dkoV7MNprWabOZMuPRSWLWqunwVZZC66mbJa2svsdYSShVBZJC+Y8uFpNsiYnnHbQ4wyUgGGEhB5sILU7XY0UfDu941nMFl2NXhxl6nQGelcYDpwcgGmMnU4cZlg6HZk6217WjHjvEqORtKEwWYYRxoaXnxAL39OdhOLK9BpjY0hqabshXAA/TG1XUhtTpxd2hr4wBj3XmA3jgH28kNUk89K4UDjHXnX6TjHGwn5+7Q1sZtMNadB+iNK3Og5iAb5Zm07QAuwVh3/kU6ztU/Zpm5BGMT8y/SZNQXUjObAgcYs1452Jpl4gBjVlced2MDzgHGrI4GaZCrA6F14UZ+szoalHE3HoBqE3CAMaujQRl3MyiB0CrRc4CR9N8m2DZD0l/mkyUzG5hBroMSCK0SWUowl0u6StJ+q1FJWg58G+gagMwso6mMu2lWV61ZU1411aAEQqtElgDzq8AJwPcl/a6kgyT9ObAW2Aj8UgH5M5u6Km64eck6yLWqthAPQLUJZFoPRtJM4C+APwS2ADOBP46IjxWTvfJ4PZghM2qLX1W5FktdepHVJR8jJs/1YJ4kBZZ9wFzgR8B3+sqdWRFaG59h/HnDhuEcLFnlWix1GIA6SN26R0iWRv5lpOqwPwX+CDgGuAv4N0nvldRTsJJ0gqTLJX1H0l5J1/d43BxJn5S0VdJ2SZ+RdGTbPp+SFB0ey3r9nDYkRq3xedTbQtybrZaylGBuB24FnhMR9zbSzpL0e8AlwG8Bz+7hPCcDZwI3AzMyvP8/Ac8AXk8qQb0PuBL4z237rQde15a2McP7WJ31Wg0yqLMfT7WaZ9RnvvZqmrWUpZH/HcBpLcEFgIj4NPDLwEM9nueqiFgcEecB3+/lAEmnACuA34+IL0bEl4HXAKdKeknb7jsi4ua2x+4e82Z1lqUhexAbn/tpqC9i5utB6iQx6iW4muq5BBMRfzXBtgeBM3o8z75e37PFGcAjEXFjy3m+Jem+xrZvTOGcNmiytKvUZfbjLCWSftuN8mwLGbQ2jbxLcO4wkIvMI/kb3ZOPl3RS+6OIDDYsI1V9tbursa3VSZIek7RH0k2STiswX1amrO0qzRvuihXV3Bizlkjq1G40aG0aeZbgPP1NbnouwUg6GLgU+H1S9+ROpueRqQ6OALZ1SN8KHN/y+nbgFuBOYAFwAfB1SadGxLfaD5a0ClgFcMwxx+ScZcvdoLWrZC2R1OnzbdsGBx0Ed9wBjz8Ohx8Ov/iLsGvXgfvW5dd+XiW4UeuBWKAsJZh3khryzwcEvJnUmP5NUiP6y/POXJtOA3bUmh4Rfx0Rl0XEDRFxBfBi4MfAWzueMGIsIpZHxPIFCxYUkmnL0aC1q9x/P/zwh7B2bbpRb98+cYmkTp9v2jS47TZ44gmYMyc933ZbSm81LL/2W9ubbr4Znnxy/+3D3AOxQFkCzCuBi4DPN15/KyI+HREvBW4Czso5b622ksbdtJtL55INABGxC7gaeG4RmbKSDdISzlu2wIMPpl//zRv0HXfAww8fWCJp3txuvTWVGnbvrv7zdRuA3Z4+aFVpnbQHyZkzUzDdvn18nzqXlGssS4BZDNwdEXuB3aRqq6bPAL+TZ8barOfAthbo3jbTrvfpCqy7OvQqqrpdpVcbNqSSR0QKGIcckm7AzfSmTje3vXvh+c+v9vNFwPOeBzNmpBvtjBnpdXuAqVO70VS1B8kTT0zpd99dfUlywGUJMA8xXoq4D/i1lm1L8spQF9cACyWd2kxoTLJ5fGNbR5IOJfUyu63g/A2/YakKKcu2bbBwYerF1LxJH344LF68f9Coawlg7lw4+OCU/1NOSc8HH3zgr/hh6B7cHiTnzEnBdM+e6kuSAy7LQMvrSYMarwI+CvyVpBOAPcBK4LO9nETSLNJAS4CnAbMlndt4fXVE7JR0D3BDRJwPEBFrJa0BPi3pjxkfaHlTRHyjcd45wNeAfwDuAeYDb2m8xyszfE7rpOiGz7o0FOeleeOdMyc9YHxusFZ1HSDYa7ffYRjg2alzxcEHwwtf6Eb9PmUJMG8j3bSJiEskCTgXOBT4G+DdPZ7nKOALbWnN18eROgwcxIE90n4X+CDwCVLJ62vA6pbte4BNwNsb77GbNLXNaRHhWSz7VeSNcNDGXPSi1xtvnXqOtep1HFG/443q8MNiGIJkTWWaTXmYeTblSeQ1W2+nG8qGDdXNBFykXm6eozbrc6s6ffY6BLoBledsyjaq8viV162k8rOfwXHH7b9vHaqJ+tXLuIy6zDhQhTqNN6nDjNBDKMtAyxmkWZR/m9SucUj7PhFxVG45s3rJ40a4YUPqIXXffeOD9448Mv2dRzXRoP4KHdWbW13bnyw3WUowlwGvBr4C/AvwRCE5svrq90Z4//3wk5+km8icOan77r33wuzZqTQE+ZeORqGqqZN+g20Zwbqu7U+WmywB5hzgjyLiI0Vlxobc44+nkeCHHppeH3poCjIR+ZSO6lLdUrV+g21ZwdqN60MvS4DZAtxfVEZsBMyenYLMrl1p4OHu3bBvX0rvt3Tk6pZx/QbbsoL1KLc/jYgsAebdwAWSboiIHUVlyIbY4sUpsDz66PjAw0WLII954FzdMq7fYFtmsB7V9qcRkWU9mL9rTMl/v6TbOHAOsIiIlXlmzobM0qWp+uW44/avEsljCg5Xt4zrN9g6WI+Ogtvaep4qRtIFwJ+QBjQeRpoOv/XhHmQ2sSInqxykiTCzmMr8b/3OylynWZ2tOCVM/9TzQEtJm0iTWr4lhnB0pgdaWu30MxBxEHqRWbVyGjyd10BLAV8bxuBiNiXtN+H582Hz5vxuylkb2/MMCsPYNuKgub8S2tqyzKb8KYqdkt+svtqrqu69d//qhU2b4DOfSc95VTdkmQrfs11PzN/PgUqYCTtLgHkQeJmkb0h6q6Q3tT3emFuuzOqk083pq19NsxI0p9l/9FE44oj0nNe0+1luAHWd9r8u/P0cqIS2tixVZB9oPP8iaSnidkEa7W82XDpVVe3dCxs3poDy+OPp7+OPT3839VvdkKVnnMcBTayf72dYq9ZKGIfUcwkmIqZN8mifXn80jI2lsRzTpqXnsbGqc1R/dVgZM4v2qqrt29PSx1ddBT/4AUyfnhYV++5306/jpn6rG7L0jBuGhb+KNNXvZ9ir1gpeITbLZJe/NsHmfcB20pLKe/rO1aAYG4PVq9M/Okg3ndWNJWpWraouX3U2iHOGtY4L2b4d7rgjzUYwZw48+WRqj5k3L13/5tQ3eY3D6bWxff788Wq7+fPTY/r00RwH1MlUx0l5CqK+ZOmmvI/917ZXh9e7gI8B/zMi9uaVyTJMqZvyokXpptJu4UJ46KF8MjZs8lpXpkz33jt+896+PeX9Rz+CJUtS3rduTdPdPPe5KX3ZsnKrUppBe+/e1Itt8+YUXF7xipTHuiq76mkq77dmTfoh1FoyjUglyhUrisvrAMmrm/JLgI8DVwNfJa0euQA4i7QE8gXAM4G3Aj8jrYA53B55JFu6DV5bwZYtcPfd6Wa0eTOsX5+CyXHHpRLMwoXphrN9e7pZLVp0YKAs+kba+it70aKUtmNHym9dA0wVJdmpdL32rAZ9yRJg3gz8XURc1Ja+RtJFwGsj4uWSDgJeyygEmKOP7lyCOfro8vMymbo0VNblP2yv30f7zVuCxx5LbW67dqV9IuCggzpXuUx0I22eP+s1ac/7Aw/Ascfuv0+dgzYMTtWTpyDqS5Zuyi8Fbuqy7d+AFzX+vhFY1E+mBsa73pWqdlrNnJnS66RODZV1mIYky/fR3sC/eHGaAXrXLjj55FQt9dOfppJCp1/f3brHrls3tWvSKe8PPHDgD526/8rOMsanSsM6BVFJsk7X/wrgGx22vaKxHWAWqcF/+DUb8i+8MFWLHX10Ci51a+Cv06/FOkzRnuX7aC9xzZmTgslPfwpPPQXPec7EJY9uVYK3356OzXpNOuV96dKUPnt2tl/ZZZZq299r2rR6lGR7MYyzGpQkS4C5GLhU0rHAVRzYBvOHjf1eBNza7SSSTiBNmvlC4FnAv0bE6ZO9uaQ5wCXA2aSS19eA1RHxaNt+ZwHvBZYCPwTeFRH/1NtHnIJVq+oXUNrVrd2j03/YMm92Wb6PTlUk06fD2Wf3lr9uVYLN8/WSh8nyvnBhKlE1f2X3ErTLbAPp9F5bt6YS7KJFrnoaYlmm6/+QpB8Dfwb8bePYp4BvA+dExJWNXf+CiZdTPpkUkG4GZmTI6z8BzwBeT+oW/T7gSuA/N3eQdCrwReDDwOrG+3xW0taIuDbDew2XurR7dFN2g2+W76PfEle3OvyTTpraNemW92OOyfYru8xSbaf3WrgwdenOEhRt4GQpwRARXwa+LGkaqfSyKSL2te2zeZLTXBURXwGQdAUwf7L3lXQKsAI4LSJubKT9GLhF0ksiollt9w7gxohoDEbhOkknA+8ERjfA1L2hsqybXbOUdP/98OCD6XtZuHDy76OfKpJuAQqmdk3yupZTLdVOpaTZ7b127XLV05DL0sj/cxGxLyIeaQ8uvR47hbc8A3ikGVwa5/kWcF9jG5JmkqrnPt927OeAUxpVbKOp7g2VZTT4tjaOH3fceLvFxo3Ffx+dRktP9ZrkdS2nMrJ9qp1FPMvAyMpUgqnQMmB9h/S7GtsAlgAHd9jvLlIgPZEJ2oaGXp0bKsuowmsvJS1alBrFqxzgOdVrkse1nKgk1K2UMtWSZpUl6Lp0zx9RUyrBVOAIDlyiGWBrYxstz+37bW3b/nOSVklaJ2ndpk2bcsimTUkZXZcHpVtsWbqVhKB7KWWq32FVJeg6dc8fUYNSgoH9p6Vpap+uptN+6pJORIwBY5Cmiuk3gzZFZXRdnqiUNKq/cjuVhJqdLTqVUvopaeZVgs5yrerUPX9EDUoJZiswt0P6XMZLLFtb0tr3gc4lIKuLgmd17VpKmj/fv3JbTVRKqXqQbNYSiUutlRuUALOe8baWVq1tM/cCT3bYbxmpW/PdheXO6q9bNc3mzV6IqtVEDfJVdxbJumiYOxdUblCqyK4B3iHp1Ii4CUDScuD4xjYiYo+k64DzgMtbjl0JrI2I0ZhdwLrrVE1Tt0GoRchSrTRZg3yVnUWyXqu6d88fAaWXYCTNknSupHOBpwELmq8lzWrsc4+kjzePiYi1wBrg05LOkXQ28BngppYxMADvAU6XdImk0yVdTBps+e6SPp4NmmH/lZu1WqnqUspEsl6rOn+WEVFFCeYo4Attac3XxwEbSflqXyHzd4EPAp+gZaqY1h0i4qZG4Hov8EbSOJlXjfQo/mHXbwP9sP/KnUpDd126tLdf2/nz09IJ0Pu1qstnGVE9Lzg27Ka04JhVq3WKmdYbTtZfqXXoRVZUHuq0YFaWz9jt2p54Ymo3G7UefzWW14JjZvWSVzfUqn/lFjkXW13mocv6Gbtd282bu1+rOvxQsP0MSi8yswPVvRtq86a6Zs3E7R5Ze0dlUXXX4qasn7F5bbdvhzvugLVr4Yc/TPPIdeJBlbXkAGODq84N9P0sagb5Bcq8G7p7DZrtsn7GuXPTImp33AFPPJHW4Xn88TRJaaf3LDJI25Q5wNjgqsuv806y3PCmTUsrXK5dm26o27fnGyjzGsTaTykh64+B5mSkEhxySJraP2I8vV3dS7MjygHGBledu6H2esPbsiUtvvXYYzBjRrppr1sHDz1Uj0DZqp9SQtYfA/PmpeWpDz88BdwZM1JvsYULOweNOpdmR5gb+W2wVd1A302vjesbNqSb5vz5cOedafmAPXvSr/a66WdQ6lTmm1u8OH0Xrd/hjh2dg8awdzcfUC7BmBWh11/srSWdffvgGc+A5z8/tTvUrZG631JC1qq6LKWeOpdmR5hLMJa4i2e+ev3F3rxpP/AAHHpoeuzaBQsWjFc/1aWEVnYpIWupp66l2RHmAGPFjsMYFEUE2F5ueM2b9qZNcNRRKbjs2gVLltRvTrQyllXo9p42kBxgzOtmVBlgmzfQhx6Cn/40lVyWLEndcru1N1TJN3zLwAHGqptROI9SQx7nqDrAzpsHZ5+9/9QozfYGN1LbAHMjv1XTxTOPkdd5jd6uwxiKqpcVzjpw0qwHDjBWzYDFPEZe5zV6uy5jKIpe1bOdp1exgjnAWDW/nvMoNeRV8qjzjABF8vQqVjC3wVhSduNtHrP85jVTcBW9o+pgFFbztEo5wFg18hhTkee4jFHsHVWXqfxtaDnAWDXyKDWUWfIYxoGodZ5eZRi/7xHkFS0bvKKldZXXypl1VMcb+TB/30PIK1qa9aPqcTJFqmPV4DB/3yOm9F5kkk6S9E1JOyX9RNK7JU3v4biTJV3bOG6zpMsk/ULbPp+SFB0ey4r7RDb0eumt5vEk+anDuCTLRaklGElHAN8A7gTOApYA7ycFurdPcNwc4F+Au4GVwJHAxcAi4Oy23dcDr2tL29h35m10TdYY7rnc8uXOB0Oj7CqyNwCHAudExGPA1yXNBi6SdHEjrZM3NY57eURsA5C0BfiKpOUR0dp4siMibi7uI9jImawx3FU6+apz5wPLpOwqsjOANW2B5HOk4HHaBMc9G1jXDC4N1wIBvCznPJrtb7KBqK7SyZfXdhkaZZdglpGqun4uIu6XtLOx7aouxx0CPNGW9hSwD3hmW/pJkh4DZgK3Am+LiBv6zbiNuIkaw/Ou0snSs6uOvcDyUMfOB5ZZ2SWYI4BtHdK3NrZ1cw/wHyUd3JL2PGA60Pq/6XbgAuDlwKsb278u6Vc6nVTSKknrJK3btGlTzx/CbD95TjWTZX4wzyVmNVfFXGSdBt6oS3rTR4EFwN9IWijpZODDwN7GI5044q8j4rKIuCEirgBeDPwYeGvHjESMRcTyiFi+YMGCKX4cG3l5VulkmR8sj7nE3PvNClR2FdlWYG6H9Dl0LtkAEBHrJa0CPgj8d1LV2BgpKD0ywXG7JF1NKtGYFadblU7WKqws84P1O5fYMPZ+G9YqwwFVdglmPamt5eckLQYOa2zrKiI+ARwN/DLwH4A3AycAvfQY83QFVr6pVGFlWTqg32UGhm02ZVcZ1k7ZAeYaYIWkw1vSVgK7gEkb4iNid0R8LyIeAV5Dyv/nu+0v6VBSz7Xb+sq12VRM5QaepT2n37afYev9NmwBcwiUXUX2EWA18CVJ7wOOBy4CPtDadVnSPcANEXF+4/Vs4G3AjaTeYy8iNeb/QURsaewzB/ga8A+kTgHzgbcATwNeWcaHsxEzWXXMRFVY3Y7NMoFnv5N9DtuARi8/UDulBpiI2Crp14EPkbokbyO1q1zUIV+t08fsBZ4D/AFpzMwdwHkRcWXLPnuATaQZAY4CdgNrgdPaBmJa3Q1CPXov7RfdbuDTpk18bJYuuv105x22AY3DFjCHgGdTbvBsyjUxKDPpNuv6W29mO3aM9yCD7p9l+nQ45JCJjy3LIATzXg3Kv50hM9Fsyl4y2eplUOrRe2m/6NZ9OaI+bR/NPK5YMfg3Ys8AUDuert/qZVDq0XutjulUheWqnOJ4BoBacQnG6qXfrrdl6acHV54j/81qzAHG6qWuN9/2Ee8w9eoYV+XYiHAjf4Mb+WukDg3PrXmYNg22boWFC914bNbGSybbYKm6Hr29C/K6dfDYYzB//njHA8h/vZc6BFazHLmKzKxde0+2p55KN/wHHhjfJ+9eX4MwzYknxrSMHGDM2rV3QT68MbPR44+Pp+Xd8aDu3bMHIQBa7TjAmLVr78m2eHEKOgcdVFzHg7rPC1b3AGi15ABj1q69J9tBB8HTnw7HH19cr6+6d8+uewC0WnIjv1m7TpNI/sZvFNvgXvd5wTw41KbAAcask7J7svU7M3LR6h4ArZYcYMzqouru2ROpewC0WnKAMbPe1DkAWi25kd/MzArhAGNmZoVwgDEzs0I4wJiZWSEcYMzMrBClBxhJJ0n6pqSdkn4i6d2Spvdw3MmSrm0ct1nSZZJ+ocN+Z0n6nqTdku6UtLKYT2JmZhMptZuypCOAbwB3AmcBS4D3kwLd2yc4bg7wL8DdwErgSOBiYBFwdst+pwJfBD4MrAbOBD4raWtEXJv/JzIbMF4SwEpU9jiYNwCHAudExGPA1yXNBi6SdHEjrZM3NY57eURsA5C0BfiKpOUR0Vwp7B3AjRGxuvH6OkknA+8EHGBstLWvc7NzZ3rthdOsIGVXkZ0BrGkLJJ8jBY/TJjju2cC6ZnBpuBYI4GUAkmYCLwI+33bs54BTGqUgs9HlGZGtZGUHmGXA+taEiLgf2NnY1s0hwBNtaU8B+4BnNl4vAQ5uPz9wF+lznji1LJsNCc+IbCUru4rsCGBbh/StjW3d3AO8StLBEfFkI+15wHSgWbZvHt9+/q1t239O0ipgFcAxxxwzSdYHhOvYrRvPiGwlq6KbcnRIU5f0po8CC4C/kbSw0a7yYWBv4zHR+dXtfSNiLCKWR8TyBQsW9JT5WvOqgzaR9nVuilg4zaxF2QFmKzC3Q/ocOpdsAIiI9aSSxn8BHgK+C3wL+DbwSMu56XD+5uuu5x8armO3iTQnq5w5s7iF08xalF1Ftp62thZJi4HDOLDtZD8R8QlJ/wgsBX4KbAYeBT7W2OVe4MnG+W9oOXQZqa3m7hzyX2/btqWSS6tZs9LNxAw8I7KVquwSzDXACkmHt6StBHaxf1DoKCJ2R8T3IuIR4DWk/H++sW0PcB1wXtthK4G1EbE9h/zXW92X3TWzkVJ2gPkIsAf4kqSXNBrZLwI+0Np1WdI9kj7e8nq2pPdJepmkFZL+D6nksjoiWhsY3gOcLukSSadLupg02PLdJXy26rmO3cxqpNQqsojYKunXgQ8BV5HaRT5ICjLt+WqdPmYv8BzgD0hjZu4AzouIK9vOf5Okc4H3Am8E7gNeNTKj+L3qYHXce8/sAIqYqPPW6Fi+fHmsW7du8h3N2rWOkG9dr94N6DYCJN0WEcs7bfNsymb9cu89s44cYMz65RHyZh05wJj1y733zDpygDHrl3vvmXXkAGPWL4+QN+uo7JH8ZsPJI+TNDuASjJmZFcIBxszMCuEAY2ZmhXCAMTOzQjjAmJlZITwXWYOkTcCPCjr9HCCv5QKmcq5ej5lsv4m2d9vWa/p80ho/VcnrGlV5fSbbJ8s16pRW5TXy/6HJ06u6Pk+PiM5LAkeEHwU/gLEqz9XrMZPtN9H2btt6TQfWDcM1qvL65HmNuqRVdo38f2jy9Kr/D3V6uIqsHFdVfK5ej5lsv4m2d9uWNb0qeeWnyusz2T5ZrsWwXp+pnsv/h6bAVWRWC5LWRZcpv60efI3qrY7XxyUYq4uxqjNgk/I1qrfaXR+XYMzMrBAuwZiZWSEcYMzMrBAOMGZmVggHGBsIkp4l6f9J2iDpq5IOrzpPNk7SZZJ+LMmNujUkabGkb0q6S9L3JV0sSUW/rwOMDYqPAG+PiKXAeuBPK86P7e+zwHOrzoR19RTwvyLimcBzgBcA5xT9pg4wVhhJJ0i6XNJ3JO2VdH2X/U5q/LraKeknkt4taXrL9qOB4yLi6kbSx4HfKf4TDLe8rg9ARNwYEY+UkvERktc1ioiHImJd4+8ngO8Ci4vOv1e0tCKdDJwJ3AzM6LSDpCOAbwB3AmcBS4D3k378vL2x2y8CD7Ycdj8l/OcYAXldHytO7tdI0pHA2cBLC8lxCwcYK9JVEfEVAElXkCbja/cG4FDgnIh4DPi6pNnARZIubqQJcN1+/vK6PlacXK+RpJnAFcAlEXFX0Zl3FZkVJiL29bDbGcCathvV50j/YU5rvH6Q/Ussx7B/icamIMfrYwXJ8xo1qsw+A9weEe/PNaNdOMBY1ZaRGu1/LiLuB3Y2thERDwMbJZ3Z2OV84EtlZnKETXp9rHK9XqPLgceBC8rKmAOMVe0IYFuH9K2NbU1vBP5c0gbgJODi4rNm9Hh9JH1M0oONvx+U9LFysmf0cI0k/Srph9ly4HZJ35a0uuiMuQ3G6qBT+8p+7S4R8V1S90orXy/X5/XlZcc6mPAaRcS/NV6XyiUYq9pWYG6H9Dl0/lVm5fL1qb/aXiMHGKvaetrq8iUtBg6jrV7ZKuHrU3+1vUYOMFa1a4AVbVO/rAR2ATdUkyVr4etTf7W9Rm6DscJImkUaJAbwNGC2pHMbr6+OiJ2kKWBWA1+S9D7geOAi4AMeY1EsX5/6G/Rr5AXHrDCSjgXu67L5uIjY2NjvJOBDwCmkOuOPARdFxN7iczm6fH3qb9CvkQOMmZkVwm0wZmZWCAcYMzMrhAOMmZkVwgHGzMwK4QBjZmaFcIAxM7NCOMCYmVkhHGDM+iDpIkmb29KmSfqMpN2SOi5LK+l0SdF4bCsls11I2tiSl9+qMi82XBxgzHIkScBHgfOA8yLi2kkOeTXwksIzNrHfBs6pOA82hDwXmVm+PgT8PrAyIq7qYf/vRsQdBedpQhFxu6StVebBhpNLMGY5kfR+4A3A70XEF6d4jsMkfUjSDyTtlHSfpL+VNHuS417bqOL6hbb0jZL+aip5MeuXSzBmOZD058BbgPMj4h/7ONUsYDrwNmATsLjx9xeAFf3m06xMDjBm/TsSeCvwwYj4ZD8niohNwBubryUdRJpN9yZJx0TE/X3l1KxEriIz699jwC3A+ZKe3e/JJP1XSbdL+hnwJHBTY9OJ/Z7brEwOMGb9exJ4GfAT4BpJx0/1RJJ+G/g0sJbUE+2FpF5eAIf0mU+zUjnAmOUgIh4FXgo8BayRdNQUT3UecEtEvCkiromIW4BeenjtbjzPaEs/Yor5MOubA4xZTiLiAeA3SW0y17Stkd6rQ4E9bWmv7uG4BxvPz2wmSHoBMGHvM7MiOcCY5Sgivg/8FulG/2VJ7SWKyXwd+DVJb5P0EkkfAH69fSdJ35T0zZakbwE/Bi6VdKak15AGfD7Wdtw7JT2VMU9mU+IAY5aziPh34JXAacDfS8ry/+xy4P3A/wC+BDwdeFWH/aY3Hs33fILUVrMPuAK4gNQbrb16bVrrcWZFUkRUnQezkSPpdOA64NnAHRGxt8K8TAeOBe4BXh4RX6sqLzZcXIIxq9a3gUcrzsO9pOBiliuXYMwq0OgA8IzGy6ci4tsV5uWXgJmNlxsiYntVebHh4gBjZmaFcBWZmZkVwgHGzMwK4QBjZmaFcIAxM7NCOMCYmVkh/j8jsrc5CNNtYwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "index = 695\n",
    "\n",
    "gs_all, Ks_all = all_params[index]\n",
    "plt.scatter(Ks_all, gs_all, color='r', alpha=0.2)\n",
    "gs_pareto, Ks_pareto = [i[0] for i in all_paretobugs[index]], [i[1] for i in all_paretobugs[index]]\n",
    "plt.scatter(Ks_pareto, gs_pareto, color='r')\n",
    "plt.xscale(\"log\")\n",
    "plt.xlabel(\"K [a.u.]\", fontsize=15)\n",
    "plt.ylabel(\"$g_{max}$ [hr$^{-1}$]\",fontsize=15)\n",
    "plt.xticks(fontsize=15)\n",
    "plt.yticks(fontsize=15)\n",
    "plt.ylabel(\"gmax\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "7126da4a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEQCAYAAABm79njAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAbKklEQVR4nO3debRkZX3u8e8DBgSUBm0QZRCwRS+YweRonCJEHPF2iASH3CTqFW2nOA/xZrG0URMXEMBEY6BFr7qiAeGSaCtIBGmMuUEFFa8SMKINQRwa6G6URojwu3/sfbQo63RXddc+Vef097NWrT57/p2iqOfs/b773akqJEnq0g6TLkCStPgZNpKkzhk2kqTOGTaSpM4ZNpKkzt1r0gVMo6VLl9aBBx446TIkaUG54oorbqqqvQYtM2wGOPDAA7n88ssnXYYkLShJrptr2aK/jJbkEUm+kuQ/knwyyX0nXZMkbW8WfdgApwPHV9VDgauBN0+4Hkna7kxl2CRZluSMJFcmuSvJmjnWOzTJxUk2JbkxyduT7Niz/AHAQVV1fjvrA8AfdP8bSJJ6TWubzWHAUcBlwE6DVkiyJ3ARcBVwNPAQ4BSaAD2+XW0/4Iaeza4H9u+mZEnSXKY1bFZX1ScAkpwLLB2wzsuAXYBjqupW4LNJdgdWJjmpnRfAwd8kacKm8jJaVd09xGrPAC5sQ2XWWTQBdHg7fQP3PJM5gHue6UiS5sFUhs2QHk7T4P9zVXU9sKldRlX9AFib5Kh2leOA8+azSEnSwg6bPYENA+avb5fNejnwF0n+AzgUOGnQzpKsSHJ5ksvXrVs37lolabs2rW02wxrUHnOPdpqq+jrwyC3uqGoVsApgZmZm69t5Vq3a6k3nxYoVk65A0nZoIZ/ZrAf2GDB/CYPPeLYoyfIkqzZu3LgNZUmS+i3ksLmatm1mVpL9gd3oa8sZVlWtrqoVS5YsGUN5kqRZCzlsLgCe1jf8zHOB24FLJ1OSJGmQqWyzSbIrzU2dAPsCuyc5tp0+v6o20QxD82rgvCQnAgcDK4FT+7pDS5ImbCrDBtgbOKdv3uz0QcDaqlqf5EjgvcBqmnaa02gCZ6skWQ4sX7Zs2dbuQpI0wFSGTVWtpelVtqX1rgKeNMbjrgZWz8zMvGRc+5QkLew2G0nSAmHY9LDrsyR1w7DpYddnSeqGYSNJ6pxhI0nqnGHTwzYbSeqGYdPDNhtJ6oZhI0nqnGEjSeqcYSNJ6pxh08MOApLUDcOmhx0EJKkbho0kqXOGjSSpc4aNJKlzho0kqXOGTQ97o0lSNwybHvZGk6RuGDaSpM4ZNpKkzhk2kqTOGTaSpM4ZNpKkzhk2Pez6LEndMGx62PVZkrph2EiSOmfYSJI6Z9hIkjpn2EiSOmfYSJI6Z9hIkjpn2EiSOmfYSJI6Z9hIkjpn2PRwuBpJ6oZh08PhaiSpG4aNJKlzho0kqXOGjSSpc4aNJKlzho0kqXOGjSSpc4aNJKlzho0kqXOGjSSpc4aNJKlziz5skvxdku8lqUnXIknbq0UfNsA/AL856SIkaXs272GTZFmSM5JcmeSuJGvmWO/QJBcn2ZTkxiRvT7LjqMerqs9X1Q+3uXBJ0la71wSOeRhwFHAZsNOgFZLsCVwEXAUcDTwEOIUmHI+fnzIlSeMyibBZXVWfAEhyLrB0wDovA3YBjqmqW4HPJtkdWJnkpHYeSb4A7Ddg+4ur6rhuypckjWrew6aq7h5itWcAF86GSuss4ETgcGB1u68njL9CSdK4TWsHgYcDV/fOqKrrgU3tsrFLsiLJ5UkuX7duXReHkKTt1rSGzZ7AhgHz17fLhpbkzCQ3tD/fkOTMQetV1aqqmqmqmb322mvUeiVJmzGJNpthDbovJnPMn3snVS8eTzmSpK01rWc264E9BsxfwuAznrFIsjzJqo0bN3Z1CEnaLk1r2FxNX9tMkv2B3ehryxmnqlpdVSuWLFnS1SEkabs0UtgkeURXhfS5AHhakvv2zHsucDtw6TzVIEkak1HPbL6e5MtJXp5kj605YJJdkxyb5FhgX2Cv2ekku7arnQ7cAZyX5MlJVgArgVP7ukOPlZfRJKkbo4bNkTR39Z8E3JjkH5I8JUlG2MfewDnt6zHAoT3TewNU1fr2WDvS3FNzAnAa8LYR6x2Jl9EkqRsj9UarqkuAS5K8guay1guBC4EbknwY+FBVXbuFfayl6VW2pWNdBTxplPokSdNpqzoIVNVtVfXBqnoicAiwFvhz4FtJLk3yrDHWOG+8jCZJ3djq3mhJDkyyEvhn4LHA+cAK4IfA2UlOG0uF88jLaJLUjVF7o+2a5PlJLgG+DfwR8H7ggKpaXlUfqKrnAC8FHAhTkgSMPoLAD2ga7c8DnlxVa+ZY78vAzdtQlyRpERk1bN4CfLSqNtuoUVXfAA7a6qomJMlyYPmyZcsmXYokLSojXUarqvdtKWgWMttsJKkbI53ZJPkgsFtVPXfAsn8AbnPgS6nPqlWTrmBuK1ZMugJtJ0btjfYU4Nw5lv0f4KnbVo4kaTEaNWz2Am6ZY9l62hEAJEnqNWrYXAc8cY5lTwRu2LZyJsubOiWpG6OGzYeAP0vyyiT3AUhyn3b4mjcDA5+CuVDYQUCSujFq1+cTgYcA7wH+JsltNM+YCbCqXS5J0j2MOhDn3cCLk5xMM0jm/Whu3vxcVX2rg/okSYvAqGc2AFTVNcA1Y65FkrRIbVXYJDkE2A+4d/+yqjp/W4uSJC0uo97UeShwNs0DzwY9k6Zoxk5bkByuRpK6MWpvtDOAnYBjgIfRjH/W+zp4rNXNM3ujSVI3Rr2M9kjgeVX1qS6KkSQtTqOe2VzLgHYaSZI2Z9SweQPw50kW9OUySdL8GvUy2ruAfYGrk6wFNvSvUFWP3vayJEmLyahh8432JUnS0EYdQeB/dlWIJGnxGrXNBoA09k/yuCS7jbuoSXHUZ0nqxshh047w/D2axw38C839NiQ5L8lrx1rdPPM+G0nqxkhhk+RNwKnA+2kG4uwdRWAN8EuPi5YkadQOAq8E3lpVJyXpH5bmGuCQ8ZQlSVpMRr2Mtg9wxRzL7sYbPiVJA4waNt8GDp9j2ROBq7atHEnSYjTqZbR3A+9Lcidwbjtv7yTHAa8HXjLG2iRJi8So99mcmWRP4K3ACe3s84FNwMqq+tiY65MkLQIjPzytqk5OcjrwOOD+wC3Av1WVN6dIkgba2sdC/xi4cMy1SJIWqVGf1PmKLa1TVe/b+nImyyd1SlI3Rj2zee9mllX774INm6paDayemZmxo4MkjdFIXZ+raof+F3A/4A+BK4FDuyhSkrSwbVWbTa+q2gCcnWQJcAZwxLbuU5K0uGzVqM9z+C4wM8b9SZIWibGETZIH0jwy+rvj2J8kaXEZtTfaOn7REWDWTsB9gZ8Cx4ypLknSIjJqm83f8sth81PgBuAzVXXzWKqSJC0qow5Xs7KjOiRJi9g4OwhIkjTQqG023+WXL6PNqaoOHrkiSdKiM2qbzbnA84Bdgc8CPwL2Bp4C3AacPdbqJEmLwqhhsx64FnhmVd02OzPJfYBPARur6p1jrG+bJdkf+BDwIJqniX4a+LOqGvoMTZK0bUZts3klcHJv0ABU1U+Av2qXT5uf0YTLfwMeCfw2dtGWpHk1atgsAR4wx7J9gPsMs5Mky5KckeTKJHclWTPHeocmuTjJpiQ3Jnl7kh1HKbiqvl9Vl7c/3wl8Hdh/lH1IkrbNqJfRPgmcnORWYHVV3ZFkZ+D3gBOB1UPu5zDgKOAymptCf0n7RNCLgKuAo4GHAKfQBOTxI9Y9u8/7A78PPHVrtpckbZ1Rw+blNO0fHwcqyY9pRg8ITRC9fMj9rK6qTwAkORdYOmCdlwG7AMdU1a3AZ5PsDqxMclI7jyRfAPYbsP3FVXXc7EQbiucC766qfx+yTknSGIx6U+dG4FlJDgMeRXNJ7QfAl6vqqhH2c/cQqz0DuHA2VFpn0ZxBHU57FlVVT9jSjtpLbx8FvlpVpwxbpyRpPLb2sdDfBL455lr6PRz4XN9xr0+yqV027CU7aB598GOawUIHSrICWAFwwAEHjFysJGluI48gkGTvJCe2DffXtGc5JHlNkseOsbY9gQ0D5q9vlw0lyeOB42gef/DVJF9L8ur+9apqVVXNVNXMXnvttZUlS5IGGXUEgUfT3My5DriU5kFpO7eLZx8zcOwY6xt0L0zmmD94B1X/2m4jSZqQUc9sTgMuAQ4BXso9v8S/BDx6THVBcwazx4D5Sxh8xrPNkixPsmrjxo1d7F6Stlujhs1vAu9rG/j7zy5uphm6Zlyupmmb+bl2NIDd2mVjV1Wrq2rFkiVLuti9JG23Rg2bjcBcDRoHAz/ctnLu4QLgaUnu2zPvucDtNJfwJEkLxKhh8wnghCS9ozlXkqXAG4HzhtlJkl2THJvkWGBfYK/Z6SS7tqudDtwBnJfkyW1vsZXAqX3docfGy2iS1I1Rw+YtwK00d/V/vp13OnANzRnHW4fcz97AOe3rMcChPdN7A1TVeuBIYEeabs4n0LQZvW3EmofmZTRJ6saoN3WuT/IY4E9oguA24BbgTOAjVXXHkPtZyxA9xNobRZ80So2SpOkzdNgkuTfNkDR/WVUfAD7QWVWSpEVl6MtoVfVTmiFqRhp1eSGxzUaSujFqm80naUZNXpRss5Gkbow6NtqFNI8YeCBwPk1X53vcb1NV54+pNknSIjFq2Px9++8xDH7aZbGIL7NJkrbOFsMmyT8Dr6qqa4CDaHqRHQl8kWYk5UUjyXJg+bJlyyZdiiQtKsO02TyZZjwyquo64D+BVcBOVXVd/6vDWjtnm40kdWPkRwy0HEVZkjS0rQ0bSZKGNmzYDHp+zNDPlJEkbd+G7Y12YZKf9c27eMA8qmqcjxmYV3YQkKRuDBM2J3RexZSoqtXA6pmZmZdMuhZJWky2GDZVtd2EjSSpG3YQkCR1zrCRJHXOsOnhqM+S1A3DpocjCEhSNwwbSVLnDBtJUucMG0lS5wwbSVLnDBtJUucMG0lS5wybHt5nI0ndMGx6eJ+NJHXDsJEkdc6wkSR1zrCRJHXOsJEkdc6wkSR1zrCRJHXOsJEkdc6wkSR1zrCRJHXOsOnhcDWS1A3DpofD1UhSNwwbSVLnDBtJUucMG0lS5wwbSVLnDBtJUucMG0lS5wwbSVLnDBtJUucMG0lS5+416QK6luRSYA8gwLeAF1XVrRMtSpK2M9vDmc3vVdWvV9WvAdcDb5p0QZK0vZn3sEmyLMkZSa5McleSNXOsd2iSi5NsSnJjkrcn2XHU41XVxnZ/OwC7AbVNv4AkaWSTuIx2GHAUcBmw06AVkuwJXARcBRwNPAQ4hSYcjx/1gEnOBx4FfBN4w1ZVLUnaapO4jLa6qvavqmfTfPkP8jJgF+CYqvpsVZ0OnAC8Psnusysl+UKStQNeH+jdWVUdBewDfAl4RSe/lSRpTvMeNlV19xCrPQO4sK8h/yyaADq8Z19PqKoDB7yOG3Dcu4APA8/fxl9BkjSiae0g8HDg6t4ZVXU9sKldNpQkeyZ5QM+sPwC+MZYKJUlDm9auz3sCGwbMX98uG2U/H0+yE03X538HXjVoxSQrgBUABxxwwCi1SpK2YFrDBgb3Gssc8wfvoOo7wMyQ664CVgHMzMzYY02SxmhaL6Otp7kRs98SBp/xSJKm2LSGzdX0tc0k2Z/mPpmrB24xBkmWJ1m1cePGrg4hSdulaQ2bC4CnJblvz7znArcDl3Z10KpaXVUrlixZ0tUhJGm7NO9tNkl2pbmpE2BfYPckx7bT51fVJuB04NXAeUlOBA4GVgKnOq6ZJC08k+ggsDdwTt+82emDgLVVtT7JkcB7gdU07TSn0QROZ5IsB5YvW7asy8NI0nZn3sOmqtbS9Crb0npXAU/qvKB7HnM1sHpmZuYl83lcSVrsprXNRpK0iBg2PeyNJkndMGx62BtNkrph2EiSOmfYSJI6Z9j0sM1Gkrph2PSwzUaSumHYSJI6Z9hIkjpn2EiSOjfND0+bd46NJqlzq1ZNuoLNW7Gik916ZtPDDgKS1A3DRpLUOcNGktQ5w0aS1DnDRpLUOcOmh8PVSFI3DJse9kaTpG4YNpKkzhk2kqTOpaomXcPUSbIOuG7SdbSWAjdNughNFT8T6jVNn4cHV9VegxYYNlMuyeVVNTPpOjQ9/Eyo10L5PHgZTZLUOcNGktQ5w2b6TfkQsZoAPxPqtSA+D7bZSJI655mNJKlzho0kqXOGzRRKsizJGUmuTHJXkjWTrkmTk+TZST6Z5HtJfpLkiiR/OOm6NBlJjk3yf5PcnOSnSa5JcnySnSZd2+b4WOjpdBhwFHAZMNUfIM2L1wPfBV5Hc/PeUcDHkiytqvdMtDJNwv2BS4CTgQ3Ao4GVwD7An06sqi2wg8AUSrJDVd3d/nwusLSqjphsVZqUNlRu6pv3MeCxVXXQhMrSFEnyF8ArgT1rSr/UvYw2hWaDRgLoD5rWV4G957sWTa2bmfKrIF5GkxamxwFXTboITU6SHYGdgd8EXg383bSe1YBhIy04SY4EjgZeNOlaNFG30YQNwEeAN02wli3yMpq0gCQ5EPgY8Imq+tBkq9GEPQ74HeANNH98vHey5WyeZzbSApHkfsAFwPXAH0+4HE1YVX2l/fELSW4CPpzklKq6dpJ1zcUzG2kBSLIr8CmaRuBnVtVtEy5J02U2eKa2d6JnNtKUS3Iv4BzgocDjq+pHEy5J0+fx7b/fnWgVm2HYTKH2r9ij2sl9gd2THNtOn19VmyZTmSbkfTSfh9cA90vymJ5lX62qOyZTliYhyWeAi4BvAnfRBM0bgLOn9RIaeFPnVGobgef6C+Wgqlo7f9Vo0pKsBR48x2I/D9uZJO8AngUcCPwM+A7wv4HTq+q/JljaZhk2kqTO2UFAktQ5w0aS1DnDRpLUOcNGktQ5w0aS1DnDRpLUOcNGC0aSY5J8LsmGJHck+VaSdyZZOqF6KsnQT0ZM8uYkR3RYz95JVrb3aW3rvp6T5IXbXtU99nlE+549Ypz71cJg2GhBSHIKzZAt3wH+BHgqcBqwHHj/hMp6bFvTsN4MHNFNKUDzMLW30dzst62eA7xwDPvp9RWa92xq73JXdxyuRlMvyXLg9cBxVfXBnkWXJllFEzzzrqoum8Rxp0mSXwHurqq7trRuVd0KbPfv2fbKMxstBK8DvtIXNABU1V1VdcHsdJKlST6c5OYkm5KsSTLTv12SFyf5Zns57rokb+5Z9pgkP0vyop55S5L8Z5K/75l3j8toSZ6Q5F+S3Nq+vpbk2e2ytcD9gbe129XsJbUkOyR5S5Jv91wefEFfvWuSnJvkf7Tr3ZrkgiT7tcsPBP5fu/ols8eY6w1Nsl+Sjyf5UZLbk1zbDoNCkg8BfwAc3lPryr46ViS5Fvgp8KAkD09yVvsebWrf29cm2aHnmL90Ga2dfk2Sv0yyrq3nb5Ps3LPOHknOTHJjkp8muT7JpM5mtZU8s9FUa/9yfhxwypCb/BOwDHgjcBPN0wsvSfLIqvp2u883AX8JnASsAX4LeEeSTVX13qq6LMnJwGlJLqqq64G/ofnj7FVz1Lk7zSMAPgG8HQjwq8Ae7SrPAi4BzgXObOfNPtb5PcAL2u2+AjwF+GCSm6vqUz2H+W3gQTSDLu4C/DWwimaQzu8DfwR8FHglvxhyfi4fafexAtgAHAw8vF32DuCAtvZXtPNu6Nn28cBDgD8DNgEbgUOAa9rj/xj4DeCE9hjv2kItbwA+R/OMnl9r17+O5r8PwKk0n4HXAT8A9geeuIV9atpUlS9fU/sC9gEKeOkQ6z69Xffwnnm7AeuAM9rp3YGfAG/r2/btNF9kO7bTOwFfpxld9+h2v0f1bVPAn7Y/z7TT991MfTcBK/vmLQPuBl7QN/8jwJd7ptfQfKnv2TPvte0xd2mnH9FOHzHEe/UTYPlmlp8LrBkwfw1wO7DPZrYNzR+yfw58p2f+EW19j+h7Dz/ft/0/AZf1TH8DeNWkP4u+tu3lZTQtFMOMGPtoYF1VXfrzjZqHjH0KeEI767E0AXROknvNvmj+sn4AsF+73Z3A82n+gj4bOLOqzt/Msa+l+QL/WJKjk+wx5O91JE3Y/GNfPRcDv5Fkx551v1xV63umZ8+M9h3yWL2+BrwryQuTHDDitldU1Q96ZyS5d5ITknwbuAP4L+AvgIPa32dz/rlv+ira/w49tb4pySuSHDJirZoSho2m3c00X17DfCE+EPjhgPk/BO7X/jzbTfqbNF+Is69L2vn792z3dZovvp1pnikzpzYEngr8CvBxYF2STyc5eAs1LwV2pDlr6a3nQzRnBw/sWXdD37Z3tv/eewvHGOS5wOU0Pfqua9uXjhxy20Hv8Yk0ly5nL+s9CnjnkPVt6Ju+s2+bP6U523krcE2S/0jyvCFr1ZQwbDTVqnk+x78CTxti9e/TdP/t9wDglvbn2X//O80XYv/ryp7tXkPTjnE18J7exu45av23qno6TVvHMTTtGB/bQs230DyT5LfnqKeTp3JW1feq6oU0nRYeS3MJ8ZNJ7j/M5gPmPRt4T1WdVFUXVdXlNL/XOGrdUFWvrqp9gF8Hvgh8NMmh49i/5odho4Xg3cBMfw8t+HlPrqe3k18E9k7yxJ7luwLPBL7Qzvo3mjaHB1XV5QNeP263exjNZaDjab5IH0XTQL1FVXV7Va0GPgj0fiH2/8UOzeW7HYElc9RzJ8Mb+Uynqu6upgv3CcCu/OIhbYNq3ZxdaM5AAWgv/4397KOqvk7T6WMHftGhQQuAvdE09apqdZJTgQ8keTxNj6+f0HzZvAxYC3ymqi5M8q/A2UneQnMJ7o00X4Qnt/va0Hbj/eskDwY+T/PFdQjwu1X1rPaL8sPAV4FTq+ruJG8D3pnk01V1dX+NSZ4JvIjmcs/1NO0oL6UJk1lXA89M81jfnwDXVNU1SU4HzkpyEs2lrXsDhwGHVNWLR3irrqcJ0hck2Qj8V3uG0V/rEuBCmk4I36K5TPgGmrObf++p9egkv0/TE+3GqrpxM8f+LPDKts3mFpoecTtvZv2hJfkC8I80HQUKeAlwG/Clcexf82TSPRR8+Rr2RXPvxyU07Rt30nxR/hU9PaOAvWi+RNfTfPFeCjxqwL7+GLiiXWc9zVnR69tl/4vmy+yhPevvSHNW9EV+0WOttzfaw2h6cP0nzV/4NwCnA/fr2cdv0dzUeBs9vcZoem+9lqYd6Q6a3nOXAs/v2XYNcG7f73AEv9y764/a9+XO5n/vge/jzjSjLlxD03X5JppOFL/as85Smi/4W9pjrJyrjnb+A9r1b6Vp0zmJJhQKuM9m6v35e9gzbyVwU8/0yTT3EP2Ypn3nEuB3Jv159DXay8dCS5I6Z5uNJKlzho0kqXOGjSSpc4aNJKlzho0kqXOGjSSpc4aNJKlzho0kqXP/H/QwsnoPnzEIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "labels = np.array([1, 2, 3])\n",
    "coex = [len(i) for i in all_paretobugs]\n",
    "plt.bar(labels, [coex.count(i)/len(coex) for i in labels], color='r', alpha=0.4, width=0.3)\n",
    "plt.xticks(labels, fontsize=15)\n",
    "plt.yticks(fontsize=15)\n",
    "plt.xlabel(\"Coexistent strains\", fontsize=15)\n",
    "plt.yscale(\"log\")\n",
    "plt.ylabel(\"Frequency\", fontsize=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df68cd27",
   "metadata": {},
   "source": [
    "### Now run the simulation for 2 resource and 2 species. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a454bcaa",
   "metadata": {},
   "source": [
    "#### First, let's fix the functions (consumption ratio of the 2 resources) of the 2 species, and generate the strains only by randomizing the g and K parameters (satellite nches). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ba53c1db",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing only_gK.py\n"
     ]
    }
   ],
   "source": [
    "%%writefile only_gK.py\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.integrate import solve_ivp\n",
    "import copy\n",
    "import pickle\n",
    "from math import *\n",
    "from matplotlib.colors import ListedColormap, LinearSegmentedColormap\n",
    "import seaborn as sns\n",
    "import scipy.stats as stats\n",
    "import time\n",
    "import multiprocessing\n",
    "from multiprocessing import Pool, freeze_support\n",
    "\n",
    "# g: maximum growth rate\n",
    "# K: Monod constant\n",
    "def OneStep(t, u, p):\n",
    "    N, R, gmax, K, x = p\n",
    "    Ni, Rj = u[:N], u[N:]\n",
    "    xmat = np.hstack([x.T, 1-x.T])\n",
    "    cmat = np.vstack([Rj for i in range(N)])\n",
    "    gmat = xmat*(gmax*(cmat/(K+cmat)))\n",
    "#     print(gmat)\n",
    "    return np.hstack([np.sum(gmat, axis=1)*Ni, -Ni@gmat])\n",
    "\n",
    "\n",
    "def Dilute(p, p_dilute, u_end):\n",
    "    N, R, gmax, K, x = p\n",
    "    c0, d, t, n, fluc_thres, elim_thres = p_dilute\n",
    "    output = u_end\n",
    "    u_add = np.zeros(N+R)\n",
    "    u_add[-R:] = np.array(c0)\n",
    "    stop, count = 1, 1\n",
    "    utransf = copy.deepcopy(u_end)\n",
    "    while stop>0 and count<n: \n",
    "        u = copy.deepcopy(utransf)\n",
    "        u = u/d + u_add\n",
    "        N_alive = u[:N]>0\n",
    "        N_temp = sum(N_alive)\n",
    "        gmax_temp = gmax[N_alive]\n",
    "        K_temp = K[N_alive]\n",
    "        u_temp = np.hstack([u[:N][N_alive], u[N:]])\n",
    "        x_temp = x[0:,N_alive]\n",
    "        p_temp = N_temp, R, gmax_temp, K_temp, x_temp\n",
    "        sol = solve_ivp(OneStep, t_span=[0, t], y0=u_temp, method=\"LSODA\", args=[p_temp], rtol=1e-6)\n",
    "        utransf_temp = sol.y[:, -1] * (sol.y[:, -1]>elim_thres*d)\n",
    "        u[:N][N_alive] = utransf_temp[:N_temp]\n",
    "        u[-R:] = utransf_temp[-R:]\n",
    "        utransf = u\n",
    "        output = np.vstack([output, utransf])\n",
    "        stop = sum(abs((output[-1, :-R]-output[-2, :-R])) > fluc_thres*output[-1, :-R])\n",
    "        count+=1\n",
    "    return output\n",
    "\n",
    "def invade(profile, p_profile, p_dilute, p_invade):\n",
    "    # load variables\n",
    "    N, R, gmax, K, x = p_profile\n",
    "    c0, D, t, n, fluc_thres, elim_thres = p_dilute\n",
    "    p_dilute_once = c0, D, t, 2, fluc_thres, elim_thres\n",
    "    gmax_inv, K_inv, x_inv = p_invade\n",
    "    p_profile_temp = N+1, R, np.vstack([gmax, gmax_inv]), np.vstack([K, K_inv]), np.hstack([x, x_inv])\n",
    "    profile_temp = Dilute(p_profile_temp, p_dilute_once, np.hstack([profile[:-R], [D*elim_thres], profile[-R:]]))\n",
    "    if(profile_temp[-1][-(R+1)]<=0):\n",
    "        return profile, p_profile\n",
    "    else:\n",
    "        profile_temp = Dilute(p_profile_temp, p_dilute, np.hstack([profile[:-R], [5*D*elim_thres], profile[-R:]]))\n",
    "        return profile_temp[-1], p_profile_temp\n",
    "\n",
    "def model(index, N_pool):\n",
    "    \n",
    "    gA1, gA2, KA1, KA2 = 1, 1, 4, 5\n",
    "    gB1, gB2, KB1, KB2 = 0.55, 0.55, 0.05, 0.03\n",
    "\n",
    "    g_diff = np.random.normal(0, 0.03, size = 2*N_pool)\n",
    "    K_exp = np.random.normal(0, 0.3, size = 4*N_pool)\n",
    "    gA1new = [gA1+i for i in g_diff[0:N_pool]]\n",
    "    gA2new = [gA2+i for i in g_diff[0:N_pool]]\n",
    "    gB1new = [gB1+i for i in g_diff[N_pool:2*N_pool]]\n",
    "    gB2new = [gB2+i for i in g_diff[N_pool:2*N_pool]]\n",
    "    KA1new = [KA1*(10**i) for i in K_exp[0:N_pool]]\n",
    "    KA2new = [KA2*(10**i) for i in K_exp[N_pool:2*N_pool]]\n",
    "    KB1new = [KB1*(10**i) for i in K_exp[2*N_pool:3*N_pool]]\n",
    "    KB2new = [KB2*(10**i) for i in K_exp[3*N_pool:4*N_pool]]\n",
    "    \n",
    "    # species A: xA = 0.2(consumes 0.2 resource 1 and 0.8 resource 2)\n",
    "    # species B: xB = 0.7(consumes 0.7 resource 1 and 0.3 resource 2)\n",
    "    # species A has g mutating around 1 and K around 4\n",
    "    # species B has g mutationg around 0.55 and K around 0.01\n",
    "    \n",
    "    xA, xB = 0.2, 0.7\n",
    "    sigma = 0\n",
    "    xAdist = stats.truncnorm(-0.2, 0.8, loc=0, scale=sigma)\n",
    "    xAnew = np.sort(xA+xAdist.rvs(N_pool))[::-1]\n",
    "    xBdist = stats.truncnorm(-0.7, 0.3, loc=0, scale=sigma)\n",
    "    xBnew = np.sort(xB+xBdist.rvs(N_pool))\n",
    "    \n",
    "    plt.scatter(xAnew, 1-xAnew, color = 'b', alpha = 0.2)\n",
    "    plt.scatter(xBnew, 1-xBnew, color = 'r', alpha = 0.2)\n",
    "\n",
    "    N, R = 2*N_pool, 2\n",
    "    gmax = np.array([[gA1new[i], gA2new[i]] for i in range(N_pool)]+[[gB1new[i], gB2new[i]] for i in range(N_pool)])\n",
    "    K = np.array([[KA1new[i], KA2new[i]] for i in range(N_pool)]+[[KB1new[i], KB2new[i]] for i in range(N_pool)])\n",
    "    x = np.array([np.hstack([xAnew, xBnew])])\n",
    "    c0, D, t, n, fluc_thres, elim_thres = [10, 10], 10, 24.0, 5000, 1e-5, 1.0e-6\n",
    "    \n",
    "    start=time.time()\n",
    "    \n",
    "    p_profile = 1, R, gmax[:1, :], K[:1, :], x[:, :1]\n",
    "    p_all = N, R, gmax, K, x\n",
    "    p_dilute = c0, D, t, n, fluc_thres, elim_thres\n",
    "    profile = Dilute(p_profile, p_dilute, np.array([1]+[0 for i in range(R)]))[-1]\n",
    "\n",
    "    survivors_0 = [0]\n",
    "    survivors = []\n",
    "    while(survivors_0 != survivors):\n",
    "        survivors_0 = survivors\n",
    "        for i in range(N):\n",
    "            _, _, _, K_profile, x_profile = p_profile\n",
    "            if(x[0, i] not in x_profile[0] or K[i] not in K_profile):\n",
    "                p_invade = gmax[i:i+1], K[i:i+1], x[:, i:i+1]\n",
    "                profile_new, p_profile_new = invade(profile, p_profile, p_dilute, p_invade)\n",
    "                N_temp, R_temp, gmax_temp, K_temp, x_temp = p_profile_new\n",
    "                survivor = profile_new[:-R]>0\n",
    "                p_profile = sum(survivor), R_temp, gmax_temp[survivor, :], K_temp[survivor, :], x_temp[:, survivor]\n",
    "                profile = np.hstack([profile_new[:-R][survivor], profile_new[-R:]])\n",
    "\n",
    "        end=time.time()\n",
    "\n",
    "        N_temp, R_temp, gmax_temp, K_temp, x_temp = p_profile\n",
    "\n",
    "        xlist = list(xAnew)+list(xBnew)\n",
    "        survivors = [list(K[:, 0]).index(i[0]) for i in K_temp]\n",
    "\n",
    "        all_survivors=survivors\n",
    "        all_params=[xAnew, xBnew], gmax_temp, K_temp\n",
    "        outputFileName = 'data/2_res_new_results_gK_separated_'+str(index)+'.pkl'\n",
    "        pickle.dump((all_survivors, all_params),open(outputFileName, \"wb\"))\n",
    "\n",
    "    return (end-start)\n",
    "\n",
    "runs = 500\n",
    "\n",
    "if __name__ == '__main__':\n",
    "    freeze_support()\n",
    "    cores = multiprocessing.cpu_count()\n",
    "    pool = multiprocessing.Pool(processes=7)\n",
    "    x1 = [i for i in range(runs)]\n",
    "    x2 = [50 for i in range(runs)]\n",
    "    tasks = [(index, N_pool) for index, N_pool in zip(x1,x2)]\n",
    "    pool.starmap(model,tasks)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "1fae0828",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAawAAAGyCAYAAAC4FPBHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAlrElEQVR4nO3deZxkZX3v8c9XjCzegGFTjCCYhIh6b6LBEdEgikZcJopxwbiBy5BcYtBIFNEoxEQUQYyJKGNUwCW4JAgorigaLxEEwY2IsokQlgFUVEQUfvePc5qpKbqrT0NVV5/pz/v16tfUWerUt6t6+tfPc57znFQVkiQtdXeZdgBJkrqwYEmSesGCJUnqBQuWJKkXLFiSpF6467QDjIHDHCWtjzLtAEtNbwtWkpXAypNPPpmVK1d2es7q1ZPNtFysWjXtBJKWI7sEJUm90NsWVlWdApwCvGTaWSRJk9fbFlaSlUlWn3LKKdOOIklaBL0tWJKk5cWCJUnqBc9hSZJ6obctLM9hSdLyYgtLktQLvW1hSZKWl94WLLsEJWl5sUtQktQLvW1hSZKWFwuWJKkXeluwPIclScuL57AkSb3Q2xaWJGl5sWBJknrBgiVJ6gULliSpFyxYkqRe6G3Bcli7JC0vvS1YkqTlpbfXYWn9snr1tBOsP1atmnYCaTJsYUmSeqG3LSxnupCk5aW3LSwHXUjS8mILS5LUC7awJEm9YAtLktQLtrAkSb1gC0uS1Au9bWFJkpYXC5YkqRd6W7A8hyVJy4vnsCRJvdDbFpYkaXmxYEmSesGCJUnqBQuWJKkXLFiSpF7obcFyWLskLS+9LViSpOXFgiVJ6gULliSpF5zpQpLUC71tYTnoQpKWF1tYkqResIUlSeoFW1iSpF6whSVJ6gVbWJKkXuhtC0uStLxYsCRJvdDbLsEkK4GVJ598MitXrpx2HGm9tnr1tBOsP1atmnaC/uptwfIcliQtL3YJSpJ6wYIlSeoFC5YkqRcsWJKkXkhVTTvDnZLk08CW084hSWN2bVXtOe0QS0nvC5YkaXmwS1CS1Au9vQ5rxg033GATcYnadNNNpx1B6rPsueee9elPf3raOaYhs620hSVJS9S111477QhLigVLktQLFixJUi9YsCRJvWDBkiT1ggVLktQLFixJUi9YsCRJvWDBkiT1ggVLktQLvZ+aSZLWW2vWwOrVs29btWpxsywBtrAkSb1gwZIk9YIFS5LUCxYsSVIvWLAkSb1gwZIk9YIFS5LUCxYsSVIvWLAkSb1gwZIk9YIFS5LUCxYsSVIvWLAkSb1gwZIk9YIFS5LUCxYsSVIvWLAkSb1gwZIk9YIFS5LUCxYsSVIvWLAkSb1gwZIk9YIFS5LUCxYsSVIvWLAkSb1gwZIk9YIFS5LUCxYsSVIv3HWuDUkOX8BxqqpeNYY8kiTNas6CBTxjAccpwIIlSZqYOQtWVe2wmEEkSRrFc1iSpF4YdQ7rAQs5UFWdf+fjSJI0u1HnsL5Nc25qPmn322AsiSRJmsWogvXoRUshSdI8Rg26+NJiBpEkaRQHXUiSesGCJUnqhVHnsBZNkh2AHYBLquqSaeeRJC09i9bCSvIbSd6e5PokP5uZ+inJO4ALgc8DFyb5YBJHHEqS1rGYLaxXAi8G3gpcD7w0yVbAU4F9gK8DjwSOAPYDjl7EbJKkJa5TCyvJTkl2GVjeOMkbk3w8yUs7vtbzgNdX1Wur6q3t8guAQ6vq/VX1nao6hqagvXCePKuSnJ3k7Pe9730dX16S1GddW1hHA2cAX22Xj6BpFf0n8OYkG1XVW+Y5xn2BswaWz2n/PWtov68AI4tgVa0GVgPccMMNXS5uliT1XNdzWA8C/guac1HAc4GXVdWewMHM0yJq/Ry4x8DyL9uvG4f224AlMhhEkrR0dC1YdwduaB/v0i7/R7v8dZrW03wuAHaeWaiqW6tq46o6b2i/BwKXdswlSVomuhasi2kKFcBewLlVdV27vCXw0w7HeCvwww77PRY4sWMuSdIy0bXr7SjgnUmeATwY2Hdg2+7AN+c7QFX9e5cXqqondMwkSVpGOhWsqnpPku8DDwUOqqrTBjZfD7xtAtkkSbpN58ENVfVl4MuzrD9knIEkSZpN55kukmyd5M1JTkvyvSQPbNcfkOThk4soSVL3C4dX0Eyf9Gc0I/h+B9iw3bwN8IpJhJMkaUbXFtZRwBeAHWmmTcrAtrOAFWPOJUnSOrqew3oI8JSqujVJhrZdB2w93liSJK2rawvrJ8BWc2y7H3D1eOJIkjS7rgXrJODQJPcbWFdJtgQOZO2sF5IkTUTXgnUQzdRM57N2aPu7aKZb+gXwuvFHkyRpra4XDv+ovb3I84A9aCayvR74V+D4qvrl5CJKkrSwC4dvBt7TfkmStKg6XzgsSdI0zdnCSnIN8PiqOjfJGmDkjRKryqHtkqSJGdUl+A7WDld/B/MULEmSJmnOglVVhw48PmRR0kiSNIeucwlum+Qhc2x7SJJtxxtLkqR1dR108U7guXNs+3Pg6PHEkSRpdl0L1i40k9/O5ovtdkmSJqZrwdqE0YMu7j6GLJIkzalrwfoW8Ow5tj0b+M544kiSNLuuM128Cfj3JBsCxwJX0ty48QU0N3X8s4mkkySp1XUuwROTvAA4jKY4Fc1NHK8AnltVH59YQkmSWNhcgu9P8gHg/sDmNDduvKCqvKBYkjRxnQsWQFVVku/SdAdeY7GSJC2WzpPfJnlikjOBm4DLgP/Trl+dZK5rtCRJGouuM108HzgZ+C6wauh53wdeNP5okiSt1bWF9RrgLVX1AuADQ9u+AzxgrKkkSRrStWDdF/jcHNtuAjYdTxxJkmbXtWD9EHjwHNt2Bi4cTxxJkmbXtWC9B3h9O7hi43ZdkuwBvBJ49yTCSZI0o+uw9jcD2wLHAbe0684ANgCOqaq3TyCbJEm36TrTRQH7JzkK2APYArge+EJVfW+C+SRJAhZ+4fCFeL5KkjQFnQtWkrsB+wAraGa6uBI4Eziuqm6eSDpJklpdLxzeieYC4XcAD6I5j/WgdvnCJF6HJUmaqK4trNXAT4A/rqrLZlYm2Q74JPAuYLfxx5MkqdG1YO0MPHuwWAFU1WVJXgd8aOzJJGm522orWLVq2imWjK7XYV0KbDTHto1oJsOVJGliuhasg4B/SPKwwZVJdgH+HnjVuINJkjSoa5fga2nmCzwjyTXANcDW7dd1wMFJDp7ZuapWjDuoJGl561qwvt1+SZI0FV1nuth30kEkSRql8x2HhyX5rSR/mGTDcQaSJGk2XS8cPjTJmwaWH0MzMvAc4OIkD5xQPkmSgO4trOcA3x1YPhL4CvCIdv1hY84lSdI6uhasewMXAyTZFvgD4PVV9VXgrcAuk4knSVKj6yjBnwKbtY8fA/yoqs5ql28CNhl3sK5OOOGEab205rH33ntPO4JG2HTTTacdQVqQrgXrS8BBSW4FDgROGti2I/DDcQeTJGlQ1y7BlwO/BE4Afgy8ZmDb84EvjzeWJEnr6nod1hU0XYGzeTxNt6AkSROzoDsOz6aqbhhHEEmSRrnDFw5LkrSYLFiSpF6wYEmSesGCJUnqhQUPukhyV+Buw+ur6saxJJIkaRZdJ7/dLMnRSa6kGcL+01m+JEmamK4trGOBRwHvBi4Ebp5UIElSa80aWL16uhlWrZru6w/oWrD2AParqn+bZBhJkubSddDFZYDnqCRJU9O1YL0SeG2S7SYZRpKkuXSdS/DUJI8FLkxyKc0EuMP7rBhvNEmS1upUsJIcAbwM+BoOupAkTUHXQRcvBl5TVYdNMowkSXPpeg7rRuCcSQaRJGmUrgXrn4BVSTLJMJIkzaVrl+CWwMOAC5Kczu0HXVRVvWqMuSRJWkfXgvV04NfAbwCPm2V7ARYsSdLEdB3WvsOkg0iSNIq3F5Ek9ULngpXkfknemeRbSa5o/z06yf0mGVCSJOh+4fAfAV+kubXIJ4CrgXsCfwY8J8mjq+rrE0spSVr2ug66OAI4F3jC4I0ak2wCnNpuf8z440mS1OjaJbgCOHz4rsLt8hE0Q94lSZqYrgXrF8AWc2zbnKarUJKkielasD4JvCnJIwdXtsuHAaeMO5gkSYO6nsP6G+Ak4EtJ1tAMuti6/ToDeMVk4kmS1Oh64fB1wCOT7Ak8FNgGuBI4s6o+O8F8kiQBHQpWko2AbwJ/XVWfBj498VSSJA2Z9xxWVd0E3AO4deJpJEmaQ9dBFx8E9p1kEEmSRuk66OIy4JlJzqa5UPhqmhnaZ1RVvXPc4SRJmtG1YB3Z/rsN8JBZthdgwZIkTUzXUYLO6i5Jmqo5C1GSi5P8Qfv4dUnuvXixJEla16iW028Dm7SPXw/cZ/JxJEma3aguwYuAv0qyBRBg1yRbzrVzVZ067nCSJM0YVbBeDRwHPJtmUMVbR+xbwAZjzCVJ0jrmLFhVdVKSzYF70wxrfxrNPbEkSVp0I0cJVtWtwOVJ9gX+s51TUJKkRdd1WPtxkw4iSdIoXl8lSeoFC5YkqRcsWJKkXrBgSZJ6oVPBSrJTkl0GljdO8sYkH0/y0snFkySp0bWFdTSwcmD5COAAYCPgzUn+dlyBkvxGku3GdTxJ0vqha8F6EPBf0BQU4LnAy6pqT+Bg4IVdDpJk/yQXJflpkjOTPG+W3R4CXNIxlyRpmehasO4O3NA+3qVd/o92+evAfec7QJK9gX8GvgocCvwPcGySjyXZeCGhJUnLT9eCdTFNoQLYCzh3YNaLLYGfdjjGgcARVfWcqjqiqvYC/gR4JPDFdpLdTpKsSnJ2krO//OUvd32aJKnHuhaso4B/SPI14K+Btw9s2x34Zodj/D6wzozuVXUaTSHcDPivJL/TJUxVra6qnatq5912263LUyRJPdepYFXVe4DHAicAj6+q9w9svh54W4fD/ISmNTZ87EuBXYFrgTOAh3bJJElaXjrNJQhQVV8Gbtf/VlWHdDzEOcBTgY/NcowfJdmj3fZ2mtuVSJJ0m84FCyDJjjR3Ht5oeFuHGzh+AHh5ks2r6vpZnv+LJH8KvBN43EJySZLWf50KVpIHAB8GHkBz9+Fh897Asao+Cnx0nn1uAVZ1ySRJWl66trCOAe5GcxPH84GbJ5ZIkqRZdC1YDwb2rqpPTDKMJElz6Tqs/SJmOW8lSdJi6VqwXgEcnOR+kwwjSdJcunYJHgb8NvDdJJcCPx7eoapWjC+WJEnr6lqwvt1+SZI0FZ0KVlXtO+kgkiSN4h2HJUm9MGcLK8nhwNur6vL28UhV9cqxJpMkacCoLsFnAB8ELm8fj1KABUuSNDFzFqyq2mG2x5IkTYPnsCRJvdB5tvYkAR4B7Mjss7UfPcZckiSto+ts7fcETqOZrb1YO2P74H2rLFiSpInp2iV4JM0dg7elKVYPA7YH/g74Pk2rS5KkienaJfgo4ADgynY5VXUZ8MYkd6FpXT1+AvkkSQK6t7DuAaypqluBG4CtB7adAew65lySJK2ja8G6BNimffwd4DkD21YCt7vlvSRJ49S1S/BU4E+AjwD/AJyU5HLgV8B2wKsmE0+SpEbXyW8PGnj8qSS7AnsBGwOfq6pPTSifJElAh4KVZEPgQOATVfUNgKo6Gzh7wtkkSbrNvOewquqXwGtoBl5IkjQVXQddnAn80SSDSJI0StdBF68EPpTkZpoBGFez7iwXVNWNY84mSdJtuhasM9t/3w780xz7bHDn40iSNLuuBeuFDLWoJElaTF2HtR874RySJI3UadBFkouT/MEc2x6U5OLxxpIkaV1dRwluD2w4x7ZNgPuMJY0kSXOYs0swyaase+3VvZJsN7TbRsDewBXjjyZJ0lqjzmG9HHg9zWCLAk6cY78ArxhzLkmS1jGqYH2IZvqlACfTTM90wdA+NwMXtPfGkiRpYuYsWFX1fZq7CZPk0cA5VfWzxQomSdKgroMu/hvYamYhjVVJ3pZk5WSiSZK0VteCdSzNOa0ZhwJHA3sCJybZZ7yxJElaV9eC9RDgCwBJ7gL8JXBwVd0f+EfgZRNJJ0lSq+vUTJsB17WP/wjYHPhgu/wFHCUoSeO31VawatW0UywZXVtYlwMPaB8/CfhuVc1ce7UZcNO4g0mSNKhrC+u9wOFJHktTsF49sG0XmkEZkiRNTNfJbw9LcgXwUOClNAVsxubAv04gmyRJt+nawqKqjgeOn2X9X4w1kSRJs+h6DkuSpKmyYEmSesGCJUnqBQuWJKkXLFiSpF7oNEowyfNHbL4VuAH4RlX9YCypFmC//fZb7JdUR1dfffW0I2iEAw44YNoRNMKmm2467QhLTtdh7cfS3MQRmvtjzRhcV0k+ATzH25BIksZtIZPfXggcBOwEbNn++2rgIuCxwAuAPwbePP6YkqTlrmsL60jg6Kr6p4F119NM13Qz8PqqelSSe9JMhLv/mHNKkpa5ri2shzP3fIH/TTNlE8A5wBZ3NpQkScMWMlv7PnNs27fdDvBbrL0NiSRJY9O1S/A1wL8leRBwCrAG2ApYSXPbkb3b/R4H/Oe4Q0qS1HW29o8muQR4FfDnwL2Aq4CvAftW1Tntfv93UkElScvbQmZrPxt4xgSzSJI0J2e6kCT1QucWVpKnA08D7gNsNLy9qlaMMZckSevoOjXTIcDrgG8A5wM3TzCTJAlgzRpYvXraKRbXqlVzburawnoR8KaqOngsgSRJWqCu57B+EzhtkkEkSRqla8E6AdhzkkEkSRqla5fgacCbk2wJfA748fAOVXXqGHNJkrSOrgXrw+2/29PMyj6sgA3GEUiSpNl0LVg7TDSFJEnz6Do106LfSViSpEFzFqwkm1TVjTOP5zvQzL6SJE3CqBbWT5M8vKrOAn5Gc55qFM9hSZImZlTBeiFw0cDj+QqWJEkTM2fBqqrjBh4fuyhpJEmaQ6cLh5Mcl+QJSTpPlitJ0jh1nenigcAngauSvDvJY5JkgrkkSVpHp4JVVTsDvwMcCewMfB64Msm/JPnjCeaTJAlYwA0cq+qSqjqsqh4M7AQcDewOnJ7ksgnlkyQJuIN3HK6qC4D3Au8DrgJ+e5yhJEkatqBBFEnuBTwDeBawC/Aj4OM0s7lLkjQxXe84/Bc0ReqRwM+Bk4A3Ap+tql9PLp4kSY2uLawjgE/QtK4+VVW/nFwkSZJur2vB2tq5AiVJ09R1WLvFSpI0VaNma78GeHxVnZtkDfPMJVhVW487nCRJM0Z1Cb4DuHrgsZPfSpKmZtTkt4cOPD5kUdJIkjSHrpPfHppkp0mHkSRpLl1nutgP+HaSbyU5OMnvTDKUJEnDuhasewOPA84AXgZ8L8nZSV6RZLtJhZMkaUbXYe23VtUXqmo/YBvgicA3gdcAlyT5ygQzSpK08Mlvq+qWqvoM8JfA/jST3z583MEkSRq00MlvfwPYk2ZewZXAxsCXgNeNP5okSWt1nfx2pkg9FdgU+ArwauCjVbVmYukkSWp1bWGdCpwFHAp8pKr+Z3KRJEm6vXkLVtsNuBtwSVVdMflIkiTdXpdBF7cApwH3n3AWSZLmNG/Bqqpbge8D95x8HEmSZtd1WPtrgNcl+d+TDCNJ0ly6Drp4LbAFcF6SK2hmcV9n9vaqWjHmbJIk3aZrwfp2+yVJ0lR0KlhVte+kgwxKshGwdVVdtpivK0lauhY8NdMieRJwybRDSJKWjq4zXXxkvn2q6pl3Po4kSbPreg5rq1nWbQ78PnAdcEGXgyT5wp14vcHjrAJWdTyWJGk90PUc1qNnW59kW+BE4KiOr7cbTXE7f579Nponz2pgdZuhRu0rSVo/LGi29mFV9cMkhwGHA6d0eMq3gQuq6lmjdkrydODDdyabJGn9Mo5BF7cA9+m475nALh32KyB3OJEkab3TddDFA2ZZfTdgJ+ANwNc6vt7hwCc77HcqsEPHY0qSloGFXDg827mi0BSrF3c5SFVdBFzUYb9fAD/omE2StAx0LVizDbq4CbjcW45IkhZD11GCX5p0EEmSRlnwKMEkmwAvork/1lXA8VVl950kaaLmLFhJjgRWVtWOA+t+k+ac1e8BPwI2A16RZEVVfW/SYSVJy9eoYe2PBj4wtO5AYEfgJVW1JXBv4FLg7yaSTpKk1qiCtT1wztC6PwPOr6r3AlTVGuBI4BETSSdJUmtUwborzUhAAJJsTnPd1fB8gJcC9xp7MkmSBowqWN8Ddh9YfnL772eG9tsauH6MmSRJup1RowT/BXh3ks2Aq4G/prlH1WeH9vsTvBuxJGnC5ixYVXVskm2A/YF7AF8H9q+qX83sk2Qr4CnAoRPOKUla5kZeh1VVhwGHjdi+Bs9fSZIWwThma5ckaeIsWJKkXrBgSZJ6wYIlSeoFC5YkqRcsWJKkXrBgSZJ6wYIlSeoFC5YkqRcsWJKkXrBgSZJ6wYIlSeoFC5YkqRcsWJKkXrBgSZJ6wYIlSeoFC5YkqRcsWJKkXrBgSZJ6wYIlSeoFC5YkqRcsWJKkXrBgSZJ6wYIlSeoFC5YkqRfuOu0AkqQ5bLUVrFo17RRLhi0sSVIvWLAkSb1gwZIk9YIFS5LUCxYsSVIvWLAkSb1gwZIk9YIFS5LUCxYsSVIvWLAkSb2Qqpp2Bg1IsqqqVk87h2bn57N0rY+fTZJPV9Wec2xbWVWnzPU4yUqAwcczy/Mdb/g4w/8OPm9cxxt+PNsxbWEtPU4ctrT5+Sxd691nM1exaq2c5/HKWR6vU2xGHG/4OLMdd9zHG358O05+K0n9dMo8j+fa3uV4sx1ntuOO+3gjj22X4BKT5Oyq2nnaOTQ7P5+ly89m/WeX4NKzXvXBr4f8fJYuP5v1nC0sSVIv2MKSJPWCBUuS1AsWrCUgye8mOSbJN5LckuT0aWdSI8kzkpyc5IokP0tyTpJnTzuXGkmenuSMJNcluSnJBUlem+Ru086m8XNY+9LwQOCJwFcB/6MtLX8DXAK8HLiW5nP6UJItq+qfp5pMAFsAXwTeAvwYWAEcAtwL+KuppdJEOOhiCUhyl6q6tX38MWDLqtp9uqkE0Bama4fWfQh4eFXtMKVYGiHJPwL7A79V/oJbr9gluATMFCstPcPFqnUusPViZ1Fn12FPxXrJLkFp4XYFzp92CK2VZANgQ+AhwF8D77R1tf6xYEkLkGQP4CnAC6edRev4OU3BAjge+NspZtGE2CUodZRke+BDwElVdex002jIrsAfA6+g+YPiX6YbR5NgC0vqIMnmwKeAy4DnTjmOhlTV19uHX0lyLXBckiOr6qJp5tJ42cKS5pFkE+ATNCfyn1RVP59yJI02U7wcxbmesYUljZDkrsBHgd8DHlFV10w5kub3iPbfS6aaQmNnwVoC2r/gn9gu/jawaZKnt8unVtWN00km4Giaz+YAYPMkuwxsO7eqfjmdWILmjrzA54HvALfQFKtXAB+2O3D944XDS0B7Mn+uvwZ3qKpLFy+NBiW5FLjvHJv9bKYsyRuAvYDtgV8DFwPvA95VVb+aYjRNgAVLktQLDrqQJPWCBUuS1AsWLElSL1iwJEm9YMGSJPWCBUuS1AsWrAlKsk97S/WfJvlRknOTvHWRMxzSzq22kOfsk6SS/K9J5VqIJDu238c9xnzcS5McMc5jjnitpyT57yQ3t9d2aRZ35Of1Dr7O7u3P+IMm/VoaHwvWhCR5NfCvwGeApwHPB04C/nSRo/wr8PgFPueTwMOBpTLDxo7A64F7jPm4ewFvH/Mxb6e9V9PxwDeAx7Svq9ndkZ9XLRNeODwhSa4APl5V+w+tz529sVz7C3CDqrr5zhynL5I8GTiFDjNLJNmoqm5alGAdJbkP8EPgKVV18oj9Nq6qXyxesuUrye7AF4H/XVXfnm4adWULa3LuAVw1vHKwWM3VLZHk9CQfG1g+NsnZSZ6a5DvATcDD2uc+cei5GyS5qp2yZp0uliR3T/LzJP93OFd7/Pe3j9fpEkyyfbv8zCTHJPlJksuTHJrkLkPHeUaS7yf5RZIvJnlw+9x9Rr1ZSV6d5MIkNyW5Osmnk9yr/cVySrvbJe2xLh3KuaJ9z35Be+O+JG9K8q0kP2uzfjDJvYZec50uwYH3+XFJvtm+V19J8sCh570oyXfa7/HaJF8a3mdg331oihXASW3eQ9ptleRvkrwtyRrgW+36LZMcl+S6JDe239vOs2VPclCSK9vP5Mg0ntjm+2mSjyf5rXne+/sk+UiSa9rv6aKZn5+h9+WpSb7bfkZfSfKAoePcpc1zYZJfJvlekhfM8np7JTmrfa3rkpya5L7tttt1CSbZvP25u7p97TOSPOyOfibqLwvW5HwdeGmSFyTZYgzH2x44HDiMZjLWS4CzgGcN7fco4J7Ah4cP0N4W4xPDz0lyP+CPZnvOkMOBnwFPBz4AvK59PHOcnYETaL73vYCTOxyTJM8HDgbeStMd9JfAhcDd22Md2O76NJquyuEutX9rv68ntv8CbA28EXgS8DLgfsAX0rROR9kOeAvwj8Cz2+N8JEnarLsB72q//yfQ3Hn4DGCzOY73yTY37ffxcJpurxl/C2wDPI/m1u4AH6d5Hw6k+azuAnwxye8OHXtvYAWwL81n8zc07+EbgL8D/oLm5+Gweb7n44FtgVXt9/SPrL1774z7Dhz7z9vv9zNJNhrY55+B1wKrad73E4H3pmkhA5DkecB/ABcBz2yzfw/YarZgSTakmdz2cTTv1VOBNcDnZ/4AuQOfifqqqvyawBfwf2gm4izgVprZpP8e2HRgn93b7Q8aeu7pwMcGlo9t9/vDof1eDvwE2HBg3THAdwaWDwGuHVjei2ZW63sPrHs1cD1wt3Z5n/b1/le7vH27fPzQ658HnDCw/FHg27Rdze26V7bP3WfEe/UvwL+P2P7k9hjbD62fyXnAPJ/FBjSz4Bew28D6S4Ejht7nXwO/N7Duqe3z7t8uHwics8CfhZn378lD64tmxvfBdXu26x81sO7uNL+kjxnKfiFN1/DMurPa/DsMrDscuHqefD8DVo7YPvPzt+vAuvu2r/UX7fLv0vycv2DouccDX2sf3wW4AviPEa81/PP6IuDmoc/krjQF7y134jPZnVn+7/m1tL9sYU1IVX0T2IlmkMXRQGj+6j07d2z03RVVdd7Quo8Av0nzS27m3k1Po2nlzOVTNL+gnjGw7lnAiTX/ObHPDi2fD9xnYPmhwCnV/kZozXnOZsB5wBPTdDGu6NAKGvbJ4RVJntB2Hf2E5hfr5e2mHec51qVV9f2B5fPbf2e+z/OAByc5KsluSe62wKzDhrOvANZU1ZdmVtTalvEjh/Y9vapuGVi+sM1/ydC6rebJeR5wWJou1u3m2OeaqjpjINMPgHPavAB70BSsE5PcdeYLOA34w/Yz/X3g3jSzqXf12PZ1Lhk4JsCXgJlu0vMY72eiJcqCNUFV9cuqOqWq/qqqHgC8mOZGgC+6A4e7epbjXwF8hbVdfHsAWzKiYFUzIOGkmeck+X3gD0Y9Z8CPh5ZvBga7hO5F0xIYNLw8m/fSdAk+EzgTuDrJGxZQuNZ5b5I8lKZQXk7T1fZwYOY+Vhsx2o+HlmeK+EYAVfV5mm6s3WhawtcmOTrJ3TtmHZmdpnvwdp91u27zDllnWxeauyXP5VnA2cBRwA+SnJdkj6F9Zrtx5TVtXmh+7jagafH/auDrWJoW0TbATNf4lSOyDNuS5rP71dDXvjTdmJP4TLREeQPHRVRV70lyOHD/dtXMaLbhXyabA8PXosw1svDDwJuSbEzzi+fcoRbCXM85pf1r+lk0ReULHb6F+VzF7c9FzHpuYlBV3Urzy/KoJNsCz6E5j3IFzbmJeQ8xtLwXzff0rJnW3sxJ/XGoquOA45JsRdOiPQq4ATjojhxuaPlKmvNmw+5J0207du0fPvukGUCzgqZb7uQk21XVde1us2Xamqarmzbbr2luoHjrLPteQ9MbAGuLXBfX0xTTv5xl2203zxzzZ6IlyhbWhCS53X/w9j/TZqz9C3qmm2qngX22pek66eqjwMY0v6T3oltL6bPAj2haNM+iOV92y+indPI1YOXMAIXWgq47q6ofVtWbaLqyZkahrdPK6WBj4FdDXZPPWUiOLqpqTVUdA/wna7PeWWcCW7cDCYDb7kj9JJrW9MRU1a1V9VXgUGAT1r1x5dZJdh3ItB3wEJrzZtD8wbMBsFlVnT3L183ABTR/hNxu5OAIp9GcH7tslmN+a5bvYRKfiZYIW1iT860kJ9EUh2to/vMfSHMx7nEAVXV5kq8Bb0hyI80fEAezgL+kq+qaJKcDR9AMpf9Ih+f8KsmJNKPKtgFuN8z9DnozzS/cE5K8j6YQv6TdNttf3QAkOYbme/4qTZfSo2m6Tl/V7nJB++9+SU4Abpztl9WAzwEvS/I2miHxuwLPvSPf0CxZD6VpAZ9O0wp+MM1IvLH8JV9Vn0ny/4APJzkIuI7m52ZjmtGLY5VkM5qL24+nGa23Ic0t5q8C/ntg12uB9yf5O+AXNAOIrqHp8qOqLkjyLprP/nCaVtFGwAOBHavqxVV1a5JXAh9M8kGa0Z1FczH1v1XV2bNEPJ5mtOPpaS5BuJima3EFcFVVHTXpz0RLhwVrcv4eeArNTAqb0/wCOIOmm2rwpPif0wxz/gBNi+uVNKP/FuIE4N3AV6v7LdtPoDmX9j80f43eaVV1dpJn0wwnfwpru3I+R9M9M5f/oils+9H8krsQeElVfbw97g+SHEgz7PulNO/T9iNynJrkVe2+L2mP/2SaX8h31tdoPp+9abq4fkDThfZPYzj2jL2AI4G30bwfZwGPqaoLx/gaM26iuf7rAJpzQjfS/OHwJ7XuRcw/oPlc30Tzx9fZwLNr3Yu096d5j19C8/N/A82glffM7FBVH0pyE/Aa4GPAz9vXm/VcZ1XdlOTR7fEOpekavYbmPZkZ0LMYn4mWAGe60EQleS7wfuB+Q4VaPZHkWJrh3zvPt680SbawNFZJ3knTovoRzTmO1wKftFhJurMsWBq3LWiuO9uC5vzLh2m6OSXpTrFLUJLUCw5rlyT1ggVLktQLFixJUi9YsCRJvWDBkiT1wv8HY6uQ9MRAt/MAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import os.path\n",
    "\n",
    "runs =500\n",
    "all_survivors = []\n",
    "all_params = []\n",
    "for index in range(runs):\n",
    "    filename = 'data/2_res_new_results_gK_separated_'+str(index)+'.pkl'\n",
    "    if(os.path.exists(filename)):\n",
    "        openfile = open(filename, 'rb')\n",
    "        survivor_pool, params_pool = pickle.load(openfile)\n",
    "        all_survivors.append(survivor_pool)\n",
    "        all_params.append(params_pool)\n",
    "        \n",
    "heatmap = np.zeros([3, 3])\n",
    "for idx, i in enumerate(all_survivors):\n",
    "    bug1 = sum([j<50 for j in i])\n",
    "    bug2 = sum([j>=50 for j in i])\n",
    "    heatmap[3-bug1, bug2-1]+=1\n",
    "    \n",
    "\n",
    "import pandas as pd\n",
    "xrange = 3\n",
    "yrange = 3\n",
    "Isurvivor = np.repeat(range(1, xrange+1), yrange)\n",
    "IIsurvivor = np.tile(range(1, yrange+1), xrange)\n",
    "df = pd.DataFrame({\"I survivor\": Isurvivor,\n",
    "                  \"II survivor\": IIsurvivor,\n",
    "                  \"percentage\": np.array([heatmap[Isurvivor[i]-1, IIsurvivor[i]-1] for i in range(9)]) })\n",
    "\n",
    "g = sns.jointplot(data=df, x='I survivor', y='II survivor', kind='hist', bins=(xrange, yrange))\n",
    "g.ax_marg_y.cla()\n",
    "g.ax_marg_x.cla()\n",
    "sns.heatmap(data=np.log(df['percentage'].to_numpy().reshape(xrange, yrange)+1), ax=g.ax_joint, cbar=False, cmap='Greys')\n",
    "# g.ax_joint.hlines(y=np.arange(0, 3), xmin=np.full(3, 0), xmax=np.full(3, 3), color=\"black\")\n",
    "# g.ax_joint.vlines(x=np.arange(0, 3), ymin=np.full(3, 0), ymax=np.full(3, 3), color=\"black\")\n",
    "\n",
    "g.ax_marg_y.barh(np.arange(0.5, yrange), df.groupby(['I survivor'])['percentage'].sum().to_numpy()/len(all_survivors), color='r', alpha = 0.4)\n",
    "g.ax_marg_x.bar(np.arange(0.5, xrange), df.groupby(['II survivor'])['percentage'].sum().to_numpy()/len(all_survivors), color='b', alpha = 0.4)\n",
    "\n",
    "g.ax_joint.set_xticks(np.arange(0.5, xrange))\n",
    "g.ax_joint.set_xticklabels([1, 2, 3], fontsize=15)\n",
    "g.ax_joint.set_yticks(np.arange(0.5, yrange))\n",
    "g.ax_joint.set_yticklabels([3, 2, 1], fontsize=15)\n",
    "g.ax_joint.set_xlabel(\"Surviving strains from species I\", fontsize=15)\n",
    "g.ax_joint.set_ylabel(\"Surviving strains from species II\", fontsize=15)\n",
    "\n",
    "g.ax_marg_x.set_yscale(\"log\")\n",
    "g.ax_marg_y.set_xscale(\"log\")\n",
    "g.ax_marg_x.set_ylim(0.001, 1)\n",
    "g.ax_marg_y.set_xlim(0.001, 1)\n",
    "g.ax_marg_x.tick_params(axis='x', bottom=False, labelbottom=False)\n",
    "g.ax_marg_y.tick_params(axis='y', left=False, labelleft=False)\n",
    "g.ax_marg_x.tick_params(axis='y', left=False, labelleft=False)\n",
    "g.ax_marg_y.tick_params(axis='x', bottom=False, labelbottom=False)\n",
    "\n",
    "g.fig.subplots_adjust(hspace=0.05, wspace=0.02)  # less spaced needed when there are no tick labels\n",
    "plt.savefig(\"heatmap_only_gK.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e9fb9d12",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.002, 0.   , 0.   ],\n",
       "       [0.038, 0.002, 0.   ],\n",
       "       [0.892, 0.064, 0.002]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "heatmap/500"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bbccc0b2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab4c947f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "055d5cd6",
   "metadata": {},
   "source": [
    "#### Then we fix the g and K (no satellite niches), and create the strains by slightly fluctuating the consumption ratios of the 2 species. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "ded6df27",
   "metadata": {},
   "outputs": [],
   "source": [
    "def invade(profile, p_profile, p_dilute, p_invade):\n",
    "    # load variables\n",
    "    N, R, gmax, K, x = p_profile\n",
    "    c0, D, t, n, fluc_thres, elim_thres = p_dilute\n",
    "    p_dilute_once = c0, D, t, 2, fluc_thres, elim_thres\n",
    "    gmax_inv, K_inv, x_inv = p_invade\n",
    "    p_profile_temp = N+1, R, np.vstack([gmax, gmax_inv]), np.vstack([K, K_inv]), np.hstack([x, x_inv])\n",
    "    profile_temp = Dilute(p_profile_temp, p_dilute_once, np.hstack([profile[:-R], [D*elim_thres], profile[-R:]]))\n",
    "    if(profile_temp[-1][-(R+1)]<=0):\n",
    "        return profile, p_profile\n",
    "    else:\n",
    "        profile_temp = Dilute(p_profile_temp, p_dilute, np.hstack([profile[:-R], [5*D*elim_thres], profile[-R:]]))\n",
    "        return profile_temp[-1], p_profile_temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9da2c119",
   "metadata": {},
   "outputs": [],
   "source": [
    "# species A: xA = 0.2(consumes 0.2 resource 1 and 0.8 resource 2)\n",
    "# species B: xB = 0.7(consumes 0.7 resource 1 and 0.3 resource 2)\n",
    "# species A has g mutating around 1 and K around 4\n",
    "# species B has g mutationg around 0.55 and K around 0.01\n",
    "\n",
    "all_survivors = []\n",
    "all_params = []\n",
    "\n",
    "for pool in range(500):\n",
    "\n",
    "    gA1, gA2, KA1, KA2 = 1, 1, 4, 5\n",
    "    gB1, gB2, KB1, KB2 = 0.55, 0.55, 0.05, 0.03\n",
    "\n",
    "    N_pool = 50\n",
    "\n",
    "    g_diff = np.random.normal(0, 0, size = 2*N_pool)\n",
    "    K_exp = np.random.normal(0, 0, size = 4*N_pool)\n",
    "    gA1new = [gA1+i for i in g_diff[0:N_pool]]\n",
    "    gA2new = [gA2+i for i in g_diff[0:N_pool]]\n",
    "    gB1new = [gB1+i for i in g_diff[N_pool:2*N_pool]]\n",
    "    gB2new = [gB2+i for i in g_diff[N_pool:2*N_pool]]\n",
    "    KA1new = [KA1*(10**i) for i in K_exp[0:N_pool]]\n",
    "    KA2new = [KA2*(10**i) for i in K_exp[N_pool:2*N_pool]]\n",
    "    KB1new = [KB1*(10**i) for i in K_exp[2*N_pool:3*N_pool]]\n",
    "    KB2new = [KB2*(10**i) for i in K_exp[3*N_pool:4*N_pool]]\n",
    "    \n",
    "    xA, xB = 0.2, 0.7\n",
    "    sigma = 0.1\n",
    "    xAdist = stats.truncnorm(-0.2, 0.8, loc=0, scale=sigma)\n",
    "    xAnew = np.sort(xA+xAdist.rvs(N_pool))[::-1]\n",
    "    xBdist = stats.truncnorm(-0.7, 0.3, loc=0, scale=sigma)\n",
    "    xBnew = np.sort(xB+xBdist.rvs(N_pool))\n",
    "\n",
    "    N, R = 2*N_pool, 2\n",
    "    gmax = np.array([[gA1new[i], gA2new[i]] for i in range(N_pool)]+[[gB1new[i], gB2new[i]] for i in range(N_pool)])\n",
    "    K = np.array([[KA1new[i], KA2new[i]] for i in range(N_pool)]+[[KB1new[i], KB2new[i]] for i in range(N_pool)])\n",
    "    x = np.array([np.hstack([xAnew, xBnew])])\n",
    "    c0, D, t, n, fluc_thres, elim_thres = [10, 10], 10, 24.0, 5000, 1e-5, 1.0e-6\n",
    "    \n",
    "    start=time.time()\n",
    "    \n",
    "    p_profile = 1, R, gmax[:1, :], K[:1, :], x[:, :1]\n",
    "    p_all = N, R, gmax, K, x\n",
    "    p_dilute = c0, D, t, n, fluc_thres, elim_thres\n",
    "    profile = Dilute(p_profile, p_dilute, np.array([1]+[0 for i in range(R)]))[-1]\n",
    "\n",
    "    for i in range(N):\n",
    "        _, _, _, _, x_profile = p_profile\n",
    "        if(x[0, i] not in x_profile[0]):\n",
    "            p_invade = gmax[i:i+1], K[i:i+1], x[:, i:i+1]\n",
    "            profile_new, p_profile_new = invade(profile, p_profile, p_dilute, p_invade)\n",
    "            N_temp, R_temp, gmax_temp, K_temp, x_temp = p_profile_new\n",
    "            survivor = profile_new[:-R]>0\n",
    "            p_profile = sum(survivor), R_temp, gmax_temp[survivor, :], K_temp[survivor, :], x_temp[:, survivor]\n",
    "            profile = np.hstack([profile_new[:-R][survivor], profile_new[-R:]])\n",
    "    \n",
    "    end=time.time()\n",
    "    print(end-start)\n",
    "    \n",
    "    N_temp, R_temp, gmax_temp, K_temp, x_temp = p_profile\n",
    "    \n",
    "    xlist = list(xAnew)+list(xBnew)\n",
    "    survivors = [xlist.index(i) for i in x_temp[0]]\n",
    "    \n",
    "    all_survivors.append(survivors)\n",
    "    all_params.append([xAnew, xBnew])\n",
    "    pickle.dump((all_survivors, all_params),open(\"2_res_new_results_x_separated.pkl\", \"wb\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "1c314663",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAawAAAGyCAYAAAC4FPBHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAlIElEQVR4nO3deZhkZX238fsrRhYTMWyKCoIxRNSYqDgibrjjQhQVl7iBy6ghxgWiiEYlRlEEUaMoYzSAS3CJiCiuKCa+KDgIbijKJkLYcUERUfi9f5zTTE1RXX0aqrr7TN+f6+pr6ix16tdVPf3t5znPeU6qCkmSlrqbLXYBkiR1YWBJknrBwJIk9YKBJUnqBQNLktQLN1/sAiag8zDHVaumWcbysXLlYlcgLQtZ7AKWmt62sJLsmmTVscceu9ilSJIWQG8DS5K0vPS2S7CqjgWOBV6w2LVIkqavty0suwQlaXnpbWBJkpYXA0uS1Auew5Ik9UJvW1iew5Kk5cUWliSpF3rbwpIkLS+9DSy7BCVpebFLUJLUC71tYUmSlhcDS5LUC70NLM9hSdLy4jksSVIv9LaFJUlaXgwsSVIvGFiSpF4wsCRJvWBgSZJ6obeB5bB2SVpeehtYkqTlxcCSJPWCgSVJ6gVnupAk9UJvW1gOupCk5cUWliSpF2xhSZJ6wRaWJKkXbGFJknrBFpYkqRd628KSJC0vBpYkqRd6G1iew5Kk5cVzWJKkXuhtC0uStLwYWJKkXjCwJEm9YGBJknqht4MutG5ZtWqxK1h3rFy52BVI09HbFpbD2iVpeeltYEmSlhcDS5LUCwaWJKkXejvowpkuJGl56W0Ly0EXkrS82MKSJPWCLSxJUi/YwpIk9YItLElSL9jCkiT1Qm9bWJKk5cXAkiT1Qm8Dy3NYkrS8eA5LktQLvW1hSZKWFwNLktQLBpYkqRcMLElSL6SqFruGmyTJF4DNFrsOSZqwy6pql8UuYinpfWBJkpYHuwQlSb3Q2+uwBthElLQuyi677FJf+MIXFruOxZBRK21hSdISddllly12CUuKgSVJ6gUDS5LUCwaWJKkXDCxJUi8YWJKkXjCwJEm9YGBJknrBwJIk9YKBJUnqBQNLktQLBpYkqRcMLElSLxhYkqReMLAkSb1gYEmSesHAkiT1goElSeoFA0uS1AsGliSpFwwsSVIvGFiSpF4wsCRJvWBgSZJ6wcCSJPWCgSVJ6gUDS5LUCwaWJKkXDCxJUi8YWJKkXjCwJEm9YGBJknrBwJIk9YKBJUnqBQNLktQLBpYkqRcMLElSLxhYkqReMLAkSb1w89k2JDlwHsepqnrVBOqRJGmkWQML2H0exynAwJIkTc2sgVVV2y5kIZIkjeM5LElSL4w7h3XX+Ryoqk6/6eVIkjTauHNYP6A5NzWXtPutN5GKJEkaYVxgPWTBqpAkaQ7jBl18fSELkSRpHAddSJJ6wcCSJPXCuHNYCybJtsC2wDlVdc5i1yNJWnoWrIWV5E+SvCvJFUl+MzP1U5L3AGcCXwHOTPKRJI44lCStZSFbWK8Eng+8HbgCeEmSzYEnAHsA3wEeABwEvBA4dAFrkyQtcZ1aWEm2T7LjwPKGSd6c5NNJXtLxtZ4FvL6qXltVb2+XnwPsX1UfqqofVtVhNIH23DnqWZlkdZLVq1at6vjykqQ+69rCOhQ4EfhWu3wQTavof4G3Jtmgqt42xzHuCJw8sHxK++/JQ/t9AxgbglW1CphJqi4XN0uSeq7rOay7A9+E5lwU8EzgZVW1C7Afc7SIWr8Fbj2w/Pv266qh/dZjiQwGkSQtHV0D65bAr9vHO7bLn2qXv0PTeprLGcAOMwtVdV1VbVhVpw3tdzfg3I51SZKWia6BdTZNUAHsBpxaVZe3y5sBV3Y4xtuBn3fY7+HA0R3rkiQtE1273g4B3ptkd+CewJ4D23YGvjfXAarqv7u8UFU9umNNkqRlpFNgVdUHkvwUuA+wb1UdP7D5CuAdU6hNkqTrdR7cUFX/A/zPiPVvmGRBkiSN0nmmiyRbJHlrkuOT/CTJ3dr1L01yv+mVKElS9wuHV9BMn/QkmhF8fwGs327eEth7GsVJkjSjawvrEOCrwHY00yZlYNvJwIoJ1yVJ0lq6nsO6F/D4qrouSYa2XQ5sMdmyJElaW9cW1q+AzWfZdifg4smUI0nSaF0D6xhg/yR3GlhXSTYD9mHNrBeSJE1F18Dal2ZqptNZM7T9fTTTLf0OeN3kS5MkaY2uFw7/or29yLOAh9FMZHsF8B/AkVX1++mVKEnS/C4cvgb4QPslSdKC6nzhsCRJi2nWFlaSS4BHVdWpSS5ljhslVpVD2yVJUzOuS/A9rBmu/h68s68kaRHNGlhVtf/A4zcsSDWSJM2i61yCWyW51yzb7pVkq8mWJUnS2roOungv8MxZtv09cOhkypEkabSugbUjzeS3o3yt3S5J0tR0DayNGD/o4pYTqEWSpFl1DazvA0+fZdvTgR9OphxJkkbrOtPFW4D/TrI+cDhwIc2NG59Dc1PHJ02lOkmSWl3nEjw6yXOAA2jCqWhu4ngB8Myq+vTUKpQkifnNJfihJB8G7gJsQnPjxjOqyguKJUlT1zmwAKqqkvyYpjvwEsNKkrRQOk9+m+QxSU4CrgbOA+7Rrl+VZLZrtCRJmoiuM108G/gM8GNg5dDzfgo8b/KlSZK0RtcW1muAt1XVc4APD237IXDXiVYlSdKQroF1R+DLs2y7GrjVZMqRJGm0roH1c+Ces2zbAThzMuVIkjRa18D6APD6dnDFhu26JHkY8Erg/dMoTpKkGV2Htb8V2Ao4Ari2XXcisB5wWFW9awq1SZJ0va4zXRSwV5JDgIcBmwJXAF+tqp9MsT5JkoD5Xzh8Jp6vkiQtgs6BleQWwB7ACpqZLi4ETgKOqKprplKdJEmtrhcOb09zgfB7gLvTnMe6e7t8ZhKvw5IkTVXXFtYq4FfAA6vqvJmVSbYGPge8D3jQ5MuTJKnRdVj7DsDrBsMKoF1+HXCfSRcmSdKgroF1LrDBLNs2oJkMV5KkqekaWPsC/5bkvoMrk+wI/CvwqkkXJknSoK7nsF5LM1/giUkuAS4Btmi/Lgf2S7LfzM5VtWLShUqSlreugfWD9kuSpEXRdaaLPaddiCRJ43S+4/CwJH+e5G+TrD/JgiRJGqXrhcP7J3nLwPJDaUYGngKcneRuU6pPkiSgewvrGcCPB5YPBr4B3L9df8CE65IkaS1dA+t2wNkASbYC/gZ4fVV9C3g7sON0ypMkqdE1sK4ENm4fPxT4RVWd3C5fDWw06cIkSRrUdVj714F9k1wH7AMcM7BtO+Dnky5MkqRBXVtYLwd+DxwF/BJ4zcC2ZwP/M9myJElaW9frsC6g6Qoc5VE03YKSJE3NvO44PEpV/XoShUiSNM6NvnBYkqSFZGBJknrBwJIk9YKBJUnqhXkPukhyc+AWw+ur6qqJVCRJ0ghdJ7/dOMmhSS6kGcJ+5YgvSZKmpmsL63DgwcD7gTOBa6ZVkCRJo3QNrIcBL6yq/5pmMZIkzabroIvzAM9RSZIWTdfAeiXw2iRbT7MYSZJm03UuweOSPBw4M8m5NBPgDu+zYrKlSZK0RqfASnIQ8DLg2zjoQpK0CLoOung+8JqqOmCaxUiSNJuu57CuAk6ZZiGSJI3TNbDeCaxMkmkWI0nSbLp2CW4G3Bc4I8kJ3HDQRVXVqyZYlyRJa+kaWE8G/gj8CfCIEdsLMLAkSVPTdVj7ttMuRJKkcby9iCSpFzoHVpI7JXlvku8nuaD999Akd5pmgZIkQfcLh+8NfI3m1iKfBS4GbgM8CXhGkodU1XemVqUkadnrOujiIOBU4NGDN2pMshFwXLv9oZMvT5KkRtcuwRXAgcN3FW6XD6IZ8i5J0tR0DazfAZvOsm0Tmq5CSZKmpmtgfQ54S5IHDK5slw8Ajp10YZIkDep6DusVwDHA15NcSjPoYov260Rg7+mUJ0lSo+uFw5cDD0iyC3AfYEvgQuCkqvrSFOuTJAnoEFhJNgC+B/xTVX0B+MLUq5Ikacic57Cq6mrg1sB1U69GkqRZdB108RFgz2kWIknSOF0HXZwHPCXJapoLhS+mmaF9RlXVeyddnCRJM7oG1sHtv1sC9xqxvQADS5I0NV1HCTqruyRpUc0aREnOTvI37ePXJbndwpUlSdLaxrWcbg9s1D5+PXCH6ZcjSdJo47oEzwL+McmmQICdkmw2285Vddyki5Mkaca4wHo1cATwdJpBFW8fs28B602wLkmS1jJrYFXVMUk2AW5HM6z9iTT3xJIkacGNHSVYVdcB5yfZE/jfdk5BSZIWXNdh7UdMuxBJksbx+ipJUi8YWJKkXjCwJEm9YGBJknqhU2Al2T7JjgPLGyZ5c5JPJ3nJ9MqTJKnRtYV1KLDrwPJBwEuBDYC3JvnnSRWU5E+SbD2p40mS1g1dA+vuwDehCRTgmcDLqmoXYD/guV0OkmSvJGcluTLJSUmeNWK3ewHndKxLkrRMdA2sWwK/bh/v2C5/ql3+DnDHuQ6Q5GnAvwPfAvYH/g84PMknk2w4n6IlSctP18A6myaoAHYDTh2Y9WIz4MoOx9gHOKiqnlFVB1XVbsAjgQcAX2sn2e0kycokq5OsXrVqVdenSZJ6rOsdhw8B3ptkd+CewJ4D23YGvtfhGH9FE1rXq6rj28Ecnwe+meTRXYqpqlXATFJVl+dIkvqtUwurqj4APBw4CnhUVX1oYPMVwDs6HOZXNK2x4WOfC+wEXAacCNynS02SpOUlVQvTQElyDHBlVT1zlu0bAp8EHg1UVXW9XYktLEnrouywww61evXqxa5jMWTUyq5dgs0Rku1o7jy8wfC2Djdw/DDw8iSbVNUVI57/uyR/B7wXeMR86pIkrfs6BVaSuwIfA+7K6OSb8waOVfUJ4BNz7HMtsLJLTZKk5aVrC+sw4BY0N3E8HbhmahVJkjRC18C6J/C0qvrsNIuRJGk2Xa/DOosR560kSVooXQNrb2C/JHeaZjGSJM2ma5fgAcDtgR8nORf45fAOVbVicmVJkrS2roH1g/ZLkqRF0SmwqmrPufeSJGl6vOOwJKkXZm1hJTkQeFdVnd8+HquqXjnRyiRJGjCuS3B34CPA+e3jcQowsCRJUzNrYFXVtqMeS5K0GDyHJUnqhc6ztScJcH9gO0bP1n7oBOuSJGktXWdrvw1wPM1s7cWaGdsH70VlYEmSpqZrl+DBNHcM3oomrO4LbAP8C/BTmlaXJElT07VL8MHAS4EL2+VU1XnAm5PcjKZ19agp1CdJEtC9hXVr4NKqug74NbDFwLYTgZ0mXJckSWvpGljnAFu2j38IPGNg267ADW55L0nSJHXtEjwOeCTwceDfgGOSnA/8AdgaeNV0ypMkqdF18tt9Bx5/PslOwG7AhsCXq+rzU6pPkiSgQ2AlWR/YB/hsVX0XoKpWA6unXJskSdeb8xxWVf0eeA3NwAtJkhZF10EXJwH3nmYhkiSN03XQxSuBjya5hmYAxsWsPcsFVXXVhGuTJOl6XQPrpPbfdwHvnGWf9W56OZIkjdY1sJ7LUItKkqSF1HVY++FTrkOSpLE6DbpIcnaSv5ll292TnD3ZsiRJWlvXUYLbAOvPsm0j4A4TqUaSpFnM2iWY5Fasfe3VbZNsPbTbBsDTgAsmX5okSWuMO4f1cuD1NIMtCjh6lv0C7D3huiRJWsu4wPoozfRLAT5DMz3TGUP7XAOc0d4bS5KkqZk1sKrqpzR3EybJQ4BTquo3C1WYJEmDug66+BGw+cxCGiuTvCPJrtMpTZKkNboG1uE057Rm7A8cCuwCHJ1kj8mWJUnS2roG1r2ArwIkuRnwYmC/qroL8CbgZVOpTpKkVtfA2hi4vH18b2AT4CPt8leBO0+4LkmS1tI1sM4H7to+fizw46qaufZqY+DqSRcmSdKgrpPffhA4MMnDaQLr1QPbdqQZlCFJ0tR0nfz2gCQXAPcBXkITYDM2Af5jCrVJknS9ri0squpI4MgR61800YokSRqh6zksSZIWlYElSeoFA0uS1AsGliSpFwwsSVIvdBolmOTZYzZfB/wa+G5V/WwiVc1DkoV+SXX03e9+d7FL0Bj3uMc9FrsEaV66Dms/nOYmjtDcH2vG4LpK8lngGd6GRJI0afOZ/PZMYF9ge2Cz9t9XA2cBDweeAzwQeOvky5QkLXddW1gHA4dW1TsH1l1BM13TNcDrq+rBSW4D7A3sNeE6JUnLXNcW1v2Yfb7AH9FM2QRwCrDpTS1KkqRh85mtfY9Ztu3Zbgf4c9bchkSSpInp2iX4GuC/ktwdOBa4FNgc2JXmtiNPa/d7BPC/ky5SkqSus7V/Isk5wKuAvwduC1wEfBvYs6pOaff7h2kVKkla3uYzW/tqYPcp1iJJ0qyc6UKS1AudW1hJngw8EbgDsMHw9qpaMcG6JElaS9epmd4AvA74LnA6cM0Ua5IkAVx6KaxatdhVLLyVK0eu7trCeh7wlqrab2IFSZI0D13PYf0ZcPw0C5EkaZyugXUUsMs0C5EkaZyuXYLHA29NshnwZeCXwztU1XETrEuSpLV0DayPtf9uQzMr+7AC1ptEQZIkjdI1sLadahWSJM2h69RMC34nYUmSBs0aWEk2qqqrZh7PdaCZfSVJmoZxLawrk9yvqk4GfkNznmocz2FJkqZmXGA9Fzhr4PFcgSVJ0tTMGlhVdcTA48MXpBpJkmbR6cLhJEckeXSSzpPlSpI0SV1nurgb8DngoiTvT/LQJJliXZIkraVTYFXVDsBfAAcDOwBfAS5M8u4kD5xifZIkAfO4gWNVnVNVB1TVPYHtgUOBnYETkpw3pfokSQJu5B2Hq+oM4IPAfwIXAbefZFGSJA2b1yCKJLcFdgeeCuwI/AL4NM1s7pIkTU3XOw6/iCakHgD8FjgGeDPwpar64/TKkySp0bWFdRDwWZrW1eer6vfTK0mSpBvqGlhbOFegJGkxdR3WblhJkhbVuNnaLwEeVVWnJrmUOeYSrKotJl2cJEkzxnUJvge4eOCxk99KkhbNuMlv9x94/IYFqUaSpFl0nfx2/yTbT7sYSZJm03WmixcCP0jy/ST7JfmLaRYlSdKwroF1O+ARwInAy4CfJFmdZO8kW0+rOEmSZnQd1n5dVX21ql4IbAk8Bvge8BrgnCTfmGKNkiTNf/Lbqrq2qr4IvBjYi2by2/tNujBJkgbNd/LbPwF2oZlXcFdgQ+DrwOsmX5okSWt0nfx2JqSeANwK+AbwauATVXXp1KqTJKnVtYV1HHAysD/w8ar6v+mVJEnSDc0ZWG034IOAc6rqgumXJEnSDXUZdHEtcDxwlynXIknSrOYMrKq6DvgpcJvplyNJ0mhdh7W/Bnhdkr+eZjGSJM2m66CL1wKbAqcluYBmFve1Zm+vqhUTrk2SpOt1DawftF+SJC2KToFVVXtOu5BBSTYAtqiq8xbydSVJS9e8p2ZaII8FzlnsIiRJS0fXmS4+Ptc+VfWUm16OJEmjdT2HtfmIdZsAfwVcDpzR5SBJvnoTXm/wOCuBlR2PJUlaB3Q9h/WQUeuTbAUcDRzS8fUeRBNup8+x3wZz1LMKWNXWUOP2lSStG+Y1W/uwqvp5kgOAA4FjOzzlB8AZVfXUcTsleTLwsZtSmyRp3TKJQRfXAnfouO9JwI4d9isgN7oiSdI6p+ugi7uOWH0LYHvgjcC3O77egcDnOux3HLBtx2NKkpaB+Vw4POpcUWjC6vldDlJVZwFnddjvd8DPOtYmSVoGugbWqEEXVwPne8sRSdJC6DpK8OvTLkSSpHHmPUowyUbA82juj3URcGRV2X0nSZqqWQMrycHArlW13cC6P6M5Z/WXwC+AjYG9k6yoqp9Mu1hJ0vI1blj7Q4APD63bB9gOeEFVbQbcDjgX+JepVCdJUmtcYG0DnDK07knA6VX1QYCquhQ4GLj/VKqTJKk1LrBuTjMSEIAkm9BcdzU8H+C5wG0nXpkkSQPGBdZPgJ0Hlh/X/vvFof22AK6YYE2SJN3AuFGC7wben2Rj4GLgn2juUfWlof0eiXcjliRN2ayBVVWHJ9kS2Au4NfAdYK+q+sPMPkk2Bx4P7D/lOiVJy9zY67Cq6gDggDHbL8XzV5KkBTCJ2dolSZo6A0uS1AsGliSpFwwsSVIvGFiSpF4wsCRJvWBgSZJ6wcCSJPWCgSVJ6gUDS5LUCwaWJKkXDCxJUi8YWJKkXjCwJEm9YGBJknrBwJIk9YKBJUnqBQNLktQLBpYkqRcMLElSLxhYkqReMLAkSb1gYEmSesHAkiT1goElSeqFmy92AZKkWWy+OaxcudhVLBm2sCRJvWBgSZJ6wcCSJPWCgSVJ6gUDS5LUCwaWJKkXDCxJUi8YWJKkXjCwJEm9YGBJknohVbXYNWhAkpVVtWqx69Bofj5L17r42ST5QlXtMsu2Xavq2NkeJ9kVYPDxzPJcxxs+zvC/g8+b1PGGH486pi2spceJw5Y2P5+la537bGYLq9auczzedcTjtcJmzPGGjzPquJM+3vDjG3DyW0nqp2PneDzb9i7HG3WcUced9PHGHtsuwSUmyeqq2mGx69Bofj5Ll5/Nus8uwaVnneqDXwf5+SxdfjbrOFtYkqResIUlSeoFA0uS1AsG1hKQ5M5JDkvy3STXJjlhsWtSI8nuST6T5IIkv0lySpKnL3ZdaiR5cpITk1ye5OokZyR5bZJbLHZtmjyHtS8NdwMeA3wL8D/a0vIK4Bzg5cBlNJ/TR5NsVlX/vqiVCWBT4GvA24BfAiuANwC3Bf5x0arSVDjoYglIcrOquq59/Elgs6raeXGrEkAbTJcNrfsocL+q2naRytIYSd4E7AX8efkLbp1il+ASMBNWWnqGw6p1KrDFQteizi7Hnop1kl2C0vztBJy+2EVojSTrAesD9wL+CXivrat1j4ElzUOShwGPB5672LVoLb+lCSyAI4F/XsRaNCV2CUodJdkG+ChwTFUdvrjVaMhOwAOBvWn+oHj34pajabCFJXWQZBPg88B5wDMXuRwNqarvtA+/keQy4IgkB1fVWYtZlybLFpY0hyQbAZ+lOZH/2Kr67SKXpPFmwstRnOsYW1jSGEluDnwC+Evg/lV1ySKXpLndv/33nEWtQhNnYC0B7V/wj2kXbw/cKsmT2+XjquqqxalMwKE0n81LgU2S7Diw7dSq+v3ilCVo7sgLfAX4IXAtTVjtDXzM7sB1jxcOLwHtyfzZ/hrctqrOXbhqNCjJucAdZ9nsZ7PIkrwR2A3YBvgjcDbwn8D7quoPi1iapsDAkiT1goMuJEm9YGBJknrBwJIk9YKBJUnqBQNLktQLBpYkqRcMrClKskd7S/Urk/wiyalJ3r7ANbyhnVttPs/ZI0kl+dNp1TUfSbZrv49bT/i45yY5aJLHHPNaj0/yoyTXtNd2aYQb8/N6I19n5/Zn/O7Tfi1NjoE1JUleDfwH8EXgicCzgWOAv1vgUv4DeNQ8n/M54H7AUplhYzvg9cCtJ3zc3YB3TfiYN9Deq+lI4LvAQ9vX1Wg35udVy4QXDk9JkguAT1fVXkPrc1NvLNf+Alyvqq65KcfpiySPA46lw8wSSTaoqqsXpLCOktwB+Dnw+Kr6zJj9Nqyq3y1cZctXkp2BrwF/XVU/WNxq1JUtrOm5NXDR8MrBsJqtWyLJCUk+ObB8eJLVSZ6Q5IfA1cB92+c+Zui56yW5qJ2yZq0uliS3TPLbJP8wXFd7/A+1j9fqEkyyTbv8lCSHJflVkvOT7J/kZkPH2T3JT5P8LsnXktyzfe4e496sJK9OcmaSq5NcnOQLSW7b/mI5tt3tnPZY5w7VuaJ9z35He+O+JG9J8v0kv2lr/UiS2w695lpdggPv8yOSfK99r76R5G5Dz3tekh+23+NlSb4+vM/AvnvQhBXAMW29b2i3VZJXJHlHkkuB77frN0tyRJLLk1zVfm87jKo9yb5JLmw/k4PTeExb35VJPp3kz+d47++Q5ONJLmm/p7Nmfn6G3pcnJPlx+xl9I8ldh45zs7aeM5P8PslPkjxnxOvtluTk9rUuT3Jckju2227QJZhkk/bn7uL2tU9Mct8b+5movwys6fkO8JIkz0my6QSOtw1wIHAAzWSs5wAnA08d2u/BwG2Ajw0foL0txmeHn5PkTsC9Rz1nyIHAb4AnAx8GXtc+njnODsBRNN/7bsBnOhyTJM8G9gPeTtMd9GLgTOCW7bH2aXd9Ik1X5XCX2n+139dj2n8BtgDeDDwWeBlwJ+CraVqn42wNvA14E/D09jgfT5K21gcB72u//0fT3Hn4RGDjWY73ubZu2u/jfjTdXjP+GdgSeBbNrd0BPk3zPuxD81ndDPhakjsPHftpwApgT5rP5hU07+EbgX8BXkTz83DAHN/zkcBWwMr2e3oTa+7eO+OOA8f++/b7/WKSDQb2+XfgtcAqmvf9aOCDaVrIACR5FvAp4CzgKW3tPwE2H1VYkvVpJrd9BM179QTgUuArM3+A3IjPRH1VVX5N4Qu4B81EnAVcRzOb9L8CtxrYZ+d2+92HnnsC8MmB5cPb/f52aL+XA78C1h9Ydxjww4HlNwCXDSzvRjOr9e0G1r0auAK4Rbu8R/t6f9oub9MuHzn0+qcBRw0sfwL4AW1Xc7vule1z9xjzXr0b+O8x2x/XHmObofUzdb50js9iPZpZ8At40MD6c4GDht7nPwJ/ObDuCe3z7tIu7wOcMs+fhZn373FD64tmxvfBdbu06x88sO6WNL+kDxuq/UyaruGZdSe39W87sO5A4OI56vsNsOuY7TM/fzsNrLtj+1ovapfvTPNz/pyh5x4JfLt9fDPgAuBTY15r+Of1ecA1Q5/JzWkC72034TPZmRH/9/xa2l+2sKakqr4HbE8zyOJQIDR/9a7OjRt9d0FVnTa07uPAn9H8kpu5d9MTaVo5s/k8zS+o3QfWPRU4uuY+J/aloeXTgTsMLN8HOLba3witWc/ZDDgNeEyaLsYVHVpBwz43vCLJo9uuo1/R/GI9v9203RzHOreqfjqwfHr778z3eRpwzySHJHlQklvMs9Zhw7WvAC6tqq/PrKg1LeMHDO17QlVdO7B8Zlv/OUPrNp+jztOAA9J0sW49yz6XVNWJAzX9DDilrRfgYTSBdXSSm898AccDf9t+pn8F3I5mNvWuHt6+zjkDxwT4OjDTTXoak/1MtEQZWFNUVb+vqmOr6h+r6q7A82luBPi8G3G4i0cc/wLgG6zp4nsYsBljAquaAQnHzDwnyV8BfzPuOQN+ObR8DTDYJXRbmpbAoOHlUT5I0yX4FOAk4OIkb5xHcK313iS5D01Qnk/T1XY/YOY+Vhsw3i+HlmdCfAOAqvoKTTfWg2hawpclOTTJLTvWOrZ2mu7BG3zW7bpNOtQ6al1o7pY8m6cCq4FDgJ8lOS3Jw4b2GXXjykvaeqH5uVuPpsX/h4Gvw2laRFsCM13jF46pZdhmNJ/dH4a+9qTpxpzGZ6Ilyhs4LqCq+kCSA4G7tKtmRrMN/zLZBBi+FmW2kYUfA96SZEOaXzynDrUQZnvOse1f00+lCZWvdvgW5nIRNzwXMfLcxKCquo7ml+UhSbYCnkFzHuUCmnMTcx5iaHk3mu/pqTOtvZmT+pNQVUcARyTZnKZFewjwa2DfG3O4oeULac6bDbsNTbftxLV/+OyRZgDNCppuuc8k2bqqLm93G1XTFjRd3bS1/ZHmBorXjdj3EpreAFgTcl1cQROmLx6x7fqbZ074M9ESZQtrSpLc4D94+59pY9b8BT3TTbX9wD5b0XSddPUJYEOaX9K70a2l9CXgFzQtmqfSnC+7dvxTOvk2sOvMAIXWvK47q6qfV9VbaLqyZkahrdXK6WBD4A9DXZPPmE8dXVTVpVV1GPC/rKn1pjoJ2KIdSABcf0fqx9K0pqemqq6rqm8B+wMbsfaNK7dIstNATVsD96I5bwbNHzzrARtX1eoRX9cAZ9D8EXKDkYNjHE9zfuy8Ecf8/ojvYRqfiZYIW1jT8/0kx9CEwyU0//n3obkY9wiAqjo/ybeBNya5iuYPiP2Yx1/SVXVJkhOAg2iG0n+8w3P+kORomlFlWwI3GOZ+I72V5hfuUUn+kyaIX9BuG/VXNwBJDqP5nr9F06X0EJqu01e1u5zR/vvCJEcBV436ZTXgy8DLkryDZkj8TsAzb8w3NKLW/WlawCfQtILvSTMSbyJ/yVfVF5P8P+BjSfYFLqf5udmQZvTiRCXZmObi9iNpRuutT3OL+YuAHw3sehnwoST/AvyOZgDRJTRdflTVGUneR/PZH0jTKtoAuBuwXVU9v6quS/JK4CNJPkIzurNoLqb+r6paPaLEI2lGO56Q5hKEs2m6FlcAF1XVIdP+TLR0GFjT86/A42lmUtiE5hfAiTTdVIMnxf+eZpjzh2laXK+kGf03H0cB7we+Vd1v2X4Uzbm0/6P5a/Qmq6rVSZ5OM5z88azpyvkyTffMbL5JE2wvpPkldybwgqr6dHvcnyXZh2bY90to3qdtxtRxXJJXtfu+oD3+42h+Id9U36b5fJ5G08X1M5outHdO4NgzdgMOBt5B836cDDy0qs6c4GvMuJrm+q+X0pwTuormD4dH1toXMf+M5nN9C80fX6uBp9faF2nvRfMev4Dm5//XNINWPjCzQ1V9NMnVwGuATwK/bV9v5LnOqro6yUPa4+1P0zV6Cc17MjOgZyE+Ey0BznShqUryTOBDwJ2Gglo9keRwmuHfO8y1rzRNtrA0UUneS9Oi+gXNOY7XAp8zrCTdVAaWJm1TmuvONqU5//Ixmm5OSbpJ7BKUJPWCw9olSb1gYEmSesHAkiT1goElSeoFA0uS1Av/H466c51iZN7eAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "file = open(\"data/2_res_new_results_x_separated.pkl\",'rb')\n",
    "all_survivors,all_params = pickle.load(file)\n",
    "\n",
    "heatmap = np.zeros([3, 3])\n",
    "for idx, i in enumerate(all_survivors):\n",
    "    bug1 = sum([j<50 for j in i])\n",
    "    bug2 = sum([j>=50 for j in i])\n",
    "    heatmap[3-bug1, bug2-1]+=1\n",
    "    \n",
    "\n",
    "import pandas as pd\n",
    "xrange = 3\n",
    "yrange = 3\n",
    "Isurvivor = np.repeat(range(1, xrange+1), yrange)\n",
    "IIsurvivor = np.tile(range(1, yrange+1), xrange)\n",
    "df = pd.DataFrame({\"I survivor\": Isurvivor,\n",
    "                  \"II survivor\": IIsurvivor,\n",
    "                  \"percentage\": np.array([heatmap[Isurvivor[i]-1, IIsurvivor[i]-1] for i in range(9)]) })\n",
    "\n",
    "g = sns.jointplot(data=df, x='I survivor', y='II survivor', kind='hist', bins=(xrange, yrange))\n",
    "g.ax_marg_y.cla()\n",
    "g.ax_marg_x.cla()\n",
    "sns.heatmap(data=np.log(df['percentage'].to_numpy().reshape(xrange, yrange)+1), ax=g.ax_joint, cbar=False, cmap='Greys')\n",
    "\n",
    "g.ax_marg_y.barh(np.arange(0.5, yrange), df.groupby(['I survivor'])['percentage'].sum().to_numpy()/len(all_survivors), color='r', alpha = 0.4)\n",
    "g.ax_marg_x.bar(np.arange(0.5, xrange), df.groupby(['II survivor'])['percentage'].sum().to_numpy()/len(all_survivors), color='b', alpha = 0.4)\n",
    "\n",
    "g.ax_joint.set_xticks(np.arange(0.5, xrange))\n",
    "g.ax_joint.set_xticklabels([1, 2, 3], fontsize=15)\n",
    "g.ax_joint.set_yticks(np.arange(0.5, yrange))\n",
    "g.ax_joint.set_yticklabels([3, 2, 1], fontsize=15)\n",
    "g.ax_joint.set_xlabel(\"Surviving strains from species I\", fontsize=15)\n",
    "g.ax_joint.set_ylabel(\"Surviving strains from species II\", fontsize=15)\n",
    "\n",
    "g.ax_marg_x.set_yscale(\"log\")\n",
    "g.ax_marg_y.set_xscale(\"log\")\n",
    "g.ax_marg_x.set_ylim(0.001, 1)\n",
    "g.ax_marg_y.set_xlim(0.001, 1)\n",
    "g.ax_marg_x.tick_params(axis='x', bottom=False, labelbottom=False)\n",
    "g.ax_marg_y.tick_params(axis='y', left=False, labelleft=False)\n",
    "g.ax_marg_x.tick_params(axis='y', left=False, labelleft=False)\n",
    "g.ax_marg_y.tick_params(axis='x', bottom=False, labelbottom=False)\n",
    "\n",
    "g.fig.subplots_adjust(hspace=0.05, wspace=0.02)  # less spaced needed when there are no tick labels\n",
    "plt.savefig(\"heatmap_only_x.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "bde23f58",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.  , 0.  , 0.  ],\n",
       "       [0.  , 0.  , 0.  ],\n",
       "       [0.99, 0.01, 0.  ]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "heatmap/500"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "490e5051",
   "metadata": {},
   "source": [
    "#### Now we modify both of those variables and look at the output:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46ce1f0e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "bd2a9ad0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting x_and_K.py\n"
     ]
    }
   ],
   "source": [
    "%%writefile x_and_K.py\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.integrate import solve_ivp\n",
    "import copy\n",
    "import pickle\n",
    "from math import *\n",
    "from matplotlib.colors import ListedColormap, LinearSegmentedColormap\n",
    "import seaborn as sns\n",
    "import scipy.stats as stats\n",
    "import time\n",
    "import multiprocessing\n",
    "from multiprocessing import Pool, freeze_support\n",
    "\n",
    "# g: maximum growth rate\n",
    "# K: Monod constant\n",
    "def OneStep(t, u, p):\n",
    "    N, R, gmax, K, x = p\n",
    "    Ni, Rj = u[:N], u[N:]\n",
    "    xmat = np.hstack([x.T, 1-x.T])\n",
    "    cmat = np.vstack([Rj for i in range(N)])\n",
    "    gmat = xmat*(gmax*(cmat/(K+cmat)))\n",
    "#     print(gmat)\n",
    "    return np.hstack([np.sum(gmat, axis=1)*Ni, -Ni@gmat])\n",
    "\n",
    "\n",
    "def Dilute(p, p_dilute, u_end):\n",
    "    N, R, gmax, K, x = p\n",
    "    c0, d, t, n, fluc_thres, elim_thres = p_dilute\n",
    "    output = u_end\n",
    "    u_add = np.zeros(N+R)\n",
    "    u_add[-R:] = np.array(c0)\n",
    "    stop, count = 1, 1\n",
    "    utransf = copy.deepcopy(u_end)\n",
    "    while stop>0 and count<n: \n",
    "        u = copy.deepcopy(utransf)\n",
    "        u = u/d + u_add\n",
    "        N_alive = u[:N]>0\n",
    "        N_temp = sum(N_alive)\n",
    "        gmax_temp = gmax[N_alive]\n",
    "        K_temp = K[N_alive]\n",
    "        u_temp = np.hstack([u[:N][N_alive], u[N:]])\n",
    "        x_temp = x[0:,N_alive]\n",
    "        p_temp = N_temp, R, gmax_temp, K_temp, x_temp\n",
    "        sol = solve_ivp(OneStep, t_span=[0, t], y0=u_temp, method=\"LSODA\", args=[p_temp], rtol=1e-6)\n",
    "        utransf_temp = sol.y[:, -1] * (sol.y[:, -1]>elim_thres*d)\n",
    "        u[:N][N_alive] = utransf_temp[:N_temp]\n",
    "        u[-R:] = utransf_temp[-R:]\n",
    "        utransf = u\n",
    "        output = np.vstack([output, utransf])\n",
    "        stop = sum(abs((output[-1, :-R]-output[-2, :-R])) > fluc_thres*output[-1, :-R])\n",
    "        count+=1\n",
    "#     print(\"c\", count)\n",
    "    return output\n",
    "\n",
    "def invade(profile, p_profile, p_dilute, p_invade):\n",
    "    # load variables\n",
    "    N, R, gmax, K, x = p_profile\n",
    "    c0, D, t, n, fluc_thres, elim_thres = p_dilute\n",
    "    p_dilute_once = c0, D, t, 2, fluc_thres, elim_thres\n",
    "    gmax_inv, K_inv, x_inv = p_invade\n",
    "    p_profile_temp = N+1, R, np.vstack([gmax, gmax_inv]), np.vstack([K, K_inv]), np.hstack([x, x_inv])\n",
    "    profile_temp = Dilute(p_profile_temp, p_dilute_once, np.hstack([profile[:-R], [D*elim_thres], profile[-R:]]))\n",
    "    if(profile_temp[-1][-(R+1)]<=0):\n",
    "        return profile, p_profile\n",
    "    else:\n",
    "        profile_temp = Dilute(p_profile_temp, p_dilute, np.hstack([profile[:-R], [5*D*elim_thres], profile[-R:]]))\n",
    "        return profile_temp[-1], p_profile_temp\n",
    "\n",
    "def model(index, N_pool):\n",
    "    \n",
    "    gA1, gA2, KA1, KA2 = 1, 1, 4, 5\n",
    "    gB1, gB2, KB1, KB2 = 0.55, 0.55, 0.05, 0.03\n",
    "\n",
    "    g_diff = np.random.normal(0, 0.05, size = 2*N_pool)\n",
    "    K_exp = np.random.normal(0, 0.3, size = 4*N_pool)\n",
    "    gA1new = [gA1+i for i in g_diff[0:N_pool]]\n",
    "    gA2new = [gA2+i for i in g_diff[0:N_pool]]\n",
    "    gB1new = [gB1+i for i in g_diff[N_pool:2*N_pool]]\n",
    "    gB2new = [gB2+i for i in g_diff[N_pool:2*N_pool]]\n",
    "    KA1new = [KA1*(10**i) for i in K_exp[0:N_pool]]\n",
    "    KA2new = [KA2*(10**i) for i in K_exp[N_pool:2*N_pool]]\n",
    "    KB1new = [KB1*(10**i) for i in K_exp[2*N_pool:3*N_pool]]\n",
    "    KB2new = [KB2*(10**i) for i in K_exp[3*N_pool:4*N_pool]]\n",
    "    \n",
    "    xA, xB = 0.2, 0.7\n",
    "    sigma = 0.1\n",
    "    xAdist = stats.truncnorm(-0.2, 0.8, loc=0, scale=sigma)\n",
    "    xAnew = np.sort(xA+xAdist.rvs(N_pool))[::-1]\n",
    "    xBdist = stats.truncnorm(-0.7, 0.3, loc=0, scale=sigma)\n",
    "    xBnew = np.sort(xB+xBdist.rvs(N_pool))\n",
    "    \n",
    "    plt.scatter(xAnew, 1-xAnew, color = 'b', alpha = 0.2)\n",
    "    plt.scatter(xBnew, 1-xBnew, color = 'r', alpha = 0.2)\n",
    "\n",
    "    N, R = 2*N_pool, 2\n",
    "    gmax = np.array([[gA1new[i], gA2new[i]] for i in range(N_pool)]+[[gB1new[i], gB2new[i]] for i in range(N_pool)])\n",
    "    K = np.array([[KA1new[i], KA2new[i]] for i in range(N_pool)]+[[KB1new[i], KB2new[i]] for i in range(N_pool)])\n",
    "    x = np.array([np.hstack([xAnew, xBnew])])\n",
    "    c0, D, t, n, fluc_thres, elim_thres = [10, 10], 10, 24.0, 5000, 1e-5, 1.0e-6\n",
    "    \n",
    "    start=time.time()\n",
    "    \n",
    "    p_profile = 1, R, gmax[:1, :], K[:1, :], x[:, :1]\n",
    "    p_all = N, R, gmax, K, x\n",
    "    p_dilute = c0, D, t, n, fluc_thres, elim_thres\n",
    "    profile = Dilute(p_profile, p_dilute, np.array([1]+[0 for i in range(R)]))[-1]\n",
    "    \n",
    "    survivors_0 = [0]\n",
    "    survivors = []\n",
    "    while(survivors_0 != survivors):\n",
    "        survivors_0 = survivors\n",
    "        for i in range(N):\n",
    "            _, _, _, _, x_profile = p_profile\n",
    "            if(x[0, i] not in x_profile[0]):\n",
    "                p_invade = gmax[i:i+1], K[i:i+1], x[:, i:i+1]\n",
    "                profile_new, p_profile_new = invade(profile, p_profile, p_dilute, p_invade)\n",
    "                N_temp, R_temp, gmax_temp, K_temp, x_temp = p_profile_new\n",
    "                survivor = profile_new[:-R]>0\n",
    "                p_profile = sum(survivor), R_temp, gmax_temp[survivor, :], K_temp[survivor, :], x_temp[:, survivor]\n",
    "                profile = np.hstack([profile_new[:-R][survivor], profile_new[-R:]])\n",
    "\n",
    "        end=time.time()\n",
    "\n",
    "        N_temp, R_temp, gmax_temp, K_temp, x_temp = p_profile\n",
    "\n",
    "        xlist = list(xAnew)+list(xBnew)\n",
    "        survivors = [xlist.index(i) for i in x_temp[0]]\n",
    "\n",
    "        all_survivors=survivors\n",
    "        all_params=[xAnew, xBnew], gmax_temp, K_temp\n",
    "        outputFileName = 'data/2_res_new_results_xandK_separated_'+str(index)+'.pkl'\n",
    "        pickle.dump((all_survivors, all_params),open(outputFileName, \"wb\"))\n",
    "\n",
    "    return (end-start)\n",
    "\n",
    "runs = 500\n",
    "\n",
    "if __name__ == '__main__':\n",
    "    freeze_support()\n",
    "    cores = multiprocessing.cpu_count()\n",
    "    pool = multiprocessing.Pool(processes=7)\n",
    "    x1 = [i for i in range(runs)]\n",
    "    x2 = [50 for i in range(runs)]\n",
    "    tasks = [(index, N_pool) for index, N_pool in zip(x1,x2)]\n",
    "    pool.starmap(model,tasks)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "b4855076",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAawAAAGyCAYAAAC4FPBHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAlaUlEQVR4nO3deZhkZX238fsrRhYTMWyKEQRNiKhJ1OCIaBTFBZeJYtxXcBmTECNGoohGIUZRBFETUcZoAJfgEhFQXFE0vkRwENyIKJsIYccNEVH4vX+c00xNUV19Gqq6+0zfn+vqa+osderXVT397ec5z3lOqgpJkpa62yx2AZIkdWFgSZJ6wcCSJPWCgSVJ6gUDS5LUC7dd7AImwGGOktZHWewClpreBlaSlcDK448/npUrV3Z6zurV061puVi1arErkLQc2SUoSeqF3rawquoE4ATgxYtdiyRp+nrbwkqyMsnqE044YbFLkSQtgN4GliRpeTGwJEm94DksSVIv9LaF5TksSVpebGFJknqhty0sSdLy0tvAsktQkpYXuwQlSb3Q2xaWJGl5MbAkSb3Q28DyHJYkLS+ew5Ik9UJvW1iSpOXFwJIk9YKBJUnqBQNLktQLBpYkqRd6G1gOa5ek5aW3gSVJWl4MLElSL/T2wmGtX1avXuwK1h+rVi12BdJ09DawnOlCkpaX3nYJOuhCkpYXW1iSpF6whSVJ6gVbWJKkXrCFJUnqBVtYkqRe6G0LS5K0vBhYkqRe6G1geQ5LkpYXz2FJknqhty0sSdLyYmBJknrBwJIk9YKBJUnqhd4OupC0cLxf2eR4v7JbrrctLIe1S9Ly0tvAkiQtLwaWJKkXDCxJUi/0dtCFM11I0vLS2xaWgy4kaXmxhSVJ6gVbWJKkXrCFJUnqBVtYkqResIUlSeqF3rawJEnLi4ElSeqF3gaW57AkaXnxHJYkqRd628KSJC0vBpYkqRcMLElSLxhYkqReSFUtdg23SpLPAlssdh2SNGFXVtXui13EUtL7wJIkLQ92CUqSeqG312ENsIkoaX2U3XffvT772c8udh2LIaNW2sKSpCXqyiuvXOwSlhQDS5LUCwaWJKkXDCxJUi8YWJKkXjCwJEm9YGBJknrBwJIk9YKBJUnqBQNLktQLBpYkqRcMLElSLxhYkqReMLAkSb1gYEmSesHAkiT1goElSeoFA0uS1AsGliSpFwwsSVIvGFiSpF4wsCRJvWBgSZJ6wcCSJPWCgSVJ6gUDS5LUCwaWJKkXDCxJUi8YWJKkXjCwJEm9YGBJknrBwJIk9YKBJUnqBQNLktQLBpYkqRcMLElSLxhYkqReMLAkSb1gYEmSeuG2s21IcvA8jlNV9aoJ1CNJ0kizBhbw1HkcpwADS5I0NbMGVlVtv5CFSJI0juewJEm9MO4c1r3mc6CqOuvWlyNJ0mjjzmF9l+bc1FzS7rfBRCqSJGmEcYH18AWrQpKkOYwbdPGVhSxEkqRxHHQhSeoFA0uS1AvjzmEtmCTbA9sD51fV+YtdjyRp6VmwFlaS30nyziRXJ7lmZuqnJO8CzgG+CJyT5ENJHHEoSVrHQrawXgm8CHgbcDXw0iRbAk8C9gS+CTwEOAR4CXD4AtYmSVriOrWwkuyYZOeB5Y2TvCnJJ5O8tONrPRd4fVW9tqre1i4/Hziwqj5QVd+rqiNoAu0Fc9SzKsmaJGtWr17d8eUlSX3WtYV1OHAK8PV2+RCaVtF/A29JslFVvXWOY9wNOG1g+fT239OG9vsaMDYEq2o1MJNUXS5uliT1XNdzWPcB/geac1HAc4B9qmp3YH/maBG1fgnccWD51+3XtUP7bcASGQwiSVo6ugbW7YGft493bpc/0S5/k6b1NJezgZ1mFqrqxqrauKrOHNrv3sAFHeuSJC0TXQPrPJqgAtgDOKOqrmqXtwB+0eEYbwN+3GG/RwLHdqxLkrRMdO16Owx4d5KnAvcD9hrYtivw7bkOUFX/1eWFquqxHWuSJC0jnQKrqt6X5IfAA4D9quqkgc1XA2+fQm2SJN2k8+CGqvoq8NUR6w+YZEGSJI3SeaaLJFsleUuSk5L8IMm92/UvS/Kg6ZUoSVL3C4dX0Eyf9Fc0I/juAWzYbt4aeMU0ipMkaUbXFtZhwJeAHWimTcrAttOAFROuS5KkdXQ9h3V/4IlVdWOSDG27CthqsmVJkrSuri2snwFbzrLt7sBlkylHkqTRugbWccCBSe4+sK6SbAHsy9pZLyRJmoqugbUfzdRMZ7F2aPt7aKZb+hXwusmXJknSWl0vHP5Je3uR5wK70UxkezXw78DRVfXr6ZUoSdL8Lhy+Hnhf+yVJ0oLqfOGwJEmLadYWVpLLgcdU1RlJrmCOGyVWlUPbJUlTM65L8F2sHa7+LryzryRpEc0aWFV14MDjAxakGkmSZtF1LsFtktx/lm33T7LNZMuSJGldXQddvBt4zizbngUcPplyJEkarWtg7Uwz+e0oX263S5I0NV0DaxPGD7q4/QRqkSRpVl0D6zvAM2fZ9kzge5MpR5Kk0brOdPFm4L+SbAgcCVxCc+PG59Pc1PGvplKdJEmtrnMJHpvk+cBBNOFUNDdxvBh4TlV9cmoVSpLE/OYS/ECSDwL3BDajuXHj2VXlBcWSpKnrHFgAVVVJvk/THXi5YSVJWiidJ79N8rgkpwLXARcCf9quX51ktmu0JEmaiK4zXTwPOB74PrBq6Hk/BF44+dIkSVqrawvrNcBbq+r5wAeHtn0PuNdEq5IkaUjXwLob8IVZtl0H3GEy5UiSNFrXwPoxcL9Ztu0EnDOZciRJGq1rYL0PeH07uGLjdl2S7Aa8EnjvNIqTJGlG12HtbwG2AY4CbmjXnQJsABxRVe+cQm2SJN2k60wXBeyd5DBgN2Bz4GrgS1X1gynWJ0kSMP8Lh8/B81WSpEXQObCS3A7YE1hBM9PFJcCpwFFVdf1UqpMkqdX1wuEdaS4QfhdwH5rzWPdpl89J4nVYkqSp6trCWg38DPiLqrpwZmWSbYFPA+8BHjr58iRJanQd1r4T8LrBsAJol18HPGDShUmSNKhrYF0AbDTLto1oJsOVJGlqugbWfsC/JHng4MokOwP/DLxq0oVJkjSo6zms19LMF3hKksuBy4Gt2q+rgP2T7D+zc1WtmHShkqTlrWtgfbf9kiRpUXSd6WKvaRciSdI4ne84PCzJ7ye5b5INJ1mQJEmjdL1w+MAkbx5YfgTNyMDTgfOS3HtK9UmSBHRvYT0b+P7A8qHA14AHt+sPmnBdkiSto2tg3QU4DyDJNsCfAa+vqq8DbwN2nk55kiQ1uo4S/AWwafv4EcBPquq0dvk6YJNJF9bVO97xjsV6ac1hr70cq7OU3eEOd1jsEqR56RpYXwH2S3IjsC9w3MC2HYAfT7owSZIGde0SfDnwa+AY4KfAawa2PQ/46mTLkiRpXV2vw7qYpitwlMfQdAtKkjQ187rj8ChV9fNJFCJJ0ji3+MJhSZIWkoElSeoFA0uS1AsGliSpF+Y96CLJbYHbDa+vqmsnUpEkSSN0nfx20ySHJ7mEZgj7L0Z8SZI0NV1bWEcCDwPeC5wDXD+tgiRJrSuugNWrF7uKxqpVi11B58DaDXhJVf3nNIuRJGk2XQddXAh4jkqStGi6BtYrgdcm2XaaxUiSNJuucwmemOSRwDlJLqCZAHd4nxWTLU2SpLU6BVaSQ4B9gG/goAtJ0iLoOujiRcBrquqgaRYjSdJsup7DuhY4fZqFSJI0TtfAegewKkmmWYwkSbPp2iW4BfBA4OwkJ3PzQRdVVa+aYF2SJK2ja2A9Bfgt8DvAo0ZsL8DAkiRNTddh7dtPuxBJksbx9iKSpF7oHFhJ7p7k3Um+k+Ti9t/Dk9x9mgVKkgTdLxz+c+DLNLcW+RRwGXAn4K+AZyd5eFV9c2pVSpKWva6DLg4BzgAeO3ijxiSbACe22x8x+fIkSWp07RJcARw8fFfhdvkQmiHvkiRNTdfA+hWw+SzbNqPpKpQkaWq6BtangTcnecjgynb5IOCESRcmSdKgruew/gE4DvhKkitoBl1s1X6dArxiOuVJktToeuHwVcBDkuwOPADYGrgEOLWqPj/F+iRJAjoEVpKNgG8Df19VnwU+O/WqJEkaMuc5rKq6DrgjcOPUq5EkaRZdB118CNhrmoVIkjRO10EXFwJPS7KG5kLhy2hmaJ9RVfXuSRcnSdKMroF1aPvv1sD9R2wvwMCSJE1N11GCzuouSVpUswZRkvOS/Fn7+HVJ7rJwZUmStK5xLac/ADZpH78euOv0y5EkabRxXYLnAn+XZHMgwC5Jtpht56o6cdLFSZI0Y1xgvRo4CngmzaCKt43Zt4ANJliXJEnrmDWwquq4JJsBd6EZ1v5kmntiSZK04MaOEqyqG4GLkuwF/Hc7p6AkSQuu67D2o6ZdiCRJ43h9lSSpFwwsSVIvGFiSpF4wsCRJvdApsJLsmGTngeWNk7wpySeTvHR65UmS1OjawjocWDmwfAjwMmAj4C1J/nFSBSX5nSTbTup4kqT1Q9fAug/wP9AECvAcYJ+q2h3YH3hBl4Mk2TvJuUl+keTUJM8dsdv9gfM71iVJWia6BtbtgZ+3j3dulz/RLn8TuNtcB0jyDOBfga8DBwL/BxyZ5ONJNp5P0ZKk5adrYJ1HE1QAewBnDMx6sQXwiw7H2Bc4pKqeXVWHVNUewKOBhwBfbifZ7STJqiRrkqw55ZRTuj5NktRjXQPrMOBfknwD+HvgnQPbdgW+3eEYfwysM6N7VZ1EE4SbAv+T5B5diqmq1VW1U1XttMsuu3R5iiSp5zoFVlW9D3gkcAzwmKr6wMDmq4G3dzjMz2haY8PHvgDYBbgSOAV4QJeaJEnLS6e5BAGq6qvAV0esP6DjIU4HngR8fMQxfpJkt3bbO2luVyJJ0k06BxZAkh1o7jy80fC2Djdw/CDw8iSbVdXVI57/qyR/CbwbeNR86pIkrf86BVaSewEfAe5Fc/fhYXPewLGqPgZ8bI59bgBWdalJkrS8dG1hHQHcjuYmjmcB10+tIkmSRugaWPcDnlFVn5pmMZIkzabrsPZzGXHeSpKkhdI1sF4B7J/k7tMsRpKk2XTtEjwI+APg+0kuAH46vENVrZhcWZIkratrYH23/ZIkaVF0Cqyq2mvahUiSNI53HJYk9cKsLawkBwPvrKqL2sdjVdUrJ1qZJEkDxnUJPhX4EHBR+3icAgwsSdLUzBpYVbX9qMeSJC0Gz2FJknqh82ztSQI8GNiB0bO1Hz7BuiRJWkfX2drvBJxEM1t7sXbG9sH7VhlYkqSp6doleCjNHYO3oQmrBwLbAf8E/JCm1SVJ0tR07RJ8GPAy4JJ2OVV1IfCmJLehaV09Zgr1SZIEdG9h3RG4oqpuBH4ObDWw7RRglwnXJUnSOroG1vnA1u3j7wHPHti2ErjZLe8lSZqkrl2CJwKPBj4K/AtwXJKLgN8A2wKvmk55kiQ1uk5+u9/A488k2QXYA9gY+EJVfWZK9UmSBHQIrCQbAvsCn6qqbwFU1RpgzZRrkyTpJnOew6qqXwOvoRl4IUnSoug66OJU4M+nWYgkSeN0HXTxSuDDSa6nGYBxGevOckFVXTvh2iRJuknXwDq1/fedwDtm2WeDW1+OJEmjdQ2sFzDUopIkaSF1HdZ+5JTrkCRprE6DLpKcl+TPZtl2nyTnTbYsSZLW1XWU4HbAhrNs2wS460SqkSRpFrN2CSa5A+tee3XnJNsO7bYR8Azg4smXJknSWuPOYb0ceD3NYIsCjp1lvwCvmHBdkiStY1xgfZhm+qUAx9NMz3T20D7XA2e398aSJGlqZg2sqvohzd2ESfJw4PSqumahCpMkaVDXQRf/C2w5s5DGqiRvT7JyOqVJkrRW18A6kuac1owDgcOB3YFjk+w52bIkSVpX18C6P/AlgCS3Af4G2L+q7gm8EdhnKtVJktTqOjXTpsBV7eM/BzYDPtQufwlHCUrS5G25JaxatdhVLBldW1gXAfdqHz8e+H5VzVx7tSlw3aQLkyRpUNcW1vuBg5M8kiawXj2wbWeaQRmSJE1N18lvD0pyMfAA4KU0ATZjM+Dfp1CbJEk36drCoqqOBo4esf6vJ1qRJEkjdD2HJUnSojKwJEm9YGBJknrBwJIk9YKBJUnqhU6jBJM8b8zmG4GfA9+qqh9NpKp52GeffRb6JdXRNdc4uf9S9vjHP36xS9AY973vfRe7hCWn67D2I2lu4gjN/bFmDK6rJJ8Cnu1tSCRJkzafyW/PAfYDdgS2aP99NXAu8Ejg+cBfAG+ZfJmSpOWuawvrUODwqnrHwLqraaZruh54fVU9LMmdaCbC3XvCdUqSlrmuLawHMft8gf9LM2UTwOnA5re2KEmShs1ntvY9Z9m2V7sd4PdZexsSSZImpmuX4GuA/0xyH+AE4ApgS2AlzW1HntHu9yjgvyddpCRJXWdr/1iS84FXAc8C7gxcCnwD2KuqTm/3+9tpFSpJWt7mM1v7GuCpU6xFkqRZOdOFJKkXOrewkjwFeDJwV2Cj4e1VtWKCdUmStI6uUzMdALwO+BZwFnD9FGuSJAFccQWsXr3YVSysVatm3dS1hfVC4M1Vtf9ECpIkaZ66nsP6PeCkaRYiSdI4XQPrGGD3aRYiSdI4XbsETwLekmQL4AvAT4d3qKoTJ1iXJEnr6BpYH2n/3Y5mVvZhBWwwiYIkSRqla2BtP9UqJEmaQ9epmRb8TsKSJA2aNbCSbFJV1848nutAM/tKkjQN41pYv0jyoKo6DbiG5jzVOJ7DkiRNzbjAegFw7sDjuQJLkqSpmTWwquqogcdHLkg1kiTNotOFw0mOSvLYJJ0ny5UkaZK6znRxb+DTwKVJ3pvkEUkyxbokSVpHp8Cqqp2AewCHAjsBXwQuSfJvSf5iivVJkgTM4waOVXV+VR1UVfcDdgQOB3YFTk5y4ZTqkyQJuIV3HK6qs4H3A/8BXAr8wSSLkiRp2LwGUSS5M/BU4OnAzsBPgE/SzOYuSdLUdL3j8F/ThNRDgF8CxwFvAj5fVb+dXnmSJDW6trAOAT5F07r6TFX9enolSZJ0c10DayvnCpQkLaauw9oNK0nSoho3W/vlwGOq6owkVzDHXIJVtdWki5Mkaca4LsF3AZcNPHbyW0nSohk3+e2BA48PWJBqJEmaRdfJbw9MsuO0i5EkaTZdZ7p4CfDdJN9Jsn+Se0yzKEmShnUNrLsAjwJOAfYBfpBkTZJXJNl2WsVJkjSj67D2G6vqS1X1EmBr4HHAt4HXAOcn+doUa5Qkaf6T31bVDVX1OeBvgL1pJr990KQLkyRp0Hwnv/0dYHeaeQVXAhsDXwFeN/nSJElaq+vktzMh9STgDsDXgFcDH6uqK6ZWnSRJra4trBOB04ADgY9W1f9NryRJkm5uzsBquwEfCpxfVRdPvyRJkm6uy6CLG4CTgHtOuRZJkmY1Z2BV1Y3AD4E7Tb8cSZJG6zqs/TXA65L8yTSLkSRpNl0HXbwW2Bw4M8nFNLO4rzN7e1WtmHBtkiTdpGtgfbf9kiRpUXQKrKraa9qFDEqyEbBVVV24kK8rSVq65j010wJ5PHD+YhchSVo6us508dG59qmqp936ciRJGq3rOawtR6zbDPhj4Crg7C4HSfKlW/F6g8dZBazqeCxJ0nqg6zmsh49an2Qb4FjgsI6v91CacDtrjv02mqOe1cDqtoYat68kaf0wr9nah1XVj5McBBwMnNDhKd8Fzq6qp4/bKclTgI/cmtokSeuXSQy6uAG4a8d9TwV27rBfAbnFFUmS1jtdB13ca8Tq2wE7Am8AvtHx9Q4GPt1hvxOB7TseU5K0DMznwuFR54pCE1Yv6nKQqjoXOLfDfr8CftSxNknSMtA1sEYNurgOuMhbjkiSFkLXUYJfmXYhkiSNM+9Rgkk2AV5Ic3+sS4Gjq8ruO0nSVM0aWEkOBVZW1Q4D636P5pzVHwE/ATYFXpFkRVX9YNrFSpKWr3HD2h8OfHBo3b7ADsCLq2oL4C7ABcA/TaU6SZJa4wJrO+D0oXV/BZxVVe8HqKorgEOBB0+lOkmSWuMC67Y0IwEBSLIZzXVXw/MBXgDceeKVSZI0YFxg/QDYdWD5Ce2/nxvabyvg6gnWJEnSzYwbJfhvwHuTbApcBvw9zT2qPj+036PxbsSSpCmbNbCq6sgkWwN7A3cEvgnsXVW/mdknyZbAE4EDp1ynJGmZG3sdVlUdBBw0ZvsVeP5KkrQAJjFbuyRJU2dgSZJ6wcCSJPWCgSVJ6gUDS5LUCwaWJKkXDCxJUi8YWJKkXjCwJEm9YGBJknrBwJIk9YKBJUnqBQNLktQLBpYkqRcMLElSLxhYkqReMLAkSb1gYEmSesHAkiT1goElSeoFA0uS1AsGliSpFwwsSVIvGFiSpF4wsCRJvXDbxS5AkjSLLbeEVasWu4olwxaWJKkXDCxJUi8YWJKkXjCwJEm9YGBJknrBwJIk9YKBJUnqBQNLktQLBpYkqRcMLElSL6SqFrsGDUiyqqpWL3YdGs3PZ+laHz+bJJ+tqt1n2bayqk6Y7XGSlQCDj2eW5zre8HGG/x183qSON/x41DFtYS09Thy2tPn5LF3r3WczW1i1Vs7xeOWIx+uEzZjjDR9n1HEnfbzhxzfj5LeS1E8nzPF4tu1djjfqOKOOO+njjT22XYJLTJI1VbXTYteh0fx8li4/m/WfXYJLz3rVB78e8vNZuvxs1nO2sCRJvWALS5LUCwaWJKkXDKwlIMkfJjkiybeS3JDk5MWuSY0kT01yfJKLk1yT5PQkz1zsutRI8pQkpyS5Ksl1Sc5O8tokt1vs2jR5DmtfGu4NPA74OuB/tKXlH4DzgZcDV9J8Th9OskVV/euiViaAzYEvA28FfgqsAA4A7gz83aJVpalw0MUSkOQ2VXVj+/jjwBZVteviViWANpiuHFr3YeBBVbX9IpWlMZK8Edgb+P3yF9x6xS7BJWAmrLT0DIdV6wxgq4WuRZ1dhT0V6yW7BKX52wU4a7GL0FpJNgA2BO4P/D3wbltX6x8DS5qHJLsBTwResNi1aB2/pAksgKOBf1zEWjQldglKHSXZDvgwcFxVHbm41WjILsBfAK+g+YPi3xa3HE2DLSypgySbAZ8BLgSes8jlaEhVfbN9+LUkVwJHJTm0qs5dzLo0WbawpDkk2QT4FM2J/MdX1S8XuSSNNxNejuJcz9jCksZIclvgY8AfAQ+uqssXuSTN7cHtv+cvahWaOANrCWj/gn9cu/gHwB2SPKVdPrGqrl2cygQcTvPZvAzYLMnOA9vOqKpfL05ZguaOvMAXge8BN9CE1SuAj9gduP7xwuEloD2ZP9tfg9tX1QULV40GJbkAuNssm/1sFlmSNwB7ANsBvwXOA/4DeE9V/WYRS9MUGFiSpF5w0IUkqRcMLElSLxhYkqReMLAkSb1gYEmSesHAkiT1goE1RUn2bG+p/oskP0lyRpK3LXANB7Rzq83nOXsmqSS/O6265iPJDu33cccJH/eCJIdM8phjXuuJSf43yfXttV0a4Zb8vN7C19m1/Rm/z7RfS5NjYE1JklcD/w58Dngy8DzgOOAvF7iUfwceM8/nfBp4ELBUZtjYAXg9cMcJH3cP4J0TPubNtPdqOhr4FvCI9nU12i35edUy4YXDU5LkYuCTVbX30Prc2hvLtb8AN6iq62/NcfoiyROAE+gws0SSjarqugUprKMkdwV+DDyxqo4fs9/GVfWrhats+UqyK/Bl4E+q6ruLW426soU1PXcELh1eORhWs3VLJDk5yccHlo9MsibJk5J8D7gOeGD73McNPXeDJJe2U9as08WS5PZJfpnkb4frao//gfbxOl2CSbZrl5+W5IgkP0tyUZIDk9xm6DhPTfLDJL9K8uUk92ufu+e4NyvJq5Ock+S6JJcl+WySO7e/WE5odzu/PdYFQ3WuaN+zX9HeuC/Jm5N8J8k1ba0fSnLnoddcp0tw4H1+VJJvt+/V15Lce+h5L0zyvfZ7vDLJV4b3Gdh3T5qwAjiurfeAdlsl+Yckb09yBfCddv0WSY5KclWSa9vvbadRtSfZL8kl7WdyaBqPa+v7RZJPJvn9Od77uyb5aJLL2+/p3Jmfn6H35UlJvt9+Rl9Lcq+h49ymreecJL9O8oMkzx/xenskOa19rauSnJjkbu22m3UJJtms/bm7rH3tU5I88JZ+JuovA2t6vgm8NMnzk2w+geNtBxwMHEQzGev5wGnA04f2exhwJ+Ajwwdob4vxqeHnJLk78OejnjPkYOAa4CnAB4HXtY9njrMTcAzN974HcHyHY5LkecD+wNtouoP+BjgHuH17rH3bXZ9M01U53KX2n+339bj2X4CtgDcBjwf2Ae4OfClN63ScbYG3Am8Entke56NJ0tb6UOA97ff/WJo7D58CbDrL8T7d1k37fTyIpttrxj8CWwPPpbm1O8Anad6HfWk+q9sAX07yh0PHfgawAtiL5rP5B5r38A3APwF/TfPzcNAc3/PRwDbAqvZ7eiNr7947424Dx35W+/1+LslGA/v8K/BaYDXN+34s8P40LWQAkjwX+ARwLvC0tvYfAFuOKizJhjST2z6K5r16EnAF8MWZP0BuwWeivqoqv6bwBfwpzUScBdxIM5v0PwN3GNhn13b7fYaeezLw8YHlI9v97ju038uBnwEbDqw7AvjewPIBwJUDy3vQzGp9l4F1rwauBm7XLu/Zvt7vtsvbtctHD73+mcAxA8sfA75L29Xcrntl+9w9x7xX/wb815jtT2iPsd3Q+pk6XzbHZ7EBzSz4BTx0YP0FwCFD7/NvgT8aWPek9nn3bJf3BU6f58/CzPv3hKH1RTPj++C63dv1DxtYd3uaX9JHDNV+Dk3X8My609r6tx9YdzBw2Rz1XQOsHLN95udvl4F1d2tf66/b5T+k+Tl//tBzjwa+0T6+DXAx8IkxrzX88/pC4Pqhz+S2NIH31lvxmezKiP97fi3tL1tYU1JV3wZ2pBlkcTgQmr961+SWjb67uKrOHFr3UeD3aH7Jzdy76ck0rZzZfIbmF9RTB9Y9HTi25j4n9vmh5bOAuw4sPwA4odrfCK1Zz9kMOBN4XJouxhUdWkHDPj28Islj266jn9H8Yr2o3bTDHMe6oKp+OLB8VvvvzPd5JnC/JIcleWiS282z1mHDta8Arqiqr8ysqLUt44cM7XtyVd0wsHxOW//5Q+u2nKPOM4GD0nSxbjvLPpdX1SkDNf0IOL2tF2A3msA6NsltZ76Ak4D7tp/pHwN3oZlNvatHtq9z/sAxAb4CzHSTnslkPxMtUQbWFFXVr6vqhKr6u6q6F/AimhsBvvAWHO6yEce/GPgaa7v4dgO2YExgVTMg4biZ5yT5Y+DPxj1nwE+Hlq8HBruE7kzTEhg0vDzK+2m6BJ8GnApcluQN8wiudd6bJA+gCcqLaLraHgTM3MdqI8b76dDyTIhvBFBVX6TpxnooTUv4yiSHJ7l9x1rH1k7TPXizz7pdt1mHWketC83dkmfzdGANcBjwoyRnJtltaJ9RN668vK0Xmp+7DWha/L8Z+DqSpkW0NTDTNX7JmFqGbUHz2f1m6Gsvmm7MaXwmWqK8geMCqqr3JTkYuGe7amY02/Avk82A4WtRZhtZ+BHgzUk2pvnFc8ZQC2G255zQ/jX9dJpQ+VKHb2Eul3LzcxEjz00MqqobaX5ZHpZkG+DZNOdRLqY5NzHnIYaW96D5np4+09qbOak/CVV1FHBUki1pWrSHAT8H9rslhxtavoTmvNmwO9F0205c+4fPnmkG0Kyg6ZY7Psm2VXVVu9uomrai6eqmre23NDdQvHHEvpfT9AbA2pDr4mqaMP2bEdtuunnmhD8TLVG2sKYkyc3+g7f/mTZl7V/QM91UOw7ssw1N10lXHwM2pvklvQfdWkqfB35C06J5Os35shvGP6WTbwArZwYotOZ13VlV/biq3kzTlTUzCm2dVk4HGwO/GeqafPZ86uiiqq6oqiOA/2ZtrbfWqcBW7UAC4KY7Uj+epjU9NVV1Y1V9HTgQ2IR1b1y5VZJdBmraFrg/zXkzaP7g2QDYtKrWjPi6Hjib5o+Qm40cHOMkmvNjF4445ndGfA/T+Ey0RNjCmp7vJDmOJhwup/nPvy/NxbhHAVTVRUm+AbwhybU0f0Dszzz+kq6qy5OcDBxCM5T+ox2e85skx9KMKtsauNkw91voLTS/cI9J8h80Qfzidtuov7oBSHIEzff8dZoupYfTdJ2+qt3l7PbflyQ5Brh21C+rAV8A9knydpoh8bsAz7kl39CIWg+kaQGfTNMKvh/NSLyJ/CVfVZ9L8v+AjyTZD7iK5udmY5rRixOVZFOai9uPphmttyHNLeYvBf53YNcrgQ8k+SfgVzQDiC6n6fKjqs5O8h6az/5gmlbRRsC9gR2q6kVVdWOSVwIfSvIhmtGdRXMx9X9W1ZoRJR5NM9rx5DSXIJxH07W4Ari0qg6b9meipcPAmp5/Bp5IM5PCZjS/AE6h6aYaPCn+LJphzh+kaXG9kmb033wcA7wX+Hp1v2X7MTTn0v6P5q/RW62q1iR5Js1w8ieytivnCzTdM7P5H5pgewnNL7lzgBdX1Sfb4/4oyb40w75fSvM+bTemjhOTvKrd98Xt8Z9A8wv51voGzefzDJourh/RdKG9YwLHnrEHcCjwdpr34zTgEVV1zgRfY8Z1NNd/vYzmnNC1NH84PLrWvYj5RzSf65tp/vhaAzyz1r1Ie2+a9/jFND//P6cZtPK+mR2q6sNJrgNeA3wc+GX7eiPPdVbVdUke3h7vQJqu0ctp3pOZAT0L8ZloCXCmC01VkucAHwDuPhTU6okkR9IM/95prn2labKFpYlK8m6aFtVPaM5xvBb4tGEl6dYysDRpm9Ncd7Y5zfmXj9B0c0rSrWKXoCSpFxzWLknqBQNLktQLBpYkqRcMLElSLxhYkqRe+P/lF4VovsxR7AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "runs = 500\n",
    "all_survivors = []\n",
    "all_params = []\n",
    "for index in range(runs):\n",
    "    openfile = open('data/2_res_new_results_xandK_separated_'+str(index)+'.pkl', 'rb')\n",
    "    survivor_pool, params_pool = pickle.load(openfile)\n",
    "    all_survivors.append(survivor_pool)\n",
    "    all_params.append(params_pool)\n",
    "    \n",
    "heatmap = np.zeros([3, 3])\n",
    "for idx, i in enumerate(all_survivors):\n",
    "    bug1 = sum([j<50 for j in i])\n",
    "    bug2 = sum([j>=50 for j in i])\n",
    "    heatmap[3-bug1, bug2-1]+=1\n",
    "    \n",
    "\n",
    "import pandas as pd\n",
    "xrange = 3\n",
    "yrange = 3\n",
    "Isurvivor = np.repeat(range(1, xrange+1), yrange)\n",
    "IIsurvivor = np.tile(range(1, yrange+1), xrange)\n",
    "df = pd.DataFrame({\"I survivor\": Isurvivor,\n",
    "                  \"II survivor\": IIsurvivor,\n",
    "                  \"percentage\": np.array([heatmap[Isurvivor[i]-1, IIsurvivor[i]-1] for i in range(9)]) })\n",
    "\n",
    "g = sns.jointplot(data=df, x='I survivor', y='II survivor', kind='hist', bins=(xrange, yrange))\n",
    "g.ax_marg_y.cla()\n",
    "g.ax_marg_x.cla()\n",
    "sns.heatmap(data=np.log(df['percentage'].to_numpy().reshape(xrange, yrange)+1), ax=g.ax_joint, cbar=False, cmap='Greys')\n",
    "# g.ax_joint.hlines(y=np.arange(0, 3), xmin=np.full(3, 0), xmax=np.full(3, 3), color=\"black\")\n",
    "# g.ax_joint.vlines(x=np.arange(0, 3), ymin=np.full(3, 0), ymax=np.full(3, 3), color=\"black\")\n",
    "\n",
    "g.ax_marg_y.barh(np.arange(0.5, yrange), df.groupby(['I survivor'])['percentage'].sum().to_numpy()/len(all_survivors), color='r', alpha = 0.4)\n",
    "g.ax_marg_x.bar(np.arange(0.5, xrange), df.groupby(['II survivor'])['percentage'].sum().to_numpy()/len(all_survivors), color='b', alpha = 0.4)\n",
    "\n",
    "g.ax_joint.set_xticks(np.arange(0.5, xrange))\n",
    "g.ax_joint.set_xticklabels([1, 2, 3], fontsize=15)\n",
    "g.ax_joint.set_yticks(np.arange(0.5, yrange))\n",
    "g.ax_joint.set_yticklabels([3, 2, 1], fontsize=15)\n",
    "g.ax_joint.set_xlabel(\"Surviving strains from species I\", fontsize=15)\n",
    "g.ax_joint.set_ylabel(\"Surviving strains from species II\", fontsize=15)\n",
    "\n",
    "g.ax_marg_x.set_yscale(\"log\")\n",
    "g.ax_marg_y.set_xscale(\"log\")\n",
    "g.ax_marg_x.set_ylim(0.001, 1)\n",
    "g.ax_marg_y.set_xlim(0.001, 1)\n",
    "g.ax_marg_x.tick_params(axis='x', bottom=False, labelbottom=False)\n",
    "g.ax_marg_y.tick_params(axis='y', left=False, labelleft=False)\n",
    "g.ax_marg_x.tick_params(axis='y', left=False, labelleft=False)\n",
    "g.ax_marg_y.tick_params(axis='x', bottom=False, labelbottom=False)\n",
    "\n",
    "g.fig.subplots_adjust(hspace=0.05, wspace=0.02)  # less spaced needed when there are no tick labels\n",
    "plt.savefig(\"heatmap_both.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "7fa8adb9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.   , 0.   , 0.   ],\n",
       "       [0.046, 0.002, 0.   ],\n",
       "       [0.882, 0.06 , 0.01 ]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "heatmap/500"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7318771",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7fd2f68c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52ba6433",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79d4287c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f573d0a8",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86a6e35f",
   "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}