unique_otus.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from __future__ import division\n",
    "\n",
    "from biom import load_table\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "from IPython.display import display"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# read in data\n",
    "full_meta = pd.read_table(\"data/797_prep_296_qiime_20150818-225202.txt\", index_col=0)\n",
    "results = pd.read_table(\"data/full_results.txt\", index_col=0)\n",
    "# sync meta with results\n",
    "meta = full_meta.loc[results.columns[:-1]]\n",
    "# read in standard deviations\n",
    "results_std = pd.read_table(\"data/full_results_std.txt\", index_col=0)\n",
    "results_var = np.square(results_std)\n",
    "# sort indicies so we are all in the same order\n",
    "meta = meta.sort_index()\n",
    "results = results.sort_index()\n",
    "results_std = results_std.sort_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1159, 503)\n",
      "502702.0\n"
     ]
    }
   ],
   "source": [
    "table = load_table(\"data/153_otu_table.s100.d1000.zzyvl.biom\")\n",
    "table = table.filter(results.columns[:-1], inplace=False)\n",
    "print table.shape\n",
    "print table.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqoAAAHcCAYAAAAeFogrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3X1QnPW9//9XgO6CMVTkLpGGamgn3CQCIXDGJtiZaG1r\nk6hd6HjbEm84NRB7fq2dEnuqNtaTaLwnN00w1RrHHiOrtannRE1b70hiCRBJxbRu1AhpAkuMgeqS\nS8j1++NMtq5Uv1zLwn6A52OGYbg+1177vl7DJq+wF1cm2bZtCwAAADBMTLQHAAAAAP4ViioAAACM\nRFEFAACAkSiqAAAAMBJFFQAAAEaiqAIAAMBIFFUAAAAYiaIKAAAAI1FUAQAAYCTHRfW9997TDTfc\noOLiYn3961/XU089FVzr6OjQkiVLVFhYqIULF6qhoSHksTt27NCiRYtUUFCgiooKtbe3D/8MAAAA\nMC45LqpLly5VV1eXNm/erJtuukmrVq3S9u3bg2tpaWnyer1avHixqqurdfjwYUnSoUOHVFVVJY/H\nI6/Xq6SkJFVVVUX2bAAAADBuTLJt2x7qzn/5y19UXl6u7du3KyMjQ5JUV1enP/zhD/qP//gPLV26\nVDt37pTb7ZYkLVmyREVFRaqurtb999+vpqYmPfLII5Kkvr4+zZs3T7/85S9VXFw8AqcGAACAsczR\nT1Tb29t1+umnB0uqJM2cOVN/+ctftHv3buXl5QVLqiQVFRVpz549kqTW1taQQhofH6/c3Fy1tLQM\n9xwAAAAwDjkqqikpKerp6dHx48eD2w4dOqT+/n4dOXJEaWlpIfsnJyers7NTktTV1TVoPSUlJbgO\nAAAAfJyjopqfn6/U1FStWLFCgUBABw4c0MMPP6xJkybp+PHjcrlcIfu7XC5ZliXp/97q/6x1AAAA\n4OMcFVWXy6UHHnhAr776qoqKinTllVfq0ksvlW3biomJGVQ6LctSfHy8JMntdn/m+lA4uJwWAAAA\nY1yc0wfMmjVL27dv15EjR5SUlKSXX35Zp59+ujIzM/XKK6+E7Nvd3a3U1FRJUnp6uvx+/6D1nJyc\nIT/3pEmT1NMT0MDACadjT0ixsTFKTEwgM4fIzTkyCw+5OUdm4SE358gsPCdzixRHRfXYsWO6/vrr\ntX79eiUnJ0uSXnjhBZWUlOjss8/Whg0bZFlW8C3+pqYmzZ07V9L/XTbQ3NwcPFYgEFBbW5uWLVvm\naOCBgRPq7+cbxgkyCw+5OUdm4SE358gsPOTmHJlFl6O3/j//+c8rEAho9erVam9v1xNPPKGnnnpK\n1113nUpKSnTGGWeopqZGPp9PGzdu1N69e1VWViZJ8ng8am5uVl1dnXw+n5YvX67MzEyVlJSMyIkB\nAABgbHN8w/97771XBw4c0OLFi/XII4/o/vvvV15enmJiYrRu3Tr5/X55PB5t3bpVa9eu1dSpUyVJ\nGRkZqq2tldfrVXl5uXp7e7VmzZqInxAAAADGB0c3/DfB0aMf8CP4IYqLi1FS0mQyc4jcnCOz8JCb\nc2QWHnJzjszCczK3SHH8E1UAAABgNFBUAQAAYCSKKgAAAIxEUQUAAICRKKoAAAAwEkUVAAAARqKo\nAgAAwEgUVQAAABiJogoAAAAjUVQBAABgJIoqAAAAjERRBQAAgJEoqgAAADASRRUAAABGoqgCAADA\nSBRVAAAAGImiCgAAACNRVAEAAGAkiioAAACMRFEFAACAkSiqAAAAMFJctAdw4k8vvKDe3uMaGDgR\n7VFG3BnTpmnGjKxojwEAABA1Y6qo3nj3/+rUpDOiPcaoOCNuux6449ZojwEAABA1Y6qonpp0hhJT\nz4z2GKPCHeiO9ggAAABRxTWqAAAAMBJFFQAAAEaiqAIAAMBIFFUAAAAYiaIKAAAAI1FUAQAAYCSK\nKgAAAIxEUQUAAICRKKoAAAAwEkUVAAAARqKoAgAAwEiOi+rhw4f1/e9/X0VFRTrvvPP061//OrjW\n0dGhJUuWqLCwUAsXLlRDQ0PIY3fs2KFFixapoKBAFRUVam9vH/4ZAAAAYFxyXFR/8IMfaPLkyXrq\nqad000036b777tP27dslSUuXLlVaWpq8Xq8WL16s6upqHT58WJJ06NAhVVVVyePxyOv1KikpSVVV\nVZE9GwAAAIwbjopqT0+PXnvtNV1//fXKzMzUeeedp9LSUu3atUu7du1SR0eHVqxYoRkzZqiyslIF\nBQWqr6+XJG3ZskWzZ89WRUWFsrKytHLlSh08eFCNjY0jcmIAAAAY2xwV1fj4eCUkJMjr9aq/v19v\nvfWWmpublZOTo9dee015eXlyu93B/YuKirRnzx5JUmtrq4qLi0OOlZubq5aWlgidCgAAAMYTR0XV\n5XLp5ptv1n//938rPz9fF154oc4991x5PB75/X6lpaWF7J+cnKzOzk5JUldX16D1lJSU4DoAAADw\ncXFOH7B//34tWLBA11xzjf72t7/ptttu0znnnKNAICCXyxWyr8vlkmVZkqS+vr7PXEeoSZOkuLjh\n3ZQhNjYm5DOGhtycI7PwkJtzZBYecnOOzMIT6bwcFdWdO3eqvr5eL730klwul3Jzc3X48GGtX79e\n55xzjt5///2Q/S3LUnx8vCTJ7XYPKqWWZSkxMXGYpzA+uVxxSkqaHJFjJSYmROQ4Ew25OUdm4SE3\n58gsPOTmHJlFl6Oi+vrrr+vMM88M+cloTk6ONmzYoPT0dL355psh+3d3dys1NVWSlJ6eLr/fP2g9\nJycn3NnHNcvq19GjHwzrGLGxMUpMTFBPT0ADAyciNNn4R27OkVl4yM05MgsPuTlHZuE5mVukOCqq\naWlpOnDggPr7+xUX938Pfeutt/SFL3xB+fn52rBhgyzLChbZpqYmzZ07V5KUn5+v5ubm4LECgYDa\n2tq0bNmySJ3LuGLbUn9/ZF4YAwMnInasiYTcnCOz8JCbc2QWHnJzjsyiy9GFBAsWLFBcXJz+8z//\nU++8847++Mc/asOGDfrud7+r4uJiTZs2TTU1NfL5fNq4caP27t2rsrIySZLH41Fzc7Pq6urk8/m0\nfPlyZWZmqqSkZERODAAAAGObo6J66qmn6uGHH5bf71d5ebnuuOMOVVVVqby8XDExMVq/fr38fr88\nHo+2bt2qtWvXaurUqZKkjIwM1dbWyuv1qry8XL29vVqzZs2InBQAAADGPse/9Z+VlaVNmzb9y7Xp\n06dr8+bNn/rY0tJSbdu2zelTAgAAYALingsAAAAwEkUVAAAARqKoAgAAwEgUVQAAABiJogoAAAAj\nUVQBAABgJIoqAAAAjERRBQAAgJEoqgAAADASRRUAAABGoqgCAADASBRVAAAAGImiCgAAACNRVAEA\nAGAkiioAAACMRFEFAACAkSiqAAAAMBJFFQAAAEaiqAIAAMBIFFUAAAAYiaIKAAAAI1FUAQAAYCSK\nKgAAAIxEUQUAAICRKKoAAAAwEkUVAAAARqKoAgAAwEgUVQAAABiJogoAAAAjUVQBAABgJIoqAAAA\njERRBQAAgJEoqgAAADASRRUAAABGoqgCAADASBRVAAAAGMlRUX3qqaeUnZ2tnJyckM+5ubmSpPb2\ndi1ZskSFhYVauHChGhoaQh6/Y8cOLVq0SAUFBaqoqFB7e3vkzgQAAADjiqOi+q1vfUsNDQ165ZVX\n1NDQoD/96U/64he/qO9973uSpKqqKqWlpcnr9Wrx4sWqrq7W4cOHJUmHDh1SVVWVPB6PvF6vkpKS\nVFVVFfkzAgAAwLjgqKi6XC4lJycHP55++mlJ0g9/+EPt3LlTHR0dWrFihWbMmKHKykoVFBSovr5e\nkrRlyxbNnj1bFRUVysrK0sqVK3Xw4EE1NjZG/qwAAAAw5oV9jeqxY8f04IMP6sYbb9TnPvc5tba2\nKi8vT263O7hPUVGR9uzZI0lqbW1VcXFxcC0+Pl65ublqaWkZxvgAAAAYr8Iuqo899pjS09P1ta99\nTZLk9/uVlpYWsk9ycrI6OzslSV1dXYPWU1JSgusAAADAx8WF+8D6+npVVlYGvw4EAnK5XCH7uFwu\nWZYlSerr6/vMdYSaNEmKixveTRliY2NCPmNoyM05MgsPuTlHZuEhN+fILDyRziusotra2qrOzk5d\neOGFwW1ut1vHjh0L2c+yLMXHxwfXP1lKLctSYmJiOCOMey5XnJKSJkfkWImJCRE5zkRDbs6RWXjI\nzTkyCw+5OUdm0RVWUX3llVdUXFysKVOmBLelp6fL5/OF7Nfd3a3U1NTgut/vH7Sek5MTzgjjnmX1\n6+jRD4Z1jNjYGCUmJqinJ6CBgRMRmmz8IzfnyCw85OYcmYWH3Jwjs/CczC1Swv6J6pw5c0K25efn\nq66uTpZlBd/ib2pq0ty5c4Przc3Nwf0DgYDa2tq0bNmycGcf12xb6u+PzAtjYOBExI41kZCbc2QW\nHnJzjszCQ27OkVl0hXUhwd/+9jdlZWWFbCspKdG0adNUU1Mjn8+njRs3au/evSorK5MkeTweNTc3\nq66uTj6fT8uXL1dmZqZKSkqGfxYAAAAYd8Iqqu+9954+//nPhx4oJkbr1q2T3++Xx+PR1q1btXbt\nWk2dOlWSlJGRodraWnm9XpWXl6u3t1dr1qwZ/hkAAABgXArrrf+T90b9pOnTp2vz5s2f+rjS0lJt\n27YtnKcEAADABMM9FwAAAGAkiioAAACMRFEFAACAkSiqAAAAMBJFFQAAAEaiqAIAAMBIFFUAAAAY\niaIKAAAAI1FUAQAAYCSKKgAAAIxEUQUAAICRKKoAAAAwEkUVAAAARqKoAgAAwEgUVQAAABiJogoA\nAAAjUVQBAABgJIoqAAAAjERRBQAAgJEoqgAAADASRRUAAABGoqgCAADASBRVAAAAGImiCgAAACNR\nVAEAAGAkiioAAACMRFEFAACAkSiqAAAAMBJFFQAAAEaiqAIAAMBIFFUAAAAYiaIKAAAAI1FUAQAA\nYCSKKgAAAIxEUQUAAICRHBdVy7L085//XCUlJZo/f77uvffe4FpHR4eWLFmiwsJCLVy4UA0NDSGP\n3bFjhxYtWqSCggJVVFSovb19+GcAAACAcclxUf3FL36hnTt36le/+pXuuusubdmyRVu2bJEkLV26\nVGlpafJ6vVq8eLGqq6t1+PBhSdKhQ4dUVVUlj8cjr9erpKQkVVVVRfZsAAAAMG7EOdn52LFjevLJ\nJ/Xwww9r1qxZkqSrr75ar732mjIzM9XR0aEnnnhCbrdblZWV2rlzp+rr61VdXa0tW7Zo9uzZqqio\nkCStXLlS8+bNU2Njo4qLiyN+YgAAABjbHBXVpqYmTZkyRXPnzg1uu+666yRJGzZsUF5entxud3Ct\nqKhIe/bskSS1traGFNL4+Hjl5uaqpaWFogoAAIBBHL31397eroyMDP32t7/VN7/5TZ1//vlat26d\nbNuW3+9XWlpayP7Jycnq7OyUJHV1dQ1aT0lJCa4DAAAAH+foJ6offvih3nnnHW3ZskWrVq2S3+/X\nzTffrISEBAUCAblcrpD9XS6XLMuSJPX19X3mOkJNmiTFxQ3vpgyxsTEhnzE05OYcmYWH3Jwjs/CQ\nm3NkFp5I5+WoqMbGxuqDDz7QPffco6lTp0qSDh48qMcee0zz58/X+++/H7K/ZVmKj4+XJLnd7kGl\n1LIsJSYmDmf+ccvlilNS0uSIHCsxMSEix5loyM05MgsPuTlHZuEhN+fILLocFdW0tDS53e5gSZWk\ns846S52dnUpPT9ebb74Zsn93d7dSU1MlSenp6fL7/YPWc3Jywp19XLOsfh09+sGwjhEbG6PExAT1\n9AQ0MHAiQpONf+TmHJmFh9ycI7PwkJtzZBaek7lFiqOimp+fr+PHj+vAgQP64he/KEnav3+/MjIy\nlJ+frw0bNsiyrOBb/E1NTcFfvMrPz1dzc3PwWIFAQG1tbVq2bFmkzmVcsW2pvz8yL4yBgRMRO9ZE\nQm7OkVl4yM05MgsPuTlHZtHl6EKCs846S1/96ldVU1Ojffv26eWXX1ZdXZ0uv/xyFRcXa9q0aaqp\nqZHP59PGjRu1d+9elZWVSZI8Ho+am5tVV1cnn8+n5cuXKzMzUyUlJSNyYgAAABjbHF/xetddd+mL\nX/yirrjiCi1fvlxXXXWVrrjiCsXExGj9+vXy+/3yeDzaunWr1q5dG7xMICMjQ7W1tfJ6vSovL1dv\nb6/WrFkT8RMCAADA+ODorX9JOvXUU7Vq1SqtWrVq0Nr06dO1efPmT31saWmptm3b5vQpAQAAMAFx\nzwUAAAAYiaIKAAAAI1FUAQAAYCSKKgAAAIxEUQUAAICRKKoAAAAwEkUVAAAARqKoAgAAwEgUVQAA\nABiJogoAAAAjUVQBAABgJIoqAAAAjERRBQAAgJEoqgAAADASRRUAAABGoqgCAADASBRVAAAAGImi\nCgAAACNRVAEAAGAkiioAAACMRFEFAACAkSiqAAAAMBJFFQAAAEaiqAIAAMBIFFUAAAAYiaIKAAAA\nI1FUAQAAYCSKKgAAAIxEUQUAAICRKKoAAAAwEkUVAAAARqKoAgAAwEgUVQAAABiJogoAAAAjUVQB\nAABgJMdFdfv27crOzlZOTk7w8w9+8ANJUkdHh5YsWaLCwkItXLhQDQ0NIY/dsWOHFi1apIKCAlVU\nVKi9vT0yZwEAAIBxx3FR9fl8WrBggRoaGtTQ0KBXXnlFt99+uyRp6dKlSktLk9fr1eLFi1VdXa3D\nhw9Lkg4dOqSqqip5PB55vV4lJSWpqqoqsmcDAACAccNxUd2/f7++/OUv6/TTT1dycrKSk5N16qmn\naufOnero6NCKFSs0Y8YMVVZWqqCgQPX19ZKkLVu2aPbs2aqoqFBWVpZWrlypgwcPqrGxMeInBQAA\ngLEvrKJ61llnDdre2tqqvLw8ud3u4LaioiLt2bMnuF5cXBxci4+PV25urlpaWsKZGwAAAOOc46L6\n9ttv6+WXX9bXv/51fe1rX9Pdd9+tjz76SH6/X2lpaSH7Jicnq7OzU5LU1dU1aD0lJSW4DgAAAHxc\nnJOd//73v6uvr09ut1v333+/Ojo6dPvtt6uvr0+BQEAulytkf5fLJcuyJEl9fX2fuY5QkyZJcXHD\nuylDbGxMyGcMDbk5R2bhITfnyCw85OYcmYUn0nk5KqpnnHGGXn31VSUmJkqSsrOzdeLECf34xz/W\nt7/9bfX09ITsb1mW4uPjJUlut3tQKbUsK3gshHK54pSUNDkix0pMTIjIcSYacnOOzMJDbs6RWXjI\nzTkyiy5HRVXSoGKZlZWl48ePKyUlRfv37w9Z6+7uVmpqqiQpPT1dfr9/0HpOTo7TESYEy+rX0aMf\nDOsYsbExSkxMUE9PQAMDJyI02fhHbs6RWXjIzTkyCw+5OUdm4TmZW6Q4KqqvvPKKfvSjH+mll14K\n/tJUW1ubkpKSNHfuXP3qV7+SZVnBt/ibmpo0d+5cSVJ+fr6am5uDxwoEAmpra9OyZcsidS7jim1L\n/f2ReWEMDJyI2LEmEnJzjszCQ27OkVl4yM05MosuRxcSFBYWKiEhQT/96U/19ttv68UXX9Tq1at1\n3XXXqbi4WNOmTVNNTY18Pp82btyovXv3qqysTJLk8XjU3Nysuro6+Xw+LV++XJmZmSopKRmREwMA\nAMDY5qioTp48WZs2bdLRo0dVVlamn/3sZ7r00kt19dVXKyYmRuvXr5ff75fH49HWrVu1du1aTZ06\nVZKUkZGh2tpaeb1elZeXq7e3V2vWrBmRkwIAAMDY5/ga1aysLG3atOlfrk2fPl2bN2/+1MeWlpZq\n27ZtTp8SAAAAExD3XAAAAICRKKoAAAAwEkUVAAAARqKoAgAAwEgUVQAAABiJogoAAAAjUVQBAABg\nJIoqAAAAjERRBQAAgJEoqgAAADASRRUAAABGoqgCAADASBRVAAAAGImiCgAAACNRVAEAAGAkiioA\nAACMRFEFAACAkSiqAAAAMBJFFQAAAEaiqAIAAMBIFFUAAAAYiaIKAAAAI1FUAQAAYCSKKgAAAIxE\nUQUAAICRKKoAAAAwEkUVAAAARqKoAgAAwEgUVQAAABiJogoAAAAjUVQBAABgJIoqAAAAjERRBQAA\ngJEoqgAAADASRRUAAABGoqgCAADASGEX1crKSi1fvjz4dUdHh5YsWaLCwkItXLhQDQ0NIfvv2LFD\nixYtUkFBgSoqKtTe3h7+1AAAABj3wiqqzzzzjF566aWQbVVVVUpLS5PX69XixYtVXV2tw4cPS5IO\nHTqkqqoqeTweeb1eJSUlqaqqavjTAwAAYNxyXFSPHTum1atX6+yzzw5u27lzp9rb27VixQrNmDFD\nlZWVKigoUH19vSRpy5Ytmj17tioqKpSVlaWVK1fq4MGDamxsjNyZAAAAYFxxXFTvuOMOXXTRRcrK\nygpua21tVV5entxud3BbUVGR9uzZE1wvLi4OrsXHxys3N1ctLS3DmR0AAADjmKOiunPnTjU1NQ16\n297v9ystLS1kW3Jysjo7OyVJXV1dg9ZTUlKC6wAAAMAnxQ11R8uydOutt+qWW26Ry+UKWQsEAoO2\nuVwuWZYlSerr6/vMdQw2aZIUFze8mzLExsaEfMbQkJtzZBYecnOOzMJDbs6RWXgindeQi2ptba1m\nzZqlr3zlK4PW3G63jh07FrLNsizFx8cH1z9ZSi3LUmJiYjgzTwguV5ySkiZH5FiJiQkROc5EQ27O\nkVl4yM05MgsPuTlHZtE15KL6P//zPzpy5IgKCwslSR999JEk6dlnn9X3v/99+Xy+kP27u7uVmpoq\nSUpPT5ff7x+0npOTM6zhxzPL6tfRox8M6xixsTFKTExQT09AAwMnIjTZ+EduzpFZeMjNOTILD7k5\nR2bhOZlbpAy5qD766KPq7+8Pfr169WpJ0o9//GMdPHhQGzdulGVZwbf4m5qaNHfuXElSfn6+mpub\ng48NBAJqa2vTsmXLInIS45FtS/39kXlhDAyciNixJhJyc47MwkNuzpFZeMjNOTKLriFfSDBt2jRN\nnz49+DF58mRNnjxZ06dPV0lJiaZNm6aamhr5fD5t3LhRe/fuVVlZmSTJ4/GoublZdXV18vl8Wr58\nuTIzM1VSUjJiJwYAAICxLSJXvMbExGjdunXy+/3yeDzaunWr1q5dq6lTp0qSMjIyVFtbK6/Xq/Ly\ncvX29mrNmjWReGoAAACMU0N+6/+TVq5cGfL19OnTtXnz5k/dv7S0VNu2bQv36QAAADDBcM8FAAAA\nGImiCgAAACNRVAEAAGAkiioAAACMRFEFAACAkSiqAAAAMBJFFQAAAEaiqAIAAMBIFFUAAAAYiaIK\nAAAAI1FUAQAAYCSKKgAAAIxEUQUAAICRKKoAAAAwEkUVAAAARqKoAgAAwEgUVQAAABiJogoAAAAj\nUVQBAABgJIoqAAAAjERRBQAAgJEoqgAAADASRRUAAABGoqgCAADASBRVAAAAGImiCgAAACNRVAEA\nAGAkiioAAACMRFEFAACAkSiqAAAAMBJFFQAAAEaiqAIAAMBIFFUAAAAYiaIKAAAAI1FUAQAAYCTH\nRfXdd9/VNddco8LCQi1YsECbNm0KrnV0dGjJkiUqLCzUwoUL1dDQEPLYHTt2aNGiRSooKFBFRYXa\n29uHfwYAAAAYlxwVVdu2VVlZqZSUFD399NO69dZbtX79ej3zzDOSpKVLlyotLU1er1eLFy9WdXW1\nDh8+LEk6dOiQqqqq5PF45PV6lZSUpKqqqsifEQAAAMYFR0W1u7tbubm5uuWWW5SZmalzzz1X55xz\njpqamrRr1y51dHRoxYoVmjFjhiorK1VQUKD6+npJ0pYtWzR79mxVVFQoKytLK1eu1MGDB9XY2Dgi\nJwYAAICxzVFRTU1N1T333KNTTjlFktTU1KTdu3erpKREr732mvLy8uR2u4P7FxUVac+ePZKk1tZW\nFRcXB9fi4+OVm5urlpaWSJwHAAAAxpmwf5lqwYIFuvLKK1VQUKALLrhAfr9faWlpIfskJyers7NT\nktTV1TVoPSUlJbgOAAAAfFxcuA+sra1Vd3e3br31Vv3Xf/2XAoGAXC5XyD4ul0uWZUmS+vr6PnMd\noSZNkuLihndThtjYmJDPGBpyc47MwkNuzpFZeMjNOTILT6TzCruo5uXlSZJqamp04403qqysTD09\nPSH7WJal+Ph4SZLb7R5USi3LUmJiYrgjjGsuV5ySkiZH5FiJiQkROc5EQ27OkVl4yM05MgsPuTlH\nZtHlqKgeOXJELS0tOv/884PbvvSlL+mjjz5Samqq9u/fH7J/d3e3UlNTJUnp6eny+/2D1nNycsKd\nfVyzrH4dPfrBsI4RGxujxMQE9fQENDBwIkKTjX/k5hyZhYfcnCOz8JCbc2QWnpO5RYqjotrR0aFl\ny5bpxRdfDF5vunfvXiUnJ6uoqEibNm2SZVnBt/ibmpo0d+5cSVJ+fr6am5uDxwoEAmpra9OyZcsi\ndS7jim1L/f2ReWEMDJyI2LEmEnJzjszCQ27OkVl4yM05MosuRxcSzJ49W7NmzdJNN92k/fv368UX\nX9Rdd92l66+/XsXFxZo2bZpqamrk8/m0ceNG7d27V2VlZZIkj8ej5uZm1dXVyefzafny5crMzFRJ\nScmInBgAAADGNkdFNSYmRuvWrdMpp5yiSy+9VD/72c/03e9+V1deeaViYmK0fv16+f1+eTwebd26\nVWvXrtXUqVMlSRkZGaqtrZXX61V5ebl6e3u1Zs2aETkpAAAAjH2Of5kqNTVVDzzwwL9cmz59ujZv\n3vypjy0tLdW2bducPiUAAAAmIO65AAAAACNRVAEAAGAkiioAAACMRFEFAACAkSiqAAAAMBJFFQAA\nAEaiqAIAAMBIFFUAAAAYiaIKAAAAI1FUAQAAYCSKKgAAAIxEUQUAAICRKKoAAAAwEkUVAAAARqKo\nAgAAwEgUVQAAABiJogoAAAAjUVQBAABgJIoqAAAAjERRBQAAgJEoqgAAADASRRUAAABGoqgCAADA\nSBRVAAAAGImiCgAAACNRVAEAAGAkiioAAACMRFEFAACAkSiqAAAAMBJFFQAAAEaiqAIAAMBIFFUA\nAAAYiaI3gA3ZAAAVr0lEQVQKAAAAI1FUAQAAYCSKKgAAAIxEUQUAAICRHBXVzs5O3XDDDfq3f/s3\nffWrX9WqVatkWZYkqaOjQ0uWLFFhYaEWLlyohoaGkMfu2LFDixYtUkFBgSoqKtTe3h65swAAAMC4\n46io3nDDDTp+/Lgee+wx3XPPPfrTn/6k+++/X5K0dOlSpaWlyev1avHixaqurtbhw4clSYcOHVJV\nVZU8Ho+8Xq+SkpJUVVUV+bMBAADAuDHkovrWW2+ptbVVK1euVFZWloqKinTDDTfo97//vXbt2qWO\njg6tWLFCM2bMUGVlpQoKClRfXy9J2rJli2bPnq2KigplZWVp5cqVOnjwoBobG0fsxAAAADC2Dbmo\npqam6sEHH9Tpp58esr23t1evvfaa8vLy5Ha7g9uLioq0Z88eSVJra6uKi4uDa/Hx8crNzVVLS8tw\n5wcAAMA4NeSiOmXKFM2bNy/4tW3bevTRR3XOOefI7/crLS0tZP/k5GR1dnZKkrq6ugatp6SkBNcB\nAACAT4oL94F33nmn3njjDdXX1+uhhx6Sy+UKWXe5XMFftOrr6/vMdQw2aZIUFze8mzLExsaEfMbQ\nkJtzZBYecnOOzMJDbs6RWXginVdYRXX16tXavHmz7rvvPn3pS1+S2+3WsWPHQvaxLEvx8fGSJLfb\nPaiUWpalxMTEMMce/1yuOCUlTY7IsRITEyJynImG3Jwjs/CQm3NkFh5yc47MostxUb3tttv0+OOP\na/Xq1Tr//PMlSenp6fL5fCH7dXd3KzU1Nbju9/sHrefk5IQ797hnWf06evSDYR0jNjZGiYkJ6ukJ\naGDgRIQmG//IzTkyCw+5OUdm4SE358gsPCdzixRHRXXNmjV6/PHHde+99+prX/tacHt+fr7q6upk\nWVbwLf6mpibNnTs3uN7c3BzcPxAIqK2tTcuWLYvEOYxLti3190fmhTEwcCJix5pIyM05MgsPuTlH\nZuEhN+fILLqGfCHB/v37tX79elVWVqqwsFDd3d3Bj5KSEk2bNk01NTXy+XzauHGj9u7dq7KyMkmS\nx+NRc3Oz6urq5PP5tHz5cmVmZqqkpGTETgwAAABj25CL6h/+8AedOHFC69evV2lpqUpLSzV//nyV\nlpYqJiZGa9euld/vl8fj0datW7V27VpNnTpVkpSRkaHa2lp5vV6Vl5ert7dXa9asGbGTAgAAwNg3\n5Lf+KysrVVlZ+anrmZmZ2rx586eul5aWatu2bc6mAwAAwITFPRcAAABgJIoqAAAAjERRBQAAgJEo\nqgAAADASRRUAAABGCuu/UMXIOjHwkbq7/q6WlqZhHWcs/a8aeXmzg/9ZBAAAgERRNVKP/x31/mOy\nbvv17miPMip6j7yrO38oFRYWRXsUAABgEIqqoaYkZ+q0qV+O9hgAAABRwzWqAAAAMBJFFQAAAEai\nqAIAAMBIFFUAAAAYiaIKAAAAI1FUAQAAYCSKKgAAAIxEUQUAAICRKKoAAAAwEkUVAAAARqKoAgAA\nwEgUVQAAABiJogoAAAAjUVQBAABgJIoqAAAAjERRBQAAgJEoqgAAADASRRUAAABGoqgCAADASBRV\nAAAAGImiCgAAACNRVAEAAGAkiioAAACMRFEFAACAkSiqAAAAMBJFFQAAAEaiqAIAAMBIYRdVy7K0\naNEiNTY2Brd1dHRoyZIlKiws1MKFC9XQ0BDymB07dmjRokUqKChQRUWF2tvbw58cAAAA41pYRdWy\nLP3whz+Uz+cL2V5VVaW0tDR5vV4tXrxY1dXVOnz4sCTp0KFDqqqqksfjkdfrVVJSkqqqqoZ/BgAA\nABiXHBfV/fv36zvf+Y46OjpCtu/cuVPt7e1asWKFZsyYocrKShUUFKi+vl6StGXLFs2ePVsVFRXK\nysrSypUrdfDgwZCfyAIAAAAnOS6qf/7zn3XOOefo8ccfl23bwe2tra3Ky8uT2+0ObisqKtKePXuC\n68XFxcG1+Ph45ebmqqWlZTjzAwAAYJyKc/qAyy677F9u9/v9SktLC9mWnJyszs5OSVJXV9eg9ZSU\nlOA6AAAA8HGOi+qnCQQCcrlcIdtcLpcsy5Ik9fX1feY6JrbY2BjFxZlxE4rY2JiQz/h/I7PwkJtz\nZBYecnOOzMIT6bwiVlTdbreOHTsWss2yLMXHxwfXP1lKLctSYmJipEbAGJaYmKCkpMnRHiNEYmJC\ntEcYc8gsPOTmHJmFh9ycI7PoilhRTU9PH3QXgO7ubqWmpgbX/X7/oPWcnJxIjYAxrKcnoKNHP4j2\nGJL+71+DiYkJ6ukJaGDgRLTHGRPILDzk5hyZhYfcnCOz8JzMLVIiVlTz8/NVV1cny7KCb/E3NTVp\n7ty5wfXm5ubg/oFAQG1tbVq2bFmkRsAYNjBwQv39Zv1BYOJMpiOz8JCbc2QWHnJzjsyiK2IXEpSU\nlGjatGmqqamRz+fTxo0btXfvXpWVlUmSPB6PmpubVVdXJ5/Pp+XLlyszM1MlJSWRGgEAAADjyLCK\n6qRJk/55oJgYrVu3Tn6/Xx6PR1u3btXatWs1depUSVJGRoZqa2vl9XpVXl6u3t5erVmzZnjTAwAA\nYNwa1lv/b7zxRsjX06dP1+bNmz91/9LSUm3btm04TwkAAIAJgnsuAAAAwEgUVQAAABiJogoAAAAj\nUVQBAABgJIoqAAAAjERRBQAAgJEoqgAAADASRRUAAABGGtYN/4FIODHQr7/+dV+0xwiKjY1RYmKC\nenoCGhiI/P/vnJc3Wy6XK+LHBQBgvKGoIuo+eP+QNj3zd03Z9Y9ojzLieo+8qzt/KBUWFkV7FAAA\njEdRhRGmJGfqtKlfjvYYAADAIFyjCgAAACNRVAEAAGAkiioAAACMRFEFAACAkSiqAAAAMBJFFQAA\nAEaiqAIAAMBIFFUAAAAYiaIKAAAAI1FUAQAAYCSKKgAAAIxEUQUAAICRKKoAAAAwEkUVAAAARqKo\nAgAAwEhx0R4AmEhODPTrr3/dF+0xIi42NkaJiQnq6QloYOBEyFpe3my5XK4oTQYAGMsoqsAo+uD9\nQ9r0zN81Zdc/oj3KqOg98q7u/KFUWFgU7VEAAGMQRRUYZVOSM3Xa1C9HewwAAIzHNaoAAAAwEkUV\nAAAARqKoAgAAwEgUVQAAABiJX6YCMGLG6+24Pg234gKAyKKoAhgxE+l2XMf8b+u6Rfs0c2b2kB/z\nWfefHQso5gBG2qgWVcuydOutt+r5559XfHy8rr76ai1ZsmQ0RwAwyibK7bh6j7Rr0zNtE6KUS9wj\nF8DoGNWiescdd6itrU2bN29WR0eHfvKTnygjI0MXXHDBaI4BACNiopRyABgto1ZUA4GA6uvrtWnT\nJmVnZys7O1vXXnutHn30UYoqAIwx0br+OFqXS3CZAxAdo1ZU9+3bp4GBARUUFAS3FRUVacOGDaM1\nAgAgQrj+eGxxUvA/+ugjSdLnPve50Rgt6vhHiNlGraj6/X6ddtppiov751MmJyfr+PHjOnr0qJKS\nkkZrFABABEyUSx0m2vXHnW816pTPp2tKcma0Rxlxn/WPkLH+y46fZqwV81F96/+TwZz82rKs0Rpj\nzOg98m60Rxg1Hx47LMmO9hijYiKdqzSxzncinas0sc73w2OHdcrn06M9BkZAX2+37nvkWZ2S2BLt\nUUbFhz1d2rjq/9OcOSP3S5CxsZG9Rf+oFVW32z2okJ78OiEhYUjHePGRH0R8LgAAAJhp1P5nqvT0\ndL3//vs6ceKfPz7v7u5WfHy8EhMTR2sMAAAAjBGjVlRzcnIUFxenPXv2BLft3r1bs2bNGq0RAAAA\nMIaMWlGNj4/XRRddpFtuuUV79+7V9u3b9dBDD+l73/veaI0AAACAMWSSbdujdjV8X1+ffv7zn+vZ\nZ5/VlClTdO211+qqq64aracHAADAGDKqRRUAAAAYqlF76x8AAABwgqIKAAAAI1FUAQAAYCSKKgAA\nAIxEUQUAAICRjCqqlmXppptuUnFxsUpLS/XQQw996r5tbW36zne+o4KCApWXl+v1118fxUnN4SSz\nk3bv3q3zzz9/FKYzl5PcXnjhBV188cUqLCzURRddpD/+8Y+jOKk5nGT2u9/9Tl//+teVn5+vyy67\nTK2traM4qVnCeY12dHSosLBQjY2NozCheZxkdv311ys7O1s5OTnBzy+++OIoTmsOJ7n99a9/1eWX\nX678/HwtXrxYr7766ihOao6hZnbVVVcpOzt70MdPf/rTUZ7YDE6+155//nldeOGFKiws1BVXXKG2\ntjZnT2YbZMWKFfZFF11kv/HGG/bzzz9vz5kzx3722WcH7ffhhx/a8+bNs++88057//799i9+8Qt7\n3rx5diAQiMLU0TXUzE7at2+fPW/ePHvBggWjOKV5hprbG2+8Yc+aNct+9NFH7Xfffdd+9NFH7by8\nPHvfvn1RmDq6hppZY2OjPXv2bHvr1q12e3u7vWrVKrukpMT+8MMPozB19Dl9jdq2bV9zzTV2dna2\n/ec//3mUpjSLk8wuuOAC+/e//73d3d0d/LAsa5QnNsNQc+vt7bXnzZtn33zzzfa7775rP/DAA/bc\nuXPtI0eORGHq6BpqZseOHQv5Htu+fbs9e/Zs+/XXX4/C1NE31NzefPNN++yzz7affvpp+91337VX\nrFhhz5s3z+7r6xvycxlTVD/88EP77LPPthsbG4Pb1q1bZ1911VWD9n3iiSfs888/P2TbBRdcYD/1\n1FMjPqdJnGRm27b9m9/8xi4sLLQvuuiiCV1UneR211132dddd13Itquvvtq+9957R3xOkzjJ7H//\n93/tX/7yl8Gve3t77ZkzZ9qtra2jMqtJnL5Gbdu2n376afuyyy6bsEXVSWbHjx+3c3Nz7XfeeWc0\nRzSSk9x+/etf2xdccEHItrKyMvvFF18c8TlNEs7r07Zte2BgwP7Wt75lP/DAAyM9opGc5PbQQw/Z\nHo8n+PU//vEPe+bMmfZf/vKXIT+fMW/979u3TwMDAyooKAhuKyoq+pdvGba2tqqoqChk25w5c9TS\n0jLic5rESWaS9Morr+jOO++c8P9trZPcLrnkEv3oRz8atP0f//jHiM5oGieZfeMb39C///u/S5KO\nHz+uhx9+WCkpKfrSl740avOawulr9OjRo7r77rt12223yZ6g/xeLk8zefvttTZo0SdOnTx/NEY3k\nJLfGxkYtWLAgZNsTTzyhc889d8TnNInT1+dJXq9Xx44d07XXXjvSIxrJSW6nnXaafD6fmpubZdu2\nvF6vpkyZoszMzCE/nzFF1e/367TTTlNcXFxwW3Jyso4fP66jR4+G7NvV1aW0tLSQbcnJyers7ByV\nWU3hJDNJWrNmzYS/NlVyltuMGTM0c+bM4Ndvvvmmdu3apXPOOWfU5jWB0+81Sdq5c6cKCwu1bt06\n3XTTTUpISBitcY3hNLdVq1bpkksuUVZW1miOaRQnme3fv1+nnnqqfvzjH2v+/PkqLy/XSy+9NNoj\nG8FJbu3t7UpKStLNN9+s+fPn69JLL1Vzc/Nojxx14fy5JkkPPvigKioqJuSfaZKz3C688EKde+65\nuvzyyzVr1iytXr1aDzzwgKZMmTLk5zOmqAYCAblcrpBtJ7+2LCtke19f37/c95P7jXdOMsM/hZvb\ne++9p2XLlqmoqEjnnXfeiM5omnAymzlzpp588kndcMMN+slPfjIhf6HKSW47duxQS0uLli5dOmrz\nmchJZm+99ZaOHz+u0tJSbdq0SV/96ld1/fXXT8hfrnWS24cffqgHH3xQaWlpevDBBzV37lxdc801\nE+6HPeH8ubZr1y51dXWpvLx8xOczlZPc3n//fXV3d+uWW27RE088oYsvvlg1NTV67733hvx8xhRV\nt9s96ARPfv3Jf7V82r7x8fEjO6RhnGSGfwont+7ubn3ve9/TpEmTdP/994/4jKYJJ7PTTz9d2dnZ\n+v73v6958+bpN7/5zYjPaZqh5nb8+HHdcsstuuWWWwb9BTDROPleq66u1ksvvaSLL75YM2fOVHV1\ntUpLS/X444+P2rymcJJbbGyscnJyVF1drezsbN14440688wz9fTTT4/avCYI58+15557TqWlpUpM\nTBzx+UzlJLe77rpLM2fO1GWXXabc3FytWLFCCQkJevLJJ4f8fMYU1fT0dL3//vs6ceJEcFt3d7fi\n4+MHfUOkp6fL7/eHbOvu7lZqauqozGoKJ5nhn5zm1tnZqSuuuEIDAwPavHmzkpKSRnNcIzjJbO/e\nvYNuP5KVlfWZb6WNV0PNrbW1VR0dHVq2bJkKCwtVWFgoSbruuut06623jvbYUeX09fnJtxCzsrLU\n1dU14nOaxkluqampmjFjRsi2M888U4cOHRqVWU0Rzt+hL7/88oS/hM5Jbq+//rqys7ODX0+aNEnZ\n2dn6+9//PuTnM6ao5uTkKC4uTnv27Alu2717t2bNmjVo3/z8/EG/ONXc3BxyYe9E4CQz/JOT3AKB\ngK699lp97nOf06OPPqqUlJTRHNUYTjKrr6/X3XffHbLt9ddfn5DXXQ41t/z8fD333HN6+umn9bvf\n/U6/+93vJEm33367brjhhlGdOdqcfK8tX75cN910U8i2ffv26ayzzhrxOU3jJLeCggLt27cvZNtb\nb72ljIyMEZ/TJE7/Dj169Kja29s1Z86c0RrRSE5yS0tLk8/nC9n29ttv6wtf+MLQn3A4tyiItJtv\nvtleuHCh3draaj///PN2UVGR/fzzz9u2bdt+vz94363e3l77K1/5in377bfbPp/Pvu222+z58+dP\nyPuoDjWzj3vyyScn9O2pbHvoud1zzz12QUGB3draavv9/uBHb29vNMePiqFm9vrrr9t5eXn2I488\nYr/zzjv2/fffb8+ZM8fu7OyM5vhRE85r1LZte+bMmRPy9lS2PfTMnnvuOXvWrFn2U089ZR84cMCu\nra21CwoK7IMHD0Zz/KgZam4HDx60CwsL7draWvvAgQP2fffdN2Ffo05en6+++qqdn58frVGNMtTc\nnnnmGTs/P9/+7W9/ax84cMBevXq1XVxc7OievUYV1UAgYNfU1NiFhYX2ueeeaz/yyCPBtZkzZ4bc\nJ7W1tdW+5JJL7Pz8fPs73/mO/cYbb0Rj5KhzktlJFNWh5/aNb3zDzs7OHvRRU1MTrdGjxsn32gsv\nvGAvWrTIzs/Pt8vKyuw9e/ZEY2QjhPMatW17wt5H1badZfbEE0/YF1xwgX322Wfb3/72t+3du3dH\nY2QjOMmtubnZvuSSS+yzzz7bvuSSSyZsbk4ye+aZZ+z58+dHY0zjOMmtvr7e/uY3v2nPmTPHvuKK\nKxz3tUm2PUFv1gcAAACjGXONKgAAAPBxFFUAAAAYiaIKAAAAI1FUAQAAYCSKKgAAAIxEUQUAAICR\nKKoAAAAwEkUVAAAARqKoAgAAwEgUVQAAABiJogoAAAAj/f++Ipfn7GKf4QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10db41090>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">10% 0.207937877481\n",
      ">25% 0.0586712683348\n",
      ">50% 0.0120793787748\n",
      ">75% 0.000862812769629\n",
      ">90% 0.0\n",
      ">95% 0.0\n"
     ]
    }
   ],
   "source": [
    "plt.hist(table.pa(inplace=False).sum(axis=\"observation\")/len(table.ids()))\n",
    "plt.show()\n",
    "print \">10%\", np.sum(table.pa(inplace=False).sum(axis=\"observation\")/len(table.ids())>.10)/table.shape[0]\n",
    "print \">25%\", np.sum(table.pa(inplace=False).sum(axis=\"observation\")/len(table.ids())>.25)/table.shape[0]\n",
    "print \">50%\", np.sum(table.pa(inplace=False).sum(axis=\"observation\")/len(table.ids())>.50)/table.shape[0]\n",
    "print \">75%\", np.sum(table.pa(inplace=False).sum(axis=\"observation\")/len(table.ids())>.75)/table.shape[0]\n",
    "print \">90%\", np.sum(table.pa(inplace=False).sum(axis=\"observation\")/len(table.ids())>.90)/table.shape[0]\n",
    "print \">95%\", np.sum(table.pa(inplace=False).sum(axis=\"observation\")/len(table.ids())>.99)/table.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(631, 168)\n",
      "(516, 171)\n",
      "(1050, 164)\n"
     ]
    }
   ],
   "source": [
    "table_forehead = table.filter([i for i in table.ids() if 'Stool' in i], inplace=False)\n",
    "table_forehead.filter(table_forehead.ids(axis=\"observation\")[table_forehead.sum(axis=\"observation\") != 0], axis=\"observation\")\n",
    "print table_forehead.shape\n",
    "table_tongue = table.filter([i for i in table.ids() if 'Tongue' in i], inplace=False)\n",
    "table_tongue.filter(table_tongue.ids(axis=\"observation\")[table_tongue.sum(axis=\"observation\") != 0], axis=\"observation\")\n",
    "print table_tongue.shape\n",
    "table_stool = table.filter([i for i in table.ids() if 'Forehead' in i], inplace=False)\n",
    "table_stool.filter(table_stool.ids(axis=\"observation\")[table_stool.sum(axis=\"observation\") != 0], axis=\"observation\")\n",
    "print table_stool.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqoAAAHcCAYAAAAeFogrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3X9UVXW+//EXPwL8xVdEUK9XVomT/BABUVpm1I1rzkwj\nqaGtmnIGy5hbiK3V1EqdZXatxkazUvxxr+Z4b1otUyq17tiPuTcrNSXBsNTugJpghqCmmOBR2N8/\nvJw6bXPYeM7hw/H5WOss1vnsfdjv/Vpgr87Z5xBkWZYlAAAAwDDB7T0AAAAAcDEUVQAAABiJogoA\nAAAjUVQBAABgJIoqAAAAjERRBQAAgJEoqgAAADASRRUAAABGoqgCAADASG0uqvn5+Zo+fbr7fnV1\ntSZNmqT09HSNHj1aW7Zs8dh/69atysnJUVpamvLy8lRVVdX2qQEAABDw2lRU3377bX344YceawUF\nBYqNjVVxcbFuu+02TZkyRd98840k6ciRIyooKFBubq6Ki4sVFRWlgoKCy58eAAAAActxUT158qTm\nzZunwYMHu9e2bdumqqoqzZ49W/3791d+fr7S0tK0bt06SdJrr72mlJQU5eXlKT4+XnPmzNHhw4dV\nUlLivTMBAABAQHFcVP/0pz9pzJgxio+Pd6+Vl5crOTlZ4eHh7rWMjAzt2rXLvX3YsGHubREREUpK\nSlJZWdnlzA4AAIAA5qiobtu2TTt37rS9bF9bW6vY2FiPtejoaNXU1EiSjh49atves2dP93YAAADg\nx1pdVF0ul5544gnNmjVLYWFhHtsaGhpsa2FhYXK5XJKkxsbGS24HAAAAfqzVRbWoqEiDBg3S9ddf\nb9sWHh5uK50ul0sRERGt2t5almU52h8AAAAdV2hrd/yv//ovHTt2TOnp6ZKkc+fOSZLeeecd/cu/\n/IsqKio89q+rq1NMTIwkqVevXqqtrbVtT0xMdDRsUFCQTp1qUFNTs6PHBaqQkGBFRnYikx8gEzsy\nsSMTT+RhRyZ2ZGJHJnYtmXhLq4vq6tWrdf78eff9efPmSZIeffRRHT58WMuWLZPL5XK/xL9z504N\nHTpUkpSamqrS0lL3YxsaGrRnzx4VFhY6HripqVnnz/PD8ENkYkcmdmRiRyaeyMOOTOzIxI5MfKfV\nL/336dNH/fr1c9+6dOmiLl26qF+/fsrMzFSfPn00bdo0VVRUaNmyZdq9e7fGjx8vScrNzVVpaamW\nL1+uiooKTZ8+XXFxccrMzPTZiQEAAKBj88qfUA0ODtaSJUtUW1ur3Nxcbdy4UYsXL1bv3r0lSX37\n9lVRUZGKi4s1YcIE1dfXa9GiRd44NAAAAAJUkNXB3qF04sR3PL3+f0JDgxUV1YVMfoBM7MjEjkw8\nkYcdmdiRiR2Z2LVk4i1eeUYVAAAA8DaKKgAAAIxEUQUAAICRKKoAAAAwEkUVAAAARqKoAgAAwEgU\nVQAAABiJogoAAAAjUVQBAABgJIoqAAAAjERRBQAAgJEoqgAAADASRRUAAABGoqgCAADASBRVAAAA\nGImiCgAAACNRVAEAAGAkiioAAACMRFEFAACAkSiqAAAAMBJFFQAAAEYKbe8B2ur06dNasGSZFBzS\n3qP4xOCkgcq59RftPQYAAEC76bBFtbb2qMqqghTZN6W9R/GJkzs+pagCAIArGi/9AwAAwEgUVQAA\nABiJogoAAAAjUVQBAABgJIoqAAAAjERRBQAAgJEoqgAAADASRRUAAABGoqgCAADASBRVAAAAGImi\nCgAAACNRVAEAAGAkiioAAACMRFEFAACAkSiqAAAAMBJFFQAAAEaiqAIAAMBIFFUAAAAYyXFRPXTo\nkO677z6lp6crOztbK1ascG976qmnlJCQoMTERPfXl19+2b1969atysnJUVpamvLy8lRVVeWdswAA\nAEDACXWys2VZys/PV2pqqtavX6+DBw/q4YcfVu/evfWrX/1K+/fv1yOPPKJx48a5H9O1a1dJ0pEj\nR1RQUKCHHnpIWVlZWrRokQoKCrRhwwbvnhEAAAACgqNnVOvq6pSUlKRZs2YpLi5ON954o4YPH66d\nO3dKkiorK5WUlKTo6Gj3LTw8XJK0du1apaSkKC8vT/Hx8ZozZ44OHz6skpIS758VAAAAOjxHRTUm\nJkbPPfecOnfuLEnauXOnSkpKdN111+n06dOqqanR1VdffdHHfvbZZxo2bJj7fkREhJKSklRWVtb2\n6QEAABCwHL30/0PZ2dk6cuSI/umf/kmjRo1SeXm5goKCtHTpUn344Yfq3r27Jk2apLFjx0qSjh49\nqtjYWI/v0bNnT9XU1FzeGQAAACAgtbmoFhUVqa6uTrNmzdLTTz+tQYMGKTg4WPHx8Zo4caJ27Nih\nmTNnqmvXrho5cqQaGxsVFhbm8T3CwsLkcrkcHTck5MKTwKGhwVJQW6c3X3BQ0IVzvISWLFq+gkwu\nhkzsyMQTediRiR2Z2JGJnbezaHNRTU5OliRNnz5djz76qB577DFlZ2crMjJSknTttdfq4MGDevXV\nVzVy5EiFh4fbSqnL5XLv31qRkZ0kScePd1ZQABfVq8JCFBXVpVX7tmSC75GJHZnYkYkn8rAjEzsy\nsSMT33FUVI8dO6aysjKNHDnSvTZgwACdO3dO3333nbp37+6xf//+/bV9+3ZJUq9evVRbW+uxva6u\nTomJiY4GPnWqQU1NzTp58owsy9FDO5RzriadOPHdJfcJCQlWZGQndyYgk4shEzsy8UQedmRiRyZ2\nZGLXkom3OCqq1dXVKiws1ObNm93Xm+7evVs9evTQSy+9pLKyMq1cudK9/969e3XNNddIklJTU1Va\nWure1tDQoD179qiwsNDRwE1NzTp//sJNAVxUmy3rwjm2Qksm+B6Z2JGJHZl4Ig87MrEjEzsy8R1H\nFxKkpKRo0KBBmjFjhiorK7V582Y9++yzeuCBB3TzzTerpKREK1euVFVVlV555RVt2LBBkydPliTl\n5uaqtLRUy5cvV0VFhaZPn664uDhlZmb65MQAAADQsTkqqsHBwVqyZIk6d+6sO++8UzNnztRvfvMb\n3XPPPUpJSdHChQv15ptvKicnRy+//LLmz5+vwYMHS5L69u2roqIiFRcXa8KECaqvr9eiRYt8clIA\nAADo+By/mSomJkYLFy686Lbs7GxlZ2f/5GOzsrK0adMmp4cEAADAFYjPUwAAAICRKKoAAAAwEkUV\nAAAARqKoAgAAwEgUVQAAABiJogoAAAAjUVQBAABgJIoqAAAAjERRBQAAgJEoqgAAADASRRUAAABG\noqgCAADASBRVAAAAGImiCgAAACNRVAEAAGAkiioAAACMRFEFAACAkSiqAAAAMBJFFQAAAEaiqAIA\nAMBIFFUAAAAYiaIKAAAAI1FUAQAAYCSKKgAAAIxEUQUAAICRKKoAAAAwEkUVAAAARqKoAgAAwEgU\nVQAAABiJogoAAAAjUVQBAABgJIoqAAAAjERRBQAAgJEoqgAAADASRRUAAABGoqgCAADASBRVAAAA\nGImiCgAAACNRVAEAAGAkiioAAACMRFEFAACAkSiqAAAAMJLjonro0CHdd999Sk9PV3Z2tlasWOHe\nVl1drUmTJik9PV2jR4/Wli1bPB67detW5eTkKC0tTXl5eaqqqrr8MwAAAEBAclRULctSfn6+evbs\nqfXr1+uJJ57Q0qVL9fbbb0uSHnzwQcXGxqq4uFi33XabpkyZom+++UaSdOTIERUUFCg3N1fFxcWK\niopSQUGB988IAAAAAcFRUa2rq1NSUpJmzZqluLg43XjjjRo+fLh27typTz75RNXV1Zo9e7b69++v\n/Px8paWlad26dZKk1157TSkpKcrLy1N8fLzmzJmjw4cPq6SkxCcnBgAAgI7NUVGNiYnRc889p86d\nO0uSdu7cqU8//VSZmZn67LPPlJycrPDwcPf+GRkZ2rVrlySpvLxcw4YNc2+LiIhQUlKSysrKvHEe\nAAAACDBtfjNVdna27rnnHqWlpWnUqFGqra1VbGysxz7R0dGqqamRJB09etS2vWfPnu7tAAAAwA+F\ntvWBRUVFqqur0xNPPKE//vGPamhoUFhYmMc+YWFhcrlckqTGxsZLbm+tkJAL3To0NFgKauv05gsO\nCrpwjpfQkkXLV5DJxZCJHZl4Ig87MrEjEzsysfN2Fm0uqsnJyZKkadOm6ZFHHtH48eN16tQpj31c\nLpciIiIkSeHh4bZS6nK5FBkZ6ei4kZGdJEnHj3dWUAAX1avCQhQV1aVV+7Zkgu+RiR2Z2JGJJ/Kw\nIxM7MrEjE99xVFSPHTumsrIyjRw50r02YMAAnTt3TjExMaqsrPTYv66uTjExMZKkXr16qba21rY9\nMTHR0cCnTjWoqalZJ0+ekWU5emiHcs7VpBMnvrvkPiEhwYqM7OTOBGRyMWRiRyaeyMOOTOzIxI5M\n7Foy8RZHRbW6ulqFhYXavHmz+3rT3bt3Kzo6WhkZGVqxYoVcLpf7Jf6dO3dq6NChkqTU1FSVlpa6\nv1dDQ4P27NmjwsJCRwM3NTXr/PkLNwVwUW22rAvn2AotmeB7ZGJHJnZk4ok87MjEjkzsyMR3HF1I\nkJKSokGDBmnGjBmqrKzU5s2b9eyzz+qBBx7QsGHD1KdPH02bNk0VFRVatmyZdu/erfHjx0uScnNz\nVVpaquXLl6uiokLTp09XXFycMjMzfXJiAAAA6NgcFdXg4GAtWbJEnTt31p133qmZM2fqN7/5je65\n5x4FBwdr6dKlqq2tVW5urjZu3KjFixerd+/ekqS+ffuqqKhIxcXFmjBhgurr67Vo0SKfnBQAAAA6\nPsdvpoqJidHChQsvuq1fv35atWrVTz42KytLmzZtcnpIAAAAXIH4PAUAAAAYiaIKAAAAI1FUAQAA\nYCSKKgAAAIxEUQUAAICRKKoAAAAwEkUVAAAARqKoAgAAwEgUVQAAABiJogoAAAAjUVQBAABgJIoq\nAAAAjERRBQAAgJEoqgAAADASRRUAAABGoqgCAADASBRVAAAAGImiCgAAACNRVAEAAGAkiioAAACM\nRFEFAACAkSiqAAAAMBJFFQAAAEaiqAIAAMBIFFUAAAAYiaIKAAAAI1FUAQAAYCSKKgAAAIxEUQUA\nAICRKKoAAAAwEkUVAAAARqKoAgAAwEgUVQAAABiJogoAAAAjUVQBAABgJIoqAAAAjERRBQAAgJEo\nqgAAADASRRUAAABGoqgCAADASBRVAAAAGMlRUa2pqdHUqVN13XXX6aabbtIzzzwjl8slSXrqqaeU\nkJCgxMRE99eXX37Z/ditW7cqJydHaWlpysvLU1VVlXfPBAAAAAEl1MnOU6dOVffu3fXKK6/o22+/\n1YwZMxQSEqJHH31U+/fv1yOPPKJx48a59+/atask6ciRIyooKNBDDz2krKwsLVq0SAUFBdqwYYN3\nzwYAAAABo9XPqO7fv1/l5eWaM2eO4uPjlZGRoalTp+qtt96SJFVWViopKUnR0dHuW3h4uCRp7dq1\nSklJUV5enuLj4zVnzhwdPnxYJSUlvjkrAAAAdHitLqoxMTF68cUX1aNHD/eaZVmqr6/X6dOnVVNT\no6uvvvqij/3ss880bNgw9/2IiAglJSWprKys7ZMDAAAgoLW6qHbr1k0jRoxw37csS6tXr9b111+v\n/fv3KygoSEuXLtVNN92kMWPG6M0333Tve/ToUcXGxnp8v549e6qmpsYLpwAAAIBA5Oga1R+aO3eu\n9u3bp3Xr1unzzz9XcHCw4uPjNXHiRO3YsUMzZ85U165dNXLkSDU2NiosLMzj8WFhYe43YjkREnKh\nW4eGBktBbZ3efMFBQRfO8RJasmj5CjK5GDKxIxNP5GFHJnZkYkcmdt7Ook1Fdd68eVq1apVeeOEF\nDRgwQAMGDFB2drYiIyMlSddee60OHjyoV199VSNHjlR4eLitlLpcLvf+TkRGdpIkHT/eWUEBXFSv\nCgtRVFSXVu3bkgm+RyZ2ZGJHJp7Iw45M7MjEjkx8x3FRffLJJ7VmzRrNmzdPI0eOdK//uHT2799f\n27dvlyT16tVLtbW1Htvr6uqUmJjoeOBTpxrU1NSskyfPyLIcP7zDOOdq0okT311yn5CQYEVGdnJn\nAjK5GDKxIxNP5GFHJnZkYkcmdi2ZeIujorpo0SKtWbNGzz//vG655Rb3+sKFC1VWVqaVK1e61/bu\n3atrrrlGkpSamqrS0lL3toaGBu3Zs0eFhYWOB25qatb58xduCuCi2mxZF86xFVoywffIxI5M7MjE\nE3nYkYkdmdiRie+0+kKCyspKLV26VPn5+UpPT1ddXZ37dvPNN6ukpEQrV65UVVWVXnnlFW3YsEGT\nJ0+WJOXm5qq0tFTLly9XRUWFpk+frri4OGVmZvrsxAAAANCxtfoZ1b/+9a9qbm7W0qVLtXTpUkkX\n3vkfFBSkvXv3auHChVqwYIEWLFigvn37av78+Ro8eLAkqW/fvioqKtLTTz+tJUuWaMiQIVq0aJFv\nzggAAAABodVFNT8/X/n5+T+5PTs7W9nZ2T+5PSsrS5s2bXI2HQAAAK5YfJ4CAAAAjERRBQAAgJEo\nqgAAADASRRUAAABGoqgCAADASBRVAAAAGImiCgAAACNRVAEAAGAkiioAAACMRFEFAACAkSiqAAAA\nMBJFFQAAAEaiqAIAAMBIFFUAAAAYiaIKAAAAI1FUAQAAYCSKKgAAAIxEUQUAAICRKKoAAAAwEkUV\nAAAARqKoAgAAwEgUVQAAABiJogoAAAAjUVQBAABgJIoqAAAAjERRBQAAgJEoqgAAADASRRUAAABG\noqgCAADASBRVAAAAGImiCgAAACNRVAEAAGAkiioAAACMRFEFAACAkSiqAAAAMBJFFQAAAEaiqAIA\nAMBIFFUAAAAYiaIKAAAAI1FUAQAAYCSKKgAAAIzkqKjW1NRo6tSpuu6663TTTTfpmWeekcvlkiRV\nV1dr0qRJSk9P1+jRo7VlyxaPx27dulU5OTlKS0tTXl6eqqqqvHcWAAAACDiOiurUqVN19uxZvfLK\nK3ruuef0P//zP1qwYIEk6cEHH1RsbKyKi4t12223acqUKfrmm28kSUeOHFFBQYFyc3NVXFysqKgo\nFRQUeP9sAAAAEDBaXVT379+v8vJyzZkzR/Hx8crIyNDUqVP11ltv6ZNPPlF1dbVmz56t/v37Kz8/\nX2lpaVq3bp0k6bXXXlNKSory8vIUHx+vOXPm6PDhwyopKfHZiQEAAKBja3VRjYmJ0YsvvqgePXp4\nrNfX1+uzzz5TcnKywsPD3esZGRnatWuXJKm8vFzDhg1zb4uIiFBSUpLKysoud34AAAAEqFYX1W7d\numnEiBHu+5ZlafXq1Ro+fLhqa2sVGxvrsX90dLRqamokSUePHrVt79mzp3s7AAAA8GOhbX3g3Llz\ntXfvXq1bt04rV65UWFiYx/awsDD3G60aGxsvud2JkJAL3To0NFgKauPwHUBwUNCFc7yElixavoJM\nLoZM7MjEE3nYkYkdmdiRiZ23s2hTUZ03b55WrVqlF154QQMGDFB4eLhOnjzpsY/L5VJERIQkKTw8\n3FZKXS6XIiMjHR87MrKTJOn48c4KCuCielVYiKKiurRq35ZM8D0ysSMTOzLxRB52ZGJHJnZk4juO\ni+qTTz6pNWvWaN68eRo5cqQkqVevXqqoqPDYr66uTjExMe7ttbW1tu2JiYmOBz51qkFNTc06efKM\nLMvxwzuMc64mnTjx3SX3CQkJVmRkJ3cmIJOLIRM7MvFEHnZkYkcmdmRi15KJtzgqqosWLdKaNWv0\n/PPP65ZbbnGvp6amavny5XK5XO6X+Hfu3KmhQ4e6t5eWlrr3b2ho0J49e1RYWOh44KamZp0/f+Gm\nAC6qzZZ14RxboSUTfI9M7MjEjkw8kYcdmdiRiR2Z+E6rLySorKzU0qVLlZ+fr/T0dNXV1blvmZmZ\n6tOnj6ZNm6aKigotW7ZMu3fv1vjx4yVJubm5Ki0t1fLly1VRUaHp06crLi5OmZmZPjsxAAAAdGyt\nLqp//etf1dzcrKVLlyorK0tZWVm64YYblJWVpeDgYC1evFi1tbXKzc3Vxo0btXjxYvXu3VuS1Ldv\nXxUVFam4uFgTJkxQfX29Fi1a5LOTAgAAQMfX6pf+8/PzlZ+f/5Pb4+LitGrVqp/cnpWVpU2bNjmb\nDgAAAFcsPk8BAAAARqKoAgAAwEgUVQAAABiJogoAAAAjUVQBAABgJIoqAAAAjERRBQAAgJEoqgAA\nADASRRUAAABGoqgCAADASBRVAAAAGImiCgAAACNRVAEAAGAkiioAAACMRFEFAACAkSiqAAAAMBJF\nFQAAAEaiqAIAAMBIFFUAAAAYiaIKAAAAI1FUAQAAYCSKKgAAAIxEUQUAAICRKKoAAAAwEkUVAAAA\nRqKoAgAAwEgUVQAAABiJogoAAAAjUVQBAABgJIoqAAAAjERRBQAAgJEoqgAAADASRRUAAABGoqgC\nAADASBRVAAAAGImiCgAAACNRVAEAAGAkiioAAACMRFEFAACAkSiqAAAAMBJFFQAAAEZqc1F1uVzK\nyclRSUmJe+2pp55SQkKCEhMT3V9ffvll9/atW7cqJydHaWlpysvLU1VV1eVNDwAAgIDVpqLqcrn0\n8MMPq6KiwmN9//79euSRR/Txxx9ry5Yt+vjjjzV+/HhJ0pEjR1RQUKDc3FwVFxcrKipKBQUFl38G\nAAAACEiOi2plZaXuuOMOVVdXX3RbUlKSoqOj3bfw8HBJ0tq1a5WSkqK8vDzFx8drzpw5Onz4sMcz\nsgAAAEALx0V1x44dGj58uNasWSPLstzrp0+fVk1Nja6++uqLPu6zzz7TsGHD3PcjIiKUlJSksrIy\n51MDAAAg4IU6fcBdd9110fX9+/crKChIS5cu1Ycffqju3btr0qRJGjt2rCTp6NGjio2N9XhMz549\nVVNT04axAQAAEOgcF9Wfsn//fgUHBys+Pl4TJ07Ujh07NHPmTHXt2lUjR45UY2OjwsLCPB4TFhYm\nl8vl6DghIReeBA4NDZaCvDW9eYKDgi6c4yW0ZNHyFWRyMWRiRyaeyMOOTOzIxI5M7LydhdeK6tix\nY5Wdna3IyEhJ0rXXXquDBw/q1Vdf1ciRIxUeHm4rpS6Xy71/a0VGdpIkHT/eWUEBXFSvCgtRVFSX\nVu3bkgm+RyZ2ZGJHJp7Iw45M7MjEjkx8x2tFVZKtdPbv31/bt2+XJPXq1Uu1tbUe2+vq6pSYmOjo\nGKdONaipqVknT57RDy6RDTjnXE06ceK7S+4TEhKsyMhO7kxAJhdDJnZk4ok87MjEjkzsyMSuJRNv\n8VpRXbhwocrKyrRy5Ur32t69e3XNNddIklJTU1VaWure1tDQoD179qiwsNDRcZqamnX+/IWbArio\nNlvWhXNshZZM8D0ysSMTOzLxRB52ZGJHJnZk4jteu5Dg5ptvVklJiVauXKmqqiq98sor2rBhgyZP\nnixJys3NVWlpqZYvX66KigpNnz5dcXFxyszM9NYIAAAACCCXVVSDfnCRaEpKihYuXKg333xTOTk5\nevnllzV//nwNHjxYktS3b18VFRWpuLhYEyZMUH19vRYtWnR50wMAACBgXdZL/3v37vW4n52drezs\n7J/cPysrS5s2bbqcQwIAAOAKwecpAAAAwEgUVQAAABiJogoAAAAjUVQBAABgJIoqAAAAjERRBQAA\ngJEoqgAAADASRRUAAABGoqgCAADASBRVAAAAGImiCgAAACNRVAEAAGAkiioAAACMRFEFAACAkSiq\nAAAAMBJFFQAAAEaiqAIAAMBIFFUAAAAYiaIKAAAAI1FUAQAAYCSKKgAAAIxEUQUAAICRKKoAAAAw\nEkUVAAAARqKoAgAAwEgUVQAAABiJogoAAAAjUVQBAABgJIoqAAAAjERRBQAAgJEoqgAAADASRRUA\nAABGoqgCAADASBRVAAAAGImiCgAAACNRVAEAAGAkiioAAACMRFEFAACAkSiqAAAAMBJFFQAAAEai\nqAIAAMBIFFUAAAAYqc1F1eVyKScnRyUlJe616upqTZo0Senp6Ro9erS2bNni8ZitW7cqJydHaWlp\nysvLU1VVVdsnBwAAQEBrU1F1uVx6+OGHVVFR4bFeUFCg2NhYFRcX67bbbtOUKVP0zTffSJKOHDmi\ngoIC5ebmqri4WFFRUSooKLj8MwAAAEBAclxUKysrdccdd6i6utpjfdu2baqqqtLs2bPVv39/5efn\nKy0tTevWrZMkvfbaa0pJSVFeXp7i4+M1Z84cHT582OMZWQAAAKCF46K6Y8cODR8+XGvWrJFlWe71\n8vJyJScnKzw83L2WkZGhXbt2ubcPGzbMvS0iIkJJSUkqKyu7nPkBAAAQoEKdPuCuu+666Hptba1i\nY2M91qKjo1VTUyNJOnr0qG17z5493dsBAACAH3JcVH9KQ0ODwsLCPNbCwsLkcrkkSY2NjZfc3loh\nIReeBA4NDZaCLmNgwwUHBV04x0toyaLlK8jkYsjEjkw8kYcdmdiRiR2Z2Hk7C68V1fDwcJ08edJj\nzeVyKSIiwr39x6XU5XIpMjLS0XEiIztJko4f76ygAC6qV4WFKCqqS6v2bckE3yMTOzKxIxNP5GFH\nJnZkYkcmvuO1otqrVy/bpwDU1dUpJibGvb22tta2PTEx0dFxTp1qUFNTs06ePKMfXCIbcM65mnTi\nxHeX3CckJFiRkZ3cmYBMLoZM7MjEE3nYkYkdmdiRiV1LJt7itaKampqq5cuXy+VyuV/i37lzp4YO\nHereXlpa6t6/oaFBe/bsUWFhoaPjNDU16/z5CzcFcFFttqwL59gKLZnge2RiRyZ2ZOKJPOzIxI5M\n7MjEd7x2IUFmZqb69OmjadOmqaKiQsuWLdPu3bs1fvx4SVJubq5KS0u1fPlyVVRUaPr06YqLi1Nm\nZqa3RgAAAEAAuayiGvSDi0SDg4O1ZMkS1dbWKjc3Vxs3btTixYvVu3dvSVLfvn1VVFSk4uJiTZgw\nQfX19Vq0aNHlTQ8AAICAdVkv/e/du9fjfr9+/bRq1aqf3D8rK0ubNm26nEMCAADgCsHnKQAAAMBI\nFFUAAADmS8hDAAATrUlEQVQYiaIKAAAAI1FUAQAAYCSKKgAAAIxEUQUAAICRKKoAAAAwEkUVAAAA\nRqKoAgAAwEiX9Zep4BvNTedUV/O1ysp2XnK/kJBgRUZ20qlTDWpqavbTdN6RnJyisLCw9h4DAAAY\njKJqoFO1B1V/uoue/M9P23sUn6g/dkhzH5bS0zPaexQAAGAwiqqhukXHqXvvn7X3GAAAAO2Ga1QB\nAABgJIoqAAAAjERRBQAAgJEoqgAAADASRRUAAABGoqgCAADASBRVAAAAGImiCgAAACNRVAEAAGAk\niioAAACMRFEFAACAkSiqAAAAMBJFFQAAAEaiqAIAAMBIFFUAAAAYiaIKAAAAI1FUAQAAYKTQ9h4A\nV57mpvP68st9PvneISHBiozspFOnGtTU1OyTY/w9yckpCgsLa5djAwAQSCiq8Lvvvj2iFW9/rW6f\nnG7vUbyu/tghzX1YSk/PaO9RAADo8CiqaBfdouPUvffP2nsMAABgMK5RBQAAgJEoqgAAADASRRUA\nAABGoqgCAADASBRVAAAAGImiCgAAACNRVAEAAGAkiioAAACMRFEFAACAkbxaVN9//30lJCQoMTHR\n/fWhhx6SJFVXV2vSpElKT0/X6NGjtWXLFm8eGgAAAAHGq39CtaKiQtnZ2XrqqadkWZYkKTw8XJL0\n4IMPKjExUcXFxXr//fc1ZcoU/eUvf1Hv3r29OQIAAAAChFeLamVlpX72s5+pR48eHuvbtm1TdXW1\n1q5dq/DwcOXn52vbtm1at26dpkyZ4s0RAAAAECC8+tJ/ZWWlrrnmGtt6eXm5kpOT3c+uSlJGRoZ2\n7drlzcMDAAAggHi1qB44cEAfffSRfv7zn+uWW27R/Pnzde7cOdXW1io2NtZj3+joaNXU1Hjz8AAA\nAAggXnvp/+uvv1ZjY6PCw8O1YMECVVdX6+mnn1ZjY6MaGhoUFhbmsX9YWJhcLpfj44SEXOjWoaHB\nUpBXRge8KiQk+MLPpyFafmdavoJMfow87MjEjkzsyMTO21l4raj+wz/8g7Zv367IyEhJUkJCgpqb\nm/Xoo4/q9ttv16lTpzz2d7lcioiIcHycyMhOkqTjxzsriKIKA0VGdlJUVJf2HsOm5XcH3yMTT+Rh\nRyZ2ZGJHJr7j1TdTtZTUFvHx8Tp79qx69uypyspKj211dXWKiYlxfIxTpxrU1NSskyfP6P8+WAAw\nyqlTDTpx4rv2HsMtJCRYkZGd3L87IJMfIw87MrEjEzsysWvJxFu8VlQ//vhj/f73v9eHH37oftPU\nnj17FBUVpaFDh+rPf/6zXC6X+xKAnTt3aujQoY6P09TUrPPnL9xEUYWBWn5GTWPqXO2JTDyRhx2Z\n2JGJHZn4jtcuJEhPT1enTp30hz/8QQcOHNDmzZs1b9483X///Ro2bJj69OmjadOmqaKiQsuWLdPu\n3bs1fvx4bx0eAAAAAcZrRbVLly5asWKFTpw4ofHjx2vmzJm68847de+99yo4OFhLly5VbW2tcnNz\ntXHjRi1evJgP+wcAAMBP8uo1qvHx8VqxYsVFt/Xr10+rVq3y5uEAAAAQwPg8BQAAABiJogoAAAAj\nefWlf+BK19x0Xl9+ua+9x/Dg7Y9PSU5Osf0BDwAAfIGiCnjRd98e0Yq3v1a3T0639yg+UX/skOY+\nLKWnZ7T3KACAKwBFFfCybtFx6t77Z+09BgAAHR7XqAIAAMBIFFUAAAAYiaIKAAAAI1FUAQAAYCSK\nKgAAAIxEUQUAAICRKKoAAAAwEkUVAAAARqKoAgAAwEgUVQAAABiJogoAAAAjUVQBAABgJIoqAAAA\njERRBQAAgJEoqgAAADASRRUAAABGoqgCAADASBRVAAAAGImiCgAAACOFtvcAADqO5qbz+vLLfe09\nxmULCQlWZGQnnTrVoKamZvd6cnKKwsLC2nEyAMAPUVQBtNp33x7Rire/VrdPTrf3KF5Xf+yQ5j4s\npadntPcoAID/Q1EF4Ei36Dh17/2z9h4DAHAF4BpVAAAAGImiCgAAACNRVAEAAGAkrlEFAAXOJxr8\nFD7RAEBHRFEFAPGJBgBgIooqAPwfPtEAAMxCUQWAAPf3Lmv4qT+A0JFwaQMQmCiqABDgAvmyBolL\nG4BARlEFgCsAlzUA6Ij4eCoAAAAYiaIKAAAAI/HSPwCgQ/PFZ+Ca9AYz3iiGKxlFFQDQoQXym8V4\noxiudBRVAECHx5vFgMDk12tUXS6XZsyYoWHDhikrK0srV6705+EBAADQgfj1GdU//elP2rNnj1at\nWqXq6mo99thj6tu3r0aNGuXPMQAA6BB8cf1tW3n7ut1z585Jkq666qrL/l7t5e9lwvXFl89vRbWh\noUHr1q3TihUrlJCQoISEBE2ePFmrV6+mqAIAcBGBfP1tzf4Sdf5/vdQtOq69R/EJri/2Dr8V1X37\n9qmpqUlpaWnutYyMDP37v/+7v0YAAKDDCdTrb+uPValbdL+APDd4j9+uUa2trVX37t0VGvp9N46O\njtbZs2d14sQJf40BAACADsKvL/3/+DqNlvsul6vV3yck5EK3Dg0NloK8N59p6o8dau8RfObMyW8k\nWe09hk8E8rlJgX1+nFvHFcjnx7l1XPXHDikkJPNCX7mCtPQ0b/FbUQ0PD7cV0pb7nTp1avX3iYy8\nsG9UVIo2v5rivQEBAABgFL/V/F69eunbb79Vc/P374qrq6tTRESEIiMj/TUGAAAAOgi/FdXExESF\nhoZq165d7rVPP/1UgwYN8tcIAAAA6ED8VlQjIiI0ZswYzZo1S7t379b777+vlStX6re//a2/RgAA\nAEAHEmRZlt+uZG5sbNS//uu/6p133lG3bt00efJkTZw40V+HBwAAQAfi16IKAAAAtNaV9ZkJAAAA\n6DAoqgAAADASRRUAAABGoqgCAADASBRVAAAAGMmooupyuTRjxgwNGzZMWVlZWrly5U/uu2fPHt1x\nxx1KS0vThAkT9MUXX/hxUv9xkkmLTz/9VCNHjvTDdO3DSSYffPCBxo4dq/T0dI0ZM0b//d//7cdJ\n/cdJJhs2bNDPf/5zpaam6q677lJ5ebkfJ/WPtvzeVFdXKz09XSUlJX6Y0P+cZPLAAw8oISFBiYmJ\n7q+bN2/247T+4SSTL7/8Ur/+9a+Vmpqq2267Tdu3b/fjpP7T2kwmTpyohIQE2+0Pf/iDnyf2PSc/\nJ++9955uvfVWpaen6+6779aePXv8OKn/OMnk448/1pgxY5Senq57771XBw4ccHYwyyCzZ8+2xowZ\nY+3du9d67733rCFDhljvvPOObb8zZ85YI0aMsObOnWtVVlZaTz31lDVixAiroaGhHab2rdZm0mLf\nvn3WiBEjrOzsbD9O6V+tzWTv3r3WoEGDrNWrV1uHDh2yVq9ebSUnJ1v79u1rh6l9q7WZlJSUWCkp\nKdbGjRutqqoq65lnnrEyMzOtM2fOtMPUvuP098ayLOu+++6zEhISrB07dvhpSv9yksmoUaOst956\ny6qrq3PfXC6Xnyf2vdZmUl9fb40YMcJ6/PHHrUOHDlkLFy60hg4dah07dqwdpvat1mZy8uRJj5+P\n999/30pJSbG++OKLdpjat1qbyd/+9jdr8ODB1vr1661Dhw5Zs2fPtkaMGGE1Nja2w9S+1dpM/vd/\n/9dKTk62ioqKrAMHDlhz5861srKyHP03x5iieubMGWvw4MFWSUmJe23JkiXWxIkTbfuuXbvWGjly\npMfaqFGjrDfeeMPnc/qTk0wsy7JeffVVKz093RozZkzAFlUnmTz77LPW/fff77F27733Ws8//7zP\n5/QnJ5n85S9/sf7t3/7Nfb++vt4aOHCgVV5e7pdZ/cHp741lWdb69eutu+66K2CLqpNMzp49ayUl\nJVkHDx7054h+5yST//zP/7RGjRrlsTZ+/Hhr8+bNPp/Tn9ryu2NZltXU1GT96le/shYuXOjrEf3O\nSSYrV660cnNz3fdPnz5tDRw40Pr888/9Mqu/OMlk9uzZ1j333OOxduutt1pr1qxp9fGMeel/3759\nampqUlpamnstIyPjoi9LlpeXKyMjw2NtyJAhKisr8/mc/uQkE+nC0+tz584N6D9L6ySTcePG6fe/\n/71t/fTp0z6d0d+cZPKLX/xCv/vd7yRJZ8+e1X/8x3+oZ8+eGjBggN/m9TWnvzcnTpzQ/Pnz9eST\nT8oK0L9/4iSTAwcOKCgoSP369fPniH7nJJOSkhJlZ2d7rK1du1Y33nijz+f0J6e/Oy2Ki4t18uRJ\nTZ482dcj+p2TTLp3766KigqVlpbKsiwVFxerW7duiouL8+fIPuckk6qqKqWmpnqsXXvttY76mjFF\ntba2Vt27d1doaKh7LTo6WmfPntWJEyc89j169KhiY2M91qKjo1VTU+OXWf3FSSaStGjRooC+NlVy\nlkn//v01cOBA9/2//e1v+uSTTzR8+HC/zesPTn9OJGnbtm1KT0/XkiVLNGPGDHXq1Mlf4/qc0zye\neeYZjRs3TvHx8f4c06+cZFJZWamuXbvq0Ucf1Q033KAJEyboww8/9PfIPuckk6qqKkVFRenxxx/X\nDTfcoDvvvFOlpaX+Htnn2vJviSS9+OKLysvLC6h/R1o4yeTWW2/VjTfeqF//+tcaNGiQ5s2bp4UL\nF6pbt27+HtunnGRysW525MiRS/48/ZgxRbWhoUFhYWEeay33XS6Xx3pjY+NF9/3xfh2dk0yuFG3N\n5Pjx4yosLFRGRob++Z//2acz+ltbMhk4cKBef/11TZ06VY899lhAvaHKSR5bt25VWVmZHnzwQb/N\n1x6cZLJ//36dPXtWWVlZWrFihW666SY98MADAfeGVSeZnDlzRi+++KJiY2P14osvaujQobrvvvsC\n7smRtvxb8sknn+jo0aOaMGGCz+drD04y+fbbb1VXV6dZs2Zp7dq1Gjt2rKZNm6bjx4/7bV5/cJLJ\nrbfeqk2bNumDDz5QU1OT3njjDX3++ec6d+5cq48X+vd38Y/w8HDbCbbc//H/pf3UvhEREb4d0s+c\nZHKlaEsmdXV1mjRpkoKCgrRgwQKfz+hvbcmkR48e6tGjhxISErRr1y69+uqrGjx4sM9n9YfW5nH2\n7FnNmjVLTzzxhO0f3UDj5GdkypQp+u1vf+t+FmjgwIH6/PPPtWbNGs2ePds/A/uBk0xCQkKUmJio\nKVOmSJISEhK0ZcsWrV+/Xvn5+f4Z2A/a8m/Ju+++q6ysLEVGRvp8vvbgJJNnn31WAwcO1F133SVJ\nmj17tn75y1/q9ddfD6jLIpxkkpWVpSlTpqiwsFDNzc267rrrNHbsWNXX17f6eMY8o9qrVy99++23\nam5udq/V1dUpIiLC9gvQq1cv1dbWeqzV1dUpJibGL7P6i5NMrhROM6mpqdHdd9+tpqYmrVq1SlFR\nUf4c1y+cZLJ7927bx6XEx8c7ehnGdK3No7y8XNXV1SosLFR6errS09MlSffff7+eeOIJf4/tU05/\nb378UmV8fLyOHj3q8zn9yUkmMTEx6t+/v8fa1VdfrSNHjvhlVn9py39zPvroo4C+5MxJJl988YUS\nEhLc94OCgpSQkKCvv/7ab/P6g9Ofk9/97ncqLS3VRx99pD//+c86ffq0+vbt2+rjGVNUExMTFRoa\nql27drnXPv30Uw0aNMi2b2pqqu1C3NLSUo8LewOBk0yuFE4yaWho0OTJk3XVVVdp9erV6tmzpz9H\n9Rsnmaxbt07z58/3WPviiy8C6vrM1uaRmpqqd999V+vXr9eGDRu0YcMGSdLTTz+tqVOn+nVmX3Py\nMzJ9+nTNmDHDY23fvn265pprfD6nPznJJC0tTfv27fNY279/v6P/2HYETv+bc+LECVVVVWnIkCH+\nGtHvnGQSGxuriooKj7UDBw7oH//xH30+pz85yeTtt9/WH//4R1111VXq0aOHGhsbtX37dl133XWt\nP2BbP57AFx5//HFr9OjRVnl5ufXee+9ZGRkZ1nvvvWdZlmXV1ta6P4usvr7euv76662nn37aqqio\nsJ588knrhhtuCMjPUW1tJj/0+uuvB+zHU1lW6zN57rnnrLS0NKu8vNyqra113+rr69tzfJ9obSZf\nfPGFlZycbL300kvWwYMHrQULFlhDhgyxampq2nN8r2vL741lWdbAgQMD8uOpLKv1mbz77rvWoEGD\nrDfeeMP66quvrKKiIistLc06fPhwe47vE63N5PDhw1Z6erpVVFRkffXVV9YLL7wQkL83luXsd2f7\n9u1Wampqe43qN63N5O2337ZSU1OtN9980/rqq6+sefPmWcOGDQvIz9ttbSaff/65NXjwYOvdd9+1\nDhw4YD344IPW7bff7uhYRhXVhoYGa9q0aVZ6erp14403Wi+99JJ728CBAz0+J7W8vNwaN26clZqa\nat1xxx3W3r1722Nkn3OSSYtAL6qtzeQXv/iFlZCQYLtNmzatvUb3GSc/Jx988IGVk5NjpaamWuPH\nj7d27drVHiP7VFt+byzLCtjPUbUsZ5msXbvWGjVqlDV48GDr9ttvtz799NP2GNnnnGRSWlpqjRs3\nzho8eLA1btw4MrEuFLMbbrihPcb0KyeZrFu3zvrlL39pDRkyxLr77rvpJtb3nSQjI8MqLCy0amtr\nHR0ryLIC9IMDAQAA0KEZc40qAAAA8EMUVQAAABiJogoAAAAjUVQBAABgJIoqAAAAjERRBQAAgJEo\nqgAAADASRRUAAABGoqgCAADASBRVAAAAGImiCgAAACP9fwgg89AQ1shhAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10e584bd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">25% 102\n",
      ">50% 25\n",
      ">75% 3\n",
      ">90% 0\n",
      ">95% 0\n"
     ]
    }
   ],
   "source": [
    "plt.hist(table_forehead.pa(inplace=False).sum(axis=\"observation\")/len(table_forehead.ids()))\n",
    "plt.show()\n",
    "print \">25%\", np.sum(table_forehead.pa(inplace=False).sum(axis=\"observation\")/len(table_forehead.ids())>.25)\n",
    "print \">50%\", np.sum(table_forehead.pa(inplace=False).sum(axis=\"observation\")/len(table_forehead.ids())>.50)\n",
    "print \">75%\", np.sum(table_forehead.pa(inplace=False).sum(axis=\"observation\")/len(table_forehead.ids())>.75)\n",
    "print \">90%\", np.sum(table_forehead.pa(inplace=False).sum(axis=\"observation\")/len(table_forehead.ids())>.90)\n",
    "print \">95%\", np.sum(table_forehead.pa(inplace=False).sum(axis=\"observation\")/len(table_forehead.ids())>.99)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqoAAAHcCAYAAAAeFogrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3X9UVPed//HXAA4YdY6IgNbKMWISQBEIQo7tkmxYYrdd\n8ceiOfnZYGLYVsSck2qr2XW1JqlJ1aQqyh6N5Wz9ka9RulHTXdsm2U02apQIig2as2BzBIPIpCZi\nOngF5vtHDrOhJFkGhuED9/k4x+Ph3jvOe847mOeZGUaH1+v1CgAAADBMSH8PAAAAAHwZQhUAAABG\nIlQBAABgJEIVAAAARiJUAQAAYCRCFQAAAEYiVAEAAGAkQhUAAABGIlQBAABgpB6HakFBgVasWOH7\nur6+XgsWLFBaWppmzpypI0eOdLr+6NGjys3NVWpqqvLz81VXV9fzqQEAADDo9ShUf/Ob3+jtt9/u\ndKywsFAxMTEqKyvTrFmztHjxYl26dEmS1NDQoMLCQuXl5amsrEyRkZEqLCzs/fQAAAAYtPwO1U8/\n/VTr1q3T1KlTfceOHTumuro6rVmzRhMnTlRBQYFSU1O1f/9+SdIrr7yi5ORk5efnKz4+XmvXrtXF\nixdVXl4euEcCAACAQcXvUH3++ec1e/ZsxcfH+45VVVVp8uTJCg8P9x1LT0/XqVOnfOczMjJ85yIi\nIpSUlKTKysrezA4AAIBBzK9QPXbsmE6ePNnlZfumpibFxMR0OhYVFaXGxkZJ0uXLl7ucHz16tO88\nAAAA8Je6HaqWZWn16tVatWqVnE5np3Mej6fLMafTKcuyJEktLS1fex4AAAD4S90O1c2bN2vKlCn6\n1re+1eVceHh4l+i0LEsRERHdOt9dXq/Xr+sBAAAwcIV198J///d/18cff6y0tDRJ0o0bNyRJv/3t\nb/WDH/xANTU1na53u92Kjo6WJMXGxqqpqanL+cTERL+GdTgcunrVo7a2dr9uh4EnNDRELtdQ9m0T\n7Nte2Le9sG976dh3oHQ7VHft2qXW1lbf1+vWrZMkLVu2TBcvXtS2bdtkWZbvJf6TJ09q2rRpkqSU\nlBRVVFT4buvxeFRdXa2ioiK/B25ra1drK/+h2wX7thf2bS/s217YN3qi2y/9jx07VuPHj/f9GjZs\nmIYNG6bx48crMzNTY8eO1fLly1VTU6Nt27bpzJkzmjdvniQpLy9PFRUV2r59u2pqarRixQrFxcUp\nMzOzzx4YAAAABraA/BOqISEh2rp1q5qampSXl6dDhw5py5YtGjNmjCRp3Lhx2rx5s8rKyjR//nw1\nNzeruLg4EHcNAACAQcrhHWA/oXTlyme8dGADYWEhiowcxr5tgn3bC/u2F/ZtLx37DpSAPKMKAAAA\nBBqhCgAAACMRqgAAADASoQoAAAAjEaoAAAAwEqEKAAAAIxGqAAAAMBKhCgAAACMRqgAAADASoQoA\nAAAjEaoAAAAwEqEKAAAAIxGqAAAAMBKhCgAAACMRqgAAADASoQoAAAAjEaoAAAAwEqEKAAAAIxGq\nAAAAMBKhCgAAACMRqgAAADASoQoAAAAjEaoAAAAwEqEKAAAAIxGqAAAAMBKhCgAAACMRqgAAADAS\noQoAAAAjEaoAAAAwEqEKAAAAIxGqAAAAMBKhCgAAACMRqgAAADASoQoAAAAjhfX3AD217Zc7dfb8\nxf4eI6i+Ee3SsicW9fcYAAAAQTFgQ/WP9Zf06fDM/h4jqNouvdvfIwAAAAQNL/0DAADASIQqAAAA\njESoAgAAwEiEKgAAAIxEqAIAAMBIfofqhQsX9NhjjyktLU3Z2dnasWOH79wzzzyjhIQEJSYm+n7f\nvXu37/zRo0eVm5ur1NRU5efnq66uLjCPAgAAAIOOXx9P5fV6VVBQoJSUFB04cEAffvihnnzySY0Z\nM0Z/93d/p/Pnz2vp0qWaO3eu7zbDhw+XJDU0NKiwsFBPPPGEsrKyVFxcrMLCQh08eDCwjwgAAACD\ngl/PqLrdbiUlJWnVqlWKi4vTnXfeqenTp+vkyZOSpNraWiUlJSkqKsr3Kzw8XJK0b98+JScnKz8/\nX/Hx8Vq7dq0uXryo8vLywD8qAAAADHh+hWp0dLReeOEF3XTTTZKkkydPqry8XHfccYeuXbumxsZG\nTZgw4Utve/r0aWVkZPi+joiIUFJSkiorK3s+PQAAAAatHv/LVNnZ2WpoaNBf//Vfa8aMGaqqqpLD\n4VBJSYnefvttjRw5UgsWLNCcOXMkSZcvX1ZMTEynP2P06NFqbGzs3SMAAADAoNTjUN28ebPcbrdW\nrVqlZ599VlOmTFFISIji4+P18MMP68SJE1q5cqWGDx+unJwctbS0yOl0dvoznE6nLMvy635DQz9/\nEtjhcPR09AHLEeJQWJg9PqihY88dv2NwY9/2wr7thX3bS6D33ONQnTx5siRpxYoVWrZsmX7yk58o\nOztbLpdLknTrrbfqww8/1Msvv6ycnByFh4d3iVLLsnzXd5fLNVSS5HSGSd6eTj8wDQkLVWTksP4e\nI6g69g17YN/2wr7thX2jJ/wK1Y8//liVlZXKycnxHZs0aZJu3Lihzz77TCNHjux0/cSJE3X8+HFJ\nUmxsrJqamjqdd7vdSkxM9Gvgq1c9amtrl2W1SkP8uumAd6O1TVeufNbfYwRFaGiIXK6hvn1jcGPf\n9sK+7YV920vHvgPFr1Ctr69XUVGR3nrrLd/7Tc+cOaNRo0bpV7/6lSorK1VaWuq7/uzZs7r55psl\nSSkpKaqoqPCd83g8qq6uVlFRkV8Dt7W1q7W1XV6vzZ5OleRt96q11V7f5B37hj2wb3th3/bCvtET\nfr2RIDk5WVOmTNFTTz2l2tpavfXWW1q/fr1++MMf6u6771Z5eblKS0tVV1enPXv26ODBg1q4cKEk\nKS8vTxUVFdq+fbtqamq0YsUKxcXFKTMzs08eGAAAAAY2v0I1JCREW7du1U033aT77rtPK1eu1Pe/\n/3099NBDSk5O1qZNm/Tqq68qNzdXu3fv1oYNGzR16lRJ0rhx47R582aVlZVp/vz5am5uVnFxcZ88\nKAAAAAx8fv8wVXR0tDZt2vSl57Kzs5Wdnf2Vt83KytLhw4f9vUsAAADYEJ8VAQAAACMRqgAAADAS\noQoAAAAjEaoAAAAwEqEKAAAAIxGqAAAAMBKhCgAAACMRqgAAADASoQoAAAAjEaoAAAAwEqEKAAAA\nIxGqAAAAMBKhCgAAACMRqgAAADASoQoAAAAjEaoAAAAwEqEKAAAAIxGqAAAAMBKhCgAAACMRqgAA\nADASoQoAAAAjEaoAAAAwEqEKAAAAIxGqAAAAMBKhCgAAACMRqgAAADASoQoAAAAjEaoAAAAwEqEK\nAAAAIxGqAAAAMBKhCgAAACMRqgAAADASoQoAAAAjEaoAAAAwEqEKAAAAIxGqAAAAMBKhCgAAACMR\nqgAAADASoQoAAAAjEaoAAAAwEqEKAAAAIxGqAAAAMJLfoXrhwgU99thjSktLU3Z2tnbs2OE7V19f\nrwULFigtLU0zZ87UkSNHOt326NGjys3NVWpqqvLz81VXV9f7RwAAAIBBya9Q9Xq9Kigo0OjRo3Xg\nwAGtXr1aJSUl+s1vfiNJWrRokWJiYlRWVqZZs2Zp8eLFunTpkiSpoaFBhYWFysvLU1lZmSIjI1VY\nWBj4RwQAAIBBwa9QdbvdSkpK0qpVqxQXF6c777xT06dP18mTJ/Xuu++qvr5ea9as0cSJE1VQUKDU\n1FTt379fkvTKK68oOTlZ+fn5io+P19q1a3Xx4kWVl5f3yQMDAADAwOZXqEZHR+uFF17QTTfdJEk6\nefKk3nvvPWVmZur06dOaPHmywsPDfdenp6fr1KlTkqSqqiplZGT4zkVERCgpKUmVlZWBeBwAAAAY\nZHr8w1TZ2dl66KGHlJqaqhkzZqipqUkxMTGdromKilJjY6Mk6fLly13Ojx492nceAAAA+KKwnt5w\n8+bNcrvdWr16tX72s5/J4/HI6XR2usbpdMqyLElSS0vL157vrtDQz9va4XD0dPQByxHiUFiYPT6o\noWPPHb9jcGPf9sK+7YV920ug99zjUJ08ebIkafny5Vq6dKnmzZunq1evdrrGsixFRERIksLDw7tE\nqWVZcrlcft2vyzVUkuR0hknenk4/MA0JC1Vk5LD+HiOoOvYNe2Df9sK+7YV9oyf8CtWPP/5YlZWV\nysnJ8R2bNGmSbty4oejoaNXW1na63u12Kzo6WpIUGxurpqamLucTExP9GvjqVY/a2tplWa3SEL9u\nOuDdaG3TlSuf9fcYQREaGiKXa6hv3xjc2Le9sG97Yd/20rHvQPErVOvr61VUVKS33nrL937TM2fO\nKCoqSunp6dqxY4csy/K9xH/y5ElNmzZNkpSSkqKKigrfn+XxeFRdXa2ioiK/Bm5ra1dra7u8Xps9\nnSrJ2+5Va6u9vsk79g17YN/2wr7thX2jJ/x6I0FycrKmTJmip556SrW1tXrrrbe0fv16/fCHP1RG\nRobGjh2r5cuXq6amRtu2bdOZM2c0b948SVJeXp4qKiq0fft21dTUaMWKFYqLi1NmZmafPDAAAAAM\nbH6FakhIiLZu3aqbbrpJ9913n1auXKnvf//7euihhxQSEqKSkhI1NTUpLy9Phw4d0pYtWzRmzBhJ\n0rhx47R582aVlZVp/vz5am5uVnFxcZ88KAAAAAx8fv8wVXR0tDZt2vSl58aPH6+dO3d+5W2zsrJ0\n+PBhf+8SAAAANsRnRQAAAMBIhCoAAACMRKgCAADASIQqAAAAjESoAgAAwEiEKgAAAIxEqAIAAMBI\nhCoAAACMRKgCAADASIQqAAAAjESoAgAAwEiEKgAAAIxEqAIAAMBIhCoAAACMRKgCAADASIQqAAAA\njESoAgAAwEiEKgAAAIxEqAIAAMBIhCoAAACMRKgCAADASIQqAAAAjESoAgAAwEiEKgAAAIxEqAIA\nAMBIhCoAAACMRKgCAADASIQqAAAAjESoAgAAwEiEKgAAAIxEqAIAAMBIhCoAAACMRKgCAADASIQq\nAAAAjESoAgAAwEiEKgAAAIxEqAIAAMBIhCoAAACMRKgCAADASIQqAAAAjESoAgAAwEh+hWpjY6OW\nLFmiO+64Q3fddZeee+45WZYlSXrmmWeUkJCgxMRE3++7d+/23fbo0aPKzc1Vamqq8vPzVVdXF9hH\nAgAAgEElzJ+LlyxZopEjR2rPnj365JNP9NRTTyk0NFTLli3T+fPntXTpUs2dO9d3/fDhwyVJDQ0N\nKiws1BNPPKGsrCwVFxersLBQBw8eDOyjAQAAwKDR7WdUz58/r6qqKq1du1bx8fFKT0/XkiVL9Npr\nr0mSamtrlZSUpKioKN+v8PBwSdK+ffuUnJys/Px8xcfHa+3atbp48aLKy8v75lEBAABgwOt2qEZH\nR+ull17SqFGjfMe8Xq+am5t17do1NTY2asKECV9629OnTysjI8P3dUREhJKSklRZWdnzyQEAADCo\ndTtUR4wYoW9/+9u+r71er3bt2qVvfetbOn/+vBwOh0pKSnTXXXdp9uzZevXVV33XXr58WTExMZ3+\nvNGjR6uxsTEADwEAAACDkV/vUf2in//85zp37pz279+vP/zhDwoJCVF8fLwefvhhnThxQitXrtTw\n4cOVk5OjlpYWOZ3OTrd3Op2+H8TyR2jo523tcDh6OvqA5QhxKCzMHh/U0LHnjt8xuLFve2Hf9sK+\n7SXQe+5RqK5bt047d+7UL37xC02aNEmTJk1Sdna2XC6XJOnWW2/Vhx9+qJdfflk5OTkKDw/vEqWW\nZfmu94fLNVSS5HSGSd6eTD9wDQkLVWTksP4eI6g69g17YN/2wr7thX2jJ/wO1aefflp79+7VunXr\nlJOT4zv+l9E5ceJEHT9+XJIUGxurpqamTufdbrcSExP9HvjqVY/a2tplWa3SEL9vPqDdaG3TlSuf\n9fcYQREaGiKXa6hv3xjc2Le9sG97Yd/20rHvQPErVIuLi7V37169+OKLuueee3zHN23apMrKSpWW\nlvqOnT17VjfffLMkKSUlRRUVFb5zHo9H1dXVKioq8nvgtrZ2tba2y+u12dOpkrztXrW22uubvGPf\nsAf2bS/s217YN3qi228kqK2tVUlJiQoKCpSWlia32+37dffdd6u8vFylpaWqq6vTnj17dPDgQS1c\nuFCSlJeXp4qKCm3fvl01NTVasWKF4uLilJmZ2WcPDAAAAANbt59RfeONN9Te3q6SkhKVlJRI+vwn\n/x0Oh86ePatNmzZp48aN2rhxo8aNG6cNGzZo6tSpkqRx48Zp8+bNevbZZ7V161bdfvvtKi4u7ptH\nBAAAgEGh26FaUFCggoKCrzyfnZ2t7OzsrzyflZWlw4cP+zcdAAAAbIvPigAAAICRCFUAAAAYiVAF\nAACAkQhVAAAAGIlQBQAAgJEIVQAAABiJUAUAAICRCFUAAAAYiVAFAACAkQhVAAAAGIlQBQAAgJEI\nVQAAABiJUAUAAICRCFUAAAAYiVAFAACAkQhVAAAAGIlQBQAAgJEIVQAAABiJUAUAAICRCFUAAAAY\niVAFAACAkQhVAAAAGIlQBQAAgJEIVQAAABiJUAUAAICRCFUAAAAYiVAFAACAkQhVAAAAGIlQBQAA\ngJEIVQAAABiJUAUAAICRCFUAAAAYiVAFAACAkQhVAAAAGIlQBQAAgJEIVQAAABiJUAUAAICRCFUA\nAAAYiVAFAACAkQhVAAAAGIlQBQAAgJH8CtXGxkYtWbJEd9xxh+666y4999xzsixLklRfX68FCxYo\nLS1NM2fO1JEjRzrd9ujRo8rNzVVqaqry8/NVV1cXuEcBAACAQcevUF2yZImuX7+uPXv26IUXXtB/\n/ud/auPGjZKkRYsWKSYmRmVlZZo1a5YWL16sS5cuSZIaGhpUWFiovLw8lZWVKTIyUoWFhYF/NAAA\nABg0uh2q58+fV1VVldauXav4+Hilp6dryZIleu211/Tuu++qvr5ea9as0cSJE1VQUKDU1FTt379f\nkvTKK68oOTlZ+fn5io+P19q1a3Xx4kWVl5f32QMDAADAwNbtUI2OjtZLL72kUaNGdTre3Nys06dP\na/LkyQoPD/cdT09P16lTpyRJVVVVysjI8J2LiIhQUlKSKisrezs/AAAABqluh+qIESP07W9/2/e1\n1+vVrl27NH36dDU1NSkmJqbT9VFRUWpsbJQkXb58ucv50aNH+84DAAAAfymspzf8+c9/rrNnz2r/\n/v0qLS2V0+nsdN7pdPp+0KqlpeVrz/sjNPTztnY4HD2cfOByhDgUFmaPD2ro2HPH7xjc2Le9sG97\nYd/2Eug99yhU161bp507d+oXv/iFJk2apPDwcH366aedrrEsSxEREZKk8PDwLlFqWZZcLpff9+1y\nDZUkOZ1hkrcn0w9cQ8JCFRk5rL/HCKqOfcMe2Le9sG97Yd/oCb9D9emnn9bevXu1bt065eTkSJJi\nY2NVU1PT6Tq3263o6Gjf+aampi7nExMT/R746lWP2traZVmt0hC/bz6g3Wht05Urn/X3GEERGhoi\nl2uob98Y3Ni3vbBve2Hf9tKx70DxK1SLi4u1d+9evfjii7rnnnt8x1NSUrR9+3ZZluV7if/kyZOa\nNm2a73xFRYXveo/Ho+rqahUVFfk9cFtbu1pb2+X12uzpVEnedq9aW+31Td6xb9gD+7YX9m0v7Bs9\n0e03EtTW1qqkpEQFBQVKS0uT2+32/crMzNTYsWO1fPly1dTUaNu2bTpz5ozmzZsnScrLy1NFRYW2\nb9+umpoarVixQnFxccrMzOyzBwYAAICBrduh+sYbb6i9vV0lJSXKyspSVlaW/uqv/kpZWVkKCQnR\nli1b1NTUpLy8PB06dEhbtmzRmDFjJEnjxo3T5s2bVVZWpvnz56u5uVnFxcV99qAAAAAw8HX7pf+C\nggIVFBR85fm4uDjt3LnzK89nZWXp8OHD/k0HAAAA2+KzIgAAAGAkQhUAAABGIlQBAABgJEIVAAAA\nRiJUAQAAYCRCFQAAAEYiVAEAAGAkQhUAAABGIlQBAABgJEIVAAAARiJUAQAAYCRCFQAAAEYiVAEA\nAGAkQhUAAABGIlQBAABgJEIVAAAARiJUAQAAYCRCFQAAAEYiVAEAAGAkQhUAAABGIlQBAABgJEIV\nAAAARiJUAQAAYCRCFQAAAEYiVAEAAGAkQhUAAABGIlQBAABgJEIVAAAARiJUAQAAYCRCFQAAAEYi\nVAEAAGAkQhUAAABGIlQBAABgJEIVAAAARiJUAQAAYCRCFQAAAEYiVAEAAGAkQhUAAABGIlQBAABg\nJEIVAAAARiJUAQAAYKQeh6plWcrNzVV5ebnv2DPPPKOEhAQlJib6ft+9e7fv/NGjR5Wbm6vU1FTl\n5+errq6ud9MDAABg0OpRqFqWpSeffFI1NTWdjp8/f15Lly7VO++8oyNHjuidd97RvHnzJEkNDQ0q\nLCxUXl6eysrKFBkZqcLCwt4/AgAAAAxKfodqbW2t7r33XtXX13/puaSkJEVFRfl+hYeHS5L27dun\n5ORk5efnKz4+XmvXrtXFixc7PSMLAAAAdPA7VE+cOKHp06dr79698nq9vuPXrl1TY2OjJkyY8KW3\nO336tDIyMnxfR0REKCkpSZWVlf5PDQAAgEEvzN8b3H///V96/Pz583I4HCopKdHbb7+tkSNHasGC\nBZozZ44k6fLly4qJiel0m9GjR6uxsbEHYwMAAGCw8ztUv8r58+cVEhKi+Ph4Pfzwwzpx4oRWrlyp\n4cOHKycnRy0tLXI6nZ1u43Q6ZVmWX/cTGvr5k8AOhyNQow8YjhCHwsLs8UENHXvu+B2DG/u2F/Zt\nL+zbXgK954CF6pw5c5SdnS2XyyVJuvXWW/Xhhx/q5ZdfVk5OjsLDw7tEqWVZvuu7y+UaKklyOsMk\n7/9x8SAzJCxUkZHD+nuMoOrYN+yBfdsL+7YX9o2eCFioSuoSnRMnTtTx48clSbGxsWpqaup03u12\nKzEx0a/7uHrVo7a2dllWqzSkd/MONDda23Tlymf9PUZQhIaGyOUa6ts3Bjf2bS/s217Yt7107DtQ\nAhaqmzZtUmVlpUpLS33Hzp49q5tvvlmSlJKSooqKCt85j8ej6upqFRUV+XU/bW3tam1t7/SDXHbh\nbfeqtdVe3+Qd+4Y9sG97Yd/2wr7REwF7I8Hdd9+t8vJylZaWqq6uTnv27NHBgwe1cOFCSVJeXp4q\nKiq0fft21dTUaMWKFYqLi1NmZmagRgAAAMAg0qtQ/eIPNCUnJ2vTpk169dVXlZubq927d2vDhg2a\nOnWqJGncuHHavHmzysrKNH/+fDU3N6u4uLh30wMAAGDQ6tVL/2fPnu30dXZ2trKzs7/y+qysLB0+\nfLg3dwkAAACb4LMiAAAAYCRCFQAAAEYiVAEAAGAkQhUAAABGIlQBAABgJEIVAAAARiJUAQAAYCRC\nFQAAAEYiVAEAAGAkQhUAAABGIlQBAABgJEIVAAAARiJUAQAAYCRCFQAAAEYiVAEAAGAkQhUAAABG\nIlQBAABgJEIVAAAARiJUAQAAYCRCFQAAAEYiVAEAAGAkQhUAAABGIlQBAABgJEIVAAAARiJUAQAA\nYCRCFQAAAEYiVAEAAGAkQhUAAABGIlQBAABgJEIVAAAARiJUAQAAYCRCFQAAAEYiVAEAAGAkQhUA\nAABGIlQBAABgJEIVAAAARiJUAQAAYCRCFQAAAEYiVAEAAGAkQhUAAABGIlQBAABgJEIVAAAARupx\nqFqWpdzcXJWXl/uO1dfXa8GCBUpLS9PMmTN15MiRTrc5evSocnNzlZqaqvz8fNXV1fV8cgAAAAxq\nPQpVy7L05JNPqqamptPxwsJCxcTEqKysTLNmzdLixYt16dIlSVJDQ4MKCwuVl5ensrIyRUZGqrCw\nsPePAAAAAIOS36FaW1ure++9V/X19Z2OHzt2THV1dVqzZo0mTpyogoICpaamav/+/ZKkV155RcnJ\nycrPz1d8fLzWrl2rixcvdnpGFgAAAOjgd6ieOHFC06dP1969e+X1en3Hq6qqNHnyZIWHh/uOpaen\n69SpU77zGRkZvnMRERFKSkpSZWVlb+YHAADAIBXm7w3uv//+Lz3e1NSkmJiYTseioqLU2NgoSbp8\n+XKX86NHj/adBwAAAL7I71D9Kh6PR06ns9Mxp9Mpy7IkSS0tLV97vrtCQz9/EtjhcPRi2oHJEeJQ\nWJg9PqihY88dv2NwY9/2wr7thX3bS6D3HLBQDQ8P16efftrpmGVZioiI8J3/yyi1LEsul8uv+3G5\nhkqSnM4wyft/XDzIDAkLVWTksP4eI6g69g17YN/2wr7thX2jJwIWqrGxsV0+BcDtdis6Otp3vqmp\nqcv5xMREv+7n6lWP2traZVmt0pDezTzQ3Ght05Urn/X3GEERGhoil2uob98Y3Ni3vbBve2Hf9tKx\n70AJWKimpKRo+/btsizL9xL/yZMnNW3aNN/5iooK3/Uej0fV1dUqKiry637a2trV2tre6Qe57MLb\n7lVrq72+yTv2DXtg3/bCvu2FfaMnAvZGgszMTI0dO1bLly9XTU2Ntm3bpjNnzmjevHmSpLy8PFVU\nVGj79u2qqanRihUrFBcXp8zMzECNAAAAgEGkV6H6xR9oCgkJ0datW9XU1KS8vDwdOnRIW7Zs0Zgx\nYyRJ48aN0+bNm1VWVqb58+erublZxcXFvZseAAAAg1avXvo/e/Zsp6/Hjx+vnTt3fuX1WVlZOnz4\ncG/uEgAAADbBZ0UAAADASIQqAAAAjESoAgAAwEiEKgAAAIxEqAIAAMBIhCoAAACMRKgCAADASIQq\nAAAAjESoAgAAwEiEKgAAAIxEqAIAAMBIhCoAAACMRKgCAADASIQqAAAAjESoAgAAwEiEKgAAAIxE\nqAIAAMBIhCoAAACMRKgCAADASIQqAAAAjESoAgAAwEiEKgAAAIxEqAIAAMBIhCoAAACMRKgCAADA\nSIQqAAAAjESoAgAAwEiEKgAAAIxEqAIAAMBIhCoAAACMRKgCAADASIQqAAAAjESoAgAAwEiEKgAA\nAIxEqAKxDEblAAAOeUlEQVQAAMBIhCoAAACMRKgCAADASIQqAAAAjESoAgAAwEiEKgAAAIxEqAIA\nAMBIAQ3V119/XQkJCUpMTPT9/sQTT0iS6uvrtWDBAqWlpWnmzJk6cuRIIO8aAAAAg0xYIP+wmpoa\nZWdn65lnnpHX65UkhYeHS5IWLVqkxMRElZWV6fXXX9fixYv1H//xHxozZkwgRwAAAMAgEdBQra2t\n1S233KJRo0Z1On7s2DHV19dr3759Cg8PV0FBgY4dO6b9+/dr8eLFgRwBAAAAg0RAX/qvra3VzTff\n3OV4VVWVJk+e7Ht2VZLS09N16tSpQN49AAAABpGAhuof//hH/fd//7e+853v6J577tGGDRt048YN\nNTU1KSYmptO1UVFRamxsDOTdAwAAYBAJ2Ev/H330kVpaWhQeHq6NGzeqvr5ezz77rFpaWuTxeOR0\nOjtd73Q6ZVmW3/cTGvp5WzscjoDMPZA4QhwKC7PHBzV07Lnjdwxu7Nte2Le9sG97CfSeAxaq3/jG\nN3T8+HG5XC5JUkJCgtrb27Vs2TL9/d//va5evdrpesuyFBER4ff9uFxDJUlOZ5jk7f3cA8mQsFBF\nRg7r7zGCqmPfsAf2bS/s217YN3oioD9M1RGpHeLj43X9+nWNHj1atbW1nc653W5FR0f7fR9Xr3rU\n1tYuy2qVhvRq3AGn5bqlN954u7/HCIqQEIeGD4/QzTffqtDQgP5nCgOFhobI5Rrq+/7G4Ma+7YV9\n20vHvgMlYAXwzjvv6Ec/+pHefvtt3w9NVVdXKzIyUtOmTdMvf/lLWZblewvAyZMnNW3aNL/vp62t\nXa2t7b6Pv7KTK+5G/Wjdfo2IiuvvUYKi+eML2rBsnqZOTevvURAkHd/fsAf2bS/sGz0RsFBNS0vT\n0KFD9Y//+I8qLCzUhQsXtG7dOj3++OPKyMjQ2LFjtXz5ci1atEhvvvmmzpw5o+eeey5Qd28bI6Li\nNHLMLf09BgAAQJ8L2Dtehw0bph07dujKlSuaN2+eVq5cqfvuu0+PPvqoQkJCVFJSoqamJuXl5enQ\noUPasmULH/YPAACArxTQN//Fx8drx44dX3pu/Pjx2rlzZyDvDgAAAIMYnxUBAAAAIxGqAAAAMBKh\nCgAAACMRqgAAADASoQoAAAAjEaoAAAAwEqEKAAAAIxGqAAAAMBKhCgAAACMRqgAAADASoQoAAAAj\nEaoAAAAwEqEKAAAAIxGqAAAAMBKhCgAAACMRqgAAADASoQoAAAAjEaoAAAAwEqEKAAAAIxGqAAAA\nMFJYfw8A4H9ZlqX33z/T32MEXUpKiqRh/T0GAMAwhCpgkPffP6Mfv/BrjYiK6+9Rgqb54wvasCxE\nsbF39vcoAADDEKqAYUZExWnkmFv6ewwAAPod71EFAACAkQhVAAAAGIlQBQAAgJF4jyqM1d7WqnPn\nzqqtrb2/RwmaDz44198jAABgDEIVxvrskwZtP/SRRkQ19/coQdN4vlyxEzP6ewwAAIxAqMJodvsJ\n+OaP6/p7BAAAjMF7VAEAAGAkQhUAAABG4qV/AECf++I/DxwaGiKXa6iuXvUM+h+WnDw5WU6ns7/H\nAAYsQhUA0Ofs+s8D//xJKS0tvb9HAQYsQhVAv+r4GDK7PMPWwY7PtNnthyNhD198tcBOgvV3GKEK\noF91fAzZ/ztqn48h45k2YPDg1YK+RagC6Hc80wZgIOPvsL5DqAIA0Afa21pt+a/N2fFtLeg7hCoA\nAH3gs08atOM3H2nEu9f6e5Sg4W0tCDRCFQCAPmK3l4S/7Fnkwf5xZHZ81jyYCFUAABAQdnwWufF8\nuWInZvT3GIMWoQoAQWbH9y7a7fHamd2eRW7+uK6/RxjUghqqlmVp9erV+v3vf6+IiAg9+uijWrBg\nQTBHAIB+x7NOANA9QQ3V559/XtXV1dq5c6fq6+v1k5/8ROPGjdOMGTOCOQYA9DuedQKA/1tIsO7I\n4/Fo//79+qd/+iclJCQoJydHCxcu1K5du4I1AgAAAAaQoIXquXPn1NbWptTUVN+x9PR0VVVVBWsE\nAAAADCBBC9WmpiaNHDlSYWH/+26DqKgoXb9+XVeuXAnWGAAAABgggvYeVY/H0+Vfquj42rKsbv85\noaGft7XD4QjccAOEI+TzD1O2iz9/ekmSt7/HCCoesz3wmO2Bx2wPdnzMzR9fUGhopsLCuj7f2dFp\ngRK0UA0PD+8SpB1fDx06tNt/jsv1+bXbN/00cMMNGLP7ewAAAICgCdpL/7Gxsfrkk0/U3v6//yqF\n2+1WRESEXC5XsMYAAADAABG0UE1MTFRYWJhOnTrlO/bee+9pypQpwRoBAAAAA0jQQjUiIkKzZ8/W\nqlWrdObMGb3++usqLS3VI488EqwRAAAAMIA4vF5v0N4B3NLSop/+9Kf67W9/qxEjRmjhwoV6+OGH\ng3X3AAAAGECCGqoAAABAdwXtpX8AAADAH4QqAAAAjESoAgAAwEiEKgAAAIxEqAIAAMBIRoWqZVl6\n6qmnlJGRoaysLJWWln7ltdXV1br33nuVmpqq+fPn6/333w/ipAgEf/b9X//1X5ozZ47S0tI0e/Zs\nvfnmm0GcFIHgz7471NfXKy0tTeXl5UGYEIHkz74/+OADPfDAA0pJSdGsWbN0/PjxIE6KQPBn37//\n/e/1ve99T2lpaXrwwQdVXV0dxEkRSJZlKTc392v/ju5trxkVqs8//7yqq6u1c+dOrVq1SsXFxfrd\n737X5TqPx6OCggJlZGTo17/+tVJTU/UP//APamlp6Yep0VPd3fe5c+dUVFSk+fPn6+DBg7r33nu1\nZMkSffDBB/0wNXqqu/v+otWrV/N9PUB1d9/Xrl3TY489pltuuUWvvfaa7rnnHi1evFh/+tOf+mFq\n9FR3911TU6OlS5fqBz/4gQ4ePKiEhAQVFBTo+vXr/TA1esOyLD355JOqqan5ymsC0mteQ/z5z3/2\nTp061VteXu47tnXrVu/DDz/c5dp9+/Z5c3JyOh2bMWOG99/+7d/6fE4Ehj/7Xr9+vffxxx/vdOzR\nRx/1vvjii30+JwLDn313OHDggPf+++/3JiQkeE+cOBGMMREg/uz7X//1X70zZszodGzevHnet956\nq8/nRGD4s+/S0lJvXl6e7+tr1655b7vtNu8f/vCHoMyKwKipqfHOnj3bO3v27K/9OzoQvWbMM6rn\nzp1TW1ubUlNTfcfS09NVVVXV5dqqqiqlp6d3Onb77bersrKyz+dEYPiz77lz5+pHP/pRl+PXrl3r\n0xkROP7sW5KuXLmiDRs26Omnn5aXf5NkwPFn3+Xl5crOzu50bN++fbrzzjv7fE4Ehj/7HjlypGpq\nalRRUSGv16uysjKNGDFCcXFxwRwZvXTixAlNnz5de/fu/dq/owPRa2E9njLAmpqaNHLkSIWF/e9I\nUVFRun79uq5cuaLIyEjf8cuXL+vWW2/tdPuoqKivffoZZvFn3xMnTux02//5n//Ru+++qwceeCBo\n86J3/Nm3JD333HOaO3eu4uPjgz0qAsCffdfV1Sk5OVn//M//rDfffFPf/OY39eMf/1i33357f4yO\nHvBn39/73vf05ptv6oEHHlBoaKhCQkK0bds2jRgxoj9GRw/df//93bouEL1mzDOqHo9HTqez07GO\nry3L6nS8paXlS6/9y+tgLn/2/UV/+tOfVFRUpPT0dP3N3/xNn86IwPFn30ePHlVlZaUWLVoUtPkQ\nWP7s+89//rNeeuklxcTE6KWXXtK0adP02GOPqbGxMWjzonf82fcnn3wit9utVatWad++fZozZ46W\nL1/Oe5IHqUD0mjGhGh4e3mXwjq+HDh3arWsjIiL6dkgEjD/77uB2u/XII4/I4XBo48aNfT4jAqe7\n+75+/bpWrVqlVatWdfnLDQOHP9/foaGhSkxM1OLFi5WQkKClS5dqwoQJOnDgQNDmRe/4s+/169fr\ntttu0/3336+kpCStWbNGQ4cO1a9//eugzYvgCUSvGROqsbGx+uSTT9Te3u475na7FRERIZfL1eXa\npqamTsfcbreio6ODMit6z599S1JjY6MefPBBtbW1aefOnV1eKobZurvvqqoq1dfXq6ioSGlpaUpL\nS5MkPf7441q9enWwx0YP+fP9HR0d3eXtPRMmTFBDQ0NQZkXv+bPv999/XwkJCb6vHQ6HEhIS9NFH\nHwVtXgRPIHrNmFBNTExUWFiYTp065Tv23nvvacqUKV2uTUlJ6fJG3IqKik5v5IbZ/Nm3x+PRwoUL\nNWTIEO3atUujR48O5qgIgO7uOyUlRb/73e904MABHTx4UAcPHpQkPfvss1qyZElQZ0bP+fP9nZqa\nqnPnznU6dv78eY0bN67P50Rg+LPvmJiYLu9P/OMf/6hvfvObfT4ngi8QvWZMqEZERGj27NlatWqV\nzpw5o9dff12lpaV65JFHJH1e4B2fs/ad73xHzc3N+tnPfqba2lo988wz8ng8+u53v9ufDwF+8Gff\n//Iv/6L6+nqtXbtW7e3tcrvdcrvd/NT/ANLdfTudTo0fP77TL+nz/7mNGjWqPx8C/ODP9/d9992n\nDz74QMXFxbpw4YI2btyo+vp6zZo1qz8fAvzgz77nz5+vffv26cCBA7pw4YLWr1+vhoYGzZkzpz8f\nAgIo4L3W80/RCjyPx+Ndvny5Ny0tzXvnnXd6f/WrX/nO3XbbbZ0+d6uqqso7d+5cb0pKivfee+/1\nnj17tj9GRi90d99/+7d/601ISOjya/ny5f01OnrAn+/vL+JzVAcmf/ZdUVHhnTt3rnfq1KneuXPn\net97773+GBm94M++9+/f7/3ud7/rvf32270PPvgg//8e4P7y7+hA95rD6+VDCgEAAGAeY176BwAA\nAL6IUAUAAICRCFUAAAAYiVAFAACAkQhVAAAAGIlQBQAAgJEIVQAAABiJUAUAAICRCFUAAAAYiVAF\nAACAkQhVAAAAGOn/A4ESAwl1zlWzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10e584c90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">25% 91\n",
      ">50% 50\n",
      ">75% 25\n",
      ">90% 10\n",
      ">95% 0\n"
     ]
    }
   ],
   "source": [
    "plt.hist(table_tongue.pa(inplace=False).sum(axis=\"observation\")/len(table_tongue.ids()))\n",
    "plt.show()\n",
    "print \">25%\", np.sum(table_tongue.pa(inplace=False).sum(axis=\"observation\")/len(table_tongue.ids())>.25)\n",
    "print \">50%\", np.sum(table_tongue.pa(inplace=False).sum(axis=\"observation\")/len(table_tongue.ids())>.50)\n",
    "print \">75%\", np.sum(table_tongue.pa(inplace=False).sum(axis=\"observation\")/len(table_tongue.ids())>.75)\n",
    "print \">90%\", np.sum(table_tongue.pa(inplace=False).sum(axis=\"observation\")/len(table_tongue.ids())>.90)\n",
    "print \">95%\", np.sum(table_tongue.pa(inplace=False).sum(axis=\"observation\")/len(table_tongue.ids())>.99)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqoAAAHcCAYAAAAeFogrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3X90VPWd//FXJnESpKTG/IJGopK6kgRMYkj2sIo9B3+t\nLgTdSTz+rMEfWSUJ3WNtN7hboViLFqql/CpEqhWPZ4WMLlL30Jba+iOgYhJMNGIdtJpgCBPkR9RJ\nLsnM9w+bWWOqX26YTD6ZeT7O4eTM/dw7877nxcCLzM0lJhAIBAQAAAAYxjHaAwAAAAB/D0UVAAAA\nRqKoAgAAwEgUVQAAABiJogoAAAAjUVQBAABgJIoqAAAAjERRBQAAgJEoqgAAADCS7aL68ccfa+HC\nhSoqKtLll1+uZ555JrjW3t6u+fPnq6CgQHPmzFF9ff2gY3fu3Km5c+cqPz9f5eXlamtrO/kzAAAA\nQESyXVQXLFiggwcPatOmTbrnnnv0wAMPaMeOHcG1tLQ0ud1ulZSUqKqqSgcOHJAkdXR0qLKyUi6X\nS263W0lJSaqsrAzt2QAAACBixAQCgcCJ7vzmm2+qrKxMO3bsUEZGhiSptrZWf/zjH/Xv//7vWrBg\ngXbt2qX4+HhJ0vz581VYWKiqqiqtXLlSDQ0NevzxxyVJPT09uuCCC/SrX/1KRUVFI3BqAAAAGMts\nfUe1ra1Np59+erCkStK5556rN998U6+//rpyc3ODJVWSCgsLtWfPHklSc3PzoEKakJCgnJwcNTU1\nnew5AAAAIALZKqopKSk6duyYent7g9s6OjrU19enQ4cOKS0tbdD+ycnJ6uzslCQdPHhwyHpKSkpw\nHQAAAPgiW0U1Ly9PqampWrp0qXw+nz744AM99thjiomJUW9vr5xO56D9nU6nLMuS9PlH/V+3DgAA\nAHyRraLqdDr1y1/+Uq+++qoKCwt144036tprr1UgEJDD4RhSOi3LUkJCgiQpPj7+a9dPhI3LaQEA\nADDGxdk9YNq0adqxY4cOHTqkpKQkvfTSSzr99NOVmZmpl19+edC+XV1dSk1NlSSlp6fL6/UOWc/O\nzj7h146JidGxYz719/vtjo0xJjbWocTEceQdJcg7upB3dCHv6DKQd6jYKqpHjx7VnXfeqXXr1ik5\nOVmS9Oc//1nFxcU677zztH79elmWFfyIv6GhQTNmzJD0+WUDjY2Nwefy+XxqbW1VdXW1rYH7+/3q\n6+M3erQg7+hC3tGFvKMLeWM4bH30/81vflM+n0/Lly9XW1ubtmzZomeeeUa33367iouL9a1vfUs1\nNTXyeDzasGGDWlpaVFpaKklyuVxqbGxUbW2tPB6PFi1apMzMTBUXF4/IiQEAAGBss33D/4cfflgf\nfPCBSkpK9Pjjj2vlypXKzc2Vw+HQ2rVr5fV65XK5tG3bNq1Zs0YTJ06UJGVkZGjVqlVyu90qKytT\nd3e3Vq9eHfITAgAAQGSwdcN/Exw+/CkfHUSBuDiHkpLGk3eUIO/oQt7Rhbyjy0DeoWL7O6oAAABA\nOFBUAQAAYCSKKgAAAIxEUQUAAICRKKoAAAAwEkUVAAAARqKoAgAAwEgUVQAAABiJogoAAAAjUVQB\nAABgJIoqAAAAjERRBQAAgJEoqgAAADASRRUAAABGoqgCAADASBRVAAAAGImiCgAAACNRVAEAAGAk\niioAAACMRFEFAACAkSiqAAAAMBJFFQAAAEaKG+0BhuvJzW693LB3tMcIu9ysSaqsuGW0xwAAABhx\nY7aofvBhm3qSZo72GGH3UWfjaI8AAAAQFnz0DwAAACNRVAEAAGAkiioAAACMRFEFAACAkSiqAAAA\nMBJFFQAAAEaiqAIAAMBIFFUAAAAYiaIKAAAAI1FUAQAAYCSKKgAAAIxku6geOHBAd9xxhwoLC3Xx\nxRfrN7/5TXCtvb1d8+fPV0FBgebMmaP6+vpBx+7cuVNz585Vfn6+ysvL1dbWdvJnAAAAgIhku6h+\n73vf0/jx4/XMM8/onnvu0S9+8Qvt2LFDkrRgwQKlpaXJ7XarpKREVVVVOnDggCSpo6NDlZWVcrlc\ncrvdSkpKUmVlZWjPBgAAABHDVlE9duyY3njjDd15553KzMzUxRdfrFmzZumVV17RK6+8ovb2di1d\nulRTpkxRRUWF8vPzVVdXJ0navHmzpk+frvLycmVlZWnZsmXav3+/du/ePSInBgAAgLHNVlFNSEjQ\nuHHj5Ha71dfXp/fee0+NjY3Kzs7WG2+8odzcXMXHxwf3Lyws1J49eyRJzc3NKioqGvRcOTk5ampq\nCtGpAAAAIJLYKqpOp1P33nuv/vu//1t5eXm68sorddFFF8nlcsnr9SotLW3Q/snJyers7JQkHTx4\ncMh6SkpKcB0AAAD4oji7B+zbt0+zZ8/Wrbfeqr/85S+67777NHPmTPl8PjmdzkH7Op1OWZYlSerp\n6fna9RMVG/t5t3Y4YuyOHhEcjhjFxUX+zRoGch74ishG3tGFvKMLeUeXUOdsq6ju2rVLdXV1evHF\nF+V0OpWTk6MDBw5o3bp1mjlzpo4cOTJof8uylJCQIEmKj48fUkoty1JiYqKtgRMTx/3t+U6R7HXc\niOB0xikpafxojxE2A3kjOpB3dCHv6ELeGA5bRfWtt97SWWedNeg7o9nZ2Vq/fr3S09P17rvvDtq/\nq6tLqampkqT09HR5vd4h69nZ2bYGPnbMp/5+v3p7j9s6LlJYVp8OH/50tMcYcbGxDiUmjgvmjchG\n3tGFvKMLeUeXgbxDxVZRTUtL0wcffKC+vj7FxX1+6HvvvaczzjhDeXl5Wr9+vSzLChbZhoYGzZgx\nQ5KUl5enxsbG4HP5fD61traqurra1sD9/X719fnl9wdsHRcp/P6A+vqi540+kDeiA3lHF/KOLuSN\n4bB1IcHs2bMVFxen//qv/9Jf//pXPf/881q/fr2++93vqqioSJMmTVJNTY08Ho82bNiglpYWlZaW\nSpJcLpcaGxtVW1srj8ejRYsWKTMzU8XFxSNyYgAAABjbbBXVb3zjG3rsscfk9XpVVlamBx98UJWV\nlSorK5PD4dC6devk9Xrlcrm0bds2rVmzRhMnTpQkZWRkaNWqVXK73SorK1N3d7dWr149IicFAACA\nsc/2T/1nZWVp48aNf3dt8uTJ2rRp01ceO2vWLG3fvt3uSwIAACAKca8IAAAAGImiCgAAACNRVAEA\nAGAkiioAAACMRFEFAACAkSiqAAAAMBJFFQAAAEaiqAIAAMBIFFUAAAAYiaIKAAAAI1FUAQAAYCSK\nKgAAAIxEUQUAAICRKKoAAAAwEkUVAAAARqKoAgAAwEgUVQAAABiJogoAAAAjUVQBAABgJIoqAAAA\njERRBQAAgJEoqgAAADASRRUAAABGoqgCAADASBRVAAAAGImiCgAAACNRVAEAAGAkiioAAACMRFEF\nAACAkSiqAAAAMBJFFQAAAEaiqAIAAMBIFFUAAAAYiaIKAAAAI1FUAQAAYCSKKgAAAIxkq6g+88wz\nmjp1qrKzswd9zcnJkSS1tbVp/vz5Kigo0Jw5c1RfXz/o+J07d2ru3LnKz89XeXm52traQncmAAAA\niCi2iuq//Mu/qL6+Xi+//LLq6+v1pz/9SWeeeaZuvvlmSVJlZaXS0tLkdrtVUlKiqqoqHThwQJLU\n0dGhyspKuVwuud1uJSUlqbKyMvRnBAAAgIhgq6g6nU4lJycHf23dulWSdNddd2nXrl1qb2/X0qVL\nNWXKFFVUVCg/P191dXWSpM2bN2v69OkqLy9XVlaWli1bpv3792v37t2hPysAAACMecO+RvXo0aN6\n5JFHdPfdd+uUU05Rc3OzcnNzFR8fH9ynsLBQe/bskSQ1NzerqKgouJaQkKCcnBw1NTWdxPgAAACI\nVMMuqk8++aTS09N16aWXSpK8Xq/S0tIG7ZOcnKzOzk5J0sGDB4esp6SkBNcBAACAL4ob7oF1dXWq\nqKgIPvb5fHI6nYP2cTqdsixLktTT0/O16ycqNvbzbu1wxAxn7DHP4YhRXFzk36xhIOeBr4hs5B1d\nyDu6kHd0CXXOwyqqzc3N6uzs1JVXXhncFh8fr6NHjw7az7IsJSQkBNe/XEoty1JiYqKt105MHPe3\n5ztFstdxI4LTGaekpPGjPUbYDOSN6EDe0YW8owt5YziGVVRffvllFRUVacKECcFt6enp8ng8g/br\n6upSampqcN3r9Q5Zz87OtvXax4751N/vV2/v8eGMPuZZVp8OH/50tMcYcbGxDiUmjgvmjchG3tGF\nvKMLeUeXgbxDZdjfUT3//PMHbcvLy1Ntba0sywp+xN/Q0KAZM2YE1xsbG4P7+3w+tba2qrq62tZr\n9/f71dfnl98fGM7oY57fH1BfX/S80QfyRnQg7+hC3tGFvDEcw7qQ4C9/+YuysrIGbSsuLtakSZNU\nU1Mjj8ejDRs2qKWlRaWlpZIkl8ulxsZG1dbWyuPxaNGiRcrMzFRxcfHJnwUAAAAizrCK6scff6xv\nfvObg5/I4dDatWvl9Xrlcrm0bds2rVmzRhMnTpQkZWRkaNWqVXK73SorK1N3d7dWr1598mcAAACA\niDSsj/4H7o36ZZMnT9amTZu+8rhZs2Zp+/btw3lJAAAARBnuFQEAAAAjUVQBAABgJIoqAAAAjERR\nBQAAgJEoqgAAADASRRUAAABGoqgCAADASBRVAAAAGImiCgAAACNRVAEAAGAkiioAAACMRFEFAACA\nkSiqAAAAMBJFFQAAAEaiqAIAAMBIFFUAAAAYiaIKAAAAI1FUAQAAYCSKKgAAAIxEUQUAAICRKKoA\nAAAwEkUVAAAARqKoAgAAwEgUVQAAABiJogoAAAAjUVQBAABgJIoqAAAAjERRBQAAgJEoqgAAADAS\nRRUAAABGoqgCAADASBRVAAAAGImiCgAAACNRVAEAAGAkiioAAACMZLuoWpalH//4xyouLtaFF16o\nhx9+OLjW3t6u+fPnq6CgQHPmzFF9ff2gY3fu3Km5c+cqPz9f5eXlamtrO/kzAAAAQESyXVR/8pOf\naNeuXfr1r3+tFStWaPPmzdq8ebMkacGCBUpLS5Pb7VZJSYmqqqp04MABSVJHR4cqKyvlcrnkdruV\nlJSkysrK0J4NAAAAIkacnZ2PHj2qp59+Wo899pimTZsmSbrlllv0xhtvKDMzU+3t7dqyZYvi4+NV\nUVGhXbt2qa6uTlVVVdq8ebOmT5+u8vJySdKyZct0wQUXaPfu3SoqKgr5iQEAAGBss1VUGxoaNGHC\nBM2YMSO47fbbb5ckrV+/Xrm5uYqPjw+uFRYWas+ePZKk5ubmQYU0ISFBOTk5ampqoqgCAABgCFsf\n/be1tSkjI0P/8z//oyuuuEKXXHKJ1q5dq0AgIK/Xq7S0tEH7Jycnq7OzU5J08ODBIespKSnBdQAA\nAOCLbH1H9bPPPtNf//pXbd68WQ888IC8Xq/uvfdejRs3Tj6fT06nc9D+TqdTlmVJknp6er52/UTF\nxn7erR2OGFvHRQqHI0ZxcZF/s4aBnAe+IrKRd3Qh7+hC3tEl1DnbKqqxsbH69NNP9dBDD2nixImS\npP379+vJJ5/UhRdeqCNHjgza37IsJSQkSJLi4+OHlFLLspSYmGhr4MTEcX97vlMkex03IjidcUpK\nGj/aY4TNQN6IDuQdXcg7upA3hsNWUU1LS1N8fHywpErS2Wefrc7OTqWnp+vdd98dtH9XV5dSU1Ml\nSenp6fJ6vUPWs7OzbQ187JhP/f1+9fYet3VcpLCsPh0+/OlojzHiYmMdSkwcF8wbkY28owt5Rxfy\nji4DeYeKraKal5en3t5effDBBzrzzDMlSfv27VNGRoby8vK0fv16WZYV/Ii/oaEh+INXeXl5amxs\nDD6Xz+dTa2urqqurbQ3c3+9XX59ffn/A1nGRwu8PqK8vet7oA3kjOpB3dCHv6ELeGA5bFxKcffbZ\n+s53vqOamhrt3btXL730kmpra3X99derqKhIkyZNUk1NjTwejzZs2KCWlhaVlpZKklwulxobG1Vb\nWyuPx6NFixYpMzNTxcXFI3JiAAAAGNtsX/G6YsUKnXnmmbrhhhu0aNEi3XTTTbrhhhvkcDi0bt06\neb1euVwubdu2TWvWrAleJpCRkaFVq1bJ7XarrKxM3d3dWr16dchPCAAAAJEhJhAIjKnP0A8f/lR9\nfX4tW/ELvdt33miPE3aTjjfq/h/dPdpjjLi4OIeSksYH80ZkI+/oQt7Rhbyjy0DeocK9IgAAAGAk\niioAAACMRFEFAACAkSiqAAAAMBJFFQAAAEaiqAIAAMBIFFUAAAAYiaIKAAAAI1FUAQAAYCSKKgAA\nAIxEUQUAAICRKKoAAAAwEkUVAAAARqKoAgAAwEgUVQAAABiJogoAAAAjUVQBAABgJIoqAAAAjERR\nBQAAgJEoqgAAADASRRUAAABGoqgCAADASBRVAAAAGImiCgAAACNRVAEAAGAkiioAAACMRFEFAACA\nkSiqAAAAMBJFFQAAAEaiqAIAAMBIFFUAAAAYiaIKAAAAI1FUAQAAYCSKKgAAAIxEUQUAAICRbBfV\nHTt2aOrUqcrOzg5+/d73vidJam9v1/z581VQUKA5c+aovr5+0LE7d+7U3LlzlZ+fr/LycrW1tYXm\nLAAAABBxbBdVj8ej2bNnq76+XvX19Xr55Zd1//33S5IWLFigtLQ0ud1ulZSUqKqqSgcOHJAkdXR0\nqLKyUi6XS263W0lJSaqsrAzt2QAAACBi2C6q+/bt0znnnKPTTz9dycnJSk5O1je+8Q3t2rVL7e3t\nWrp0qaZMmaKKigrl5+errq5OkrR582ZNnz5d5eXlysrK0rJly7R//37t3r075CcFAACAsW9YRfXs\ns88esr25uVm5ubmKj48PbissLNSePXuC60VFRcG1hIQE5eTkqKmpaThzAwAAIMLZLqrvv/++Xnrp\nJV1++eW69NJL9fOf/1zHjx+X1+tVWlraoH2Tk5PV2dkpSTp48OCQ9ZSUlOA6AAAA8EVxdnb+6KOP\n1NPTo/j4eK1cuVLt7e26//771dPTI5/PJ6fTOWh/p9Mpy7IkST09PV+7fqJiYz/v1g5HjK3jIoXD\nEaO4uMi/WcNAzgNfEdnIO7qQd3Qh7+gS6pxtFdVvfetbevXVV5WYmChJmjp1qvx+v37wgx/oX//1\nX3Xs2LFB+1uWpYSEBElSfHz8kFJqWVbwuU5UYuK4vz3fKZK9jhsRnM44JSWNH+0xwmYgb0QH8o4u\n5B1dyBvDYauoShpSLLOystTb26uUlBTt27dv0FpXV5dSU1MlSenp6fJ6vUPWs7Ozbb3+sWM+9ff7\n1dt73O7oEcGy+nT48KejPcaIi411KDFxXDBvRDbyji7kHV3IO7oM5B0qtorqyy+/rO9///t68cUX\ngz801draqqSkJM2YMUO//vWvZVlW8CP+hoYGzZgxQ5KUl5enxsbG4HP5fD61traqurra1sD9/X71\n9fnl9wdsHRcp/P6A+vqi540+kDeiA3lHF/KOLuSN4bB1IUFBQYHGjRun//zP/9T777+vF154QcuX\nL9ftt9+uoqIiTZo0STU1NfJ4PNqwYYNaWlpUWloqSXK5XGpsbFRtba08Ho8WLVqkzMxMFRcXj8iJ\nAQAAYGyzVVTHjx+vjRs36vDhwyotLdWPfvQjXXvttbrlllvkcDi0bt06eb1euVwubdu2TWvWrNHE\niRMlSRkZGVq1apXcbrfKysrU3d2t1atXj8hJAQAAYOyzfY1qVlaWNm7c+HfXJk+erE2bNn3lsbNm\nzdL27dvtviQAAACiEPeKAAAAgJEoqgAAADASRRUAAABGoqgCAADASBRVAAAAGImiCgAAACNRVAEA\nAGAkiioAAACMRFEFAACAkSiqAAAAMBJFFQAAAEaiqAIAAMBIFFUAAAAYiaIKAAAAI1FUAQAAYCSK\nKgAAAIxEUQUAAICRKKoAAAAwEkUVAAAARqKoAgAAwEgUVQAAABiJogoAAAAjUVQBAABgJIoqAAAA\njERRBQAAgJEoqgAAADASRRUAAABGoqgCAADASBRVAAAAGImiCgAAACNRVAEAAGAkiioAAACMRFEF\nAACAkSiqAAAAMBJFFQAAAEaiqAIAAMBIwy6qFRUVWrRoUfBxe3u75s+fr4KCAs2ZM0f19fWD9t+5\nc6fmzp2r/Px8lZeXq62tbfhTAwAAIOINq6g+99xzevHFFwdtq6ysVFpamtxut0pKSlRVVaUDBw5I\nkjo6OlRZWSmXyyW3262kpCRVVlae/PQAAACIWLaL6tGjR7V8+XKdd955wW27du1SW1ubli5dqilT\npqiiokL5+fmqq6uTJG3evFnTp09XeXm5srKytGzZMu3fv1+7d+8O3ZkAAAAgotguqg8++KDmzZun\nrKys4Lbm5mbl5uYqPj4+uK2wsFB79uwJrhcVFQXXEhISlJOTo6amppOZHQAAABHMVlHdtWuXGhoa\nhnxs7/V6lZaWNmhbcnKyOjs7JUkHDx4csp6SkhJcBwAAAL4s7kR3tCxLS5Ys0eLFi+V0Oget+Xy+\nIducTqcsy5Ik9fT0fO26HbGxn3drhyPG9rGRwOGIUVxc5N+sYSDnga+IbOQdXcg7upB3dAl1zidc\nVFetWqVp06bpn/7pn4asxcfH6+jRo4O2WZalhISE4PqXS6llWUpMTLQ9cGLiuL895ymS/Z475jmd\ncUpKGj/aY4TNQN6IDuQdXcg7upA3huOEi+r//u//6tChQyooKJAkHT9+XJL0u9/9TnfccYc8Hs+g\n/bu6upSamipJSk9Pl9frHbKenZ1te+Bjx3zq7/ert/e47WMjgWX16fDhT0d7jBEXG+tQYuK4YN6I\nbOQdXcg7upB3dBnIO1ROuKg+8cQT6uvrCz5evny5JOkHP/iB9u/frw0bNsiyrOBH/A0NDZoxY4Yk\nKS8vT42NjcFjfT6fWltbVV1dbXvg/n6/+vr88vsDto+NBH5/QH190fNGH8gb0YG8owt5RxfyxnCc\n8IUEkyZN0uTJk4O/xo8fr/Hjx2vy5MkqLi7WpEmTVFNTI4/How0bNqilpUWlpaWSJJfLpcbGRtXW\n1srj8WjRokXKzMxUcXHxiJ0YAAAAxraQXPHqcDi0du1aeb1euVwubdu2TWvWrNHEiRMlSRkZGVq1\napXcbrfKysrU3d2t1atXh+KlAQAAEKFO+KP/L1u2bNmgx5MnT9amTZu+cv9Zs2Zp+/btw305AAAA\nRBnuFQEAAAAjUVQBAABgJIoqAAAAjERRBQAAgJEoqgAAADASRRUAAABGoqgCAADASBRVAAAAGImi\nCgAAACNRVAEAAGAkiioAAACMRFEFAACAkSiqAAAAMBJFFQAAAEaiqAIAAMBIFFUAAAAYiaIKAAAA\nI1FUAQAAYCSKKgAAAIxEUQUAAICRKKoAAAAwEkUVAAAARqKoAgAAwEgUVQAAABiJogoAAAAjUVQB\nAABgJIoqAAAAjERRBQAAgJEoqgAAADASRRUAAABGoqgCAADASBRVAAAAGImiCgAAACNRVAEAAGAk\niioAAACMZLuofvjhh7r11ltVUFCg2bNna+PGjcG19vZ2zZ8/XwUFBZozZ47q6+sHHbtz507NnTtX\n+fn5Ki8vV1tb28mfAQAAACKSraIaCARUUVGhlJQUbd26VUuWLNG6dev03HPPSZIWLFigtLQ0ud1u\nlZSUqKqqSgcOHJAkdXR0qLKyUi6XS263W0lJSaqsrAz9GQEAACAi2CqqXV1dysnJ0eLFi5WZmamL\nLrpIM2fOVENDg1555RW1t7dr6dKlmjJliioqKpSfn6+6ujpJ0ubNmzV9+nSVl5crKytLy5Yt0/79\n+7V79+4ROTEAAACMbbaKampqqh566CGdeuqpkqSGhga9/vrrKi4u1htvvKHc3FzFx8cH9y8sLNSe\nPXskSc3NzSoqKgquJSQkKCcnR01NTaE4DwAAAESYYf8w1ezZs3XjjTcqPz9fl112mbxer9LS0gbt\nk5ycrM7OTknSwYMHh6ynpKQE1wEAAIAvihvugatWrVJXV5eWLFmin/70p/L5fHI6nYP2cTqdsixL\nktTT0/O16ycqNvbzbu1wxAx39DHN4YhRXFzk36xhIOeBr4hs5B1dyDu6kHd0CXXOwy6qubm5kqSa\nmhrdfffdKi0t1bFjxwbtY1mWEhISJEnx8fFDSqllWUpMTLT1uomJ4/72fKdI9jpuRHA645SUNH60\nxwibgbwRHcg7upB3dCFvDIetonro0CE1NTXpkksuCW779re/rePHjys1NVX79u0btH9XV5dSU1Ml\nSenp6fJ6vUPWs7OzbQ187JhP/f1+9fYet3VcpLCsPh0+/OlojzHiYmMdSkwcF8wbkY28owt5Rxfy\nji4DeYeKraLa3t6u6upqvfDCC8HrTVtaWpScnKzCwkJt3LhRlmUFP+JvaGjQjBkzJEl5eXlqbGwM\nPpfP51Nra6uqq6ttDdzf71dfn19+f8DWcZHC7w+ory963ugDeSM6kHd0Ie/oQt4YDlsXEkyfPl3T\npk3TPffco3379umFF17QihUrdOedd6qoqEiTJk1STU2NPB6PNmzYoJaWFpWWlkqSXC6XGhsbVVtb\nK4/Ho0WLFikzM1PFxcUjcmIAAAAY22wVVYfDobVr1+rUU0/Vtddeqx/96Ef67ne/qxtvvFEOh0Pr\n1q2T1+uVy+XStm3btGbNGk2cOFGSlJGRoVWrVsntdqusrEzd3d1avXr1iJwUAAAAxj7bP0yVmpqq\nX/7yl393bfLkydq0adNXHjtr1ixt377d7ksCAAAgCnGvCAAAABiJogoAAAAjUVQBAABgJIoqAAAA\njERRBQAAgJEoqgAAADASRRUAAABGoqgCAADASBRVAAAAGImiCgAAACNRVAEAAGAkiioAAACMRFEF\nAACAkSiqAAAAMBJFFQAAAEaiqAIAAMBIFFUAAAAYiaIKAAAAI1FUAQAAYCSKKgAAAIxEUQUAAICR\nKKoAAAAwEkUVAAAARqKoAgAAwEgUVQAAABiJogoAAAAjUVQBAABgJIoqAAAAjERRBQAAgJEoqgAA\nADASRRVmqRFlAAATFUlEQVQAAABGoqgCAADASBRVAAAAGImiCgAAACNRVAEAAGAkiioAAACMZKuo\ndnZ2auHChfrHf/xHfec739EDDzwgy7IkSe3t7Zo/f74KCgo0Z84c1dfXDzp2586dmjt3rvLz81Ve\nXq62trbQnQUAAAAijq2iunDhQvX29urJJ5/UQw89pD/96U9auXKlJGnBggVKS0uT2+1WSUmJqqqq\ndODAAUlSR0eHKisr5XK55Ha7lZSUpMrKytCfDQAAACLGCRfV9957T83NzVq2bJmysrJUWFiohQsX\n6re//a1eeeUVtbe3a+nSpZoyZYoqKiqUn5+vuro6SdLmzZs1ffp0lZeXKysrS8uWLdP+/fu1e/fu\nETsxAAAAjG0nXFRTU1P1yCOP6PTTTx+0vbu7W2+88YZyc3MVHx8f3F5YWKg9e/ZIkpqbm1VUVBRc\nS0hIUE5Ojpqamk52fgAAAESoEy6qEyZM0AUXXBB8HAgE9MQTT2jmzJnyer1KS0sbtH9ycrI6Ozsl\nSQcPHhyynpKSElwHAAAAvixuuAf+7Gc/09tvv626ujo9+uijcjqdg9adTmfwB616enq+dt2O2NjP\nu7XDETPMycc2hyNGcXGRf7OGgZwHviKykXd0Ie/oQt7RJdQ5D6uoLl++XJs2bdIvfvELffvb31Z8\nfLyOHj06aB/LspSQkCBJio+PH1JKLctSYmKi7ddOTBz3t+c8RbLfc8c8pzNOSUnjR3uMsBnIG9GB\nvKMLeUcX8sZw2C6q9913n5566iktX75cl1xyiSQpPT1dHo9n0H5dXV1KTU0Nrnu93iHr2dnZtgc+\ndsyn/n6/enuP2z42ElhWnw4f/nS0xxhxsbEOJSaOC+aNyEbe0YW8owt5R5eBvEPFVlFdvXq1nnrq\nKT388MO69NJLg9vz8vJUW1sry7KCH/E3NDRoxowZwfXGxsbg/j6fT62traqurrY9cH+/X319fvn9\nAdvHRgK/P6C+vuh5ow/kjehA3tGFvKMLeWM4TvhCgn379mndunWqqKhQQUGBurq6gr+Ki4s1adIk\n1dTUyOPxaMOGDWppaVFpaakkyeVyqbGxUbW1tfJ4PFq0aJEyMzNVXFw8YicGAACAse2Ei+of//hH\n+f1+rVu3TrNmzdKsWbN04YUXatasWXI4HFqzZo28Xq9cLpe2bdumNWvWaOLEiZKkjIwMrVq1Sm63\nW2VlZeru7tbq1atH7KQAAAAw9p3wR/8VFRWqqKj4yvXMzExt2rTpK9dnzZql7du325sOAAAAUYt7\nRQAAAMBIFFUAAAAYiaIKAAAAI1FUAQAAYCSKKgAAAIxEUQUAAICRKKoAAAAwEkUVAAAARqKoAgAA\nwEgUVQAAABiJogoAAAAjUVQBAABgJIoqAAAAjERRBQAAgJEoqgAAADASRRUAAABGoqgCAADASBRV\nAAAAGImiCgAAACNRVAEAAGAkiioAAACMRFEFAACAkSiqAAAAMBJFFQAAAEaiqAIAAMBIFFUAAAAY\niaIKAAAAI1FUAQAAYCSKKgAAAIxEUQUAAICRKKoAAAAwEkUVAAAARqKoAgAAwEgUVQAAABiJogoA\nAAAjDbuoWpaluXPnavfu3cFt7e3tmj9/vgoKCjRnzhzV19cPOmbnzp2aO3eu8vPzVV5erra2tuFP\nDgAAgIg2rKJqWZbuuusueTyeQdsrKyuVlpYmt9utkpISVVVV6cCBA5Kkjo4OVVZWyuVyye12Kykp\nSZWVlSd/BgAAAIhItovqvn37dM0116i9vX3Q9l27dqmtrU1Lly7VlClTVFFRofz8fNXV1UmSNm/e\nrOnTp6u8vFxZWVlatmyZ9u/fP+g7sgAAAMAA20X1tdde08yZM/XUU08pEAgEtzc3Nys3N1fx8fHB\nbYWFhdqzZ09wvaioKLiWkJCgnJwcNTU1ncz8AAAAiFBxdg+47rrr/u52r9ertLS0QduSk5PV2dkp\nSTp48OCQ9ZSUlOA6AAAA8EW2i+pX8fl8cjqdg7Y5nU5ZliVJ6unp+dr1ExUb+/k3gR2OmJOYduxy\nOGIUFxf5N2sYyHngKyIbeUcX8o4u5B1dQp1zyIpqfHy8jh49OmibZVlKSEgIrn+5lFqWpcTERFuv\nk5g47m/Pd4pkr+NGBKczTklJ40d7jLAZyBvRgbyjC3lHF/LGcISsqKanpw+5C0BXV5dSU1OD616v\nd8h6dna2rdc5dsyn/n6/enuPn9zAY5Rl9enw4U9He4wRFxvrUGLiuGDeiGzkHV3IO7qQd3QZyDtU\nQlZU8/LyVFtbK8uygh/xNzQ0aMaMGcH1xsbG4P4+n0+tra2qrq629Tr9/X719fnl9wf+/ztHIL8/\noL6+6HmjD+SN6EDe0YW8owt5YzhCdiFBcXGxJk2apJqaGnk8Hm3YsEEtLS0qLS2VJLlcLjU2Nqq2\ntlYej0eLFi1SZmamiouLQzUCAAAAIshJFdWYmP/7gSaHw6G1a9fK6/XK5XJp27ZtWrNmjSZOnChJ\nysjI0KpVq+R2u1VWVqbu7m6tXr365KYHAABAxDqpj/7ffvvtQY8nT56sTZs2feX+s2bN0vbt20/m\nJaOav/+4vJ371dTUMNqjjLgvX9OUmzt9yF0jAABAZAvZNaoYece8f1X30VN1329eH+1Rwqr70If6\n2V1SQUHhaI8CAADCiKI6xkxIztRpE88Z7TEAAABGHHffBQAAgJEoqgAAADASRRUAAABGoqgCAADA\nSBRVAAAAGImiCgAAACNRVAEAAGAkiioAAACMRFEFAACAkSiqAAAAMBJFFQAAAEaiqAIAAMBIFFUA\nAAAYiaIKAAAAI8WN9gDA/4+/v0/vvLN3tMcYFbm50+V0Okd7DAAARgVFFcb79EiHNj73kSa88slo\njxJW3Yc+1M/ukgoKCkd7FAAARgVFFWPChORMnTbxnNEeAwAAhBHXqAIAAMBIFFUAAAAYiaIKAAAA\nI1FUAQAAYCSKKgAAAIxEUQUAAICRKKoAAAAwEkUVAAAARqKoAgAAwEj8z1SAofz9fXrnnb2jPUZY\nxMY6lJg4TseO+TR1aq6cTudojwQAMABFFTDUp0c6tPG5jzThlU9Ge5Sw6T70oX52l18FBYWjPQoA\nwAAUVcBgE5IzddrEc0Z7DAAARgXXqAIAAMBIFFUAAAAYiY/+AWCUWZalt95qGe0xwmbgh+cyM78t\nh4O/hgB8Nf6EAIBR9tZbLfrhQ09rQnLmaI8SNt2HPtTPf1Cq884rGO1RABgsrEXVsiwtWbJEf/jD\nH5SQkKBbbrlF8+fPD+cIAGAkfnAOAIYKa1F98MEH1draqk2bNqm9vV3/8R//oYyMDF122WXhHAOA\noaLp3rFfFI3nHK2i7TIPiUs9cHLC9jvG5/Oprq5OGzdu1NSpUzV16lTddttteuKJJyiqACRF571j\nJanzvd1Kn1I02mOElb+/T3v3vq3+fv9ojxJW77yzVxufa42qyzwkLvXA8IWtqO7du1f9/f3Kz88P\nbissLNT69evDNQKAMSAaPwLvPtQ22iOE3adHOlS77SNNSO4e7VHCauAfJdH2ezxa/2GSmzud/2nv\nJIWtqHq9Xp122mmKi/u/l0xOTlZvb68OHz6spKSkcI0CADAA/yiJHtH4D5Oj3vd1+9y9OvfcqaM9\nSljFxjp08cUXhez5wvrR/5f/VTHw2LKsE36e2NjPb/3qcMSEbrgxpPvQh6M9Qth9dvSApMBojxF2\n0Xje0XjOUnSedzSesxTd533qN9NHe4yw6unu0i8e/51OTWwa7VHC6rNjB7Vv9xgsqvHx8UMK6cDj\ncePGnfDzJCZ+vu+KZf8VuuEAAABgnLD9z1Tp6ek6cuSI/P7/uz6lq6tLCQkJSkxMDNcYAAAAGCPC\nVlSzs7MVFxenPXv2BLe9/vrrmjZtWrhGAAAAwBgStqKakJCgefPmafHixWppadGOHTv06KOP6uab\nbw7XCAAAABhDYgKBQNiu6u7p6dGPf/xj/e53v9OECRN022236aabbgrXywMAAGAMCWtRBQAAAE5U\n2D76BwAAAOygqAIAAMBIFFUAAAAYiaIKAAAAI1FUAQAAYCSjiqplWbrnnntUVFSkWbNm6dFHH/3K\nfVtbW3XNNdcoPz9fZWVleuutt8I4KULBTt5//vOfddVVV6mgoEDz5s3T888/H8ZJEQp28h7Q3t6u\ngoIC7d69OwwTIpTs5P3OO+/o+uuvV15enkpKSvTqq6+GcVKEgp28//CHP+jKK69UQUGBbrjhBrW2\ntoZxUoSSZVmaO3fu1/4ZfbJ9zaii+uCDD6q1tVWbNm3S4sWLtXr1av3+978fsp/P51NFRYWKior0\n9NNPKz8/X//2b/+mnp6eUZgaw3Wiee/du1fV1dUqKyvTs88+q2uuuUYLFy7UO++8MwpTY7hONO8v\nWrJkCe/rMepE8/7kk09066236pxzztFvf/tbXXrppaqqqtLHH388ClNjuE40b4/Ho7vvvlt33HGH\nnn32WU2dOlUVFRXq7e0dhalxMizL0l133SWPx/OV+4SkrwUM8dlnnwXOO++8wO7du4Pb1q5dG7jp\nppuG7Ltly5bAJZdcMmjbZZddFnjmmWdGfE6Ehp28V6xYEbj99tsHbbvlllsCDz/88IjPidCwk/eA\nrVu3Bq677rrA1KlTA6+99lo4xkSI2Mn7N7/5TeCyyy4btK20tDTwwgsvjPicCA07eT/66KMBl8sV\nfPzJJ58Ezj333MCbb74ZllkRGh6PJzBv3rzAvHnzvvbP6FD0NWO+o7p371719/crPz8/uK2wsFDN\nzc1D9m1ublZhYeGgbeeff76amppGfE6Ehp28r776an3/+98fsv2TTz4Z0RkROnbylqTDhw/r5z//\nue677z4F+D9Jxhw7ee/evVuzZ88etG3Lli266KKLRnxOhIadvE877TR5PB41NjYqEAjI7XZrwoQJ\nyszMDOfIOEmvvfaaZs6cqaeeeupr/4wORV+LG/aUIeb1enXaaacpLu7/RkpOTlZvb68OHz6spKSk\n4PaDBw/qH/7hHwYdn5yc/LXffoZZ7OQ9ZcqUQce+++67euWVV3T99deHbV6cHDt5S9IDDzygq6++\nWllZWeEeFSFgJ++2tjZNnz5d9957r55//nmdccYZ+uEPf6jzzz9/NEbHMNjJ+8orr9Tzzz+v66+/\nXrGxsXI4HNqwYYMmTJgwGqNjmK677roT2i8Ufc2Y76j6fD45nc5B2wYeW5Y1aHtPT8/f3ffL+8Fc\ndvL+oo8//ljV1dUqLCzUxRdfPKIzInTs5L1z5041NTVpwYIFYZsPoWUn788++0yPPPKI0tLS9Mgj\nj2jGjBm69dZb1dnZGbZ5cXLs5H3kyBF1dXVp8eLF2rJli6666irV1NRwTXKECkVfM6aoxsfHDxl8\n4PG4ceNOaN+EhISRHRIhYyfvAV1dXbr55psVExOjlStXjviMCJ0Tzbu3t1eLFy/W4sWLh/zhhrHD\nzvs7NjZW2dnZqqqq0tSpU3X33XfrrLPO0tatW8M2L06OnbxXrFihc889V9ddd51ycnK0dOlSjRs3\nTk8//XTY5kX4hKKvGVNU09PTdeTIEfn9/uC2rq4uJSQkKDExcci+Xq930Lauri6lpqaGZVacPDt5\nS1JnZ6duuOEG9ff3a9OmTUM+KobZTjTv5uZmtbe3q7q6WgUFBSooKJAk3X777VqyZEm4x8Yw2Xl/\np6amDrm856yzzlJHR0dYZsXJs5P3W2+9palTpwYfx8TEaOrUqfroo4/CNi/CJxR9zZiimp2drbi4\nOO3Zsye47fXXX9e0adOG7JuXlzfkQtzGxsZBF3LDbHby9vl8uu2223TKKafoiSeeUEpKSjhHRQic\naN55eXn6/e9/r61bt+rZZ5/Vs88+K0m6//77tXDhwrDOjOGz8/7Oz8/X3r17B2177733lJGRMeJz\nIjTs5J2Wljbk+sT3339fZ5xxxojPifALRV8zpqgmJCRo3rx5Wrx4sVpaWrRjxw49+uijuvnmmyV9\n3sAH7rN2+eWXq7u7Wz/96U+1b98+/eQnP5HP59MVV1wxmqcAG+zk/atf/Urt7e1atmyZ/H6/urq6\n1NXVxU/9jyEnmrfT6dTkyZMH/ZI+/8vt9NNPH81TgA123t/XXnut3nnnHa1evVoffvihVq5cqfb2\ndpWUlIzmKcAGO3mXlZVpy5Yt2rp1qz788EOtWLFCHR0duuqqq0bzFBBCIe9rw7+LVuj5fL5ATU1N\noKCgIHDRRRcFHn/88eDaueeeO+i+W83NzYGrr746kJeXF7jmmmsCb7/99miMjJNwonn/8z//c2Dq\n1KlDftXU1IzW6BgGO+/vL+I+qmOTnbwbGxsDV199deC8884LXH311YHXX399NEbGSbCTd11dXeCK\nK64InH/++YEbbriBv7/HuC//GR3qvhYTCHCTQgAAAJjHmI/+AQAAgC+iqAIAAMBIFFUAAAAYiaIK\nAAAAI1FUAQAAYCSKKgAAAIxEUQUAAICRKKoAAAAwEkUVAAAARqKoAgAAwEgUVQAAABjp/wHnKh5W\nOPNiaQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10e49e050>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">25% 92\n",
      ">50% 36\n",
      ">75% 12\n",
      ">90% 2\n",
      ">95% 1\n"
     ]
    }
   ],
   "source": [
    "plt.hist(table_stool.pa(inplace=False).sum(axis=\"observation\")/len(table_stool.ids()))\n",
    "plt.show()\n",
    "print \">25%\", np.sum(table_stool.pa(inplace=False).sum(axis=\"observation\")/len(table_stool.ids())>.25)\n",
    "print \">50%\", np.sum(table_stool.pa(inplace=False).sum(axis=\"observation\")/len(table_stool.ids())>.50)\n",
    "print \">75%\", np.sum(table_stool.pa(inplace=False).sum(axis=\"observation\")/len(table_stool.ids())>.75)\n",
    "print \">90%\", np.sum(table_stool.pa(inplace=False).sum(axis=\"observation\")/len(table_stool.ids())>.90)\n",
    "print \">95%\", np.sum(table_stool.pa(inplace=False).sum(axis=\"observation\")/len(table_stool.ids())>.99)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Unique per family"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Family OTUs</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>[310018]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>57</th>\n",
       "      <td>[2936049, 996116]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>[4426982]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>[1094701, 570755]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>[4325163]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>[1056382, 968121]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>[966091]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>[938243]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>[4318710]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>26</th>\n",
       "      <td>[4353325]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>27</th>\n",
       "      <td>[91024]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1053</th>\n",
       "      <td>[890316]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>30</th>\n",
       "      <td>[188500, 543684, 834883, 1035401, 166076]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>31</th>\n",
       "      <td>[971266]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35</th>\n",
       "      <td>[589597, 548694, 3242243]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>40</th>\n",
       "      <td>[869044]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>41</th>\n",
       "      <td>[1026308, 322234, 102388, 907206, 4312339]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>44</th>\n",
       "      <td>[321764]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29</th>\n",
       "      <td>[1808779, 540826, 278560, 138177, 832652, 4444...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>49</th>\n",
       "      <td>[113760]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50</th>\n",
       "      <td>[539866, 1667532, 4328665, 946945, 1087966, 43...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>53</th>\n",
       "      <td>[778762, 138602]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>55</th>\n",
       "      <td>[4452721]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>52</th>\n",
       "      <td>[151649]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>61</th>\n",
       "      <td>[890460]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>62</th>\n",
       "      <td>[1116372, 114587]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>66</th>\n",
       "      <td>[4365144]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>69</th>\n",
       "      <td>[145053]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>390</th>\n",
       "      <td>[652471]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75</th>\n",
       "      <td>[996126, 350938]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>77</th>\n",
       "      <td>[1009173, 326069]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>78</th>\n",
       "      <td>[4441792]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>37</th>\n",
       "      <td>[489457]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1020</th>\n",
       "      <td>[529839, 4471279, 4296486]</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                            Family OTUs\n",
       "2                                              [310018]\n",
       "57                                    [2936049, 996116]\n",
       "12                                            [4426982]\n",
       "13                                    [1094701, 570755]\n",
       "17                                            [4325163]\n",
       "18                                    [1056382, 968121]\n",
       "19                                             [966091]\n",
       "20                                             [938243]\n",
       "22                                            [4318710]\n",
       "26                                            [4353325]\n",
       "27                                              [91024]\n",
       "1053                                           [890316]\n",
       "30            [188500, 543684, 834883, 1035401, 166076]\n",
       "31                                             [971266]\n",
       "35                            [589597, 548694, 3242243]\n",
       "40                                             [869044]\n",
       "41           [1026308, 322234, 102388, 907206, 4312339]\n",
       "44                                             [321764]\n",
       "29    [1808779, 540826, 278560, 138177, 832652, 4444...\n",
       "49                                             [113760]\n",
       "50    [539866, 1667532, 4328665, 946945, 1087966, 43...\n",
       "53                                     [778762, 138602]\n",
       "55                                            [4452721]\n",
       "52                                             [151649]\n",
       "61                                             [890460]\n",
       "62                                    [1116372, 114587]\n",
       "66                                            [4365144]\n",
       "69                                             [145053]\n",
       "390                                            [652471]\n",
       "75                                     [996126, 350938]\n",
       "77                                    [1009173, 326069]\n",
       "78                                            [4441792]\n",
       "37                                             [489457]\n",
       "1020                         [529839, 4471279, 4296486]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1, 34)\n"
     ]
    }
   ],
   "source": [
    "unique_family_site = pd.DataFrame(index=[\"Family OTUs\"])\n",
    "\n",
    "for family in set(meta.family):\n",
    "    infamily = meta.loc[meta.family==family].index\n",
    "    outfamily = meta.loc[meta.family!=family].index\n",
    "    infamily_table = table.filter(infamily, inplace=False)\n",
    "    outfamily_table = table.filter(outfamily, inplace=False)\n",
    "    infamily_present = infamily_table.sum(axis=\"observation\") > 0\n",
    "    outfamily_absent = outfamily_table.sum(axis=\"observation\") == 0\n",
    "    unique_family = infamily_table.ids(axis=\"observation\")[np.logical_and(infamily_present, outfamily_absent)]\n",
    "    if len(unique_family) > 0:\n",
    "        unique_family_site[family] = [unique_family]\n",
    "\n",
    "display(unique_family_site.transpose())\n",
    "print unique_family_site.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Unique per family per site"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Site</th>\n",
       "      <th>Family OTUs</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>30</th>\n",
       "      <td>forehead</td>\n",
       "      <td>[543684]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29</th>\n",
       "      <td>forehead</td>\n",
       "      <td>[832652]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50</th>\n",
       "      <td>forehead</td>\n",
       "      <td>[1667532, 4328665, 140863, 4305723]</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        Site                          Family OTUs\n",
       "30  forehead                             [543684]\n",
       "29  forehead                             [832652]\n",
       "50  forehead  [1667532, 4328665, 140863, 4305723]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sites = ['forehead', 'tongue', 'stool']\n",
    "unique_family_site = pd.DataFrame(index=[\"Site\", \"Family OTUs\"])\n",
    "\n",
    "for family in set(meta.family):\n",
    "    meta_infamily = meta.loc[meta.family==family]\n",
    "    meta_outfamily = meta.loc[meta.family!=family]\n",
    "    for site in sites:\n",
    "        infamily_site = meta_infamily.loc[meta_infamily.env==site].index\n",
    "        infamily_table = table.filter(infamily_site, inplace=False)\n",
    "        outfamily_site = meta_outfamily.loc[meta_outfamily.env==site].index\n",
    "        outfamily_table = table.filter(meta_outfamily.index, inplace=False)\n",
    "        infamily_present = infamily_table.pa(inplace=False).sum(axis=\"observation\") > 1\n",
    "        outfamily_absent = outfamily_table.pa(inplace=False).sum(axis=\"observation\") == 0\n",
    "        unique_faimly = infamily_table.ids(axis=\"observation\")[np.logical_and(infamily_present, outfamily_absent)]\n",
    "        if len(unique_faimly) > 0:\n",
    "            unique_family_site[int(family)] = site, unique_faimly\n",
    "display(unique_family_site.transpose())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['797.50.HB.Forehead', '797.50.HA.Forehead'], dtype=object)"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "table.ids()[table.data('4305723', axis=\"observation\")>0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:sourcetracker]",
   "language": "python",
   "name": "conda-env-sourcetracker-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}