power_and_sample_size.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as st\n",
    "import scipy.optimize as op\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "figsize = [3, 2.25]\n",
    "alpha = 0.05\n",
    "n_iter = 1000000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculation functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_lift_dist(reach, true_lift, conv_ctrl, ctrl_scale_factor, n_iter=1000000):\n",
    "    \n",
    "    lam_C_c = conv_ctrl\n",
    "    lam_C_t = ctrl_scale_factor * lam_C_c * (1 + true_lift * reach)\n",
    "\n",
    "    C_t = np.random.poisson(lam_C_t, n_iter)\n",
    "    C_c = np.random.poisson(lam_C_c, n_iter)\n",
    "\n",
    "    R_c = C_c * reach\n",
    "    \n",
    "    L = (C_t - ctrl_scale_factor * C_c) / (ctrl_scale_factor * R_c)\n",
    "    \n",
    "    return L\n",
    "\n",
    "def get_crit_val(reach, conv_ctrl, ctrl_scale_factor, alpha=0.05, n_iter=1000000):\n",
    "    \n",
    "    # Get the distribution of L under H0.  True lift = 0\n",
    "    L = get_lift_dist(reach, 0, conv_ctrl, ctrl_scale_factor, n_iter)\n",
    "    \n",
    "    crit_val = np.percentile(L, 100*(1 - alpha))\n",
    "    \n",
    "    return crit_val\n",
    "\n",
    "def get_lift_study_power(reach, min_lift, conv_ctrl, ctrl_scale_factor, alpha=0.05, n_iter=1000000):\n",
    "    \n",
    "    crit_val = get_crit_val(reach, conv_ctrl, ctrl_scale_factor, alpha, n_iter)\n",
    "    \n",
    "    # Get the distribution of L under H1\n",
    "    L = get_lift_dist(reach, min_lift, conv_ctrl, ctrl_scale_factor, n_iter)\n",
    "    \n",
    "    power = (100 - st.percentileofscore(L, crit_val)) / 100\n",
    "    \n",
    "    return power\n",
    "\n",
    "def lift_study_power_helper(conv_ctrl, *args):\n",
    "    reach, min_lift, alpha, min_power, ctrl_scale_factor, n_iter = args\n",
    "    power = get_lift_study_power(reach, min_lift, conv_ctrl, ctrl_scale_factor, alpha, n_iter)\n",
    "    \n",
    "    return power - min_power\n",
    "\n",
    "def get_min_conv_ctrl(reach, min_lift, ctrl_scale_factor, alpha=0.05, power=0.8, n_iter=1000000):\n",
    "    \n",
    "    a = 100\n",
    "    b = 50e6\n",
    "    xtol = 10\n",
    "    rtol = 0.01\n",
    "    \n",
    "    conv_ctrl = op.bisect(lift_study_power_helper, a, b, \n",
    "                          args=(reach, min_lift, alpha, power, ctrl_scale_factor, n_iter), \n",
    "                          xtol=xtol, rtol=rtol)\n",
    "    \n",
    "    return conv_ctrl\n",
    "\n",
    "def get_lift_diff_power(reach, lift_a, min_lift_uplift, conv_ctrl_a, ctrl_scale_factor, alpha=0.05, n_iter=1000000):\n",
    "    \n",
    "    # Calculate the critical value under H0\n",
    "    L_a = get_lift_dist(reach, lift_a, conv_ctrl_a, ctrl_scale_factor, n_iter)\n",
    "    L_b = get_lift_dist(reach, lift_a, conv_ctrl_a, ctrl_scale_factor, n_iter)\n",
    "    \n",
    "    L_ab = L_b - L_a\n",
    "    \n",
    "    crit_val = np.percentile(L_ab, 100*(1-alpha))\n",
    "    \n",
    "    # Calculate the test power under H1\n",
    "    # NOTE: WE ARE ASSUMING THE CONTROLS IN A AND B HAVE THE SAME CONVERSION RATE WHICH IS NOT THE CASE\n",
    "    # - HOWEVER IT'S THE ONLY ESTIMATE WE HAVE\n",
    "    L_b = get_lift_dist(reach, lift_a + min_lift_uplift, conv_ctrl_a, ctrl_scale_factor, n_iter)\n",
    "    \n",
    "    L_ab = L_b - L_a\n",
    "    \n",
    "    power = (100 - st.percentileofscore(L_ab, crit_val)) / 100\n",
    "    \n",
    "    return power\n",
    "\n",
    "def lift_diff_power_helper(conv_ctrl_a, *args):\n",
    "    reach, lift_a, min_lift_uplift, ctrl_scale_factor, min_power, alpha, n_iter = args\n",
    "    power = get_lift_diff_power(reach, lift_a, min_lift_uplift, conv_ctrl_a, ctrl_scale_factor, alpha, n_iter)\n",
    "    \n",
    "    return power - min_power\n",
    "\n",
    "def get_min_conv_ctrl_a(reach, lift_a, min_lift_uplift, ctrl_scale_factor, alpha=0.05, power=0.8, n_iter=1000000):\n",
    "    \n",
    "    a = 100\n",
    "    b = 50e8\n",
    "    xtol = 10\n",
    "    rtol = 0.01\n",
    "    \n",
    "    conv_ctrl_a = op.bisect(lift_diff_power_helper, a, b, \n",
    "                          args=(reach, lift_a, min_lift_uplift, ctrl_scale_factor, power, alpha, n_iter), \n",
    "                          xtol=xtol, rtol=rtol)\n",
    "    \n",
    "    return conv_ctrl_a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "########################\n",
    "# Single cell lift study\n",
    "\n",
    "def plot_power_by_sample_size(ax, reach, min_lifts, convs_ctrl, ctrl_scale_factor, alpha=0.05, n_iter=1000000):\n",
    "    \n",
    "    for min_lift in min_lifts:\n",
    "        powers = []\n",
    "        for conv in convs_ctrl:\n",
    "            powers.append(100 * get_lift_study_power(reach, min_lift, conv, ctrl_scale_factor, alpha, n_iter))\n",
    "\n",
    "        ax.plot(convs_ctrl, powers, label='{:.0f}%'.format(100 * min_lift))\n",
    "\n",
    "    ax.set_xlim([0,max(convs_ctrl)])\n",
    "    ax.set_ylim([0, 100])\n",
    "    ax.get_xaxis().set_major_formatter(matplotlib.ticker.FuncFormatter(lambda x, p: thousands_tick_formatter(x, p)))\n",
    "    ax.set_xlabel('Conversions in control')\n",
    "    ax.set_ylabel('Power (%)')\n",
    "    ax.legend(loc='lower right', title='$L_m$')\n",
    "    \n",
    "def plot_power_by_reach(ax, reaches, min_lifts, conv_ctrl, ctrl_scale_factor, alpha=0.05, n_iter=1000000):\n",
    "    \"\"\"Vary the reach holding the total audience size constant\"\"\"\n",
    "    \n",
    "    for min_lift in min_lifts:\n",
    "        powers = []\n",
    "        for reach in reaches:\n",
    "            powers.append(100 * get_lift_study_power(reach, min_lift, conv_ctrl, ctrl_scale_factor, \n",
    "                                                     alpha, n_iter))\n",
    "\n",
    "        ax.plot(100 * reaches, powers, label='{:.0f}%'.format(100 * min_lift))\n",
    "\n",
    "    ax.set_xlim([0,100])\n",
    "    ax.set_ylim([0, 100])\n",
    "    ax.set_xlabel('Reach (%)')\n",
    "    ax.set_ylabel('Power (%)')\n",
    "    #ax.legend(loc='upper right', title='$L_m$')\n",
    "    \n",
    "def plot_power_by_ctrl_scale_factor(ax, reach, min_lifts, conv_ctrl_at_50, ctrl_splits, \n",
    "                                    alpha=0.05, n_iter=1000000):\n",
    "    \"\"\"Vary the fraction of audience in the control group\"\"\"\n",
    "    \n",
    "    for min_lift in min_lifts:\n",
    "        powers = []\n",
    "        for split in ctrl_splits:\n",
    "            s = (1 - split) / split\n",
    "            conv = split * conv_ctrl_at_50 / 0.5\n",
    "            powers.append(100 * get_lift_study_power(reach, min_lift, conv, s, \n",
    "                                                     alpha, n_iter))\n",
    "\n",
    "        ax.plot(ctrl_splits, powers, label='{:.0f}%'.format(100 * min_lift))\n",
    "\n",
    "    ax.set_xlim([0,1])\n",
    "    ax.set_ylim([0, 100])\n",
    "    ax.set_xlabel('Fraction of audience in control group')\n",
    "    ax.set_ylabel('Power (%)')\n",
    "    #ax.legend(loc='upper right', title='$L_m$')\n",
    "\n",
    "#######################\n",
    "# Multi-cell lift study\n",
    "\n",
    "def plot_power_by_sample_size_mc(ax, reach, lift_a, min_lift_uplifts, convs_ctrl_a, ctrl_scale_factor, \n",
    "                                 alpha=0.05, n_iter=1000000):\n",
    "    \n",
    "    for min_lift_uplift in min_lift_uplifts:\n",
    "        powers = []\n",
    "        for conv in convs_ctrl_a:\n",
    "            powers.append(100 * get_lift_diff_power(reach, lift_a, min_lift_uplift, conv, ctrl_scale_factor, \n",
    "                                                    alpha, n_iter))\n",
    "            \n",
    "\n",
    "        ax.plot(convs_ctrl_a, powers, label='{:.0f}%'.format(100 * min_lift_uplift))\n",
    "\n",
    "    ax.set_xlim([0,max(convs_ctrl_a)])\n",
    "    ax.set_ylim([0, 100])\n",
    "    ax.get_xaxis().set_major_formatter(matplotlib.ticker.FuncFormatter(lambda x, p: thousands_tick_formatter(x, p)))\n",
    "    ax.set_xlabel('Conversions in cell $A$ control')\n",
    "    ax.set_ylabel('Power (%)')\n",
    "    ax.legend(loc='lower right', title='$D_m$')\n",
    "    \n",
    "def thousands_tick_formatter(x, p):\n",
    "    if x >= 1000:\n",
    "        x = str(int(x / 1000)) + 'k'\n",
    "    else:\n",
    "        x = int(x)\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Single cell lift study\n",
    "\n",
    "### Power visualisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOgAAACrCAYAAABsZdC6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFCtJREFUeJztnXu8VWWZx78/EeegYCWgGZcAR2moMSJqtGbmg6Z9EoJs\ntNtHTRonM6K8ZIWXyjELa3KwpotiCt6q0RGL0BwPJjlNZB0EVLxAMprHTPBYGQji5Zk/3nVwc9x7\nn/fsvddea+39fD+f9VnXd63f3uc8+33X+z7v88jMcBwnn+yWtQDHcSrjBuo4OcYN1HFyjBuo4+QY\nN1DHyTFuoI6TY9xAHSfHuIE6To5xA3WcHLN71gLqYcSIETZu3LisZbQvDz4Y1hMnZqsjh6xatepJ\nMxtZ730KbaDjxo2jq6sraxnty7RpYb1iRZYqcomkRxpxH2/iOk6OSd1AJQ2StFrSsmR/H0mdkjYk\n61clx2dL+lbaepwGMnlyWJzUaEYNeipwf8n+POA2MzsQuC3Zd4rIxReHxUmNVN9BJY0GZgBfBs5I\nDr8HmJZsXwmsAD7Xp9wM4Fxgppk9OZBnPvfcc3R3d7N9+/bahWdIR0cHo0ePZvDgwVlLcXJA2p1E\nFwOfBYaVHNvPzB5Ptv8A7FdaQNJ7CcY83cz+ONAHdnd3M2zYMMaNG4ekGmVng5nR09NDd3c348eP\nz1pO/xx/fFhfc03znnn7fDjsrOY9L2NSa+JKejewycxWVbrGwmzx0hnjhxNq0xmVjFPSyZK6JHVt\n3rz5Zee3b9/O8OHDC2ecAJIYPnx4cWr/7u6wNIvb5++6bgPSfAd9OzBL0sPAD4HDJV0DPCFpf4Bk\nvamkzEOE2vagSjc1s4VmNtXMpo4cWX6YqYjG2UuRtadKGxllKakZqJmdZWajzWwc8EHgZ2Z2PLAU\nODG57ETgxyXFHgGOAa6S9Pq0tDlOUchiHPRC4EhJG4Ajkv2dmNkDwHHA9ZIOyECf4+SGpngSmdkK\nQm8tZtYDvKPMNYuBxcn2amBSM7Q5dXDooc15TrnmbZt0FhXa1S+GBZ3rG3q/04+s+Hq8C0OHDmXL\nli079xcvXkxXVxff+lb+fDEWdK6P/ly7ML893wubibv6tTmN/gFrKm3QceQG6uxkwMZ6zDFhcVKj\n5Zu4WbFt2zYml/ipPvXUU8yaNStDRdXpNc4BNXd7elJU5ECEgUoannTsOANgyJAhrFmzZud+7zvo\n1q1bmTNnDnvssQfTpk3juOOOy1Clk3dimri/knS9pOnyUfS6WbJkCcceeyyXXXYZS5cuzVRLod8/\n24QYAz0IWAicAGyQ9BVJNXT5ORB8hceMGQPAoEGDMlZTnlwZbht0BFWj3yZu4i/bCXRKOgy4Bpgj\naS0wz8xWpqyxLmoaPkiR0aNH093dzeTJk3nxxRezllMf73jZcLbTYKLeQYHjCTXoE8AnCe56k4Hr\ngQJMu2g+pWOgALNnz2b27Nls3bqVuXPnctNNNzFz5syM1DWolvz85+u/h1OVmF7clcDVwNFmVjp1\noUvSJenIal322msvFi1alLUMpyDEvIOea2ZfKjVOSe8DMLOvpqbMyT9HHRUWJzViDLRcSJLWd4Js\nc6KawNu2hSUtYjqIWrwTqWITV9JRwHRglKRvlpzaG3g+bWGO41R/B/090AXMAkqjIvwFOD1NUY7j\nBCoaqJmtBdZKutbMvMZ0nAyo1sS9zszeD6yWVBo3SITh0YNTV+fkm3e/O2sFLU+1Ju6pydr/Ck55\nzjwzawWBFp68XbEXtyQ05pPAo2b2CPBXwBsJ76eOkx4t3jsbS8wwyx1Ah6RRwK0Ej6LFaYpqBYYO\nHbrL/uLFi5k7d25GanalYb6206a9lEDJSYUYTyKZ2TOSTgK+Y2Zfk7Sm31J5odG/xC3alHLySUwN\nKkmHEiLt3ZQc63cahqQOSb+WtFbSOkn/mhz35EmOE0lMDXoqwXPoRjNbJ2kCcHtEuWeBw81si6TB\nwC8k/RT4J0LypAslzSN4Kn2u2o2KSNEiKpSj5mBiTsPotwY1szvMbFav362ZbTSzT0WUMzPrndIx\nOFmMkDzpyuT4lcDRfctKmiFppaQRkZ8jd/RGVOhdzj//fAA2btzISSedxLHHHpuxwhajRTuVYqab\nHQScCYwrvd7MDo8oO4jghfTXwLfN7E5JqSZPyjsTJkzg8ssvbw0Dff/707lvixpbLcQ0ca8HLgG+\nB7wwkJub2QvAZEmvBG6U9IY+562PE8ThwFTgnWb2dLl7SjoZOBlg7NixA5Hj0OBoCXPmNO5eTlli\nDPR5M/tuPQ8xsz9Juh14F0nyJDN7vELypAmEMCtdFe61kBCChalTp1q5a5wm8cwzYb3nntnqaGFi\nDPQnkuYANxI6fgAws6eqFZI0EnguMc4hwJHAV3kpedKFlE+e9BlgiaT3mdm6gXyYsmQ0LFIpokJP\nTw/nnHMOq1evZv78+Zx1Vr6Hbap2FE2fHtYrVjRNT7sRY6C9mcg+U3LMCDVdNfYHrkzeQ3cDrjOz\nZZJWAtcl46qPALu8yJjZA5J6kyfNNLOHYj5IURg+fDiXXOKBKJw4YoKG1RRzyMzuBt5U5rgnT3Kc\nSPodZpG0p6RzJS1M9g9Msmc7jpMyMZ5Ei4AdwNuS/ceAC1JT5KRGruLdpkELDs/EvIMeYGYfkPQh\ngMQv1yPMOzB7duPv2YJGVg8xNeiOpBfWAJKs189WL5ItIdZ2Mcmj9oo17+zZ6Rips5MYA/0icAsw\nRtK1wG3AZ1NVVQcdHR309PTk8h+9P8yMnp4eOjo6spYSx5NPhsVJjZhe3E5JdwGHEMKdnGpmuf2r\n9KZW2Lx5c9ZSaqKjo4PRo0dnLSOOXndFHwdNjWoxiab0OdTrPztW0lgzuys9WbUzePBgxo/3bBRO\na1CtBr0oWXcQ/GPXEmrQgwlueIemK81pJC3fg9uiVItJdJiZHUaoOaeY2VQzezPB+eCxZgl0nAHR\nYr3AMZ1EE83snt4dM7sX+Jv0JDltS4sZVyOIGQe9W9L3CHlBIYQ+uTs9SU5h+PjHs1bQ8sTUoB8B\n1hFCn5wK3Jccc9qIsu+wH/hAWJzUiBlm2Q4sSBbHeYlHHw3rMWOy1dHCxDRxHac8J5wQ1j4Omhox\nTVzHcTIiZrrZ3zZDiOM0jBbqDY6pQb+TBKCeI+kVqStyHGcnMXFx/4EwtDIGWCXp+5KOTF2Z0160\nUK3XSKI6icxsg6RzCS5+3wTelMwJPdvMlqQp0Mkxn/501gpanpjA1QcTxj1nAJ3ATDO7S9JrgJWA\nG2i7MnNm1gpanph30P8A7gLeaGaf6J3FYma/B86tVEjSGEm3S7ovSZ50anLckyc1mdQc5R98MCxO\nasQY6Azg+2a2DUDSbpL2BDCzq6uUex74tJlNIswl/YSkSYRkSbeZ2YGEyd/z6vkAToZ87GNhcVIj\nxkCXA0NK9vdMjlXFzB4vqW3/AtwPjKJNkic5GdMinU4xBtpRkqWMZHtAsf4ljSNMU7sTiEmeNI+Q\nPCm3kRvaEZ9T2nxiDHRraXQFSW8GtsU+QNJQ4AbgtL4JkSwEDuqbPOlzwIxKmc0knSypS1JXUcOa\nOE4sMQZ6GiENw/9I+gXwn8DcmJsniXtvAK4tGY55IkmaRIXkScMIyZPKYmYLk8njU0eOHBkjw3EK\nS8xslt9Ieh0wMTn0oJk911+5ZJz0cuB+M/v3klPNTZ7kpMe5FTvxnQYRO5vlLbyUwHeKJMzsqn7K\nvB04AbhH0prk2NkEw2zb5EnNJtX3xiOOSO/eDhDnqHA1cACwhpcS+BpQ1UDN7BeEIGPl8ORJrcCa\n5Hd38uRsdbQwMTXoVGCSFTEStJMup50W1j4fNDViOonuBV6dthCnGKTSZG6RMcs0iKlBRwD3Sfo1\nu2bYnpWaKsdxgDgDPS9tEY7jlCdmmOXnkl4LHGhmyxM/3EHpS3PqxT1/ik9MyJOPAv8FXJocGgX8\nKE1RTkH4ylfCklda4N02pon7CeCtBD/a3snb+6aqyikGb3tb/9c4dRHTi/usme3o3ZG0O7v6zzpt\nxs6m8y9/GRYnNWJq0J9LOhsYksQimgP8JF1ZTiE4++ywrmcctAWaoWkSU4POAzYD9wAfA26mSiQF\nx3EaR0wv7ovAZcniOE4TienF/T9JG/suzRDn1I4PsSQUvAkd64vbSwfwPmCfdOQ4RWFB53pOz1pE\nGxDTxO3pc+hiSauAL6QjySkMF1+ctYKWJ2a62ZSS3d0INapnRXPqn2ZW8OZnM4gxtItKtp8HHqbP\nJGunTVmeBHf0idupEdPEPawZQpwCcsEFYe0GmhoxTdwzqp3vE2/IcZwGEuOoMBX4OMFJfhRwCjCF\nEH1vWHrSHKdBFPhdN+YddDQwJYkOj6TzgJvM7Pg0hTmOE1eD7gfsKNnfQZ9o8OWQdIWkTZLuLTnm\niZOcQIFrtWYSY6BXAb+WdF5Se97JS7lVqrEYeFefY544qQk0zYvo0kvD4qRGTIbtLxPyg/4xWT5i\nZv3O0jWzO4Cn+hz2xEktxILfCSZO7P/CPFDQGjumBoWQLOlpM/sG0C1pfI3P88RJLcSElT+Dn/jM\nwzSJGWb5IqEndyKwCBgMXEOIHF8zZmaS+iZOmgq8s2+SpT56TgZOBhg7dmw9Epw6mXLDIvjZEM+0\nnSIxNeh7gVnAVtiZWbvW4ZW6Eiclz/fkSVUoxCyWgjY3syDGQHeUpgmUtFcdz+tNnATlEycdA1wl\n6fV1PMNpIo/+MToTpVMDMQZ6naRLgVcmEf6WEzF5W9IPgJXAREndSbKkC4EjJW0Ajkj2d2JmDwC9\niZMOGNhHcZx+KGDNHeOL+/UkFtHThPfQL5hZZ0S5D1U45YmTHCeSqgYqaRCwPHGY79confbils9+\nDYCPDqRQAWuxLKlqoGb2gqQXJb3CzP7cLFFObTS7g2jLvvs39XntSIwv7hZCEt5Okp5cADP7VGqq\nnEJw0IqbAVgAnH5k1c73/HD7fDjsrKxVRBNjoEuSxXF24eBlPwBg/bTpGStpXSoaqKSxZvY7M4vx\nu3UcJwWqDbPsTJAk6YYmaHHqIGsHhajnewfRgKlmoCrZnpC2EMdpGgX6oahmoFZh23GcJlGtk+iN\nkp4m1KRDkm2SfTOzvVNX5+SaZZ//xi77CzrXF6c3tyBUNFAz8yzaTlW2v8ITDKRN7HxQx3kZk25d\nwqRbfQQuTdxAW4CsenAn3Xojk269Me7iAnXM5Ak3UKc9KcgPhhtowcl6/LMvedNTdNxAnfTJa22V\nV10luIEWmLzWVrvoKoAR5BlPI+jUzI8uWJi1hPrJ+ewWN9CCkofa8/mOIRXPLehcz+m7uwt3vXgT\n16mZg5dey8FLr614fuXGvsnZc0qOm+FuoAUkD7UnwEF33MJBd9xS9twhvytY8zenRporA5X0LkkP\nSvqtJM/bUlBKjbMwtSjk0khzY6BJgLJvA0cRovp9SJJH9+tDXmrPSpSrOd1IaydPnURvBX5rZhsB\nJP2QkGzpvkxV5YC8G2Uv1Zq1vUZ66IThzZJTOznq2c2TgY4CHi3Z7wb+LiMtmVMEo9x7+++Bgb1v\n9q1Nc2uw5WrSDIw2TwYaRWnyJODZ0gTBOWYEkPdsbQPWeEZvuX+5qJmfLcPv8uyBXNyQvIx5MtDH\ngDEl+6OTY7tgZguBhQCSusxsanPk1U4RdNaqsdmfrQjfJQSdjbhPbjqJgN8AB0oaL2kP4IOEZEuO\n07bkpgY1s+clzQX+GxgEXGFm6zKW5TiZkhsDBTCzm4GbB1CkKKPhRdBZq8Zmf7YifJfQIJ0KqT8d\nx8kjeXoHdRynD7k2UEn7SOqUtCFZv6rCdVdI2tR3yCW2fBN1lnVllHSepMckrUmWhiY76c+FUoFv\nJufvljQltmyOdD4s6Z7k+2tID2odOl8naaWkZyWdOZCyL8PMcrsAXwPmJdvzgK9WuO4fgSnAvbWU\nb4ZOQsfXQ4Qo/XsAa4FJybnzgDNT0lbxuSXXTAd+Soh5fAhwZ2zZPOhMzj0MjGjC/2SMzn2BtwBf\nLv271vJ95roGJbj69SZvuhI4utxFZnYH8FSt5RtAzHN2ujKa2Q6g15UxbWKe+x7gKgv8CnilpP37\nKyvpNZJukLRa0gOS3pqRzmbSr04z22RmvwGeG2jZvuTdQPczs8eT7T8A+zW5fCOfU86VcVTJ/ieT\nZtsVDW6K9/fcatdULCtpd0JttsjM3kRowdyfkU4I6UmWS1qVeJulRYzOhpXNfJhF0nLg1WVOnVO6\nY2YmqeYu53rLp6zzu8CXCP9kXwIuAv65Fp1N5GjgfjNbBmBmz2Ss5+/N7DFJ+wKdkh5IWlaFJnMD\nNbMjKp2T9ISk/c3s8aQps2mAt6+3fCN1VnRlNLMnSu51GbCsVp0DeW7ENYOrlJ0M/ConOjGz3vUm\nSTcSmpNpGGiUS2qjyua9ibsUODHZPhH4cZPLN/I5FV0Z+7xHvRdo5ASAGBfKpcCHk17SQ4A/J032\namX/ALy+9waSRmalU9JekoYlOvYC3kljv8OB6mxc2bR7versMRsO3AZsAJYD+yTHXwPcXHLdD4DH\nCS/l3cBJ1cpnqHM6sJ7Qk3dOyfGrgXuAu5M/2P4N1vey5wKnAKck2yJMln8o0TE1QvPQROs6YA0w\nKyudhF7RtcmyrlRnSn/v/nS+Ovk/fBr4U7K9d7Xvs9LinkSOk2Py3sR1nLbGDdRxcowbqOPkGDdQ\nx8kxbqCOk2PcQNsISVvKHDtF0oeT7dcls0FWS3qzpDnNV+mU4sMsbYSkLWY2tMr5ecDuZnaBpHHA\nMjN7Q7P0OS8nc1c/J1sknQdsIQQIPw14QdI7gCeAAyStATrN7DPZqWxf3EAdIMSDknQJsMXMvp7U\noG8ws8nZKmtv/B3UcXKMG6jj5Bg3UKcSfwGGZS2i3XEDbS/2lNRdspxR6UIz6wH+V9K9kv6tiRqd\nEnyYxXFyjNegjpNj3EAdJ8e4gTpOjnEDdZwc4wbqODnGDdRxcowbqOPkGDdQx8kx/w+uwgC4vjkZ\ndQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x121faa898>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "reach = 1\n",
    "true_lift = 0.05\n",
    "conv_ctrl = 20000\n",
    "split = 0.5\n",
    "ctrl_scale_factor = (1 - split) / split\n",
    "\n",
    "h0_dist = get_lift_dist(reach, 0, conv_ctrl, ctrl_scale_factor)\n",
    "h1_dist = get_lift_dist(reach, true_lift, conv_ctrl, ctrl_scale_factor)\n",
    "\n",
    "crit_val = get_crit_val(reach, conv_ctrl, ctrl_scale_factor, alpha=alpha)\n",
    "\n",
    "power = (100 - st.percentileofscore(h1_dist, crit_val)) / 100\n",
    "\n",
    "fig = plt.figure(figsize=figsize)\n",
    "ax = fig.gca()\n",
    "\n",
    "ax.hist(h0_dist, bins=100, label='$\\mathrm{H_0}$', alpha=0.5)\n",
    "ax.hist(h1_dist, bins=100, label='$\\mathrm{H_1}$', alpha=0.5)\n",
    "plt.axvline(crit_val, color='r', ls='--')\n",
    "ax.set_xlim([-0.1,0.1])\n",
    "plt.draw()\n",
    "\n",
    "xlabels = [x.get_text() for x in ax.get_xticklabels()]\n",
    "xlocs = list(ax.get_xticks())\n",
    "\n",
    "xlabels.append(r'$c$')\n",
    "xlocs.append(crit_val)\n",
    "\n",
    "ax.set_xticklabels(xlabels)\n",
    "ax.set_xticks(xlocs)\n",
    "ax.get_yaxis().set_major_formatter(matplotlib.ticker.FuncFormatter(lambda x, p: thousands_tick_formatter(x, p)))\n",
    "ax.legend()\n",
    "ax.set_xlabel('Lift')\n",
    "ax.set_ylabel('Frequency density')\n",
    "\n",
    "plt.savefig('img/single_cell_dist.pdf', bbox_inches='tight');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.999455"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "power"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Power vs audience size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOoAAACrCAYAAABokACHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXd8VGX2/99nJpMGSaiB0KvSCcWKIgo2LGBF1y6o+1vr\n7vq17bLq6rrqqguuri7qIgKLdS2LiiBSRAUBwdCbdAIJpJI65fz+eC4wQMokZFpy36/XvObOc+/c\n58xkPnnaOecRVcXGxiaycYTbABsbm+qxhWpjEwXYQrWxiQJsodrYRAG2UG1sogBbqDY2UUDQhCoi\n/xaRLBFZ7VfWTETmiMgm67mp37lHRWSziGwQkQuDZZeNTTQSzBb1beCiY8oeAeaqandgrvUaEekF\nXAf0tt7zTxFxBtE2G5uoImhCVdWFQM4xxaOAKdbxFGC0X/m7qlqmqluBzcCpwbLNxibaiAlxfa1U\nNdM63gu0so7bAov9rttllR2HiNwJ3AnQqFGjQT169DjuGk9WFp79B8Dnqyu7bWxqzJqy0v2q2rIu\n7hVqoR5GVVVEauy/qKqTgEkAgwcP1mXLlh055/Wy5+FHKJg5k6Sxo2g0dCgxLVviiIsDpxMROXIj\n/2Og1FPGvB3fsCjzO3YV7gIg3hlHq8TWtExoQUpcE5JiG5MQk0hcTCyxjjhczhic4sQhDhzixIkD\nBBziQJDj6gBMeQ2p+TuiHAXcxVBWYB7lReAuAU+peXaXgqfEXOMuNeXlReZZa/nP2ekCh/VwOkGc\nIA5wxIDDCRIDTgfgAIfDOu8058UBDgHE731OBj7w3va6+kpCLdR9IpKmqpkikgZkWeW7gfZ+17Wz\nympE0Q+LKZg5kxa/+Q0t7r3naGFWgqry8eaPmfjTRHJKc0jvns7ozjdxWuvT6NKkCw6xJ8brlJI8\nyN8JBXvMc94O8yjYAyW5UJxjntVb8fudcRCfDIkpkNAU4ltAXBIkNIP4FHAlQmyieXYlgDMWYhtb\nZQkQm2SOYxKMOJ0uiImv8J/qCfPAe3V2q1AL9TPgFuBZ6/lTv/L/iMhLQBugO/BjTW9evGQxxMTQ\nfNzYgESaX5bPHxb9gQW7FjAwdSATz51Iemp6Tau1OZai/ZC1DnJ+gdyt1vM28yjNP/pahwuatIeU\ndpDay4gvsRk0SoXGqZDUGhJbGHHGpxixNUCCJlQRmQEMA1qIyC7gcYxA3xeRscB24FoAVV0jIu8D\nawEPcLdqZf9SK6do8RIS+vfHkZhY7bWlnlLunns3aw+s5ZFTH+H6HtfbrWdN8ZRD1hrI/BkyM4w4\nD2yCouwj1zhioElHaNYZ2p0KTTsaUSa3g5S20Li16UraVEnQhKqq11dyangl1/8F+Ett6/MWFFC6\nZg0tfv3rQGzjkW8fISM7gxeHvcj5Hc+vbbUNB58X9q2G7T/A7uWwNwMObAafx5yPS4FWveCki6Bl\nD0jtCc27GkE6wzYVUm+oN99g8bJl4PORePpp1V770aaPmLtjLg8OftAWaWV43bDzR9g02zzvXQXl\nheZcUhtI6w8nj4TWfaFNOjTtHJxxng1Qj4RatHgxEhdHQnrVY8x9Rft4cdmLnNr6VG7udXOIrIsC\nVE3XdetC2PYt/LLACNPhgjYDoP8YaH8adDzTdF1tQkq9EWrxj0tJGDgAR2xsldc9v/R5PD4Pj5/x\neEATTvUaVdj9E2S8C2s/g4N7TXmTDtD3Kug6HLoMMxM5NmGl3gjVnZlJ4qBBVV6zvWA7c7bPYWzf\nsXRI7hAiyyIMrwd2/AAbZ8G6z8zSiDMOTroQul8AXc4xQrWJKOqFUNXnw1dQgDOl6v/8U9dOJcYR\nw696/CpElkUQBZmwYhosnwwFu836Yudz4JyHoedlZunDJmKpF0L1FRaCKs6Uyn9s+WX5fLr5U0Z2\nHknLxDrx6op8VGH7d7D4NdjwpXEi6HIuXPgMdBsBcY3DbaFNgNQLoXrzzSK6I7lyoX6y+RNKvaXc\n3LsBTCB5ymDVB7DkdTNbm9AUzrgbBt1qlkxsoo56ItQCgCpb1NnbZ9OzWU9OanpSqMwKPeVFsOzf\n8MOrUJgJLXvCZROh77XGbc4maqkfQi0wLWplY9S9RXvJyM7g3gH3htKs0FFeDEvfgO9ehuL90Ols\nGPUqdD3PXtusJ9QLofryDwm14hb1mx3fADCiw4iQ2RQSvB74aQoseA4O7jPjz2GPQofqnT5soot6\nIVRvgen6OpIrblG/3vE1XVK60KVJl1CaFVy2fguf/x72b4AOZ8A1bxtnBJt6Sf0Qal7lLWpeaR7L\n9y1nXN9xoTYrOOTvhrl/Nk4KTTrCmOnQ4xK7i1vPqR9CLShA4uNNgPgx/JT1Ez71MaTNkDBYVof4\nvPDDKzDvr2aZ5ezfw9kP2pNEDYT6IdT8PJyVdHtXZq3E5XDRu0XvEFtVhxzYAh//Gnb9CCdfAhc9\nA007hdsqmxBSL4RqvJIqnkhakbWCXs17Eec8vrWNeHw+WPYWzPmT8SS68k3oe7XdzW2A1AuhevPy\ncVSwNFPmLWPNgTXc0POGMFh1guTtgE9+YyJZuo2Ay1+B5LRwW2UTJsIiVBH5LTAOk8ZqFXAbkAi8\nB3QCtgHXqmpuIPfzFhTgand86NXaA2tx+9zRl15l1zKYcZ3xMLr8FRhwo92KNnBCngNDRNoC9wGD\nVbUP4MQk364wOXcgeAsKKhyjrshaAUB6yygRqqrxLHr7EohtBHd8AwNvskVqE7a9Z2KABBGJwbSk\ne6g8OXe1ePPzKxyjrshaQcfkjjRPaH7iFgebskJ470aY+VuzLjr2a2jRPdxW2UQIIReqqu4GXgB2\nAJlAvqrOpvLk3EchIneKyDIRWZadnY2Wl6PFxRW6D649sJa+LfoG54PUJXk74a0LTYTLBU/Djf+F\nxg0kwscmIMLR9W2KaT07Y1KDNhKRG/2vUVXFjF+PQ1UnqepgVR3csmXLI15Jx7So+WX5ZBVnRb4T\n/r418Nb5JsftDR/AmffaWflsjiMcv4gRwFZVzVZVN/Bf4Eys5NwAxyTnrpJDQnUeE+K2MXcjQGQL\ndeePMPlic3z7LOhWYYJGG5uwCHUHcLqIJIpJWjQcWMeR5NxwdHLuKqnMfXBT7iYAujeN0HHetu9g\n6hWQ2BzGzoZWUeyQYRN0Qr48o6pLRORD4CdMsu0VmL1kGlNBcu7qqCzEbWPuRlLiUmiZEIFjvVUf\nwqd3m9xEt/zPZIO3samCsKyjqurjmMz5/pRRSXLuqqgsxG1T3ia6N+keeZkGF02Arx+HjkPg2neg\nUYtwW2QTBUT9rMWh7A7+IW4+9bE5d3PkjU/XfmZE2vtKuOkTW6Q2AVMPhGq1qH5C3X1wN8We4sga\nn+5ZYRzr2w6G0a9BTNX5h21s/Il+oRYW4GjcGHE6D5dF3ETS7uXwzigzcTRmGrjiw22RTZApKHXX\n6f2i3ilfy8qR+KN/+JvzNgPQrUm3cJh0NHtWwjtXQEITuHWm7VhfjzlY5mHp1hz+l7GHL1ZlVv+G\nGhD9QnW7EZfrqLKdhTtpmdCSRq5GYbLKImu9WYKJT4ZbPzf7gNrUGw6WeVi5I4+fd+WxYEM2y3fk\n4vUpSXExXDmwHX+tw7rqpVB3H9xNm8ZtwmSRReE+mHaV2dH65k8brEjdbje7du2itLQ03KacMKpK\nqduHVxW3x0eJ24tXldgSN6iHu4Z2YUi3Fgzq2JR4l9MWqj8VCXXPwT30b9k/TBYB7hJ493ooyYHb\nvmzQSa937dpFUlISnTp1irylsgBxe3zklbjZf7CMGK/vUEQJHRJdJMfH0L8wnxE9D9K5c+eg2VA/\nhOq3g5vH52Fv0V5Gdh4ZJoMUPr3HTCCNmWb2Dm3AlJaWRp1IfaoUlXkoKPVQVOah1O0FIDE2hrZN\nE0h0OXE65PBnSopvwfoD+4NqU7VCFZF2mHjRszFO9CXAauBz4EtV9QXVwmpQd/lRLereor141Uu7\npDDt4fndBFj9IZw33my+ZBMVInV7fOSXuiks9VBc5sGrikOExFgnrZPjSU5wEe9yVvjeUHy+KoUq\nIpOBtsBM4DmMo3w8cBJwEfAHEXlEVRcG29DK0PKju767D+4GoG3jtqE3ZsMs+PpJ49Bw9u9DX79N\nwKgq5R4fZR4fOUXlh5dT4mKcNEl00TjeRVJcDA5HZPyTqa5FfVFVV1dQvhr4r4jEAmHdTPPYMWrY\nhJq1Dj4aB2n9zHYSUdCKNETKPF5yi8rJLXbj9prOYIzDQWpSPE0SK281w02VQq1IpCLSFUhU1VWq\nWg5sDpZxgaBuN47EI7ltdxXuwilOWjcKoaN7eRG8e4PJsXvdDDvXbgShqnh9SmGph9zicg6WeRAg\nKd5Fq+R44mIcJLicEdNyVkaNJpNE5DGgG+ATkThVvSk4ZgXOsZNJe4r20LpRa2IcIZwn++oxyPnF\nRMKkhKHL3YD517/+RUZGBq+++upR5WUeL9kFZeSXuPGqyUEQG+OgVXI8TRNjiY2JLqe86sao9wGv\nqqrXKuqvqmOscxnBNi4Qjp1M2l24O7Td3vWfw/K3Ycj90Pns0NVrA8CqVavo29ek2yn3+Cgu95BX\n7Kaw1A0iNE1wEedykhhrHtEwsVUR1TU7B4BZIvIPVf0MmC0iszA+wl8F3boAqGgyaUjbEG1fkbfT\n5N5N6w/n/jE0ddocxc8ZGVwy+mo2ZRVSUm7akxiHg5ZJcTRvHIfLGV0tZ2VUN0adLiIfAQ+KyDjg\nT8AMwKWq+bWtVESaAG8CfTC5kW4HNlCLvL7+k0mlnlKyS7JD06J6PfDRWLMnzNWT7WiYEKJqxpwH\nisrJyFhFcpvO+HyQlpJAozgn8S4njihtOSsjkIFcV+B9jLCessrGA7UWKjARmKWqV1szx4nAY5i8\nvs+KyCOYvL4PV3cjf6FmFhlH6JAIdf4zsHMJXPVWg/Y8CiWlbi/ZhWUUlrrx+JT9e/eQlJREete2\nUd2tDYTqxqhvA26MkHar6h0iMgB4Q0SWquqfa1qhiKQAQ4FbAayZ43IRGQUMsy6bAswnUKFak0n7\nivcBBH/G95cF8O1LJoN936uDW1cDp9Tt5WCZh+IyL/kl5ThESE5wkZzgYufPW0nv349GcUd+xnv3\n7mXMmDFccsklrFmzhjPPPJM5c+bwxBNP0KdPnzB+khOjuhZ1gKr2BxCRFQCqugK4zBJWbegMZAOT\nRaQ/sBy4nxrk9QXuBOjQoQPavMXhFjW7OBuA1MTUWpoWAKX5JgC8eTe4+Png1dOAKXV7KSh1U1Di\nobjcA5hxZ/PGcaQmxRFjjTv9J5IOsXLlSm644QbuvPNORo8ezR133EGTJk3Yvn17VAu1upH2LBH5\nSkS+Af7jf0JVA8oSWAExwEDgNVUdABRxzPYVNcnrq+VHZn2zS4xQg5rQ7OsnoTATrviX2XbCpk5Q\nVUrcXnYcKGbjvkL25pfiUyUtJZ4erZPomZZEmyYJh0UKRqiTJk2iU6dOdOrUiTPOOIOVK1cydOhQ\n3G43zZs3x+FwsHr16uMEHW1UN5n0sIgkAz5VPVhHde4CdqnqEuv1hxih7hORNFXNrEleXzyeo1rU\nRq5GJLqC5HCw/QezDeLpv4F2g4JTRwPChI15yS9xk1/ipszjwyFCalI8zRvHVjtjO3369OPKxo4d\ny0knnURGRgY9e/YEYNu2bXToEFYHuhOmujHqjcB/KnO8t7yU0lR1UaAVqupeEdkpIier6gZM5sG1\n1uMW4FkCzetrLWT7t6hBa01LC+Dju0yKz3P/EJw6Gghur4/sQuOM4Pb6EIRGcU5aNI4jOcF1Qksq\nb731FgDp6emkp5vIpalTp9aJ3eGkujFqc2CFiCzHjCWzMU753YBzgP3UYNc1P+4Fplszvr9gtl10\nUNO8vscKtTiblolBEuqXD5ttJ277EuIaB6eOesyhoOucojJyi92oQnJCDEnx8STHxxzVpbU5nuq6\nvhNF5BXgPGAI0A8T5rYOuElVd9SmUlVdCQyu4FSN8vrqIaFas77ZJdn0aRGECYNNX8PP/4Gh/wcd\nTq/7+9djPF4fmfkl5Je4Kff4EBGaJLhITYojLkId4CORatdRLffBOdYjsvAdaVFVlezibFIT6njG\n1+eFOeOhWRcY+lDd3rses3FfIX/7agPXdnNAYTmN42No2TiOlASX3XrWgujO8ODX9T3oPkipt7Tu\nu74rp0PWWrhmiu19VA2lbi9Lt+UwbfF2Zq/dR+PYGMb1TaNHWlK9ceULF1EuVDPHJS7X4TXUOp1M\nKiuEb/4C7U6FXrVdNq7/ZBWU8sLsDXy6cg9lHh9NEl38ZlhXxp3Vhb07ttgirQMCScXiAK5W1fdD\nYE/N8GtRs0rMak6dtqgL/wYH98J10+1A8GPweH3MzMhk9tq9zN+Qjdvr49rB7TmvRypndm1BQqwZ\nf+4Ns531hUDGqD4ReQjj7xtR+E8m1XmLun8T/PBPSL8R2lU079UwKS738NWavbw6bwubsw7SOjme\nUeltuGtoVzq1sB1AgkWgXd+vReRBTHRL0aFCVc0JilWB4teiHvZKqosWVRW+fAhcCTDi2E3nGiYH\nyzy8Nn8zk7/bRnG5ly4tG/H6jQO5sHfriHeGv/3225k5cyapqamsXm2Sljz88MN8+eWXpKen8847\n7wAwbdo09u/fzwMPPBBOcyskUKGOsZ7v9itToEvdmlND/IVanE1iTGLdZMdf9z/Y8g1c9Cw0DqLf\ncBSwK7eYGT/uYMaPO8kpKufy/m248fSODO7YNOLTlxzi1ltv5Z577uHmm28GID8/n59++omMjAzG\njRvHqlWr6NatG5MnT2bWrFlhtrZiAhKqqgYvs/CJoAoIEusiuzC7bpzxy4tNapXU3nDKHSd+vyhE\nVflmfRaTFv7Ckq05OASG92zF3ed2I719k1rf98n/rWHtnoI6tBR6tUnm8cuq3q196NChbNu27fBr\nh8OB2+1GVSkuLsblcvHCCy9w77334jommXukEJBQRSQR+B3QQVXvFJHuwMmqOjOo1lWHKogcblHr\npNu7+J9HPJCc0T0pXlNUlUWb9/OPbzbz49YcOjRL5MELTmL0gLa0a1p/ErYlJSUxcuRIBgwYwPDh\nw0lJSWHJkiWMHz8+3KZVSqC/xMkYF8Izrde7gQ8w+X7DhqqCWJNJdeGVVJIH378MJ10MHc+s/vp6\nwi/ZB3nnh+0s2JjN1v1FpCbF8dSo3lx3aoc6XVqpruULJQ899BAPPWQcWMaNG8ef//xn3nzzTWbP\nnk2/fv344x8jK7VOoELtqqpjROR6AFUtlkiYQfAbo+aU5tA8vvmJ3e+HV0y86bmP1YFxkU9RmYdp\ni7fz0pyNAJzepTn/75yujBrQhriYhuHet2LFClSVk08+mUcffZSvvvqK2267jU2bNtG9e4Tsr0vg\nQi0XkQSsGFEraqYsaFYFiiVUt0MpchfRLL5Z7e91MAsWvwa9Rpsk2vWYXbnFvDpvC5+u3E1xuZcR\nPVvxzBV9SE1ueBssjx8/nkmTJuF2u/F6TXI0h8NBcXFxmC07mkCF+gQwC2gvItMxDvq3BsmmwLGE\nWuAzX2rT+Ka1v9eC58BTavaMqafkFZcz4etNTF+yHUEYPaAN1wxuz+COTSN+ieVEuP7665k/fz77\n9++nXbt2PPnkk4wdO5ZPPvmEwYMH06aN2aIzPT2dvn370q9fP/r3D+NugBUQ6KzvbCvU7XRAgPtV\nNbjbVwXAIYeHPJ+Jaa+1UPdvgmWTYfDt0CICdimvY3YcKGb6ku28u3QnhaVuxpzSnvuGdyctJSHc\npoWEGTNmVFg+evRoRo8effj1Cy+8wAsvvBAqs2pEoLO+04AFwLequj64JtWAQ0L1GqHWeow690nj\n3HBOtbnUooZDM7hvLdrKgo3ZOEQY0TOV355/Ej1aJ4fbPJsaEmjX9y3Mtov/sManK4CFqjoxaJYF\ngiXUHG8hUMsWdccS4+Bw7h+gcRBzLYUIVWXBxmwmfL2JlTvzSE2K477zunP9qR1ondLwxqD1hUC7\nvvNEZCFwCnAu8GugNyY/b60QESewDJOG9FIRaUZNE3BbQj3gMSmGayxUVZjzJ2jcCs64u/rrI5yf\nduTyl8/XsXx7Lm2bJPDMFX25alDbBjODW58JtOs7F2gE/AB8C5yiqoElH6uc+zGZIg71wx6hpgm4\nVcHhINedT4wjhiRXUs0s2DgLdi6GS/8e1RkFN+wt5OW5m/h8VSYtk+J4enQfrh3cPuo2QrKpnEC7\nvhnAIMwWFPlAnoj8oKoltanU2sX8EuAvGI8ngJon4FZFXC5yy3JpFtesZjOXqjD/WWjaCQbcXBPz\nIwJVZfn2XCYt/IU56/aR6HJy33nduOucrkclpLapHwTa9f0tgIgkYZZlJgOtgbha1jsBeAjwbwJr\nnIC7e7NmSOs0ckpyat7t3TQbMlfC5a9ElaugqvK/jEwmLdzC6t0FpCS4uOfcbtw+pDNNG9kZKOor\nAfWNROQeEXkPM4k0Cvg3cHFtKhSRS4EsVV1e2TWBJuBOiIszXkllNRSqqlk3bdIB+l9X048QNlbt\nyueq177nvhkrKHX7eGpUb75/5Dx+f8HJtkgrYefOnZx77rn06tWL3r17M3GimVZ5+OGH6dev3+GI\nGjBhbhMmTAiXqVUSaFMSD7wELFdVzwnWOQS4XERGWvdNtpZ/ap6A+1DXtzSXti1qsDHU1gWwezlc\nOgGckRkt4c+u3GJeX7CF6Ut20LxRLM9f1Y+rB7WLmjCzcBITE8OLL77IwIEDKSwsZNCgQQwbNqze\nhrm9YO0T82trHPitqv5cmwpV9VHgUQARGQY8qKo3isjfqEUC7kNCrdEa6ncTzUxv/+tr/gFCyIa9\nhUycu5FZq/ciItx8ekd+d8HJpCRE/j+XCvnyEdi7qm7v2bovXPxspafT0tJIS0sDTNRMz5492bFj\nR70Nc7sPMy78r1U0TUQmqeo/6tCWZ6lNAm5XDAfdBwPv+mb+bILChz8OrshcV9y2v4i/f72Rz37e\nQ+PYGO46pys3nd6RNk0ahidRsNi2bRsrVqzgnHPOYd26dVEV5oaqVvvAzPo28nvdCMgI5L3BfPRt\n2VI3XDJS+7zdR9/f8L4GxAe3qf6lrWpxbmDXh5A9ecX6yEcZ2uXRz7XHH7/UZ79cp7lFZeE264RY\nu3ZtuE1QVdXCwkIdOHCgfvTRR8edGzt2rC5fvlzfeOMNveaaa/Spp56q8f0r+pzAMq2j33qgY1QB\nvH6vvVZZeFHFG2PMaBYXQOTM/s2w5mM44x5IqH2mgrqm1O3ljYW/8Or8zXh9yo2ndeDu87qRmhSZ\nLX604Xa7ueqqq7jhhhu48sorjzpX38LcJgNLRORj6/VojFtheFHFbTndNEsIQKiL/g7OWDjz3uDa\nFSBen/LJit28NGcju/NKuLhPax4b2ZP2zepPNoVwo6qMHTuWnj178rvf/e648/UqzE1VXxKR+cBZ\nVtFtajY0Di+quB1mFadpXDVj1NztkPEunDIu7AnLfD5l7vosXvhqAxv2FdK3bQp/u7ofZ3ZrEVa7\n6iPfffcdU6dOpW/fvod3d3vmmWcYOXJk/QlzE5F4jF9vN2AV8E898eWZusNfqNVNJi1+DRA4877g\n21UJqsrnqzJ5ac5GfskuolPzRF751QBG9kmzl1qCxFlnnXU4HPJY6lOY2xTAjfHvvRjoCURM0lNV\npdzhI0ZiSI6tInSrtABWTIPeV0BKDdZb65DdeSU888U6Ps/IpGdaMhOvS2dk3zR7uwebgKhOqL1U\ntS+AiLwF/Bh8k2qAKqVOL8lxyVX7+f48A8oL4bRfh842i+zCMp6ftZ6PV+xGBP7vwpO5a2gXe0cz\nmxpRnVDdhw5U1RNx6TpUKRNv1a2pzwdL/gXtToF2g0Jmmsfr471lO3l+1gZKyr3ceHpH7hjahbb2\nWqhNLahOqP1F5FDGZAESrNeCcckNb6qAQ0KNq8KMTV9BzpaQZhactyGLp2euZUt2Ead2asYzV/ah\nW2oNQ/BsbPyobsfxyI44VqVUPCTFViGC716GlPYmu2BQTTFhZ68v2MLX67Lo0qIR/7ppEBf0alWv\nE4fZhIboie+qALWEWmnXd9cy2PE9XPjXoIay7ckr4bfvrWTJ1hyS4mN45OIe3D6ksx24bVNnRLVQ\n8SlFUl65UL//B8SlwMCbglK9qvLZz3t4/LM1uD0+nry8N9cMbkdibHR/rfWNTp06kZSUhNPpJCYm\nhmXLltXb3dwiE1VKcFcs1LydJmnZGXdDXN2PD5dvz+GZL9azfHsu/dqlMGFMOl1aNq7zemzqhnnz\n5tGihXEoqbe7uUUs1jpqq4rGqMveAhROrdsd2XKLyhn/6WpmZpj8RHZsaOA89+NzrM+p22yzPZr1\n4OFTa5bmNRp3c4v6QZTHKce3qO4SWD4FelxisjjUAarK3HX7uHjit3y1Zi8PjOjOgv8bxrWntLdF\nGuGICCNGjGDQoEFMmjTpqN3c0tLSDoe5+XspRRrR3aICHifHz/qu+hBKcuDUu074/qrK4l9ymPD1\nRpZszaFry0a8ecsQ+rRNOeF7NzRq2vLVFYsWLaJt27ZkZWVx/vnn06NHj6jbzS3kLaqItBeReSKy\nVkTWiMj9VnkzEZkjIpus54AiwT1Ojl9HXfomtOwJnc6q+E0Bsq+glOvfWMz1byxmS3YRT43qzawH\nhtoijTLatjVuo6mpqVxxxRX8+OMRBzv/MLcPPviA999/ny1btrBp06ZwmVsh4ej6eoDfq2ovzF42\nd4tIL47k9e0OzLVeV3+zY1vU3ctNdsFTxsIJrF8u2JjNpf9YRMaufJ68vDeLHj6Xm87oZPvmRhlF\nRUUUFhYePp49ezZ9+hzZR3f8+PE89dRT9SPMrS5RkxI00zouFJF1QFtqk9cXq0X1H6Mu/Te4GkG/\nMbWyb3PWQZ74bA2LNu+nS4tGTBt7Gie3tr2KopV9+/ZxxRVXAODxePjVr37FRRddBBBVYW5SWQhQ\nSCoX6QQsxCT23qGqTaxyAXIPvT7mPYfz+vaOix809N4u/PXp70mJS4GSXHixJ/QfA5fVbLcNr095\n49tfeGnK2y21AAANK0lEQVTORhJcTu4f3p0bTu9gbwdxgqxbt46ePXuG24ygU9HnFJHlqjq4Lu4f\ntskkEWkMfAQ8oKoF/m52qqoiUmleX2ASQJ/4BPU4obHLWr9c9SF4SmDQbTWy5Zfsgzz4wc/8tCOP\nC3u34qnRfew0KDYRRViEKiIujEinq+qhzIY1z+sLxMTG43RYrd6KqSZ9ZJv0gOwodXv51wKTqyg+\nxsHE69K5vH8b2zfXJuIIuVCtbu1bwDpVfcnv1GfUNK8vEBtvbe6UmWFSgV78t2rfcyjTwl+/WM/u\nvBIu7ZfGny7tRWqy3YraRCbhaFGHADcBq0RkpVX2GLXJ6wvExlmJwFZMA2cc9L26yutL3V7Gf7Ka\nD5bvomdaMs9f3Y8hdq4imwgnHLO+i6g81ejwmt4vLqExeN2w6gPjiZRYeTbCxb8c4I+frGZz1kHu\nG96d+4d3x2l7FdlEAVHtmZTVKpbiTq1gyzzjiVTJkozH6+PZL9fz5qKttGuawDu3n8rQk6J/d3Gb\nhkNUr96XxiiJKc1g9YcQ3wS6nnfcNblF5dw6eSlvLtrKzWd0ZM5vz7FF2sC4/fbbSU1NPcrRIScn\nh/PPP5/u3btz/vnnk5trNrf/7rvv6NevH4MHDz7snZSXl8cFF1yAz+cLi/0Q5UL1qpckZwKs/xx6\njYKYo7ceXLUrn0v/sYgft+bw/NX9+POoPiTE2uuiDY1bb731uPC1Z599luHDh7Np0yaGDx/Os8+a\njaZefPFFvvjiCyZMmMDrr78OwNNPP81jjz2GwxE+uUR119enPpIL90H5waMmkco9Pl6dt5lX522m\nZVIcH/z6DPq3j5wtLBoqe595hrJ1dRvmFtezB60fqzof1tChQ9m2bdtRZZ9++inz588H4JZbbmHY\nsGE899xzuFwuiouLD4e/bdmyhZ07dzJs2LA6tbumRLVQAZL2b4FGLaHjEAByisq5a+oylm7LZXR6\nGx6/rLe9ya/Ncezbt+/wdoytW7dm3759ADz66KPcfPPNJCQkMHXqVB588EGefvrpcJoK1AOhJu/f\nDF2GgcPJ5qyD3P72UvYVlPLy9QO4vH+bcJtn40d1LV+4EJHDTi7p6eksXrwYgIULF5KWloaqMmbM\nGFwuFy+++CKtWrUKuY1RPUYFSC7Ohc5D+X7Lfq7853cUl3t4987TbZHaVEmrVq3IzMwEIDMzk9TU\no/cjUlWefvppxo8fz5NPPsnzzz/PHXfcwcsvvxwOc+uBUH0+fnL047bJS2mVHM/HvxnCgA4Bbmps\n02C5/PLLmTJlCgBTpkxh1KhRR51/5513GDlyJM2aNaO4uBiHwxHe8Le62mg1HI/+nZtqwQu9tdf4\nL3XEi/P1wMHo3vS3PhIJGxlfd9112rp1a42JidG2bdvqm2++qfv379fzzjtPu3XrpsOHD9cDBw4c\nvr6oqEiHDRum5eXlqqq6cOFC7dOnjw4cOFDXr19fYR2RspFxRBJTfpCFJSfTJDGWqWNPo5k9aWRT\nATNmzKiwfO7cuRWWJyYmMm/evMOvzz77bFatWhUU2wIluru+Pi9fFZ/MC9f0p3WK7VBvU3+JaqHm\nkEy7gRdyRtfm4TbFxiaoRLVQSxLbcPdlQ8Jthk01aBiziISCUHy+qBZq2yYJNIqL6mF2vSc+Pp4D\nBw7UW7GqKgcOHCA+PrhDL/tXbhNU2rVrx65du8jOzg63KUEjPj6edu3aBbUOW6g2QcXlctG5c+dw\nmxH1RFzXV0QuEpENIrJZRALK7WtjU9+JKKGKiBN4FbgY6AVcbyXntrFp0ESUUIFTgc2q+ouqlgPv\nYhJz29g0aCJtjNoW2On3ehdwmv8F/gm4gTIRWR0i26qiBbDftgGIDDsiwQaAk+vqRpEm1GpRvwTc\nIrJM6ygT+YkQCXZEgg2RYkck2HDIjrq6V6R1fXcD7f1et7PKbGwaNJEm1KVAdxHpLCKxwHWYxNw2\nNg2aiOr6qqpHRO4BvgKcwL9VdU0Vb5kUGsuqJRLsiAQbIDLsiAQboA7tCOtubjY2NoERaV1fGxub\nCrCFamMTBUStUEPlaigi7UVknoisFZE1InK/Vd5MROaIyCbrualVfquIvBIkW5wiskJEZobRhiYi\n8qGIrBeRdSJyRpjs+K3191gtIjNEJD7YdojIv0Uky3/tvrI6rXOPWr/PDSJyoV/5wZrWHZVCDbGr\noQf4var2Ak4H7rbqegSYq6rdgbnW62BzP7DO73U4bJgIzFLVHkB/y56Q2iEibYH7gMGq2gcz8Xhd\nCOx4G7jomLIK67R+I9cBva33/NP63daKqBQqIXQ1VNVMVf3JOi7E/DDbWvVNsS6bAow+9r0icomI\n/CAiJ7yvo4i0Ay4B3vQrDrUNKcBQzP62qGq5quaF2g6LGCBBRGKARGBPsO1Q1YVAzjHFldU5CnhX\nVctUdSuwGfO79belhWXLJdXVHa1CrcjVsG2wKxWRTsAAYAnQSlUzrVN7gVbHXHsF5r/rSFWtC3e2\nCcBDgP9ORaG2oTOQDUy2uuBvikijUNuhqruBF4AdQCaQr6qzQ22HRWV1VvkbFZFWwOfAn1T18+oq\niah11EhGRBoDHwEPqGrBoczqAKqqIuK/znUeMBi4QFUL6qDuS4EsVV0uIsMquibYNljEAAOBe1V1\niYhM5JjuZSjssMaBozD/OPKAD0TkxlDbcSwV1FkZLkw3+W5VXRDIvaO1RQ2pq6GIuDAina6q/7WK\n94lImnU+Dcjye8sWIAk4qY5MGAJcLiLbMN3880RkWohtANMq7FLVJdbrDzHCDbUdI4Ctqpqtqm7g\nv8CZYbCDKuqs6jfqAZYDFxIg0SrUkLkaimk63wLWqepLfqc+A26xjm8BPvU7tx24CnhHRHqfqA2q\n+qiqtlPVTpjP+o2q3hhKGyw79gI7ReRQVMhwYG2o7cB0eU8XkUTr7zMcM3cQajuoos7PgOtEJE5E\nOgPdgR+tcwrcDvQQkYcDqqWuMnmH+gGMBDZi/lP+IYj1nGV9sRnASusxEmiO6b5sAr4GmlnX3wq8\nYh0PwPyQu9ahPcOAmdZxyG0A0oFl1vfxCdA0THY8CawHVgNTgbhg2wHMwIyJ3ZjexdjK6rSu/4P1\n+9wAXOxXftB6jsO4y/6murptF0IbmyggWru+NjYNCluoNjZRgC1UG5sowBaqjU0UYAvVxiYKsIVa\nQ0SktYi8KyJbRGS5iHwhInW9iF5rRORyCW400ZtBDICoqt7RtalXRJ4QkQeDYVMosV0Ia4C1uP4x\nMEVVr7PK+mP8OzeG2BanqnqPLVfVzwhinilVHRese1fDaGAmZg30KEQkRlU9oTcpdNgtas04F3Cr\n6uuHClT1Z1X9Vgx/s+IjV4nIGAARGSYi8+VIDOd069qLROSDQ/exrjsUZ3qBFVXxk4h8YPkZIyLb\nROQ5EfkJuEZE7hMTJ5shIu9a1xyOvRSRTiLyjXV+roh0sMrfFpGXReR7EflFRK62ytNEZKGIrLQ+\nx9nHfgHWZxlsHR8Ukb+IyM8isthyND/2+sYiMtn6TjJE5Cqr/HqrbLWIPOd3/XH3FJEzgcuBv1m2\ndbXsmCAmJef9lX3WekOwPHrq4wMTA/n3Ss5dBczBxEa2wri5pWE8ifIxvp4O4AeMt1OMdU0j6/2v\nATdikkcv9Ct/GBNhAbANeMivzj1AnHXcpAIvnP8Bt1jHtwOfWMdvAx9Y9vTChAwC/B7Ly8v6HEkV\nfM75mDhQMB5bl1nHzwN/rOD654AJfq+bAm2sz97S+h6+AUZXdU/L5quPseOffq8r+6xPAA+G+7dz\nog+7Ra07zgJmqKpXVfcBC4BTrHM/quouVfVhXBA7qemqzQIuExNTeQnGT/R0jHi+E5GVGP/Rjn71\nvOd3nAFMFxM5UlHX7wzgP9bxVMvGQ3yiqj5VXcuR0KylwG0i8gTQV038bVWUY7qjYJzMO1VwzQhM\nkD8AqpqL+V7mq3Gq9wDTMXGugd7zEP7fRVWfNeqxhVoz1gCDavG+Mr9jL0fmBt4FrsWEYC2zhCHA\nHFVNtx69VHWs3/uL/I4vwYhgILDUEnxtbBI4HBg9FBPl8baI3FzNPdxqNVvHfK4ToSb3LKriXL3C\nFmrN+AaIE7P/DQAi0s8ay30LjBGT16gl5gf/YyX3OcQCjMjuwIgWYDEwRES6WfdvVNGssog4gPaq\nOg/TPU4BGh9z2feYaBuAGywbK0VEOgL7VPUNTCaJgdXYHwhzgLv96miK+V7OEZPhwAlcj/kuqqIQ\nE6ZWGTX6rNGGLdQaYP2nvwIYYS3PrAH+ions/xjTFf0ZI+iH1ISFVXU/L6abd7H1jKpmY8aZM0Qk\nAzOm7VHB253ANBFZBawAXlaTFsWfezFd2QzgJkzOpaoYBvwsIiuAMZj8SCfK00BTa9LoZ+BcNRkR\nHgHmYb6v5ar6aVU3wfwj+z8xmSW6VnC+pp81qrCjZ2xsogC7RbWxiQJsodrYRAG2UG1sogBbqDY2\nUYAtVBubKMAWqo1NFGAL1cYmCvj/grgIwbUosewAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1216915c0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "reach = 1\n",
    "min_lifts = [0.01, 0.02, 0.05, 0.1]\n",
    "convs_ctrl = np.linspace(100, 100000, 100)\n",
    "ctrl_scale_factor = 1\n",
    "\n",
    "fig = plt.figure(figsize=figsize)\n",
    "ax = fig.gca()\n",
    "\n",
    "plot_power_by_sample_size(ax, reach, min_lifts, convs_ctrl, ctrl_scale_factor, alpha, n_iter)\n",
    "\n",
    "plt.savefig('img/single_cell_power_vs_conv.pdf', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Power vs reach"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOcAAACrCAYAAACdbos3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXd4lFXah+9n0nsjhFRC702KIIgUpSnq7vop7qeiy9pA\nURFRV1wrCrqifnakuZZdEUUpC4IoKC41UgRCINSEAOllMkkmM3O+P97BjS6QkGTmnUne+7rmmpm3\nnd8k88x7ynN+R5RSGBgYeB4mvQUYGBicGyM4DQw8FCM4DQw8FCM4DQw8FCM4DQw8FCM4DQw8FJcF\np4gsFJFcEdlbY1u0iKwTkUPO56ga+x4XkUwRyRCR0a7SZWDgLbjyzrkYGPObbY8B65VSHYD1zveI\nSFdgAtDNec7bIuLjQm0GBh6Py4JTKfU9UPibzdcBHzhffwBcX2P7P5VSVUqpo0AmMMBV2gwMvAFf\nN5cXp5Q65Xx9Gohzvk4EttQ4Ltu57b8QkbuAuwBCQkL6du7c2UVSvQylUFVVOKxW1NmHzYay2cBm\nQ9nt4HCAkRDmNvZVVeYrpWLre767g/MXlFJKRC76q6KUmgfMA+jXr5/asWNHo2vzBhxVVZRv3oxl\n8xYsaWlUZWSgbHYw+UBgED7xCfjGxeEbE4NPTDQ+4RGYQkMwBQVjCgpE/P0RP3/nsx/i6wMmH8TH\nBCLaa5Nor0VATCC/1iAi5xZ3Iepzzm9QSlFWXcaZ8lzOlJ/mTMUZ8i355FcUUFRVSFFlETZlP+e5\nPvgQ7BdEkG8wQb6B+PsEEOATgL+PH/4+/viZ/PAz+ePn44uv+OFr8sXX5IOv+OJr8kXEhI/44ON8\nNplMCCZMJsGEDwIIAgJDr7nneEM+p7uD84yIxCulTolIPJDr3H4SSK5xXJJzm0ENlFJYtm2neOlS\nzOvX47BYEH9/gnr1InribQR06ox/2zb4t26NT2io3nIbhWpHNQeLDrIrdxfpBekcKj7EidITmKvN\nvxxjEhNxEXHEJ8QTH9qR9kEtaRHUgpigGKICo4gJjCEiIIJw/3CCfIPq96NSL+5p0NnuDs7lwERg\ntvP5qxrbPxGRuUAC0AHY5mZtHouy2ShZsZLChQupOnQIU1gY4VdfTdioqwgeMABTQIDeEhsNq93K\n7rzd/HjyR7ad3kZGYQZWhxWAmMAYOkZ1ZHy78aSEpZASnkJKWAqJYYn4mfx0Vt74uCw4ReQfwDCg\nhYhkA0+hBeUSEZkEHAduBFBK7RORJcB+wAZMUeo89ZJmRvmWrZyZNYuqQ4cI6NiR+FmzCB83FlNQ\nkN7SGo1cSy4bsjawIWsD209vp9Jeia/40jO2J3/s8ke6xnSlT8s+tApppbdUtyLePGWsKbc5HRUV\n5L78N4o++QS/pCRaTp9O2OhRbqySuZbCykLWHF3DqqOr2JO3B4Ck0CSGJg1lYPxA+rfqT6i/d1fN\nRSRNKdWvvufr1iFkcH6sx4+Tfd/9VB06RPTE24h96CFMgYF6y2owNoeN77K+4/ODn7Pl1Bbsyk6n\nqE5M7TOV4cnDaRfZrsn8+DQGRnB6GJa0NLLunYyIkPz++4RePkRvSQ3GbDXzxaEv+Dj9Y3LKc4gP\nieeO7ncwrs04OkR10Fuex2IEpwdRvmULWfdOxi8ujuT57+OflKS3pAZRai3lo/0f8dH+jyirLuOS\nlpcwY8AMhiUNw8dkJIDVhhGcHoJ5049kT5mCf0oyKQsX4htb77Fr3bHarSzJWMK7e96lpKqEkSkj\nubPHnXRr0U1vaV6FEZweQPmWLVpgpqaSsmghvtHRekuqFzaHjWWZy3hv93ucsZzh0vhLebjvw3SJ\n6aK3NK/ECE6dsWZnk/3Ag9odc/EifKOiaj/JA9lxegezts4isziTXrG9eHbwswyKH2R08DQAIzh1\nxFFZSfbUqeBwkPTWW14ZmFX2Kt7e9TaL9i4iITSB14a9xoiUEUZQNgJGcOrImRdepGp/OknvvI1/\nSoreci6aH7J/4MVtL5JVlsUNHW/gkX6PEOwXrLesJoMRnDpRsmIlxUuWEHPnnwkbPlxvORdFha2C\nv23/G0sOLiE1PJV5V81jUMIgvWU1OYzg1AFrVhann3qKoL59iX3gAb3lXBSHiw8zfeN0Moszub3b\n7dzf5378ffz1ltUkMYLTzSi7nZzHHgeTicSXX0J8vedfsObYGv76418J8g3ivSvf47LEy/SW1KTx\nnm9GE6Hwg79TkZZGwpzZ+CUk6C2nTiileG/Pe7y16y36tOzDK1e8Qmyw947DegtGcLoRW34++W++\nSejw4YRfe63ecupEtb2aZzY/w1eHv2J82/E8fdnTRjXWTRjB6Uby/u8NHFYrcY/O8IqhhjJrGdM2\nTGPLqS1M7j2Ze3re4xW6mwpGcLqJyoMHKV66lOhbb8E/NVVvObVypvwM966/l6PFR3lu8HNc3/76\n2k8yaFSM4HQT+W+8iSk4mBb33qu3lFo5aT7JpK8nUVxVzNtXvm0Mk+iELo7vIvKQiOwTkb0i8g8R\nCbyQ4bS3U5meTtm6dUTfdhs+kZF6y7kgR0qOMHH1RMqsZcwfNd8ITB1xe3CKSCIwFeinlOoO+KAZ\nSp/TcLopkPfWW5jCwoi+faLeUi7I3vy9TFw9EZvDxsLRC+neorvekpo1eq2V4gsEiYgvEAzkcH7D\naa+mMj0d8zfriZ44EZ/wcL3lnJeMwgzuXHsnIX4hfDj2QzpFd9JbUrPH7cGplDoJ/A04AZwCSpRS\nazm/4fSvEJG7RGSHiOzIy8tzi+aGkP/Ou5hCQ4m+7Va9pZyXrNIs7l53NyF+ISwavYjk8OTaTzJw\nOXpUa6PQ7pJt0GwwQ0TklprHKM117JzOY0qpeUqpfkqpfrEePiG58uBBytauJfq2Wz32rpljzmHS\n2knYlI15V80jPjReb0kGTvSo1l4JHFVK5SmlqoEvgMtwGk4D/MZw2mspePc9TMHBRN92m95Szkmu\nJZc/r/0z5moz866aR9vItnpLMqiBHsF5AhgoIsGijWiPBNL5j+E0/Npw2iuxZp+kdM0aIidM8Mge\n2pKqEu5edzcFFQW8e+W7dI3pqrckg9/g9nFOpdRWEVkK/IRmIL0Tbe2TUM5hOO2tFH34dzCZPLKt\naam2MPmbyZwoPcE7V75Dz9ieeksyOAe6JCEopZ5Cc4CvSRXaXdTrsZeWUvzZUsLHjsWvlWe5lFvt\nVh787kH2Fexj7rC5DIg3Vlr0VIwMIRdQ/NlnOCwWYu64XW8pv8KhHMzcNJPNpzbz/ODnGZEyQm9J\nBhdAr3HOJouy2Sj86GOCBwwgsKtnteNe/+l1Vh9bzUN9H+K69tfpLcegFozgbGTK1q7FduqUx2UD\nfZX5FQv3LuSmTjdxR7c79JZjUAeM4GxkCj74AL/WKYQOG6a3lF/IKMzguS3PMaDVAB4b8Jgx7ctL\nMIKzEancv5/K3XuIvuVWxOQZf9riymIe/O5BIvwjeGnoS/iajG4Gb8H4TzUixcu+RPz8iBh/jd5S\nAM3FYNrGaZyxnGHRmEXEBMXoLcngIjCCs5FQViulK1cSOnKkxyQdvLjtRbaf3s4LQ16gV2wvveUY\nXCSeUfdqAph/+AF7URER13tGL+iSjCV8dvAz/tT9T4xvN15vOQb1oNY7p4gkoc23vBwtUb0C2Aus\nAlYrpRwuVeglFC9bhk+LFoQO0X89zbQzaby49UWGJg1lap+pessxqCcXvHOKyCJgIWAF5gA3A5OB\nb4AxwCYRGepqkZ6OLS8P84aNRFx7re4+tGfKz/DwhodJCkti9uWzjXUwvZjavkmvKKX2nmP7XuAL\nEfEHvG+Rj0ameNmXYLMRecMNuuqodlTz8MaHsdgsLBi9gDD/MF31GDSMCwbnuQJTRNoBwUqpn5VS\nViDTVeK8AeVwULx0KcH9+hHQto2uWubvmc/uvN28NPQl2kW201WLQcO5qDqYiPwFaA84RCRAKeV5\nUy7cjGXbNqpPnCD2vim66vg572fe2/Me17S9hrFtxuqqxaBxuGBwishU4C2llN25qZdS6ibnvj2u\nFucNlCz7ElNoKGGjRummwWw189gPj9EiqAWPX/q4bjoMnCgFB1Y2+DK13TkLgDUi8oZSajmwVkTW\noHUkfd3g0r0cR0UFZevWETZ2DKbAQF00KKV4dvOznDSfZOHohYT7e6YdSrOh4DCsngGZ3zT4Uhfs\nrVVKfQyMB3qKyHIgDfg98D9KqUfqW6iIRIrIUhE5ICLpIjLIG31ry779FofFQsR4/dY9+TLzS1Yf\nW82U3lO4JO4S3XQ0e8oL4NtZ8PZAyNoGo19s8CXrkoTQDlgC3AVMAV4HghpY7uvAGqVUZ6AXmk2J\n1/nWlq5YiW9cHMED+utSfnZZNrO3zaZ/q/5M6jFJFw3NnjP7YNE4eLkdfP8SdL0e7tsOgyY3+NK1\ntTkXA9Vo3rInlVJ3ikgf4H0R2a6UevZiCxSRCGAocDuAs8fXKiLXAcOch30AbAAevdjruwtbURHm\nTZuIvu02XZLc7Q47T2x6ApOYmDV4FiYxkr3cis0KW9/R7paBETDsceg0FuIbz/KltjZnH6VULwAR\n2QmglNoJjHcGU31oA+QBi0SkF1pV+QEuwrcW7S5OSop+Q6ylK1aCzUbEdfpUaRftW8RPuT8xa8gs\nw87S3RxaB2seg4JM6HQ1XPt/ENKi0Yup7ed2jYh8LSLfAp/U3KGUqq87ni9wCfCOUqoPUM5vqrDe\n4FtbvGwZgV27EtjJ/c7o+/L38dbOtxidOprxbY28WbdReAT+8Uf42Jls8sfP4OZPXBKYUHsSwqMi\nEg44lFLmRiozG8hWSm11vl+KFpxnRCReKXXK031rK9PTqUpPJ27mTLeXbbaamfH9DGKCYnhy4JPG\nxGl3UFkC378MW94FH38Y+RQMmgK+AS4ttrY25y3AJ+dLbndmC8UrpTbVtUCl1GkRyRKRTkqpDDTH\nvf3Ox0RgNh7uW1u8bBni50f41ePcWm7NYZMFoxcQERDh1vKbHeUFsPsT+PF1KM+D3v8LI56EcPc0\nI2prc8YAO0UkDa1tmAcEomUJXQHkU79e1fuBj525uUeAO9Cq2B7vW6tsNkpXriJ0+HB8o9w72vPF\noS9YfWw1U/tMpW9cX7eW3ayw2+Df/wcb54CtElIvhz8ugUT3DlXVVq19XUTeBEYAg4GeaFPG0oFb\nlVIn6lOoUmoX0O8cuzzet7Z861bshYWEX3O1W8s9XnqcOdvncGn8pcawiSvJ2QUrH4ScndDlWq0X\nNk4fF8Vac2udqXvrnI9mT+nq1ZiCgwkd6r6ZcjaHjb/88Bf8TH48P/h5Y9jEFeQegH+/oVVjg1vA\n/yyGbr/TVZJhU3IRKKuVsnXfEDpypFvT9RbvW8ye/D28fMXLtArxLAd5r6eyFL55GnYsAN8gGHCX\ndrcM0t9qxgjOi6B882YcJSWEj3XfrI9DRYd4e9fbjGo9ijGpY9xWbrPg+L/hi7ugJBsGTobLp0OI\n55ig1cWmxATcoJRa4gY9Hk3pv/6FKSyMkCGD3VKe1W7liU1PEOoXyl8u/YtbymwW5OyCtMXw0wcQ\nmQKT1kKy560ZU5c2p0NEZqDl1zZbHJWVlK37hrBxYzH5+7ulzLlpc0kvTOe14a8ZtpaNQWUJrJoO\nPy8BnwC4ZCKMeg4CPNMxoq7V2m9EZDrwKVpGDwBKqUKXqPJAzBs2aDNQrnGPJ+364+v5OP1jbuly\nCyNTPL4T27Ox22DPp7DhRSjNgSsehYH3QpBnT3yqa3De5HyuOd1fAc1mKeTSVavwjY0luL/rZ6Dk\nWnJ5avNTdI3pyrS+01xeXpMmazusmAq5+yG+F/xhAaRcqreqOlGn4FRK6WuOozP20lLMGzYS9ceb\nER/XutmdXaavylbF7Mtn4+fj59Lymizl+fDt81rbMjwBbvwQuowHL0p3rFNwikgwMA1IUUrdJSId\ngE5KqYZ7MXgBZeu+QVVXE3616xMPPsv4jM2nNjPz0pm0iWjWv4n1w1IIW96Bre+CtRwuvQdGPOGx\n7coLUddq7SK09L3LnO9PAp8BzSI4S1evxi8picAePVxaTo45h7lpcxkYP5AbO3lk9qLnUl2pza/8\n4VWoKtHukiOehFj3zxpqLOoanO2UUjeJyM0ASimLNJPpELbCQso3bybmT39y6QwQu8POX//9VxSK\npy972phtcjEc/g5WTdOmdHUcCyOfhLhueqtqMHUNTquIBOGcY+mcjVLlMlUeRNnadWC3Ez7OtYkH\nr6a9ytZTW3l60NMkhia6tKwmQ5UZ1j6htSuj28FtX0HbYTqLajzqGpxPA2uAZBH5GC0J/nYXafIo\nSv/1L/zbtCGgc2eXlbHyyEo+2P8BN3e+mT90/IPLymkyWAq1dLvtC6DsNFw2FYY/AX76OCC6irr2\n1q51ThsbCAjwgFIq36XKPIDq06exbN9Oi8mTXVbNLKkqYc62OfSO7c2M/jNcUkaTwW6DHQvhu1lQ\nWQxth8MNi6D1IL2VuYS69tZ+BGwEflBKHXCtJM+hdNUqUMqli+G+ufNNSq2lzBw401h1+kIc/UHz\ng83dD22ugNEvQKvueqtyKXX9NixAWwLwDWd7cyfwvVLqdZcp8wBKlq8gsFdP/FNTXXL9PXl7WHJw\nCTd2vJFO0d7bq+hSik/Auqdg3xcQkQI3fQSdr/Gq8cr6Utdq7Xci8j3QHxgO3AN0Q/OfrRci4gPs\nQLPcvEZEotHSA1OBY8CNSqmi+l6/oVRmZFCVkeEyn6D8inwe+u4h4kPiua/PfS4pw6spztLS7fZ8\nCiZfuOIxGPIg+DXUMtl7qGu1dj0QAmwGfgD6K6UaasD1AJqjwtn1A86aSs8Wkcec73XzrS1Zvhx8\nfV3SS2tz2Ji2YRpl1WV8eOWHhhdQTezVWgLBdy+CskP/P8Nl90NEkt7K3E5dq7V7gL5Ad6AEKBaR\nzUqpivoU6lwt+2pgFlrmEYDHmEorh4PSVf8idPBgfKOjG/3683+ez87cncy5fI5RnT2LzQqZ62D9\nc5CXDh1Gw7iXIaq13sp0o67V2ocARCQMbQhlEdAKqK834GvADKBmTpXHmEpbtu/Advo04dOnN/q1\n9+Xv493d7zKuzTjGtXWve59HohSkLYJvntF6YCNSYMIn0Glcs2hXXoi6VmvvQ+sQ6ovWHlyIVr29\naETkGiBXKZUmIsPOdYxSSonIeU2lgXkA/fr1O+cxDaV05QokOJiwEcMb9brl1eU8+sOjtAhqYUye\nBig5CWtnap09ba6AQfdpSQS+7pkv6+nUtVobCMwF0pRStgaWORi4VkTGOa8b7hyq8QhTaYfVSunX\nawm7ciSm4OBGvfasLbPIKsti/qj5zbudWWXWFv3Z8i4oB4yYCUMeBh3WnPFk6vTXUEr9DagE7hGR\n+5xrnNQLpdTjSqkkpVQqMAH4Vil1C7AczUwadDSVNm/ciKO0lIjxjbvMwfLDy1lxZAX39LyH/q30\nWZVMd6wW2LEI3hqgGTV3/wPcnwZDHzEC8xzUtVo7Fa2d94Vz00ciMk8p9UYjapmNB5hKl65YiU9M\nDCGDGi/r5GjJUZ7f8jz94vpxV8+7Gu26XoNS2pDI138BSwEk9NGsJz3Qt8eTqGu19s/ApUqpcgAR\nmYM2rNKg4FRKbUDrlUUpVYDOptL2sjLMGzYQedNNiG/jZOtU2iqZ8f0MAnwCmH35bHxMrp2s7VEo\npTncff8SHNkAyZdqSQQpg5p9Z09dqOs3UAB7jfd257YmRdnatSirtdHS9c6ubZJRmMGbI98kLuSc\nHdBNk8z1sPElyNqimTSPmQMD7oTm9OPUQC5msvVWEVnmfH89Wkpfk6Jk+Qr8Wqc02qTqTw58wooj\nK5jSewpDk9znEK8rFUWw6mHY+zmEJ8LYl6DPreDfuJ1rzYG6jnPOFZENwBDnpjuci+g2GazZ2Vi2\nbaPFlCmNMgMlozCDV3a8wrCkYc2jnXnyJ9j5Iez9AqxmGD4TBk91+TJ5TZnalgAMRMujbQ/8DLzd\nCEMpHknxp5+CyUTkDQ2fT1lpq+TR7x8lMiCSZwc/27TXNjn8HXz3AmRv05Yz6HKNNl6Z0FtvZV5P\nbXfOD4BqtISDsUAX4EFXi3I3DquV4qWfEzp8GH6tGr4Wyextszlccpj3rnyPqEDP9katNyfT4NtZ\ncHi9ltUzZg70vhkCm/H4bSNTW3B2VUr1ABCRBcA210tyP2Vfr8VeVETUhJsbfK3PD37O54c+584e\nd3JZ4mW1n+Bt5KbD+mch418QFA2jntcW/zGqr41ObcFZffaFUsrWVE2niv75T/xapxByWcPGNtML\n0nlh6wsMih/ElN5Taj/BmzDnaj6wOz8E/zDN2e7Su73SctJbqC04e4lIqfO1AEHO94KWAht+/lO9\ng8qMDCrS0mg5YwbSgCyV8upypm+cTmRgJHOGzmk645lFxzQf2J8+BHuV5gM79BEIbvzZOga/praV\nrZvIN+z8FP3jH0hAABG/u77e17A77Dz545OcNJ9kwegFTaOdeTINNr0GB1aC+ECPG7Ql8lq011tZ\ns6FZm9bYzWZKlq8gfNw4fKPqF1BKKZ7f+jzrjq9jer/p9I3r28gq3UzOTq2jJ3Od1rkz+AGtTRme\noLeyZkezDs6SL79CWSxE3Tyh3td49adXWXpwKXf2uJOJ3SbWfoIn4nBova47FkHGKm31rZFPaS4E\ngV7fcnErFVY7h/PMZBdZGnytZhucyuGg6OOPCezRg6CePet1jcV7F7No7yJu6nQT9/e5v5EVugGH\nA/Yv09Ls8g5oaXZDZ8Bl9xlDInXE4VAcyS9n29FCfjiUx4aMPCqq7bWfWAeabXCW/3sz1qNHSZgz\nu17nf33sa15Je4UxqWN4fMDj3rV8gsMBB9doyQNnfobYzvD7+dD1OmOicy1YrDa2Hi1k29FCfs4u\nYXd2MWWVWl5OXHgAN/RN4rJ2MSRHB9NjTsPKarbBWfTRR/jExBA29uINvPYX7Gfmppn0admHWUNm\neU/PbHkB7PpYM2YuOgpRqfD797V5ld7yGdyIUorjBRa2HStkY0Yee3NKyCq04FDgaxI6x4dxTc8E\n+iRH0r9NNKkxwY36I90sg7PqyFHMGzcSc8/dF72EfI45h/vX309UYBSvDnsVfx8Pv9M4HHBis+bT\ns/8rsFsheaC22E+Xa8FY/xPQAnFXVjH/PlxAZq6ZI3lmjuSVU1al3RVbhQfSt3UU1/VKoF9qNAPa\nRBPo59ofNLcHp4gkA39HM/BSwDyl1Ovu9K0tmD8f8fcn+pZbLuq8osoi7l53NxX2ChZftZiYoBhX\nyGscqitg2/uw/X3NmDkgHPreDn3vgLiueqvTlcpqO4fOmDmcZ+ZEoYUjeWZ2Z5dwNL8cgPiIQNrF\nhnJ9n0S6JoTTOzmSzq3C3N500ePOaQMeVkr95HTzSxORdWiufi73ra3OyaFk+XKiJkzAN6buwWWp\ntjBl/RRyzDnMGzWPjlEdG1ta41B0DPYs0VbeKj0JqZdri/x0GQ/+IXqr04W8sioOnC4lM9fM+vRc\nNh8pwO74jzdcYmQQHeJCufeKdozu3oqIIM+oTbg9OJ32l6ecr8tEJB1IxE2+tQULFoIIMZP+VOdz\nqu3VTNswjX0F+3h12KueN5Zpt8G+ZVrV9fiP2rY2Q+F370Gby/XV5mYcDsWu7GI2ZOSRmVvGgVNl\nHHHeEQHatgjhz5e3oVdSJB1ahpIcHezy6ml90bXNKSKpQB9gK27wra0+eZLiJUuIvP46/OLj63RO\nlb2KaRum8WPOjzw96GlGpIy4qDJdyqk9WgL6nk+1hWNj2ms5rz1vgshkvdW5HJvdQfqpMnYcL2T7\nsUL255RyqqSSKpsDk0DrmBDatwxlwoBkeiRG0qZFCHHhAV7Ts65bcIpIKPA58KBSqrTmH8xVvrV5\nb78NQIvJk+t0fJm1jIc2PMTWU1t5cuCTnrF2ptWiJQxseReObwJE8+a56lnodHWTdbFzOBQHTpeR\ndryQw3nlpJ8qZU92yS9jiomRQfRKjmBUt1Z0SwhnWMeWRAR7RvW0vugSnCLihxaYHyulzjr6udS3\nturIUUqWfUn0rbfU6a6ZVZbFfevv40TpCWYNmcW17a5tTDkXT34mbH5Da09WWyAsAUbNgl4TIKSF\nvtoaEaUUZ0qrOFlcwfGCcrYfK2JfTglH8soxO3tOQ/x9aN8ylJv6J9MnJZL+qdEkRDa9BY706K0V\nNP+hdKXU3Bq7zvrWzsYFvrW5c1/BFBhIzF21W4YcKznGpK8nUWmvZN6oefr5zJaegiPfwe5/wtGN\n4BMAPW/UxiVTh3j9MEhhuZUDp0spMFvJKrKwMSOPn0+WYLH+J8MmLMCXnskR/P6SRHolRXJp22gS\nI4O8pmraEPS4cw4GbgV+FpFdzm1/wYW+teVbt2H+Zj2xDz5Yaw/twaKD3LPuHuzKzuIxi+kQ1aGx\nZNSNyhI4uFZbVv3EZm1bRIrW49r3dght6V49jcSZ0ko2Hy5gd3Yx6adKOZxXTl5Z1a+O6RIfzo39\nkmnXMpSkqCASI4NoFxuKj6npB+K50KO31tlQOieN7lurHA5y58zBNz6e6NsvnJi+MWsjM76fQYhf\nCAtGLaB9lJumR9msWsL5jkVwbJO29F1UGxj5V2h/FcR195q2pFKKw3nl7DxRRFahhSP55ezPKf2l\nxzTIz4fO8WEM6xhLx7gwusSHExceQExoANEhHp7Q4WaafIZQybJlVO7fT8LLL2MKDDznMUopFu1b\nxGtpr9E5ujNvjHjD9R6z1RXa8MeeT+HEVrBVaHfIwQ9Ah6u0LB4PDUilFIXlVk6VVLL1aCFpxwux\n2RXFFdVknC6jpEIz0DAJJEQG0TU+nBv7JzOkfQs6twrD18czP5en0aSD015aSu4rcwnq04fwa64+\n5zHl1eU8s/kZVh9dzejU0Tw3+DmCfF3UuVCSraXQZX4DJ7ZoHTvR7bTqavuR0G6ER+W4ajMuzGTm\nlnOyuIKM06UcOF1GZq75V+3ClOhggv19CA3w5eqe8fRKiqBv62haxwTjZwRivWnSwZn35pvYi4qI\ne3/eOTuPLnHPAAAL4klEQVQQMgozmLZhGtnmbB645AEmdZ/U+B0NpTlasnn6SjjlbGLHdoE+t2i5\nralDdF2aQClFblkVh86YySmu4HC++ZfXWYUWymsEYXSIv3YX7JdM65hg4sID6ZEYQXK0YRjtCpps\ncJo3/UjRhx8RdfMEgrp1+6/9646v44lNTxDqF8rC0QsbL+vHVqWtD3J0I2Tv0F4rOyQN0NqQXa+H\nmHaNU9ZFUGCu4qcTxZRVVpNvriLjtJljBeUcyy+noNz6y3H+PibaxoaQFBXMwLYxdEsIp3OrcBIi\nA4kO8W8WvaSeQpMMzurcXHIefZSA9u1p+cgjv9pXUlXCy9tf5qvDX9GzRU9eHf4qLYMb0ANaXaFZ\ne5zZpy3Wc/g7qC4Hky/EddMmLve9HaLbNugz1YbN7uBMWRVZhRYOOHtDzVU2iixWThZVkJlnRtVI\n2WgZFkDb2BCu7BJH5/gwOrUKIzkqmFYRgUZV1ENocsGpbDZypj+Co7ycxA8WYwr6T/tx08lNPPXj\nUxRUFnBnjzu5p9c99ZvyZc6Dg6shfQUc/R5sldr28ETodRN0GK3ltDZyonmxxUpWYQW5ZZXklFRy\noqCcg2fMZOaaOV1a+atk7vBAXyKC/YgI8qN1TAjjeyUwuH0MMSEBRAb7ERls9Ix6Ok0uOPPeegvL\ntm3Ev/giAe21oZCiyiLmps3ly8wvaR/ZnjdGvkHXmDpOm3I4oCBTWy3r+GZt7LHoqLYvKlW7K7Yd\nBq16aMHZgGqf1eYgu8jCsYJy9p4s5VCumdMlFeSbrRSWW3/pBT2Lv6+J9rGh9E+NIikqmETn2GDH\nuDCvyiE1ODdNKjhL13xNwbvvEfGH3xP5u+upsFXw6YFPmb93PuXWciZ1n8S9ve8lwOcC7uRKaUnk\nRzdq1dSjP0BFobYvOEZbW7LfHdDmCojvVadgrLY7KLZUU1JhJd9sJavQwsEzZeSbrRRbrOSZqzhd\nUkW++deD8snRWrB1T4wgMsiPlOhgUpwdMfERgcSGBmBqpgP0zYEmE5zlW7eR88gjBPXuTdzMmaw6\nsoq5aXPJteQyOGEwD/d7+L+zfZSC4uPaWpKn92iTkk/tAUu+tj88ETqOgdaDtOTyFh3/KxiVUpRb\n7ViqbJwotJB+uozc0kpyS6vIKrJwotBCTnEFjt+k6Af4mmgZHkB4oB+xYQF0i48gPjKQpKhgUmOC\n6dQqjLBA707PM2gYTSI4LWlpZE+ejF9KCrsfuZqH1v6RzOJMusZ0Zfbls7XcWIcD8g9B1jatA+f0\nHs1xrrJEu0hQNESmQMfRkHgJlUlDyPNPpsLmIL+sipwTleTvP0JhuZV8cxX5Zit5ZVVkF1p+sbI4\ni0kgOiSApKggLkmJ4vd9EokNCyA8yI8WoQHERwTSOiak2aalGdQNUeqiZl15FP369VMb33yTrHsn\nUxEdzAv/688Bnzw6Rrbn9qQrubraF1PuXshNR53Zj1RrKWR23xDKorpwKqANmaTwk09PjpOAAgrK\nreQUV/5XFfMsgX4mYkICaBHqT4tQLQDjI4MIDfAlPiLQmY4WaASeASKSppTqV+/zvTk4OyW3Ul+E\nRZMTDc/cbCIxyJ+7is2MKD3zS/JusYRzhGT22pP52Z7CTkd7DqsEFNpwQXigLy3CAgj290EQokL8\nSYgIJCkqiJZhgQT5+/xyt4sNCyAkoElUNgzcQLMOzu6BQeqpK1LZObqKcZVWgipbclQlclzFke+f\nTFF0L0yRSYQGakMH0SH+JEQGER3sT5C/ieToYGJDjV5NA9fQ0OD06ttAVWQgcdOfZnxMRwJiWhMR\n7M/QQF+ig/2N5GoDr8ergzMiqR1Dr7pJbxkGBi7BuL0YGHgoHhecIjJGRDJEJNPpX2tg0CzxqOAU\nER/gLWAs0BW4WUSatz25QbPFo4ITGABkKqWOKKWswD/RzKYNDJodntYhlAhk1XifDVxa84CaptJA\nlYjsdZO2C9ECyDc0AJ6hwxM0AHRqyMmeFpy1UtNUWkR2NGQcqbHwBB2eoMFTdHiChrM6GnK+p1Vr\nTwI11xFIcm4zMGh2eFpwbgc6iEgbEfEHJqCZTRsYNDs8qlqrlLKJyH3A14APsFApte8Cp8xzj7Ja\n8QQdnqABPEOHJ2iABurw6txaA4OmjKdVaw0MDJwYwWlg4KF4bXDqkeYnIski8p2I7BeRfSLygHN7\ntIisE5FDzucoN2jxEZGdIrJSRw2RIrJURA6ISLqIDHK3DhF5yPm/2Csi/xCRQHdoEJGFIpJbc5z9\nQuWKyOPO72qGiIyuSxleGZw6pvnZgIeVUl2BgcAUZ7mPAeuVUh2A9c73ruYBIL3Gez00vA6sUUp1\nBno59bhNh4gkAlOBfkqp7midiBPcpGExMOY3285ZrvM7MgHo5jznbed3+MIopbzuAQwCvq7x/nHg\ncR10fAVcBWQA8c5t8UCGi8tNcv7zRwArndvcrSECOIqzU7HGdrfp4D8ZZdFoIw8rgVHu0gCkAntr\n++y//X6ijUYMqu36Xnnn5NxpfonuFCAiqUAfYCsQp5Q65dx1GnDxEmW8BswAHDW2uVtDGyAPWOSs\nXs8XkRB36lBKnQT+BpwATgElSqm17tTwG85Xbr2+r94anLoiIqHA58CDSqnSmvuU9tPosvEpEbkG\nyFVKpZ3vGFdrcOILXAK8o5TqA5Tzm+qjG/4WUWgTI9oACUCIiNziTg3nozHK9dbg1C3NT0T80ALz\nY6XUF87NZ0Qk3rk/Hsh1oYTBwLUicgxt1s4IEfnIzRpA+/XPVkptdb5fihas7tRxJXBUKZWnlKoG\nvgAuc7OGmpyv3Hp9X701OHVJ8xPNCWwBkK6Umltj13Lg7LLZE9Haoi5BKfW4UipJKZWK9rm/VUrd\n4k4NTh2ngSwROTvzYiSw3806TgADRSTY+b8ZidYp5da/RQ3OV+5yYIKIBIhIG6ADsK3Wq7my08CV\nD2AccBA4DDzhpjKHoFVV9gC7nI9xQAxaB80h4Bsg2k16hvGfDiG3awB6Azucf48vgSh36wCeAQ4A\ne4EPgQB3aAD+gdbOrUarRUy6ULnAE87vagYwti5lGOl7BgYeirdWaw0MmjxGcBoYeChGcBoYeChG\ncBoYeChGcBoYeChGcHoJImIXkV3O2RcrRCSyka//tIhMr8NxIiLfiki4iMSKyCanputrHPOViCTU\neP83ERnRmHqbA0Zweg8VSqneSpt9UQhM0UnHOGC30tIWbwbeRfMbfhBARMYDO5VSOTXOeQP3zJJp\nUhjB6Z1spkbitIg8IiLbRWSPiDxTY/uXIpLmnO94V43tY0TkJxHZLSLra1y3q4hsEJEjIjL1PGX/\nL//JfKkGgtEG/u0i4osWpC/VPEEpdRyIEZFWDfjMzQ93ZLIYj0bJSDE7n32Az4Axzvej0IykBO3H\ndiUw1Lkv2vkchJZBEwPEos2QaPObY54G/o0WaC2AAsDvHDqOA2HO1xHAKrQsoZFocytvP4/+94E/\n6P139KaHR7nvGVyQIBHZhXbHTAfWObePcj52Ot+HouVufg9MFZHfObcnO7fHAt8rpY4CKKUKa5Sx\nSilVheakn4s25Sn7NzqilVJlznNLgKvhlxkijwG/E5H30VL5XlFKbXael4s2c8SgjhjVWu+hQinV\nG2iNdpc82+YU4EWltUd7K6XaK6UWiMgwtFkbg5RSvdCCN7CWMqpqvLZzbutUm4ic63vzJDALrR26\nCS3x++ka+wOBilrKN6iBEZxehlLKglZ9fNjZxvsa+JNzjikikigiLdGqnEVKKYuIdEazVQHYAgx1\nzo5ARKIvUkIG0LbmBhHpACQppTagtUEdaBMEgmoc1hGtam1QR4zg9EKUUjvRZoLcrLSZ/58Am0Xk\nZ7R5lWHAGsBXRNKB2WhBiVIqD20hqC9EZDfw6UUWvwptNkxNZqHNugBttsa9aNP6Xodf5sC2R2ub\nGtQRY1aKwUXhnET8d6XUVRdxzu+AS5RST7pOWdPDuHMaXBRK88h5X0TCL+I0X+AVF0lqshh3TgMD\nD8W4cxoYeChGcBoYeChGcBoYeChGcBoYeChGcBoYeCj/Dy3kdcJmNTkuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x121890828>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "reaches = np.linspace(0.001, 1, 100)\n",
    "min_lifts = [0.01, 0.02, 0.05, 0.1]\n",
    "conv_ctrl = 20000\n",
    "ctrl_scale_factor = 1\n",
    "\n",
    "fig = plt.figure(figsize=figsize)\n",
    "ax = fig.gca()\n",
    "\n",
    "plot_power_by_reach(ax, reaches, min_lifts, conv_ctrl, ctrl_scale_factor, alpha, n_iter)\n",
    "\n",
    "plt.savefig('img/single_cell_power_vs_reach.pdf', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Power vs control/test split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOYAAACrCAYAAAByrOAJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXt4VNW1wH9rJjN5P0kIIRCe4SEoIBEVRB5iFbWi1qpt\nre/6rNa+bq3eVq/aXnqr1l5bW21V0CpqW1T0WkVBREUUKC8FeUMCBBJCyIs8JjPr/rHPhJDnJDNJ\nZsL5fd98c+acffZec85ZZ++99tpri6piY2MTXjh6WgAbG5vm2IppYxOG2IppYxOG2IppYxOG2Ipp\nYxOG2IppYxOGdJliisizIlIkIl802pcmIu+JyDbrO7XRsZ+LyHYR2SIi53WVXDY2kUBX1pjzgPOb\n7LsHWKKqucAS6zcichJwFTDGOudJEXF2oWw2NmFNlymmqi4HDjfZPQeYb23PBy5ptP9lVa1V1V3A\ndmBSV8lmYxPuRHVzeZmqWmhtHwAyre1sYGWjdHutfc0QkZuBmwHi4+Mnjho1quWSVKkvKsJbXo7W\n1oVAdBsbiygnzoQEojIzEZer1WRr1qw5pKoZnSqi08IFiaqqiHTYH1BVnwaeBsjLy9PVq1c3S+M5\ncIC9d95FTdVR4i+4gLiJebiy+uFISMARHQ0OJzgEETn+xKa/W8Cnyp7y3Wwr3cbOsl3sq9xLYdUB\nPD7PcekcOEh0JxDniic2Kpa4qFhinDG4nG7cThdR4sLljMIpTqIcLpzixOlw4JQonA4nThyIOHCI\n4BAHgiDWN4BDBBDzLVj7zTFpSNXor/mPiWkkKdpon//a0uyY4r9FLdwqNce1IZ02nNOQl/rwWT98\n6rWOKT58+HyKFy8+nw+fevHho95bj0fr8fm81KuXWm8tHm8dtd466ny11NbXUeutocpTRaWnihpv\nTTOxkt1JZCdkMzhpMENThjEidQTJ0Ult39gWXFNV1fwlnxf1ePAdPUp9UTE1mzZR/q9/IU4n/R95\nhMTp01vMUkT2tF1o63S3Yh4UkSxVLRSRLKDI2r8PGNgo3QBrX6co+p/fUrtjBwP+8ASJs2YFIa7B\npz4+K/yMt3e9zbKCZRypPQJA38S+jBg0glOTz2Fw8mD6J/QnPTadjNgMkqOTcYht9O5qauprKKkp\n4XD1YQ4ePUh+RT67ynaxpXQbb5YuwVPyDpTA6LTRnD/kfC4edjHpselBl5t+x+3svf12Dj70MAlT\nprRZc3YG6UondhEZDLylqmOt378FSlR1rojcA6Sp6n+IyBjgJUy/sj/GMJSrqt628m+pxvQcPMj2\nc2aRdvXVZN7zs6Dkr6mv4dUtr7LgqwXsrdxLvCue6QOnMzV7KnmZeWTGZ7afiU2P4fF62HR4E6sO\nrGJp/lI2HtqIy+Fi9pDZ3Dj2RoamDA0q/4ply9h7621kPfwQKZdf3uy4iKxR1bxOZa6qXfIBFgCF\ngAfTZ7wR6INRum3A+xjF9Ke/D9gBbAFmB1LGxIkTtSkHH39cN40arbX5+c2OdYS3d76tM1+ZqWPn\njdVr3r5G3975ttbU1wSVp03PsvPITv31yl/raX87TU+ed7Le99F9Wlpd2un8fD6f7vzG5bpt5jnq\nq6trdhxYrZ3Vn86eGA6fporpra3VLWdO1vzbbu/A5T2eqroqvfeje3XsvLF61ZtX6eeFn3c6L5vw\npKS6RB9d9aiOnz9epy6Yqh/v/bjTeZUtXqybRo7Sio+b5xGMYvaqTlDN+vV4Dx8m5fJvdOr8wzWH\nuWnxTby5401uHXcrL1zwAqf1Oy3EUtr0NGkxafwo70e88vVXyIjL4I4ld/Dqllc7lVfcxIkA1G7b\nFkoRe5di1u7eDUD0iJEdPvdwzWGu/de1bC3dyu9m/I47xt9BlKPHjNY23cCI1BE8P/t5pmRP4aGV\nD3VKOaPS0nD26WMrZlvU7dqNuN24svp16Lya+hruXHonhVWFPHXuU5yTc04XSWgTbsS74nl8xuNM\nGzCNh1c+zLu73+1wHtHDh1O7fXtI5epdirl7N+5BOYizY958D376IBuLNzJ36lwmZk7sIulswhWX\nw8Uj0x5hfN/x/OKTX7DzyM4OnR+dm0vdtu1+I2ZI6H2KOXhIh85Zvnc5b+58k1vG3cKsQcGPedpE\nJjFRMTw67VFio2L58Yc/pqa+ueNCa0QPH26cD/bvD5k8vUYxtb6euoIC3IMHB3xOlaeKh1Y+xPCU\n4dx88s1dJ5xNRJARl8F/n/XfbD+ynT+t/1PA50WPyAWgJoT9zF6jmJ79+8Hj6ZBiPrnuSQ5WHeT+\nM+/H5Qyt54ZNZDI5ezJzhs3hhU0vUFBRENA50cOHA1AXwn5mr1HMOssi6x4yOKD0O4/s5KXNL3FZ\n7mWM7zu+y+SyiTzuOvUuohxRPLb6sYDSO5OSiMrMDKlltvco5q5dAAHVmKrK3M/nEuuK5a5T7+pi\nyWwijb5xfblh7A28n/8+Xxz6ov0TMAag2m12jdmM2t27cSQn40xNbTft+uL1fFr4KbeecitpMWnd\nIJ1NpHH16KtJcCUw/8v57ScG3EOGUJefH7Lye41iGovsoOZTuVrg5S0vk+BK4PIRzR2PbWwAEtzm\n+Xhvz3vsr2zf2upMTsZXWYl625x3ETC9RjE9e/Jx5wxqN11JdQmLdy/m4mEXE+eK6wbJbCKV74z+\nDoLwwqYX2k3rTEoEwFdZGZKye41iesvLcaamtJtu4baFeHwerhx1ZTdIZRPJ9Ivvx7mDzmXRjkXU\neduOguFINBOxvRUVISm7Vyim+nz4qqpwJiS0nU6V17e/zqR+kxiaHNxcPJsTgznD51BeV87yvcvb\nTOevMb1lZSEpt1copq+qClRxJCS2mW5L6RbyK/KZPWR2N0lmE+mcnnU66bHpvLnjzTbTORKtpqxd\nYx7DfzEciW3XmIt3L8YpTmbmzOwOsWx6AVGOKC4YcgHL9y3nSM2RVtM5k6ymbHl5SMrtEcUUkR+K\nyJci8oWILBCRmLaCQbeH1+pwOxNbrzFVlcV7FpPXL88eIrHpEF8f9nXqffVtzjxxRnqNKSLZwF1A\nnppYQE5MsOcWg0EHgt8S1lZTdmvpVvaU7+Frg77WeeFtTkhGpo5kYOJAlu9rvZ/paKgxI1QxLaKA\nWBGJAuKA/bQeDLpd/G8pZ0J8q2mW5i9FELsZa9NhRITJ/Sez6sCqVq2zjvh4EMFXEaFNWVXdBzwC\n5GOCdZWp6mJaDwZ9HCJys4isFpHVxcXFAHgrrBqzjabsx/s/Zmz62JCELrQ58Tgr+yyq66tZW7S2\nxePicOBISorcGtPqO84BhmBCVcaLyNWN01iBjFqcdaqqT6tqnqrmZWSYINftNWXLasv44tAXTMme\nEqq/YXOCManfJKIcUXyy75NW0zgTEyO3xgRmAbtUtVhVPcBCYDJWMGiAJsGg28VXaTVlW7HKflr4\nKT71MaW/rZg2nSPOFceEvhP4ZH/riulISozcGhPThD1DROLEOLaeA2wGFgHXWmmuBd4INENvRSU4\nnUhsbIvHP9n3CYnuRMamjw1OcpsTmin9p7C1dCtFR1uuM5yJSXgjtcZU1c+AfwD/BjZaMjwNzAXO\nFZFtmFp1bqB5+ioqcCQktOjArqqs2LeCM7POtKPe2QTFpH5mAbrW+pnOpER8Iaoxe+RJVdX7gfub\n7K7F1J4dxldV2ao73q7yXRRVFzG5/+TOZG1j08CotFFEO6NZX7ye8wY3X1vZkZgU2Q4GocZbUdmq\nRXZD8QYAJvSd0J0i2fRCXE4XY/qMYX3R+haPOxMT8dmKeQzTlG15DHNj8UYSXAkMTh7cvULZ9ErG\nZYxj0+FN1Hprmx1zJCXiO3oUra8PupxeoZjeykqcrQyVbDy0kTHpY+wl8WxCwri+46j31bO5ZHOz\nY84QTv3qFU+rr7Llpmx1fTVbS7dySvopPSCVTW9kXMY4ANYVrWt2zJEUOn/Z3qGYFRUtjmFuLtmM\nV72cnH5yD0hl0xtJj00nOyGb9cXN+5nOEPrLRrxiqireykoc8c0Vc+OhjQCcnGErpk3oGJcxjg2H\nNjTbf2yGSfAGoMhXzJoaqK9vsSm7oXgD/eP72/6xNiFlZNpIio4WUVZ7fLSCUM4wiXjF9DXMxWxe\nY24q2cSY9DHdLZJNLyc3xSyJsK30+ADP3VpjisgAEfmJiLwhIqtEZLmIPCkiF4r0vKmzYWZJE6ts\ndX01+yr3NVxEG5tQkZtqKeaR4xUzlDVmm54/IvIckA28BfwG41geA4wAzgfuE5F7VLXtSEVdiN+B\nvek45s6ynSjKsJRhPSGWTS8mMy6TRHdisxrTERcHDkdI/GXbc8l7VFVbihH/BbBQRNxATtBSBIF/\nzKhpWJEdR3YAMDxleLfLZNO7ERFyU3KbKaY4HDgSE/GVdXFTtiWlFJFhInKydbxOVUO7lG4H8bUy\nSXr7ke1EOaIYmDSwJ8Sy6eXkpuay/UjzxWqd8fEmamOQdKiPKCL3AvcBPxCR9sNTdwO+Ksv408SJ\nfceRHQxOGozLYS+vZxN6RqSOoNJTSWFV4XH7xe1GPW0Hhw6E9vqYdwF/VFX/ggzjVPVK61jzgZwe\nwN+UdbSgmLZjQQvU10HJNijZDtGJkDwQDu+CmjJIGwoCVB+BpGxIHQyuGFAFbx1ERfe09GFDgwGo\ndBv9E/o37Be3G19dFysmUAK8IyJPqOoiYLGIvIOpaVuP5deNNDRlGynmUc9R9lXu45LhAcfz6l1U\nH4F9a6C6FKqKoWwvHNkDJTvh0BbwdcDJ2hkN6gOfBxL6QeZJRoHTR0DmGOg/AdytB0HrrfhtF9uO\nbGPawGkN+8XtRrtaMVX1RRH5J/ATEbkJ+CWwAHCpaqdjwYtICvBXYCwmts8NwBbgFWAwsBu4QlVL\n28vLV1mJIy4OcTob9u0s2wmcIIaf+jrYtxoKPoPirVC6C/auOl75omIgJcco1IjzoO9JkD4caiuM\n0qYOgdgUOLwTEIhJhvJ9cCTf1KQOJ7ji4fAOKNpklL7Guv2OKMgaDxmjTJ4ZoyHnDJNfLybRnUjf\nuL7sKtt13H6jmJ6g8w9kovQw4FWMIj1k7fsFEMwiDb8H3lHVyy3LbhxwLyau7FwRuQcTV/Zn7WXk\nq6lpFlJk+xFjj+q1QyWqsOtDWDMftr8PtZYVMLE/pAyEyXfC0BmQ2A/i0iEuDQJYnpC+owMvv7II\nDmyA3R+bF8H292Dd38xxhwuGToPsPBiQB4OnmiZxL2NQ0iDyy49fE7NbakwRmQd4MIqzT1W/JyIT\ngL+IyCpVfbCjBYpIMnA2cB0Yyy5QJyJzgOlWsvnAMgJQTK2rQ9zu4/btKttlLLKJvcgiW10KX75u\nFHHfGqgohNg0GHMJ5H4Ncs6E+G5yPRSBxExIPBdyz20k4xE4+CVseRu2LYbtSwAFVxwMmgJDzoZT\nrjAvjF5ATmIOHxR8cNw+cbtCshRfezXmBFUdByAiawFUdS3wdUuROsMQoBh4TkTGAWuAH9CBuLLA\nzQA5OTmWYh5veS2oKGBAwoDIj/FzeBesfgZ2LTcPvK/eGGQGTYHhs2DMpeFVE8WmwOAp5nPer6Cu\nCvI/hS3vmP/w3i9gyYPmZXL6bTBgYk9LHBQ5STkcrjlMRV0FiW4zXCduN1rbfBJ1R2nvyX1HRN4F\nXMBLjQ+oasBR7Foo81TgTlX9TER+T5PlEFRVRaTVuLKY4F3k5eWp1tXhaFJj5pfnk5PUo34PwVG4\nHj56FDYtMv27QVNgyg9g9NdNfy6QZmk44I43L5Dhs8zvkh2w6q+w9m+w8e+QPRHybjQvGHfkLSI8\nKNEslJxfkc+YPsYn29FNxp+fiUgS4FPV0CyVC3uBvVa0PDAR8+7BiiurqoUdiSvrq6tFXMcUU1XJ\nr8jntH6nhUjcbqRwPSyba5qC0Ulw1t0w6RZIyuppyUJDn2Fw/n/DjHth3QKjpG/cDu/8HMZ/C06/\nFdKG9LSUAeN3XskvP6aY4nLj64ZxzKuBl1TV18rxYUCWqn4caIGqekBECkRkpKpuwUTG22R9rsWE\nrQw4rqzW1SHRx8bXSmpKqK6vjqz+5a7l8OH/wO6PjEV0xn1w+i1muzcSnQin3wyTvgd7PoHVz8Gq\nZ+Dzp2HUhXDaTTBkWti3DPzP2J7yPQ37ussq2wdYKyJrMH3BYowT+3BgGnCIDqzK1Yg7gRcti+xO\n4HrM2OirInIjsAe4IpCMtM5znPHHbyWLiKZs1SF4/wFY+4IZ0D/3QTj12l4/1NCACAw+y3zKH4bP\nnzKW5s1vGkvu7N+YsdIwJTYqlsy4TAoqChr2ddc45u9F5A/ATGAKcApQjYmc/l1VzW/r/DbyXQfk\ntXCow3Flta7OePVb+C9STmIYK2ZtJXz6B1jxBHiqYcrdMP0ecLUcSf6EICkLZj0A0+4xwy5LH4Y/\nT4Uzb4fpPw9bJ4acpJwWasyu9/zBcsd7z/qEHU2HS/Ir8nGKk6yEMOyX1R01ho+PHoHKg8aYM/MX\nkDGypyULH1wxpik75jLTmljxBGz8h2nej/+2MYaFETmJOSzNX9rwO1SK2eMTnYOl6XBJQXkB/RP6\nh5fzuqrpRz0+Fv71U+OBc+P7cOXfbKVsjbg0uPh/4YbFkDwAFn0f/nAarH/ZXM8wYVDSIEprSymv\nM04eEu0Gny/o2LKRr5i1tccNl+RX5IeX4afqEDw/B96627jCXf8OXP8vGBiBVuOeIOd0uPE98xJz\nx8Frt8DCm8FT09OSAce6TAXlpgvlfxaDrTUDCS3iEJGADDE9gc/TvCkbNop5cBP8ZYbxY73ocbj2\nTRh0ZthbG8MOEdPsv+UjmPGfsPFVeO5844DRw2QnZgOwv2o/QMOz2OWKaQ2V/EdQpXQhxiprhkvK\nasuoqKvoecOP1wPLfwtPT4P6Wrj+bci73lbIYBGBaT+FK180M2WeOhu++GePipQVb2wZhZXGac2v\nmMFO/Qq0Kfu+FZBroIik+T9BlRwiGht//BbZHq0xq0vhhUuNVXHUheYtnx3Zrmdhx+iL4NaPTP/8\nHzfAa7ca63YPkOROIi4qrmHCtN/ZJdixzECdSa+0vu9otE+BoUGVHgIaK+b+StOcaDxxtVs5uAn+\nfp2ZPnXJn403i03XkDrI9NU//B/TOinbC99aYJwXuhERISs+65hihqgpG5Biqmp4+kmpgtfbYJX1\nX5xuHypRhc//Aov/E2KS4LuvwZCp3SvDiYjTBTPvg/RcU2s+ez7M+YOZvN2N9Evo11wxg3TLC6gp\nay3L/p8i8rT1O1dELgqq5FBgmc39F6OwqpB4VzyJrm58a3o98OYPzDDIsBlw26e2UnY3p1wB337F\nRGv4y0wzg8Xnbf+8EJEVn8WBqgMADZVElxt/LJ4D6gD/ssz7gIeDKjkUWIrpN1EXVhaSFZ/V4pLv\nXUJNObx4Ofx7Pkz9MVy1ABIyuqdsm+PJPRe+v8o4IXz0qOnnVx/plqKz4rM4XHOYmvqa7rPKWgxT\n1f/BTJpGVY9iwjb1KP7QgX4n9sKqwgYrWZdTXgjPXWBm8M95Es75JTgiflg4solJhjl/hIufgD0r\n4PmL4ejhLi/W/8wdqDrQfeOYFnUiEosx+PhnlQQ/GzRYfJZiuo41ZbtFMYu3wDPnmvg6334VJnyn\n68u0CZxTr4GrXoKir2DehVBxoEuL6xdvIjIUVhV2e435APAOMFBEXgSWEA5jm9ZsNHG7Oeo5ypHa\nI11v+Nm7Bp49z4xPXvd/MLzDfvc23cGIr8F3XoXSPeZ+le7usqIa15jdOo6pqouByzBxehYAeaq6\nLKiSQ0Ej44+/892lNeauj2D+180k5hvfhf7ju64sm+AZOh2uXWT6mi9cClUlXVJMZlwmgnR/jSki\nf8Mo5g5VfUtVDwVVaoho6GO6XQ3m6i4bw8z/DF660kShu+Fd44huE/4MyDPdjbJ98PK3usTH1uV0\nkRGb0UQxg3MwCLQp+wyQBTwhIjtF5J8i8oOgSg4FjWpMv69il9SYhevhxW+ayHDXLOo9oT5OFHJO\nh8ueMj7Lr15jYvGGGP9YZoNiBhmQK9Cm7AfArzDxZP+CmeR8WzAFi4hTRNaKyFvW7zQReU9Etlnf\nqe1m4jN9TEd0NIWVhTjFGfrVo4s2w/OXGMeBaxYZ5bSJPMZcChc+BtvehYU3gTe4aVlN8Y9ldndT\ndgnwCcY1bwtwmqqOCqpkE7Jyc6Pf92ACPudijEvthyxpVGMWVhWSGZcZ2pCVZfuMUjrdcM0bphlr\nE7mcdiOc92vY9Aa8fmtInRD6xfUzdg6X5WDQHZ4/wAaMg8FYTHiRsdbwSacQkQHAhZjo7n7mYAI9\nY323u/CINlFMv9k6JNRVwYKrzPd3F5oIbzaRz5l3wDn3m/CZb/80ZNlmxGVQ662lSkwTtltqTFX9\noaqejTEAlWA8gYJxq3gcM9zSOPpewAGfRWS1iKyuKDOrNPitsiFTTK/HzFo4+AVc/mxYB4Sy6QRT\nfwST7zLBtDe8GpIsM2KNx9ehOqMW3TJcIiLfF5FXgLWYmu1ZYHZnCrR8bItUdU1radRUha0GfFbV\nPFXNS0ywfGJdLoqOFpEZH4L+n88Hb3wftr4DF/zWjIfZ9D7Oud8E0n7zbmNHCJKMOKOYxTXFIYn7\nE2hTNgZ4DBilqrNU9b9UdWl7J7XCFOBiEdkNvAzMtIZjDlqBngk44LPlYFCpNXh8noa3VlAs/y1s\neNnMlD/tpuDzswlPnFHwjWcgOgFevCJo7yC/0bH4aDESHd09wyWq+ghQA9xq1Z7jOlugqv5cVQeo\n6mDgKmCpql4NLMIEeoYAAz77+5glPrN4rf+t1Wl2LIVl/w2nXAln/yS4vGzCn6QsMyvlaAm8dIUJ\nK9pJGpqy1Ye6r8a0VpZ+Eehrff4mIncGVXJz5gLnisg2YJb1u20sxTxUb/qaQdWYFQfgnzeZdR4v\n+p0dBuREof8E+OZzcGAjLPxepy218a54YqNiKTpa1K1N2ZuA01X1l6r6S+AM4HtBlQyo6jJVvcja\nLlHVc1Q112outz8twFLMYq9Z37ZvbN/OCgKv327ivl4xP2yDC9t0ESPOg/N/Y9aMWdLhlSUBE8kg\nIzbDqjFd3aaYAjR+lXgJg2lf+HzgcFBcZ3Q4Pa6TzgWfPQU7lsB5D9txXk9UTr/ZLE/xye8hf2Wn\nskiPTae4ujgkK351ZKL0ZyLygIg8AKzEuOn1KKqKuN0UHy0mwZVAbFQnhlb3rIDF98GI2WZJOJsT\nl/N+bZxI3rijU8G9+sb1NTWmq5sUU1Ufwyz8c9j6XK+qjwdVcijwK2Z1cecMP2X7jO9k6mC49M92\nv/JEJzrBTLIu2W7iN3WQ9Nj0Y33MID1/2luGLwa4FbO610bgSVUNrZNhMKgi0abG7LDhR9WszVh3\n1MyrPFFW2LJpm6HT4czvm0Wfcs6Eky8P+NSMuAyq66vxuZzQxTXmfIzD+kaMQ8EjQZUWanyKw9XJ\nGnP1s7BzGXztIbtfaXM8sx6AgWfAors6FO3dXzl4nF0/7eskVb1aVZ8CLgfODqq0UKO+hj5mh2rM\n0t2w+BcwdAbk3dBl4tlEKE4XXP4MiMNEQAxwESN/5VDn9HV5H7NB7cOqCWuhqvhcTup8dYErpt/l\nzuE0MUjtfqVNSyQPgK89CLs+hH8/H9Ap/mewRrpeMceJSLn1qQBO8W+LSHlQJYcCVbxR5i8E3JRd\n9VezpPp5vzYX38amNU69zqxsvfgXUL6/3eR+t7waR33XKqaqOlU1yfokqmpUo+2koEoOBap4nKaZ\nEVCNWbID3r8fhp8LE67uYuFsIh6Hw6zR6a0zzu7tNGmT3ElEO6OpFk+3jWOGJz6lzlpguN0a0+c1\n3j0Ol7nYdhPWJhDShpqYwdvebXeKmIiQHptOldR1z9ol4Yqqj1qHcUhqt8Zc9VcoWGkW+0nqoUWH\nuhFVpfSoh7JqDw6B2nof5dUe9pZWU+3xkhrnBpRqjxenw4HX5+NQRR11Xh9OhxDlEFxOBy6ng7R4\nF1nJsXhVqfF4iXKYfQNS44hxhdfS613C6beY5f7evddMA4xtPepNSnQKR+XQia2YqFLj8BLviifO\nFdd6uqpD8MGvjBV23FXdJ18XUHbUQ0HpUSpq6qnz+thXWs3O4ko2HyjncJUHVaWqrp7SKg+VtV1r\nrxOBYRkJDM9IoKquHqdDGJAaS1ZyLOkJbpJjXSTFuEiOc9EnPpo+CW5czghspDmccNFj8PR040t7\n0e9aTWoUsxCfp3uW4QtPVKl2eEmLaWepzqUPmyk958+NiCasqrKn5Cj7jxi3sD2Hj7K+4Aif7ixh\nT8nRZumjoxyM7JdIdkosIhDvdpIS52ZgWhypcS5UIdrlIDHGRXZKDHHuKA5X1eEQIdbtxOvzmWZY\nQjTRUQ58qni8isfrw+P1UVJZR2FZDVFOMcd9UFJVy47iKr7cV8a2ogoSY1zU+3yszT9CWXXLD6XT\nIZwyIJm8QakMSU9gZL9ExvRPioxaN2scTLoFPvuzsU+0suZpSkwKldTCCa+Y4iE1po3IBfvXwpp5\ncPqt0DfY+GGhx+dTdhRXsr2oksKyGj7fdZgVOw5RXnN8bZcUE8WkIX341qQcBqXFkRLnxuUUslJi\n6ZcUg9PRsRdO/5TA/YqzkmMZm50ccPoaj5eSqjrKqz0csZrTJVW1FByu5rNdJcxfsYc6rxVFXyDB\nHUVmcgznjO7LWcPTGdkvkYyE6O5bHCpQZtwLXy40sYJufL/FtWpSoi3FDJKIVkz1+TgqHlKjW2nz\ne+vNAHFCX5jeftC9rkZV2XWoitW7S1mzp5RtRRVsK6qkopES9k+OYfbYLCbkpDCoTzwikJUcQ05a\nXPg9qK0Q43KSnRJLdivK7/Up+49Us6mwnM2F5ZRVe9heVMkzH+3iqQ93ApAQHcXg9DhOzk5mTP9k\nRmclcnJ2Cu6oHmwKxyTBuQ/BazfD2hdg4rXNkqREp7Bbgo9b2+2KKSIDgecxwbYUeFpVf28tHf8K\nMBjYDVyhqqVtZqZKldSRGtOKYq76iwnWfPlzPeILq6rsO1LNuoIjLNtSzIdbiymuMG/T1DgXo/ol\ncfG4/oyqDkf3AAARS0lEQVQfmMJJ/ZPITIqhT7w7YhSwszgdwsC0OAamxXHemGMB1MqqPXy5r4wt\nByvYU3KU7UWVvL3xAAs+LwAgMSaKGSP7clZuOqcNTiMnLa7DLYWgOeUKWPMcLPkvOGlOs+cqJToF\nTwha5j1RY9YDP1bVf4tIIrBGRN7DrIuyRFXnisg9mLiyP2szJ1UqqW1ZMSuLYOmvYPgsE+y3m6is\nrWfF9kN8uLWYpV8VUVhmQvInx7o4e0QGZw7tw6QhaQzLiO/1CthRkmNdTB6ezuThx+bVqir7y2rY\nuLeMpV8dZOlXRSxabwb73VEOzs7N4Iq8AUwZnk58dDc8ziLGVvH0dPjoEfja8cvEpsSk4AmBGN2u\nmFaIykJru0JENgPZmOh7061k84FlBKCYtQ4vA6JbMP4sfRjqq7vc4FPv9bF+bxlr80tZtqWYlTtL\nqPcpcW4nU3PTuW36MMYNSGFM/ySiItEi2cOISEOz+Pyx/VBVthysYMPeMjYXlvPm+kLe33yQKIcw\nfmAK00dmcPaIDMb0T+662rT/eLP04so/w8Trj4s5nBqdSn2E1pgNiMhgYALwGR2IKwvcDDAmOgZP\nlDSvMQs3GP/GM26D9Nwukb3g8FEWrd/Piyv3sN+qFYdmxHPjWUOYPrIvEwel9mx/qJciIozql8So\nfsbx7N4LRrNyZwmf7ijh4+2HeGTxVh5ZvJWUOBfTR2Rw/th+nDM6M/TDNDN/CV++boZPrpjfsDuS\nm7IAiEgC8E/gblUtb9ysU1UVkVbjygJPA4yNiVWPk+aKueRB0/afFrolPCtqPLy+dh8fbj3EzuJK\ndh6qAmDysD78/ILRnD4kjb5JMSErzyYwXE4HU3MzmJqbwX8AxRW1rNhhuhIffFXE6+v2k5EYzYUn\nZzE1N52puRmheWEmZhpL/0ePwMEvG4KCp0RHaFMWQERcGKV8UVUXWrsPikiWqhYGHFcWqHdy/Dhm\nweew/T0zp64ND41A2X2oinkrdvP31QVU1XkZkh7PiMwEvjUph/PG9COnTxuODTbdTkZiNHPGZzNn\nfDb1Xh/LtxXz0mcFvLwqn3krdpOeEM0VeQO4dEI2uZmJwRV25h3w+dOwbC5c+QJg+pgR2ZQVUzU+\nA2y2Qpb48ceVnUuAcWUBPFFNaswPfgVx6XBa54P4HSirYclXB1m0bj+f7TqMyylcdEp/rjlzEOMH\npthGmwghyulg5qhMZo7KpMbjZcWOQ7y4Mp8/f7iDJ5ftYGx2Et+alMOc8dkkdMZwFJcGZ9wOH841\n3aesU4h2RuNwRwMdjxl0nOxBnd05pgDfBTaKyDpr370YhXxVRG4E9gBXBJKZx8mxccz8lVZUgodN\n/JYOoKps2FvGn5bt4J0vTVTuoenx/PS8kXxz4gC7mRrhxLicDUpaVFHD/20o5JVVBdz32hf86v82\nM2d8f66dPLih7xowZ9wGK/8EH/4GrnrRlBWbSMQppqp+TOuhL8/pcIZu17HoeMt/a2rLDkS7K66o\n5ZmPd7Fo3T72l9WQGBPFHTOGccn4bIb3TbBrx15I38QYrp8yhOsmD2ZdwREWfJ7Pa2v3seDzAiYP\n68PNZw/l7NwMHIFYdWNT4MzbTQR/q9aMjU0iwJ5Yq0S05w9AdKylPPvXwvb3zWIx7vb7fV8dKGf+\nit0s/Pc+PF4fM0dl8oNZuZw/NovkWFc3SG7T04gIE3JSmZCTys9nj+blVQXMW7GL655bxYDUWC6b\nkM0NZw0hJc7ddkan3wqfPtlQa3oG9uWZu51w21edli3iFTMmxmqyfvQoxCS3uxBQweGjPLp4C6+v\n20+My8GlE7K5ZdowhqTb0ddPZFLj3dw2fRg3njWEtzcWsnDtPp74YDvzVuzmlmnD+PakHFLjW1HQ\n2BQTMHr5I3B4JwlJffiib2HLaQMk4hUzLi7ZRCbY/JZZCCim5T7Cih2H+NOyHXy8/RAup4M7Zgzj\ne1OHtv82tDmhcEc5uGRCNpdMyOarA+XM/ddX/PbdLfzvkm18M28At5w9jIFpLbTI8m6Ajx6D1c+R\nkpzCkdpglo+NcMWsj4LopFT49I9mOfZJNzdLsza/lP9dso0PthTTLymGO2fmctVpAzs0u8LmxGRU\nvyTmXT+Jrw6UM++T3byyqoBXVhVw+/Th3DFj+PHjoUn9YdSFsPZvpJ77IyrqKoIqO6IVs6CvA9/w\n/rD0j8a5OOHYokJF5TX87J8b+GBLMcmxLn4+exTXTh4cGXP/bMKKUf2SmPuNU/jBrFx+86+v+P2S\nbfzfxkLuu2A000dmHDMQnnYjbF5EyuE9QZcZ0YrpUx9ph3ZAfY2Jnm3x3qaD3PPPDVTV1XPP7FF8\n94xB3ePgbNOryUqO5fGrJnDx+P48+OYmrp+3ipmj+vLQJWPNFLch0yBtKCn717WfWTtEvDNn6sHN\nJmp231GUVtVx98tr+d7zq+mbFMOb3z+LW6cNs5XSJqTMHJXJ4h9O474LRvPpjhLOfexDFv57r5ks\nMeoiUg5sCrqMyFfM0nwYfRHvfHGAc3/3IW9tKOTuWbm8cceU4F2ubGxawR3l4HtnD2XxD8/m5Oxk\nfvTqen6+cAOe4eeTUh+BE6VDTarXy7OHTuLBRWs4OTuZF248ndFZPR/y1ubEYGBaHC/edDqPvbeV\nJ5ftIP9QCn91JrIw9hRG8EWn841oxRytbtLpz4MravjmxAH86tKT7alWNt1OlNPBf5w/imEZCdyz\ncAOfJZ7KtJ0fBZdniGTrERx1VbxWNZ4fnTuCO2cOt93nbHqUb0wcgAgs+MfJTHe/H1ReEV29KDBo\n8hXcdU6urZQ2YcFlpw4g/ZTz+Ju3427fjYloxaxLP4k5sy/oaTFsbI7jPy+bxLpT7g8qj4hWzGi3\nOyICONucWMS6nTzyzXFB5RHRimlj01uxFdPGJgwJO8UUkfNFZIuIbLfiy9rYnHCElWKKiBP4IzAb\nOAn4loic1LNS2dh0P2GlmMAkYLuq7lTVOuBlTCBoG5sTinBzMMgGChr93guc3jhB44DPQK2IdN7v\nKXSkA4dsGYDwkCMcZAAY2dkTw00x26VxwGcRWa2qeT0sUljIEQ4yhIsc4SCDX47OnhtuTdl9wMBG\nvwdY+2xsTijCTTFXAbkiMkRE3MBVmEDQNjYnFGHVlFXVehH5PvAu4ASeVdUv2zjl6e6RrF3CQY5w\nkAHCQ45wkAGCkENUW1y7x8bGpgcJt6asjY0NtmLa2IQlEaGY7bnpieF/reMbROTUHpDhO1bZG0Vk\nhYgEN72gk3I0SneaiNSLyOU9IYOITBeRdSLypYh8GGoZApFDRJJF5E0RWW/JcX0XyPCsiBS1Np7e\n6WdTVcP6gzEC7QCGAm5gPXBSkzQXAP/CLFZ0BvBZD8gwGUi1tmeHWoZA5WiUbinwNnB5D1yLFGAT\nkGP97ttDz8W9wG+s7QzgMOAOsRxnA6cCX7RyvFPPZiTUmIG46c0BnlfDSiDFWvy222RQ1RWqWmr9\nXIkZgw01gbos3olZGDi4Jac6L8O3gYWqmg+gqj0lhwKJ1pqsCRjFrA+lEKq63Mq3NTr1bEaCYrbk\nppfdiTRdLUNjbsS8JUNNu3KISDZwKfCnLig/IBmAEUCqiCwTkTUick0PyfEHYDSwH9gI/EBVfV0g\nS1t06tkMq3HM3oCIzMAo5lk9JMLjwM9U1deDcZCigImY9U5jgU9FZKWqbu1mOc4D1gEzgWHAeyLy\nkaqWd7McHSYSFDMQN72uduULKH8ROQX4KzBbVUtCWH5H5MgDXraUMh24QETqVfX1bpRhL1CiqlVA\nlYgsB8YBoVTMQOS4HpirprO3XUR2AaOAz0MoR3t07tkMdae8Czr5UcBOYAjHOvljmqS5kOM72J/3\ngAw5wHZgck9eiybp5xF6408g12I0sMRKGwd8AYztATn+BDxgbWdaCpHeBfdlMK0bfzr1bIZ9jamt\nuOmJyK3W8T9jrI8XYBTjKOZN2d0y/BLoAzxp1Vb1GuIZDgHK0aUEIoOqbhaRd4ANgA/4q6qGdHpe\ngNfiIWCeiGzEKMbPVDWk08FEZAEwHUgXkb3A/YCrkQydejZtlzwbmzAkEqyyNjYnHLZi2tiEIbZi\n2tiEIbZi2tiEIbZi2tiEIWGvmCLitWYp+D+DQ5DnJY3j1YrIgyIyK9h82ylzgTW74IchzvcBEfmJ\ntd3l/6MVGW7tIre79sodLyIdXlXKmvnyVlfIFCrCfhwTqFbV8a0dFJEoVe2oY/IlwFuYGRCo6i+D\nkK9dRKQfcJqqDu/Kcrr6f7RRbpePn7bCeIyn09tND3TyuQgIEXGqqrcr8m4g1F4QXeBVUdnCvusw\nQbqWAh9iZg4sAf6NcVae0yjtNZiB7vXAC5jpWYeBXRg/ymE08pDB+HeutfJ5Foi29u8G/qtRGaNa\nkCsGeM46vhaYYe3fAFRb5U1tcs7Xgc+s9O8Dmdb+B4CfNEr3BTDY2r4P4972MbDAn67J/5hoXZs1\nmEH4LGv/MuA3GLe0rX55MIP0j1jlbADubCufJv+hQdbW8m/hnJ9Z12k9xm0OjKKttMp/jWPT6Jrl\nifH2yQeKret6pSXHC8An1nVp7X5MB95qQaY44FXMC/s1677k+Z9D4FFL3rNo+zlJt7bzgGWNrtEL\nwKfANuB7bT73Pa14ASim17rw64DXGinmXiDN+h0FJFnb6RgvCwHGWDfSf6H86efRyFXN/9u6kQXA\nCGv/88DdjS64/2G9HePN0lTWH2M8UMD4ZOZbeQ6mdZetVI45etwEPNqWYmIUZaP1ECVZ//UnTf6H\nC1gBZFj7r2wk17JGZVwAvG9t3wb8A4jyX6u28mlHMZvl3yT9bCvfuCb3ZQMwzdp+EHi8HZmvA/7Q\nRI41QGw792M6LSvmT4CnrO2xmClifsVU4Apru73npDXFXI9x6k+3zu/f2nMfyU3Z91TVPw9OgF+L\nyNkYF7BsjG/kTODvarlhNUrfGiOBXXpsFsR84A7MjA2Ahdb3GuCyFs4/C3jCKusrEdmDmQLV1myG\nAcAr1hw9N6Ymb4upmBfUUQARaSm850jMg/We5R7oBAobHW/8PwZb27OAP6vV/FPVwyIytp18WqOl\n/BszC3jO/x+sspKBFFX1RzuYD/y9A3n6WaSq1dZ2a/ejNc4Cfm+l/0JENjQ65sXMcYX2n5PWeMOS\nrVpEPsDMKW1xckEkKGZrVDXa/g5mhvpEVfWIyG7MWy3U1FrfXkJ37Z4AHlPVRSIyHfNmBfO2bmyc\n68j/EeBLVT2zleOB/o/28mmNrrhOgeZZ1caxYKjRwPqVje9b03vW1P+1VX/YsLfKBkgyUGQp5Qxg\nkLV/KfBNEekDICJp1v4KILGFfLYAg0XEb6T5LqZ/FSgfYV4SiMgIzIyTLQHI7p8GdG2j/bsxISuw\n4sQMsfYvBy4RkVgRScT0UVv6HxkicqZ1vktExrQjx3vALSISZZ2T1sl8AuE94HoRifOXpaplQKmI\nTLXSBHLtW7uPfjp6Pz4BrrDSnwSc3Eq6tp6T3ZjuBsA3mpw3R0RirOdxOibAeYv0FsV8EcizZhFc\nA3wFoCZY9K+AD0VkPfCYlf5l4KcislZEhvkzUdUajPf/3628fEBHLI5PAg7r3FeA61S1tp1zHrDK\nW8PxC+H8E0gTkS+B72PNZVTVf1t5r8dMJ2p2c9WE2rgc+I31v9dhjF5t8VdMH2yDdc63O5lPu6jq\nOxjj3WoRWYfp24F5Mf3WakKOx/Qz2+ID4CRrGO3KFo539H48iXkRbQIeBr4EylqQv63n5L+A34tZ\nt6RpDbvBknkl8JCq7m9NEHt2iY2NhZj1WV2qWmO9sN8HRlovqGDzfgAzwvBIIOkjuY9pYxNq4oAP\nRMSF6V/fHgql7Ax2jWljE4b0lj6mjU2vwlZMG5swxFZMG5swxFZMG5swxFZMG5sw5P8BzljbYrmQ\nZyUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x127e76fd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "reach = 1\n",
    "min_lifts = [0.01, 0.02, 0.05, 0.1]\n",
    "conv_ctrl = 20000\n",
    "ctrl_splits = np.linspace(0.001, 0.999, 100)\n",
    "\n",
    "fig = plt.figure(figsize=figsize)\n",
    "ax = fig.gca()\n",
    "\n",
    "plot_power_by_ctrl_scale_factor(ax, reach, min_lifts, conv_ctrl, ctrl_splits, alpha, n_iter)\n",
    "\n",
    "plt.savefig('img/single_cell_power_vs_split.pdf', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Minimum sample size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10%: 1352, 54068\n",
      "5%: 5107, 204271\n",
      "2%: 31190, 1247589\n",
      "1%: 124459, 4978355\n"
     ]
    }
   ],
   "source": [
    "reach = 1\n",
    "ctrl_scale_factor = 1\n",
    "\n",
    "for min_lift in [0.1, 0.05, 0.02, 0.01]:\n",
    "    min_conv_ctrl = get_min_conv_ctrl(reach, min_lift, ctrl_scale_factor)\n",
    "    N = 2 * min_conv_ctrl / 0.05\n",
    "    print('{:.0f}%: {:.0f}, {:.0f}'.format(100 * min_lift, min_conv_ctrl, N))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multi-cell lift study\n",
    "\n",
    "### Power visualisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOgAAACrCAYAAABsZdC6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFdxJREFUeJztnXm0HUW1h78fIXgDCQJJiIGAITIJiDFEBRFXosJjjCiD\nusBHFJwQRRE0DAo+5QVcjwfyUDEIhNEnoyLgk4BBhBWGBJJAZI6AF5HABYLEhDDs90f1SQ4355xb\nZ+jT3efsb61e3V3dVbVP5+5U1a5dtWVmOI6TT9bKWgDHcarjCuo4OcYV1HFyjCuo4+QYV1DHyTGu\noI6TY1xBHSfHuII6To5xBXWcHLN21gI0w4gRI2zs2LFZi+E0wsMPh/M222QrR0rMmzfveTMb2Ww5\nqSuopEHAXOBpM9tX0kbAr4GxwBPAwWb2oqSpwEQzOyq27LFjxzJ37tzWC+2kz6RJ4XzrrVlKkRqS\nnmxFOe3o4h4NPFh2Pw24xcy2Am5J7h3HqUCqCippDLAP8Muy5E8AFyXXFwH7V8i3j6Q5kkakKZ+T\nIePHh8OpSdpd3LOA7wDDytJGmdkzyfU/gFHlGSR9EjgG2NvMXuxfoKQvAV8C2HzzzdOQ2WkHZ52V\ntQSFIDUFlbQvsMTM5kmaVOkdMzNJ5evdPgpMBPYws5er5JkBzACYOHHiGmvlXnvtNXp7e1mxYkWz\nPyETenp6GDNmDIMHD85aFCcHpNmC7gpMkbQ30AOsL+lS4FlJo83sGUmjgSVleR4HxgFbEwxLddPb\n28uwYcMYO3Yskpr8Ce3FzOjr66O3t5ctttiibfWeOesRAL61+9Ztq5NDDw3nSy9tX50FJLUxqJkd\nb2ZjzGws8Bngj2Z2KHAdcFjy2mHAb8uyPQkcAFwsaftG6l2xYgXDhw8vnHICSGL48OGFbf3rorc3\nHE5NspgHPQ24QtLhBIU8uPyhmT0k6RDgSkn7mdnj9VZQROUskaXsZ856pL2taD3Mnv7W+8nHZyNH\nm2mLgprZrcCtyXUf8LEK78wEZibX9wHbtUO2bqbUta2UlltFLTF7elcoqbv6dSmVlLP/84HeaRv9\nW88uotCufk6B2WWXgd8ZSDFLzzu4Je14BW11KxDb9Rs6dCivvPLKqvuZM2cyd+5czjnnnJbKU1im\nd2+rWA8dr6DOW8lNt9WJwsegTjYccEA4nJp4C5oSy5cvZ3yZr+kLL7zAlClTMpSoMVKbeunrq/6s\nXqNQB49FB1RQScOTqRGnDoYMGcL8+fNX3ZfGoMuWLePII49knXXWYdKkSRxyyCEZSunknZgu7p2S\nrpS0t4rsAZATrrnmGg488EDOO+88rrvuurbVm6tpk7SYPb3jpmRiFHRrgnP654BHJf2npJzPYueX\n3t5eNttsMwAGDRqUsTRO3hmwi2shutIsYJakycClwJGSFgDTzGxOyjI2Rd48YsaMGUNvby/jx4/n\nzTffzFqcKFLxLvrYGs5kHdf6tYKoMShwKKEFfRb4OsHhfTxwJdC+ZRcFonwOFGDq1KlMnTqVZcuW\ncdRRR3HDDTew3377ZSRdDvje97KWoBDEWHHnAJcA+5tZ+fKDuZLOTUeszmW99dbjwgsvzFoMpyDE\njEFPMrMfliunpIMAzOz01CRzOpu99gqHU5MYBa20qVfnTTh1KLm13i5fHo406KCxbNUurqS9gL2B\nTSWdXfZofeD1tAVz8keqS9E6SKlaSa0x6N8J245MAeaVpf8T+FaaQjmOE6iqoGa2AFgg6TIz8xbT\ncTKgVhf3CjM7GLiv3857IkyP7pi6dE5T5HLsWWLffbOWoBDU6uIenZz9Szqt59hj0y2/Qxzoq1px\nyzaXfh74m5k9CbwNeC9hfOo4TsrETLPcBvRI2hS4ieBRNHOgTJJ6JN0taYGkRZJ+kKRvJGmWpEeT\n84ZJ+lRJHbPdwNChQ99yP3PmTI46KjouVK5pydTNpEnh6EAH91YS40kkM/tXsk3mz8zsx5LmD5gL\nXgU+amavSBoM3C7p98CnCMGTTpM0jTDP+t2Gf8FAtPofv+BdJqdYxLSgkrQLcAhwQ5I24DIMC5Qc\nUgcnh+HBkxwnmpgW9GiC59C1ZrZI0jhgdkzhSWzQecCWwE/N7C5JTQVPKgpZ7qiQa+utUxcxy81u\nI4xDS/eLgW/EFG5mbwDjJW0AXCtph37P6w6eVJToZtV2VFi8eDGnnnoqS5cu5aqrrspQwi6h4Nbc\nAbu4kraWNEPSTZL+WDrqqcTMXiK0unuSBE9Kyq4UPGkYYZF4tbJmmNlEM5s4cmTTEcbbzrhx4zj/\n/POzFiN7Dj44HE5NYrq4VwLnEoLwvhFbsKSRwGtm9pKkIcDuwOmsDp50GpWDJx0HXCPpIDNbFFuf\nUzDevTRrCQpBjIK+bmY/b6Ds0cBFyTh0LeAKM7te0hxSDp7ktIemdvxbsTKce9ZpnUAdSIyC/k7S\nkcC1hKkTAMzshVqZzGwh8L4K6e0NnpTR2KPajgp9fX2ceOKJ3HfffUyfPp3jjy/m2Khpjk8M+Wd+\nMVs5ck6MgpZieR5XlmaEQLtOnQwfPpxzz/WNKJw4Yqy4vudQQfDplRoU1JobY8VdV9JJkmYk91tJ\ncgd6p3HctS+aGE+iC4GVwIeS+6eBH6UmkVMocrulSocQMwZ9l5l9WtJnARK/XN9h3mmOf5uQtQSF\nIEZBVybzmAYg6V2UWXPziJlR1P9Dwj7hXcCeO2UtQSGI6eKeDPwfsJmky4BbgO+kKlUT9PT00NfX\nV8g/dDOjr6+Pnp6erEVJn6XLwuHUJMaKO0vSvcDOhO1Ojjaz51OXrEFKoRWee+65rEVpiJ6eHsaM\nGZO1GOlzyuXh7POgNam1J1H/QUJpBcrmkjY3s3vTE6txBg8ezBZbdN/MkBtqIpk9vVBTLbVa0DOS\ncw9hhckCQgu6I2E7zl3SFc3pOHx6pW5q7Uk02cwmE1rOCckKkp0I7ntPt0tApxj4dEs6xBiJtjGz\n+0s3ZvYA8O70RHIcp0TMNMtCSb8kxAWFsPXJwvREcrqCKR/MWoJCEKOgnwe+yup9cm8DGll+5jir\nmez7nscQM82yAjgzOZycUdhx35KXwnnjDbKVI+fEjEEdp/VMvzIcWVCgvXhdQZ2WUtgWPacM2MWV\n9J5yK67j1E1BWqs8EtOC/iwJ4XCkpLenLpHjOKsYUEHNbDfC1MpmwDxJl0vaPXXJHMeJmmbBzB6V\ndBLBxe9s4H3JmtATzOyaNAV0OpSDPpy1BIUgZgy6I2EudB9gFrCfmd0raRNgDlBRQSVtBlxMCO1g\nwAwz+4mkjYBfA2OBJ4CDzexFSVOBiWbWGSHAUqbwxpgPuTNaDDEt6P8QNq0+wcyWlxLN7O9Jq1qN\n14FvJ8o8jNA9ngVMpZ3RzZzsqGUceipZDrh5htEBCrCRWIyRaB/g8pJySlpL0roAZnZJtUxm9kxp\nSZqZ/RN4ENgUj27W8UQ5zp/5m3A4NYlR0JuBIWX36yZp0UgaS1gFcxcQE91sGiG6WW4XhjtOO4jp\n4vaUxfkkCci7bmwFkoYCVwPfNLOXy/cK6uToZo7TCmJa0GXluytI2glYXuP9VSSRta8GLiuz9nZ1\ndDPHqYcYBf0mIZDRnyXdTrDADmhpTaZhzgceNLP/LntUim4GlaObHQBcLGn7CNkcp6OJWc1yj6Rt\ngW2SpIfN7LWIsncFPgfcL6kUyfYEQthBj27WBIWfYgE4dHLWEhSCKEcF4P2Eecu1gQmSMLOLa2Uw\ns9sJexhVon3RzZx8stOWWUuwmhxvJBbjqHAJ8C5gPqsD+BrBCcFxKrLzUzNg3PDqLzz293DecpP2\nCFRQYlrQicB2VsSdoJ1MmbO4D4BdKinqT28IZ98XtyYxRqIHgHekLYjjOGsS04KOAP4i6W7eGmF7\nSmpSOY4DxCnoKWkL4ThOZWKmWf4k6Z3AVmZ2c+JFNCh90Zz+dMT0Sl7JqeN8jBX3iwTXuo0I1txN\ngXOpMFXiODs/NSPuxcP3SFeQDiHGSPQ1gtPByxAWbwMbpymU01mUrLlvYYd3hsOpSYyCvmpmK0s3\nktYmCebrOA3zwJPhcGoSo6B/knQCMCTZi+hK4HfpiuV0POffFA6nJjFW3GnA4cD9wJeBGwk7LDht\nwo1D3UuMFfdN4LzkcJzOJmfW3Bgr7l+pMOY0s3GpSOQUkoGstzXd/pyqxPrilugBDiJMuTiOkzIx\nXdz+NvKzJM0Dvp+OSE5X8LV9spagEMR0cSeU3a5FaFFj15E6TmV8mVkUMYp2Rtn16ySbTacijbMG\nHWvBnfdYOOdp4XYOieni+t4UTuu5dHY451VBc2LNjeniHlPreb8NwZwuJNr/ljJrblrCdBixVtz3\nE3bjA9gPuBt4NC2hHMcJxCjoGGBCEr4BSacAN5jZobUySboA2BdYYmY7JGkeOMlx6iDGF3cUsLLs\nfiX9wjVUYSawZ7+0aYTASVsBtyT3TgWi4ps46ZNxdPCYFvRi4G5J1yb3+7M6+FFVzOy2JCZLOZ8A\nJiXXFwG30i+ymaR9gJMIYQ49NkuOqWfsuQbfWiNmllOBGCvuqZJ+D+yWJH0+2be2EWICJx1DCJz0\nYoN1OAVgzuuh8+bGotrEOhysC7xsZhdKGilpCzP7azMVNxI4CTx4Uqew4fwkYID75tZkwDGopJMJ\n3dDShNBg4NIG62sqcBJ48KROYfRNcxl909ysxcg9MUaiTwJTgGUQImsTFKkRPHDSALhxyCknpou7\nsrw7Kmm9mIIl/YpgEBohqRc4GQ+c1BE0ZRwqIhl6FcUo6BWSfgFskOzw9wUiFm+b2WerPPLASc4q\n5izu8zWiNYix4v5XshfRy4QQhN83s1mpS+Y4Tm0FlTQIuDlxmHeldFrGY0fslbUIhaCmgprZG5Le\nlPR2M1vaLqG6kW4zDK3caP2sRaifDMaiMWPQVwhRsmeRWHIBzOwbqUnl5JJWGoeG3/0QAH0f2Nb3\nK6pBjIJekxyO0zJG3boACArqVKeqgkra3MyeMrMB/W4dp6uYPb1t3dxajgq/KV1IuroNsnQtRRh/\ndt3cZ06opaAqu/Y9cJ3UmbO4r3KgpS6mloJalWvHcdpELSPReyW9TGhJhyTXJPdmZgW0k+eLbu7a\nPvLV/VIpt220acqlqoKamUfRdlLj9WHrVn3m7n+riVnN4jgtZ+QdDzDyjgeyFqN5Zk9PdVsU3yE+\nA7q5a1ti5B2LAHhu1x0qPnfnhYC3oI6TY7wFbTN5bz19vrNBUjIaeQvq5Jpunxv1FrRN5L3ldFpE\ni90AXUGdVbSze/vQ0Z+q6/1unXpxBU2ZIrScWYw733zb4LbX2TZaOO3iY9AUKYJyZsWo2fMZNXt+\nXXm6cTzqLWgXk6XFdvg9DwPw7OTxdeftpjnSXCmopD2BnwCDgF+a2WkZi9QQRWg5O2E6pRsUNTcK\nmmxQ9lNgd6AXuEfSdWb2l2wlG5giKGSJTlDM/pR3eztNWXOjoMAHgMfMbDGApP8lREPLpYIWQSk7\nURkHotNa1Twp6KbA38rue4EPZiQL4EpYZAYyJhVFgfOkoFGURzcDXpVUhCURI4C8xzrNQsYRHHFG\nvXUW4VtC2OS9afKkoE8Dm5Xdj0nS3oKZzQBmAEiaa2YT2yNe4xRBzixkbKTOInxLCHK2opw8zYPe\nA2wlaQtJ6wCfIURDc5yuJTctqJm9Luko4A+EaZYLzGxRxmI5TqbkRkEBzOxG4MY6shTFQlIEObOQ\nsZE6i/AtoUVyysw37HOcvJKnMajjOP3ItYJK2kjSLEmPJucNq7x3gaQl/adcYvO3Uc49JT0s6TFJ\n08rST5H0tKT5ybF3i+WrWG/Zc0k6O3m+UNKE2Lw5kvMJSfcn368lFtQm5NxW0hxJr0o6tp68a2Bm\nuT2AHwPTkutpwOlV3vsIMAF4oJH87ZCTYPh6nLBL/zrAAmC75NkpwLEpyVa13rJ39gZ+T9jzeGfg\nrti8eZAzefYEMKINf5Mxcm4MvB84tfzftZHvmesWlODqVwredBGwf6WXzOw24IVG87eAmHpWuTKa\n2Uqg5MqYNjH1fgK42AJ3AhtIGp2WzJI2kXS1pPskPSTpA03K2U4GlNPMlpjZPcBr9ebtT94VdJSZ\nPZNc/wMY1eb8raynkivjpmX3X0+6bRe0uCs+UL213onJWxeS1ia0ghea2fsIPZ8Hm5QTQniSmyXN\nS7zN0qKZb1J33synWSTdDLyjwqMTy2/MzCQ1bHJuNn/Kcv4c+CHhj+yHwBnAFxqRswDsDzxoZtcD\nmNm/ACTVzBTBh83saUkbA7MkPZT0rApN5gpqZh+v9kzSs5JGm9kzSVdmSZ3FN5u/lXJWdWU0s2fL\nyjoPuL5ROeupN+KdwRF562U8cGeL5cTMSuclkq4ldCfTUNAol9RW5c17F/c64LDk+jDgt23O38p6\nqroy9htHfRJo5QKAGBfK64B/T6ykOwNLky57Gu6X/wC2L91IGtmsnJLWkzQsKW89YA9a+w3Laeab\n1J83batXkxaz4cAtwKPAzcBGSfomwI1l7/0KeIYwKO8FDq+VP0M59wYeIVjyTixLvwS4H1iY/ION\nbrF8a9QLfAX4SnItwmL5xxM5Jg4kcxOyDE1+4yJgPjClWTkJVtEFybGoFXI2+T3fkfwdvgy8lFyv\n38j3dE8ix8kxee/iOk5X4wrqODnGFdRxcowrqOPkGFdQx8kxrqBtRNIbyWqLRZIWSPq2pLWSZxMl\nnZ1cv03Szcm7n5a0W5JnvqQh2f6Kykj6D0lVnTkqvD9J0vXJ9ZTSyg5JIyXdlfjp7ibphLRkLgI+\nzdJGJL1iZkOT642By4E7zOzkfu/tDPzIEu8lSecCt5vZpZH1iPBv+2ZLf0ALkTSJsNJj337pnwE+\nbmZHJPervllXkuaErh9rTHC/0u9+HNBHmICfRHDx2xh4DFhKmMj/MmGlzl+By5J8xxG8UhYCP0jS\nxgIPAxcTJuvfSfComQPcC1wJDE3efQL4QZJ+P7Btkj4UuJDVThMHJOkVy+n3W2YCB9Yqv9/7k4Dr\nk+upwDkEN8CngOeS33468EZyfVnW/35ZHN7FzRALu+gPIihlKW0JcATwZzMbb2a/IHjeHGdmh0ja\nA9iK4Gs6HthJ0keS7FsBPzOz7YFlwEmE1mgCMBc4pqz655P0nwOlRcXfI7jPvcfMdgT+KGnEAOVU\no1L5A32P+cD3gV8nv/27wPLk+pCYMjqNzJ3lnbrZIznuS+6HEhTzKeBJC+skISxo3g64I1kpsg6h\nFSxxTXKeB5Si6X6c4B8KgJm9KGnfAcqpRqXynTpxBc0QSeMIXbglwLtjswHTk5a1vKyxhFaz/L1Z\nZvbZKuW8mpzfoPbfwUDlVCO2fKcG3sXNiGQVx7nAOZYMxCL5A/AFSSVj06aJwak/dwK7StoyeW89\nSVsPUPYs4GtlMm7YYDmt5DVJHRyOuzauoO1lSGmahbDq5SaCMSUaM7uJYP2dI+l+4CpgWIX3niMY\nX34laSGhW7rtAMX/CNhQ0gOSFgCTGyynlcwAFkq6rI115gafZnGcHOMtqOPkGFdQx8kxrqCOk2Nc\nQR0nx7iCOk6OcQV1nBzjCuo4OcYV1HFyzP8D7LQYZJWXPrgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x12768fe10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "reach = 1\n",
    "lift_a = 0.05\n",
    "min_lift_uplift = 0.05 # NOTE: This is the absolute uplift\n",
    "conv_ctrl_a = 10000\n",
    "split = 0.5\n",
    "ctrl_scale_factor = (1 - split) / split\n",
    "\n",
    "L_a = get_lift_dist(reach, lift_a, conv_ctrl_a, ctrl_scale_factor)\n",
    "L_b = get_lift_dist(reach, lift_a, conv_ctrl_a, ctrl_scale_factor)\n",
    "\n",
    "L_ab_h0 = (L_b - L_a)\n",
    "\n",
    "crit_val = np.percentile(L_ab_h0, 100*(1-alpha))\n",
    "\n",
    "# Calculate the test power under H1\n",
    "# NOTE: WE ARE ASSUMING THE CONTROLS IN A AND B HAVE THE SAME CONVERSION RATE WHICH IS NOT THE CASE\n",
    "# - HOWEVER IT'S THE ONLY ESTIMATE WE HAVE\n",
    "L_b = get_lift_dist(reach, lift_a + min_lift_uplift, conv_ctrl_a, ctrl_scale_factor)\n",
    "\n",
    "L_ab_h1 = (L_b - L_a)\n",
    "\n",
    "power = (100 - st.percentileofscore(L_ab_h1, crit_val)) / 100\n",
    "\n",
    "fig = plt.figure(figsize=figsize)\n",
    "ax = fig.gca()\n",
    "\n",
    "ax.hist(L_ab_h0, bins=100, label='$\\mathrm{H_0}$', alpha=0.5)\n",
    "ax.hist(L_ab_h1, bins=100, label='$\\mathrm{H_1}$', alpha=0.5)\n",
    "plt.axvline(crit_val, color='r', ls='--')\n",
    "ax.set_xlim([-0.1,0.1])\n",
    "plt.draw()\n",
    "\n",
    "labels = [x.get_text() for x in ax.get_xticklabels()]\n",
    "locs = list(ax.get_xticks())\n",
    "\n",
    "labels.append(r'$c$')\n",
    "locs.append(crit_val)\n",
    "ax.set_xticklabels(labels)\n",
    "ax.set_xticks(locs)\n",
    "ax.get_yaxis().set_major_formatter(matplotlib.ticker.FuncFormatter(lambda x, p: thousands_tick_formatter(x, p)))\n",
    "ax.legend()\n",
    "ax.set_xlabel('Difference in lift')\n",
    "ax.set_ylabel('Frequency density')\n",
    "\n",
    "plt.savefig('img/multi_cell_dist.pdf', bbox_inches='tight');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.773379"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "power"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Power vs audience size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOoAAACtCAYAAAC+yeOaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXl4lOXV/z9nJnsIARICISxB1iA7uKEiglrFfcel7lr7\nutRaq2Jrq75at9rq2/5s3YqoWEVr3UUQQRAVBYkECBCWAGHJBmRfZjm/P+4HDJhlksxkljyf68qV\nebb7Ps/MfOdennPOLaqKjY1NaOMItgE2NjYtYwvVxiYMsIVqYxMG2EK1sQkDbKHa2IQBtlBtbMKA\ngAlVRP4lIkUisqbBvh4iskBE8qz/3Rscmykim0Rkg4j8LFB22diEI4FsUV8GTj9s373AQlUdAiy0\nthGREcAM4EjrmmdFxBlA22xswoqACVVVlwB7D9t9LjDbej0bOK/B/jdUtU5VtwKbgKMDZZuNTbgR\n1cH19VLV3dbrPUAv63UG8E2D8wqsfT9BRG4CbgJITEycMHz4cJ8r91ZX49m/H291NVpXD7ZXlk0A\nWVtXW6KqPf1RVkcL9SCqqiLSaqWo6vPA8wATJ07UFStWtHiNZ/9+dv3u91QuXIik9SLhmKOJGzqU\nqLReOJO74khIQGJjEacTHE4QATHXikjjhTa1/zD215Wxcd9GtpVvY0fFDgqr9lBcXYIHzyHnJTjj\nSYxJJCEqkfjoeOKj4oh1xhLrNP9jnDHEOKKJdkTjdDhxiJMoh5MoicIhThziwOEQnDgRcSCIdRuC\niODA7OOg2XLw3uTHPT7f18Hz/YgvpamP54UC43925TZ/ldXRQi0UkXRV3S0i6UCRtX8n0K/BeX2t\nfe3GXVxM/qUzcBUXk/bbu+h26QycXRL9UXSjqCqrS1azcPtClhYsZdP+TQBESRRHpB9BZtexHNe1\nP/2S+pHRJYPeib3pldCLuKi4gNlkEyyu9FtJHS3U94Grgces/+812P+6iPwF6AMMAb5tb2Wqyp6H\n/hd3SQmZr71K/Jgx7S2ySXZV7uKdvHf4YPMH7KraRZQjigm9JnD2oLOZ0GsCWT2yiHHGBKx+m8gm\nYEIVkX8DU4BUESkA/ogR6FwRuR7YBlwCoKprRWQusA5wA7eoqqfRgltBxaefUrFgAT1/c2fARFpY\nVcg/V/+Td/PexaMeJvWZxK3jbmVKvykkxSQFpE6bzoeEc5hbc2NUVWXzqafh6JrEwLlzkSj//ibV\ne+p5ee3LvJjzIi6vi4uGXMR1I68jvUu6X+uxCV9EZKWqTvRHWUGbTAo0NdnZuAoKSH/0Ub+LdHXx\nau5fdj9byrZw6oBT+fWEX9MvqV/LF9rYtJGIFWr5Rx8jMTEknXqK38pUVebkzuGpFU+RmpDKs9Oe\n5cS+J/qtfBubpohIoarbTfm8eXQ56SScXbr4pUy3182DXz/Iu5ve5eR+J/PICY/YY1CbDiMihVr9\n7bd4SkroeuaZfimvzlPHXYvvYnHBYn4x+hf8z9j/wSF2PINNI7hqoXwn1Jb5tdiIFGrl0i+RmBi6\nTDmp3WV5vB5mLp3J4oLF/O6Y3zFj+Aw/WGgTEahC6SbYtgx2rTJ/hWvB6/Z7VREp1LoN64kdMgRH\nXPudCB779jEWbFvA3UfdbYu0M1NVCqV5sHcr5H8Ju76HsgKoKzfH47pB+hiYdDukDoX4bvDgdL9V\nH3FCVVVq12+gy8lT2l3We5ve440Nb3D1iKv5+Yift984m/DA64E9ObBjuWkhd66EwjU/Ho/rBv2P\nhcwTIW04ZE6GlEGtcr9sLREnVHdxMZ69e4kb5ruzfmNs3r+ZR5Y/wsReE7ljwh1+ss4mJKktN8Lc\n9Bls/9q8rq80x+J7QO9RMPV+02Im94PUIcYnvAOJOKHWbdgIQOywYW0uw+V1ce/Se4mPiueJyU8Q\n5Yi4t6nz4vVAUa5pIbd9BVsWwf7t5pg4IWMCjL0c+h4NAyZB1z4BbSl9JeK+gXUb1gMQN7ztQn0p\n5yXW713P0yc/Tc8Ev0Qp2QQLVdN13fy56cru+PbHcWVsVxg4GSZcCz2HG2HGdwuuvU0QcUKtXb+B\nqPR0nMnJbbp+075NPLf6Oc4YeAbT+k/zs3U2AUcVqkrMxM+mhbDuPfMaMWIcdRH0Pw56j4aUweAM\nDwmEh5WtoG7DeuKGDm3z9X9Z+RfinfHMPHqmH62yCSiqUJIHa/8Lq16DsgZd2f7HwfG/gqyzIL57\n8+WEMBElVG99PXVbttJlattawuW7l7N051LunHAn3ePC90PtFNTsgw3zYOsS87jkgDiPOBmO/SX0\nGAj9joGEHsG1009ElFDrt24Fj4fYoUNafa1XvTy14inSE9O5POvyAFhn027KdppWM+9TMxHkdUNC\nimk1T7gDhpwG3SIzOCKihOresweAmIxG0y01y8LtC8ndm8sjJzxCrDPW36bZtJWyAtjwiRlr5n8J\nKKSNgEm3QdbZ0Gd8SMzKBpqIEqqryGR2iUpLa9V1Hq+HZ7OfZWDyQM4c6B//YJs2omoenaz9L+R+\nCCUbzP7UoTDlXhh9CfQ4Irg2BoGgCFVEfg3cgMlVlQNcCyQAbwKZQD5wiarua0257gNCTU1tlT3z\n8uexaf8mnpz8JM4OfpBtY7F3K6x5G3LehuL1ZiIo8wQYfxUMPsV4AHViOlyoIpIB3A6MUNUaKwXL\nDGAEJjn3YyJyLyY59z2tKdtdVIyzRw8kxvfcRKrKS2teYnC3wZyWeVprqrNpL14P5C+Fr/4OmxaY\nff2PgzOfghHnQWLrfnAjmWB1faOAeBFxYVrSXcBMTI4lMMm5F9NqoRa1utu7onAFefvyeHDSg3bo\nWkegapwOsufA+g+huhQSUuHk38GYyyJ2Mqi9dLhQVXWniPwZ2A7UAPNVdb6INJWc+xAaJuDu37//\nIceMUFvnSTQndw7JsclMH+i/SAebRnDVmgmhr/9mfGmjE2H4dBh2BgybDtHxwbYwpAlG17c7ZgmL\ngcB+4C0ROSQBanPJuQ9PwN3wmLuoiNgs38cyuyp3sWjHIq498lo7r24gUDXhYKvnwuo3zbPP1KFw\n9jMw8kKItTNk+Eowur6nAFtVtRhARN4BJtF0cm6fULcbd2kp0a3o+r6T9w6qyqXDLm1NVTYt4fUY\n39rFj8HOFeCMMa3m+KuMQ4LDHmK0lmAIdTtwrIgkYLq+04AVQBWNJ+f2CXfpXvB6fR6jetXLh1s+\n5Nj0Y+0Un/6ifDd8+RdY8w5Ul0ByfzMxNPKikHV2DxeCMUZdLiJvA99jkm2vwnRlu9BIcm5fcbfy\nGerKwpXsrNzJreNubU01NoejalKQZL9u/Gy9buOIkHUWDD8bouzVAfxBUGZ9VfWPmMz5DanDtK5t\nwl1sCbWnb0J9f/P7JEQlMLXf1LZW2bnxemH9B/DV36DgO3DGwpHnG6eEHgODbV3EETGeSe6iYsC3\nFrXWXcv8/PmclnkaCdEJgTYtsnDXQc5bsOwZKNkI3QfCGU8ajyG7exswIkioReBwEJXScrTE8t3L\nqXZXc3rm4Qui2zSJuw5WzDJj0MpC6DUKLpoFI87t8LQknZHIEWpxEVEpKT4tX7FoxyISoxM5qvdR\nHWBZmONxm0crXzxmUpZkngjn/9PM3nYCZ/hQIWKE6vLRK8mrXr4o+IITMk6wl0FsDlXY8DF89qBx\njE8fC2f9FQZNswUaBCJGqO6iYqJ7927xvDUlayipKWFKvymBNyocUTXZ+BY/anINpQyBS16BrHNs\ngQaRiBGqp6SE+JFHtnje4h2LcYqTEzPsxZ1+wq5VMP9+4yif3N94EI29sl15hVwuFwUFBdTW1vrR\n0NAiLi6Ovn37Eh0dHbA6IkeoFRU4krq2eN4XBV8wLm0cybFtS34WkZTvgs8fNs9CE1Jg+p9h/NV+\neQZaUFBAUlISmZmZSAS2yKpKaWkpBQUFDBwYuMdSESFUra9Ha2txdm3ed7S4upiN+zZyx3g7oTYA\nNfth6VPw7fOgXjj+djjxLohr+QfPV2prayNWpAAiQkpKCsXFxQGtp0WhikhfTLzoiUAfjNvfGuAj\n4BNV9QbUQh/wVJqs5o6k5oX61a6vADg+4/iA2xTSeD2w8mVY9AhU74UxM4yjQvfMgFQXqSI9QEfc\nX7NCFZFZQAbwIfA4xlE+DhgKnA78TkTuVdUlgTa0ObzlJqGys2vzLcGyXctIiUthaPe2pxMNe7Yu\nhXn3mnQnA46H0x81SzXYhDQttahPqeqaRvavAd4RkRigfyPHOxRPRQXQfIvqVS/f7PqGEzJO6JwB\n4uW7YN5MWPeumSi6eLZxVojw1i5SaFaojYlURAYBCaqao6r1wKZAGecrXkuozmaEmluay766fUzK\nmNRRZoUGHjesegUWPACeOphynxmL2oHaYUWrJpNE5D5gMOAVkVhVDYm1CD3lLbeoX+/+GoDj0o/r\nEJtCgtwPYcH9sHeL8Sg6+xmzPGCE8Nxzz/HAAw/Qq1cvVJWhQ4fyxBNPBHT2NVg02wcUkdtFpKEj\n5xhVvU5VbwBCZmDjqWh5jPrt7m8Z3G0wKfEpHWVW8KgqgblXw5tXQFQczPg3XP1BRIkUICcnh4ce\neojs7Gx++OEHpk2bxgUXXIBqo8lBwpqWBmulwDwROcfani8i80RkPvBpYE3zHW95811fl8dFdnE2\nR/c+uiPNCg6bFsI/Jhn3v6n3wy+WmNxEETgWXb16NSNHjjy4ffPNN7Nnzx527NgRRKsCQ7NCVdU5\nwNnAaBF5H1gJXABcrKq/bWulItJNRN4WkfUikisix4lIDxFZICJ51n+fF3/xVJSD04kkNB6yllOS\nQ427JrKFWlcJH94Jr11gFkO68XOYfBc4A+ctE2zWrFnDkUce6o0WHx/Pvn2tSgcdFvgy/TkImIvJ\n/HcL8AzQ3pmIZ4B5qjoc04XOxeTxXaiqQ4CF1rZPeCsqcXbp0uTzrG/3fIsgTOw9sZ1mhyi7V8Nz\nk2HFv+C4W+GmxWaV7Ahmx44ddOnSha4Nhjsul4vdu3dzxBGRl0m/peeoLwMHcu/uVNUbRWQc8IKI\nfKeqD7W2QhFJBiYD1wBYM8f1InIubczr66kox9HM+PS7Pd8xrMewyHMbdNfD1383ScQSUuCajyCz\nczhz5OTkHNLtBZg1axZTp06lqqqKs846izPPPJO1a9cyadIkFixYwAMPPPCTa8KFlmZ9x6nqGAAR\nWQWgqquAsy1htYWBQDEwS0TGYLrTvwLanNfXW17R5Pi0zlNHdlE2lw6PsEyDhWvhrWtNCNrws8yM\nbifKLH/4+HT+/Pk8+uijfPzxx2RnZ3PFFVdw0003cd5553HjjTfSrVs3tm3bFrZCbanrO09EPhWR\nz4HXGx5Q1VZlCWxAFDAe+IeqjsNkHzykm6tm2q7JvL6qOlFVJ/bsaZJtG4f8xoW6rnQd9d56JvaK\nkG6vKnz/CrwwDWr3w+VvwYw5nUqkYFrUOXPmMGHCBMaPH8/s2bOZN28eWVlZZGdnM3nyZFwuFykp\nKTgcDtasWcOoUeE7HGjJ4eEeEekKeFW10k91FgAFqrrc2n4bI9Q25/X1VlQQM6BxB6nsomwAxqaN\nbZfRIUFZAbx3K2xZZJ6LXvgSJDXa8Yh45syZ0+SxvLw8hg4dyurVq8nKygIgPz//JysrhBMtjVGv\nBF5vyvHe8lJKV9Uvfa1QVfeIyA4RGaaqGzCZB9dZf23K69tciFt2UTb9k/rTIy7MV57etQpenwH1\nlSYMbeL1diLrJnjppZcAGDt2LGPHmh/oV199NZgmtZuWxqgpwCoRWYkZSxZjnPIHAycBJbRidrYB\ntwFzLF/hLZhlFx20Ma+vt7y80TGqqpJdnM0JGSe0wcQQweuFFS+ZgO7EnnD9Aug1IthW2XQwLXV9\nnxGRvwNTgeOB0Zgwt1zg56q6vS2Vqmo20NigsdV5fdXtxltV1egYtaCigL21exnTM2ScqFpH2U54\n95ew9QuTq+j8f0KX1q1WZxMZtOjrq6oeYIH1F3J4rVjUxoLGs4vDeHya8zZ8dKdxqj/raZhwTUR6\nF9n4RthnePgxaPynY9TsomwSoxMZlBxGPq71VfDBr0yS675HwfnPRZyPrk3rCXuh/hg0/tMWdXXJ\nakaljsIZLgmiq0pgzsWwO9ss7HvCne1KLGYTObQ4bSgiDhFp1YJNHcnBELcuhwq1zlPHpn2bGJka\nJg+4ty+HF6ZC0Tq4dA6cdLctUpuDtChU69HM3R1gS5v4McTtUKFu2rcJt7rJ6pEVDLN8RxWWPAmz\nTgcUrv7QRLvY2DTA1wdxn4nIXSLSz4py6SEiIfFg0lvR+Bh13d51AGSlhLBQXbXwzk0mVefIC+Hm\nZdDPXmbD31x33XWkpaUd4j54zz33MHr0aK666qqD+1577TWefvrpYJjYIr4K9VJM5MwSzPPUlZjF\nh4OOt4kWNbc0l6ToJPp26RsMs1qmqhReORdy5pq40Qte8GuaTpsfueaaa5g3b97B7bKyMr7//ntW\nr15NTEwMOTk51NTUMGvWLG655ZYgWto0Pg2CVDVkc1t4yitABEeXLofszy3NJSslKzRTVRathzcu\nM89JL5oFIy8ItkUdwoMfrGXdrnK/ljmiT1f+eHbzKyRMnjyZ/Pz8g9sOhwOXy4WqUl1dTXR0NH/+\n85+57bbbAprtvj341KKKSIKI/F5Enre2h4jIWYE1zTe8lRU4EhORBu50Lq+Ljfs2hub4dP3H8OIp\nUFcBV7/faUQaSiQlJTF9+nTGjRtHeno6ycnJLF++nPPOOy/YpjWJr9OKszDd3QMp/HYCb2Hy/QYV\nb3UNjvhD49i37N9Cvbc+tManXi988bhZvrDPODOzm5wRbKs6lJZavo7k7rvv5u67zRzpDTfcwEMP\nPcSLL77I/PnzGT16NL///e+DbOGh+DpGHaSqT2CCyFHVaiAk+pTe2lrkMKHm7s0FQmgiyV0H79xg\nRDr2Crh2XqcTaaiyatUqVJVhw4bx1ltvMXfuXDZv3kxeXl6wTTsEX1vUehGJx4oRtaJm6gJmVSvQ\n2hoccXGH7Fu/dz3xUfEMSBoQJKsaUF8Nb1xuQtNOeQCOv8N2BQwh7r//fp5//nlcLhcejwcwY9jq\n6uo2l+nyePkuf6+/TAR8F+oDwDygn4jMwTjoX+NXS9qIt6YWiT9UqHn78hjcbXDwPZLqq+H1S2Db\nMjj3WRh3RXDt6aRcdtllLF68mJKSEvr27cuDDz7I9ddfz7vvvsvEiRPp06cPYMLiRo0axejRoxkz\nxrdADlWlut5DZZ2b+99dw4591Xy/bR/ltW6/3oP4mgNVRFKAYzFd3m9UtcSvlrSBiRMn6tvDhyPO\nKAbMfhkwb9xJb57Eyf1P5sFJDwbPuKpSePNK2P618dcdE2GpYHwkNzf3YPB2pFDn8lBR56am3kNF\nrRu310vh9i3cOb+Y/j0SyErvyilZvThjVPpKVfVLahGfWlQReQ34Aliqquv9UbG/0JpanKk/piEp\nrS1lX90+hnQbEjyjitablrRiD1z4Ioy6KHi22LQZVcXtUarr3VTUuqn3eHF5lDq36SJHORwkxjpJ\njo/DURbH6j+eFrDHgb52fV/CLLv4N2t8ugpYoqrPBMSqVnD4ZNLGfRsBGNI9SELN+wzevtas7XLt\nx9A3QnI1dQI8Xi91Li8ur1JR46Ks1oXHa3qcTocQF+UkNspBj8QYkuOjiHY6Dgpzt0MC+szeV4eH\nRSKyBDgKOBm4GTgSk5+3TVhLZazApCE9y3JJfBPIBPKBS1S1xUzKWnPoZFLePjNbFxShrphlYkjT\njoTL34DkEPWKsgHA61Wq6t1U1XmoqnNTXe9BrZx6DhGS46OJj3ESH+0kIcYZVOcZX7u+C4FE4Gtg\nKXCUqvqcfKwJfoXJFHHAb+5AAu7HRORea7vFvL6mRf1RqBv3bSQ1PrVjcySpmkWBlzwJg0+Fi1+G\n2C4tXmbTsagqNS4z8VNZ66aq3oOqIkBctJPUpBgSY6KIcprW0+EIndl5X7u+q4EJwEigDNgvIl+r\nak1bKrVWMT8TeAS409rdpgTc3tpaHHE/dn3z9uV17PjUXWcyA+bMhXE/h7P+GtHLSIQLbo+Xeo+X\nereXGpeHmnoPNS7Pwa5sXLSTlMQYusRFkRgThTOERNkYvnZ9fw0gIkmYxzKzgN5AbBvrfRoTOtfQ\nk75NCbg1IRGH1aJ6vB62lG3h0mEdNMPqrjMzu3nzYdofTKC3/Yw0KHi9ZpKn1u2lrNpFRa3rYGJo\nESEuykFyfDSJsVF0iTXjy3DCV1/fW0XkTcwk0rnAv4Az2lKh5SNcpKormzrH5wTcqamgilgt6vaK\n7dR56jpmfOquN5nq8+abLPUn/sYWaQfi9SqVtS52l9WQV1jB2l3l5BVVsmNvNdX1HlKTYhmQkkh8\n3X5uv/I8zj35GH52wlG88sI/iHY6wi7MzdeubxzwF2Clqrb3Se7xwDkiMt0qt6v1+Kf1Cbi9Jt3w\ngcmkLfu3ADC42+B2mtgC9VXw5s9h80I440mTeMwmoLg9XmpdHqrrzTPM6gPjSxESY5z0TIohLtpJ\nXLSTmCgHDutHszohlqeeeorx48dTUVHBhAkTmDJlysEwtxtuuIGcnBwGDx7MrFmzDgmHCyV87fr+\n2Von5mZr5mupqv7QlgpVdSYwE0BEpgB3qeqVIvIkrUzArZZQD0wmbSkzQh2YHMCovPpqeO1C2LEc\nzvkbjL+q5WtsDJ/cC3tyfDpVUTxe60/1wG8yCUAXh3lc4hTB0Wc0jjMeb7Kc9PR00tPTARM1k5WV\nxfbt2yM2zO12YA6QZv29JiK3+dmWx4BTRSQPOMXabp6DLarp+m4p20KvhF4kRif62TQLjwvmXmVE\neuFLtkj9iKK4vV7qPB5q3ablrHV5cXsVQYiJchAXbRwMEqKjiHU6iXI4cLQiNiQ/P59Vq1Zx0kkn\nRWyY2w3AMapaBSAij2Me1fytPZWr6mLM7C6qWkprE3Bb7o+OBi3qEckBWhvT44L/3ACbFpgxqR1H\n2nrO+PG31+tVymtdlNW4qHN7qXN7D3Zlo51CYkwU3RKiSYiNOtiNbQ+VlZVceOGFPP3003Tt2jXs\nwtx8FaoAngbbHkIgzO1g1zcuHq962Vq2lQuGBEBA7np46xrY8BGc9og9Jm0DPzoXGAeDapcZY0Y7\nHcRHO0mKjaJrfHRAHAtcLhcXXnghV1xxBRdccOj3o2GY28yZM/n000+59tprycvLY8iQILqhHkZr\nAseXi8h/re3zMG6FwcXrBQRHfByFVYXUuGv836J63PCf641Ip/8Zjr7Rv+VHIKrKhsIKNhZWkuZx\nsaW4kup6D1413dj4GCep1jPMLrFRAfX4UVWuv/56srKyuPPOO39yPBBhboHA18mkv4jIYuDAakvX\nWgsaBxdVEEHi4thcthnw80SSxw3v3gy578Ppj9sibYFNRRV8krOH/67ayZaSKgBePCedFK/SIzGG\npLgoEjrYuWDZsmW8+uqrjBo16uDKbn/605+YPn26X8LcOopmw9xEJA7j1zsYyAFe8sPjGb8xbsgQ\nfd0ZxREff8SbNV/y5IonWXzJYlLiU9pfuKsW3r7OtKTT/ggn/vTXuDNT7/by+fpCvt++n3W7ylm/\np5ySynoAjh7YgwvGZTCuf3c8e3cwYkTkrz7XWDifiHRYmNtsTPqVpRgHhyzgDn9U7Be8Ck7zHHXL\nni0kxyb7x8fXXW+yMmxeaHd3Lbxe5estpSxaX8Se8lq+3lxKaVU9MVEOhvdOYurwNEb37capI3rR\nq+uPvte5+4I+lRERtCTUEao6CkBEXgK+DbxJvnNgfWWJi2Nr2VaOSD6i/eMdrwfeudGItBM/J/V6\nzTjzi43FrMjfxw8F+ymuqCM2ykF6chxHZfZgxtH9OH5wati544UjLQnVdeCFqrpDLkduA8+k/PJ8\nTu53cvvLe/82WPeumd3tZCKtqfcwf90e/vP9TlZt30eFlU5kUM9Ejh+UwtSsXpw2ohdx0WGy6FYE\n0ZJQx4jIgYzJAsRb24JxyQ1uandLqOVSx97avQzo2o5kZqrwyW8hew5MmQmTbvWTkaFLrcvDV5tL\neHfVLr7aXEpJpclX179HAueM6cPYft04cUhPeifHtVCSTaBpacXx0P7pVEViY9lRWQDQPqEuegS+\nexEm3Q4ntRhdF7bsKatlSV4x89cWsmxTCTUuD90Sopk2vBcDUhIY3787kwalhFQspk2Yr4+qXu/B\nbi9AZtfMthW0/DkT9D3u53DqQxETBaOq5JdW8+WmElbk72Xltn0U7DMhxH2S47hoQl+mDk/j+MGp\nxETZ48xQJqyFiteLxMezvWI7DnHQN6kNqU9yP4BP7oFhZ8JZT4e9SKvq3CzNK2HR+iK+3FTCzv1G\nmGlJsUzM7M61xw/kmIE9OLJP19BclycAZGZmkpSUhNPpJCoqihUrVnDPPffwySefMHbsWF555RXA\nhLmVlJRwxx2h82DjAGEuVMURF8e2sm2kJ6YT44xp3fXblxv/3YwJJltgmC4cXFRey/Kte/l07R4+\nyy2k1uUlKS6K4welcvNJR3DCkJ5kpiR0GmE2xqJFi0i1slU2XM0tosLcQhVV06Lml+e3vtu7ezXM\nuRi69oHL3oCYhIDY6G9UldzdFcxft4clG4vZXFxFWY2ZnO+eEM1FE/oyfVQ6R2X2CLnHJo9/+zjr\n9/o32+zwHsO55+jWzSmE42puYS1UrDHqtvLNjEsb5/t1pZvh1fMhNgmueg+69AycjX7A5fGyIn8f\nizcUsWBdIVtKqhCBsf26cc6YPmSmJjJxQHeO7NOVqBATZyggIpxyyik4nU5+8YtfcNNNNx0Mc5s2\nbdrBMLf7778/2KY2SZgLVXHHOKl2V/s+41tVAnMuAvXCVe9Ct/6BtbGNbCutYsnGYr7aXMqXeSVU\n1LmJdgpHD+zB9ScO5LQRvemZ1NaUVcGhtS2fv/jyyy/JyMigqKiIU089leHDh0dsmJvfEJF+wCuY\n5GUKPK+qz7Qlr6+ql5oo8yzVJ6G664xrYPkuuPoDSA2dMCaAzcWVvJe9i49W72JzsXFqz+gWz/RR\n6UzNMrMiFjYlAAAPOklEQVSzXWLD+7c1GGRkmJXz0tLSOP/88/n222+ZPHkyEHlhbv7EDfxGVb+3\nshquFJEFmOyGrcvr6/VS5TDjsxaFqgof/cZkZ7j4Zeh3dHvvo93sKatlaV4xm4oq+WpzKTk7yxCB\nYwemcOWxAzh5WBoDOvkkUHupqqrC6/WSlJREVVUV8+fP5w9/+MPB4xEV5uZPrJSgu63XFSKSC2TQ\nlry+Xi/ljjqiHdGkJ6Y3X/E3z8KqV2Hyb+HI89tzC23G7fGSu7uCLzaaseYPBWUAxDgdZKUn8fsz\nszh7TJ9DnNpt2kdhYSHnn28+b7fbzeWXX87pp58OEDlhbgGvXCQTWIJJ7L1dVbtZ+wXYd2D7sGsO\n5vUdER8/4cGrJjPrNCcfXfBR0xWt+Y8JWcs6By6eDY6Om3CpdXlYtL6ID3N2s3h9EVX15ld7TL9u\nnDaiF9Oy0hjcs0vETgJF4mpujRHsMLeAISJdgP8Ad6hqecPunaqqiDSZ1xd4HmBkQoKWSS29Ewc1\nXdG2r+G/N0P/4+CCFzpEpPur6/l8fRFfbS5lwbpCympcpHaJ4ZyxGUwalMIxA3uQZreaNq0gKEIV\nkWiMSOeo6jvW7jbk9VX2Sw29E3s3fnzfNnjzCkjuBzNeh+jAiaOsxsXiDUV8nLObReuLqfd46ZYQ\nzUlDe3LhhL6cMDg15JdNsAldgjHrK5h8S7mq+pcGh96nlXl9AfZRRa+ERla/cNXCG1eA1w2XvwkJ\n/l00qt7tZeW2fSzbVMKXm0pYXbAfrxpXvZ8fN4Bzx/ZhZJ9k27ndxi8Eo0U9Hvg5kCMi2da++zAC\nnSsi1wPbgEt8Kaw2ShtvUT/7IxTmwOVz/fYYxuNVlm0q4c0VOw6ON50OYUzfZG6dOoSThvZkXL9u\ntjht/E4wZn2/pOlUo63L6wvURfNToW74BJb/E465GYb+rPVGNqCi1sVXm0v5YqMJDSuprKN7QjTn\njstgytCeHDcohaS40HQ7s4kcwv7peV3UYUItyoX/3Ai9R8MpD7apzFqXh8Ubinn/h518lltEvdtL\nQoyTKcN6cuaoPpwyIo3YqNAO1bWJLMJeqPUNW9SaffD6pcbB/rI3WjV5VOvysDC3iI/X7GbR+iKq\n6z2kJMZw+dH9OX1kb8b3727HbIYp1113HR9++CFpaWmsWbMGgL1793LppZeSn59PZmYmc+fOpXv3\n7ixbtoxf/vKXxMTE8O9//5shQ4awf/9+LrnkEubNm4ejAx/tNSTsv3kSG0dStLXM6se/hfKdZoY3\nOcOn69ftKufe/6xmwv8u4JbXv2f5llLOHZvBa9cfw/L7pvHAOUdy7BEptkjDmGuuueYn4WuPPfYY\n06ZNIy8vj2nTpvHYY2a5jaeeeoqPP/6Yp59+mn/+858APPzww9x3331BEymEeYvqiRISuqUaF7u1\n70LOWzDlPujb/DPm7aXVvP/DTj5cvZv1eyqIi3Zwzpg+nDcug2MGptiPUQLEnj/9ibpc/4a5xWYN\np/d99zV7zuTJk8nPzz9k33vvvcfixYsBuPrqq5kyZQqPP/440dHRVFdXHwx/27x5Mzt27GDKlCl+\ntbu1hLVQC9PjcA/PhOq9xo83fWyTibLr3B7mry3kla/z+S7f+PpPGNCdP5w1ggvH9yU5wZ4Q6kwU\nFhYeXI6xd+/eFBYWAjBz5kyuuuoq4uPjefXVV7nrrrt4+OGHg2kqEOZCdXldZnw6/36o3W9iS52H\nCi57x35mLdvKwtwiKuvc9O+RwN2nD+PcsRlkdIsPkuWdk5ZavmAhIgcDH8aOHcs333wDwJIlS0hP\nT0dVufTSS4mOjuapp56iV69GntsHmLAWqsfroXd9LWS/Bif8GnqPtPYrC3ML+deyrXyzZS9d46I4\na3Q6Z4xK58TBqfZzTht69erF7t27SU9PZ/fu3aSlpR1yXFV5+OGHeeONN7jtttt44oknyM/P5//+\n7/945JFHOtzesBaqovTe+o1xETzpHqrq3Pz72+3M/jqfHXtr6JMcx33Th3P5MQPsOE6bQzjnnHOY\nPXs29957L7Nnz+bcc8895Pgrr7zC9OnT6dGjB9XV1TgcjuCGv6lq2P7FZcbpskfT1LXkr/qPxZt0\n7IOf6oB7PtSL//GVfrR6l7rcHrUJLuvWrQu2CTpjxgzt3bu3RkVFaUZGhr744otaUlKiU6dO1cGD\nB+u0adO0tLT04PlVVVU6ZcoUra+vV1XVJUuW6MiRI3X8+PG6fv36Ruto7D6BFeqn73pQw9zaS+ag\nbrr8ujhuin2BH0ocnDysJ7dPG8K4/t2DbZqNhR3mFuZhbv4gtbaKL10nU+JM5OVrRzJlWFrLF9nY\nhCFhLVRVeD96OnNvPs6ewbWJaMJaqFucmbx542W2SEMcVY3ovE8dMXwMa7+4Qb2SGZzWJdhm2DRD\nXFwcpaWlHfJlDgaqSmlpKXFxgc3YEdYtqk3o07dvXwoKCiguLg62KQEjLi6Ovn3bsO5RK7CFahNQ\noqOjGThwYLDNCHtCrusrIqeLyAYR2WTl97Wx6fSElFBFxAn8P+AMYARwmYiMCK5VNjbBJ6SEChwN\nbFLVLapaD7yBScxtY9OpCbUxagawo8F2AXBMwxMaJuAG6kRkTQfZ1hypQIltAxAadoSCDQDD/FVQ\nqAm1RbRBAm4RWeEvF632EAp2hIINoWJHKNhwwA5/lRVqXd+dQL8G232tfTY2nZpQE+p3wBARGSgi\nMcAMTGJuG5tOTUh1fVXVLSK3Ap8CTuBfqrq2mUue7xjLWiQU7AgFGyA07AgFG8CPdoR1mJuNTWch\n1Lq+NjY2jWAL1cYmDAhboXaUq6GI9BORRSKyTkTWisivrP09RGSBiORZ/7tb+68Rkb8HyBaniKwS\nkQ+DaEM3EXlbRNaLSK6IHBckO35tfR5rROTfIhIXaDtE5F8iUtTw2X1TdVrHZlrfzw0i8rMG+ytb\nW3dYCrWDXQ3dwG9UdQRwLHCLVde9wEJVHQIstLYDza+A3AbbwbDhGWCeqg4Hxlj2dKgdIpIB3A5M\nVNWRmInHGR1gx8vA6Yfta7RO6zsyAzjSuuZZ63vbJsJSqHSgq6Gq7lbV763XFZgvZoZV32zrtNnA\neYdfKyJnisjXIpLaXjtEpC9wJvBig90dbUMyMBmzvi2qWq+q+zvaDosoIF5EooAEYFeg7VDVJcDe\nw3Y3Vee5wBuqWqeqW4FNmO9tQ1tSLVvObKnucBVqY66Gvi020w5EJBMYBywHeqnqbuvQHqDXYeee\nj/l1na6q/nBnexq4G/A22NfRNgwEioFZVhf8RRFJ7Gg7VHUn8GdgO7AbKFPV+R1th0VTdTb7HRWR\nXsBHwB9U9aOWKgmp56ihjIh0Af4D3KGq5Q1Ti6iqikjD51xTgYnAaapa7oe6zwKKVHWliExp7JxA\n22ARBYwHblPV5SLyDId1LzvCDmsceC7mh2M/8JaIXNnRdhxOI3U2RTSmm3yLqn7hS9nh2qJ2qKuh\niERjRDpHVd+xdheKSLp1PB0oanDJZiAJGOonE44HzhGRfEw3f6qIvNbBNoBpFQpUdbm1/TZGuB1t\nxynAVlUtVlUX8A4wKQh20EydzX1H3cBKwOdVtsNVqB3maiim6XwJyFXVvzQ49D5wtfX6auC9Bse2\nARcCr4jIke21QVVnqmpfVc3E3OvnqnplR9pg2bEH2CEiB6JCpgHrOtoOTJf3WBFJsD6faZi5g462\ng2bqfB+YISKxIjIQGAJ8ax1T4DpguIjc41Mt/srk3dF/wHRgI+aX8ncBrOcE641dDWRbf9OBFEz3\nJQ/4DOhhnX8N8Hfr9TjMF3mQH+2ZAnxove5wG4CxwArr/XgX6B4kOx4E1gNrgFeB2EDbAfwbMyZ2\nYXoX1zdVp3X+76zv5wbgjAb7K63/sRh32f9pqW7bhdDGJgwI166vjU2nwhaqjU0YYAvVxiYMsIVq\nYxMG2EK1sQkDbKHa2IQBtlBtbMKATilUEektIm+IyGYRWSkiH4uIv13L2oyIfBXm5Vc29rqJc/8m\nItsCbE83EfmfNl7b6tjRQNDphGq5nP0XWKyqg1R1AjCTwyItOsIOEWn0/VfVSYGsO9Dl+4oVjXQy\nECMiSQGsqhvQqFCb+xxCikC53oXqHyaKYkkTx+7EuKStwUTJHNififElfQFYC8wH4oHHMBEQB857\nALjLen0lxrczG3gOE9yciXEne8UqZwAm1OkHq85LG7qYNWVTU/ZYxxIbK/Ow+6xsqZzDzr8K4zL4\nA/Bqg/0/ucdG7K9s7L22js3G+C4vBI5t4pyf1N2G9+QNoMay88kmPoemPvsm7e/Q722wDejwGzaZ\nAf7ayP4JQI71Re9ifYDjGnwJ3MBYa3uu9SUdB3zRoIx1mIiJLOADINra/6z1hcvExJMea+2/EHih\nwfXJDb8cTdnUlD3NlXnYvVY2d1+HnXskxqc61do+4D/b6D0e/uVu6otulbsSEODvwA1NnHNI3W18\nTzKBNQ3KPfxzaO6zDwmhhn6T33GcAPxXVatUtRITOnVig+NbVTXber0SyFTVVUCaiPQRkTHAPlXd\ngYnmmAB8JyLZ1vYR1rXbVPUb63UOcKqIPC4iJ6pqWSts+ok9PpZ5OE2Vc4CpwFtqBVqr6oEMB83d\noy88jAmaVkxL2FhES2N1t+U9aYyGn0NLn33Q6YyB42uBi9pwXV2D1x5M1xfgLau83sCb1j4BZqvq\nzIYFWGOyqgPbqrpRRMZjonEeFpGFqvpQe+xpQ5lN3VdLNHqPPl0ocgwmj9A4Efl/QBzmB6a9tOZe\nqpo5FnJ0xhb1cyBWzKpwAIjIaMz45TwrxjEROB9Y6kN5b2LGWRdhRAtmzHWRiKRZ5fcQkQGHXygi\nfYBqVX0NM3Yaf9gpS1trkw9ltpbPgYtFJOXAvVj7fbrHJvgTcLaqZqqJsR1D4y1qY3W3+j0BKjBB\n403RljI7lE7XoqqqWrlznraCdmuBfOAOTJa5A8G9L1pd25bKW2vNWO5UK3eOqq4Tkd8D860ZRRdw\nCyanTkNGAU+KiNc655eHlf29iPzEJqtlbopmy2wt1v09AnwhIh5gFXBNM/fY7KMWETkFiFHVzxrU\nUSgiXUSkR4OudaN1q+o1rX1PVLVURJaJSfP5CSaDZcPjjb7PLb03HYkdj2pjEwZ0xq6vjU3YYQvV\nxiYMsIVqYxMG2EK1sQkDbKHa2IQBtlBtbMIAW6g2NmHA/weCu79B1FcuFAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1216376d8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "reach = 1\n",
    "lift_a = 0.05\n",
    "min_lift_uplifts = [0.01, 0.02, 0.05, 0.1]\n",
    "convs_ctrl_a = np.linspace(100, 100000, 100)\n",
    "ctrl_scale_factor = 1\n",
    "\n",
    "fig = plt.figure(figsize=figsize)\n",
    "ax = fig.gca()\n",
    "\n",
    "plot_power_by_sample_size_mc(ax, reach, lift_a, min_lift_uplifts, convs_ctrl_a, ctrl_scale_factor, alpha, n_iter)\n",
    "\n",
    "plt.savefig('img/multi_cell_power_vs_conv.pdf', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Minimum sample size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10%: 2745, 219596\n",
      "5%: 10754, 860346\n",
      "2%: 67453, 5396260\n",
      "1%: 264745, 21179569\n"
     ]
    }
   ],
   "source": [
    "reach = 1\n",
    "lift_a = 0.05\n",
    "ctrl_scale_factor = 1\n",
    "\n",
    "for min_lift_uplift in [0.1, 0.05, 0.02, 0.01]:\n",
    "    min_conv_ctrl_a = get_min_conv_ctrl_a(reach, lift_a, min_lift_uplift, ctrl_scale_factor)\n",
    "    N = 4 * min_conv_ctrl_a / 0.05\n",
    "    print('{:.0f}%: {:.0f}, {:.0f}'.format(100 * min_lift_uplift, min_conv_ctrl_a, N))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:expt_dashboard]",
   "language": "python",
   "name": "conda-env-expt_dashboard-py"
  },
  "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.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}