thanos.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# **Thanos snaps**\n",
    "\n",
    "## [Riddler Express, Apr 26, 2019](https://fivethirtyeight.com/features/how-many-earthlings-would-survive-63-thanos-snaps/)\n",
    "\n",
    "### solution by [Laurent Lessard](https://laurentlessard.com)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 351,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import binom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Obtain analytic distribution of the number of snaps that will occur"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 324,
   "metadata": {},
   "outputs": [],
   "source": [
    "# P(s,N) is the probability there will be s more snaps given there are N thanoses currently\n",
    "Nmax = 250\n",
    "p = 0.5  # probability of survival in one snap\n",
    "P = np.zeros((Nmax+1,Nmax+1))\n",
    "P[0,0] = 1.0\n",
    "    \n",
    "for s in range(1,Nmax+1):\n",
    "    for N in range(1,Nmax+1):\n",
    "        P[s,N] = np.dot( P[s-1,range(N)],  binom(N-1,p).pmf(range(N)) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 321,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEuCAYAAADx63eqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XucVXW9//HXW9BEQbxzElC88POEoikTWpoNaidMD9ZJS1OT0jj8DK2kX1Ka5aVfdrHrsVNkipY6GeUJlVJTscy84Y3QVERExHuKoniZ/Jw/1nd0sdkzs2dmr9kzrPfz8diPvdftuz7fvdfan7W+66aIwMzMrCzWaXQAZmZmvcmJz8zMSsWJz8zMSsWJz8zMSsWJz8zMSsWJz8zMSsWJrwEk/UTSV+pU1taSVkoakLrnSTq2HmWn8n4v6eh6ldeF+Z4p6RlJT/T2vCviaJa0rIHz/7CkR9NvvFsXp50l6cyiYjOQNEpSSBrY6Fisdk58dSZpiaRVkl6U9LykmyRNlfTmdx0RUyPijBrL2r+jcSJiaUQMjoh/1iH2r0n6ZUX5B0TEBT0tu4txjASmA2Mi4l96c9590HeAaek3vrOtZ26Dp+0Vkl7Kdb+3gTFbAdKGTEgan+u3g6RuXYwtaX9Jd6Tl5lFJH039N5f0F0nPpv+wv0raq1716Auc+Irx7xExBNgGOAs4Cfh5vWeyFm9lbgM8GxFPNTqQeurm77UNsLCyZ26DZ3BEDE69d831+3OPgrW+6h9Aj/fiJY0BLgZOBoYC7wTmp8ErgU8BWwCbAN8ELl+b/m+c+AoUESsiYg7wMeBoSTvD6k1QaevqirRl9Q9Jf5a0jqRfAFuTLXArJX0x16xyjKSlwHXtNLVsL+lWSSsk/U7SpmleazTbte1VSpoIfBn4WJrf3Wn4m02nKa5TJD0i6SlJF0oamoa1xXG0pKWpmfLk9r4bSUPT9E+n8k5J5e8PXANsleKYVWXaZknLJE1PcTwu6ZO54as190qaLOnGXHdIOk7Sg2nP/AxJ26ct2xckXSppvYp5fjnVaYmkI3L93ybpO6nOTyprxh5UEedJyppsz69Sl6rfaSp3JTAAuFvSQ+19l53YRNKVqZ63SNo+N+8fpC39FyTNz+8lKtv7vzTF86KkhZKacsPfkb7n59OwSblhH5R0b5ruMUlfqBaYsr2VG9Jy+oykX+WGhbKWkgclPSfpHElKw7aXdJ2yPZJnJF0kaePctEskfSnF8Jyk8yWtn4ZVXd+qxHaapB+lz+sq2yv6VuoeJOkVSZvkJjmi2nKfft8Zkh5K8V6qt9bHLq0zyQXALpLe18l4nTkF+GlE/D4iWiPi2Yh4CCAiXomI+yPiDUDAP8kS4KY9nGffERF+1fEFLAH2r9J/KfB/0+dZwJnp8zeAnwDrptd7AVUrCxgFBHAhsCEwKNdvYBpnHvAYsHMa5zfAL9OwZmBZe/ECX2sbNzd8HnBs+vwpYBGwHTAY+C3wi4rYfpbi2hV4FXhHO9/ThcDvgCFp2geAY9qLs2LaZqAVOD19Zx8EXgY2qYw5dU8Gbsx1BzAH2AjYKcV5barXUOBe4OiKeX0XeBvwPuAlYMc0/PuprE1TXS4HvlEx7TfTtIOq1KXd7zQX6w41LHdrjEe2nP0DGA8MBC4CWnLDjwQ2S8OmA08A6+eWhVfSdzuAbDm9OQ1bN8X8ZWA9YF/gxdx38jjw3vR5E2D3dmK+hGyPYx1gfWDvivpcAWxMtgH4NDAxDdsBeH/6TrcA/gR8v2KZ/hswMv0uf6GG9a0itn2BBenze4CHgFtyw+6uZbkHPgfcDIxI8f4UuKSb68wssr29E0jLc/ouIjfOj4Hn23ndkxtvMXAGsCD9Xr8ENq2Y3z3Aa20xNuo/tYhXwwNY2160n/huBk5On2flVsTTyRLAGn9ulWXlVpTtqvTLJ76zcsPHpIV3AD1PfNcCx+WG7Qi8TvbH2RbHiNzwW4HDqtRrQFrBx+T6/ScwL31eI86K6ZuBVW11Tv2eAvasjDl1T2bNxLdXrns+cFKu+2zSHylvJa8Nc8MvBb5CtjX8ErB9bti7gYdz075GSibt1KXd7zQXa08S37m57g8Cf++gjOfImkvbloU/VixHq9Ln95IlyXVywy8BvpY+L02/50adxHwhMDO/zFTUJ58ILwVmtFPOh4A7K5bpqRX1fqiz9a2izEFkiX8zYAZZkl9GtnFyGvDDivWv6nIP3Afslxv2drqxzuR+zzPJEuhS4AAqEl+tr7RcLgH+T6rTb4CLqoy3PnA4aUNwbXm5qbP3DCfb+q70bbKt56slLZY0o4ayHu3C8EfItmw3rynKjm2VysuXPRAYluuXPwvzZbKVqtLmZHsKlWUN70Isz0ZEaw3zas+Tuc+rqnTny3ouIl7KdT9C9l1sAWwAzE9NZ88Df0j92zwdEa90EEct32lPtPt7pKbi+1JT4/Nke7ubdzDt+sqa1LcCHo2sKSwfd9vv9xGyZPNIasp8dzuxfZFs4+HW1Fz6qVpil7SlpJbUjPoC2d5K5fJduQ5slT7XtL5FxCrgdrI9/H2AG4CbgL1SvxtqiZXsGO1lueXjPrKmw66uM/nYXiXbWzuD7PvrjlXA+RHxQESsBP4/2W9WOa9XIuISYIakXbs5rz7Hia8XSHoX2Z/CjZXDIuLFiJgeEdsB/w6cKGm/tsHtFNle/zYjc5+3JtvCfIZs72SDXFwDWP1PurNyl5OtyPmyW1k9adTimRRTZVmPdbGc9qxWT6CnZ4ZuImnDXPfWZN/FM2R/IDtFxMbpNTTeOtkEeu877ZJ0PO8k4KNkTcQbAyuo7Y90OTCy4tjYm79fRNwWEQcDWwL/Q7a3toaIeCIiPh0RW5HtIf5Y0g41zP8bZN/rLhGxEVmTbWXclevA8jTPjta3SjeQNWvuBtyWuj9A1nT8pxrihCwBH5BbPjaOiPUjoqfL+vlkGyofzvdUdox5ZTuv/ElS99D5spm3Lllz/FrBia9AkjaSdBDQQtaEuKDKOAelg/wCXiDbGmy7NOFJurewHSlpjKQNyJp2Zkd2ucMDZFvtB0pal+wA99ty0z0JjKp2sD+5BPi8pG0lDSbbSvxVxZ5Xp1IslwJflzRE0jbAiWRb7vVwF/AfkjZIf6TH1KHM0yStlxLGQcCv0x7Pz4DvSdoSQNJwSR/oQrl1+U67YQhZgn0aGCjpVLJjnrW4hWzj4ovpxI9msiTSkr6jIyQNjYjXeWuZXoOkQyWNSJ3Pkf0R13JZzhCyMw+flzQc+H9VxvmMpBHpRJIvA79K8+xofat0A/AJ4N6IeI3UhE7WlP10DXFCdjzx62kZR9IWkg6ucdp2peXja2QbL/n+UyN3tm/Fa6fcqOcDn5S0XfqfOInsmCqS9pS0d/otB0k6iWwP9Zaext1XOPEV43JJL5Jt7Z1MdmLEJ9sZdzTwR7IV+a/AjyNiXhr2DeCU1ExS9cy4dvyC7HjAE2Rt9CdAdpYpcBxwLtnW+Utkxy3a/Dq9PyvpjirlnpfK/hPwMNkxkOO7EFfe8Wn+i8n2hC9O5dfD98iOYTxJdhbcRT0s7wmyP+blqaypEfH3NOwksqazm1Oz2x/JjtPVqp7faVdcBfyebGPokTTfzprQAUhJYBLZMaZnyE6o+ETuOzkKWJK+j6lke2TVvAu4RdnZq3OAz0bEwzWEcBqwO9ke6pVkJwRVuhi4mmz5WsxblwB0tL5VuonsWF/b3t29ZN9TrXt7AD8gq9vV6T/hZmCPLkzfkUvITkzpsog4j+wY6y1kv/+rpP8Jso3hc4Bnyf4nPggcGBHLexpwX9F29qCZ2VpB0hKyk5v+2OhYrG/yHp+ZmZVKoYlP0kRJ90taVO3sKWUXFj8t6a70qts9Js3MzKoprKkznTH4ANmFpsvIzoo6PCLuzY0zGWiKiGmFBGFmZlahyD2+8cCiiFicDoa3AD0+m8nMzKwnikx8w1n9LLFlVL9A+SOS7pE0W9ld+c3MzApT5N22q10IW9muejnZfetelTSV7NTzfdcoSJoCTAEYNGjQuJEji82Pb7zxBuus07/P++nvdejv8YPr0Fe4Dn1Db9ThgQceeCYituh0xKLuhUZ2z8Krct1fAr7UwfgDgBWdlTtu3Lgo2vXXX1/4PIrW3+vQ3+OPcB36Ctehb+iNOgC3R4Pv1XkbMDrdkWI94DCyCznfJOntuc5JZPexMzMzK0xhTZ0R0SppGtkdIgYA50XEQkmnk2XlOcAJyp7j1Up2A+fJRcVjZmYGxR7jIyLmAnMr+p2a+/wlsiZQMzOzXtG/j5aamZl1kROfmZmVihOfmZmVihOfmZmVihOfmZmVihOfmZmVihOfmZmVihOfmZmVihOfmZmVihOfmZmVihOfmZmVihOfmZmVihOfmZmVihOfmZmVihOfmZmVihOfmZmVihOfmZmVihOfmZmVihOfmZmVihOfmZmVihOfmZmVihOfmZmVihOfmZmVihOfmZmVihOfmZmVihOfmZmVihOfmZmVihOfmZmVihOfmZmVysBGB2DWl4yacWXdypo+tpXJdSxvyVkH1q0sszLzHp+ZmZWKE5+ZmZWKE5+ZmZWKE5+ZmZWKE5+ZmZWKE5+ZmZWKE5+ZmZWKE5+ZmZWKE5+ZmZVKoYlP0kRJ90taJGlGB+MdIikkNRUZj5mZWWGJT9IA4BzgAGAMcLikMVXGGwKcANxSVCxmZmZtitzjGw8siojFEfEa0AIcXGW8M4BvAa8UGIuZmRlQbOIbDjya616W+r1J0m7AyIi4osA4zMzM3qSIKKZg6VDgAxFxbOo+ChgfEcen7nWA64DJEbFE0jzgCxFxe5WypgBTAIYNGzaupaWlkJjbrFy5ksGDBxc6j6L19zo0Kv4Fj62oW1nDBsGTq+pWHGOHD61fYTXq78sRuA59RW/UYcKECfMjotNzRYp8LNEyYGSuewSwPNc9BNgZmCcJ4F+AOZImVSa/iJgJzARoamqK5ubmAsOGefPmUfQ8itbf69Co+Ov5GKHpY1s5e0H9VrElRzTXraxa9fflCFyHvqIv1aHIps7bgNGStpW0HnAYMKdtYESsiIjNI2JURIwCbgbWSHpmZmb1VFjii4hWYBpwFXAfcGlELJR0uqRJRc3XzMysI4U+gT0i5gJzK/qd2s64zUXGYmZmBr5zi5mZlYwTn5mZlYoTn5mZlYoTn5mZlYoTn5mZlYoTn5mZlYoTn5mZlYoTn5mZlYoTn5mZlYoTn5mZlYoTn5mZlYoTn5mZlYoTn5mZlYoTn5mZlYoTn5mZlYoTn5mZlYoTn5mZlYoTn5mZlYoTn5mZlYoTn5mZlYoTn5mZlYoTn5mZlcrARgdga49RM66sW1nTx7YyuY7lLTnrwLqVZWb9m/f4zMysVJz4zMysVJz4zMysVJz4zMysVJz4zMysVJz4zMysVJz4zMysVJz4zMysVJz4zMysVJz4zMysVJz4zMysVJz4zMysVJz4zMysVJz4zMysVJz4zMysVJz4zMysVJz4zMysVApNfJImSrpf0iJJM6oMnyppgaS7JN0oaUyR8ZiZmRWW+CQNAM4BDgDGAIdXSWwXR8TYiHgn8C3gu0XFY2ZmBjUmvpTEumo8sCgiFkfEa0ALcHB+hIh4Ide5IRDdmI+ZmVnNFNF5rpH0MDAbOD8i7q2pYOkQYGJEHJu6jwL2iIhpFeN9BjgRWA/YNyIerFLWFGAKwLBhw8a1tLTUEkK3rVy5ksGDBxc6j6I1og4LHltRt7KGDYInV9WtOMYOH1rTeGtDHerJ60Lf4DrUZsKECfMjoqmz8QbWWN4uwGHAuZLWAc4DWir22CqpSr81smxEnAOcI+njwCnA0VXGmQnMBGhqaorm5uYaw+6eefPmUfQ8itaIOkyecWXdypo+tpWzF9S6eHZuyRHNNY23NtShnrwu9A2uQ33V1NQZES9GxM8i4j3AF4GvAo9LukDSDu1MtgwYmeseASzvYDYtwIdqicfMzKy7aj7GJ2mSpMuAHwBnA9sBlwNz25nsNmC0pG0lrUe2xzinotzRuc4DgTWaOc3MzOqp1naYB4HrgW9HxE25/rMl7VNtgoholTQNuAoYAJwXEQslnQ7cHhFzgGmS9gdeB56jSjOnmZlZPdWa+D4RETfme0jaKyL+EhEntDdRRMylYo8wIk7Nff5sV4I1MzPrqVqv4/thlX4/qmcgZmZmvaHDPT5J7wbeA2wh6cTcoI3Imi/NzMz6lc6aOtcDBqfxhuT6vwAcUlRQZmZmRekw8UXEDcANkmZFxCO9FJOZmVlhOmvq/H5EfA74L0nVLj6fVFhkZmZmBeisqfMX6f07RQdiZmbWGzpr6pyf3m/onXDMzMyK1VlT5wI6eGJCROxS94jMzMwK1FlT50G9EoWZmVkv6ayp02dympnZWqXDO7dIujG9vyjphcr33gnRzMysfjrb49s7vQ/paDwzM7P+ouanZEraHdib7GSXGyPizsKiMjMzK0itz+M7FbgA2AzYHJgl6ZQiAzMzMytCrXt8hwO7RcQrAJLOAu4AziwqMDMzsyLU+liiJcD6ue63AQ/VPRozM7OCdXYB+4/Ijum9CiyUdE3qfj9wY0fTmpmZ9UWdNXXent7nA5fl+s8rJBozM7OCdXY5wwW9FYiZmVlvqOnkFkmjgW8AY8gd64uI7QqKy8zMrBC1ntxyPvDfQCswAbiQtx5ZZGZm1m/UmvgGRcS1gCLikYj4GrBvcWGZmZkVo9br+F6RtA7woKRpwGPAlsWFZWZmVoxa9/g+B2wAnACMA44Cji4qKDMzs6LUtMcXEbcBpL2+EyLixUKjMjMzK0it9+psSk9jvwdYIOluSeOKDc3MzKz+aj3Gdx5wXET8GUDS3mRneu5SVGBmZmZFqPUY34ttSQ8gIm4E3NxpZmb9Tmf36tw9fbxV0k+BS8ju1fkxfNsyMzPrhzpr6jy7ovuruc9R51jMzMwK19m9Oif0ViBmZma9odazOodK+q6k29PrbElDiw7OzMys3mo9ueU8spNZPppeL5Cd1WlmZtav1Ho5w/YR8ZFc92mS7ioiIDMzsyLVuse3Kl27B4CkvYBVxYRkZmZWnFr3+KYCF+aO6z2H79VpZmb9UKeJL92fc8eI2FXSRgAR8ULhkZmZmRWg06bOiHgDmJY+v+CkZ2Zm/Vmtx/iukfQFSSMlbdr2KjQyMzOzAtSa+D4FHAfcANyee3VI0kRJ90taJGlGleEnSrpX0j2SrpW0TVeCNzMz66paE98Y4BzgbuAu4EfATh1NIGlAmuaANP3hksZUjHYn0BQRuwCzgW/VHrqZmVnX1Zr4LgDeAfyQLOm9I/XryHhgUUQsjojXgBbg4PwIEXF9RLycOm8GRtQauJmZWXcoovN7TUu6OyJ27axfxfBDgIkRcWzqPgrYIyKmtTP+fwFPRMSZVYZNAaYADBs2bFxLS0unMffEypUrGTx4cKHzKFoj6rDgsRV1K2vYIHiyjleKjh1e2x321oY61JPXhb7BdajNhAkT5kdEU2fj1Xod352S9oyImwEk7QH8pZNpVKVf1Swr6UigCXhfteERMROYCdDU1BTNzc01ht098+bNo+h5FK0RdZg848q6lTV9bCtnL6h18ezckiOaaxpvbahDPXld6Btch/qqda3cA/iEpKWpe2vgPkkLgEjH6CotA0bmukcAyytHkrQ/cDLwvoh4tebIzczMuqHWxDexG2XfBoyWtC3wGHAY8PH8CJJ2A35K1iT6VDfmYWZm1iU1Jb6IeKSrBUdEq6RpwFXAAOC8iFgo6XTg9oiYA3wbGAz8WhLA0oiY1NV5mZmZ1ap+ByCqiIi5wNyKfqfmPu9f5PzNzMwq1Xo5g5mZ2VrBic/MzErFic/MzErFic/MzErFic/MzErFic/MzErFic/MzErFic/MzErFic/MzErFic/MzErFic/MzErFic/MzErFic/MzErFic/MzErFic/MzErFic/MzErFic/MzEql0Cewm1nvGzXjyrqVNX1sK5PrWN6Ssw6sW1lm3eU9PjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzK5VCE5+kiZLul7RI0owqw/eRdIekVkmHFBmLmZkZFJj4JA0AzgEOAMYAh0saUzHaUmAycHFRcZiZmeUNLLDs8cCiiFgMIKkFOBi4t22EiFiShr1RYBxmZmZvKrKpczjwaK57WepnZmbWMIqIYgqWDgU+EBHHpu6jgPERcXyVcWcBV0TE7HbKmgJMARg2bNi4lpaWQmJus3LlSgYPHlzoPIrWiDoseGxF3coaNgieXFW34hg7fGhN47kOq2tUHerJ63Pf0Bt1mDBhwvyIaOpsvCKbOpcBI3PdI4Dl3SkoImYCMwGampqiubm5x8F1ZN68eRQ9j6I1og6TZ1xZt7Kmj23l7AX1WzyXHNFc03iuw+oaVYd68vrcN/SlOhTZ1HkbMFrStpLWAw4D5hQ4PzMzs04VlvgiohWYBlwF3AdcGhELJZ0uaRKApHdJWgYcCvxU0sKi4jEzM4NimzqJiLnA3Ip+p+Y+30bWBGpmZtYrfOcWMzMrFSc+MzMrFSc+MzMrFSc+MzMrFSc+MzMrFSc+MzMrFSc+MzMrFSc+MzMrFSc+MzMrFSc+MzMrFSc+MzMrFSc+MzMrFSc+MzMrFSc+MzMrFSc+MzMrFSc+MzMrFSc+MzMrlUKfwG61GzXjyrqWN31sK5PrVOaSsw6sSzlmZn2B9/jMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxUnPjMzKxU/CBaM+tz6vlg5no+lBn8YOa1gff4zMysVJz4zMysVJz4zMysVJz4zMysVApNfJImSrpf0iJJM6oMf5ukX6Xht0gaVWQ8ZmZmhSU+SQOAc4ADgDHA4ZLGVIx2DPBcROwAfA/4ZlHxmJmZQbF7fOOBRRGxOCJeA1qAgyvGORi4IH2eDewnSQXGZGZmJVfkdXzDgUdz3cuAPdobJyJaJa0ANgOe6cqM6nnND9T3uh9f82NWTr4Wse9SRBRTsHQo8IGIODZ1HwWMj4jjc+MsTOMsS90PpXGerShrCjAlde4I3F9I0G/ZnC4m3z6ov9ehv8cPrkNf4Tr0Db1Rh20iYovORipyj28ZMDLXPQJY3s44yyQNBIYC/6gsKCJmAjMLinMNkm6PiKbeml8R+nsd+nv84Dr0Fa5D39CX6lDkMb7bgNGStpW0HnAYMKdinDnA0enzIcB1UdQuqJmZGQXu8aVjdtOAq4ABwHkRsVDS6cDtETEH+DnwC0mLyPb0DisqHjMzMyj4JtURMReYW9Hv1NznV4BDi4yhm3qtWbVA/b0O/T1+cB36Ctehb+gzdSjs5BYzM7O+yLcsMzOzUnHiSyStL+lWSXdLWijptEbH1F2SBki6U9IVjY6lOyQtkbRA0l2Sbm90PN0haWNJsyX9XdJ9kt7d6Ji6QtKO6ftve70g6XONjqurJH0+rc9/k3SJpPUbHVNXSfpsin9hf/kNJJ0n6SlJf8v121TSNZIeTO+bNCo+J763vArsGxG7Au8EJkras8ExdddngfsaHUQPTYiId/aV05+74QfAHyLiX4Fd6We/R0Tcn77/dwLjgJeByxocVpdIGg6cADRFxM5kJ9n1qxPoJO0MfJrsTli7AgdJGt3YqGoyC5hY0W8GcG1EjAauTd0N4cSXRGZl6lw3vfrdAVBJI4ADgXMbHUtZSdoI2IfsrGUi4rWIeL6xUfXIfsBDEfFIowPphoHAoHSd8AaseS1xX/cO4OaIeDkiWoEbgA83OKZORcSfWPOa7PwtKi8APtSrQeU48eWkJsK7gKeAayLilkbH1A3fB74IvNHoQHoggKslzU937elvtgOeBs5PTc7nStqw0UH1wGHAJY0Ooqsi4jHgO8BS4HFgRURc3diouuxvwD6SNpO0AfBBVr8xSH8yLCIeB0jvWzYqECe+nIj4Z2raGQGMT80M/Yakg4CnImJ+o2Ppob0iYneyJ3t8RtI+jQ6oiwYCuwP/HRG7AS/RwGadnkg3n5gE/LrRsXRVOoZ0MLAtsBWwoaQjGxtV10TEfWRPrbkG+ANwN9Da0KDWAk58VaRmqXms2Ubd1+0FTJK0hOxpGPtK+mVjQ+q6iFie3p8iO640vrERddkyYFmuxWA2WSLsjw4A7oiIJxsdSDfsDzwcEU9HxOvAb4H3NDimLouIn0fE7hGxD1nz4YONjqmbnpT0doD0/lSjAnHiSyRtIWnj9HkQ2Urz98ZG1TUR8aWIGBERo8iap66LiH61hStpQ0lD2j4D/0bW3NNvRMQTwKOSdky99gPubWBIPXE4/bCZM1kK7Clpg/S4s/3oZycZAUjaMr1vDfwH/ff3yN+i8mjgd40KpNA7t/QzbwcuSA/QXQe4NCL65eUA/dww4LL0WMaBwMUR8YfGhtQtxwMXpabCxcAnGxxPl6VjSu8H/rPRsXRHRNwiaTZwB1nz4J30obuHdMFvJG0GvA58JiKea3RAnZF0CdAMbC5pGfBV4CzgUknHkG2UNOyuXb5zi5mZlYqbOs3MrFSc+MzMrFSc+MzMrFSc+MzMrFSc+MzMrFSc+Mx6maR5kgq/+bakE9KTIS4qel5m/Ymv4zPrRyQNTDcrrsVxwAER8XCRMZn1N97jM6tC0qi0t/Sz9By0q9MdfVbbY5O0ebpFHJImS/ofSZdLeljSNEknphtV3yxp09wsjpR0U3rO2vg0/YbpOWa3pWkOzpX7a0mXA2vcZDnN42/p9bnU7ydkN8ueI+nzFePvpOzZk3dJukfS6E7q++kU092SfpMubEfSLEk/kfRnSQ+ke8VWLb+OP41ZjznxmbVvNHBOROwEPA98pIZpdgY+TnZ/0a8DL6cbVf8V+ERuvA0j4j1ke2XnpX4uNVb9AAACMUlEQVQnk91m7l3ABODbuac6vBs4OiL2zc9M0jiyu8LsAewJfFrSbhExlewRPBMi4nsVMU4FfpBuyN5Edm/Rjur724h4V3pW5X3AMbmyRgHvI3sU1k+UPei1vfLN+gQ3dZq17+GIuCt9nk/2J9+Z6yPiReBFSSuAy1P/BcAuufEugey5ZZI2SveJ/Teym4x/IY2zPrB1+nxNRFQ+3wxgb+CyiHgJQNJvgfeS3Z6rPX8FTk7PbvxtRDyYbhHXXn13lnQmsDEwGLgqV9alEfEG8KCkxcC/Viu/g1jMep33+Mza92ru8z95a0OxlbfWnfU7mOaNXPcbrL6hWXmvwAAEfKTtyecRsXV6LA1kjzaqRh1XYU0RcTHZo4ZWAVdJatuLbK++s4BpETEWOI3V67xGPToo36xPcOIz67olwLj0+ZBulvExAEl7kz0gdQXZntTx6UkCSNqthnL+BHwoPYFgQ7Knc/+5owkkbQcsjogfkt0xf5eOxgeGAI9LWhc4omLYoZLWkbQ92THF+7tRvlmvclOnWdd9h+wu80cB13WzjOck3QRsBHwq9TsD+D5wT0p+S4CDOiokIu6QNAu4NfU6NyI6auaELOkeKel14Ang9BRHe74C3AI8QtZkOyQ37H7gBrKnakyNiFckVSvfrM/w0xnMrFtSwr0iImY3OhazrnBTp5mZlYr3+MzMrFS8x2dmZqXixGdmZqXixGdmZqXixGdmZqXixGdmZqXixGdmZqXyv3VnOTGI0ucVAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.bar( range(64), P[range(64),63] )\n",
    "plt.xlim(2.5,10.5)\n",
    "plt.ylim(0,0.5)\n",
    "plt.xlabel('number of snaps')\n",
    "plt.ylabel('probability')\n",
    "plt.title('Distribution of number of Thanos snaps when N=63')\n",
    "plt.grid()\n",
    "plt.tight_layout(pad=0)\n",
    "plt.savefig('thanos1.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute expected survival\n",
    "\n",
    "If there are `s` snaps, then there will be `binom(pop,p**s)` survivors. Exected value of this is `pop*p**s`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 325,
   "metadata": {},
   "outputs": [],
   "source": [
    "# expected survival rate\n",
    "x = np.arange(Nmax+1).reshape((-1,1))\n",
    "exprate = np.dot( P.T, p ** x )\n",
    "var = np.dot( P.T, ( p ** x - exprate ) ** 2 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 347,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfMAAAEuCAYAAACar/u6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XecXGXZ//HPtb1vymY3yaZ3kgABNgGlJYBKFRQVsKOiKCo86mN7fgr6WFEeAUWlIwakqhQVkJJQpCSBkNASSO+9bfruXr8/ztlkMuzszmbP7pT9vl+vee3MnDNnrrn2nLnmPuc+9zF3R0RERDJXTqoDEBERkY5RMRcREclwKuYiIiIZTsVcREQkw6mYi4iIZDgVcxERkQynYp5CZvZZM3v2IF97hZlNTXLexWZ2ykG+zzQz+8LBvDYqZvYJM3ssguUcdL47g5n90cx+EPP4y2a2xszqzax3F7y/m9mIzn4faR8zO97M5iU57wHrdLb9Tzv43fW6mU1ux/zFZvaQmW0xs3uj+t7pKhlVzMN/7M7wy6759rsUxpPyQtcduPsd7v7+VMcRz8xuM7OfHOzr3f1id//fcFn5wP8B73f3MnffEPdeQ8Iv6n/EPT/VzK5oIbbvx2wju8ysMebx6wcbs0QvvgC7+zPuPjqVMaVCR7eneO4+zt2nteMlHwFqgN7u/tF0/d5JJKOKeeis8Muu+fbVVAckHWNmeamOIV4KYqoBioC2Cu0xZnZsWwtz9581byPAxcDzMdvMuAjiFclYCbbvwcB8d2/o6nhiHex3TyYW8xaZ2R/M7L6Yx780sycsMNnMloetlfVhC/8TMfMWmtmvzWxpuJvzj2ZWHDP9bDObbWZbzWyBmZ1qZj8Fjgd+F7uHwMzGmNm/zWyjmc0zs4/FLKe3mT0YLuclYHgbn+nTZrbEzDaY2Q9a2+VkZh8MdyttDvcYHBI3y0Qze8PMNpnZrWZWFL6up5k9bGbrwmkPm9mAJHM+ycxmhp9njZn9X/j8ZDNbHjfvvtgtOERwX9iq3Ap8P9zj0itm/iPC/1V+7K7E8H/z67hlP2Bm3wjvfzf8H20LP++HkvwszS3fz5vZUuDJ8Pl7zWy1BbvenjazceHzXwQ+AXw7/P8/FD7f38zuD/O5yMy+3sp73mZmPzGzUUDzbtXNZvZkK6FeCUTWegFOMbO3w//9dWZmYWzDzezJcN1bb2Z3mFmPmNgXm9m3zGxOmJu7m9epcPpFZvZOuB08aGb9w+fNzH5jZmvD180xs/HhtITboZlVhevm5nCZz5jZu76/Dnb54fT/NrNVZrbSzD5nMS1mi9sLZ+/evd3adn9bmNt/hOvli2Y2PJz2dDjbq+F6dF789nMw67SZTQw/Y17Mc+ea2ewE858eLnubma0ws2+Fzzd/d34zzOkqM7sw5nVnmNkrFnwHLLOYvUS2f5v6YpjTVWb2zQTv3+L2FJrQynp2pgXfzZvN7D9mdljMtNa+cz4b9/4/An4InBe+/+db+B+7mV1sLW8vuWZ2lQXbyiIz+2o4f144vdLMbg5zsMKC7T43nPZZM3suXG83AleY2Qgzmx5+5vVmdneL/+hY7p4xN2AxcEqCaSXA/PCfdDywHhgQTpsMNBDsxiwETgS2A6PD6VcDDwK9gHLgIeDn4bRJwBbgfQQ/fmqBMeG0acAXYmIoBZYBFwJ5wJFhHOPC6XcB94TzjQdWAM8m+DxjgXrgOKAA+DWwt/nzA1cAU8P7o8LP8z4gH/g28A5QEJO314CB4Wd8DvhJOK03cG6Yv3LgXuDvMXEc8BnjYnwe+FR4vww4JibfyxP978LY9wLnhDktJiieF8XM/yvgj+H9zzbnCTghzLGFj3sCO4H+4eOPAv3D5Z4X5qVf/HJa+CxDAAduD/8/xeHznwvzUkiwnsyOec1tzXkMH+cAswi+FAqAYcBC4AMJ3nPf62PeP6+N+MoI1pvmXE4Frmhju2nxc4fLexjoAQwC1gGnhtNGEKxPhUAf4Gng6rj/50thrnsBbwIXh9NOIljvjwxf/1vg6XDaB8Ic9QAMOCTm/9Padvhz4I8E63c+wTZuLXymg13+qcAagu2yFLgzzM+IBNv6vpzS9nZ/G7CR4LskD7gDuCvu/zAi5vFkYrYf2rFOx8X8BnBazLS/Ad9MsI6sAo6P2aaOjPvu/HGY99OBHUDPmOmHhrEdFubwnLh19i9hjg4lWMcSfYffRsz2lMR6diSwFjgayAU+E85fmMx3TgvvfwXhd2oruU20vVwc5ntAmL/Hidmegb8D14d5qA4/05di3qcB+BrB+lEc5ux/wliLgONa28bdPSOLeT2wOeYWWwAmEWw0S4AL4jaOBqA05rl7gB8QbPDbgeEx094DLArvXw/8JkE80zhwAz8PeCZunuuBy8OVbS/hD4Fw2s9IXFx+CPwl5nEJsIeWi/kPgHti5s0h+MKfHJO3i2Omnw4sSPC+E4BNiT5j3LxPAz8CquKen0zbxfzpuOlfAJ4M7xvBl+MJ8RtVOG1pzLSLml+XIMbZwNktbZxx8w0h2PiGtbKsHuE8leHj2ziwmB8NLI17zfeAWxMsb9/rSb6Y5wFfAV4In+9oMT8u5vE9wHcTLOMc4JW4/+cnYx5fyf4fXzcDV8ZMKyNY94cQFPr5wDFATsw8bW2HPwYeIKboJYjzYJd/C/CLmGmjSL6YJ9zuY/7PN8Vtf2/F/R8SFvP2rNNxMX8HuCO834ugCPdLsMylwJeAiha25Z2x6yVBAT0mwXKuJvy+ZP86G/uddyVwc1vbQ5Lr2R+A/42bfx5wYsxrE37ntPD+V9B2MW9xeyFojHwpZtop7N9ea4DdxPyAAC4Anop5n/jvjduBGwgbpMncMnE3+znu3iPmdmPzBHd/iaAlZASJjrXJ3bfHPF5C8GuvD0GhnBXuqtkMPBI+D0FrdkGSsQ0Gjm5eTrisTwB9w+XlERSp2BgS6R87r7vvADa0Mu+SmHmbwtfWxswT/77Nuz1LzOx6C3bnbyUo0D2adwG14fMEX3pvmdkMMzszide0FA/AfcB7LNgdewLBhvBM/Is8WNPvItgYAD5O0NIh/DyfjtnttpmgpVV1MHGFu85+Ee7i3Erw5UAryxsM9I/7/3+fYGOO0o1AjZmdFcGyVsfc30FQeDGzajO7K9wluJXgR0P8527xtbx7fawnWHdr3f1J4HfAdcAaM7vBzCpoezv8FcHepsfMbKGZfbelD9OB5R+wvdH6thmvte2+WaJctakD6/RU4CwzKwM+RvCDY1WCec8l+JGxJNy9+56YaRv8wOPIsevJ0Wb2lAWHlbYQtFDjY2vxu6cdEuVuMPDNuLwPbGX58d85B6O1dT52+bH3BxPs1VgVE+f1BC30RLF9m6COvWTB4dPPtRVYJhbzhMzsEoLdeisJkhGrp5mVxjweFM63nuCX57iYHwiVHnQcgiDJiY5te9zjZcD0uB8bZe7+ZYJdMg0EK1tsDImsIthl0/zZigl2ibdkJcEK0zyvhe+zImae+PddGd7/JjAaONrdKwgKKQQrUqvc/W13v4BgpfwlcF+Y4+0EX5zN8eSy/0tz38vjlrUZeIzgS+fjBHsl4vPb7C/AR8xsMEFr+P7wfQYTFLqvEvRI7UFweKHNz5Igro8DZxP8yq4kaGkQs7yW/v+L4v7/5e5+ejvev+0A3fcS7BH5X9r32drj5wSf77BwvfhkO94rfn0sJVh3VwC4+7XufhQwjuDH4H/Txnbo7tvc/ZvuPgw4C/iGmZ3c0psfzPIJtrfWts0D1mkOLNStbfcd0pF12t1XEBwK+xDwKeDPrcw7w93PJtiW/867G0OJ3Elw6GKgu1cSHAqJjy3Rd8+7wkjyPZstA34al/cSd/9LRMtvjwO+rznwMy8jaJlXxcRZ4Qd2RI3/Plzt7he5e3+CPSa/tzZOOcyaYm5BJ6KfEHzpfIqgI8WEuNl+ZGYFZnY8cCZwb9iKvRH4jZlVh8uqNbMPhK+5GbjQzE42s5xw2phw2hqC46LNHgZGmdmnLOi4lW9BR5RD3L0R+CtB54YSMxtLcIwnkfsIflW/18wKCL68E23A9wBnhDHmExTo3cB/Yua5xMwGWNDJ7PtAc4eKcoIvuc3htMtbiekAZvZJM+sT5nBz+HQjwW7OIgs6x+QD/4/gR1Zb7gQ+TdBKuDPRTO7+CsGPo5uAR8MfAhAcj/JwGhZ01Bmf7OdpQTlBHjcQfJH/LG56/P//JWCrmX3HgnNWc81svJlN7EAMifyZIKendsKyIfjs9QTrRS1BQUzWnQTbzAQzKyTI24vuvjjcHo4O14vtwC6gsa3t0IKOTiPCH6pbCdazxvg3PtjlE2xDnzWzsWZWwru3g9nAh8NtdwTBXqlmCbf7JPMVvx7F6ug6fTtBw+ZQgmPm7xJ+J37CzCrDH4rN+U1GObDR3XeZ2SSCH8DxfhDmbRxBv4JEnblay0NLbgQuDv/fZmal4XdOeTuWEZV7gEvDdaoHwSEOAMK9IY8BV5lZRVhHhpvZiYkWZmYftf0dkTcRrAOt/k8ysZg/ZAeeZ/43C3oMTgV+6e6vuvvbBAXrz+GXCQS7RzYR/Cq8g+AY8lvhtO8Q7MJ7wYJdio8TtFabd91fCPyGoCPcdPa3Oq4haCFuMrNr3X0b8H7g/PB9VhO0WJtj+CrBbpnVBMeHbk30Id39dYIOEXcR/OrbRnCsancL884j+BHzW4IWyFkEp/DtiZntToIVamF4a+4RfTVBh4v1wAsEux6TdSrwupnVE+TifHff5e5bCI7r3kTQGtsOLE+8mH0eBEYCa9z91Tbm/QtBi3lf0Xf3N4CrCFojawi+wJ5rx+eJdzvBbsEVBJ1bXoibfjMwNtx19vfwB9tZBP0OFhHk9CaCVn2kwve6nOBYaGf4EUEHoy3APwh+iCbF3Z8g6MdxP8G6O5xgmwCoIPgS3kSQ2w0EnTuhle2QYL14nOAHxvPA773lc4gPavnu/i+CbeHJcJ74Mwp+Q9BnZQ3wJ2IO7SSx3bflCuBP4Xr0sdgJEazTfyP4vvpb3GHGeJ8CFod5uZjg+yQZXwF+bGbbCPr5tNSin06Q0yeAX7t7ooFYDtie2npjd59J0GfmdwT/73eI66XehW4k+H6dA7wC/JNgT2xzAf40QafYNwhivQ/o18ryJgIvht+tDwKXuvui1gJo7hGc1SwYBWiquyd1ylU6suC412ZgZFv/VBHpODNzgu3tnVTH0hFmtoCgc9bjXfy+Qwh+1OZ7is/d7mpmdhpBR73Bbc4ckUxsmXcbZnZWuHuqlKB1MZf9nbBERFplZucS7KJtbewC6aDwsNrpZpYXHpa6nASHNTqLinl6O5tgt91Kgt2M57fSKUxEZB8zm0Zw+tYlYZ8B6TxGcGhqE8Fu9jcJDjt0XQCqDSIiIplNLXMREZEMp2IuIiKS4dLualVRqKqq8iFDhiScvn37dkpLSxNOl/ZRPqOlfEZL+YyW8hm9WbNmrXf3+IG12iUri/mQIUOYOXNmwunTpk1j8uTJXRdQllM+o6V8Rkv5jJbyGT0za8/wwS3SbnYREZEMp2IuIiKS4VTMRUREMpyKuYiISIZTMRcREclwad+bPRyX/PcEVyya5u53tPESERGRbiUlLXMzu8XM1prZa3HPn2pm88zsHTP7bvj0h4H73P0i4INdHqyIiEiaS9Vu9tsIroW9j5nlAtcBpwFjgQvMbCwwAFgWztbqxdmT1aTx6EVEJIukpJi7+9PAxrinJwHvuPtCd98D3EVw1bDlBAUdIoj36395hf95dmdHFyMiIpI20umYeS37W+AQFPGjgWuB35nZGcBDiV5sZl8EvghQU1PDtGnTWpxv0/rd7NjblHC6tF99fb3yGSHlM1rKZ7SUz/SUTsXcWnjO3X07cGFbL3b3G4AbAOrq6jzRcIPP1L/BS6sXaTjCCGl4x2gpn9FSPqOlfKandDo1bTkwMObxAGBl1G9SXpTH7kZoaGyKetEiIiIpkU7FfAYw0syGmlkBcD7wYNRvUlYY7IzYvjuSvnQiIiIpl5Ld7Gb2F2AyUGVmy4HL3f1mM/sq8CiQC9zi7q+3c7lnAWfV1tYmPKazcvleAB6b9gx9StLpt0zm0jG0aCmf0VI+o6V8pifzLDxNq66uzhNdAvVfc1fx5Tte5l+XHs8h/Sq6OLLspGNo0VI+o6V8Rkv5jJ6ZzXL3uo4so9s1TcuKgp0R23Y1pDgSERGRaHS7Yl5elA/Atl17UxyJiIhINLphMQ9a5vW71TIXEZHskE7nmXdYMh3gNu8KTkmbNecNKje/3XXBZTF1iImW8hkt5TNaymd6yqpi7u4PAQ/V1dVdlKiDxs49jTDtEfoOGsrkySO6NL5spQ4x0VI+o6V8Rkv5TE/dbjd7UX4OOQb16gAnIiJZotsVczOjOE+92UVEJHt0u2IOUJJn6s0uIiJZo1sW86I8U292ERHJGlk1AlxMb/aLpk6dmnC+n/ynntzcXL53dHHXBZfF6uvrKSsrS3UYWUP5jJbyGS3lM3pTpkzp8AhwWVXMm7U2nCvAOVc9wp68Uv556fFdGFX2Uu/WaCmf0VI+o6V8Rk/DuR6k4jzYtlvHzEVEJDt0z2Keb+rNLiIiWaNbFvOSPKN+VwPZeIhBRES6n25ZzIvyoKHJ2bW3KdWhiIiIdFhWDeeazNjsALmNewDj0Sen06OoW/6eiZTGao6W8hkt5TNaymd6yqpinszY7AD/Wfk4sJtDj5rE8D46xaKj1Ls1WspntJTPaCmf6albNktLwp8w6gQnIiLZoFsW8+I8A9CQriIikhW6aTEP/urKaSIikg26aTFvbpmrmIuISObrlsW8JD8s5rrYioiIZIFuWcyLcoO/OmYuIiLZIKtOTUv2PPOdO7ZTmm/MmbeIaXkruyy+bKXzTqOlfEZL+YyW8pmeuuVV06ZNm8ZVc/PoUZLPnz9/dBdGlp103mm0lM9oKZ/RUj6jp6umdcDg3iUs2bAj1WGIiIh0WLcu5is272Rvo8ZnFxGRzNZ9i3mvUhqbnJWbd6Y6FBERkQ7pvsW8dwmAdrWLiEjG68bFvBSAJRu2pzgSERGRjum2xby6vJDCvBy1zEVEJON122Kek2MM6lXCko0q5iIiktm6bTGHYFf7UrXMRUQkw2XVoDExI8BdNHXq1ITz1dfXU1ZWxl/e3M1Tyxq4/n0lmFnXBZplmvMp0VA+o6V8Rkv5jN6UKVM6PGhMVg3n6u4PAQ/V1dVd1NoIRc0jGC0tXMyjS15n7FHvoaaiqKvCzDoaESpayme0lM9oKZ/pqVvvZh8S9mhfuE492kVEJHN162I+sibYVfT22m0pjkREROTgdeti3reiiPKiPOavUTEXEZHM1a2LuZkxqqac+WvqUx2KiIjIQevWxRxgVE0Z89dsI5t69YuISPeiYl5TzuYde1lXvzvVoYiIiBwUFfOacgDe1q52ERHJUCrmYTGft1qd4EREJDN1+2JeVVZAz5J8nZ4mIiIZq9sX8+Ye7WqZi4hIpurWY7M3+/Mbu3luRQO/P6WEHI3R3m4aqzlayme0lM9oKZ/R09jscdo7Nnuz9eXLeeLeVxk0ro4R1eWdHme20VjN0VI+o6V8Rkv5TE/dfjc7wOEDKgF4ddmWFEciIiLSfirmwLA+ZZQW5DJn+eZUhyIiItJubRZzMysxsx+Y2Y3h45Fmdmbnh9Z1cnOM8bWVvLpcLXMREck8ybTMbwV2A+8JHy8HftJpEaXI4QN78MaqrexpaEp1KCIiIu2STDEf7u5XAnsB3H0nkHVdvg8bUMmehiZdQU1ERDJOMsV8j5kVAw5gZsMJWupZ5fABPQB4VcfNRUQkwyRTzK8AHgEGmtkdwBPAdzozqFQY0LOYXqUFvLJUxVxERDJLm+eZu/tjZjYLOIZg9/ql7r6+0yPrYmbGxCE9eXHRhlSHIiIi0i7J9GZ/wt03uPs/3P1hd19vZk90RXBd7eihvVm2cScrNu9MdSgiIiJJS1jMzazIzHoBVWbW08x6hbchQP+uCrArHTOsNwAvLlTrXEREMkdrLfMvAbOAMeHf5tsDwHWdH1rXG9O3nMrifF5cuDHVoYiIiCQt4TFzd78GuMbMvubuv+3CmFImJ8eYOKSXjpuLiEhGSaYD3G/NbDwwFiiKef72zgwsVY4Z1ovH31zD6i276FtZ1PYLREREUiyZDnCXA78Nb1OAK4EPdnJcKfOe4cFx82feXpfiSERERJKTzHnmHwFOBla7+4XA4UBhp0aVQmP7VVBdXsi0+SrmIiKSGczdW5/B7CV3nxSeaz4F2Aa85u7juiLA9jCzs4CzamtrL5o6dWrC+err6ykrK0s4/ea5u5m5poHfnVRCbk7WjVwbubbyKe2jfEZL+YyW8hm9KVOmzHL3uo4so81j5sBMM+sB3EjQm70eeKkjb9pZ3P0h4KG6urqLJk+enHC+adOm0dr0nb1X8cwdL1M+9HAmDe0VeZzZpq18Svson9FSPqOlfKanVou5mRnwc3ffDPzRzB4BKtx9TpdElyLHjqwiL8d4at5aFXMREUl7rR4z92Af/N9jHi/O9kIOUFGUT92Qnjz11tpUhyIiItKmZDrAvWBmEzs9kjRzyiE1vLV6G0s2bE91KCIiIq1KpphPAZ43swVmNsfM5ppZ1rfOTx3fF4B/zF2V4khERERal0wHuNM6PYo0NKBnCYcP7MG/5q7mK5NHpDocERGRhNpsmbv7kpZuXRFcqp1xaF/mrtjC0g07Uh2KiIhIQsnsZu+2ThvfD4B/vqZd7SIikr5UzFsxsFcJEwb24O+vrEh1KCIiIgmpmLfh3CNreWv1Nl5fuSXVoYiIiLQoYTE3s21mtrWF2zYz29qVQabSWYf3Jz/XuH+WWuciIpKeEhZzdy9394oWbuXuXtGVQaZSj5ICTh5TwwOzV7C3sSnV4YiIiLxL0rvZzazazAY13zozqHRz7lED2LB9D09qRDgREUlDyVzP/INm9jawCJgOLAb+1clxpZUpo/tQU1HIHS8uTXUoIiIi75JMy/x/gWOA+e4+lODa5s91alRpJi83hwsmDeLp+es0vKuIiKSdZIr5XnffAOSYWY67PwVM6OS40s4FkwaRm2Pcqda5iIikmWSK+WYzKwOeBu4ws2uAhs4NK/3UVBTxgXE13D1zGTv3NKY6HBERkX2SKeZnAzuA/wIeARYAZ3VmUOnqs+8dyuYde7n/5eWpDkVERGSfZIr5F4H+7t7g7n9y92vD3e7dzsQhPTl8YA9ufnYRTU2e6nBERESA5Ip5BfComT1jZpeYWU1nB5WuzIyLjh/KovXbefzNNakOR0REBEjuqmk/cvdxwCVAf2C6mT3e6ZGlqVPH9WVAz2L+MH0B7mqdi4hI6rVnbPa1wGpgA1DdOeGkv7zcHL48eTivLN3Mc+90y6MNIiKSZpIZNObLZjYNeAKoAi5y98M6O7B09pGjBtCvsohrnpiv1rmIiKRcMi3zwcBl7j7O3S939zc6O6h0V5iXy8UnDmfG4k08v0CtcxERSa3WrprWfDGVK4GlZtYr9tY14aWv8yYOpH9lEb945C31bBcRkZRqrWV+Z/h3FjAz/Dsr5nGXMLNhZnazmd3XVe+ZjKL8XL71gdHMWb6Fh+asTHU4IiLSjbV2CdQzw79D3X1Y+Lf5NiyZhZvZLWa21sxei3v+VDObZ2bvmNl3W1uGuy90988n835d7ZwJtYztV8GvHp3H7gaNCiciIqmRTAe4B8zsAjMrOYjl3wacGre8XOA64DRgLHCBmY01s0PN7OG4W1r3ms/JMb5/+iEs37ST2/+zJNXhiIhIN2Vt9cY2sxOB84AzgJeAu4GH3X1XUm9gNiScf3z4+D3AFe7+gfDx9wDc/edtLOc+d/9IK9O/SDBaHTU1NUfdddddCZdVX19PWVlZMuEn5aqZu1iwuZErTyihrMAiW26miDqf3Z3yGS3lM1rKZ/SmTJkyy93rOrKMvLZmcPfpBAPF5AInARcBtxCMDHcwaoFlMY+XA0cnmtnMegM/BY4ws+8lKvrufgNwA0BdXZ1Pnjw5YQDTpk2jtent1XfMVk6/5hle3l3DD98/NrLlZoqo89ndKZ/RUj6jpXympzaLOYCZFRNcXOU84EjgTx14z5aargl3D4TjwF/cgffrdGP6VnDexEH86fnFnHtULeP6V6Y6JBER6UaSOWZ+N/AmQav8OmC4u3+tA++5HBgY83gAkPHdwb976hh6luTzvb/OpVGnqomISBdqtZibWQ4wl6CAX+zuT7p7UwffcwYw0syGmlkBcD7wYAeXmXKVJfn84MyxzFm+hdufX5zqcEREpBtJpgPc8+7+noNauNlfgMkEw8CuAS5395vN7HTgaiAXuMXdf3owy2/h/c4Czqqtrb1o6tSpCefrrA4c7s5Vs3bzzqZGfnpcMb2L2zP0feZSh5hoKZ/RUj6jpXxGL4oOcMkU8x8Bc4C/eoYMRF5XV+czZyYe16YzO3As27iD9/1mOseN6MONnz4Ks+zv3a4OMdFSPqOlfEZL+YyemXW4mCfTdPwGcC+w28y2mtk2M9vakTfNZgN7lfBfp4zi8TfX8LdXVqQ6HBER6QaSuZ55ubvnuHuBu1eEjw/2tLRu4QvHD2PikJ788IHXWbphR6rDERGRLJdMb/YTWrp1RXCZKjfH+M15EzDgsrtfoaGxo30GRUREEkvmmPlDMQ+LgEnALHc/qTMDOxip7gAX74WVDfxxzm4+NCKfs0cUdPr7pYo6xERL+YyW8hkt5TN6XdIB7l0vMBsIXOnuF3TkjTtTKjvAxfuvu2fz4KsruedLx3DU4Oy8cqw6xERL+YyW8hkt5TN6XdUBLt5yYHxH3rQ7+dHZ46jtUcxX7niZddt2pzocERHJQskcM/+tmV0b3n4HPAO82vmhZYeKonz++Mmj2LJzL5fc+TJ7dfxcREQilkzLfCYwK7w9D3zH3T/ZqVFlmbH9K/jFhw/jpUUb+fk/30p1OCIikmWSuWravouqmFlPDhxXXZJ0zhG1vLp8M7c8t4jDB1Zy9oTaVIckIiJZIpne7NOADxIU/tnAOmC6u3+j06Nrp3TrzR6vocn51YxdLNhAGM5HAAAgAElEQVTSxHcmFjGyZ26Xx9AZ1Ls1WspntJTPaCmf0euq4VxfcfcjzOwLwEB3v9zM5rj7YR15486UTr3Z423cvocP//45tuzcy9++cixDqkpTEkeU1Ls1WspntJTPaCmf0euq3ux5ZtYP+BjwcEfeTKBXaQG3XTgJgM/e+hIbt+9JcUQiIpLpkinmPwYeBd5x9xlmNgx4u3PDym5Dqkq56TN1rNyyi4tun8nOPY2pDklERDJYMmOz3+vuh7n7V8LHC9393M4PLbsdNbgXV583gZeXbuLLd8xiT4NOWRMRkYPTPS64naZOP7QfP/vQoUybt45L79IY7iIicnDaPZxrOkv33uyJPLp4L395aw/H9s/j84cWkJNh10BPt3xmOuUzWspntJTP6KVkbPZMkM692RO55vG3+c3j87lg0iB+es54cnIyp6CnYz4zmfIZLeUzWspn9KLozZ5w0Bgza/U8cnf/v468sRzo6yePYE9jI9c9tYC9jU388tzDyM2ggi4iIqnT2ghw5eHf0cBE4MHw8VnA050ZVHdkZnzr/aMpyM3lN4/PZ29jE1d99HDyctWtQUREWpewmLv7jwDM7DHgSHffFj6+Ari3S6LrZsyMS08ZSX6eceUj89jT0MTV50+gMC87RooTEZHOkUyzbxAQO7LJHmBIp0QjAHxl8gh+cOZY/vXaai68dQbbdu1NdUgiIpLGkinmfwZeMrMrzOxy4EXg9s4NSz5/3FD+72OH89KijZx3/Qus3bYr1SGJiEiaSmbQmJ8CFwKbgM3Ahe7+s84OTODDRw7gps/UsXjDds79w394Z219qkMSEZE0lNSpaWZ2HDDS3W81sz5Ambsv6vTo2ilTzzNvy8LNjVz98i72NsGXDy/ksD5tXrm2S2VaPtOd8hkt5TNaymf0uuqqaZcDdcBodx9lZv2Be9392I68cWfKxPPM27J80w4uun0W81Zv5X/OGMvnjh2CpcngMpmYz3SmfEZL+YyW8hm9rrpq2ocIrme+HcDdV7L/tDXpIgN6lnDfxe/hfWNr+N+H3+B7f52r8dxFRARIrpjv8aD57gBmlvkX4M5QpYV5/OETR/G1k0Zw14xlfPKmF1m3bXeqwxIRkRRLppjfY2bXAz3M7CLgceCmzg1LEsnJMb75/tFcc/4EXl2+mTOufYYXFm5IdVgiIpJCyfRm/zVwH3A/wWhwP3T3azs7MGnd2RNq+fslx1JWmMfHb3yB6556h6am7BtnX0RE2tZmMTezX7r7v939v939W+7+bzP7ZVcEJ607pF8FD37tOM44rD+/enQeF942g43b97T9QhERySrJ7GZ/XwvPnRZ1IHJwygrzuPb8CfzknPE8v2ADp1/zDM8v0G53EZHuJGExN7Mvm9lcYIyZzYm5LQLmdl2I0hYz45PHDOavX3kvxQW5fPymF/jpP95g197GVIcmIiJdIOF55mZWCfQEfg58N2bSNnff2AWxtVu2DhrTHrsbnLvn7eHJZQ0MKDMuOqyQwRWde6GWbM5nKiif0VI+o6V8Rq+rBo05Bng95qpp5cBYd3+xI2/cmbJx0Jj2emreWr593xw279jDZaeM4uITh3fa9dG7Qz67kvIZLeUzWspn9Lpq0Jg/ALGDgm8Pn5M0NmV0NY9ddgLvG1vDrx6dx4d//xxvrtqa6rBERKQTJFPMzWOa7+7eRCvXQZf00bO0gOs+fiTXXnAEyzft5KzfPsuvHn1Lx9JFRLJMMsV8oZl93czyw9ulwMLODkyiYWZ88PD+PP6NEznniFque2oBp1+jgWZERLJJMsX8YuC9wApgOXA08MXODEqi17O0gF9/9HCmfv5oGpqc8294gW/d+yrr6zUcrIhIpmtzd7m7rwXO74JYpAscN7KKRy87gWuffJubnlnIo6+v5hvvG8WnjhlMXm4yv+1ERCTdJDMC3Cgze8LMXgsfH2Zm/6/zQ5POUlyQy3dOHcOjl53AEYN68qOH3uDM3z6rXe8iIhkqmabYjcD3gL0A7j4HtdSzwrA+Zfzpwolc/6mj2LargfNveIGv3DGLxeu3pzo0ERFph2R6pZe4+0tmB5yj3NBJ8UgXMzM+MK4vJ47qw/XTF3L90wt47PU1fPKYwXz95JH0Ki1IdYgiItKGZFrm681sOPuvZ/4RYFWnRiVdrig/l0tPGcm0/57MxyYO5M8vLOHEK5/iuqfe0alsIiJpLplifglwPcEY7SuAywh6uEsWqi4v4mcfOpRHLzueo4f15lePzmPKr6dx78xlNOoSqyIiaSmZ65kvdPdTgD7AGHc/zt2XdH5okkojqsu56TN13PXFY6guL+S/75vD+34znQdmr1BRFxFJM8mMzd4buBw4jmBX+7PAj9097bo+60IrncPdmbmmkQfe2cPyeqd/qXH2iAIm9s0lx0z5jJjyGS3lM1rKZ/S66kIr/waeBpqr4yeAyWFrPS3pQiudo6nJ+ddrq7n68fm8vbaeUTVlXHbKKIrWv8VJU6akOrysofUzWspntJTP6EVxoZVkerP3cvf/jXn8EzM7pyNvKpkpJ8c447B+nDq+L/+Yu4prHp/PV+54mQFlxtYeKzjzsH4aeEZEJAWS+eZ9yszON7Oc8PYx4B+dHZikr9ycYLz3x/7rRK4+bwJNwGV3z2byr6dx+/OL2blHvd9FRLpSMsX8S8CdwO7wdhfwDTPbZma6pmY3lptjnHNELT85tpibPl1HdXkhP3zgdY775ZP89om32bJjb6pDFBHpFpIZm728KwKRzJVjxkljazj5kGpmLN7EH6a9w1X/ns8fpy/gvImD+Ox7hzCod0mqwxQRyVrJjM3++bjHuWZ2eeeFJJnKzJg0tBe3XjiJf116PO8bW8Ptzy/mxF8/xUW3z+T5BRtoq8OliIi0XzK72U82s3+aWT8zOxR4AVBrXVp1SL8Krj7/CJ79zklcMnkEs5Zs4oIbX+C0a57hnhnLNKqciEiEkhk05uPAn4C5BB3fLnP3b3V2YJId+lYW8a0PjOY/3z2JX557KADfvn8O7/3Fk/zykbdYumFHiiMUEcl8bR4zN7ORwKXA/cAhwKfM7BV317ewJK0oP5fzJg7iY3UDeX7hBm57bjHXT1/AH6Yt4PiRVXzi6MGcfEg1+Tq1TUSk3ZI5z/wh4BJ3f8KCS6d9A5gBjOvUyCQrmRnvHV7Fe4dXsWrLTu6esYy7Zyzj4qmzqC4v5PyJAzlv0iBqexSnOlQRkYyRTDGf5O5bATzovXSVmT3YuWFJd9CvspjLThnFV6eMYNq8ddzx4hJ++9Q7/O6pd5gyuprzJw1i8ug+aq2LiLQhYTE3s2+7+5XuvtXMPuru98ZMvhD4fueHJ91BXm4Op4yt4ZSxNSzbuCNorc9cxhO3z6SqrIBzJtRy7lEDOKRfRapDFRFJS601ec6Puf+9uGmndkIsIgzsVbKvw9xNn66jbnAv/vT8Yk675hnOuPYZbn1uERu370l1mCIiaaW13eyW4H5Lj0UilR/TWt+4fQ8Pzl7B/S+v4EcPvcHP/vkmU0ZX85GjBjB5dDUFedoNLyLdW2vF3BPcb+mxSKfpVVrAZ48dymePHcpbq7dy/6zl/O2VlTz2xhoqi/M5bXxfPnh4f44e1pvcHP3OFJHup7Vifng49roBxTHjsBtQ1OmRibRgTN8K/ueMsXzn1DE88/Z6Hnx1JQ+9upK7ZiyjuryQMw/rzwcn9OfwAZUEJ1+IiGS/Nq9nnknM7CzgrNra2oumTp2acL76+nrKysq6LrAsl+p87m50Xl3byAurGpizrpEGh+oS4+i+eUzql8eAMsuowp7qfGYb5TNaymf0pkyZ0uHrmWdVMW9WV1fnM2fOTDh92rRpTJ48uesCynLplM8tO/fy6OureXD2Sv6zYD1NDkOrSjl1fF9OG9+XQ2vTv8WeTvnMBspntJTP6JlZh4t5MueZi2SMyuJ8PlY3kI/VDWTdtt089sZqHnltNTc8vZA/TFtAbY9iTh3fl1PH9+WoQT3J0TF2EckCKuaStfqUF/KJowfziaMHs3nHHv79xhoeeW01f35+CTc/u4g+5YV8YFwNp43vx9FDe5GnwWlEJEOpmEu30KOkgI/WDeSjdQPZtmsvT761lkdfX839s1Yw9YWlVBbnM3l0H04+pIYTR/Whsjg/1SGLiCRNxVy6nfKifM6eUMvZE2rZuaeR6fPX8fiba3jqrbU8MHsleTnGxCG9OPmQak45pIYhVaWpDllEpFUq5tKtFRfk7juG3tjkzF62mSfeXMMTb67lJ/94k5/8402G9ynl5ENqOHlMNUcN7qnd8SKSdlTMRUK5OcZRg3ty1OCefPvUMSzbuIMn3lzD42+u5dbnFnHD0wupKMrjuJFVnDiqDyeM6kO/Sl3dTURST8VcJIGBvUr2jTy3bddennl7PdPmrWX6/HX8c+5qAEbVlHHCyD6cOLoPE4f0oig/N8VRi0h3pGIukoTyonxOP7Qfpx/aD3dn/pp6np6/junz13H780u46dlFFOXncPTQ3vta7cP7lKb9Oe0ikh1UzEXaycwY3bec0X3LueiEYezc08gLizYwfd46nn57HT9++A0AansUc/zIKt4zvDfvGd6b6nKNgiwinUPFXKSDigtymTK6mimjqwFYtnEHT7+9jqfnr+Ofc1dx14xlQLBL/r3Dg+J+zNDeVJbo9DcRiYaKuUjEBvYq2TdYTWOT88bKrTy3YD3/WbCBu2cs47b/LCbHYHxtJe8Z3pv3Dq9i4pCelBRocxSRg6NvD5FOlJtjHDqgkkMHVHLxicPZ09DE7GWb+U9Y3G95dhHXT19Ifq4xYWAPjh7am6JtDdTtbqCsUJuniCRH3xYiXaggL4dJQ3sxaWgvLjsFduxpYObiTTy3YD3PL9jAH6YvoLHJ+c3LjzGufwUThwTzThzSi16lBakOX0TSlIq5SAqVFORxQtj7HaB+dwO3PTSN3RUDeXHRRv78QjCOPMDI6rJ9PwQmDulF/x46x11EAirmImmkrDCP8VV5TJ48GoDdDY3MWb6FlxZt5KVFG3lg9krueHEpAAN6Fu8r7EcO6snI6jJdBU6km1IxF0ljhXm5TBwSFOxLpkBjk/Pmqq28tGgjMxZvZPq8dfz15RUAlBfmMWFQD44c1JMjB/dkwsAeumCMSDehYi6SQXJzjPG1lYyvreRzxw3F3Vm8YQcvL9nEy0s3MWvJJn775Ns0OZjBiD5lHDkoGKL2yME9GFal1rtINlIxF8lgZsbQqlKGVpVy7lEDgOC4+6vLNvPykk3MWrqJR15fzd0zg3PdK4ryOGJQT44c1JPDB1Zy+IAe9FTHOpGMp2IukmXKCvM4dkQVx46oAqCpyVm4fjsvL920rwV/9RPrcA/mH9SrhMMGBIX9sAFBq79Up8WJZBRtsSJZLifHGFFdxojqMj5WNxCArbv28tryLby6fAtzlm/mlaWbeXjOqmB+g5HV5Rw2oJLDBvZgwoAejO5bTkGeLv0qkq5UzEW6oYqifN47oor3hq13gHXbdjN3xWZmLwsK/BNvreXeWcsBKMjN4ZD+FRw+oJLDBvTg0NpKhvcp1bXdRdKEirmIANCnvJCTxtRw0pgaANyd5Zt2Midsvc9etpn7Zy3n9ueXAFCYl8OYfhWM71/BuP6VjK+tYFRNuS4DK5ICKuYi0iIzY2CvEgb2KuGMw/oBwalxC9fV8/rKrby2Yguvr9zKg6/uP/c9L9yl31zcx/WvZGz/Cg1NK9LJ0n4LM7NzgDOAauA6d38sxSGJdFu5OcbImnJG1pRzzhG1wP4WfHNxf23lFqbPX8f9Lwe76M1gSO9SxvXfX9wP6VtOn/JCXe9dJCKdWszN7BbgTGCtu4+Pef5U4BogF7jJ3X+RaBnu/nfg72bWE/g1oGIukkZiW/CnHdpv3/Nrt+46oAU/e9n+TnYAvUoLOKRfOWP6VjCmbzmH9KtgRHWZdtOLHITObpnfBvwOuL35CTPLBa4D3gcsB2aY2YMEhf3nca//nLuvDe//v/B1IpIBqiuKqK4oYsqY6n3Pbd6xh7dWb+OtVVt5a/U23ly9jTteXMKuvU1A0PIfWlXKIf2aC3xQ7PtVFqkVL9KKTi3m7v60mQ2Je3oS8I67LwQws7uAs9395wSt+ANYsAX/AviXu7/cmfGKSOfqUVLAMcN6c8yw3vuea2xylmzYvq/Iv7l6G68s3cRDr67cN09FUR5j+lUwtl8Fo/uWM6qmnJE1ZVQUabhaEQDz5pEjOusNgmL+cPNudjP7CHCqu38hfPwp4Gh3/2qC138d+AwwA5jt7n9MMN8XgS8C1NTUHHXXXXcljKm+vp6ysrKD/UgSR/mMlvIZ2LHXWVHfxLJt+2/LtzWxq3H/PL2KjP5lOdSWGbVlOdSW5dC/LIfivP2teOUzWspn9KZMmTLL3es6soxUdIBraV9Zwl8U7n4tcG1bC3X3G4AbAOrq6nzy5MkJ5502bRqtTZf2UT6jpXwm1tQUdLabv2Yb89du4+019cxfs41py+vZ3dCwb77aHsWMqC5jVE0ZTZv3ctbhExhZXaaR7SKg9TM9pWLNXg4MjHk8AFiZYF4RkX1ycoxBvUsY1LuEU8bW7Hu+sclZtnEH89ds4+21QYGfv6ae5xduYE9DEze/9hwQFPlRNWXhbvpyRlaXMaxPKeXaXS8ZLhXFfAYw0syGAiuA84GPpyAOEckSuTnGkKpShlSV8v5x+59vbHLu/ddT9BwylrfDAj9/zTaee2cDexqb9s1XU1HI8D5l4a2U4dXBfXW8k0zR2aem/QWYDFSZ2XLgcne/2cy+CjxK0IP9Fnd/vTPjEJHuKTfH6Fuaw+RxffnAuL77nm9obGLJxh28s7aeBevqWbB2OwvW1fP32SvYtmv/7vqSglyG9SlleJ8yhlWVMbw6uD+0qlSn0Ela6fQOcF3JzM4Czqqtrb1o6tSpCedTB45oKZ/RUj6j1Z58ujtb9jir6p1V25vCm7N6exPrd+7/rjSgqtjoV5pDv1KjX1kO/Upz6FuaQ0UBWd2a1/oZvSg6wGVVMW9WV1fnM2fOTDhdHTiipXxGS/mMVlT53LmnkUXrgxZ8cNvOgrX1LFxfv+88eQguQTukqoQhvYPrzA/pHez+H1pVSs+S/Iwv9Fo/o2dmGdmbXUQk4xQX5DK2fwVj+1cc8HxTk7Nyy04WrNvOwnX1LF6/nUUbdjBn+Rb+OXcVTTHtpYqivKDAh0W++f7Q3qVUlqgTnhw8FXMRkQ7IyTEG9CxhQM8SThzV54BpexqaWLZpR1Dg129n8YbtLF6/g5mLN/HgqyuJ3THasyR/X2EfElPkh1SVqLe9tEnFXESkkxTk5ezrJR9v195Glm7cERT5sNAvWr+d/yzYwF9fWXHAvFVlBQzsVcLgXiUMCsfBH9y7lEG9SqguLyQnJ7N33UvHZdUxc3WASw3lM1rKZ7QyMZ+7G521O4KOd2t2NLF2h7Mu/Ltxlx8wylZ+DvQpNvqU5NCn2KguyaFPSfi32CjIjbbQZ2I+0506wCWgDnBdS/mMlvIZrWzL556GJlZs3snSjTuC24bt4f2dLN2wne17Gg+Yv7q8kMG9g9b8oF4lDO4d/B3Uq5SqsoJ2d8jLtnymA3WAExHpZgrychga9o6P5+5s3L4nptDvYEl4//kFG/jrywfuvi/Ozw132xeHx/2Lw1twv7I483vfdxcq5iIiWcLM6F1WSO+yQo4Y1PNd03ftbWT5pp0s27iDJRu2B635jdtZvmknLyzcSP3uhgPmLyvMo7ZH8QFFfvPqBnov38KAnsX0yIJT7bKFirmISDdRlJ/LiOoyRlS/+5i3u7N1ZwPLNu1g+aadLN/3dycrNu/kpUUb2RYW++tmPwtAaUHuu1r0tTH3s+G8+kyRVcfM1QEuNZTPaCmf0VI+o7N9r7Nsw3a2WxHrdzrrdzaFf4P7Ow9s2FOYG4yUV1WcQ+9io1eR0bto//0ehUaeeuKrA1wi6gDXtZTPaCmf0VI+o9VaPrfs3MuKuFb98k07WLZpJ6u27GTzjr0HzJ9jUFNRRL/KIvr3KKa2R/G++8237tC6Vwc4ERFJG5XF+VQW579rlLxm23c3sGrLTlZu3sXKzTuD25bg/msrtvDYG2vY09B0wGuK8nOCwl5ZTP8eRTH39z/WRW9UzEVEpIuUFuYxorqcEdXlLU53dzZs3xMW+v0Ff9WWXazYvJPp89exdttu4nco9yotoH+PIvpVBq37vpVBa7+51V9TUZT1BV/FXERE0oKZUVVWSFVZIYcNaHmePQ1NrNkaFvq4Vv7SDTt4YcGGfR31YvUsyadvZTF9KwrDv2Ghjyn8FUV5GbtLX8VcREQyRkFeDgPDIW0Tqd/dwOotu4Lb1l2s3rIz/Bs8nrtiC+vr97zrdSUFufStKKJvZdH+v3H3q0rTc/hcFXMREckqZYV5CU/Ba7a7oZG1W3fvL/L7Cn/w98VFG1mzdRcNTQfu08/LMWriC35FEdUVhfStCFr4NRVFFBd07W79rOrNrlPTUkP5jJbyGS3lM1rdKZ9N7mzd42zaFdw27nI27w7+btrVFDy324kbQReA4jzoWWj0KDJ6FObE3LcD7uflWCSnpmVVy9zdHwIeqquru6i1U1F0qkq0lM9oKZ/RUj6jpXweyN3ZuquBtVt3sWbrbtZs3cWabbtY23x/6y6WbN3NjDW72Nv47sZzr9KCSOLIqmIuIiLSlcxs3yl5I2ta7qUP0NTkbNqxJyj423btK/6rt+7ilQjiUDEXERHpZDk5+8fNH8uB5+H/PIrlR7AMERERSSEVcxERkQynYi4iIpLhVMxFREQynM4zlw5TPqOlfEZL+YyW8hk9XQI1AV0CtWspn9FSPqOlfEZL+YxeFJdA1W52ERGRDKdiLiIikuFUzEVERDJcVh4zN7N1wJJWZqkC1ndRON2B8hkt5TNayme0lM/ojXb3xGPBJiErh3N19z6tTTezmR3tbCD7KZ/RUj6jpXxGS/mMnpkl7rGdJO1mFxERyXAq5iIiIhmuuxbzG1IdQJZRPqOlfEZL+YyW8hm9Duc0KzvAiYiIdCfdtWUuIiKSNbpdMTezU81snpm9Y2bfTXU8mcjMFpvZXDOb3dwL08x6mdm/zezt8G/PVMeZrszsFjNba2avxTzXYv4scG24vs4xsyNTF3l6SpDPK8xsRbiOzjaz02OmfS/M5zwz+0Bqok5fZjbQzJ4yszfN7HUzuzR8XuvoQWgln5Guo92qmJtZLnAdcBowFrjAzMamNqqMNcXdJ8ScovJd4Al3Hwk8ET6Wlt0GnBr3XKL8nQaMDG9fBP7QRTFmktt4dz4BfhOuoxPc/Z8A4fZ+PjAufM3vw+8F2a8B+Ka7HwIcA1wS5k3r6MFJlE+IcB3tVsUcmAS84+4L3X0PcBdwdopjyhZnA38K7/8JOCeFsaQ1d38a2Bj3dKL8nQ3c7oEXgB5m1q9rIs0MCfKZyNnAXe6+290XAe8QfC9IyN1XufvL4f1twJtALVpHD0or+UzkoNbR7lbMa4FlMY+X03pSpWUOPGZms8zsi+FzNe6+CoKVF6hOWXSZKVH+tM4evK+Gu31viTnso3y2g5kNAY4AXkTraIfF5RMiXEe7WzG3Fp5Td/72O9bdjyTYvXaJmZ2Q6oCymNbZg/MHYDgwAVgFXBU+r3wmyczKgPuBy9x9a2uztvCcchqnhXxGuo52t2K+HBgY83gAsDJFsWQsd18Z/l0L/I1gF9Ca5l1r4d+1qYswIyXKn9bZg+Dua9y90d2bgBvZv5tS+UyCmeUTFJ473P2v4dNaRw9SS/mMeh3tbsV8BjDSzIaaWQFBJ4MHUxxTRjGzUjMrb74PvB94jSCPnwln+wzwQGoizFiJ8vcg8Omwx/AxwJbmXZ2SWNwx2w8RrKMQ5PN8Mys0s6EEnbZe6ur40pmZGXAz8Ka7/1/MJK2jByFRPqNeR7PyQiuJuHuDmX0VeBTIBW5x99dTHFamqQH+Fqyf5AF3uvsjZjYDuMfMPg8sBT6awhjTmpn9BZgMVJnZcuBy4Be0nL9/AqcTdILZAVzY5QGnuQT5nGxmEwh2Ty4GvgTg7q+b2T3AGwS9jC9x98ZUxJ3GjgU+Bcw1s9nhc99H6+jBSpTPC6JcRzUCnIiISIbrbrvZRUREso6KuYiISIZTMRcREclwKuYiIiIZTsVcREQkw6mYi3QiM3Mzuyrm8bfM7IqIln2bmX0kimW18T4fDa/49FTMc4fGXO1po5ktCu8/bmaTzezhzo5LRPZTMRfpXLuBD5tZVaoDidXOK4V9HviKu09pfsLd5zZf7YlgkIv/Dh+fEnWsItI2FXORztUA3AD8V/yE+Ja1mdWHfyeb2XQzu8fM5pvZL8zsE2b2kgXXkR8es5hTzOyZcL4zw9fnmtmvzGxGeBGHL8Us9ykzuxOY20I8F4TLf83Mfhk+90PgOOCPZvardnzuMjO7z8zeMrM7wlGwMLMfhnG9ZmY3xDw/zcx+GX7G+WZ2fPh8kZndGsb1iplNCZ8fF847O/yMI8PnPxnz/PVhLnLDXL8WLudd/wuRTNetRoATSZHrgDlmdmU7XnM4cAjBpT0XAje5+yQzuxT4GnBZON8Q4ESCCzY8ZWYjgE8TDKk50cwKgefM7LFw/knA+PDSivuYWX/gl8BRwCaCq+Kd4+4/NrOTgG+5+8x2xH8EwfWYVwLPEYyC9SzwO3f/cfiefwbOBB4KX5MXfsbTCUZxOwW4BMDdDzWzMWFco4CLgWvc/Q4LhmbONbNDgPMILgS018x+D3wCeB2odffx4fv2aMfnEMkIapmLdLLwCkm3A19vx8tmhNdB3g0sAJqL8VyCAt7sHndvcve3CYr+GILx8j8dDh35ItCbYHxngJfiC3loIjDN3de5ewNwB9CRq+G95O7Lw4tIzI6JeYqZvWhmc4GTCAp+s+YLesyKme/MQm4AAAH0SURBVP844M8A7v4WsAQYBTwPfN/MvgMMdvedwMkEP0ZmhJ/9ZGAYQV6GmdlvzexUoLUrgIlkJLXMRbrG1cDLwK0xzzUQ/qAOdzcXxEzbHXO/KeZxEwdut/HjMTvBJRS/5u6Pxk4ws8nA9gTxtXTZxY6Ijb8RyDOzIuD3QJ27Lws7Aha18JpG9n/GFuNy9zvN7EXgDOBRM/tCOO+f3P178fOb2eHABwha+h8DPnewH0wkHallLtIF3H0jcA9BZ7JmiwlakgBnA/kHseiPmllOeBx9GDCP4EJCX7bgsouY2SgLrnDXmheBE82sKuwcdwEw/SDiaU1z4V5vwbWdk+mJ/zTBrnLC3euDgHlmNgxY6O7XEnTAOwx4AviImVWH8/cys8Fh58Mcd78f+AFwZJQfSiQdqGUu0nWuAr4a8/hG4AEze4mgECVqNbdmHkHRrQEudvddZnYTwW7ql8MW/zrgnNYW4u6rzOx7wFMELdx/unukl7F1981mdiPBoYLFBJckbsvvCTrfzSXYk/FZd99tZucBnzSzvcBq4MfuvtHM/h/BcfUcYC9BS3wncGv4HMC7Wu4imU5XTRMREclw2s0uIiKS4VTMRUREMpyKuYiISIZTMRcREclwKuYiIiIZTsVcREQkw6mYi4iIZDgVcxERkQz3/wEcIOcOBWvw2gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 504x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(7,4))\n",
    "plt.semilogy( exprate, label='true rate' )\n",
    "# plt.semilogy( 1.42 * p ** (np.log2(x+1)), label='approximate rate')\n",
    "plt.grid(which='both')\n",
    "plt.xlabel('Number of Thanoses')\n",
    "plt.ylabel('Expected survival rate')\n",
    "plt.title(\"Expected global survival rate if N Thanoses sequentially snap their fingers\")\n",
    "plt.xlim((-5,Nmax))\n",
    "# plt.legend();\n",
    "plt.tight_layout(pad=0)\n",
    "plt.savefig('thanos2.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 348,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.02219897]\n",
      "[1.66492281e+08]\n"
     ]
    }
   ],
   "source": [
    "# print global survival rate and expected number of survivors\n",
    "print( exprate[63] )\n",
    "print( exprate[63] * 7.5e9 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}