{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "ImportError",
     "evalue": "No module named pymc",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mImportError\u001b[0m                               Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-1-17bed78f2f88>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[1;32mimport\u001b[0m \u001b[0mpymc\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      2\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mnumpy\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mmatplotlib\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpyplot\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mscipy\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mstats\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0mstats\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[0mget_ipython\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmagic\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34mu'matplotlib inline'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mImportError\u001b[0m: No module named pymc"
     ]
    }
   ],
   "source": [
    "import pymc\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as stats\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 282,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "num_obs = 100\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 283,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_obs = np.random.poisson(1./50., size=num_obs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 284,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma_r = np.random.normal(loc=.01, scale=.003, size=num_obs)\n",
    "sigma_f = np.random.normal(loc=1./.05, scale=5, size=num_obs)\n",
    "sigma_m = np.random.normal(loc=.001, scale=.005, size=num_obs)\n",
    "\n",
    "m_r = np.random.normal(loc=.03, scale=.02, size=num_obs)\n",
    "m_f = np.random.normal(loc=1.7, scale=.2, size=num_obs)\n",
    "m_m = np.random.normal(loc=1.7, scale=.3, size=num_obs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 285,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs4AAAGICAYAAABRBWbAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XmcJHV9+P9XL7vLtTu7bNDlWh1AUSQooBwRkQHB4BHE\nIyiKEjxiokE0RAXMTwZNFG+MRL8/5XBVRBGUgNeXRRnFoASU+xBRkEN2QVnYIXIs7Hz/eFczNb3d\n01Xd1Ud1v56PRz2mu7rq05+ufnfNuz71qU+BJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEnr\nOR1YBVybmrcEWAHcDFwILE69dhzwG+Am4MVdqqOUxzLgYuB64DrgXcn8ceBO4MpkOqgXlZNmUW9/\nvAfwP0TMXg7s3oN6Sc0Yuxoa+wC7MjPYPw68L3n8fuCk5PGzgKuAecAocAswpyu1lLLbAtglebwA\n+DWwI3AC8M+9qpSUQb398QTw18njlxAHhVK/MXY1MJoltpcAq2vmHQwsTx4vBw5JHr8COAtYC9xG\nJM57FFJLqTgriQM8gAeBG4Gtk+eVntRIyqbe/vhuYFHyeDFwV1drJGVj7GqojDLzKDEd/JXU888B\nb0i9dirw6o7WTGrPKPB7ouX5BOKA72rgNGZ2QZL6xSgz98dPBe4Abie6Gi3rQZ2kLEYxdjUA2u1K\nMZVMs70+w/bbb19dx8mpnanaatyqBcA5wNFEy/MXgG2Jbhx3A5+qXcHYdSpouoXinEb0038K8B6i\nL+l6jF2nAqYi4xYyxK5x61TQVHTsNjXKzKPEm4h+ogBbJs8Bjk2mqh8Ce9Ypb6pMTjjhhFK9T6vl\n5F0vy/LtLjPba8SPoVXzgP8LvLvB66PMjPlCY7fImOp13FhOfm3Gbm1srkk9rgAPNFivsPrX0+n9\nZDf2w2X/DJ0uv824bTV2O/qZpqbK970AU4ynpn1rnqenPvyf1YvyC4jdGVppcT4fOCJ5fARwXmr+\n64D5RMvd04krZkttbGysVO/Tajl518uyfLvLdGjbV4iWjhuAk1Pzt0w9fiX1E2ep39wC7Js83p8Y\n7UgqA2NXpTS3yetnEYG9OdEX6YPEKBpnA28h+oQemix7QzL/BuAx4B0UnOX3golz68v3aeK8N3A4\ncA0xDBLA8cBhRDeNKeBW4O2deHOpDfX2x38P/CewIfBQ8lzqN8auBkazxPmwBvMPaDD/I8kk9auf\nUf9Myw+6VYEiDwh6fcBlOV3VaH9cr0tcV3V6+3Rj+5f9M/RJjDZi7HbKaGeLhwHYRgXrxfBbSZcT\nqXWVSgW6H7/Grtpm7KqMjNv+UKlU4nZdWYyD26/42PUGJZIkSVIGJs6SJElSBibOkqS+tWRkhEql\nkmman3G5VpZfMjLS600h5TOHXL+HkcXGeBbNLg6UJKlnVk9OZh6eqUK+oZzyLF+ZnMxRstQH1pG9\nPzQwOW6MZ2GLsyRJkpSBibMkSZKUgYmzJJXD6cAq1r+r5VHAjcB1wMe6XSmpCeNWA8U+zpJUDmcA\nnwO+kpq3H3Aw8GxgLfCkHtRLmo1xq4Fii7MklcMlwOqaef8IfJRIPgDu7WqNpOaMWw0UE2dJKq+n\nAy8EfgFMAM/raW2kbIxblZZdNSSpvOYCmwF7AbsDZwPb1VtwfHz8icdjY2OMjY11vnYqrYmJCSYm\nJjpVvHGrjulw7Hb9vvPgvedVgKLvPZ+Rsau2tRm7o8AFwM7J8x8AJwE/SZ7fAuwJ/KlmvdLGbqVS\n6Y9xnIGybsMiGLf9oVKpZB+beZxc4zgzPpgxXnS+YFcNSYUbGVmS645VlUqFkZElva52GZ0H7J88\n3gGYz/rJh9RvjFuVll01JBVucnI1+dr+YHKyFyfASuUsYF/gL4A7gA8SQ32dTgz19Sjwpp7VTqrP\nuNVAMXGWpHI4rMH8N3a1FlI+xq0Gil01JEmSpAxMnCVJkqQMTJwlSZKkDEycJUlqYi5kHiFmychI\nr6srqUPaSZyPA64nror9OrAhsARYAdwMXAgsbreCUsGWARcTsXsd8K5kvrErqaHHiHFiskyrJyd7\nVEtJndZq4jwKvA3YjRjQfAPgdcCxRPKxA/Cj5LnUT9YC7wF2Iu5a9U5gR4xdSZLURKuJ8xoiAdmE\nOIO1CfAH4GBgebLMcuCQdisoFWwlcFXy+EHgRmBrjF1JktREq4nzfcCngNuJhPl+orVuKbAqWWZV\n8lzqV6PArsBlGLuSJKmJVm+Asj3wbiLxeAD4FnB4zTLV7l7rGR8ff+Lx2NgYY2NjLVZDw2JiYoKJ\niYkii1wAnAscDdR2SDR2VZgCY/d04GXAPUQXubRjgE8AmxMNG1I/MXY1MFq9x+1rgQOBtybP30j0\nF90f2I84Hb4lcRHWM2vWnZqayncrXqlWpVKB1uN3HvBd4AfAycm8m4AxjN1CxPeTd1tVGIbt20bs\n7kN0L/oKM5OPZcCXgGcAz6V+8lHa2K1UKpkjKW/U5Vk+97Il3d6NtLnPbTV2Sxu3nVKpVGA848Lj\nZF82WX4Qt3ebsbueVrtq3EQkyhsnlTkAuAG4ADgiWeYI4Lx2KygVrAKcRsTryan552Psqr9dAqyu\nM//TwPu6XBcpD2NXA6PVrhpXE0eOVwDrgF8BXwQWAmcDbwFuAw5tv4pSofYmuhVdA1yZzDsOOAlj\nV+XzCuBOIp6lMjF2VUqtJs4AH0+mtPuI1mepX/2MxmdajF2VySbA8US3uaqGpyPtn688OnBdSVrm\n2B30uB1ZPMLkA477XaQOx25xfT5ysM+S2lZ0n6WMjN2M7OPcWJuxO0p0ids5mS4C/py8tg1wF7AH\ncRFWWmlj1z7O/aGAfe4o+WO3tHGbVa4+y5Cv33KeZZPlB3F7F50vtNPiLEkFmlvdweWycOFmrFkz\nlBfjX8vMYRNvpfHFgVI/MXZVWu3ccluSCpTnpsbT0+RkvWuOBtJZwKXE3S3vAI6seX3wmoo0KIxd\nDQxbnCWpHA5r8vp2XamFlJ+xq4Fhi7MkSZKUgYmzJEmSlIGJsyRJkpSBibMkSQWaSwyBlXVaMjLS\n6ypLMCd73I4sHt6Y9eJASZIKVB0fJqvKpDfAUB9YR+ZxnyfHhzdmbXGWJEmSMjBxliRJkjIwcZak\ncjgdWEXcda3qE8CNwNXAt4FFPaiX1Iyxq4Fh4ixJ5XAGcFDNvAuBnYDnADcDx3W7UlIGxq4Ghomz\npFmNjCzJNUJApVLpdZUH1SVA7f3FVxCX9ABcBmzT1RpJ2Ri7GhgmzpJmNTm5mhgjIM+kHngz8P1e\nV0JqgbGr0nA4Okkqvw8AjwJfb7TA+Pj4E4/HxsYYGxvreKVUXhMTE0xMTHTjrWaNXeNWeXU6dntx\nTnVqasoWKbUn6Q7Q7fgdytiNbZ33c3drnVivTN9Lm7E7ClwA7Jya93fA24AXAQ83WK+0sVupVDJH\nRd4IyrN8p5Z9Yvk+/34K2OeOkj92Sxu3WVUqlcxjJwOxbNbl8yzbQtll+W6KzhdscZak8joIeC+w\nL42TZqkfGbsqJfs4S1I5nAVcCjwDuIPoF/o5YAFxodWVwOd7VjupMWNXA8MWZ0kqh8PqzDu967WQ\n8jN2NTDaaXFeDJxDDGB+A7AnsIQ4eryZGKNxcbsVlApWbyD+ceBOotXjStYfb1SSJKmtxPmzxPAx\nOwLPBm4CjiUS5x2AHyXPpX5SbyD+KeDTwK7J9MNuV0qSJPW/VhPnRcA+TJ9qeQx4ADgYWJ7MWw4c\n0lbtpOLVG4gfejPCjCRJKpFWE+dtgXuJ1rtfAV8CNgWWEqfBSf4ubbeCUpccBVwNnIZdjCRJUh2t\nXhw4F9gN+CfgcuBk1u+W0fAWYg5orrw6PKD5F4APJY8/DHwKeEu9BY1d5dXFG0lIkjqs1dPTWwA/\nJ1qeAV4AHAdsB+wHrAS2BC4Gnlmz7sAPaK7O68BNJLK8NpSx6w1QiuXNe/LxBij9wbjtDG+A0nlF\nx26rXTVWEmMx7pA8PwC4nkg4jkjmHQGc11btpO7YMvX4lcwccUOSJAlobxzno4AzgfnAb4EjgQ2A\ns4nT3LcBh7ZZP6loZxF3qtqcOPg7ARgDdiEaiW4F3t6rykmzOB14GXAP02dElgDfBJ7K9D73/l5U\nTpqFsauB0c5wdFcDuwPPAV5FjKpxH9H6vAPwYvwRqP8cBmxFHPAtI3bobyKGVHwOMRLMqoZrS71T\nbyhFhwBVGRi7GhjecluSyqHeUIoOAaoyMHY1MEycJam8HAJUZWXsqpTa6eMsSeofDYcABYdSVD5d\nHkax74evHVk8wuQDk5mWXbhoIWvuX9PhGqmRTsduL+6WNvDDy6jzHBqpexyOrlgFD6V4E3Fx62xD\ngEKJY9fh6PpDAfvcUfLHbt/Eba5h48azf58OR9d5/TIcnSSp987HIUBVTsauSsnEWZLK4SzgUuAZ\nxFCKRwInAQcCNwP7J8+lfmPsamDYx1mSyuGwBvMP6GotpPyMXQ0MW5wlSZKkDEycJUmSpAxMnCVJ\nkqQMTJwlSZKkDEycJUmSpAxMnCVJkqQMTJylITEysoRKpZJ7kiRJwcRZGhKTk6uJGwHnnVQCxwHX\nA9cCXwc27G11pEyMW5WOibMkldso8DZgN2BnYAPgdb2skJTBKMatSsg7B0pSua0B1gKbAI8nf+/q\naY2k5oxblZItzpJUbvcBnwJuB/4A3A9c1NMaSc0ZtyolW5wlqdy2B95NnPp+APgW8AbgzPRC4+Pj\nTzweGxtjbGysW/VTCU1MTDAxMdHJtxjcuJ2DF1b3UKdjtxff7NTUlBccqT3JTqnb8Vvq2I1t1kr9\nW1mvW+vEemX6XjoQu68FDgTemjx/I7AX8M7UMqWN3Uqlkjkq8kZQnuU7tewTy/f59zPscVupVGA8\n48LjdGbZPiu7X76bZoqOXbtqSFK53UQkHBsT/xwOAG7oaY2k5oxblVK7ifMGwJXABcnzJcAK4Gbg\nQmBxm+VLRTsdWEUMf1Rl3KrMrga+AlwBXJPM+2LvqiNlYtyqlNpNnI8mjhCr7fXHEgnIDsCPkudS\nPzkDOKhmnnGrsvs4sBMxrNcRxGgFUr8zblU67STO2wAvBU5luu/IwcDy5PFy4JA2ypc64RJgdc08\n41aSJDXVTuL8GeC9wLrUvKXEaXCSv0vbKF/qFuNWkiQ11epwdC8H7iH6N481WKbh/XpLObyMeqoL\nQyNVzXqfaWO3H83NPfTTwoWbsWbNfR2qz0xdjF1JGaxcuZKf/vSnmZdfvNjLXjSt1eE5PkIMHfMY\nsBEwAnwb2J1IpFcCWwIXA8+sWbdvhpdRebU5vMwocUHrzsnzm2get1Dy2B3k4ehaea9efZcOpZiP\nw9H1h0GK249+9KOc+NkTmb/V/EzLT145GQ/GM77BeIeW7aeyP8TM/gZNLFy0kDX3r8lRmeIUHbut\ntjgfn0wA+wL/QiTSHyc6+H8s+XteuxWUuuB8jFtJGgpTU1Os3WEtj7zokWwrXEVrx/SDbB25kvLJ\n8clO1aTrihrHuRpSJxEDmt8M7J88l/rJWcClwDOAO4AjMW4lSVIGRdxy+yfJBHHv+QMKKFPqlMMa\nzDduJUnSrLxzoCSV32LgHOBGYmz9vXpbHSkT41alU0SLsySptz4LfB94DbFf37S31ZEyMW5VOibO\nklRui4B9iAtbIUY7eqB31ZEyMW5VSibOklRu2wL3EreTfw7wS+Bo4M+tFvjoo4/yi1/8gnXrso03\nNTIywm677dbq22k4FR63UjeYOEtSuc0FdgP+CbgcOBk4FvhgeqE8N+/58Y9/zOsPOYRnb7RRpgpc\nsmYNDz/yCPPmzcu0/JKREVZPDs7wVIOoCzfuKTxuJeh87Jo4S1K53ZlMlyfPzyESkBnSCUgzjz/+\nOH+10UZ874FsZ87nz5mT6wYeqycnc91MRN1Xm6SeeOKJRb9F4XErQedj11E1JKncVhJjku+QPD8A\nuL531ZEyMW5VSrY4S1L5HQWcCcwHfkvc2Efqd8atSsfEWZLK72pg915XQsrJuFXp2FVDkiRJysDE\nWZIkScrAxFmSJEnKwMRZ6rGRkSVUKpVc08jIkl5XW1JB5kLm3/6SkZFeV1fKb072GB9Z3N8x7sWB\nUo9NTq6GzKPaVtdxdFtpUDxG9j1AxRvHqIzWAePZFp0c7+8Yt8VZkiRJysDEWZIkScrAxFmSBsMG\nwJXABb2uiJSDcatSMXGWpMFwNHADeTvMS71l3KpUTJwlqfy2AV4KnAp45ajKwrhV6Zg4S1L5fQZ4\nL3HtulQWxq1Kp9Xh6JYBXwGeTJxe+SLwH8AS4JvAU4HbgEOB+9uupdQdtwFrgMeBtcAePa2NlM3L\ngXuIfqJjjRYaHx9/4vHY2BhjYw0XlZiYmGBiYqKTb2HcqiM6HbutJs5rgfcAVwELgF8CK4Ajk78f\nB94PHJtMUhlMETvw+3pcDymP5wMHE6e8NwJGiIaNN6UXSicgUjO1SeqJJ55Y9FsYt+qITsduq101\nVhJJM8CDwI3A1sSPYHkyfzlwSFu1k7rPfnYqm+OJs4DbAq8DfkxN8iH1IeNWpVREH+dRYFfgMmAp\nsCqZvyp5LpXFFHARcAXwth7XRWqVoxOojIxblUK7t9xeAJxLDCdTe4/EKRr8EOyzpLy60N8OYG/g\nbuBJRJejm4BL0gv0T+zOpVKxcbybRkaWJLdHz2fjjRfwvvcd04Ea1fWTZJLKxLhVabSTOM8jkuav\nAucl81YBWxBdObYkOv6vxz5LyqsL/e0gkmaAe4HvEBcHNkyce+sx8jfQmGi3I5Lm/I1iDz1UmRE3\nHYpdSVIXtNpVowKcRgxafnJq/vnAEcnjI5hOqKV+twmwMHm8KfBi4NreVUeSJPWbVluc9wYOB64h\nhpIBOA44CTgbeAvTw9FJZbCUaGWG+F2cCVzYu+pIkqR+02ri/DMat1Yf0GKZUi/dCuzS60pIkqT+\n5Z0DJUkqiblApVLJNC0ZGel1daWB0+6oGpIkqUvyXBZcmawd7EpSu2xxliRJkjIwcZYkSZIyMHGW\npHJbBlwMXA9cB7yrt9WRMjN2VTr2cZakclsLvAe4irib6y+JO1/e2MtKSRkYuyodW5wlqdxWEokH\nwINE0rFV76ojZWbsqnRMnCVpcIwCuwKX9bgeUl6jGLsqAbtqSNJgWACcAxxNtN7NMD4+/sTjsbEx\nxsbGulUvldDExAQTExPderuGsWvcKq9Ox66JsySV3zzgXOBrwHn1FkgnIFIztUnqiSee2Km3mjV2\njVvl1enYtauGJJVbBTgNuAE4ucd1kfIwdlU6Js6SVG57A4cD+wFXJtNBPa2RlI2xq9Kxq4akITSX\nSqXS60oU5WfYCKJyMnZVOibOkobQY8BUznUGJtGWJLXIIz1JkiQpAxNnSZIkKQMTZ0lS27bYfHMq\nlUqmSZIamkPmfUmlUmFk8UhXq2cfZ6kADz74IOPj/8Yjj6zNtd6GG87rUI2k7lo9OZm517ips6SG\n1gHj2RefHJ/sVE3qMnGWCnDrrbdyyilf5pFH3ptrvQ03/HiHaiRJkopm4qy+NTIywuRkd48k27Hh\nhk/ikUeOybXO/Pmn8cgj93SoRpIkqUid6ON8EHAT8Bvg/R0oX0OiB0lzl2J3og/LspzulNMxPd3v\nrpvKO7RfPhMdLb0771H28juk9/nCrZbf8/foxmcoUNGJ8wbAKcSP4VnAYcCOBb+H1AldjN2JPizL\ncrpTTkf0fL87ZeI88OV3QM/jFoDbLL/n79Hp8gtWdOK8B3ALsRnWAt8AXlHwe3TVxMREqd6n1XLy\nrpdl+aKW6ZKBi10NDWNXZWTcqpSK7uO8NXBH6vmdwJ4Fv0dXTUxMMDY2Vpr3abWcvOtlWb6oZbqk\nrdidM2cODz10K4sWvbjpsg8//Fs22uhSAP7859/nrKa0nsL3u3PmzOGyhx/mxYsWZVp+3QMPtPN2\nGk49yxfmzJnDvBvnsdGfNuLhex5mo3s3mnX5NaxhKvedRqVsXg18KfX8cOBzNcvcQtzr1smpnekq\nimXsOnVruoViGbtO3ZiMW6eyToXGbtEtzncBy1LPlxFHkWlPK/g9pSIYuyorY1dlZNxKRCL+W2AU\nmE+0CnpxoMrA2FVZGbsqI+NWSrwE+DXRNH5cj+si5WHsqqyMXZWRcStJkiRJkgZLlsH3/yN5/Wpg\n1wbrnpZ6/BPgxmT5bwOLWiynWp9jgHXAkjbKOSqp03XAx1os53PA/wBXApcDuzcp53RgFXBtTbm/\nTZa/GbgQWJxaJ0tZ6frl3daNysm7rWcrJ8+2rt1GVXtQf1u34jbgmqSs/2mjnLR69V4CrKD+91pE\n+eNE/9crk+mgNspfBlwMXE98T+9K5hf1GRqVP04xn2Ej4DKia8MNwEeT+UV+B43eY5xiPkPWujba\nR/8tsX0fB3ZLzR8FHiJifg1wf8HlQ7SM/xb4X2JbtFJ+o88/CjwCPJxMP29Qbpb9d6P3+2Xyt93/\ne+l1x5mOi1uI0VJaKb/RPrGo+jcqP13/dvcv0sDagPiBjwLzqN+/7qXA95PHewK/qLPuhsSO7kVJ\nObcQg/kDnJRMrZRzFbAf8EPivkpL2qjPpcljgCe1WM6DwFuT5V5CJAaNygHYh9hhXVtT7ieBPxDb\n+v3J9iFHWenvLM+2nq2cPNu6WX2ybuvactImgL9OHle3dauqn6dI9er9ceB9yeP091pU+ScA/9xG\nmWlbALskjxcQXQV2pLjP0Kj8Ij/DJsnfuURMvYBiv4NG71HUZ8hS19n20c8EdiB+G7WJ87UdLP9Z\nyXKfJA4mbgGObaH8RvXbjtjv1lunKsv+u9H7bQD8EfhCweVX46Kd/6vQeJ9YRP1nKz9XXHfiltut\n2hY4FfhWryuiUtoUWA58EXh9huWzDL5/cFImROvLYuKfYnrd3ZK/uyflnJqsV11nmxbL+QbwWaZ3\nrq3W589Ey8naZL17WyznFiK5JFn+rlnKAbgEWJ08Tpf7cuKo/xXJuoc0+WzpsjZh5neWZ1vPVk6e\nbT1bOXm2de02Srub6dbz6rZuR6XN9WvVq3f6c6a/16LKh+I+x0qmh7N8kDhDsDXFfYZG5UNxn+HP\nyd/5RDKxmmK/g0bvAcV8hix1nW0ffRPR6tjt8l8BnEXsx05O1r+mhfIb1e85wKMN1qn32Rrtvxu9\n3x5JffcruHyIuGjn/yo0/u0XUf/Zyq/WP5N+SpxvZbpFS8rrVcDZwN8znUzNpt7g+1tnXGar1Pyt\niRbUrWuWAXgzceTbSjlPJk4TXtNmfUaAZxNH3RPA81os5wxiu94OfII4XZllG9bWeynxT2lr4pTZ\n0iafLW3uLMs029azlZNnW89WTp5tPZtjgU8xc1u3agq4CLgCeFsb5TSzlPg+Yeb3WqSjiFOvp9Fe\nN4S0UaIF6jI68xmq5VdbvYr6DHOI5HwV091Ciq5/vfeAYj5Dlrq28tuBaIR7GvBNopW8yPK3Spar\n1v9OYOMWym/0+bcguslcSexDNqlTpyz770bvtzXRzWRpnWXaKR8iLs4iYn5xg2VmK382RdS/mcxx\n3YnEuVEfkiz9SaW0PLGU/rE8nqHsqYx1aHYU2qicDxAtB19voZyNiST1sjr1yFufDYhuF3sB7yUO\nLlop523AT4GnAO8hvpt65dTbHo220VTNa1nKqqfZtp6tnE3Ivq2b1Sfrtm5WzmlEv9jabd2KvYnE\n7SXAO4lTlZ1W+70W4QtEQrQL0SL/qQLKXACcCxwNTNa8VsRnWACck5T/IMV+hnVJOdsAL2T6bFBV\nEfWvfY8x8n2GFcS+u3aqbdhoVNdm9V9BnBH7Zqrs7xL7qkni1PvXie+hqPIPY/0+z1nrX5lluer8\n+5I675rU/x1Md/2qLauZeu83VfN3tnXzqMbFMUQf82ax3eq+vtX6N1s+12+zE4nzGazfsXoD4JRk\n/rOI4NsReCPwGeJIRqqVJ5buZHow/SxxnWXw/dpltkmWSc+/i0jaq+suI46GXwq8ocVytk/KOJQ4\nE7MNcUHEn1qoz5+Zbu26nPhn+McWynkakThDJAN7NPhc9boVpJdbxfT3tSVwT51lGpW1lvW/syzb\nerZy8mzrZvXJuq2bdb3YA/hO8ri6rVt1d/L33qTMdsqazSqmT4emv9ei3MN0gnEq7X+OeUTS/FXg\nvGRekZ+hWv7XUuUX/RkAHgC+BzyXzn0H1fd4Hvk+w4HAznWm8zPWtdk++kDid/baVNl/SbR6riLO\nmP0W+KsCy/8k8Vuq1n8b4iK+LOWnf/uNPv/vmW5N/RWx7R/JUG7t/rvR+91F7Mur79fq/73adatx\ncRdxUeYedZaZrfxm+8R269+s/Fy/zU4kzvX6kDTq9/JVokXlD8QFLP+HyPhtkRbki6VvE7dw/Tyx\nY27mCuDpTA++/9o6650PvCl5vBdxhfaqmnWvBp6azJtPdDfaNanTwy2WczOxw9+bOAq+k2jlOLuF\n+mzCdAv8Dsm8b7VQzjridwqwf1LHRp+rVrrc7wJHJusewXRSkaWsh5j5neXZ1o3K+TX5tvVs9cmz\nrWdzC7Bv8ri6rVuxCbAwebwp8GLWP3tTlOr3CTO/16JsmXr8Str7HBWiVf8Gop9qVVGfoVH5RX2G\nzZk+lbwxkeBdSbHfQaP32CK1TDufIUtds+yjYWbr4uZE48r5RG7xdOD5BZZ/PvA6Yj/27mT9Z7dQ\nfqPPf2tqnWcQCeCZNeVm+X/S6P2uSOp78SyfuZXyq7F9BbAT8LsWyp9NEfWfTZH7l5aN1rzxa2h+\nT3qpnlE6F0v1Bt9/ezJVnZK8fjUzT9Gl1z0j9fiPRKvB7cn0+RbLSdfnXqZHRshbzgeIA9S7knqN\ntVjOKUR3hjuIRHPXJuWcRSTajyTrfDop63fJ8jcTicW7U+tkKeteohW1lW3dqJy827pROXm3de02\nOjKZ/zymhwL7OTOHU8pj26SMq4hh0Yq6wUS13o8yXe8lRF/qIoZCqy3/zcBXiD7oVxNJRjv9d19A\nHAhexczhp4r6DPXKfwnFfYadidbIq5Ly3pvML/I7aPQeRX2GRnXdimjdrmp0g5RXErHxEHEx5g+S\n+a8mYr1RwFTYAAAgAElEQVQ6HN1dBZcPcDyxH6s3HF3W8ht9/lcRjULV4ei+nMxv5f9So/f7JdHN\nsd3/e+ny03FxafJ6K+U32icWVf9G5Re5f2nZKDOTnVdj4qzWjGIsSZKkPtCtUTWy9CeVsjCWJEnS\nQBllZivhXOL07ijRL6XeoNVSPaMYS5IkaUA16kPSqG+M1IixJEmSJEmSJEmSJEmSJEmSJEmSJEmS\nJEmSJEkaPvOBE4i7Na6rmR6mvTvoSZIkzcpERGUxH/gRcUvkvwL2Am4E/g/wFOAvelc1yX2pJA06\nExGVyUnAL4ANUvP+Cfhlb6ojPcF9qSQNARMRlcUi4M/AS2vmHwNcM8t6hwAfBb4EHAbcR9yJNIt3\nAMfXmb8bcDnw7ozlaPC5L5WkAddqItKqvwT+Gzi6A2Vr8L2cuKvo/Jr53yZa9Rq5CdgOOChZ9zc5\n3nMpsGmD184nEmipk/vSRgdvnfR68h1gStJQaDURacfXiARayutvgbtq5m0NTAI7NlhnW+CS1PN9\ngLMKqMsGwO1ApYCyVH6d3JfOdvDWKRuS7wBTKoQXCajftZKItOtXHSpXg+/JRCvYkuT5fOCHwDsb\nLL8XcDbwE+A4YA7wr8A/JK9vBLwd+AYwN5l3KrATsHHy2jnAvOS1hcAngNcB/w58r4DPpMHQi31p\nJxV1gCll5kUCKoOsichORLJxBfCfwGXAa4hWiROAI4mEZJ/UOhsB48lrxwLLgO2BbwIjRFLyJ+AF\nxX4kDbgXA18HTiSS3Jc3Wf544O9Tzy8k4hng1UQM3whsksz7LfE7eA3T3Tqqr50PvDB5/GHg/a1+\nCA2crPvSpxNd1T4JvBJ4ExHPzwHeAHwQOCJZdmMidqsHb08n+tR/nOi3/0bgq8lyVZsRv43XEwd3\nm9N4/70A+BjwWuADwMmpctIHmFJXeJGAyiJLIrI/sZO9hzjoW0r8YzgXeFGyzCbM7Mv3HWC/1OO/\nJpLoU4iEBSKZnlPQ55DqORPYO3k8j+heUbUJsCvw3eT5KNPdOjYFnk/ELsDuwA2pdc8lGkSkqiz7\n0jcR+89VTDegnUX0hQbYGbggeZw+eNs0te7dTCfo3yJiGCLf+G9gh+T5O4iuSo323+kDwX9jZpKf\nPsCUOq6IiwQ+CKzBpEL940Dgy6nnexMtGFXbEC3IEK3Il6de2zr5+2WiRe/LRKuz1GlXMt0t7vlE\nspBuoftXplukXw98hOnY/DzwquT5Mcnzqt8RifjCjtRag2pToqX5lNS8q5hOhF9PdAeqLps+eKuu\n+9nUujcyHYN/Q+QYhwJvA56XWq52/117IPhtphPw2gPMoWci1nn7EEd+F9XM3xu4NGMZ5xMtH+sK\nrJfUjn2YeZHV84GLU8/3J/qSQnRPmki9Vu37twfR328d8Bai9WOzDtRVgui3vAi4P3m+J9Ei99rU\nMlsS16EAHEAc8FXPiPw10Y/5UOAB4N5k/r7AHcCziRY9Kav/Zea+dDMiUb0vef43RMyNJMseTnTH\nGEmt+7Nk2acTLdePEF3jngn8X6Jf/5eY2bBRu/9+ITP30bsAVyfvszuRzKcPMIeaiXPnbQz8EXg0\nNW9r4ojvs3XXWN++zExKpF7bm5k73lXAQ8njDYkWjvcmz1emXgN4KrGjXk38Ln6T/H0xcXZG6oQd\ngF+nnl9BJLqrUvPOBF5GtPRdRyQNv0te+znR3/QC4gLCJxNJ9wjwB+I30YlhGzXY0vvSfYiDOYgG\nt72IxPh1ybz0wVt13Wri/AqiNfo1RPJ9E7A29T67Md1to3b/nT4Q3I9o3NgJeBr1DzCH2tzmi6hN\nPyGS5yXEUeR84DTiAqkba5bdMJl/O5Fc/4QI7n2IU9qHEi177yYOev6RSMqfCnyBuGJX6oaNgZtT\nz79G9Iv7O+Kiv3cQMQuRjJxEJB1TwINEa8kPk9d/QCTaFxKtJVIn7MzMVrVLmJk8QJwFbHQm8PCa\n5/+YenwBUn5ziOT2D8nzHYHvJ48fB/6H2Dd+M5lXPXj7r9S6dyev3Uyc3buVyAUuIPKFNxHDJN7N\n9ChGtfvvbzA9QsxDybL7E417mxIXKnpQqK7KeuV3o4urfk20bkD8gPYkjjzfl8z7LOUc4kaSOm0v\n4pT1ScTpa0nSAGh0cdUOzGzN+DGRXD+XaME+n+lkW5I0087ERX7HNFtQkpqxq0b/aHRx1c7EUHYQ\np1eeTZxuWU2M9fhqouP/dl2rqSSVx7XJJElt8+LA/tHo4qpJouM+xMgDXySS6R8SV3KfzPTFBJIk\nSdLAm0OcTvw74i5UOyfzNyAuuvoH4s5rEN04jibuFHQM8Iwu1rPsNiLulHQVMW7lR5P5S4AVxAUT\nF+It0NV/lhFnpa4nRnx4VzJ/HLiTGKP4SuCgXlROmsXpRONQuuV/D+LityuJYf9270G9JEkZVG+Z\nO5foBvMC4ral1Yst309cSCT1ky2I8VUh7vz1a+Ki4BOAf+5VpaQM9iFuqJFOnCeIi9wBXoJDrqok\n7KqhYVQdK3g+0aK/GjgYWJ7MXw4c0oN6SbNZSZwpgRjS70am78JY6UmNpGwuIfazaXcTN6SBOMN3\nF5KkvjSHSEAmiZZmmLlTr7D+Tl7qJ6PEHe4WEC3OtxF3+joNuxmpP40ys8X5qcR1OrcTXY2W9aBO\nkqQcFhFdNfZj/UT5vtqFt99++yniBh5OTu1Mt9CeBcTQldWzIk8mDvYqxPUQp9WuYOw6FTC1G7ej\nzEycLwJemTz+W+IakxmMW6eCpnZjV1LK/wf8C3F70i2SeVsmz2tNFe2EE04YyjLLUMdOlUnsyFs1\nj7iZx7sbvD5K/aHXCv8cgxoTwNTULNMJNc+L2Lb98tln02bc1ovNNanHFaZHj+pI3Ba1PbpRDjDF\neMZp3+Rvm9uqyHjpt21dQOzOYB9nDZvNmT6VvTFwIHFV9/nErUxJ/p7X/apJs6oQrck3EMNQVm2Z\nevxKHLNY5XALsG/yeH9m3gJa6lveAEXDZkvi4r85yfRV4EdE8nw2MVb2bcChPaqf1MjewOHANUS8\nAhwPHEaMtjEF3Aq8vSe1kxo7i0iSNyf6NX8Q+HvgP4n7FjyUPJf6nomzhs21wG515t8HHNDlujA2\nNjaUZZahjp0qsw0/o/5Zwh90uyIwxDFReInl+extOKzB/D27VYGitke/lcNoMcUUGS99t40K5hBG\nUnZJdympdZVKBbq/7zV2M6pUKrk6RFaAYdi2xm33VCqV6dudZTU+HHHYiqJj1z7OkiRJUgYmzpIk\nSVIGJs6SJElSBibOkiRJZTYn+vLmmUYWj/S61qXkqBqSJEllto7cFxROjk92oiYDzxZnSZIkKQMT\nZ0mSJCkDE2dJktQppwOrWP9W8EcBNwLXAR/rdqWkVtnHWZIkdcoZwOeAr6Tm7QccDDwbWAs8qQf1\nklpii7MkSeqUS4DVNfP+EfgokTQD3NvVGkltMHGWJEnd9HTghcAvgAngeT2tjZSDXTUkSVI3zQU2\nA/YCdgfOBrart+D4+PgTj8fGxhgbG+t87VRqExMTTExMdKz8SsdKlgbP1NTUVK/roJKrVCrQ/X2v\nsZtRpVIhz5aqAMOwbduM21HgAmDn5PkPgJOAnyTPbwH2BP5Us95Qxm2lUsk9JjPjtLTOMGzfove5\ndtWQJEnddB6wf/J4B2A+6yfNUl8ycZaU28jIkvy3dx1Z0utqS+q+s4BLiQT5DuBIYoi67Ygh6s4C\n3tSz2kk52cdZw2YZMSzSk4Ep4IvAfxAnud7K9NXdxwE/7EH9SmFycjXkOqENk5P2DJOG0GEN5r+x\nq7WQCmLirGGzFngPcBWwAPglsILIAj+dTJIkSesxcdawWZlMAA8Sd67aOnluk6gkSWrIPs4aZqPA\nrsRYohC3gL0aOA1Y3KM6SZKkPmWLs4bVAuAc4Gii5fkLwIeS1z4MfAp4S+1KjimqvDo9pqik/jWy\neITJByZ7XY3CtPJ5Fi5ayJr713SoRt3nqWkNo3nAd4mxRE+u8/ooM8ccrRrKMUXriXEx826LylCM\nGdpMG2OKNrqwdQnwTeCpwG3AocD9Nesauxk5jnN9jj/emm6OydyNcZxb/Ty9/B4dx1lqT4XoinED\nM5PmLVOPX0kMkyT1k+qFrTsRd1x7J7AjcCxxgesOwI+S55KkDrCrhobN3sDhwDXAlcm844khk3Yh\nWvJuBd7ek9pJjTW6sPVgYN9k/nJgApNnSeoIE2cNm59R/0zLD7pdEakNo8SFrZcBS4FVyfxVyXOp\nn5wOvAy4h/W7wB0DfALYHLivy/WScjNxltQlc6t9zTJZuHAz1qzx/2gdC4BziQtba6/SmaJB53Mv\nbFUeBV/UegbwOaKPftoy4EDg90W9kdRpJs6SuuQx8lxQ6J0G65pHJM1fBc5L5q0CtiC6cWxJtOqt\nJ504S83UHlydeOKJ7RR3CXGWpNangfcB/9VO4VI3eXGgJJVDowtbzweOSB4fwXRCLfWzVwB3Eteb\nSKVhi7MklUO9C1uPA04CzibGHb+NGI5O6mebEBdlH5iaV/cUk12MlFenx873XKiUXenHFC1Kq+M4\n51tnMMd9djzc/uY4zvUVELejTI+PvzNwEfDn5LVtgLuAPZjZ1aj0ces4zq29T5GK3ufa4ixJkrrp\nWmaO/nIr8FwcVUMlYB9nSZLUSWcBlxI36bkDOLLm9XI3K2uo2OIsSZI66bAmr2/XlVpIBbDFWZIk\nScrAxFmSJEnKwMRZkiRJysA+zpIkScNmzhNDtSkHE2dJkqRhs47WxosecnbVkCRJkjIwcZYkSZIy\nMHHWsFkGXAxcD1wHvCuZvwRYAdwMXAgs7kntJGnwnA6sIu4YWPUJ4EbgauDbwKIe1EvKzcRZw2Yt\n8B5gJ2Av4J3AjsCxROK8A/Cj5LkkqX1nAAfVzLuQ2A8/h2iwOK7blZJaYeKsYbMSuCp5/CDR4rE1\ncDCwPJm/HDik+1WT1AlLRkaoVCqZJnXEJcDqmnkriMvTAC4DtulqjaQWmThrmI0CuxI77aXEqUSS\nv0t7VCdJBVs9OckUZJrUE28Gvt/rSkhZOBydhtUC4FzgaGCy5rWG/0PHx8efeDw2NsbY2FhnaqeB\nMTExwcTERK+rIfWrDwCPAl+v96L7XOXV6X2u56U0jOYB3wV+AJyczLsJGCO6cmxJXED4zJr1pqam\nbJOC6qD5ebdF3nUqDOL2TroDdHvfO9SxW6lUMkde/ihlIOO0VgFxOwpcAOycmvd3wNuAFwEP11mn\n9HFbqVRaGyt5wNbp5fdY9D7XrhoaNhXgNOAGppNmgPOBI5LHRwDndblekjRMDgLeC7yC+kmz1JdM\nnDVs9gYOB/YDrkymg4CTgAOJq7v3T55Lktp3FnAp8AzgDqJP8+eILnMriP3w53tWOykH+zhr2PyM\nxgeMB3SzIpKmLRkZYfVk7eUG9W22cCH3rVnT4RqpQIfVmXd612shFcDEWZLUc9WRL7KoZEywJalo\ndtWQpHKod/e1ceBOZnY7kiR1iImzJJVDvbuvTQGfJsYj3xX4YbcrJUnDxMRZksqh3t3XwGFFJalr\nTJwlqdyOAq4mhllc3OO6SNJA8+JASSqvLwAfSh5/GPgU8JZ6C3oHNuXhHS+l+kycJam87kk9PpW4\nM1td6cRZaqb24OrEE0/sXWWkPmJXDUkqry1Tj1/JzBE3pH5Rb0SYJcTNT24GLsRuRioJE2dJKod6\nd1/7GHAN0cd5X+A9Paud1Fi9EWGOJRLnHYAfJc+lvmdXDUkqB+++prK6BBitmXcwcbAHsByYwORZ\nJWDiLEkqlblApeIofCW3lOi+QfJ3aQ/rImVm4ixJKpXHIPPtucGBrktgigZfqaPBKK9Ojwhj4ixJ\nkrptFbAFsJK4yPWeegs5Gozy6vSIMF4cKEmSuu184Ijk8RHAeT2si5SZibM05EZGllCpVHJNkpRD\n7YgwRwInAQcSw9HtnzyX+p5dNTRsTgdeRpwW3DmZNw68Fbg3eX4c8MOu16xHJidXk6/HKNhrVFIO\n9UaEATigq7WQCmCLs4ZNvfFEp4BPA7sm09AkzZIkKTsTZw2bS4DVdebbhCpJkmZl4iyFo4i7r52G\nt36VJEl12MdZgi8AH0oefxj4FPCWegs6pqjy6vSYopKk7jFxlmaOH3oqcEGjBR1TVHl1ekxRSVL3\n2FVDisH3q14JXNurikiSpP5li7OGzVnAvsDmxHiiJwBjwC7E6Bq3Am/vVeUkSVL/MnHWsKk3nujp\nXa+FJOk44HBgHXGm70jgkZ7WSGrCrhrSAPEugJJKYhR4G7AbcTOqDYDX9bJCUha2OEsDxLsASiqJ\nNcBaYBPg8eTvXT2tkZSBLc6SJKnb7iOG/rwd+ANwP3BRT2skZWCLsyRJ6rbtgXcTXTYeAL4FvAE4\nM72QY+crr06PnW/iLEmSuu15wKXAn5Ln3waezyyJs5RFp8fOt6uGpD41N/eFjiMjS3pdaUnZ3ATs\nBWxMXGhxAHBDT2skZWDiLKlPPUZc6Jh9iosjB9bpwCpm3qBnCbACuBm4EFjcg3pJrbga+ApwBXBN\nMu+LvauOlI2JsySVwxnAQTXzjiUS5x2AHyXPpbL4OLATMRzdEcQoG1JfM3GWpHK4BKhtUj8YWJ48\nXg4c0tUaSdKQMXGWpPJaSnTfIPm7tId1kaSB56gakjQYqp2963JYr86YC5nvwLnZwoXct2ZNZytU\nkE4P6SWVlYmzJJXXKmALYCWwJXBPowUd1qszqpewZlGZnOxkVQrV6SG9NETmZD+4fMIGxP0kc1i4\naCFr7u/8gamJsySV1/nERVUfS/6e19vqSFKNdcB4znXG868zOd6dA1P7OEtSOZxF3DDiGcAdwJHA\nScCBxHB0+yfPJUkdYouzJJXDYQ3mH9DVWkjSELPFWcPGm0hIUn9YDJwD3EjcNXCv3lZHas7EWcPG\nm0hIUn/4LPB9YEfg2UQCLfU1E2cNG28iIUm9twjYhzgLCDFAyQO9q46UjYmz5E0kJKnbtgXuJc4C\n/gr4ErBJT2skZeDFgdJM3kRChfJGElJdc4HdgH8CLgdOJrrJfTC9kPtc5dXpfW7OEamlgTAKXADs\nnDy/CRhj+iYSFwPPrLPe1NRU1lsd9EYMMp+3jv26TmvvUY7vqOv73lLEbuYbidC5SOp42X3+PTTS\ngbjdAvg50fIM8AIicX55apm+j9tmKpVKV8Ywdp1Yvl68FB27dtWQpm8iAd5EQpK6YSUxHvkOyfMD\ngOt7Vx0pG7tqaNicBewLbE7stD9I3DTibOAtwG3Aob2qnCQNkaOAM4H5wG+Jm/pIfc3EWcPGm0hI\nUn+4Gti915WQ8rCrhiRJkpSBibMkSZKUgYmzJEmSlIGJsyRJXTCXGBor67RkZKTXVZZUw4sDJUnq\ngsfIOUb05GSnqiKpRbY4S5IkSRmYOEuSJEkZmDhLkiRJGZg4S5KkXtkAuBK4oNcVkbIwcZYkSb1y\nNHAD+a6blHrGxFmSJPXCNsBLgVOBSo/rImXicHSSVH63AWuAx4G1wB49rY2UzWeA9wIOWK3SMHGW\npPKbAsaA+3pcDymrlwP3EP2bxxotND4+/sTjsbExxsYaLtpxI4tHmHzAsbX73cTEBBMTEx0r38RZ\nkgaDp7pVJs8HDia6amxEtDp/BXhTeqF04txrkw9MwnjOlfIur7bVHmCdeOKJhZZvH2dJKr8p4CLg\nCuBtPa6LlMXxwDJgW+B1wI+pSZqlfmSLsySV397A3cCTgBXATcAl6QW6fcp7ycgIq71ldGl1+nR3\nHY6qoVIwcZam3YYXWKmc7k7+3gt8h4jdholzN6yenMyVCdnPpL90+nR3jZ8kk9T37KohTateYLUr\nJs0qj02AhcnjTYEXA9f2rjqSNLhscZZmsuFLZbOUaGWG2KefCVzYu+pI0uAycZamVS+wehz4/4Ev\n9bY6Uia3Arv0uhKSNAxMnKVpfXeBlcqvBxdZSZI6xMRZmtZ3F1ip/Lp8kZUkqYO8OFAKXmAlSZJm\nZYuzFLzASpIkzcrEWQpeYCVJ3bWMuM32k4mLs78I/EdPayQ1YeIsSZJ6YS3wHuAqYAHwS+LC7Bt7\nWSlpNvZxliRJvbCSSJoBHiQS5q16Vx2pOVucJUnqQ3OBSiXbPZk2W7iQ+9as6WyFOmuUuGvrZXlW\n+pdj/4Xrbrwu1xttvOHGnPHFM1i8eHGu9SQwcZYkqS89RnT8zaIyOdnJqnTaAuAc4Gii5fkJzcbO\nP/3Lp7P6uatjXKSMNrpoI7Z5yjb87+T/tlxh9a9Oj51v4ixpgMzN3EJXtXDhZqxZc1+H6iOpiXnA\nucDXgPNqX8w0dv4zgEU53vDn85i8axIyFD2zMjmXV090eux8E2dJAyRPG12YnMyXaEsqTAU4DbgB\nOLnHdZEy8eJASZLUC3sDhwP7AVcm00E9rZHUhC3OkiSVXEkvJPwZNuCpZEycJUkquSG6kFDqKY/0\nJEmSpAxMnCVJkqQMTJwlSZKkDOzjLElq6ne/+x3XXHNNr6shST1l4ix12J/+9CdWrFiRe70XvvCF\nbLXVVh2okZTfh449lmu++12eMm9e02XXrFvXhRpJUveZOEsd9p3vfIejjvoI8+btkXmdRx65Brid\nRx/1lrCdl/9ug3Gzs7W51ij7HQqnHn+cox96iCMeeqjpstcDf9n5KklS19nHWZp2EHAT8Bvg/UUV\nOjU1xZw5L2Jy8ht1pn+oO//RR1+dJM1TOSeAiaKqnlJ0mUWX106Z1YG86k0XN5i/dpZ16k+Tk6tb\nrF8mHYnd2Ux0403aNGGZ/a57cXur5XSlnCLLKrJOBTJxlsIGwCnEjvxZwGHAjp1/24khLbPo8spU\nZuF6ErsTnX6DAkxYZj/rbtzeZjldKafIsooqp2AmzlLYA7iF+KmuBb4BvKKXFZIyMnZVRsatSsk+\nzlLYGrgj9fxOYM+iCl+79lfAh+q88pMG839a1Ftr8HU0dtPOA36fPG4UuQD3dOLNNWiKi9tfABs2\nWeY2nmiGf2TNIy29jSRp2quBL6WeHw58rmaZW8jf6djJqXa6hWIZu07dmIxbp7JOhcauLc5SuAtY\nlnq+jGgBSXta96ojZWbsqoyMW0kqsbnAb4FRYD5wFV25OFBqm7GrMjJuJankXgL8mjitc1yP6yLl\nYeyqjIxbSZIkSZI0+JYAK4CbgQuBxQ2WazZw/zHAuqS8dsv8BHAj8Dtgkji92ehmAf+RrH81sGuG\nsrPcgKBemcuIu4NcD1wHvKuAMqs2AK4ELiiozMXAOcQ2vAHYq83yjiM+97XA15m+nr9Zmc8Efg48\nTMRHWhE3gjgdWJXUq5HZtnsrZb4hKesa4L+BZxdQR4DdiTvTvKqAOgKMETF1HdmGOG5W5ubAD4nu\nBdcBf5ehzNl+M2lZv6Ms5eX9frLWEfJ9R800i/+s27uo30BRcV9UvBcV40XFdZGxXFQc9yp2paH3\nceB9yeP3AyfVWWYD4tTiKHHv5dq+ecuInc+tRNLcbpkHEv0BbwG+kJRXrz/gS4HvJ4/3JAZpmq3s\nZp9jtjK3AHZJHi8gTre2W2bVPwNnAuc3qX/WMpcDb04ezwU2a6O8UeIApposfxM4ImMdnwQ8D/g3\nZibOWdbNYh/in1Ojf4rNtnsrZf4VsCh5fFCGMpuVB7E9fgx8lxh5od06Lib+mW6TPN+8gDLHgY+m\nyvsTzS+2b/SbScvzHWUpL+/3k6VMyP8dzSZL/I+TbXsX9RsoKu6LiveiYryouC4ylouK467FrjdA\nkWY6mEi0SP4eUmeZZgP3f5rpRLmIMlcQR8i3EC3WW9V5z9r3uYzYmW4xS9lZbkBQr8ylwEriHxzA\ng0Rr7lZtlgmx438pcCpQybBtmpW5iPhncXry2mNEy2+r5a1J1tmE+IeyCTE6QJY63gtckbyeVtSN\nIC4BZruv92zbvdUyfw48kCpzm1mWzVIewFHEGYJ7myyXtczXA+cyPWLDHwso825gJHk8QiQYjzUp\ns9FvJi3Pd5SlvLzfT5YyIf93NJss8Z91exf1Gygq7ouK96JivKi4LjKWi4rjrsWuibM001LiVBbJ\n33o/9HoD92+dPH5F8vyaAstMv/5m4ii+9vXZytiqwfxm79mozNqd1ijRinFZG2VWl/kM8F6im0uW\n5ZvVc1tiB3kG8Cti3Nht26jjfcCngNuBPwD3AxdlrGMj7aybR5bvsh1vYbqFqVVbE7+hLyTPp9os\nD+DpxJmfi4kDlzcWUOaXgJ2IGLgaODrn+qNM/2bSWv2OGpWXlvf7aVRm0d9Rlvhvd3vP9l7t/gba\nifuitmVRMd7Kdh6luFhuVFZalu09W53a3t6O46xhtIJoia31gZrn1cHTa9Wb9wpgP2A7oovGtcSP\n9CVtlFn7+i5EwvZ14mYB9VQazG9UZha1ZabXW0AcvR9NHOW3WmYFeDlx07krif56RdRzLrAb8E/A\n5cDJRIv/AzRXb1tuD7yb2DE/AHyL6H/3cMY61lNEcpjVbN9lO/YjDur2brOck4FjiXpVyBfPjcwj\nYuBFxBmCnxOnen/TRpnHE61bY0RMrACeQ1yD0Eztb6ZW3u+oWXmQ//uZrcyiv6MsMdjO9q5V5G+g\n3bgvalsWFeN5t3ORsVxUHHc8dk2cNYwOnOW1VURSvRLYkvp3D643cP9/At8DfkTsuCB+lP9OnDZr\npcz0zQB2IW4GsE+D1+uVsU2yzLwGZWe5AUG9Mu9KHs8jTg9+jbgbc5bPMVuZryZO770U2Ig4VfgV\n4PNtlFlJlr08mX8O0cd8JLVsnvLGgEuJU5gA3waeD3w1Qx0bybLNijDbd9mOZxMtVQfR/LR0M88l\nTtVD9LF8CXH6/vyGazR3B/EbfCiZfkokA+0kzs8nftsQF+veCjyDaO2bTb3fTFre76hZeZD/+2lW\nZtHfUZb4b3V7N3uvdn4DRcR9UduyqBjPs52LjOWi4rjbsSuJuPCuelX3sdS/kC/rwP3piwPbKfMg\n4gYkvocAAAf4SURBVMKPW5u8Z/pijL2YvoCiUdlZPkejMitEUvuZHJ+jWZlp+zI9qka7Zf4U2CF5\nPE58H62WtwtxxfbGxDZYDrwzYx2rxpl5cWCRN4IYJduFUY22e94yn0L0T90rY1nNyks7g+xXvc9W\n5jOJ7jQbEAe11wLParPMTwMnJI+XEonekiblNfrNpOX5jrKUl/f7yVJmWp7vqJEs8Z9ne49SzG9g\ntnLybNfZyklrti1nKydPjM9WTtbtXGQsFxXHvYhdScRO4iLWHzpuK6JFuSrLwP2/Y3o4unbK/A3w\n+2T+w0QXgerrb0+mqlOS5a4mTt01K7ve/CxlvoDoh3wV0bXiSiLBb6fMtH2Z2QrQTpnPIVqcryZa\niBe1Wd77mB6ObjnRypGljlsQLUMPEC0mtxOnFRutm9dZRN/ER5nuD593u+ct81Si9b0aA/9TQB2r\nsv5jy1LmvzD9nc02TFXWMjcnDuyuTsp8fYYy6/1mXlKnrlm/oyzl5f1+staxqqjko9lvJ+v2Luo3\nUFTcFxXvRcV4UXFdZCwXFce9il1JkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJ\nyq7S6wpIkpTTu4CdiDugnQC8lbhr2N7Enc/+i7hb2Dpgd+Ac4q5pktSWub2ugKSeMwlRmewI3Ad8\nDrgGeAT4QPL3ZuBrwPZELD8G/A1wOsasOs99qSQNuB2Bw4G/JHbmnwQ2TF47GFgD/DvTB9l/A9zb\n5TpKaUcAmxFx+0dgUeq1NwL3A4tT896AMavOc186JOb0ugKSemoP4HvALkQr3oeJljuIhGQd8Ami\n5Q5gpNsVlGosB1YD+wErgAdSr+0D/IRInqv2AyYalHUI8FHgS8BhxG9gNGM93gEc3+C13YDLgXdn\nLEvl575UkobIaax/yvCLxKnFtFOBb3Xg/f8S+G/g6A6UrcH0O+K0d9pvmBlD84lW6Vc3KOMmYDvg\noGTZ3+R4/6XAprO8fj6RQGu4dGtfOtuBW6e8nnwHl5I0sIpIQtr1NSKBlpp5KtGCt0Nq3jbJvJ1T\n814DrCROj78KeFrqtW2BS1LP96G4/qYbALfjBfjDqFv70mYHbp2wIfkOLiVpIBWRhBThVwWXp8F1\nJHBXzbw3AvfUzPss8Hkigf1Yav5ewNlEt47jiG6L/wr8Q/L6RkTy8w2m+6SeSlz4tXHy2jnAvFSZ\nC4lT8a8j+rJ+r6VPpjLrl31ppxR5cFla9nGWtD9wNzEiQdV+RIvItal5+wDfBh4H9iRaH/4BuAL4\nT+Ay4G+TZTckrio/kkhI9knmbwSMJ/OPBZYl87cnWjJGiKTkT8ALivl4GkDbAWfWmfe1mnnfJGLs\n34nT5VW/AK5Kyvgokdi8kOkW6JcBXwaeQ7QOQvwmfpO8dkbyWjpxPpMYOeEbSXk/beWDqdRa3ZdC\ntEh/Engl8Cbg60SMvQH4IHFRLMSB298z88Dt6UR/+o8T/fbfCHw1WRbiYtoTia4W/w5sThwENtp/\nLyAONF9LjFhzcjJ/X+JgU5KG2oeJHW7aCcCna+Y9n0gMPkK0kOxP7GDvAf4C2ILpJONc4EXJ402I\nYcMAvkP8I6k+fnHy+EjgFKZPXS7DA3t11pnEMGEQCcjtqdc2AXYFvps8H2U6qd6U+C18J7X87sAN\nqefnEq3aGi6t7kvfRHS9WEXsSyFado9JHu+cLA/RWl3tj79JMq+6/t3AkmTet4gLFTcgrh+ptoK/\ng+imNNv++3ziQBLg34B3Jo8vJBJuSVKLDiRa5tL2JloxqrYhWpD3JkYaSM+v+jLw2+SvV5urG65k\neti65xPJwsap1/+VaNmDaKn7CNOx+XniFHv1+THJvKrfEcn4wsJrrUG0KdHSfEpq3lVMJ8GvJ7oB\nVZetPXCrrv/Z1Lwbifj7G6Lh4lDgbcDzUsvU23/XHgR+mziIrD24HFr/r727Z40qiOIw/igBsVlI\nExAEK0URRRQ0KlEMQiqxU9EUgpX2fgOLfAdrIXai4FtjoSiIioqFWhhEiaQRhMRgxJfiP5e9d3fB\nwC4E2ecHgTs3c+EWu8OZM2dmzehI6scEzQ1WkEH9Qa09SZb3DtJc5vtcu95PzkH9DVwgGZDRQb+s\nVIyQI8KqY+sOkKzc6VqfTcDHcn2cTPqqFZEpUsN8qrS/0T6T9yjwCdhNMnvSvyzRHEtHSaD6tbRP\nkM9bq/SdJqUYrY7nH5X2VlJD/YOc7HKP1PRfpZnU6DV+H6F5fOMe4BUZ11/SnFwOJQNnSf04TPfA\nuwAsl+sNJMtxudz/Xuu3hQTLm8m5vCtk+XGFlHDU+0qDtA14V2s/I0HuQu3eNVLPfBZ4QzJxH8r/\nnpCa06qUYxYYI4F3C5gn343XSKtTH0snyEQOUmoxToLiM+Ve58Ster4KnE8CN0hZxwvgZ63fXtpl\nG73G7/ok8BjZhLsT2Ef35HIo+ZPbkvqxkeZGGMgGrSvAebLp7xIpw5gDZkjA8QdYJMuAU8Dd8uwd\nEmjfp/3jAdKg7aKZVXtIdwDxuPz1Mt3RXgQu1tq3kFZvPQlu50t7B3C7XP8CnpJx8Xq5V03cbnY8\n/6W035MVvjnyWTxE6qDXlT7VCUa9xu9Z2qfDLJf+k2TF5RxOBiVJGhrjZNl6Bti+xu8i6T9kxlmS\nNCyWgOdks+rbNX4XSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZI0TP4Cx61S\nWejT998AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1be678b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(1, [12, 6], dpi=300)\n",
    "fig.subplots_adjust(hspace=.3)\n",
    "ax1 = fig.add_subplot(2,4,1)\n",
    "counts, bins, ignored = ax1.hist(data_obs, normed=False, color='black')\n",
    "ax1.set_xlabel(r'$\\sigma_{obs}$', fontdict={'size':18})\n",
    "ax1.set_xscale('log')\n",
    "\n",
    "ax2 = fig.add_subplot(2,4,2)\n",
    "counts, bins, ignored = ax2.hist(sigma_r, normed=False, color='blue')\n",
    "ax2.set_xlabel(r'$\\sigma_{rock}$', fontdict={'size':18})\n",
    "\n",
    "ax3 = fig.add_subplot(2,4,3)\n",
    "counts, bins, ignored = ax3.hist(sigma_f, normed=False, color='red')\n",
    "ax3.set_xlabel(r'$\\sigma_{fluid}$', fontdict={'size':18})\n",
    "\n",
    "ax4 = fig.add_subplot(2,4,4)\n",
    "counts, bins, ignored = ax4.hist(sigma_m, normed=False, color='green')\n",
    "ax4.set_xlabel(r'$\\sigma_{mineral}$', fontdict={'size':18})\n",
    "\n",
    "ax6 = fig.add_subplot(2,4,6)\n",
    "counts, bins, ignored = ax6.hist(m_r, normed=False, color='blue')\n",
    "ax6.set_xlabel(r'$m_{rock}$', fontdict={'size':18})\n",
    "\n",
    "ax7 = fig.add_subplot(2,4,7)\n",
    "counts, bins, ignored = ax7.hist(m_f, normed=False, color='red')\n",
    "ax7.set_xlabel(r'$m_{fluid}$', fontdict={'size':18})\n",
    "\n",
    "ax8 = fig.add_subplot(2,4,8)\n",
    "counts, bins, ignored = ax8.hist(m_m, normed=False, color='green')\n",
    "ax8.set_xlabel(r'$m_{mineral}$', fontdict={'size':18})\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 286,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Python27\\lib\\site-packages\\ipykernel\\__main__.py:17: RuntimeWarning: invalid value encountered in power\n"
     ]
    }
   ],
   "source": [
    "\n",
    "    \n",
    "# prior information of unknown model parameters\n",
    "phi_r = pymc.Normal('phi_r', .7, .2)\n",
    "phi_f = pymc.Normal('phi_f', .5, .3)\n",
    "phi_m = pymc.Normal('phi_m', .5, .3)\n",
    "\n",
    "#m_r = pymc.Normal('m_r', .03, .015)\n",
    "#m_f = pymc.Normal('m_f', 2, .3)\n",
    "#m_m = pymc.Normal('m_m', 3.5, .3)\n",
    "\n",
    "# expected outcome is governed by Archies equation of Glove 2010\n",
    "@pymc.deterministic\n",
    "def archies(s_r=sigma_r, s_f=sigma_f, s_m=sigma_m, \n",
    "            p_r=phi_r, p_f=phi_f, p_m=phi_m, \n",
    "            m_r=m_r, m_f=m_f, m_m=m_m):\n",
    "    archie_return = (s_r*p_r**m_r)+(s_f*p_f**m_f)+(s_m*p_m**m_m)\n",
    "    return archie_return\n",
    "\n",
    "\n",
    "# likelihood of the observations is a normal distribution\n",
    "# and is the data itself\n",
    "sigma_obs = pymc.Normal('sigma_obs', mu=archies, value=data_obs, observed=True)\n",
    "\n",
    "archie_model = pymc.Model([archies, sigma_r, sigma_f, sigma_m, \n",
    "                           phi_r, phi_f, phi_m, \n",
    "                           m_r, m_f, m_m, \n",
    "                           sigma_obs, \n",
    "                           data_obs])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 287,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " [-----------------100%-----------------] 10000 of 10000 complete in 3.0 sec"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Python27\\lib\\site-packages\\ipykernel\\__main__.py:17: RuntimeWarning: invalid value encountered in power\n"
     ]
    }
   ],
   "source": [
    "model_fit = pymc.MCMC(archie_model)\n",
    "model_fit.sample(iter=10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 288,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "             Rock   Fluid    Mineral\n",
      "Mean Percent   0.61   0.12   0.52\n",
      "Median Percent 0.56   0.05   0.60\n",
      "Mode Percent   0.54   0.05   0.63\n",
      "Min Percent    0.18   0.00   0.00\n",
      "Max Percent    1.00   0.85   1.00\n",
      "STD            0.24   0.22   0.30\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAGNCAYAAAAM1oc0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+YXNV52PHvgJCNzA7ylkQISbVUVzTIxQETBAl2GCeY\nyHWDcJKCcE2UoqSkymNct02MaBtGdoJJ/NgxiR9oa/NDODxyqN1Q2cEE7DAOpAHFlB8CmSBRyUEK\nCCcW0tgOiQTbP85Z9uzu/NqZOzN3Zr6f57nP3Dn3zsy7Wr173zn33HNBkiRJkiRJkiRJkiRJkiRJ\nkiRJkiRJkiRJkiRJktRDrwceBh4DdgIfi+3jwH3AM8C9wMLkNZuAXcDTwIU9i1QaHfXysgzsAx6N\ny7uT19TLy7OAHXHbDd0MWtIsWeaypB5bEB/nAQ8Bbwd+G/i12P5h4Pq4voqQ6McBy4HdwDG9ClQa\nIbXy8lrgP9TYt1ZeFuK27cDquH43sKY74Uqqo9Nc9hgrdUErifX9+DgfOBY4CFwEbIntW4CL4/pa\nYCtwBNhLSN7Jg6+k7NTKS5gqfFO18vIcYDEwRiiSAW5nKpcl9UanuewxVuqCVgrkYwjfWA8A9wNP\nAYvic+Ljorh+CuG00KR9wJJMIpWUqpWXAB8AHgduZmroU728nNm+H/NV6rUscllSxlopkF8FzgCW\nAj8OvHPG9om41DNr25vf/ObJ17i4uIRlN3MzMy9LwE3Aitj+PPCJOb5nTeari8usZa752kinuTyR\nPjFfXVxmLW3l61zGLh0C/ohwUc8B4OTYvhh4Ma7vB5Ylr1ka26Z59tlnmZiYyOVy7bXX9j2GQY3P\n2NpfgDfPMXdn5uWPxDyc/IPwWaZOvdbKy32xfekg52vef7fGNpzxdZCvjbSby9Ny1nwdztjyHl+e\nY2s3X5sVyCcxdWrneOBdhCtqtwHrY/t64K64vg1YRxhLtQJYydT4RknZqJeXJyf7vJcwOwXUz8sX\ngMOE8cgF4HKmcllS92WVy5IyNq/J9sWEi/COicvngK8REvhOYAPhQoFL4v47Y/tO4CiwkfANWFJ2\n6uXl7YRTshPAHuDKuH+jvNwI3EY4ON8N3NOLH0ASkG0uS8pQswJ5B/C2Gu3fAS6o85rr4jKQSqVS\nv0NoKM/xGVvP1MvLn2/wmnp5+QhwehZB9Uuef7fG1r68x5eRLHN5IOT595rn2CDf8eU5tnbVmkam\nFybiuBBJQKFQgP7lYzPmq5QwX6XB0W6+OsG4JEmSlLBAliRJkhIWyJIkSVLCAlmSJElKWCBLkiRJ\nCQtkSZIkKWGBLEmSJCUskCVJkqSEBbIkSZKUsECWJEmSEhbIolgcp1Ao1F2KxfF+hyhJ0sgqLiw2\nPE6/drxeWOx3qEOjX/eS917xORLuU97o91HA31d3tXuv+B4xX6WE+apeKxQKUG5hxzIer2doN1/t\nQZYkSZISFsiSJElSYl6/A5CUf29905vqbivMm8d/v+MOzj333B5GJElS91ggS2rq9//qr+pu+/CC\nBezZs8cCWZI0NCyQJTX11gbbFh57bM/ikCSpFxyDLEmSJCUskCVJkqSEBbIkSZKUsECWJEmSEhbI\nGnnealuSJKWcxUIjr1o9SKNbbVereb2jrCRJ6gZ7kCV17MorrmjYCz9eLPY7REmSWmYPsqSOVV9+\nuUEfPBSq1Z7FIklSp+xBliRJkhIWyJIkSVLCAlmSJElKWCBLkiRJCQtkSZIkKWGBLEmSJCUskCVJ\nkqSEBbIkSZKUsECWJEmSEhbIkiRJUsICWZIkSUpYIEuSJEkJC2RJkiQpYYEsDZ7XAw8DjwE7gY/F\n9nHgPuAZ4F5gYfKaTcAu4GngwqT9LGBH3HZDV6OWNFOWuSwpQxbI0uB5GXgncAbw1rj+duBqwkH1\nVOBr8TnAKuDS+LgGuBEoxG03ARuAlXFZ05OfQBJkk8sex6UuMLGkwfT9+DgfOBY4CFwEbIntW4CL\n4/paYCtwBNgL7AbOARYDY8D2uN/tyWsk9Uanuby6V4FKo8QCWRpMxxBOyx4A7geeAhbF58THRXH9\nFGBf8tp9wJIa7ftju6TeySKXJWVsXr8DkNSWVwmnZU8E/phwajY1EZdMlJP1UlykUVGpVKhUKt16\n+05zeda2crn82nqpVKJUKnUaozQwssrXQvNdumJiYiKzY7c6VCgUaPz3t8Aw/77y8POHGNrOx/8K\n/B3wi4Ta9QXC8In7gR9iavzi9fHxHuBa4Ftxn9Ni+2XA+cAvz3j/hj/9ZWNjfL5abfIvyFD/H9Jo\n6TBfG2k3lx9O3sPj6xAqFArTeyrqKfu3dqZ289UhFtLgOYmpq9qPB94FPApsA9bH9vXAXXF9G7CO\nMMZxBeFivO2Eg+9hwnjkAnB58hpJ3ZdVLkvKmEMspMGzmHDhzjFx+RzhSvdHgTsJs1LsBS6J+++M\n7TuBo8BGprrMNwK3EQ7OdxN6pCT1Rpa5LClDzbqclxGubP9BQhL+D+B3CR39vwh8O+53DfCVuL4J\nuAJ4BbiKMIfjTJ4CypE8DDHopzz8/F08ZZsFh1hIidznq7k2dBxi0b5287VZD/IR4EOEK2xPAB4h\nzM04AXwyLql0jsYlwFcJ8zi+OtfAJEmSpH5oNgb5BUJxDPBd4JtMTSlTqxp3jkZJkiQNtLlcpLcc\nOBN4KD7/APA4cDNTFxk4R6MkSZIGWqsX6Z0AfAH4IKEn+SbgI3HbR4FPEC4mqKXmYBjnadQo6/K8\nqpIkqQOtDFo+Dvgy4SK8T9XYvhz4EnA6rc3RCF5EkCt5uEitn/Lw8+f+op8GG71IT6Mm9/lqrg0d\nL9JrX7fmQS4QhlDsZHpxvDhZfy+wI647R6MkSZIGWrMhFucB7weeIMzLCGFKt8sIt8acAPYAV8Zt\nztEoSZKkgdasQH6Q2r3MX6nRNum6uEiSJEkDx1tNS5IkSQkLZEmSJClhgSxJkiQlLJAlSZKkhAWy\nJEmSlLBAliRJUluKC4sUCoWmS3Fhsd+hzkmrt5qWJEmSpqkeqrZ0l79qudr1WLJkD7IkSZKUsECW\nJEmSEhbIkiRJUsICWZIkSUpYIEuSJEkJC2RJkiQpYYEsSZIkJSyQJUmSpIQFsiRJkpSwQJYkSZIS\nFsiSJElSwgJZkiRJSlggS5IkSQkLZEmSJClhgSxJkiQlLJAlSZKkhAWyJEmSlLBAliRJkhIWyJIk\nSVLCAlmSJElKWCBLkiRJCQtkSZIkKWGBLEmSJCUskCVJkqSEBbI0eJYB9wNPAU8CV8X2MrAPeDQu\n705eswnYBTwNXJi0nwXsiNtu6GbQkmbJMpclZWhevwOQNGdHgA8BjwEnAI8A9wETwCfjkloFXBof\nlwBfBVbG/W8CNgDbgbuBNcA9Xf8JJEE2uXwq8GqP4pVGhj3I0uB5gXBABfgu8E3CwRKgUGP/tcBW\nwsF4L7AbOAdYDIwRimOA24GLuxKxpFqyyOXV3Q1RGk0WyNJgWw6cCTwUn38AeBy4GVgY204hnK6d\ntI9wEJ7Zvp+pg7Ok3lpO+7ksKWMOsZAG1wnAF4APEnqfbgI+Erd9FPgEYfhEx8rJeiku0qioVCpU\nKpVufkQnuTwxs6FcLr+2XiqVKJVK2UUq5VxW+VrrFE4vTExMzMpp9UmhUKDG39h0D4b595WHnz/E\nMKd8PA74MvAV4FM1ti8HvgScDlwd266Pj/cA1wLfIlwgdFpsvww4H/jlGe/V8Ke/bGyMz1erTf4F\nGer/QxotbeRrI1nk8sPJ/h5fh1ChUJjeU1FPufd/a/McG7Sfrw6xkAZPgXDadSfTD6iLk/X3Eman\nANgGrAPmAysIF+htJ4x/PEwYj1wALgfu6mbgkqbJKpclZcwhFtLgOQ94P/AEYQoogGsIPcBnELrD\n9wBXxm07gTvj41FgI1Nd5huB24DjCbNYOIOF1DtZ5vKcFBcWqR6qNt1v7MQxDr90uJ2PkAaaQyyU\niyEG/ZSHnz/jU7ZZc4iFlMh9vraQa3k/La7p8vz7ynNs4BALSZIkKRMWyJIkSVLCAlmSJElKWCBL\nkiRJCQtkSZIkKWGBLEmSJCUskCVJkqREswJ5GeFWtE8BTwJXxfZx4D7gGeBeYGHymk3ALuBp4MIs\ng5UkSZK6rVmBfAT4EPAW4FzgV4DTCPeDvw84FfgaU/eHXwVcGh/XADe28BmSJElSbjQrXl8AHovr\n3wW+CSwBLgK2xPYtwMVxfS2wlVBY7wV2A6uzC1eSJEnqrrn07i4HzgQeBhYBB2L7gfgc4BRgX/Ka\nfYSCWpIkSRoI81rc7wTgi8AHgeqMbRNxqafmtnK5/Np6qVSiVCq1GIo0+CqVCpVKpd9hSJKkGlop\nkI8jFMefA+6KbQeAkwlDMBYDL8b2/YQL+yYtjW2zpAWyNGpmfincvHlz/4KRJEnTNBtiUQBuBnYC\nn0ratwHr4/p6pgrnbcA6YD6wAlgJbM8qWEmSJKnbmvUgnwe8H3gCeDS2bQKuB+4ENhAuxrskbtsZ\n23cCR4GNNB5+IUmSJOVKswL5Qer3Ml9Qp/26uEiSJEkDxzmKJUmSpIQFsiRJkpSwQJYkSZISFsiS\nJElSwgJZkiRJSlggS5IkSQkLZEmSJClhgSxJkiQlLJAlSZKkhAWyJEmSlLBAliRJkhIWyJIkSVLC\nAlmSJElKWCBLkiRJCQtkSZIkKWGBLEmS1KLiwiKFQqHpUlxY7Heo6sC8fgeg5orFcarVg3W3j429\nkcOHv9PDiCRJGk3VQ1Uot7Bfudr1WNQ9FsgDIBTHEw22F3oXjCRJ0pBziIUkSZKUsECWJEmSEhbI\nkiRJUsICWZIkSUpYIEuSJEkJC2Rp8CwD7geeAp4Erort48B9wDPAvcDC5DWbgF3A08CFSftZwI64\n7YauRi1ppixzWVKGLJClwXME+BDwFuBc4FeA04CrCQfVU4GvxecAq4BL4+Ma4EZgcm7Am4ANwMq4\nrOnJTyAJssllj+NSF5hY0uB5AXgsrn8X+CawBLgI2BLbtwAXx/W1wFbCwXgvsBs4B1gMjAHb4363\nJ6+R1H1Z5PLqHsUqjRQLZGmwLQfOBB4GFgEHYvuB+BzgFGBf8pp9hIPwzPb9sV1S7y2n/VyWlDHv\npCcNrhOALwIfBGbe03SCRrdfnKNysl6KizQqKpUKlUqlmx/RSS7P2lYul19bL5VKlEqljgOUeqW4\nsBhu593E2IljHH7p8Kz2rPLVAlkaTMcRDqifA+6KbQeAkwmnbRcDL8b2/YSLgSYtJfQ87Y/rafv+\nWh9WzihoaRDNLDI3b96c5dt3msuzcjYtkKVBUz1UbemgUy3XLqKzyleHWEiDpwDcDOwEPpW0bwPW\nx/X1TB1stwHrgPnACsLFeNsJB9/DhPHIBeDy5DWSui+rXJaUMXuQpcFzHvB+4Ang0di2CbgeuJMw\nK8Ve4JK4bWds3wkcBTYydVp2I3AbcDxwN3BPt4OX9Josc1lShiyQpcHzIPXP/lxQp/26uMz0CHB6\nFkFJmrMsc1lShhxiIUmSJCUskCVJkqSEBbIkSZKUsECWJEmSEhbIkiRJUsICWZIkSUpYIEuSJEkJ\nC2RJkiQpYYEsSZIkJSyQJUmSpIQFsiRJkpSwQJYkSZISFsiSJElSwgJZkiRJSlggS5IkSYlWCuRb\ngAPAjqStDOwDHo3Lu5Ntm4BdwNPAhZlEKUmSJPVIKwXyrcCaGW0TwCeBM+Pyldi+Crg0Pq4Bbmzx\nMyRJkqRcaKV4fQA4WKO9UKNtLbAVOALsBXYDq9sNTpIkSeq1Tnp3PwA8DtwMLIxtpxCGXkzaByzp\n4DMkSZKknprX5utuAj4S1z8KfALYUGffiVqN5XL5tfVSqUSpVGozFGnwVCoVKpVKv8OQJEk1tFsg\nv5isfxb4UlzfDyxLti2NbbOkBbI0amZ+Kdy8eXP/gpEkSdO0O8RicbL+XqZmuNgGrAPmAyuAlcD2\ntqOTJEmSeqyVHuStwPnAScBzwLVACTiDMHxiD3Bl3HcncGd8PApspM4QC0mSJCmPWimQL6vRdkuD\n/a+LiyRJkjRwnKNYkiRJSlggS5IkSQkLZEmSJClhgSxJkiQlLJAlSZKkhAWyJEmSlGj3TnqSJCmH\nPv7xj/c7BGngWSBLkjRErvniNXBc/e0TL3v/LqkZC+QcKBbHqVYP9jsMSdIQOHr+UVjQYIeDwDd6\nFU17iguLVA9VG+4zduIYh1863KOINGoskHMgFMeNvtEXehWKJEl9Vz1UhXKTfcqNC2ipE16kJ0mS\nJCUskCVJkqSEBbLUoWJxnEKh0HApFsf7HaYkSWqRY5ClDjUfQw7VquPIJUkaFPYgS5IkSQkLZGkw\n3QIcAHYkbWVgH/BoXN6dbNsE7AKeBi5M2s+K77ELuKF74UqqIas8lpQxC2RpMN0KrJnRNgF8Ejgz\nLl+J7auAS+PjGuBGpuYOvAnYAKyMy8z3lNQ9neaxx3CpS0wuaTA9QJjuf6Zag53XAluBI8BeYDdw\nDrAYGAO2x/1uBy7OOlBJdXWax6u7Fpk04iyQpeHyAeBx4GZgYWw7hXDKdtI+YEmN9v2xXVJ/zSWP\nJXWBs1hIw+Mm4CNx/aPAJwjDJzpWTtZLcZFGRaVSoVKp9Orj5pLHtafPeRA4Lq4vB1ZkF5yUd1nl\nqwWyNDxeTNY/C3wpru8HliXblhJ6n/bH9bR9f603LmcWojR4SqUSpVLpteebN2/u5sfNJY9r5itv\nBxZ0IzQp/7LKV4dYSMNjcbL+XqaujN8GrAPmE/qSVhLGHb8AHCaMRy4AlwN39SpYSTXNNY8ldYE9\nyNJg2gqcD5wEPAdcSxj5cAbhtOse4Mq4707gzvh4FNjI1KnZjcBtwPHA3cA9vQheEpBdHg+M4sIi\n1UPVfochNWWBLA2my2q03dJg/+viMtMjwOmZRCRprrLK44FRPVRtbcxWK/tIXeQQC0mSJClhgSxJ\nkiQlLJAlSZKkhAWyJEmSlLBAliRJkhIWyJIkSVLCAlmSJElKWCBLkiRJCQtkSZIkKWGBLEmSNEKK\nC4sUCoWmS3Fhsd+h9o23mpYkSRohrd7yu1qudj2WvLIHWZIkSUpYIEuSJEkJC2RJkiQpYYEsSZIk\nJSyQJUmSpIQFsiRJkpSwQJYkSZISFsiSJElSwgJZkiRJSlggS5IkSQkLZEmSJCkxr4V9bgHeA7wI\nnB7bxoE/AN4E7AUuAV6K2zYBVwCvAFcB92YXbraq1Sq33norr7zySt19Tj31VN7znvf0MCpJkiT1\nUysF8q3A7wG3J21XA/cBvw18OD6/GlgFXBoflwBfBU4FXs0u5Ow89thj/NqvXc/ExKU1t09M/A2L\nF9/Mt75lgSxJkjQqWimQHwCWz2i7CDg/rm8BKoQCeS2wFThC6FneDawGHuo40i55/evfzKFDv1Nn\n65NMTKzraTySJEnqr3bHIC8CDsT1A/E5wCnAvmS/fYSeZEmSJGkgtNKD3MxEXBptn6VcLr+2XiqV\nKJVKGYQiDYZKpUKlUul3GJIkqYZ2C+QDwMnAC8BiwgV8APuBZcl+S2PbLGmBLI2amV8KN2/e3L9g\nJEnSNO0OsdgGrI/r64G7kvZ1wHxgBbAS2N5JgJIkSVIvtdKDvJVwQd5JwHPArwPXA3cCG5ia5g1g\nZ2zfCRwFNtJ4+IUkSZKUK60UyJfVab+gTvt1cZEkSZIGjnfSU+4Vi+MUCoW6S7E43u8QJUnSEMli\nFgupq6rVgzQaqVOtFnoXjCRJ6r9joFBofPwfO3Gs7be3QJYkSdJgeRUoN96lWq62/fYOsZAkSZIS\nFsiSJElSwgJZkiRJSlggS5IkSQkLZGkw3UK45fuOpG0cuA94BrgXWJhs2wTsAp4GLkzaz4rvsQu4\noYvxSpotqzyWlDELZGkw3QqsmdF2NeHAeirwtfgcYBVwaXxcA9wITM6NcxPhjpgr4zLzPSV1T6d5\n7DFc6hKTS2pqXsMblfTJA8DBGW0XAVvi+hbg4ri+lnDL+COEW8PvBs4BFgNjwPa43+3JayR1X6d5\nvLr7IUqjyQJZauoo4UYl9ZbcWEQ4XUt8XBTXTwH2JfvtA5bUaN8f2yX1z1zzWFIXeKMQaThlWr2X\nk/VSXKRRUalUqFQq/fjoZnlce9uDwHFxfTmwItOYpHzbQzjH0iELZGl4HABOBl4gDJ94MbbvB5Yl\n+y0l9D7tj+tp+/5ab1zOOFBpkJRKJUql0mvPN2/e3M2Pm0se18xX3g4s6F6AUq6tYOpL4dfbfxuH\nWEjDYxuwPq6vB+5K2tcB8wl/NlYSxh2/ABwmjEcuAJcnr5HUH3PNY0ldYA+y1BPzsr6gbytwPnAS\n8Bzw68D1wJ2EWSn2ApfEfXfG9p2EAdUbmTo1uxG4DTgeuBu4J8sgJTWUVR5LypgFstQTkxf61TPn\n4vmyOu0X1Gm/Li4zPQKcPtcPl5SJrPJY6o5j6OdsTX1lgSxJkqTZXqX5RSjNtg8oxyBLkiRJCXuQ\nJUmShsEID4nImgVyE/v27Wr4n21s7I0cPvydHkYkSZJUQytDImhxnxFngdzExMQ/0OjiqmrVb2qS\nJEnDxDHIkiRJUsICWZIkSUpYII+AYnGcQqFQdxl2o/7zS5KkuXEM8gioVg+S8U0qBsqo//ySJGlu\n7EGWJEkdKS4sNjxT5xk7DRp7kCVJUkeqh6pOL6ahYg+yuq7ZGOBicbzfIUqSJL3GHmR1XbMxwM4l\nLUmS8sQeZLVgnj3AkiRpZNiDrBYcxR5gSZI0KuxBliRJkhIWyJIkSVLCAlmSlLnxYvN5cceLxX6H\nKUk1OQZZkpS5g9Vqw/tXAhSq1Z7EIklzZQ+yJEmSlLBAliRJytoxNB1mVFzoMKO8coiFJElS1l6l\n6a21q2WHGeWVPciSJElSwgJZkiRJSlggS5IkSQnHIGsIzKNQ8HbXkiQpG/Ygd2xe4ytUi+P9DnAE\nHAUmGiySJEmtswe5Y5PFWW3Vqj2bkiRlLk6j1tSxwCvNdxs7cYzDLx3uOCwNBwtkSZI0eFqYRg3i\nPi3s55RrSjnEQpIkSUrYgzwUvEhNkqSB0+owEfVcpwXyXuAwYXTPEWA1MA78AfCmuP0S4KUOP0cN\nNR4HDSafJEm5M5dhIuqpTodYTAAl4ExCcQxwNXAfcCrwtfhckiRJGghZjEGe2T15EbAlrm8BLs7g\nMyRJkqSeyKIH+avAN4Bfim2LgANx/UB8LkmSJA2ETscgnwc8D/wAYVjF0zO2171TQ7lcfm29VCpR\nKpU6DEUaJJW4SJKkvOm0QH4+Pn4b+EPCOOQDwMnAC8Bi4MVaL0wLZGn0lOIyaXN/wpAkSbN0MsRi\nATAW198AXAjsALYB62P7euCuDj5DkiRJ6qlOepAXEXqNJ9/nDuBewnjkO4ENTE3zJkmSJA2ETgrk\nPcAZNdq/A1zQwftK6sxe5jY/+Sbgirj/VYQvupL6ay/eZ0DqG281LQ2fucxPvgq4ND6uAW7EvwtS\nHnifAQ2XeNfAZkteeKtpaTjVmp/8/Li+hTCFxtXAWmAroYdqL7CbcDB+qBdBSmqo1TyW8m/A7hpo\nT5E0fOYyP/kpwL7ktfuAJT2IUVJj3mdA6iN7kKXh0/b85Mn2acrJeonpE9RJw65SqVCpVHr9se3n\n8YPAcXF9ObCiK/FJ+bSHcD60QxbI0vCZy/zk+4FlyWuXxrZpyl0KVBoEM29mtXlzT+Ytb/s+A7yd\nMBGrNIpWMPWl8Ovtv41DLKThMtf5ybcB64D5hD8pK4HtvQpWUk3eZ0DqM3uQpeEy1/nJd8b2ncBR\nYCONh19I6j7vMyD1mQWyMjAvV1OzjLh25ie/Li6S8sH7DEh9ZoGsDBylcaejxbMkSRocjkGWJEmS\nEhbIkiRJUsICWZIkSUpYIEuSJEkJC2RJkiQpYYEsSZIkJSyQJUmSpIQFsiRJkpSwQJYkSZISFsiS\nJElSwgJZkiRJSlggS5IkSQkLZEmSJClhgSxJkiQlLJAldd08oFAoNFzGi8V+hylJEhCOW5LUVUeB\niSb7FKrVXoQiSVJTFshdN49CodDvICRJktQiC+Sua6nvrBeBSJIkqQWOQZYkSZISFsiSJElSwiEW\nygHHaUuSpPywQFYONBunbfEsSZJ6xyEWkiRJUsICWZIkSUpYIEuSJEkJC2RJkiQpYYEsSZIkJSyQ\nJUmSpITTvEmSpNqOwXnqNZIskCVJUm2vAuUW9mtlH2mAWCBLkiTZW66EBbIkSZK95Up4kZ4kSZKU\nsECWJEmSEhbIkiRJUsICWZIkSUpYIEuSJEkJC2RJkiQp0a0CeQ3wNLAL+HCXPqNLKv0OoIlKvwNo\noNLvABqo9DuAPBvgfIVKpdLvEOoa5NjGi0UKhULdZbxY7Gt8I2yg81Ud2NPvABrIc2xt6kaBfCzw\naUISrwIuA07rwud0SaXfATRR6XcADVT6HUADlX4HkFcDnq/5LqQGObaD1SoTUHc5WK32Nb4RNfD5\nqg7s7XcADeztdwDZ60aBvBrYTfjnOgJ8Hljbhc+R1DnzVRoc5qvUI924k94S4Lnk+T7gnC58TiZe\nfvlbwEeSlq8nzw/0PiCpt1rK14/MbEg8+fd/n3VMA228WJzWu7p58+Zp2984NsZ3Dh/udVgaDq0d\nX/8MOK7Bu7ycbVCSWvOzwGeS5+8Hfm/GPrupf+bOxWUUl930h/nq4jL3xXx1cRmcpa187UYP8n5g\nWfJ8GeFbbuqfduFzJc2d+SoNDvNVGmDzgGeB5cB84DG8iEDKK/NVGhzmqzTg3g38JaFbe1OfY5HU\nmPkqDQ7zVZIkSZI0XFqZ0Px34/bHgTN7FBc0j+1fx5ieIFwT/NbehdbyRPBnA0eBn+lFUIlW4isB\njwJP0ttJiJvFdhJwD+HU5JPAL/QsMriFMDXKjgb79CsfwHztRJ5z1nxtj/namTznrPnavrzmbN7z\ndZpjCaeAlhMmnKk1VupfAHfH9XOAh3IU248CJ8b1NTmLbXK/PwG+TLiyuVdaiW8h8BSwND4/KUex\nlYGPJXH9Ld25WLWWdxCSsl4C9ysfwHztRJ5z1nxtn/navjznrPnavjznbOb52q1bTUNrE5pfBGyJ\n6w8TfvEBJrRyAAAIZ0lEQVSLuhjTXGL7c+BQEttSeqPVieA/AHwB+HaP4prUSnzvA77I1NXVf5Oj\n2J4HJu+RWyQk79EexfcAcLDB9n7lA5ivnchzzpqv7TNfuxufx9jZ8pyvkO+czTxfu1kg15rQfEkL\n+/QiSVqJLbWBqW8e3dbqv9ta4Kb4fKIHcaWf3Sy+lcA4cD/wDeDy3oTWUmyfAd4C/DXhNMsHexNa\nS/qVD/U+23xtTZ5z1nztHvO1vjznrPnavkHO2TnnQze7vVv9D1Vo83WdmMtnvBO4AjivS7HM1Eps\nnwKujvsWmP1v2E2txHcc8DbgJ4EFhJ6Chwhjf7qpldiuIZwWKgFvBu4DfhioNnhNL/UjH+byOebr\nbHnOWfO1u8zX2vKcs+Zr+wY9Z+eUD90skFuZ0HzmPktjW7e1EhuEiwY+Qxgf1ajrPkutxHYW4dQG\nhDE+7yac7tjW9ehai+85wmmfv4vLnxISpNsJ3EpsPwb8Zlx/FtgD/DPCN/F+61c+1Pps87V1ec5Z\n87V7zNf68pyz5mt348trzvYzH2ZpZULzdND0ufRukH4rsf1jwlibc3sU06S5TgR/K729wraV+H4I\n+CphQP8CwqD5VTmJ7ZPAtXF9ESG5x3sQ26TltHYRQS/zAczXTuQ5Z83XzizHfG1HnnPWfO1ufP3M\n2eXkM19rqjWh+ZVxmfTpuP1xwmmDvMT2WcLg8kfjsj1HsaV6XSBDa/H9J8KVtjuAq3IU20nAlwj/\n33YQLnjola2EcVn/QOgFuIL85AOYr92ML9XrnDVf22O+dibPOWu+ti+vOZv3fJUkSZIkSZIkSZIk\nSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZKk0bG53wFIaso8lfLDfNQ0x/Q7AHXFcf0OQFJT5qmU\nH+ajprFAHj5jwKF+ByGpIfNUyg/zUbNYIA+fs4GH+x2EpIbMUyk/zEfNYoE8fM4G/qLfQUhqyDyV\n8sN81Czz+h2AMrES+HfAd4F3EMZSLQR+A3ipj3FJmmKeSvlhPkpD7nLgT4BFwHygHNvXAH+KZwmk\nPDBPpfwwH6Uh91PAc8B4fF4iJPikF4Gf7HFMkqYzT6X8MB/VEr8lDbZPAJ8GvhOfvx14INn+OqCY\n0Wf9c+DPgA9m9H7SqGgnTy8GPgZ8BrgsvnZ5i5+3Ebimzra3EcZa/vsW30saNt06bjbKu255H3P7\n2yCNhJOAVwkXF0z6rWR9NfAPhFNIWfl9QqEsqTXt5unTwD8h9GzNB3bN4TMXAW9osH0boVCWRk03\nj5vN8q4bXsfc/jZoDrxIb3B9DzgSHyHM41iN6wXCXYGuBw5k+JmrgCczfD9p2LWTpyuAbwP/Ly7v\nAL4xh89slPPHAmcAj87h/aRh0c3jZpbH2latZm5/GzQHx/Y7ALXtKOHb8JnAHwPvAnYDLwC/R0jW\nDwNvAX4W+G/ADwPXEiZE30n49vlfCD1V/xKYAP4KeH1sX0HowfoW8APxs+4FfiE+/p+4v6TaWs3T\nSecCv07oiVpMyLHLgW8SDoSvBzYAvwrcRegN+yzwLOFq/A3AJuAP4zYIRcBvAv8IWEcYWndHF35W\nKe9azceVhLz7V4RcXB3bdwMXAO8lDGt4HDgeuILpebcS+HngZ+Lrzybk7D0xhjfG/U8h5OQThONt\nvWP1CcBHCTl8cdzvnvgZk38bJCWOBf4z8DvAfYRvv58CfjzZ5ycIyfUiIblOJpyyBfgiUxcjLCAk\nKYQkf2ey/lPAvyGM2/rZ2L4Mx7BLrWglT1PXAP82eX4v4YsuhPx7HeGguCC2PUvI6Z9jajjGguT1\n25LP+ijTC3Jp1LSSjz9PGDJxgHDcBNgK/Me4fjrwpbheK+8mX/88UxcD/k/C2ZtjCdfznBrbNxKK\n40bH6jSHfwP4lbie/m2QVMfHGmx7F3DbjLbzmP6tcynwt7E9nTB9SXy8jXAgvo3sLvyTRk2jPJ10\nByEPIczNmp6lWUDo/fpyfL6cqQuM3gD8GOFL7aSzCT1Qk75I6KWWVD8f30DoJf500vYYU8Xu+4CP\nJ/vOzLvJ19+QtH2TcDbnpwmdUZcAvwT8SLJPrWP1zBz+X4S/ATP/Nihj9gAOh/mECwvqeQfTr9KF\nkND3J89/Avg68KPxcdL++LgaOI1w+mgD4RvuG9sPWRo5zfJ00irgqbh+NuHAfHx8/n3gPYQeJQh5\n/ADhS+v3gPcDn2PqS+yPA5Xkvc8EHiEcqKVR1igfv8f04+YbCQXp5MwXPw38EfXzbvL1D8bnKwnD\nOP6ecIHsHwN3EmapSTuqah2rZ+bwGYShHT/G9L8NypgF8nA4g5Ao9ZzH7KQ7APxdXH8d4Zvsr8b2\n7yf7vYmQoAcJf0x2xccLZ+wnqbFmeQrhwukTmbqT1zmE07GXJvssJlwXAGE85F8wNfTppwgH7kvi\n80OEC/4AzifM//pWwildaZTN5bj5DkIeQhgicS6h+F0X22bm3eTrJwvktYRrBn4O+L+ECwUnvY2p\n4Ra1jtVpDr+T0Gn1FuAsZv9tUIYskIfDGcDDDbYfDzwzo+33CQfjXyBckLeRMITiDsKFQOsJ46jO\niq+/J77uK4Qeru8Tvg1Lak2zPIVwoPzL5Pk3CMVseoX8HYRe5PcRZpU5mzDbBcCfE3J3cgjG54Ef\nJBxEi8BfEw7CTyCNtkb5eAyhiP3r+Pw04O64/gqwndCp9IXYNjPvJl//fHz+DGFM8XcJY5cLhOPr\nesJFg5PH51rH6skcXkfI4ecJZ3z/gtl/GyRJGkqX4kV0ktR39iBLUv+dSxiXeCbwv/sciySNPG8U\nIkn99z3CxXN/S7iLniRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJ0nD7/6o3\nf0pPBaqxAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x29469c70>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pr_arr = model_fit.trace('phi_r')[:]\n",
    "pf_arr = model_fit.trace('phi_f')[:]\n",
    "pm_arr = model_fit.trace('phi_m')[:]\n",
    "\n",
    "pr_arr = pr_arr[np.where((pr_arr>=0)& (pr_arr<=1))]\n",
    "pf_arr = pf_arr[np.where((pf_arr>=0)& (pf_arr<=1))]\n",
    "pm_arr = pm_arr[np.where((pm_arr>=0)& (pm_arr<=1))]\n",
    "\n",
    "# make sure that the rock percentage is confined to values where the \n",
    "# total proportion is equal to 1\n",
    "pr_arr = pr_arr[np.where(pr_arr > 1-(stats.mode(pf_arr)[0][0]+stats.mode(pm_arr)[0][0])*1.2)]\n",
    "\n",
    "fig_2 = plt.figure(2, [12, 6], dpi=300)\n",
    "ax1 = fig_2.add_subplot(1, 3, 1)\n",
    "h1 = ax1.hist(pr_arr, bins=20, color='blue')\n",
    "ax1.set_xlim(0, 1)\n",
    "ax1.set_xlabel(r'$\\phi_{rock}$', fontdict={'size':18})\n",
    "\n",
    "ax2 = fig_2.add_subplot(1, 3, 2)\n",
    "h2 = ax2.hist(pf_arr, bins=20, color='red')\n",
    "ax2.set_xlim(0, 1)\n",
    "ax2.set_xlabel(r'$\\phi_{fluid}$', fontdict={'size':18})\n",
    "\n",
    "\n",
    "ax3 = fig_2.add_subplot(1, 3, 3)\n",
    "h3 = ax3.hist(pm_arr, bins=20, color='green')\n",
    "ax3.set_xlim(0, 1)\n",
    "ax3.set_xlabel(r'$\\phi_{mineral}$', fontdict={'size':18})\n",
    "\n",
    "\n",
    "print '             Rock   Fluid    Mineral'\n",
    "print 'Mean Percent   {0:.2f}   {1:.2f}   {2:.2f}'.format(pr_arr.mean(),\n",
    "                                                          pf_arr.mean(),\n",
    "                                                          pm_arr.mean())\n",
    "print 'Median Percent {0:.2f}   {1:.2f}   {2:.2f}'.format(np.median(pr_arr),\n",
    "                                                          np.median(pf_arr),\n",
    "                                                          np.median(pm_arr))\n",
    "print 'Mode Percent   {0:.2f}   {1:.2f}   {2:.2f}'.format(stats.mode(pr_arr)[0][0],\n",
    "                                                          stats.mode(pf_arr)[0][0],\n",
    "                                                          stats.mode(pm_arr)[0][0])\n",
    "print 'Min Percent    {0:.2f}   {1:.2f}   {2:.2f}'.format(pr_arr.min(),\n",
    "                                                          pf_arr.min(),\n",
    "                                                          pm_arr.min())\n",
    "print 'Max Percent    {0:.2f}   {1:.2f}   {2:.2f}'.format(pr_arr.max(),\n",
    "                                                          pf_arr.max(),\n",
    "                                                          pm_arr.max())\n",
    "print 'STD            {0:.2f}   {1:.2f}   {2:.2f}'.format(pr_arr.std(),\n",
    "                                                          pf_arr.std(),\n",
    "                                                          pm_arr.std()) "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}