docs/Notebooks/Example1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Example 1: Finding the habitable zone\n",
    "\n",
    "In Section 6 of [Bixel & Apai (2021)](https://ui.adsabs.harvard.edu/abs/2021AJ....161..228B/abstract), we propose that the concept of a \"habitable zone\" could be validated by searching for a region of space where planets with atmospheric water vapor are statistically more frequent.\n",
    "\n",
    "In this example, you will use Bioverse to determine whether a [LUVOIR](https://ui.adsabs.harvard.edu/abs/2019arXiv191206219T/abstract)-like imaging survey could test this hypothesis.\n",
    "\n",
    "## Setup\n",
    "First, we'll import the Bioverse code:"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "source": [
    "# Import numpy\n",
    "import numpy as np\n",
    "\n",
    "# Import the relevant modules\n",
    "from bioverse.survey import ImagingSurvey\n",
    "from bioverse.generator import Generator\n",
    "from bioverse.hypothesis import Hypothesis\n",
    "from bioverse import analysis\n",
    "\n",
    "# Import pyplot (for making plots later) and adjust some of its settings\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.rcParams['font.size'] = 20."
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "For this example, we will use the LUVOIR-like imaging survey and host star catalog."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "source": [
    "generator = Generator('imaging')\n",
    "survey = ImagingSurvey('default')"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Injecting the statistical effect\n",
    "\n",
    "The first step is to inject the statistical effect we are searching for into the simulated planet population. Specifically, we will simulate the likelihood that a planet has atmospheric water vapor as follows:"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "source": [
    "def habitable_zone_water(d, f_water_habitable=0.75, f_water_nonhabitable=0.01):\n",
    "    d['has_H2O'] = np.zeros(len(d),dtype=bool)\n",
    "\n",
    "    # Non-habitable planets with atmospheres\n",
    "    m1 = d['R'] > 0.8*d['S']**0.25\n",
    "    d['has_H2O'][m1] = np.random.uniform(0,1,size=m1.sum()) < f_water_nonhabitable\n",
    "\n",
    "    # exo-Earth candidates\n",
    "    m2 = d['EEC']\n",
    "    d['has_H2O'][m2] = np.random.uniform(0,1,size=m2.sum()) < f_water_habitable\n",
    "\n",
    "    return d\n",
    "\n",
    "generator.insert_step(habitable_zone_water)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Next, let's make a simulated dataset using this modified Generator object and the imaging Survey. Let's start with a relatively optimistic assumption that 75% of EECs are habitable and only 1% of non-EECs have \"false positive\" water vapor."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "source": [
    "sample, detected, data = survey.quickrun(generator, f_water_habitable=0.75, f_water_nonhabitable=0.01)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Let's take a look at the simulated data set by plotting which planets have H2O versus their insolation. You might notice that planets within the habitable zone (approx 0.3 < S < 1.1) are more likely to have water-rich atmospheres."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "source": [
    "x, y = data['a_eff'], data['has_H2O']\n",
    "\n",
    "# Now plot the water-rich/water-poor planets versus a_eff (in log-space)\n",
    "fig, ax = plt.subplots(figsize=(16,2))\n",
    "ax.scatter(x,y)\n",
    "ax.set_xscale('log')\n",
    "ax.set_yticks([0,1])\n",
    "ax.set_yticklabels(['no H$_2$O','has H$_2$O'],fontsize=24)\n",
    "ax.set_xlabel('Semi-major axis (scaled) (AU)',fontsize=24)\n",
    "\n",
    "# Highlight the boundaries of the habitable zone (water-rich planets should be more common here)\n",
    "ax.axvline(0.95,linestyle='dotted',c='black',lw=3)\n",
    "ax.axvline(1.7,linestyle='dotted',c='black',lw=3)\n",
    "\n",
    "plt.show()"
   ],
   "outputs": [
    {
     "output_type": "display_data",
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA/UAAAC5CAYAAABpyR9zAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA0PklEQVR4nO3debwcVZ338e83NxcS1hA2SViCjmyOGwZxZERQEB5cwAUVRGVmBB8VHVzQAVFhQBkHxwV3HDUioA4ucRk1oI9xl5GIgoqBUYMQUFmyQiAh9/f8cU7ddDrV3dV9e7nL5/161atvV51z6neqqrvu6Tp1yhEhAAAAAAAw8UwbdAAAAAAAAKAzNOoBAAAAAJigaNQDAAAAADBB0agHAAAAAGCColEPAAAAAMAERaMeAAAAAIAJikY9AAAAAAATVOVGve1ltsP2Eb0LBwAAAAAAVDXlrtTbXpB/nFjczbRtrH+27bNt/9D2n22vt/0X2z+yfY7tnbu1LgAAAADA5DZ90AFMJbZPlvRhSbPyrBFJqyTtImk3SYdJOsv2ayLiyoEECQAAAACYMKbclfpBsf1KSZcrNeiXSDpO0syImC1phqRjJf08L788pwcAAAAAoCGu1PeB7cdLukSSJX1V0okRsaFYnv9eZPu7kq6SdIKkS2xfGxG/rLKOXXbZJebNm9flyAEAAAAAg7ZkyZK7I2LXsmUdNeptz5Z0rqTnSpoj6W5J35T09oi4s0H6Z0t6jqRHSZqr1EvgVknflvSeiLijJN9Wkl4l6UWSDpK0raR7Jf1Z0o8kXR4RP+2kDn12oaStJN0h6WW1DfpaEfGQ7VMl3SRpD0kXKG23lubNm6frrruuO9ECAAAAAMYN27c2WtZJo35PSQsk7SPpfkmh1LB/haSjbB8cESvq8pwj6Y0171dLminpwDydYvuoiLihJujpkq6W9NQ8K5TuP99Z6f7zx+S/x3Wj3vaekv5PfvuhiFjdLH1ErLL9IUnvlPRM23tGxO29jhMAMDgLr1+uixct1R0r12nOrJk665j9dcLj5zacP5YyG6VbvnKdhmxtjNDcHq2rk3xlaSRtEfOsmcOypZX3b2gZQ6vt3e626LQeVbdvNwx6/Rg79iGARjpp1H9Q0jJJJ0XET3Pj+zhJn5E0T9LZkt5cl2e5pH+T9F+SbomItbaHJD1OqfF6jKQrbT86IiLnOVmpQX+/pFdK+mJEPJDzzVW6gr1DB/H321OVut1L0sKKeRYqbRdLOlwSg+YBwDi33377jf598803V8638PrlOvvLN2rdho2SpOUr1+nsL9+o6269V19asnyL+ZJa/iPfqMz6vPXpNuZTcC/W1Um+sjRnXfUrydKGjbFZzCvXbeoE1yyGqtu76rZYeP1yvfCoQ1X8+xKnX1qpHlW3bzcMev0YO/YhgGY6GSjvQUlHFd3eI+KhiPiaUhdzSXpBfYaIeF9EnB0R10fE2jxvY0QskXS8pN8qdcs/vCbbk/LrZRFxeUQ8UJPvTxHx4Yi4qIP4C0/Oj5RrOCl1+x+rg/Lrg5KWVszzO0nr898HdiEGAECP3XLLLaNTOy5etHT0H/XCug0b9blrbyudf/Gi1qeSRmXW5y1L16t1dZKvLM2GkRht0DfTKIZ2tnersoryNty7XA+tuEMPrbijcj2qbt9uGPT6MXbsQwDNdNKovzQi7imZvzC/7mt726qFRcSDkq7Jbw+rWVR0U9+j7QirGZa0e4tpRhfWMzu/roiIkSoZcrriFoaGz623fbrt62xfd9ddd40xTADAINyxcl3p/I1R3nBtlL5Kmvr5rcrq5ro6yVdl/e2uo93t3SrfWOox1vpVNej1Y+zYhwCa6aRR//MG85fX/D2rfqHtA2x/yPYNtlfbHrEdtkPSP+dkc2qyfCu/Hm/7a7afZ7thA7cD348IN5uUbinYgu1DbF9i+0bba23fYfsbtueXJe8wvpb5IuLSiJgfEfN33bV0IEQAQJ8sXbp0dGrHnFkzS+cPufw00Ch9lTT181uV1c11dZKvyvrbXUe727tVvjmzZmrOaR8fncrSd7qNumXQ68fYsQ8BNNNJo35N2cyie3w2XLvM9osl3SDpNZIerTSK/SpJf8nTfTnp6BX+iPi+pLdLekjp/vkvSbrb9k2232P7kR3E3i1vkfRCSYslvV7S+yTtL+la28fVpS16Nexku9L2zulm5bf3jjVYAEDv7bfffqNTO846Zn/NHB7abN7M4SGddOhepfOLgeI6KbM+b1m6Xq2rk3xlaYanWcNDrX8vbxRDO9u7VVlFeTvsvreGZ8/V8Oy5letRdft2w6DXj7FjHwJopufPqbe9q6RPKDX0vyDpYkk31D7WzfYFSo/I2+wsHREX2L5c6d72IyT9naQD8vTPtv8pIi7rdR1KvFfSyRFR3Pcu2/8p6TdKYwt8sybtTfl1a6WG/01q7QClR+BJabwBAMAkVQxyVTaq9fx9Znc02nWzMhul63T0+6rr6iRfozRlMVcd/b7K9m5nW4ylHv0a4GzQ68fYsQ8BNONocQ/ZaEJ7mdJj7I6MiMUN0hSF7RsRy/K8lyl1Y/+tpEeX3Vdu+1JJp0n6TESc2iSGIUlPkXS+0qB690l6eET8tVIlUhkLJL1cqfv9Ed1Km9NfIen5ETGjZt6ekv6k9IPFOVUG97N9jtLo9yFp7yqPtJs/f37wnHoAAAAAmHxsL4mIstu9O+p+36498+sNDRr0lvS0KgXlke8XS3qWpA1K3fVLKzYgcyTdXTsjN8iL8QHOsN30MXx5+Rn57Td5Rj0AAAAAoJF+NOpX5de/zQ34eqdJekT9TNtblaQtrJdUPNdj67GF1x22D1N6Jv0XSha/XelHiDmSLrM9XJJGtqcr9WrYI6d/e2+iBQB025w5c0YnTC3sewDAIPWjUf8dpW7kfyvpEtuzpHRF2vZZkj6sTYPJ1brM9qdtH2N7+2Km7XlKDd8ZktZJ+mFvw2/N9m6SrpR0q9KtAZuJiCVKA+pJ0vGSfmL72KJxb3u67WdI+rGkE3K6MyPiF72OHQDQHXfeeefohKmFfQ8AGKSeD5QXEUttv1+pUXuGUhf0FZJ2kDQkaZGk6yS9tS7rDKUB8k6VFLZXKQ0et01evlHSKyPibg1Q/sHhm5K2l3R4RKwuSxcRH7a9WtIHlW4Z+JakEdsrJe2otC0kabWkMyLis72OHQAAAAAwsfW8US9JEfEG2zdJepWkg/J6fynps5I+JOltJdn+RenK9dMkPVKpS/qQpN9L+oGk90fEDT0PvgnbMyV9XdKBko6OiF83Sx8Rn7X9TUmvlHScpP2UHl13r6RblBr6Hxv0DxUAgPYtX7580CFgQNj3AIBBqjz6PTaXu84vlHS0pOdExLcHGQ+j3wMAAADA5NRs9Pu+XKmfbGxPk3SFpGMlvWTQDXoAAAAAwNREo74z75F0oqRrJE23fUrtwoi4fCBRAQAAAACmFBr1nTk4vx6dp3o06gEAAAAAPUejvgMRccSgYwAAjC/bbz/69FWtWbNmgJGg39j3AIBBolEPAEAXrF27dtAhYEDY9wCAQZo26AAAAAAAAEBnuFIPAEAXrF69etAhYEDY9wCAQaJRDwBAF9TeV42phX0PABgkut8DAAAAADBB0agHAAAAAGCColEPAAAAAMAExT31AAB0ge3RvyNigJGg39j3AIBB4ko9AAAAAAATFI16AAAAAAAmKLrfAwDQBXS7nrrY9wCAQeJKPQAAAAAAExSNegAAAAAAJiga9QAAAAAATFDcUw8AQBesWbNm9O/tt99+gJGg39j3AIBBolEPAEAX7LDDDqN/M3Da1MK+BwAMEt3vAQAAAACYoLhSDwBAF2y33XaDDgEDwr4HAAwSjXoAALqg9r5qTC3sewDAINH9HgAAAACACYpGPQAAAAAAExSNegAAAAAAJijuqQcAoAvuuOOO0b/nzJkzwEjQb+x7AMAg0agHAKAL5s6dO/o3zyqfWtj3AIBBovs9AAAAAAATFFfqAQDogj322GPQIWBA2PcAgEGiUQ8AQBfU3leNqYV9DwAYJLrfAwAAAAAwQdGoBwAAAABggqJRDwAAAADABMU99QAAdMHNN988+vd+++03wEjQb+x7AMAg0agHAKAL9t9//9G/eVb51MK+BwAMEt3vAQAAAACYoLhSDwBAFzzykY8cdAgYEPY9AGCQaNQDANAFtfdVY2ph3wMABonu9wAAAAAATFA06gEAAAAAmKBo1AMAAAAAMEFxTz0AAF2wZMmS0b+f8IQnDDAS9Bv7HgAwSDTqAQDogvnz54/+zbPKpxb2PQBgkOh+DwAAAADABMWVegAAuuDggw8edAgYEPY9AGCQaNQDANAFtfdVY2ph3wMABonu9wAAAAAATFA06gEAAAAAmKBo1AMAAAAAMEFxTz0AAF2wePHi0b+POOKIgcWB/mPfAwAGiUY9AABdcOSRR47+zbPKpxb2PQBgkOh+DwAAAADABMWVegAAuuCpT33qoEPAgLDvAQCDRKMeAIAuqL2vGlML+x4AMEh0vwcAAAAAYIKiUQ8AAAAAwARFox4AAAAAgAmKe+oBAOiCr3/966N/P/vZzx5gJOg39j0AYJDM81Qnh/nz58d111036DAAYMqyPfo359aphX0PAOg120siYn7ZMrrfAwAAAAAwQdH9HgCALnjWs5416BAwIOx7AMAg0f1+kqD7PQAAAABMTnS/BwAAAABgEqJRDwAAAADABEWjHgAAAACACYqB8gAA6IIrr7xy9O+TTz55gJGg39j3AIBBYqC8SYKB8gBgsHhW+dTFvgcA9BoD5QEAAAAAMAnR/R4AgC446aSTBh0CBoR9DwAYJLrfTxJ0vwcAAACAyYnu9wAAAAAATEI06gEAAAAAmKBo1AMAAAAAMEExUB4AAF1w6aWXjv59+umnDzAS9Bv7HgAwSAyUN0kwUB4ADBbPKp+62PcAgF5joDwAAAAAACYhut8DANAFp5122qBDwICw7wEAg0T3+0mC7vcAAAAAMDnR/R4AAAAAgEloyjXqbS+wHbYXdzNtmzHMtn227R/a/rPt9bb/YvtHts+xvXM31wcAAAAAmJy4p77PbJ8s6cOSZuVZI5JWSdpF0m6SDpN0lu3XRMSVAwkSAAAAADAh0KjvI9uvlPRRSZa0RNLbJH03ItbbHpb0NEkXSDpE0uW2t4+Ijw8sYABAZf/xH/8x+vcb3/jGAUaCfmPfAwAGacoNlGd7gaSXS/p+RBzRrbQV1vt4ST+TtJWkr0o6MSI2lKSbLukqSSdIWi/p0Ij4ZavyGSgPAAaLZ5VPXex7AECvDWygPNvL8j3pR+T7yN9r+4+2H7S93PYnbO/Roozn2f627btyvtttX2H74F7G3gMXKjXo75D0srIGvSRFxEOSTpV0Z05/Qb8CBAAAAABMLP3qfr+npAWS9pF0v6SQNEfSKyQdZfvgiFhRm8H2NEmflvSyPGujpDWS5ko6WdKLbZ8RER/tSw3GwPaekv5PfvuhiFjdLH1ErLL9IUnvlPRM23tGxO29jrMXFl6/XBcvWqo7Vq7TnFkzddYx++uEx8+tvLxbZYxn3Yq/rBxJOu9rv9HKdek3pJ22GdY7nv2o0fLL8lx367363LW3aWOEhmyddOheuvCER4+pLs3q2M/9t/D65U23R5W6lKVZvnLd6DxL2marId2/fmNpvnMX3thw+7YTXyd1r7Kd69MdecCu+t7v7mqar75OT3r4Tlp2z7q219Xuvm/nuGpVj1Z1mLfzTP3kD/equAi7zfA0vet5j9msjDe84Q1jrmNt+hnD0/TgQyMaCZV+Ftv5btxx5rBsaeX9Gyrv127so16p/ewN2doYobkNvnN2nDmsDRtHdN/6jZKkmcPTNGN4aHRbdFKn2vJnbTOs2Yc+VxtG0sHxuPOv1nnPaf69Uht3Wfz166j9Tj//67/Rivu3vC4wa+bw6Hrr43tww0bdv2Fki3Rl9anyObel4WnW+o2b90qwpZnTp2ndhpHNYq763VP2/dcof9XP/6xthhUhrVq3+bFfuw9mNThGVty/oel+anYctqp37fYsNDoGlq9cp2mWRmo29zRLf/fw2Zt9T1X9XPdC1XNm7T6ute1WQxoemja6n5r9DyNV37ZV/4/Z7PjW5ufxIw/YVf99w52jn7uyz1CVbdPo+6oshmmWtp4+TQ9sGBmNJVR+LuimTr/ze3muqC175vA0rXtoRNHgvFi1rPp9McjPTrf0tPu97WVKDfmVkpZJenVE/DR3MT9O0meUBoy7OCLeXJf3XyRdpPQDwNslfSAi1tieK+l9kk5UGmTuyIj4QRsxLVCfu9/bfomky/PbgyLipgp5DpL0m/z2Ja0GzRuP3e8XXr9cZ3/5Rq3bsHF03szhIV30vEePngybLe9WGeNZt+IvK2d4yNq4MTRSl3Z4yLr4BY+VpC3yTJO2SC9Jpzxp75Zfmo3q8vwnzNWXliwvrWNZDL3afwuvX66zrvrV6D/ehWJ71P+D28mxWaY237kLb9TlP/vTFmlOedLemr/P7MrxtavqcValTvX5GtWpWZ52YuqkTtKWx1WzmKrUocw0S+994eMaxttuHats/+Kz2Ml3YzO92Ee90qxujb5zmmm3TlW27fA06+ITW3+vNIpF2vIYHp5mjUjaONL4/7bhadaLnrhXy/rXxldlP3f6GRkeshTa7Hut0bFW9v03zdLQNGvDxs3zt3te6bZWx1mrejfbnt2sR78+r1XPmWX7uJFG/8MMTbOmqf1t2+z/mE6O77LPeJlW31ednoeq/F/Wrk6/83t5rmjnvNiNsgrj4VxXZjw8p/5BSUdFxE+l1MU8Ir6m1CVdkl5Qm9j2tpLOzm/fHREXRsSanHe5pJMk/Ugp/gvVmSfnx8k1nCS9qMOy6x2UXx+UtLRint8p3VMvSQd2KY6+unjR0i0+OOs2bNTFi5ZWWt6tMsazbsVfVs6GkpNhMf/iRUtL85Sll6TPXXtbRzGs27BRn7v2toZ17Of+u3jR0tJ/JortUZ+2k2OzTG2+Rtvxc9fe1lZ87aq6navUqT5flWOj6rra2ffN8rdbjyp1KDMSahpvu3WsEncRayffjc30Yh/1SrO6NfrOaabdOlXZthtGqn2vNIql9Dt9JJo26Is0VepfG1+V/dzpZ2TDxtjie63RsVb2/TcS2qxBX+Rv97zSba2Os1b1brY9u1mPfn1eq54zqzbopcb/w2wc6WzbNtvmnRzfZZ/xMq2+rzo9D3X6mWym0+/8Xp4r2jkvdqOswng417WrX93vL42Ie0rmL5T0Hkn72t42Iu7L858haQelRu2/12eKiI22L5C0SNJTbD8sIv7cZkzDknZvM0+nZufXFRHRqN20mYgYsb1CKcbS59bbPl3S6ZK09957dyPOrrqjpkty2fxWy7tVxnjWrfh7nX5jhR49jcpslLdZDL3Yf+2sbyzHZrN8jbbFxoiebo+qx1nV9dSmq3JstLOuTmLoJH9t2qp1aHd97cZYJfYi1k6/G5vp9j7qlVbr72R/dnLctJturMd2Ve1+Jqvs57F8RlqV3SyGRjo5r3TbWI6zVnm7WY9+bJNunzPHEkOz82wjnR7fVerUKk2n56FufyZrY6k6f6z5xhJTrU7/D+l2+kHr15X6nzeYv7zm71k1fxeD4P2q/l77Gj+Q9FBd+nZ8PyLcbFK6PWALtg+xfYntG22vtX2H7W/YLu0OoXSbbSea5ouISyNifkTM33XXXTtcRe/MmTWz6fxWy7tVxnjWrfg7Sd9OniG3PoQbldcob7MYerH/mpVZv2wsx2azfI22xZDdVnztqrqdq66nNl2VY6OddXUSQ/38dsuoWod24mi2bCx1L2Lt9LuxmW7vo15ptf5O9mcnn+d207WTbyzbuN3PZJX9PJbPSKuym8XQSCfnlW4by3HWKm8369GP7dHtc+ZYYmh2nm2k0+O7Sp1apen0PNTtz2RtLFXnjzXfWGKq1en/Id1OP2j9atSvKZsZEQ/UvB2u+btoodY2+svyFlf/+92ifYukF0paLOn1Svf47y/pWtvHlaQv4twpDwDYUk43K7+9dyzBDspZx+yvmcNDm82bOTw0OvhJq+XdKmM861b8ZeUMD7n0Az48ZJ11zP6leRodnCcduldHMcwcHtJJh+7VsI793H9nHbO/hqdt+cVfbI/6tJ0cm2Vq8zXajicduldb8bWr6nauUqf6fFWOjarramffN8vfbj2q1KHMNGuzeM8777zRqVWMZarEXcTayXdjM73YR73SrG6NvnOaabdOZetf+aMrRicp3W87ls9X6Xf6NGuo5DuiPk2V+tfGV2U/d/oZGR7yFt9rjY61su+/ac73p9flb/e80m2tjrNW9W62PbtZj359XqueM8v2cSON/ocZmtbZtm22zTs5vss+42VafV91eh7q9DPZTKff+b08V7RzXuxGWYXxcK5rV7+633dq60EH0MB7JZ0cEcU977L9n0oD210o6Zt16YuB8bZWavy3HChP0gFKj7STpN+OKdoBKQaXaDQaZqvl3SpjPOtW/I3KkZqPfl+Wp9PR75vVZf4+s5vWsR/7ryizyujy7R6bVUe/L7Zjs+3bi9Hvqx5nZelajQhbVqcqo9+P9divkr9qParUocro9+eff/7o3+edd17bdaxP32z0+3a/GzsZ/X68fr/Wf/bKRpOu/c7p9uj39dtl1jbDuvXHnxtdPu/oU0tHxm4Ud7PRsDsd/b62/q1Gv6+yn+s/I90e/b7Z93Oj/FXPK70c/b6Iod3R7+u3Z6HRMTDeR79v55zZ69Hvq5xn621xfKt7o99X+b4qi2EQo993+p3fy3NFfdljGf2+2b5g9PtWhW8a/f7IiFjcIE0RwL4RsSzPu0DSuZJ+GxGPapBvhlIPgOmSnhkR9Q3pRjEtUI9Gv7d9haTnR8SMuvl7SvqT0v/750TERRXKOkfpkXYhae9Wj7Qbj6PfA8BU4pougL08t2L8Yd8DAHqt2ej34/VK/S/y6yNtz80j3tc7XJvi/0XJ8kGYI+nu+pkRcbvtbyk9xu8M2x9u9qx62ztIOiO//WaVZ9QvWbLkbtu3dhDzjpJWdZAP1U3lbTwZ6j6e6zBeYhtEHP1aZ0frcbV7/HZRyTkDE1vFfd9r4+W7YTKbytt4MtR9PNdhvMTGub0zvTy379NowXht1F8tabXSCPhnSTqzdqHtIUlvy29/2MHI911n+zBJT1W6v77M2yUdrdTwv8z2iRGxRf8j29OVBujbQ9KGnK+liOhoXAHbl0bE6Z3kRTVTeRtPhrqP5zqMl9gGEUe/1tnL9di+rtEv7sBYjJfvhslsKm/jyVD38VyH8RIb5/aOyx7Iub1fA+W1JT/a7l357etsv9X2dpJke66kz0n6e6XHap87mCg3sb2bpCsl3Srp/LI0EbFEaVA9STpe0k9sH2t7OJcx3fYzJP1Y0gk53ZkR0eteCF/vcfmY2tt4MtR9PNdhvMQ2iDj6tc7xso2BdnDc9t5U3saToe7juQ7jJTbO7RPIuLynPs8fkvQpSS/LszYqXb2fpXRv+oik10bER9qMaYG6eE+97e0lfU/SwyUdHhG/blHmSyV9UKnbh5TqsTK/L4ZkXC3pjIj4bKv6AAAmJq7UAwAwuXClvk5EbIyIl0t6gVJ3/JWStpN0p9KV+ie226DvNtszlX7pOVDSs1o16CUpN9QfIemtSlfl75G0vdJj636idFvBI2jQA8Ckd+mgAwAAAF01kHN7T6/UT2a52/xCpfvknxMR3x5sRAAAAACAqWa8DpQ3rtmeJukKScdKegkNegAAAADAINCo78x7JJ0o6RpJ022fUrswIi4fSFQAAAAAgCmF7vcdsL1Y6fF1pSJiXDykFgAw8dl+tdLjXfeQ9BulJ6P8cLBRAQCATtg+XNKbJD1B6XHn/xARC8ZS5rgdKG88i4gjIsKNpkHHBwCYHGy/SNIHlB7z+nilAVW/ZXvvgQYGAAA6tZ2kX0v6Z0nrulEgV+oBABinbF8r6YaIOK1m3i2SvhgRZw8uMgAAMFa21yo9ynzBWMrhSj0AAB2y/QLbH7T9Q9urbYftpuOq2N7T9qds32H7QdvLbL/f9k516bZS6pp3dV0RV0t6cndrAgAApN6e23uFgfIAAOjcuZIeK2mtpNslHdAsse1HKHWh303SVyX9TtITlbrgHWv7sIi4JyffRdKQpL/UFfMXSUd1qwIAAGAzvTy39wRX6gEA6NzrJe0naQdJr6qQ/iNKJ/3XRcQJEfEvEfE0Se+TtL+kd5bkqb9PziXzAABAd/Tj3N5VNOoBAOhQRHwvIm6JCgPU2H64pGdIWibpw3WL3yHpPkkvtb1tnne3pI2SHlaXdjdtefUeAAB0QY/P7T1Box4AgP54Wn69OiJGahdExBpJP5a0jaQn5XnrJS2RdHRdOUcrdfMDAACD1da5vVdo1AMA0B/759ebGyy/Jb/uVzPvvZJOtf0K2wfa/oDSM20/1qMYAQBAdW2f221vZ/txth+n1B7fO7/v+HG1DJQHAEB/7JhfVzVYXsyfVcyIiC/Y3llp0J49lJ5re1xE3NqrIAEAQGVtn9slzZf0vZr35+fpM5JO7SQIGvUAAIwPzq+b3cMXER9RGoQHAABMLFuc2yNicc38rqD7PQAA/VH8Wr9jg+U71KUDAADj27g4t9OoBwCgP5bm1/0aLH9kfm10Xx4AABhfxsW5nUY9AAD9Udw/9wzbm51/bW8v6TBJ6yT9rN+BAQCAjoyLczuNegAA+iAifi/paknzJL2mbvH5kraVdFlE3Nfn0AAAQAfGy7ndEdE6FQAA2ILtEySdkN8+TNIxkv4g6Yd53t0R8aaa9I9Qesb8bpK+KukmSYdKOlKpa96TI+KefsQOAAC2NBHP7TTqAQDokO3zJL2jSZJbI2JeXZ69JP2rpGMl7SzpTkkLJZ0fEff2JFAAAFDJRDy306gHAAAAAGCC4p56AAAAAAAmKBr1AAAAAABMUDTqAQAAAACYoGjUAwAAAAAwQdGoBwAAAABggqJRDwAAAADABEWjHgAAAACACYpGPQBgUrB9qu2wvXjQsXSL7Xm5TjHoWMYD20fk7bGsz+t9Z17vP/Rzvd1ge3GO/dQ+r3dBXu95dfOfnOdf1aX1HGD7Idvf70Z53WZ7O9v32L7N9oxBxwNgcqJRDwBTlO3puSH8bdt32l5ve4Xtm2z/t+232D5k0HECg2T7YZLOlLRM0mcHGswkEBE/kbRY0gtsz+9Cke+SNCTpglYJbZ9Z/Ehm+9IK6YsfJhZ3mjYi1kq6RNKeks5oVQ4AdIJGPQBMQbZ3lfRTSZ+WdIykh0naKMmS9pd0nKR/k3TNoGLswCpJSyX9adCBdNEGpTotHXQg48T9Stvi931c5zskbSPp3RHxUB/XO5ldmF8vGkshtg+V9FxJ10bEdypkeXnN3y/s45XzD0haI+ls2zv0aZ0AphAa9QAwNV0uab7SP5pvlrRHRMyMiFmSdpR0tKSPSFo5qADbFRFfiYgDIuJlg46lWyJiea7TAYOOZTyIiP/J2+Pp/Vif7Z0knSppndJnBl0QEd+V9EdJR9l+9BiKekN+rXLV/TGSHqfU4+Iape+5E8aw7soiYqWkL0marc1/WACArqBRDwBTjO0DJD0jv/3HiLg4Iv5cLI+INRHxnYh4jaQDBxIkMD6cImmGpG/kbtTons/n13/qJLPtnSUdL2m9pC9XyFI0pj+nTT/Q9LOB/bn82lF9AaAZGvUAMPXUXhn7RrOEEbGu0bI8ANQ5tn9ue5XtB2zfYvsS23s1yDM6aJftHWz/u+3f215n+w+2/7W2S6ztp9teZPtu2/fZ/oHtpzQou+OB8uoHYLN9jO3v2L7X9krb19j+u5r0O+bB027Osd9m+922ZzYofw/br8pjFdxi+37bq21fb/t827Ma5Gs5UJ7tx9u+PMfwYN5Wi2w/v0meZbncI2zPtf2RvP0ftP3LiptNtg+1fZHtn9lensdl+Gsep+EFDfIca3skT89okOacHN8q2/Nq5jccKM/2tHwMfM9pYLINtu+y/Rvbn7J9bNV61SgGxvuvRglsP9b2ZXmbPmh7Td6W3873cG/TIN+Btj+Wj6H78nF2Y/78PKEu7Va2n2n7E7Z/lffxA7ZvtX1Fffp25O320nyM35X34R22v+DUvb1Z3kNtfz1/Ttba/qXtf7Zd5f/LYpueYnu4g9BfImlrSdfkK+HN4pye00vSFZK+otT74mjbe3Sw7k58V9I9kh5r+/F9WieAqSIimJiYmJim0CTpREmRp0d0WMaBSt1Yi3I2SFpb8/5eSYeV5Fucl79e0k3577VKV9uKvF/LaV8taUTpXv9VNcsfbFD2qXn54g7qc0TOu6zJetdJ+ntJu0q6sSb2B2vSfKNB+V+sSROSVuTyi/f/K2nPknzzijQNyj29rpwVkh6qef9ZSUMl+Yp9d7qku/Lf9+X6/LLiNtuurk7rJa2um/fxBnk/mJcvlzS7btnja46HUxvtp5Iyr6hb98q6ffOzNo+JXWry7tEgzXF1x+4DdcdMSDqgJN9r6/bTWqXxAor3i+vSP6uuzPvy8Vj7+XtpgxgXl23LvGx7pa7oRTkjdfFvlHRGg3JfXFeHFTmOUDreP5P/Pq9B/mk1x8uhHXxmv5Hznl0h7TNz2l/VzPt8nndWk3wLyvZHp2kl/XdO85Z268vExMTUbOJKPQBMPUtq/v6w06B5ldneUdI3Je0jaaGkgyXNjIjtJO2r1JDcSdKXGl2BVhp8zJKekvNtJ+k0pUbCs22/TdL7lQbr2zkidlRq4P5U0laS3tdOzG3YNZd9Uc16983rnSHpvZI+JmlY0lOUGkXbS3pFjv2Zto8rKfcWSedKepTSttopl3eEpJ9LeoSkj7cTqO0nS/qoUuPoi5L2yuXOkvRWpcbDKZLOblLMf0i6U+lHkm3zvii9wl5iROk4OEnSXEkzImIHpX3/WqWG6um2TyzJ+2ZJv5M0R2l7FnWaodQ1eljSlyNiQZVAbB8u6eQc0+sl7RBpfIgZeR2nSvpRxXoVDsuvd0bEnQ3SfDDH+g1J+0fEjHzM7CjpcEmfUGro18Z6otJo6ENK++2gvN23zbGeos0/o1Lalp+W9HRJu+R9NVPpM/h+SdMlXWp77zbreJmkoyTdoNTw3TbHv5Okc5SO6Q/YPqw2k+1H5HiGJF2t9OPgTrneb1S6V/34ZiuOiBFJ1+e3pb1vGrFtSUXPmfptVebl+fWKmnlX1C3rh+vya1v1BYCWBv2rAhMTExNT/ydtuopWXPn+jtKI1MdL2rVF3gtzvoWS3CBNcUXqTXXzF2vTlcW/Kcn3yZq4PlWyfB+lhltI2rtu2aka+5X6kPTpkuV716x3fYvYt4i7xbpnS/przrtv3bJ5RVwl+b6bl/1I5Vfj35WXr1Fq5NYuW6ZNV1d379Ex9tK8ju81WH6wNl3lfmme9/78/k6lxmuj/bSsbv6b8/xvdTH+d+Yyv91g+W41x0ylbaj0A8BtOc+VXYy1OPbeUbKs+MydWjf/qDz/j6rrLVGyXb/RYH2/U/oxpz7fuTXb5rwmcV+S01zVZn0fWVP+w1qk3Unph5URpR++avfF3bmMJzTIu0DdvVL/vJzmrm7teyYmJqYIrtQDwFR1mtJV5/VKV76frnR1d6Gkv9r+H9svyVfE6hVXtt4XEdGg/GJQqKMbLL8qIv63ZH7tY6m2eNxVRNyq1FVdkv62QdljVbbePyldbZcax/7dTuKKiHsl/SS//btmaQu2Z0s6Mr+9KCI2liR7t1JjZjulbuJlLouIv7QRbju+nl+fZHuofmFE/EKpx4Ykfcj2P0p6XX7/jxFxdxvrWp1fd6t4P3cVxb3WjeJYo9RQrE3bytOVnle+UdJZnYe2hWJbH9Y01eaKz/GCfAyWuTK/Hlnsw/yd8Lw8/30R8UBJvvcr3U7QSrFt272vvTb9PS3Svljp3vsfRsRtxcyI2CDpqvy2X1fri/ru0uE4AgBQikY9AExBEbE+It4oaS9J/1epEX6L0lUkSTpEqRv0F2obSU4D4O2Z315l+89lk9IVOOXyy9zYYP5f8+sD2tR4r1c0QndqXMOOPaBNjfd6RWy/brC8aVy2n5gHbPtdHlQsvGkQvKKr8pyKcT5e6faFkPT9sgQRsUqbuiYf3KCcn1ZcXynb023/Ux4U7s48UFxRpxU52Qw13lfvVuppsIPS1V9L+mhEfKvNUL6j9APVwZIW2z7FdtVt2cgu+XVF2cJIg0gW236R7XNtP67sB4waT8qvv4qI5e0EY3u27bfZ/kkeCPChmm39lZysnTo/Ob++vsnnuOguvo2knfPfD1e6xUNqfOytVbVu8cW23aVpqi0V6dfmxnkzRYP9ypJlRRf8k/rUyK49lnZumAoA2jR90AEAAAYnIv6qdC/3xyXJ9u6Sni3p7UoN8hMl/VjSB3KW2itkVe7FLx35W6l7dZniivNfmvQCKNJU+ifc9s9V/uPCeyLiPXXzqqy3VexbxGX7TZL+XanRWqRdodQQldK9yDOU7quuotj2q6L5o9Zur0tf766K69uC7e0kLdKmxqGUBm+7S5uuYO+eX7dVyRXviBix/QqlbtxSui3gTe3GEhH/a/tVkj6kdL/yU3KMyyR9W9KlEXF94xJKbZ1f1zdJ8wql++kPlHRBntba/oHSD2Wfj4iHatIX2+NP7QRi+yBJ/68mv5R6ChSD5W2l9MNJ1eNH2vRZLsYAaKX4LNceS3c0SV/lR4viKn/pUyOaqLJvZHt/SYcq3e5zVUmSH0u6Vem2nuMkfbXNONpV26uh3ToDQENcqQcAjIqIv0TEfypd8SyuPP9jTZLa88aOEeEW07x+xd7ErkqNofppu36s3PajlK5IW6nR+ShJW0fE7Ih4WEQ8TGnANGlTo7+qrVsnaaqs235Vb1Nq0N+tdDV094jYJiJ2y3WaW5O2Wb3+oebvPZQGDWxbRHxKaVDDM5UaZ/cojUnwfyUtsX1Om0UWXdJnNVnnHyQ9RtJzJV2q9ESH4naHz0q6Nv/4UWh3/xY+rXTM/kLSsZK2j4gdImL3vK2LwQjbKb/4LB9f4XPsiFjWZsxVYil6cLTqQl+v2Dc7NrhFqFBcpR+WdE9t75jcw2FEqUFfm7ZWOz86FD96NHwMqDbvsdJunQGgIRr1AIAt5PuZi6tW+9Usqr3/+qD+RdS5iJjXoJFyXp9CeL7S+XZRRLw2In5bcg/87iX5mimusM9s8fSC4laJjq/IN1E0JF8bEZflXh+1WtbJ9t9r073lv1b6keJy21t1ElD+UeoDEXGC0o85T1Tqmm5JF9h+TBvFFT0Lmt7mEREPRcTCiHhlRByk9MPEWUoNwoO1adwASfpzft1HFeUR7Z+o9APMcyJiUUnvjHaPH2nTZ7ndz3HtsdSsu3+V++SLbdvO+Am16YeUnj6xhXzb0CltlPks2/Vd4ouGd5W6FGmaNdaL+q6PiNVN0gFAW2jUAwAauS+/jnZxjYg/alNj4Hlb5ECZomFd2v3b9rbadK91Vddr0/gHR5YlyI8efEJ++4s2y6+iab2URldvyPb2Slezp0n6lKSnKY1b8BilJyyMSSQ/V/rx4fa8nr9vo4il+XXfNtf753xbx/vzrKfWLP5Zfn2M7bmqZvSHmSb34Tfd1g0U4yk8v818f5C0Mv99eFmCfEzPr1DWvPz6u2aJStyiTbd4NNo/T1O67eYBpR8fdmoy/Urpav5JdWUUx/ZetvdUA7ZnSnpcXZ4y8/Lr0iZpAKBtNOoBYIqxvW9+znSzNNsoPWtakn5Zt3hBfn217QOblOHcsJzqVuXXRzdY/lY1uNrYSB6t/Hv57VsajPj+FqX79NcqPU++2xrWK3c5f2uL/B9UauT8UdKZEXGX0j3qkvTG/Oz5Sppd2c+9IorB1Nq5XeHH+fVvbc8oWedwi67fRTfs2nV+V+le8yFJF1eMo9jOu9verSSOR0s6uWJZtRbk1/m2X9Ysoe3R3gp5zIkv5bdn2i7bpq9T4/E0ah2SX39YIe2ofJW7GLCy0Y8HRXf6RRFxZ0SsbDRpU33qu+BfrTR2gSS9sUlIr1S67WJE6QkijXRUXwBohUY9AEw9j5K01PaXbb/Q9mjXUtvb2n620j+dxRWwD9Tl/zelq3XbSvq+7ZfX3jdsey/bpymNfv3cXlZkgrgmvz7T9jn5BxPZ3tX2xZLOVmf3175NqRFxsKTPF1cSbW+X7x//l5zu33rU1beo13ttP7Vo4No+RKnx2nBEc9vPU2pAjUh6WUSskaSI+LrSKPjTJF1me4eKsbzL9hdtn+D0uL9iPbvbvkTpWI6amKu4QelRedO16SpsrUdJ+rXtM23vV1P/YdvPl/SGnG5RkSGP1F40Dk+y/V+2D6iJdw/bp+WYCzcp9TSw0tMo/qZmPc/LdWo2WGKpiPi2pC/nt5+yfX7dd8FOto+3/VWlx1/WukjpCviBkhba3jfnmWn7TKUBA1epifxDQfHj4o/ajb8mzyH1C/L3UdGT6Mv1y0sUaebnQQkljf548O789kzb73V6Akixnt3yZ634geaTEfH7JuuhUQ+gN8oeXs/ExMTENHknSccoNXBqp/uVutTWzntI0jkNyvgbSb+tSbtRqWF6f10ZL6/LtzjPP7VBuUfk5cuaxF9ahqRT8/zFHWyTjtdbpQylK4HFNhlRGuhrJL//pNJV05B0Xl2+eUW+But8Zd72teU+VLOuyyUNleRblpcfMYbj6OFK91cX61qn1Lgsjqdn1CybV5PvYTX5LiopdztJv8/LF1TZxkpd3WuPu1VKDfLaeaXHcos6firnvaBk2ePqyn8gfwY21sz7uaQdSvK+oS7dGm3+2Vlcl/65delXS3ow/32r0r3jjY69hset0g9zX6mrx8q8/Wrnfbok74vrjrUVSj0iiiv5n1HJMV2Tv4j5+x0ef4fX1N91y/4hL1svaaeK5f0u53l33fxpkj5atz1Wa8vvy4WSZjYpf0+lz+hapYEOO/rcMTExMZVNXKkHgCkmIhZJ2l/p0WELtel58Nsp/aP6C6VG0mMj4l0NyvhfpWelv1qpG/i9Ss8af0jpCucHle4l/mxvajHhvEjpyvlNSg0fK3XvfnlE/FOnhUbEx5Wu/l2p9Ki97ZQaZNdIOjEiToktB+Xrikgjvz9R6YeDvyp1KV+p9OzvQyLi6gZZP6l0Ff+X2nwQuaLctZJeptQAenm+Gt3K+5S6fH9V0s1K23drSbdJ+oKkwxsdyy18Mr++qGTZTZJeIOljSvdRr1T6DKxWuor8WkmHRUkviYh4r9Ln59NKP7AMK/0ocINSz5jX16X/itI94tco/QAwrNSYfU8u53Z1ICLui4jnSnqW0tXq5UojvW+l9L1wZa7jq0vyfl7SYZL+O9d9K6Uf+s5UGscgWqz+xfn1k01TNY79B0r7em9t/lhFaVM3+u9FxIqKRRZd8E+xPVSznpGIeJXS99kVSreLTFe6teX2nO94Sc+NiGYj379Y6bj8r8g9UwCgWxzR6jsXAAAMQu5qfYvSaNljfXwdOmD710pd7Q+JiOsGHc9kkEeZv1PpqvWeEXF/h+W8Sanr+4cj4owuhth1tpco3Srz5Ij4aav0ANAOrtQDADB+FY8Mq39cHPqn6E3wpoFGMbm8Tqm3wb932qDPPqr0mMBTSx5HN27YPlKpQb+IBj2AXqBRDwDAOJQfk/Wa/PZ/BhnLVBYRX5J0raQTbe836Hgmuvwow9cqXam/pEXypiLiPkn/qjQ2wOtbJB+kc5VuRzh70IEAmJymDzoAAACwOdvfVbqHd0jp3vL6JxCgv16ldN/0HKX7uNG5fZQa8z8e41X6wick7awOngDQD3kk/h9I+nxEXD/oeABMTtxTDwDAOGP7l5L2k3SjpAsjPeoNAABgCzTqAQAAAACYoLinHgAAAACACYpGPQAAAAAAExSNegAAAAAAJiga9QAAAAAATFA06gEAAAAAmKBo1AMAAAAAMEH9f0zxEHwQsvvUAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1152x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The effect injected by `habitable_zone_water` is apparent in the dataset. But is this effect of high enough statistical significance to confirm the habitable zone hypothesis?"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Defining the hypothesis\n",
    "\n",
    "Next, we'll create a new Hypothesis object representing the hypothesis that planets within one region of space are more likely to have atmospheric water vapor than those outside of it. Since we do not know the planets' sizes, the only independent variable is the distance from the star modulated by the stellar luminosity, called `a_eff` (this is already calculated). Likewise, the dependent variable is the presence or absence of water vapor, called `has_H2O`, which is either 0 or 1. The relationship between these values is parameterized by `a_inner` (the inner edge of the habitable zone in AU), `delta_a` (the width in AU), `f_HZ` (the fraction of HZ planets with H2O), and `df_notHZ` (the fraction of non-HZ planets with H2O *divided by* `f_HZ`). By defining the four parameters in this way, we can easily avoid parameter combinations inconsistent with our hypothesis (such as the fraction of non-HZ planets with H2O being higher than `f_HZ` planets)."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "source": [
    "# Define the hypothesis in functional form\n",
    "def f(theta, X):\n",
    "    a_inner, delta_a, f_HZ, df_notHZ = theta\n",
    "    in_HZ = (X > a_inner) & (X < (a_inner + delta_a))\n",
    "    return in_HZ * f_HZ + (~in_HZ) * f_HZ*df_notHZ\n",
    "\n",
    "# Specify the names of the parameters (theta), features (X), and labels (Y)\n",
    "params = ('a_inner', 'delta_a', 'f_HZ', 'df_notHZ')\n",
    "features = ('a_eff',)\n",
    "labels = ('has_H2O',)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "In addition, we must consider the prior probability distribution of these parameters. Conservatively, we suppose the inner edge might extend far inward, or be slightly farther from the Sun than Earth, and that the HZ could be very narrow or very wide. Similarly, the fraction of planets with water vapor in the HZ and outside of it could span many orders of magnitude. We therefore choose to impose log-uniform prior distributions on these parameters across the following ranges:\n",
    "\n",
    "0.1 < `a_inner` < 2 AU\n",
    "\n",
    "0.01 < `delta_a` < 10 AU\n",
    "\n",
    "0 < `f_HZ` < 1\n",
    "\n",
    "0 < `f_notHZ` < 1\n",
    "\n",
    "After deciding on the bounds, we can initialize the Hypothesis object."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "source": [
    "bounds = np.array([[0.1, 2], [0.01, 10], [0.001, 1.0], [0.001, 1.0]])\n",
    "h_HZ = Hypothesis(f, bounds, params=params, features=features, labels=labels, log=(True, True, True, True))"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We also need to define the null hypothesis against which `h_HZ` is to be compared. The null hypothesis says that the fraction of planets with water vapor is independent of their orbits - this is a one parameter hypothesis where `f_H2O`, the fraction of planets with water vapor, is the only parameter. Again, the prior distribution on `f_H2O` should be broad (0 to 1) and log-uniform in shape. We can define this null hypothesis and attach it to `h_HZ` as follows:\n"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "source": [
    "def f_null(theta, X):\n",
    "    shape = (np.shape(X)[0], 1)\n",
    "    return np.full(shape, theta)\n",
    "bounds_null = np.array([[0.001, 1.0]])\n",
    "h_HZ.h_null = Hypothesis(f_null, bounds_null, params=('f_H2O',), features=features, labels=labels, log=(True,))"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "\n",
    "Now that the Hypothesis has been formed, we can formally test it using our simulated dataset. The `fit()` method will automatically extract the appropriate variables from the data and estimate the Bayesian evidence for both the hypothesis and null hypothesis using nested sampling. It will return the difference between the two i.e. the evidence in favor of the habitable zone hypothesis."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "source": [
    "results = h_HZ.fit(data)\n",
    "print(\"The evidence in favor of the hypothesis is: dlnZ = {:.1f} (corresponds to p = {:.1E})\".format(results['dlnZ'], np.exp(-results['dlnZ'])))"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "The evidence in favor of the hypothesis is: dlnZ = 13.1 (corresponds to p = 2.1E-06)\n"
     ]
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Generally speaking, a result where dlnZ > 3 is considered significant, so this was a succcessful test of the hypothesis.\n",
    "\n",
    "**Note for iPython/Jupyter users:** We will need to reload the `generator` and `h_HZ` objects here to replace the ones you have defined above. These will produce the same results, but the next code block uses the `multiprocessing` module which is incompatible with functions defined in iPython."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "source": [
    "# Reload `generator` and `h_HZ`\n",
    "generator = Generator('imaging')\n",
    "from bioverse.hypothesis import h_HZ"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "\n",
    "## Computing statistical power\n",
    "\n",
    "We have not yet tackled the most uncertain part of this analysis: namely, just how common are habitable worlds? The fraction of Earth-sized planets in the habitable zone with atmospheric water vapor could far smaller than the assumed 75%, in which case it might be impossible to test the habitable zone hypothesis. To quantify the importance of this assumption, we will need to repeat the previous analysis several times with different values of `f_water_habitable`.\n",
    "\n",
    "The `analysis` module enables this through its `test_hypothesis_grid()` function, which loops the planet simulation, survey simulation, and hypothesis test routines over a grid of input values. Let's use it to iterate over values of `f_water_habitable` ranging from 1% to 100%. For each value, we will repeat the analysis 20 times to average over Poisson noise.\n",
    "\n",
    "This may take a few minutes. To speed things up, we will run 8 processes in parallel (you may need to change this number for an older CPU)."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "source": [
    "f_water_habitable = np.logspace(-2, 0, 10)\n",
    "results = analysis.test_hypothesis_grid(h_HZ, generator, survey, f_water_habitable=f_water_habitable, t_total=10*365.25, processes=8, N=20)"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stderr",
     "text": [
      "100%|██████████| 200/200 [01:53<00:00,  1.76it/s]\n"
     ]
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now, let's plot the average Bayesian evidence for the hypothesis as a function of `f_water_habitable`."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "source": [
    "plt.plot(f_water_habitable, results['dlnZ'].mean(axis=-1), lw=5)\n",
    "plt.xlabel('Fraction of EECs with H2O', fontsize=20)\n",
    "plt.ylabel('$\\Delta$lnZ', fontsize=20)\n",
    "plt.axhline(3, lw=5, c='black', linestyle='dashed')\n",
    "plt.xscale('log')"
   ],
   "outputs": [
    {
     "output_type": "display_data",
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEbCAYAAAABNllnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAwE0lEQVR4nO3deZxddX3/8ddn9qwzWWYSkiEZEkhICLIkAgGBkAAiLkAjrbZaQSu12lpt+2vd0AS12tZSRGsVq2KxtbaggBU0EMImICQxAoGEkGRC9mWyJ7PP5/fHOZPcudvcO3P3eT8fj/s4c7/ne77ne+89cz/3+z3f8z3m7oiIiORCWb4rICIiQ4eCjoiI5IyCjoiI5IyCjoiI5IyCjoiI5ExFvitQiMaPH+9NTU35roaISFFZtWrVPnevT5ZHQSeOpqYmVq5cme9qiIgUFTPb0l8eda+JiEjOKOiIiEjOKOiIiEjOKOiIiEjOKOiIiEjOKOiIiAxxa3cc4tM/fYkV6/fQ3tWd1X1pyLSIyBD30Es7+fHzb/Dj599gRFU5l8+s530XTuXi08dnfF9q6YiIDHGPvLL7xN/HOrp56KVdbNx3LCv7UtARERnCtrQc47XdR2PSr5o1ISv7U9ARERnCIls5vd7UWMvE2pqs7E9BR0RkCFu2NjboXD07O60cUNARERmy9h/rYOWW/THpV82emLV9KuiIiAxRy1/dTY/3TZsydjgzJozM2j4VdEREhqh453Oumj0BM8vaPhV0RESGoNaObp7csDcmPZvnc0BBR0RkSHr69X20dfb0SRszvJK5U8dkdb95DzpmdpOZeT+PlOZlMLPmJGXsyvZrEREpFo+8EvuVuPDMCVSUZzcsFMI0OGuApQnWXQosBB5Oo7xDwB1x0mOvfhIRGYK6e5zlr+6JSb8qy11rUABBx93XEASeGGb2bPjnXWkUedDdlwyuViIipWv1GwdoOdbRJ626oozLZmR+rrVoee9eS8TM5gAXAduBX+S5OiIiJSPeqLVLzxjP8Krst0Py3tJJ4k/D5ffcPZ25tqvN7H3AFOAY8CLwZJpliIiUJHdPOFQ6Fwoy6JjZMOB9QA/w72luPhG4Jypts5nd7O5PZKJ+IiLFauPeo2yOmkHaLBhEkAuF2r32+0Ad8LC7b01jux8AiwgCzwjgbOA7QBPwsJmdk2hDM7vFzFaa2cq9e2PHrouIlIJfxZlrbe6UMdSPqs7J/gs16NwSLr+TzkbuvtTdH3P33e5+3N1fdvePALcDw4AlSba9y93nufu8+vr6AVdcRKSQ5bNrDQow6JjZbOBiYBvwUIaK/Xa4vCxD5YmIFJ09h9tYs/VgTPqQDjoMfABBMr0D0kdkqDwRkaLzaJxrc6bXj2BaffYm+IxWUEHHzGqA9xMMIPheBoueHy43ZbBMEZGisizOLARXn5W92xjEU1BBB7gRGAM8lGgAgZlVmtmZZjY9Kv0sMxsbJ/9U4Jvh0x9lusIiIsXgaHsXz7zeEpOey641KLwh070DCJLNQDAZeBXYQjAqrdeNwKfMbAWwGTgCTAfeDtQQnB/6WobrKyJSFJ58bS8d3X0n+KwfVc25jXU5rUfBBB0zmwW8hYEPIFgBzATOI+hOGwEcBJ4muG7nHnf3hFuLiJSweKPWrpzVQFlZ9u6dE0/BBB13fxXo99W7e3O8fOGFn7r4U0QkSmd3D8tfjQ06V2fxttSJFNo5HRERybAXNu/ncFtXn7ThVeXMnz4u53VR0BERKXHL4nStXT6jnprK8pzXRUFHRKSE5XuCz2gKOiIiJeyVnYfZfrC1T1p5mbHwzIa81EdBR0SkhMVr5VzQNJa64VV5qI2CjohISSukrjVQ0BERKVnbD7aydsfhmHQFHRERybhH1sbOtTbrlNGcOnZ4HmoTUNARESlRj8S5IDSfrRxQ0BERKUmHWjv5zab9MelXK+iIiEimPb5+D109faebnFRbw1mTRuepRgEFHRGRErRsbfyuNbPcTvAZTUFHRKTEtHd18/j62LuEXpWHCT6jKeiIiJSYZze2cKyju0/aqJoKLpwWc5/LnFPQEREpMfEuCF14ZgOV5fn/ys9/DUREJGN6egprgs9oCjoiIiXkxe2H2HOkvU9aZblx+Yz6PNWoLwUdEZES8sgrsbMQzJ8+nlE1lXmoTSwFHRGREhKvay3fF4RGUtARESkRzfuO8druozHphXI+BxR0RERKRrxWzjmNtUwYXZOH2sSnoCMiUiIKedRaLwUdEZES0HK0nZVb4kzweVb+ZyGIpKAjIlIClq/bQ9T8nkwdN5wzGkbmp0IJKOiIiJSAuF1rs/I/wWc0BR0RkSLX2tHNUxv2xqQX2vkcUNARESl6T7++j7bOnj5pY0dUMXfqmDzVKDEFHRGRIrdsbewsBAvPbKCiACb4jFZ4NRIRkZR19ziPrYt375zC61oDBR0RkaK2+o0DtBzr6JNWXVHGpWeMz1ONklPQEREpYvFGrV16Rj3DqyryUJv+KeiIiBQpd497PqeQJviMpqAjIlKkXt9zlOaW433SzGDhrIY81ah/CjoiIkVqWZyutblTxjB+ZHUeapMaBR0RkSIV9945ZxVu1xoo6IiIFKXdh9tYs/VgTPpVswtrgs9oCjoiIkXo0VdjWzmnN4zktPEj8lCb1CnoiIgUoWK4d048CjoiIkXmaHsXz7zeEpNeyEOlexVE0DGzZjPzBI/YQejJy2o0s++b2Q4zaw/LvsPMCm/mOxGRAXhi/V46uvtO8NkwqppzGuvyU6E0FNIlq4eAO+KkH021ADObDjwDNAAPAOuAC4C/BK4xs0vcPfbngYhIEXnkldjf4otmTaCsrLDunRNPIQWdg+6+ZJBlfIsg4Hzc3b/Rm2hmtwOfBL4MfGSQ+xARyZvO7p64E3wWQ9caFEj3WiaY2TTgaqAZ+Neo1V8AjgHvN7PCHtohIpLEC5v3c7itq0/aiKpy5k8fl6capaeQWjrVZvY+YApBgHgReNLdu1PcfmG4XObufTo73f2Imf2aIChdBCzPUJ1FRHIq3iwEl8+sp6ayPA+1SV8hBZ2JwD1RaZvN7GZ3fyKF7WeGy9cSrN9AEHRmoKAjIkXI3Yt2qHSvQule+wGwiCDwjADOBr4DNAEPm9k5KZRRGy4PJVjfm14Xb6WZ3WJmK81s5d69sfcaFxHJt1d2Hmb7wdY+aeVlxhUzC3eCz2gFEXTcfam7P+buu939uLu/7O4fAW4HhgFLMrCb3mEdnqAOd7n7PHefV19fn4HdiYhkVrxWzoWnjaVueFUeajMwBRF0kvh2uLwshby9LZnaBOtHR+UTESkqy9YWd9caFH7Q6R0XmMqIs/XhckaC9WeEy0TnfERECta2A8d5ZefhmHQFncyaHy43pZB3Rbi82sz6vC4zGwVcArQCz2WueiIiufFonK61WaeMpnHM8DzUZuDyHnTM7CwzGxsnfSrwzfDpjyLSK83szHD2gRPcfSOwjGDwwceiiltK0Fr6D3c/lsHqi4jkxCNxZpUulgtCIxXCkOkbgU+Z2QpgM3AEmA68HagBHgK+FpF/MvAqsIUgwET6KME0OHea2aIw34XAFQTdap/N2qsQEcmSQ8c7eW7T/pj0Yutag8IIOisIrrE5j6A7bQRwEHia4Lqde9w97oizaO6+0czmAbcB1wDXAjuBO4Gl7h77qYmIFLgV6/fQ3dP3a3By3TDOmjQ6wRaFK+9BJ7zwM5WLP3vzN3Ny+HO89VuBmwdfMxGRwhBvqPSVsxowK/wJPqPl/ZyOiIgk1t7VzePr40zweVZh35Y6EQUdEZEC9szGFo519J2CcnRNBRecFjP+qigkDTrhzdA+mWphZvYBM3ts8NUSERGI37V2xZkNVJYXZ5uhv3M6NwFuZhcBH3D3tn7yNwGXZ6BeIiJDXk+Px70+pxhHrfVKJVT2EAxrfsrMTslyfUREJPTi9kPsOdLeJ62qvIzLZxTv/JCpBJ1/Bu4F5gIvhEOSRUQky5atjb0t9fzp4xhVU5mH2mRGKkHnuLv/PsGtnicBT5jZe7JbLRERKfZ758ST8pkod78V+EOCa2T+08y+mLVaiYgMcc37jrFhz9GY9CETdADc/b+BBcBu4DNmdp+ZFddscyIiRSBeK+ecU+uYMLomD7XJnLTH3Ln788CbgTXADcDTZtaY4XqJiAxpy16JPZ9TjBN8RhvQQG933w68BbgPOBd43szmJ91IRERS0nK0nVVbDsSkF3vXGgxiRgJ3b3X3G4EvAROBx4C3ZapiIiJD1fJ1e4ia35Op44ZzRsPI/FQogwZ9Sau7f55ggIETdLuJiMggxDufc/XsCUU5wWe0/oLOFoLbDCQVDjC4nGCAgYiIDFBrRzdPbdgbk37V7OKc4DNa0mlw3P20VAty9xfMbAYwbtC1EhEZop7asJe2zp4+aWNHVDF36pg81SizMno/HXc/CsQOLBcRkZTE61pbeGYD5WXF37UGaZ7TMbPLzez/zGyPmXWaWXecR1e2KisiUsq6e5zl6+LcO6cERq31SrmlY2ZvB+4HyoE3gPWAAoyISIas2nKA/cc6+qTVVJZx6RnFO8FntHS615YAncDb3X1ZdqojIjJ0PRLngtC3nF7PsKryPNQmO9LpXpsD/EQBR0Qk89w9/lDps0qnaw3SCzpHgf3ZqoiIyFD2+p6jNLcc75NWZrDozIY81Sg70gk6ywFNdSMikgXL4rRy5k4dw7iR1XmoTfakE3T+DphuZp+zUrgsVkSkgMQLOqUw11q0dAYSfAFYCywFPmhma4g/W4G7+4cGXzURkaFh9+E2frf1YEx6qcxCECmdoHNTxN9N4SMeBxR0RERS9Oirsa2cMxpGctr4EXmoTXalE3RSnhJHRERSt2zt0OhagzSCjrtvyWZFRESGoqPtXTy7sSUmvVSDzqBvbSAiIgP3xPq9dHT3neCzYVQ15zTW5adCWaagIyKSR/FmIbhy9gTKSmSCz2gJu9fM7PsDLFOj10REUtDZ3cNjcSb4LNWuNUh+TuemAZap0WsiIil4fvN+Drf1nTd5RFU5F08v3duSJQs6Gq0mIpJF//fijpi0BTMbqK4onQk+oyU8p+PuW9J9AFuBc3NW+yzZsWMHZpaRR38ytZ8lS5Yk3c+SJUv0mvSa9JoK7DV9dfE5MWVHd60V02tKRUbuHGpmU4E/AW4GJmaqXBGRoWR4VTlXlNgEn9EGHBzMrBy4DrgFuBKw8OGZqZqIyNDysStOp3ZYZb6rkVVpBx0zm0bQqrkJmEAQaDYA9wAXAG/PYP1ERIaE0xtG8uFLp+W7GlmXUtAxswrgBoJWzRUE54IOAt8Ffujuz4b5vpudaoqIlLYvXT+HqorSv3TS3BP3hpnZGcCHgQ8A44EeYBnwQ+ABd2+Pyv9d4IPuXtRDL+bNm+crV67MdzVEpMQ8u7GF9373uZj0xec38s+/HzuooNiY2Sp3n5csT38tnfUE52heBv4B+E93j52ZTkREkuro6uHWB16OSa8dVslnrj0zDzXKj1TacoeAnwI/y0bAMbNxZvYnZvYzM3vdzFrN7JCZPW1mHzKzlNubZtZsZp7gETvXhIhIjnz3qU28vudoTPrfXXNmyd0dNJn+Wjq3Ah8kuIHb583sSYKutXvdPfbdG5gbgX8DdgIrgDcIBij8HvDvwNvM7EZP1g/Y1yHgjjjpmaqviEhatu4/zp3LN8Sknzeljve8+dQ81Ch/kgYdd/8y8GUzeyvBuZ13ApcB3zSznxIMIlg+yDq8BrwL+IW7n5hq1cw+AzwPLCYIQPelWN5Bd18yyDqJiGSEu/OFB9fS3tV3JunyMuPL159dshN7JpJS15W7/8rd3w2cCnwW2AW8D1hmZm+Y2ZfMbMZAKuDuj7n7zyMDTpi+C/h2+HTBQMoWEcm3X63dHXdSz5submL2pNF5qFF+pXWdjrvvAb4KfNXMFhEMob4O+AzwaaA1w/XrDJddSXP1VW1m7wOmAMeAF4En3b07w3UTEUnqWHsXS3++NiZ94ugaPnnVgH6nF70Bz0gQdqstN7PxBBeK/gkwgwzNSBBeG/TH4dNfprHpRIILVSNtNrOb3f2JTNRNRCQVdzz6GjsPtcWkf+GdsxlZPTRnCxv0lUjuvs/dv+buZwILgR8PvlpA0KKaAzzk7r9KcZsfAIsIAs8I4GzgO0AT8LCZJRwIb2a3mNlKM1u5d+/eQVVcROTVnYf5/q+bY9IXzKznmjkTc1+hApH04tB8MbOPA18H1gGXuPv+QZb3NeCvgfvd/Yb+8uviUBEZjJ4e593ffobVbxzsk15dUcYjn7ycKeOG56diWZbKxaEFN+eCmX2MIOC8Alwx2IAT6h2QcFkGyhIRSep/Vm6NCTgAf7Hw9JINOKkqqKBjZp8AvkkwA8IV4Qi2TOgdOjIiQ+WJiMTVcrSdr/5yXUz6tPoRfPiy0p/Qsz8FE3TM7O+AfwHWEASc2DGGAzc/XG7KYJkiIjG+8vA6Dh7vjEn/0vVzSvqOoKkqiKBjZrcSDBxYBSxy931J8laa2ZlmNj0q/SwzGxsn/1SC1hPAjzJYbRGRPp7fvJ97V22LSb/hvMlcPH18HmpUePI+Zs/MPgDcBnQDTwEfj3Pb02Z3vzv8ezLwKrCFYFRarxuBT5nZCmAzcASYTnB/nxrgIeBrWXkRIjLkdXT18Ln7X4pJH11TwWeunZWHGhWmvAcd4LRwWQ58IkGeJ4C7+ylnBTATOI+gO20EwT1/nia4bueeNOZvExFJy/ee3sxru2OnePzba86kftTQmdCzP3kPOuE8aUvSyN9McLfS6PQnCIKTiEhObd1/nK8vfy0m/ZxT6/jDC6bkoUaFqyDO6YiIFLOlP19LW2ffCT3LDL58/ZwhN6FnfxR0REQGYdnaXTz6auxg2w9c3MScybV5qFFhU9ARERmgY+1dLHkwdkLPCaOr+ashOqFnfxR0REQG6M7lG9gRZ0LPz7/jLEbVVOahRoVPQUdEZADW7zrC957eHJN+2Yx6rj176E7o2R8FHRGRNPX0OJ+7/yW6evpehVFVUcZt7zqLONcaSkhBR0QkTfeu2sYLzQdi0v/8itNpGq8pHpNR0BERScOBYx185eFXY9KnjR/Bn16uCT37o6AjIpKGrz68jgNxJvT8oib0TImCjohIilY27+cnK7fGpF937iQuOV0TeqZCQUdEJAWd3T189mcvx6SPqqngs2/XhJ6pUtAREUnBD369mfW7j8Sk/7+3zqRhVE0ealScFHRERPqx/WAr//LIhpj0NzXW8kcXTs1DjYqXgo6ISD+WPriW1s7uPmnBhJ5nU64JPdOioCMiksSjr+xm2Su7Y9Lff9FUzm7UhJ7pUtAREUngeEcXX4gzoWf9qGr++q0z81Cj4qegIyKSwDcee53tB1tj0m99x2xGa0LPAVHQERGJ47XdR/juk5ti0t9y+nje+aZT8lCj0qCgIyISxd353P0vx53Q84vXz9GEnoOgoCMiEuW+1dt5fvP+mPQ/u3w6p2lCz0FR0BERiXDgWAd//1DshJ5N44bzZwum56FGpUVBR0Qkwj/+ah37j3XEpN923RxqKjWh52Ap6IiIhFZt2c+Pn4+d0PMdbzqFy2bU56FGpUdBR0QE6Eo0oWd1Bbe+Y3YealSaFHRERIC7n2lm3a7YCT3/+uoZTBitCT0zRUFHRIa8HQdbuf2R12LS50wezfvnN+W+QiVMQUdEhrzbfv4Kxzv6TuhpmtAzKxR0RGRIW7FuD79cuysm/X0XTuWcU+tyX6ESp6AjIkNWa0c3n38wdvDA+JHV/I0m9MwKBR0RGbK+uWIDW/fHm9BzFrXDNKFnNijoiMiQ9PqeI9wVZ0LPS04fx7vOmZSHGg0NCjoiMuSs33WEv733RTq7oyb0LC/jtus0oWc2VeS7AiIiubD/WAcPrtnOvau38fL2w3HzfOTyaUyvH5njmg0tCjoiUrI6unp4fP0e7l21jRXr98S0bCJNGTucj15xeg5rNzQp6IhISXF31u44zL2rtvHg73bEnbwzntuuO0sTeuaAgo6IlIQ9R9p44Lc7uG/1trjT2SRSVV7GX189gwUzG7JYO+mloCMiRauts5vlr+7hvtXbeOK1vXT3JO4+i3b25FoWnz+Zd507mbEjqrJYS4mkoCMiRcXdWbP1IPet3saDa3ZwuK0r5W3rR1Xze+dNZvHcRmZMGJXFWkoiCjoiUhR2HmrlZ7/dzr2rtrFp77GUt6uqKOPq2RNYPLeRS08fT0W5rhTJp4IJOmbWCNwGXAOMA3YC9wNL3f1ArssRkfxr7ejmV2t3cd/qbTz9+j489d4zzp9Sx+K5jbzj7EnUDtfsAoWiIIKOmU0HngEagAeAdcAFwF8C15jZJe7ekqtyRCR/3J0Xmg9w36pt/OKlnRxtT737bFJtDTecP5nfO79R19sUqIIIOsC3CALFx939G72JZnY78Engy8BHcliOiOTY1v3H+enq7dy3ehtv7D+e8nbDKst525yJLJ7byPxp4yjTrQgKmnk67dVsVMBsGrARaAamu3tPxLpRBN1jBjS4e8KO3EyVAzBv3jxfuXLlQF+SiKToWHsXD720k3tXbeM3m/ente2Fp41l8dxGrj37FEZWF8rv56HNzFa5+7xkeQrhk1oYLpdFBgoAdz9iZr8GrgYuApbnoBwRyaKeHue5TS3cu3obD7+0i9bO7v43Cp06dhiLz29k8fmNnDp2eBZrKdlSCEGn96YVsfeKDWwgCBYzSB4sMlWOiGTBlpZj3LdqG/et3s72g7G3E0hkZHUF1549kcXnN/LmprHqPityhRB0asPloQTre9PrslmOmd0C3AIwZcqUfnYlIqk42t7FQy8G3WfPN6fefWYGl0wfz7vnNvLWsyYyrErT05SKQgg6/en9WTPYk09Jy3H3u4C7IDinM8h9iQxZJ7rPVm3j4ZfT6z6bVj+Cxec3csN5k5lUNyyLtZR8KYSg09sCqU2wfnRUvmyXIyIDMNDus9E1FbzznEksntvIeafW6V42Ja4Qgs76cDkjwfozwmWiczWZLkdEUjTQ7rMyg8tm1PPuuY1cOWuCZnceQgoh6KwIl1ebWVmcoc6XAK3AczkqR0SSGEz32RkNI3n33KD7rGF0TRZrKYUq70HH3Tea2TKCkWUfA74RsXopMAL4Tu+1NWZWCUwHOt1940DLEZH0DLT7rHZYJe86ZxLvntvImxpr1X02xOU96IQ+SjB9zZ1mtgh4FbgQuIKgO+yzEXknh+u3AE2DKEdE+jGY7rPLZ9Tz7rmnsmhWg7rP5ISCCDphK2UeJyfqvJZgBoE7CSbqTOloz1Q5IkOZus8kmwoi6AC4+1bg5hTyNXNy+POAyxGRvpr3HeO+1dv4qbrPJIsKJuiISO4dbe/iFy/u4N5V23ihOfU7f6j7TAZKQacEuTvuaLoQiUvdZ5JPCjoloLO7h5e2H+LZjS08t6mF375xkKPtXVSVl1FdWUZ1RTnVFWXUhH/3Lqsry6gJl8H6IF9knsi81XHyxtumstzUxZIFPT1Oe1cPbZ3dtHV109YZ/t0Z/t3VTXtnVHpv/jCttaObp1/fp+4zyRsFnSLU1d3DyzsO89ymFp7d2MLK5v0c64j9tdrR3UNHdw9HSP0mWJlQZpwIUGNHVDG5bhiNY4bTOGbYicfkuuE0jKou6dZYe1c3uw+1s/NQK7sOt7H3SHufANA3cPTQ3tUdd317GDw6unr632mGqPtMskVBpwh09ziv7DjMs5v28ezGFl5oPpDW3RRzrcehtbOb1s5uDh7vTHg/+6ryMibV1TB5zDAa64YHyzFBgJo8ZhgTR9dQXqBBqbWjOwgmh9rYeaiNXYfb+j4/1EbLsY58VzNtZzSM5MZ5jVx/rrrPJDsUdApQT4/zys6gJfPcphZ+s3k/R9oKN8gMVEd3D80tx2luOQ7E3kW8osw4pa7mREspWJ5sNU2sraGyvCzj9Tra3sWuQ63sjAggwfJk2qHWzozvN196u89unNfI2ZPVfSbZpaBTAHp6nHW7jgTdZZtaeH7z/pL6Uhuorh5n6/5Wtu5vBWIvsSozmDi65kQQmhzRddc4Zhin1NVQXXGyW8jdOdzaxc7DrXGDya7wcaSAW5GZou4zyRcFnTzo6XE27DnKsxv38WzYkjl4fHBBZnRNBRdOG8f8aeO4aNo4Zk4cRWd3D+3huYLeE9CRy/auk+cS2jsj1yfbJuLviGVbxLK7Jzd3huhx2HGojR2H2ni+OXa9GTSMqmZi7TCOtHay81BbWiO1ClFVRRk14QCO4BH+3TswpDf9RJ6TadUVZVRXllM3rJILTxur7jPJCwWdHHB3Xt9zlGfD7rLnNu1n/yD7+0fVVHDhaWO5KAwys04ZHXP+o7ysPPwFWzmofaWrqzsITsc7utl9uI1tB1rZduA42w60sv1g64nn2e4ydIfdh9vZfbg9q/tJxAzqR1ZzSm0NE2trmDC6hhHVFdRURAaD3hGAkX9HBZWKk0GjlAdeyNCgoJMF7s7GvcdOdJf9ZlML+44OLsiMrK7gzU1jmD99HPOnjWf2pNggUygqysuoKC9jRHUF9aOqmTM5/i2ODrV2sj0MQL3BaPuBVrYdDALUYFt/2VReZkwYVc3E2hpOqR0WLmsilsNoGFWdlXNOIsVMQScD3J3mluMnrpN5blMLe44M7tf18Kpy3tw0lvnTg5bMnEmjqSixL7DaYZXUDqtk9qTRcdcfbe+KH5TCVlO2RodVlhsTRp8MHqfU1jBxdGRQGcb4kVUl93mI5IKCTgZ88O4XWLF+76DKGFZZzrymMVw0bRzzp4/j7Mm1Q/5X8sjqCmZOHMXMiaPirj/e0cWOg61sPRGMervvgqC0N07gr64o6xM8TrRMRp98Pm5ElbqxRLJEQScDzpgwKu2gU11RxrymMcw/EWTqqKoY2kEmXcOrKji9YRSnN8QPSm2d3ew42MqeI+2MqqlgUu0w6oZXakiwSB4p6GTA/GnjuOvJTUnzVFWUMXfKmBPdZeecWttnOK9kXk1lOdPqRzKtfmS+qyIiIQWdDJjXNIbyMuszVLiqvIzzptSd6C4799Q6XQshIkOegk4GjKqp5PwpdQAnrpM5f+oYBRkRkSgKOhnyk1vm6+SziEg/dOY6QxRwRET6p6AjIiI5o6AjIiI5o6AjIiI5o6AjIiI5o6AjIiI5Y+65ufdJMTGzQ8CGJFlqgUMJ1o0H9mW8UtmX7DUV8r4GU1a626aaP5V8yfLo+CqcfZXi8dXf+sEcY1PdvT5pDnfXI+oB3DXQ9cDKfNc/G6+5UPc1mLLS3TbV/Knk6+cY0vFVIPsqxeOrv/XZPsbUvRbfzwe5vhjl8jVlcl+DKSvdbVPNn0q+ZHl0fBXOvkrx+EpnXxmn7rUMM7OV7j4v3/WQ0qTjS7It28eYWjqZd1e+KyAlTceXZFtWjzG1dEREJGfU0hERkZxR0BERkZxR0CkwZvZpM3vBzA6b2V4z+7mZzcl3vaQ0mNllZvagmW03Mzezm/JdJylOZvZRM9tsZm1mtsrMLk1lOwWdwrMA+BZwMbAQ6AIeNbOx+ayUlIyRwMvAXwKtea6LFCkz+wPg68DfA+cBzwAPm9mUfrfVQILCZmYjCa4cvt7dS/H6DckTMzsK/Lm7353vukhxMbPfAC+6+4cj0jYA97r7p5NtO+RbOmb2bjP7hpk9FXZpuZn9qJ9tGs3s+2a2w8zazazZzO4wszFZqOIogs/pQBbKliwrguNLSlA2jzszqwLmAsuiilhG0EOTlG5XDZ8DzgGOAtuAM5NlNrPpBE3JBuABYB1wAUF3xTVmdom7t2Swfl8H1gDPZrBMyZ1CP76kNGXzuBsPlAO7o4rZDVzZX8WGfEsH+CQwAxgN/FkK+b9F8MF83N2vd/dPuftC4F+AmcCXIzOb2ZfCXxnJHgvi7cjMbgfeAix29+4Bv0LJp4I9vqSkZfW4C0Wfm7E4aTF0TidC+M+5AvhPd39fnPXTgI1AMzDd3Xsi1o0CdhK88Q3ufixMH0/wyyCZN9z9eNS+/gV4D3CFu68b4EuSAlJIx1e4rc7pDAGZPu7C7rXjwHvd/X8j8v4rMMfdL09WH3WvpWdhuFwW+cEAuPsRM/s1cDVwEbA8TN9HmtOEm9nXCQLOAgWcISUnx5dIlLSOO3fvMLNVwFXA/0Zkvwq4r7+dqXstPTPD5WsJ1vfeg2fGQHcQ/lq4GXgvcMDMJoaPkQMtU4pGLo6vkWZ2rpmdS/D/PyV83u9QVylZAznubgduMrM/MbNZ4Q/lScC3+9uZgk56asNlopsf9abXDWIfHyUYsbacoFnb+/ibQZQpxSEXx9c84LfhYxiwNPz7tkGUKcUt7ePO3X8CfIJgwMIagnPP17r7lv52pu61zLJwOeATZe5u/eeSISoTx9fjEeWIpCLucefu3yIYgJAWtXTS0xvxaxOsHx2VTyQdOr4kH3J63CnopGd9uEzUp35GuEzUNyqSjI4vyYecHncKOulZES6vNrM+7104tPASgvmsnst1xaQk6PiSfMjpcaegkwZ330gw1UMT8LGo1UuBEcB/9F5DIZIOHV+SD7k+7ob8xaFmdj1wffh0IvBWYBPwVJi2z93/JiJ/9HQRrwIXAlcQND8v1jQl0kvHl+RDIR93CjpmS4AvJMmyxd2borY5lWCI6TXAOIIhzfcDS919f1YqKkVJx5fkQyEfd0M+6IiISO7onI6IiOSMgo6IiOSMgo6IiOSMgo6IiOSMgo6IiOSMgo6IiOSMgo6IiOSMgo6IiOSMgo7knJktMTMP791etMxstJndaWbNZtYVvqZz812vYjaQY8PMmsJt7s5axSRjFHSKTPjPlexxUwHU8aZCqUuW/SPwF8BLwFcIJkfclWwDM3s8hc9wSUT+phTyu5k1xdlXlZl9yMx+YWY7zazdzI6Y2Rozu8PM3pTRdyOLwtf4eI721ZzoPY3I0/s5LohIqzSzG8zse2b2spkdNrPjZvaSmd0WzticbL/zzOwHZrbJzFrD7V8ys38ys8kZe4F5pjuHFq+lCdLX5LISA/RN4L+BN/JdkUF6B/Cau79zANv+EGhOsO7xOGmHgDuSlHcw8omZzSCYN2sWsA94hOD9rgJmAx8BPm5m17v7gynXOvuK+diYDvwUOEZwu4BfACMJJtu8FfgDM7vE3fdFbmRmBnwV+Fugi+Cz+l+Cz+piglvVf9TMPuDu9+botWSNgk6Rcvcl+a7DQIX/dPv6zVj4JgFPDnDbu8NbR6fqYKqfuZlNAJYDjQSB6jPu3hqVp4FgQsgxadQh64r82DhCcGuAH0beBsDMqgiC0dsJ3vO/iNruVoKA0wy8w93XRq40s8XAj4D/NrOr3H0Fxczd9SiiB8F9yr2fPE1hvrsJ7gb4E2AP0AMsCPPMBb4O/A7YD7QBG4B/BsYkKfsPCL7QerdpBn4MzAvXP95bxziPpjDPkvD5gjjlLwJ+GVH+awS/Amvj5O3dVwXwmbD+7cBW4B+AqjTf21OAfw1fUwewl+DLYm6C/UY/Hk9hH73bxrz2fj7L5jRex3fDbf4rhbzVEX9PAL5GcCfJYwStp/XhcTQthbKeDd//EVHpT4b1+V5U+uww/T8i0vocG8BNSY6nJXGO9yaCltK+8PhZSfBFns5x0Bx5vGboc7w4zP9SnM+3Mzzezk6y/UfC7dcBZem8nkJ7qKVT2qYDvyH44v5PYBhwOFz3YeAG4AngUaAcOB/4K+BtZnahux/pLSjsAvgB8AGCf+ifEnwpNxLcc2M9wT/43QRfVtcR3JdjTUR9DiarrJn9KfBvBF94/0sQKBcAfwe8M+yaiFfGfwGXAg+Hr+9agl+ODcDNyfYZse/TgKcJWi+PEQTSU4Ebgbeb2WJ3/78w+90EXzpfALaEzyFxd1nOmNkw4P3h00RdsCe4e3u43XDg1wTHzCPAzwEDphJ8lvcS3I8lmeXARQSfxS8jyr0wXL8oKv/CiO0SWRO+juj3GmK7IacCz4f1vAcYS/Aj6QEzu9Lz20LoDJddUek3E/xo+h93fynJ9v9O0CKaCVzOybt9Fh0FnSIVebI5QrO73x3x/C3AV9z9M3HyfgX4mLt3R5X7IYID/KMErYVeHyYIOC8AV7n7oYhtygm+4HH3u4P4xHXA/VH1SfZ6pgJ3AkeBC9x9XcS6bwF/RnDi/pY4m08HzvLwnh9m9lmCFtwfm9mn3T3pyf3QtwkCzufc/ctR+34S+KGZTXX3o72vycy+QPCeL0nlNUa5KckIrW/HqXNdgs8cYJe7fzv8ex5QDWx39/Vp1GcRwft4h7t/MnJF2D1UnUIZjwGf5WRrFYIAVEUQyK4ys+ke3Kmyd5+928Xl7muANSm+1wsIWj8ngq2Z/VdYl/9H+l/UnzCzgwnWNaVZ1gfD5S+j0t8SLh9NtrG7d4UDKf6Q4PbRRRt08t7U0iO9B4m7Gk5073Cyu2EXEd0nKZZvBCetH4tKfyks87wUyrgpzHtTgvVLiOqaIPiycuDv4+QfQ9CCaaVvd9Dj4TZXxtlmabiu364VgtaaE/ySroyz/p5w/R/H+SweT/P97a1zsse5EfmbUsi/JiL/74dpz6VZr3cmev/TKKMm/IxWR6T9I8Gv/Plh+beE6WXAAYKBGEmPjf7e64j3qBkoj7N+C8GdMlN9Hc0pvOcer54JynsXQdf2VqK6roFXwnKuSaGcr4Z5vzXQz6gQHhoyXaTc3eI8FkRl+52H3SfRwuGdf25mT5vZfjPrNjMn+OcYDUyOyDsCmAPsdvffZuklnR8uY371uvsB4LcEX2pnxtl2ZZy0reEylRPl54XLp9y9M876x6LyZcIVCT5D8+DXfbQtSfKfG5HPwqWnWZ8ngO3Ap8zsl2b2cTObG7ZiU+LubQS3PD7XzMaFyQuBF9z9WWA3J1s35wN1JGnlDMAaj2q5h7YysAETpyV6zwner36Z2cUE3b/HgMXhsdwnS7hM5fMa6GdbUBR0SluybqWfAN8gOHn+AMEv0qXh4xB9u1PqwuX2zFfxhNpwuTPB+t70uugVHv88T2/feSpfmgPedwHaES4b09nI3Q8TnI/5AScHmawEdpnZUjOrTLGo5QRfjleY2RiCQN17zuYxYGF4fnBRRP5MOZggvYs8fNeZ2XyC84w9BC2Z5+Nk6z22pqRQZO9nmug4LQo6p1Pa4v4iMrN5BIMIHgWujfx1b2ZlBCfhIx0Ml9m8QK33HNFEYG2c9adE5cvWvuPJ5r4zbSXBCLJGM5vpaZzXcfdtwIfCoDCboJXyMeDzBF/at6ZQTG/L5UqCL9sy+gad9wLnEAQdp5jPTSRhZpcSXKfTA7zV3Z9LkPVpgoE4VxKMOkxUXjnBOSsIBnwULbV0hqbTw+WDcbqTLiAY5XaCB9ccvAxMMLNUuph6uzhS7poh6D6Dk/9YJ5hZHXAuwRDYV9MoM919v8XM4v0QuyJcrs7CvjPKg+tx7gmf9hskzCxmgIAH1rr7N4CrwuTrU6zCCwTn3xYRBK1WgqHUcDL4XEtwMvxFj7pQMoke0jue8sbMFhK0cLoIBt0kCjgQjMbrBm4ws7OS5PsgwUCX9aTYtVeoFHSGpuZwuSAyMbxg8F8TbHNnuPyOmdVGrjCzMjM7JSKpJVym0mXQ60cEJ5z/wsxOj1r3RYLzTD9KdI5qMMJf+I8QnJD+ROQ6M7uQYMTQAeBnmd53lnwO2Ab8UTiFyrDoDGY23szuBN4TPp+TYNqXCeHyeCo7Ds+pPEnww+ZG4Onez8zdNxMce38JDCe98zktBEPYC5qZXQ38H8EPpEXu/kKy/O6+Cfh7oBJ40MxmxynzeoLuzm7go+7ek+l655K614amFwia6L9nZs8QNPEnAG8j+CW1I842/04wvPOPgQ1m9gDBdTqTCH7Rfp9g5BEEv2yPEww5HUtwAhngGx4x1DqSuzeb2ScIgt5qM/ufsPzLCUY+rSO4XidbPkLwnvxT+MWxkpPX6fQAN3vEdUsZkGzI9Bp3vz8qLdmQaQhmOGgGcPfdZraIYBqcvwE+YGaR0+DMIvjBUc3JFsyVwO3h8bCO4BqpRoKh7z3AP6X6wghaNO8gGEYffc5mOfChiL/TKfM9ZvZzYBVBK+JJdx/ojBAZZ2YzCc6P1gAPAdeZ2XXR+Tx22PcSYATBNXK/M7NfEXQxVxJcVHohQYvxve6eyYEX+ZHv4XN6pPcgzRkJkuQZC3yL4JdnG7CR4BfX8DCtOcF2f0TQvD8UbreZ4MLT86PyXUMQfI5ycnhpU7huCYlnJLgaWEbQsmgHXicY5FAXJ+/jid4L+hm2nWCbyQQXp24huEJ8H8EX95uTfBaPp/n5PR7xfiR63B2RvymF/IneyyqCL/iHCE4+dxBM1fISQcv17Ii8s4DbCYLt3vC9bya4KPTiNF/j2RH1enPUuveG6Z3AqDjbxj02CALYfxH8gOkmwYwEyd7zNOrfHHm89vM5LohIW5DKZ5WkzAsI5uTbTBBkjhJ0a38NaEznMyjkh4UvVkREJOt0TkdERHJGQUdERHJGQUdERHJGQUdERHJGQUdERHJGQUdERHJGQUdERHJGQUdERHJGQUdERHLm/wOUJSWiWJd+jwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "A more useful metric than the average Bayesian evidence is the survey's \"statistical power\", which is the likelihood that the survey could successfully test the hypothesis under a certain set of assumptions. Run this cell to plot the statistical power versus the fraction of EECs with water vapor."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "source": [
    "power = analysis.compute_statistical_power(results, method='dlnZ', threshold=3)\n",
    "plt.plot(f_water_habitable, power, lw=5)\n",
    "plt.xlabel('Fraction of EECs with H2O', fontsize=20)\n",
    "plt.ylabel('Statistical power', fontsize=20)\n",
    "plt.axhline(0.8, lw=5, c='black', linestyle='dashed')\n",
    "plt.xscale('log')"
   ],
   "outputs": [
    {
     "output_type": "display_data",
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaoAAAEbCAYAAACLGcAmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA5dUlEQVR4nO3deZxcVZn/8c+TfSEr2SAhNAlZQUGMLGHLQgK4gaIOqGw6oIKAOvpzHxNn3B1FUERHEQUc3GYARTE7ASRAEBDJSkgCISEr2ddOP78/zq2kuupWdd3qqq6q7u/79arXTZ977r2num76qbPcc8zdERERqVbtKl0AERGRfBSoRESkqilQiYhIVVOgEhGRqqZAJSIiVa1DpQvQWvTr18/r6uoqXQwRkZry9NNPb3L3/vnyKFCVSF1dHQsXLqx0MUREaoqZrW4qj5r+RESkqilQiYhIVVOgEhGRqqZAJSIiVU2BSkREqppG/YmISGK/W/gKr27dw+TRAznh6J60a2dlu5YClYiIJHbPEy/z7CtbuXnWcgb06MzEUQP42ITh1PXrXvJrqelPREQS2bRzH8+t2Xro5w079vGbha9gZapUKVCJiEgi85ZuJHMpw+H9u3PskaWvTYEClYiIJDRnyfqstEmjB5TtegpUIiJSsP31DcxftikrfdLogWW7pgKViIgUbOGqLezcV98orUeXDoyr61O2aypQiYhIwWYv2ZCVds7I/nRsX75wokAlIiIFmxMTqCaXsX8KFKhERKRAL23cycpNuxqlmcGEUQpUIiJSBeJqU286pjd9u3cq63UVqEREpCCxzX5jyjfaL0WBSkREmrRj7wGeXLklK72cz0+lKFCJiEiTHlm+ifqGxtNRHNWrC6MH9Sj7tRWoRESkSbMXZzf7TRo9ACvXBH9pFKhERCSvhgZn3tK4/qnyN/uBApWIiDThuTVb2bxrf6O0zh3accawfi1yfQUqERHJK26035nH96Nrp/Ytcn0FKhERySsuULXEaL8UBSoREcnptW17eWHt9qz0idUYqMzsCjM7v5yFERGR6hJXmxo9qAeDe3dtsTIkqVHdAVxQroKIiEj1iZ+NouVqU5AsUL2WML+IiNSwvQcO8tiLcYskVm+gegiYaGYKViIibcCClzaz58DBRml9unXk5GPKt0hinCRB54tAD+DnZtYyg+dFRKRi4pr9JowaQPt25Z+NIl2HBHn/B9gGXAFcamarCM2BnpHP3X1yaYonIiKV4O45p01qaUkC1YS0f3cGRkWvTJmBS0REaszyDTt5deueRmnt2xnnjOzf4mUpOFC5u/qmRETaiLja1Lhj+9Cra8cWL4uCj4iIZJmzZH1WWksPS09RoBIRkUa27t7P06tfz0qfNLr8q/nGSRSozKydmd1gZgvMbJuZ1afte5OZ3WZmI0tfTBERaSkPL9tIxhqJDO3bjeH9u1ekPEmmUOoEzARuBoYDO4D0MYorgQ8BHyhh+UREpIVVcpHEOElqVJ8BJgLTgYHAz9J3uvtWYD6g+QBFRGpU/cEGHl62MSu9Uv1TkCxQfQB4zN2/6u4NxA9DXwkMLUnJaszatWsxs5K8mlKq60ybNi3vdaZNm6b3pPek99TG3lPHDu35x7TzWf2tt7P10XsA6NapPace17cs76kQSQLVccCCJvJsAbLfjYiI1KyzR/Sjc4eWWSQxTpJAtQfo3USeocDWYgsjIiLVZ3KFRvulJAlUzwJTLQyqyGJmvQj9U0+WoFwiIlIlJoxu+dko0iUJVP8NHAPcY2Y903eYWW/gTqAPcHupCiciIpX1xiG9GNCjS0XLYO6FT81nZj8HrgYOAK8D/YFngBMI8//9yN1vKEM5q964ceN84cKFlS6GiEjRrrzjyawRf584bwSfOK98j8ea2dPuPi5fnkQP/Lr7hwnPSi0iBCkDTgFeBD7cnCBlZkPM7A4zW2tm+8xslZndbGYFLXxiZleZmTfxOphxTF0T+e8t9v2IiNSSXfvqeXzF5qz0SsyWninJ7OkAuPudwJ1m1pXQ1LfN3Xc1pxBmNhz4GzAAuB9YApwK3ARcYGZnunv2b7CxZwnPeMU5G5gE/CXH/ueA+2LS/9nENUVEWoXHXtzE/oMNjdL69+jMiUf3qlCJDkscqFLcfQ9hJGAp3EYIUje6+62pRDP7HvBJ4GvAR5soz7OEYJXFzB6P/vnTHIc/6+7TEpVYRKQVmbs0ZjaKUQNo18KLJMZJMoXSPDP7kpmdaWZFB7iY8w4DpgKrgB9l7P4KsAu43MyKmmTKzE4ETgdeBR4svqQiIq1TrkUSJ1ZBsx8kq1GdCZxDmJFil5k9AswB5rj7M80ow6RoOyOa8eIQd99hZo8RAtnpwOwizv+RaPtzdz+YI8/RZvYR4EhgM/C4u/+jiGuJiNScF9ZuZ8OOfY3SOrVvx1kj+lWoRI0lCVR9gHOB8wjB5QLgQsDN7HVgHiGQzHH3pQnOm1oleFmO/csJgWokCQNV1I/2QaCBjLkJM0yJXunHzgOudPeXk1xTRKTWzFmSXZs6bVhfjuhcssazZkmywu9OQtPZgwBm1g+YTAhak4B3RS9Pcl4g1VO3Lcf+VHrvBOdMeV903IPu/krM/t3AfxAGUrwUpb0RmAZMBGab2cm5BouY2bXAtQBDh7bJKQ5FpBWYHROoqmG0X0rRCye6+ybgz4TA9RdgI2G4eql73lLnK/yBr8OujbY/idvp7hvc/d/d/e/uvjV6zSfU4J4Ajgf+NdfJ3f2n7j7O3cf171/ZJ7dFRIqxccc+/rFma1Z6NQWqRPW6aPqkMwk1qcnAm4H2wH7gceBWkvcjpWpMucZA9szIV2hZxwLjgTWEgFowd683s58BpxH65X6Q5HgRkVoxb+kGMud9GN6/O8ceWZlFEuMUHKjMbCbhD38XQp/P08B3CAMqHnP3vUWWIdWflevR5xHRNlcfVi6FDKLIJ/V4dvV8WiIiJRbXPzV5TGUnoc2UpEY1mdD8NhP4d3cv1eSzc6PtVDNrlz7yz8x6EGpwe2h6iZFDzKwLcDkhoP68yHKdHm1fyptLRKRG7a9v4JHlm7LSq6nZD5L1Ud1NeBZpKvC4ma0xszvN7HIzG1JsAdx9BTADqAOuz9g9nVCj+VVqQIOZdTSz0dFsFrm8lzBK8c85BlEQneu0uNngzWwS4UFjCO9bRKTVeWrVFnbuq2+U1qNLB958bEEz17WYJKP+rgAws5EcHqL+NuAKwhD15Rwenv6HhOW4jjCF0i1mNhlYTOgfmkho8vtiWt7B0f7VhOAWJzWIItdMFCnfAk6IhqKvidLeyOFnu77s7n8r+F2IiNSQuId8zx3Zn47tix5nVxbFzPW3jBA8brOwjvDJhGeqbiJMc/SRpOd19xVmNg74anSutwLrgFuA6e6+pdBzmdkY4CwKG0RxF2FI/VsIz4R1BNYDvwV+6O6PJHkfIiK1JG7apMljqqvZD5ox15+ZnUTotzqPMOlratBBMcPIiZrori4g3yryDIF398X59mfk/TnF92GJiNSslzbuZOWmxo+ImsG5I2s4UEV9Qqlh6RMJ0w2lAsIiwui/2YQZKkREpIrFjfY7ZWgf+naPXcS9opLUqJYTaktGmED2AQ73Sa0vfdFERKRc4gJVtY32S0kSqH4DzCIEppVlKo+IiJTZ9r0HeHJldtd/zQcqd7+snAUREZGW8ciyTdQ3NB5OcHSvLowe1KNCJcqvqMEUZtYRGE2Y8HUbsNjdD5SwXCIiUiaxzX5jBhAGclefRIPlzaynmd0ObCWspjsPeAbYama3m1nvEpdPRERK6GCDMy9uNd8qbfaDZKP+egKPAScAO4BHCM86HUV4lupa4CwzG+/u20tfVBERaa7n1mxl8679jdK6dGzH+OHVsUhinCQ1qs8TgtSPgWPdfYK7X+buE4BjCcvIj43yiYhIFZob0+x35vB+dOnYvgKlKUySQPVuYIG7X+/uW9N3uPs2d7+BsNTHJSUsn4iIlFDctEkTq7jZD5IFqqE0/TDvw8AxRZdGRETKZt22PSxal90zU839U5AsUO0Gmno3/aN8IiJSZeYu2ZiVNnpQD47u3bUCpSlckkD1FPBeMxsRtzOaYul9UT4REakyc5ZkTyJUjZPQZkryHNV3COtGPWVmtxIWPFwHDAImADcARwDfLXEZRUSkmfYeOMijL8Ytklhdq/nGSTIzxWwzuw74AfCF6JViwAHg4+4+q7RFFBGR5nr8pc3sPdDQKK1v906cfEzvyhQogaTrRv3EzP5CWOb9TUAvwswUzwB3u/vq0hdRRESaa07MaL8JI/vTvl11zkaRrpiFE18GvlaGsoiISBm4e85pk2pBda03LCIiJbds/U5e3bqnUVr7dsbZI/pXqETJJA5UZvYBM5ttZlvMrD7azjazD5SjgCIi0jyzY0b7vaWuD726dqxAaZJLMtdfR+D3wNsJgyfqgY2ElX4nAhPM7H3AezSTuohI9YibNmlyDYz2S0k61987gCcIgamLux8FdAEmAU8SgthnS11IEREpzuu79vP06tez0qt92qR0SQLVFcCLwAR3f9jdGwDcvcHd5xGepXoJuKrEZRQRkSI9vGwjGWskcuyR3Rjev3tlClSEJIFqCHC/u++P2+nu+4D7gcGlKJiIiDRf7Gi/0dW7SGKcJIFqLdBUz1vHKJ+IiFRY/cGGmlskMU6SQPVr4D3RAopZotV93wPcU4JyiYhIMz29+nW2761vlNa9U3tOPa5vhUpUnCSB6qvAQuBJM3u/mQ0xs47R9gPAAsKAiv8oR0FFRCSZOTG1qbNH9Kdzh+pdJDFOkpkpUk+LGXBXzH4DRgB7M9o+3d0Tz4AhIiLNEzdtUq01+0GyQPUI4E3mEhGRintly26Wb9iZlT5hdG3MRpEuyezpE8pYDhERKaG40X4nDenFgB5dKlCa5tFcfyIirdDsmEBVSw/5plOgEhFpZXbtq2fBis1Z6bU0bVI6BSoRkVbmsRc3sf9g40USB/TozAlHxz5dVPUUqEREWpm4/qmJowbQrgYWSYyjQCUi0orU+iKJcRSoRERakRfWbmfDjn2N0jq1b8dZx/erUImaT4FKRKQVmR3zkO9pw/rSvXPtzrugQCUi0orMiVnNd3KNDktPyRlizWxosSd195eLPVZERIqzccc+nluzLSt9Uo0OS0/JVxdcRXFTJnkT541lZkMIE99eQFjefh1wHzDd3bOXp4w/xyrg2By717v7oBzHjQe+BJxOWLH4ReAO4FZ3P1j4uxARqZy5MZPQHj/gCIYe2a0CpSmdfAHlV7TQ3H5mNhz4GzCAsPjiEuBU4CbgAjM7092zn16Ltw24OSY9e9KrcO2LgD8Ae4HfAFuAdwDfB84E3lvwGxERqaC5MaP9ar3ZD/IEKne/qgXLcRshSN3o7remEs3se8Anga8BHy3wXFvdfVohGaO1tf4bOAhMcPeFUfqXgTmE9bcudfd7C30jIiKVsL++gfnLNmal1+q0SekqPpjCzIYBUwlNjT/K2P0VYBdwuZl1L8Pl3wP0B+5NBSkAd99LaAoE+FgZrisiUlJPrtzCrv2Neyp6dunAm4/tU6ESlU41jFecFG1nuHujOT/cfYeZPUYIZKcDsws4X2cz+yAwlBDk/gHMz9HXlLr2QzH75gO7gfFm1tnd98XkERGpCnEP+Z47agAd21e8PtJsxQx6eAtwPjAY6ByTxd39wwlOOSraLsuxfzkhUI2ksEA1iOyFHVea2dXu/nCh13b3ejNbCZwADAMWF3BtEZEW5+7MjhmWPqkG156KU3CgsrBs753ABwmr+Xq0TfG09CSBqle0zR5T2Ti9dwHn+gVhgccXgB2EAPNx4FrgL2Z2hrs/V6prm9m10bkZOrTo0fwiIs3y0qZdrN68u1FaO4NzR9Z+/xQk66P6OHA5obYyjhCUbgbGA18gBIZ7CcGhlFLBsMkRiO4+3d3nuPt6d9/t7v90948C3wO6AtNKeW13/6m7j3P3cf37t45vLiJSe+JG+50ytA99u3eqQGlKL0mguhJY6u5Xufvfo7St7r7A3b8JTAQu4XC/T6FStZZeOfb3zMhXjNuj7TkVuLaISFnFTZvUGkb7pSQJVKMIQ7bTHWo6dPdngD8B1yUsw9JoOzLH/hHRNlcfViFSn2LmyMGc1zazDsBxQD3wUjOuLSJSNtv3HuCpVVuy0ifX8GzpmZIEKqNxzWIX0Dcjz3JgdMIyzI22U82sUXnMrAfhods9wIKE5013RrTNDDipwHtBzDHnAN2Av2nEn4hUq0eWbaK+oXHvxODeXRk1sEeFSlR6SQLVq4SRfikvAW/OyDOCEMAK5u4rgBlAHXB9xu7phFrQr9x9F4CZdTSz0dFsFoeY2Qlmlhk4MbNjgR9GP96dsfv3wCbgUjMbl3ZMF+A/ox9/nOT9iIi0pLjRfhNH9yeMf2sdkgxPf5LGgekvwGeiWRz+F5gAXERo/kvqOsIUSreY2WTCUPDTCP1ey4AvpuUdHO1fTQhuKe8FPmdmc4GVhMEdw4G3Eebv+zPw3fSLuvt2M7uGELDmmdm9hCmU3klo6vw9YVolEZGqc7DBmbc0ezaKyTU+CW2mJDWqPwDtzey46OdvE4LFdMJDtbcCW4HPJS1EVKsaRxj+fhrwb4QgcwtwRoHz/M0F/o/Qr/R+4FPAucCjhIEgb3f3/THXvi/KN58wGOQG4EB0/KXu3iLzHYqIJPXcmq1s2dX4z1qXju04Y/iRFSpReRRco4r+oN+X9vMWM3sTcA0hqKwiNNGtK6Yg7v4KcHUB+VbR+PmtVPrDQOYDvYVe+zHgrcUcKyJSKXNiRvudObwfXTq2r0BpyqdZUyi5+zYymtNERKRlzI55fmpSKxrtl1L7k0CJiLRB67btYfG67VnpE0e14UBlZh81sxVmdnSO/YOj/UmmTxIRkSLETUI75qieHN27awVKU15JalTvB9a5+9q4ne7+KrCGMBegiIiUUVz/VGtYJDFO0pkpnmsizz9I/sCviIgksPfAQR5bsSkrvTVNm5QuSaDqRRh+ns92oPZX6RIRqWKPr9jM3gONlu+jb/dOnHxM78oUqMySBKp1wBubyPNGIPvpMxERKZm42SgmjOpP+3atZzaKdEkC1VzgAjM7K26nmZ0NXEhhixuKiEgR3J25S7LrA5NaabMfJAtU3wL2A7PM7HtmNjWaX2+qmX0fmAnsi/KJiEgZLF2/g1e37mmU1qGdcfaI1rsmXpKZKZaa2fuAXwOfAG5K222E/qn3u7uWbBcRKZO4YelvqetLr64dK1CalpFoZgp3f9DMhgFXEebk600YYLEA+GWBc/KJiEiR4hZJbM3NflDEFEpRMPqvMpRFRETyWLhqC0+vfj0rvTVOm5ROUyiJiNQAd+cbf1mSlT6sf3eG9ctcvLx1yVmjMrNzon8+6e57035ukrvPb3bJRETkkBmL1sfWpj527vBWtUhinHxNf/MAB8YQFi9M/VyI1jXHvIhIBdUfbODbD2XXpkYN7MG7TxlSgRK1rHyB6quEwLQp42cREWlBv3t6DSs27spK/+yFo1rtQ77pcgYqd5+W72cRESm/3fvr+f7MZVnppx3Xt1Uu6REnyTIfQ82sZxN5epjZ0OYXS0REAO54dCUbduzLSv/chaNbfd9USpJRfytp/JBvnBujfCIi0kxbdu3n9odfykp/6xsG8aahbWf+7ySByqKXiIi0gFvnLGfnvvpGaR3aGZ85v22tplTq56gGAtk9fiIiksjLm3dz94LVWemXnTqU41r5c1OZ8s5MYWZXZCSdHJMGYTj6UOBy4PkSlU1EpM36r5lLOXCw8UDrbp3ac+PkERUqUeU0NYXSnRweku7ARdErU6pJcDcwvSQlExFpo/756jbuf3ZtVvo1Zw+jf4/OFShRZTUVqK6OtgbcAdwH3B+T7yCwGXjc3beWqnAiIm3RN2OmSup3RCeuOWdYBUpTeXkDlbv/MvVvM7sSuM/df1X2UomItFGPLN/Ioy9uykq/cfIIjuiceB7xViHJelQTy1kQEZG2rqHBY2tTdUd247JT2+4jqkke+O1jZmPNrHNG+tVmdr+Z/drMTi19EUVE2oYHnlvLC2u3Z6V/5vzRdGzfdhe7SFKP/DrwQeDQnB1mdgNwM4cHU1xsZuPcfVHJSigi0gbsqz/Id2cszUo/aUgv3vqGQRUoUfVIEqLPBGa7+560tE8DrwLnAO+L0j5VorKJiLQZdy94mTWv78lK/9yFY9rMVEm5JKlRDQZmp34ws7HAMcBn3f3RKO29hKAlIiIF2r73AD+cszwrfeKo/pwx/MgKlKi6JKlRdQX2pv18JuHZqllpaSsIAU1ERAp0+7wVvL77QKM0M/h/F7StqZJySRKoXgXSf2vnA9uB59LS+gDZdVcREYn12ra93PFY9lze737TEMYclXfBijYjSdPfXOBKM/s4oWb1TuAP7t6Qlud44JUSlk9EpFW7edYy9h5oaJTWqUM7PjV1ZIVKVH2S1Ki+AewEfgD8lBCspqV2mtkA4FzgbyUsn4hIq7V8/Q5+uzD7u/1V4+sY3LtrBUpUnZI88LvSzE4A3hMlPeDuL6dlORb4EfDrEpZPRKTV+vZfl9LQeN5ZenbpwHUThlemQFUq0Xwc7v4a8MMc+54CnipFoUREWruFq7Ywc9H6rPTrJx5P726dKlCi6tV2H3UWEakQd+frf16clX50ry5cOb6u5QtU5XLWqMzs3wnDz3/k7luinwvh7v4fSQtiZkOArwIXAEcC6wiztU9399cLOP5I4F3A24A3EIbJ7yesj/UL4BcZAz8wszoge7jNYb9x90uTvhcRkXxmLFrP31/empX+ySkj6dKxfcsXqMrla/qbRghUvwG2kDZwogkOJApUZjacMAhjAGEZkSXAqcBNwAVmdqa7b27iNO8FfkwIcHOBlwkrDr8b+BlwoZm919095tjnCEEx0z+TvA8RkabUH2zg2w9lTzw7elAP3n3KkAqUqPrlC1Sp2dJfzvi5HG4jBKkb3f3WVKKZfQ/4JPA14KNNnGMZYcj8g+k1JzP7AvAkcAkhaP0h5thn3X1ac96AiEghfrtwDSs27spK/+wFo2nfrm1PlZRLzkDl7g/n+7lUzGwYMBVYRRg1mO4rwLXA5Wb2b+6e/ekeLt+cHOmvmdnthGA3gfhAJSJSdrv313PzrGVZ6acd15cJo/pXoES1IckyH1eY2RubyHOimV2RsAyTou2MzD4kd98BPAZ0A05PeN50qblJ6nPsP9rMPmJmX4i2ed+niEgx7nh0JRt27MtK//xbNfFsPklG/d0JXNxEnosIAxeSGBVts79mBKmZGot6TNvMOgCp4PlQjmxTgFSt63bgOTOba2Ztd6UyESmpzTv3cfvDL2Wlv+0NR3HyMb1bvkA1pNTD09sTBlMk0SvabsuxP5Xeu5gCAd8ETgT+7O5/zdi3mzDw482EeQr7EGbXmEtoJpxtZt1zndjMrjWzhWa2cOPGjUUWT0TaglvnvMjOfY0bdTq0Mz59/qgcR0hKqQPVSKDJoeQJperDSQMgZnYj8G+EUYSXZ+539w3u/u/u/nd33xq95hP6zJ4gzF34r7nO7+4/dfdx7j6uf3+1L4tIvJc37+aeJ1ZnpV926lCO65fzu7BE8s5MYWZ3ZCRdHD17lKk9MBQ4G3gwYRlSNaZeOfb3zMhXEDO7njAv4SJgsrtvKfRYd683s58BpxHW1/pBkmuLiKT77oylHDjY+Lt2t07tuXHyiAqVqLY0NYXSVWn/duDk6BXHCbWQTyYsQ2rt5Vx9UKlPMlcfVhYz+wTwfcJzUJPdfUPCMgGk2vL0dUdEivb8mm088NzarPRrzh5G/x6dK1Ci2tNUoDou2hrwEnAz8bWLg8Dr+YaP5zE32k41s3YZz0D1ICzQuAdYUMjJzOyzhH6pZ4Ep7r6piDLB4VGG2b2fIiIF+lbMw739jujENecMq0BpalPeQOXuhxpVzWw6MDc9rRTcfYWZzSD0C10P3Jq2ezqhRvOTVBA0s47AcOCAu69IP5eZfZkwDdPTwNSmmvvM7DTgGXffn5E+icM1w7uLfW8i0rbNX7aRR1/M/q580+QRHNE50ZzgbVqSZT6ml7Ec1xGmULrFzCYDiwn9QxMJTX5fTMs7ONq/GqhLJZrZlYQgdRB4BLgx5rmEVe5+Z9rP3wJOMLN5wJoo7Y0cfrbry+6u9bVEJLGGBuebf8muTdUd2Y1LT9WTL0kUFdLNrD3QD4htYM1Yp6pJUa1qHIcnpX0rYc6+WwiT0hYyECLVTNke+ESOPA8TngdLuYswke1bgAuBjsB64LfAD939kSTvQ0Qk5YHn1rJo3fas9M+cP5qO7bVwRRIWP0drjsxmbyD0/0wkR5AizJ7e5uq048aN84ULF1a6GCJSBfbVH2Tyfz3Mmtf3NEo/6Zje3HfdeM1CkcbMnnb3cfnyFBxQzGw0h5eZnwm8gzDr+HrgFEINKzVruYhIm3XX46uzghTA5y4YrSBVhCT1zy8TmsbGu/tFUdr/ufsFhGa3XwBjgULXrRIRaXW27z3AD+e+mJU+cVR/zhh+ZAVKVPuSBKoJwJ/c/fm0NAOIRuR9hDArReJFE0VEWovb561g6+4DjdLM4LMXjq5QiWpfkkDVj8MTxEKYibxb6gd3ryc0/U0tTdFERGrLa9v2csdj2YuGX3LKEEYP6hlzhBQiSaDaAhyR9vMmwrRJ6faTeyokEZFW7eZZy9h7oNFqRXTq0I5PTilq8QeJJAlUK0h7bonwUO0UMxsAEM0yfhGQ/XVCRKSVW75+B79d+EpW+tXj6xjcu2sFStR6JAlUM4CJacte3A70BZ4xs98BzwPHAj8rbRFFRKrftx5aSkPG0z49u3TgugnHV6ZArUiSQPXfwIeBrgDu/iDhwdquwCXAAMJMD7eUtogiItXtqVVbmLV4fVb69ROPp1e3jhUoUeuSZAqldcBvMtJuMbMfEQZabPAkTw+LiLQC7s43/rw4K/3oXl24cnxdyxeoFSq4RmVmQ80sa9iKux909/Xu7mbWQ8u3i0hb8tcX1vP3l7dmpX9q6ii6dGzf8gVqhZI0/a0Ebmoiz41oMIWItBH1Bxv49l+zJ54dPagH73rT4AqUqHVKEqiMw8vCi4i0eb9duIaXNmYvw/fZC0bTvp3+XJZKqafwHQgUs3iiiEhN2b2/nu/Pyl54/PRhfZkwqn8FStR65R1MYWZXZCSdHJMGYWmNocDlhGHqIiKt2s8fWcnGHfuy0j934RhNPFtiTY36uxNIjeRzwgO9F8XkS30quwmr8oqItFqbd+7jJ/Nfykp/2xuO4uRjerd8gVq5pgLV1dHWgDuA+4D7Y/IdBDYDj7v71lIVTkSkGt0650V27qtvlNahnfHp80dVqEStW95A5e6/TP07Wur9Pnf/VdlLJSJSpV7evJt7nlidlX7ZqUM5rl/3mCOkuZI88DuxnAUREakF352xlAMHG89t0L1Te26cPKJCJWr9mrVkvJm9E5hEaBqc7+5/KEmpRESq0PNrtvHAc2uz0q85Zxj9e3SuQInahrzD083sHWY238zOjdn3C+D/CA/53gD81swUqESkVXJ3vvlQ9lRJ/Y7oxDVnD6tAidqOpmpU7wROAZ5ITzSztwNXEp6Z+j6wA7gWuNjMLnP3/ylDWUVEWtzLm3czY9FrzFi0nidXbsnaf9PkEXTv3KzGKWlCU7/dUwkj+fZmpH+IMFz9anf/PYCZ3UVYs+oDgAKViNSkhgbnH69uY+ai15i5aD3L1u/Mmfe4ft259FRNb1puTQWqQcDjMennAFuBQ0197v6amT0InFmy0omItIC9Bw7y+IrNzFi0ntmL17Mh5kHeOJ85fxQd25d6gh/J1FSg6kNYgv6QaHb0vsAfY5b1WEloLhQRqWqv79rPnCUbmLV4PQ8v28ju/QcTHX/m8Udy4YmDylQ6SddUoNoBDMlIe3O0fSbHMZnNhCIiVWH15l3MXLSemYvWs3D16xzMXJK3AMf07cpFJw3mI+cO01RJLaSpQPU88DYzO8LdUw217yL0Tz0ak/84YF0JyyciUrQk/U35nDSkF+eNGciUEwYyamAPBagW1lSgugf4CfCwmf0SGEkYLPEaMDc9o4VP7izi+7RERFpEsf1N6Tq1b8cZw49kytiBnDdmIIN6dSlDSaVQTQWqnwPvBs4HTiY82HsAuMndMxt0JxMGX8wqcRlFRPJK9TfNXLSe+cuT9zcB9OrakUmjBzBl7EDOGdmfIzTkvGo0Nddfg5m9DbgMGE+YePZ/3f3ZmOz9gB8AD5S6kCIimVL9TTMWrWfhqi0U0d3EMX27MmXMIM4bO4C31PXVCL4q1eRXBndvIDQB3tNEvnuBe0tULhGRRhoanOfWbGXW4vXN6m9645BeTFF/U01R3VYk8sqW3cxavJ65SzeyY+8BBvfuypSxA5kwagC9unasdPHapL0HDvK3FZuYuWiD+pvaMAUqabPcnX++up2Z0fQ4S17b0Wj/My9v5U//WEeHdsapx/VlytiBTBk7kCF9ulWoxG1DKfubzhszkHNHqb+p1ln2M7tSjHHjxvnChQsrXQxpwr76gyx4aQszF73GrEUbeG178sf+xhzVkyljBjBl7CBOHNxTTUclUIr+piF9uh76MqH+ptphZk+7+7h8efQ1Q1q9bbsPMHdp+Ib+8LKNWSuzJrV43XYWr9vOLXNe5KheXThvzEDOGzuQM4YdSacO+uNYiFR/U+rh2+Ub1N8kuSlQSav0ypbdzFy0nlmL1/PEyi1FzUBQiHXb9nLXgtXctWA1R3TuwLmj+jN17EAmjBxAr27q10qn/iYplgKVtAruzvOvbmNW1HyU2d9UqBMH92TKmEGMOaoHj724iZmL1rN2W2HNgzv31fPgP9bxYEa/1nljBnJM37bZr7UlNZ+enm+SZlAfVYmoj6rllaK/qWN74/Rhh7+hH927a6P97s4La7cfGhL9wtrtRZV19KAeTB07sE30a63atItZi/V8kxSmkD6qqglUZjYE+CpwAXAkYc7A+4Dp7v56Oc9jZuOBLwGnA12AF4E7gFtjZuCIpUDVMkrR39SjSwcmjgrf0M8d1Z+eXQpvont16x5mRU2Kj6/YTH0Rf4VbW79WqfqbNJ9e21QzgcrMhgN/AwYA9wNLCIs2TgSWAme6++ZynMfMLiKsq7UX+A1hWZN3AKOA37v7ewt5DwpU5VOK/qbUM1FTxg7k1ONK8w19+94DzFu6kZmL1jNvyQZ2FBE0a7Vf63B/03pmLd7ARvU3SZFqKVD9FZgK3Ojut6alfw/4JPATd/9oqc9jZj0JtadehCC2MErvAswBzgAui2bdyEuBqnRK3d903tgBjD2qvM1t++sbeHLllkOzdBfar5Wu2vu11N8k5VATgcrMhhGWsF8FDI+mbErt60FoujNggLvvKuV5zOxDhIl3f+XuV2acbxIwG5jv7uc29T6SBqq5SzfwzzXbCs7fVry2fS+zF5env6mllLJf69xR/TmiU2X/mB842MCClVua3d80ZexAxtX1UX+TNFIrz1FNirYz0oMLgLvvMLPHCLWk0wmBo5TnSR3zUMz55gO7gfFm1tndk7dt5DFn8QbuWrC6lKdsk5rT31QuZsaJg3tx4uBefOK8kUX3ay15bUfRtclKO2lIr/CFYaz6m6T5qiFQjYq2y3LsX04IMCPJH6iKOU/OY9y93sxWAicAw4DFea4tLagc/U3lNLh3V64cX8eV4+tK0q9VjTq1b8f4448Mg0TU3yQlVg2Bqle0zdUOlkrvXYbzNOvaZnYtcC3A0KFDmyieNEdL9jeVU88uHXnnSUfzzpOOLkm/ViWpv0laSi3cWam/SM3tTCvmPHmPcfefAj+F0EdVfNEkU7X0N5VTpw7tOGtEP84a0Y9p7zyhJP1a5Zbe3/SWuj50qPLarLQO1RCoUrWWXjn298zIV8rzlOraiZ07sr+WjojRrp0xcuARnDOyOvqbWkpcv9bDSzeyduueShcNgN7dOnLWiH7qb5KKqIZAtTTajsyxf0S0zdX31JzzLAXGRcc8nZ7ZzDoAxwH1wEtNXDux86KOZpE4g3t35f2nqTlZBKAa6u1zo+1UM2tUnmhY+ZnAHmBBGc4zJ9peEHO+c4BuwN9KPeJPREQKV/FA5e4rgBlAHXB9xu7pQHfCc06pZ586mtnoaBaKos8T+T2wCbjUzA6N448e+P3P6McfF/3mRESk2Sr+wC/ETn20GDiNMPXRMmB8auojM6sDVgKr3b2u2POkHXMxIWDtBe4lTKH0TqIplID3eQG/JM1MISKSXCEP/Fa8RgWHakPjgDsJgeXfgOHALcAZhczzV+x53P0+4FzCA76XADcAB4BPAZcWEqRERKR8qqJG1RqoRiUiklzN1KhERERyUY2qRMxsG2Gaplx6kft5rH6EQR21Jt97quZrFXuupMclyd9U3ubs1/3VstdqqfsryTGF5MuXp5z317Hu3j9vDnfXqwQv4KfF7gcWVrr85XjP1XqtYs+V9Lgk+Ztz/zS1X/dXy16rpe6vJMcUkq+Je6ii95ea/krnj83cX4ta8j2V8lrFnivpcUnyN/f+0f1VPddqqfsryTGF5MuXp6L3l5r+qoCZLfQmOhNFiqX7S8qpJe4v1aiqw08rXQBp1XR/STmV/f5SjUpERKqaalQiIlLVFKhERKSqKVDVODP7vJk9ZWbbzWyjmf3RzE6sdLmk9TCzc8zsATN71czczK6qdJmkNpnZdWa20sz2mtnTZnZ2IccpUNW+CcBtwHhgEmH9rFlm1reShZJW5Qjgn8BNhKVyRBIzs38BfgB8HXgTYQLxv5hZkwuvaTBFK2NmRxCeIL/Y3VvjszVSQWa2E/i4u99Z6bJIbTGzJ4B/uPs1aWnLgd+7++fzHasaVUJm9h4zu9XMHoma29zM7m7imCFmdoeZrTWzfWa2ysxuNrM+ZShiD8Ln+noZzi0toAbuMWmFynnfmVkn4M2ENQPTzSC0BuVVDUvR15ovAScBO4E1wOh8mWPWyFoCnEpoRrnAzM70ApcxKdAPgGeBx0t4TmlZ1X6PSetUzvuuH9AeWJ9xmvXAeU0VTDWq5D4JjAR6Ah8rIP9thA/yRne/2N0/5+6TgO8TFmf8WnpmM/vP6JtMvteEuAuZ2feAs4BL3P1g0e9QKq1q7zFp1cp630Uy+5osJi2L+qiaIfrPPBe4x90/GLN/GLACWAUMd/eGtH09gHWED2qAu++K0vsRvn3k87K778641veBS4GJ7r6kyLckVaaa7rHoWPVRtQGlvu+ipr/dwGXu/ru0vD8CTnT3c/OVR01/5TUp2s5I/yAB3H2HmT0GTAVOB2ZH6ZtIOGW+mf2AEKQmKEi1OS1yj4lkSHTfuft+M3samAL8Li37FOAPTV1MTX/lNSraLsuxP7V+1chiLxB9I7kauAx43cwGRa8jij2n1JSWuMeOMLOTzexkwt+ModHPTQ4rllarmPvue8BVZvavZjYm+oJ9NHB7UxdToCqvXtE214JjqfTezbjGdYSRfrMJ1e3U69PNOKfUjpa4x8YBz0SvrsD06N9fbcY5pbYlvu/c/TfAJwiDNp4l9Ke/1d1XN3UxNf1VlkXbojsK3d2aziVtWCnusXlp5xEpROx95+63EQZhJKIaVXmlvlX0yrG/Z0Y+kaR0j0kltOh9p0BVXkujba7+gRHRNlc7r0hTdI9JJbTofadAVV5zo+1UM2v0u46GcJ5JmDttQUsXTFoN3WNSCS163ylQlZG7ryBMEVIHXJ+xezrQHfhV6vkWkaR0j0kltPR9pwd+EzKzi4GLox8HAecDLwGPRGmb3P3TafkzpxlZDJwGTCRUi8drehtJp3tMKqGa7zsFqoTMbBrwlTxZVrt7XcYxxxCG8l4AHEkYPn4fMN3dt5SloFKzdI9JJVTzfadAJSIiVU19VCIiUtUUqEREpKopUImISFVToBIRkaqmQCUiIlVNgUpERKqaApWIiFQ1BSoREalqClRSE8xsmpm5mU2odFmaw8x6mtktZrbKzOqj93RypctVy4q5N8ysLjrmzrIVTEpGgaoNiP5D5ntdVQVlvKpaylJm3wZuAJ4HvkGYwPO1fAeY2bwCPsNpafnrCsjvZlYXc61OZvZhM3vQzNaZ2T4z22Fmz5rZzWb2xpL+Nsooeo/zWuhaq3L9TtPypD7HCWlpHc3sXWb2czP7p5ltN7PdZva8mX01mok833XHmdkvzOwlM9sTHf+8mX3HzAaX7A1WmFb4bVum50h/tiULUaQfAvcCL1e6IM30dmCZu7+jiGN/CazKsW9eTNo24OY859ua/oOZjSTM0zYG2ATMJPy+OwFjgY8CN5rZxe7+QMGlLr9avjeGA/8L7CIsnfEgcARhQtgvA/9iZme6+6b0g8zMgG8C/w+oJ3xWvyN8VuOBTwPXmdmV7v77FnovZaNA1Ya4+7RKl6FY0X/UTU1mrH5HA/OLPPbOaFn4Qm0t9DM3s4HAbGAIIbh9wd33ZOQZQJi0tE+CMpRdjd8bOwjLZPwyfUkMM+tECGBvI/zOb8g47suEILUKeLu7v5C+08wuAe4G7jWzKe4+l1rm7nq18hfg4aPOm6cuyncnYdXO3wAbgAZgQpTnzcAPgOeALcBeYDnwX0CfPOf+F8IfwdQxq4D/AcZF++elyhjzqovyTIt+nhBz/snAQ2nnX0b4ttkrJm/qWh2AL0Tl3we8AnwL6JTwd3sU8KPoPe0HNhL+wLw5x3UzX/MKuEbq2Kz33sRnuSrB+/jv6JhfF5C3c9q/BwLfJaz4uotQS1sa3UfDCjjX49Hvv3tG+vyoPD/PSB8bpf8qLa3RvQFcled+mhZzv9cRamSbovtnIeGPf5L7YFX6/Vqiz3F8lP/5mM/3QHS/vSHP8R+Njl8CtEvyfqrtpRqVZBoOPEH4Y38P0BXYHu27BngX8DAwC2gPnAJ8CrjQzE5z9x2pE0XNE78AriT8Efhfwh/yIYQ1a5YS/ijcSfgDdxFhXZtn08qzNV9hzewjwI8JfyR/RwiuE4DPAu+Imk3izvFr4GzgL9H7eyvhG+oA4Op810y79nHAo4Ra0hxC8D0GeC/wNjO7xN3/FGW/k/CH6ivA6uhnyN2U12LMrCtwefRjrubhQ9x9X3RcN+Axwj0zE/gjYMCxhM/y94T1jPKZDZxO+CweSjvvadH+yRn5J6Udl8uz0fvI/F1DdhPpscCTUTnvAvoSvljdb2bneWVrIgeibX1G+tWEL1q/dffn8xz/M0LNaxRwLodX5a05ClRtSHqHe5pV7n5n2s9nAd9w9y/E5P0GcL27H8w474cJ/ymuI9RKUq4hBKmngCnuvi3tmPaEoIC73xliGhcB92WUJ9/7ORa4BdgJnOruS9L23QZ8jDB44dqYw4cDJ3i0Zo6ZfZFQU7zCzD7v7nkHOERuJwSpL7n71zKuPR/4pZkd6+47U+/JzL5C+J1PK+Q9Zrgqz8i222PK3DvHZw7wmrvfHv17HNAZeNXdlyYoz2TC7/Fmd/9k+o6o6apzAeeYA3yRw7ViCEGrEyH4TTGz4R5WlE1dM3VcLHd/Fni2wN/1BEIt61CANrNfR2X5DMn/uH/CzLbm2FeX8FwfirYPZaSfFW1n5TvY3eujwSTvJywNX7OBquJVOr3K/yJ3M8ihpicON4W8RlrTToHnN0LH/ZyM9Oejc76pgHNcFeW9Ksf+aWQ0mxD+wDnw9Zj8fQg1pT00bqqaFx1zXswx06N9TTb7EGqFTvjG3jFm/13R/itiPot5CX+/qTLne52clr+ugPzPpuV/X5S2IGG53pHr95/gHF2iz+jvaWnfJtQmzojOf22U3g54nTAYJe+90dTvOu13tApoH7N/NWFF20Lfx6oCfuceV84c53snodn9FTKa1YFF0XkuKOA834zy3lbsZ1QNLw1Pb0Pc3WJeEzKyPedR006maCjtx83sUTPbYmYHzcwJ/6F6AoPT8nYHTgTWu/szZXpLp0TbrG/X7v468AzhD+HomGMXxqS9Em0LGSzwpmj7iLsfiNk/JyNfKUzM8Rmah1pEptV58p+cls+irScsz8PAq8DnzOwhM7vRzN4c1ZYL4u57CcuZn2xmR0bJk4Cn3P1xYD2Ha1GnAL3JU5sqwrOe0UIQeYXiBo0cl+t3Tvh9NcnMxhOapncBl0T3cqMs0baQz6vYz7aqKFBJpnxNXr8BbiUMILif8M13evTaRuOmnt7R9tXSF/GQXtF2XY79qfTemTs8vt8q1RdQyB/aoq9dhdZG2yFJDnL37YT+pV9weKDNQuA1M5tuZh0LPNVswh/UiWbWhxDcU31Qc4BJUX/n5LT8pbI1R3o9Ffj7aGZnEPpNGwg1pidjsqXuraEFnDL1mea6T2uC+qgkU+w3LzMbRxhIMQt4a3otwszaEQYipNsabcv50GGqz2sQ8ELM/qMy8pXr2nHKee1SW0gYeTfEzEZ5gn4qd18DfDgKJGMJtaHrgX8n/KH/cgGnSdWQziP8gW5H40B1GXASIVA5tdzXkoeZnU14jqoBON/dF+TI+ihhMNJ5hNGauc7XntAHB2HQS81SjUoKdXy0fSCmqetUwujAQzw8E/JPYKCZFdL8lWp+KbjZiNC0B4f/Mx5iZr2BkwnDjRcnOGfSa59lZnFf+CZG27+X4dol5eF5qbuiH5sMLGaWNUjCgxfc/VZgSpR8cYFFeIrQnziZEOj2EIatw+GA9VbCgIB/eMbDr3k0kOx+qhgzm0SoSdUTBh7lClIQRjEeBN5lZifkyfchwmCfpRTY7FitFKikUKui7YT0xOgh0B/lOOaWaPsTM+uVvsPM2pnZUWlJm6NtIc0ZKXcTOt1vMLPjM/b9B6Hf7O5cfW7NEdUkZhI65T+Rvs/MTiOMtHod+L9SX7tMvgSsAT4QTb/TNTODmfUzs1uAS6OfT8wxZdDAaLu7kAtHfUTzCV+G3gs8mvrM3H0l4d67CehGsv6pzYTHBaqamU0F/kT4UjXZ3Z/Kl9/dXwK+DnQEHjCzsTHnvJjQFHsQuM7dG0pd7pakpj8p1FOE5oN3m9nfCM0PA4ELCd/Y1sYc8zPCUNorgOVmdj/hOaqjCd+c7yCM2ILwDXo3YXhvX0InOsCtnjasPZ27rzKzTxAC5d/N7LfR+c8ljBhbQnieqlw+SvidfCf6Y7OQw89RNQBXe9pzZSWQb3j6s+5+X0ZavuHpEGa6WAXg7uvNbDJhCqVPA1eaWfoUSmMIX1I6c7imdB7wveh+WEJ4hm0I4TGDBuA7hb4xQs3p7YRHFjL7oGYDH077d5JzXmpmfwSeJtRW5rt7sTODlJyZjSL093YB/gxcZGYXZebz7CH204DuhGcYnzOzvxKavzsSHhQ+jVAzvczdSzn4pDIqPexQr/K/SDgzRZ48fYHbCN9w9wIrCN/sukVpq3Ic9wFC08O26LiVhIeJT8nIdwEhYO3k8FDeumjfNHLPTDEVmEGowewDXiQM9Ogdk3dert8FTQyRz3HMYMIDx6sJMwVsIvyxf0uez2Jews9vXtrvI9frzrT8dQXkz/W77EQICn8mdMDvJ0zz8zyhhvyGtLxjgO8RAvTG6He/ivCg7/iE7/ENaeV6S8a+y6L0A0CPmGNj7w1C0Ps14UvPQXLMTJHvd56g/KvS79cmPscJaWkTCvms8pzzVMIckCsJgWknocn9u8CQJJ9BNb8serMiIiJVSX1UIiJS1RSoRESkqilQiYhIVVOgEhGRqqZAJSIiVU2BSkREqpoClYiIVDUFKhERqWoKVCIiUtX+P71ZmRHcqLVzAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can see that a LUVOIR-like survey will be able to determine the existence of the habitable zone - but only if habitable planets are relatively common."
   ],
   "metadata": {}
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "9805f964417e1fb5bf47cad3ce563e418ccf3a95b9919f79d65dd2cee73654e5"
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3.8.8 64-bit ('base': conda)"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}