Error_Propagation_With_Sympy.ipynb

{
 "worksheets": [
  {
   "cells": [
    {
     "metadata": {},
     "cell_type": "heading",
     "source": "Error Propagation with Sympy",
     "level": 1
    },
    {
     "metadata": {},
     "cell_type": "markdown",
     "source": "Based on modified version of sympy error propagation code by Jerry Terrell (see [here](http://arengr.wordpress.com/2013/12/12/python-error-propagation-with-sympy/) )\n\nThe second example makes use of functions from Matthew Rocklin's excellent sympy.stats module, which can also be used for sampling-based uncertainty propagation (see [here](http://people.cs.uchicago.edu/~mrocklin/tempspace/scipy2012-sympystats-paper.pdf)). The second example below shows how to go about this. \n\n\nBig thanks to Jerry and Matthew!  \n\n\n\n\n<div align=\"right\"> John Griffiths 17/01/2014\n"
    },
    {
     "metadata": {},
     "cell_type": "heading",
     "source": "Importage",
     "level": 2
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "%matplotlib inline\nfrom matplotlib import pyplot as plt\n\nfrom sympy import *\nfrom sympy.stats import Gamma, E, density, sample_iter\n\nfrom sympy.interactive import printing\nprinting.init_printing()\n\nfrom uncertainties import ufloat \n\nimport numpy as np\n\nimport pandas as pd",
     "prompt_number": 322,
     "outputs": [],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "def plot_dist(ax, rv, urange, la):\n    D = density(rv)\n    f = D.args[1]\n    u = np.linspace(urange[0], urange[1], 500)\n    v = np.array([re(f.subs(D.args[0][0],i).evalf()) for i in u], dtype=np.float)\n    m = np.amax(v)\n    col = list(np.random.rand(3))\n\n    ax.fill_between(u,np.zeros(u.shape),v, color=col+[0.5])\n    ax.plot(u,v, color=col)\n    ax.grid(True)\n    ax.set_xlim(urange)\n    ax.set_ylim([0.0, 1.1*m])\n    ax.set_ylabel(la)",
     "prompt_number": 323,
     "outputs": [],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "heading",
     "source": "Error propagation function",
     "level": 2
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "def run_error_propagation(variables, variables_values, uncertainties, function_definition, function_name):\n    \n  #==========================================\n  # DO NOT MODIFY BELOW THIS LINE\n  #===========================================\n  partial_terms = [];\n  # create the partial derivative terms\n  for i in range (0, len(variables)):#size(variables)):\n    partial = diff(function_definition, variables[i])\n    partial_terms.append(partial);\n    \n  # make symbolic substitutions for the partial terms\n  partial_term_values = [];\n  for i in range(0,len(partial_terms)):\n    ptsv = partial_terms[i];\n    for k in range(0,len(variables)): #size(variables)):\n      #ptsv = ptsv.subs(variables[k], variables_values[k]);\n      ptsv = ptsv.subs({variables[k]: variables_values[k]});        \n    partial_term_values.append(ptsv);\n\n  product_terms = []\n  for t in range (0, len(variables)):\n    product_terms.append(partial_term_values[t]**2 * uncertainties[t]**2)\n    #print(product_terms[t])\n  uval = sqrt(sum(product_terms))\n    \n  # determine the value of the function_definitiontion\n  fval = function_definition; # fval = function_definitiontion value\n  for i in range(0, len(variables)): #size(variables)):\n    fval = fval.subs(variables[i], variables_values[i]);\n    \n  # put everythign into a dataframe\n  dfvs_partial_term_values = ['%f' %p.evalf() for p in partial_term_values]\n  dfvs_product_term_values = ['%f' %p.evalf() for p in product_terms]\n      \n  #dfvs_fvals = ['%f' %p.evalf() for p in fval]\n  #dfvs_uvals = ['%f' %p.evalf() for p in uval]  \n  dfvs_variable_names = ['%s' %p.evalf() for p in variables]    \n  dfvs_variable_values = ['%s' %p for p in variables_values]\n  dfvs_uncertainties = ['%s' %p for p in uncertainties]\n    \n  #dfvs = pd.DataFrame([dfvs_partial_term_values, dfvs_variable_values,dfvs_uncertainties]) #index=dfvs_variable_names)\n  df = pd.DataFrame([ dfvs_variable_values,dfvs_uncertainties, \n                      dfvs_partial_term_values, dfvs_product_term_values], \n                      columns=dfvs_variable_names, \n                      index=['variable_values', 'uncertainties', 'partial_term_values', 'product_term_values'])\n    \n  returndict = {'fval': fval, 'uval': uval, 'partial_terms': partial_terms,\n                'partial_term_values': partial_term_values,   'product_terms': product_terms,\n                'df': df}\n\n  return returndict\n",
     "prompt_number": 324,
     "outputs": [],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "heading",
     "source": "Example application",
     "level": 2
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "# define yout independent variables here\nh,r = symbols('h r')\nvariables = [h,r]; # make sure this matches the variables above\n\n# these are the initial guesses for the independent variables\n# in the order they are listed above\nvariables_values = [10,6];\n\n# this is the uncertainty in the variables above\nuncertainties = [0.1,0.1];\n\n# this is your function definition\nfunction_definition = h*pi*r**2;\nfunction_name = 'V'; # name of the function\n\nres = run_error_propagation(variables, variables_values, uncertainties, function_definition, function_name)",
     "prompt_number": 325,
     "outputs": [],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "markdown",
     "source": "Partial derivative terms"
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "res['partial_terms']",
     "prompt_number": 326,
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 326,
       "text": "⎡   2         ⎤\n⎣π⋅r , 2⋅π⋅h⋅r⎦",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAHEAAAAaBAMAAAByY656AAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAu90idpmJMlRmqxBE\n781pshE+AAABsklEQVQ4EcWUv0sDMRiG3/vB1XpXtYqTy3GDq4f+AwUH1zooukgH28GpY3FyEdza\nwc2lODjoYF1FaKmgaxFxEwoidGtRUKGF83rfJdfYi4OLGZo3+Z4nTUJaQGz5+Zo4IR/NOu5I0XAn\n2yNDaWw6OdhCVXUT78KEbKBVf5iphjaQwcL8mAlM9AUC2ka6yGauHlhCjKlmeDUI1zCiXbzyWoyZ\n50UKF8Azn/riadxMZnmRgr/BS7ZdKzrJuHmHe1Ht1iLTsHmNTN3zW9bau9nNJDpHC7zKwlKNCExt\nHc8h4MIbKpRWj1dgKB2llfK8HhNYb32ACNRfcE4cmYkDw10DSmrVYgdiUtCrMyGBchHpkAvPqcC/\nwVq9IQh8sO6ngMAO8BhyobkfvKQyZ4WgV/0hEYvQPoGAC81tDK/bXzGuHUJzERDmAEbLJI5Ms4/h\nM3mK85DIQneJSPaQyp0QR6ZuD0fRU6+3R9ZYLuRPQYT/A5xu5IgjU6ngDLBsxqvRUwG6nvcGIvxP\nYzPkyGQG7295kgaJWZEKvBBvmn82k3xleYj/TjkfVf7JdIT/22g7v6emU/0Gj858aApwfb4AAAAA\nSUVORK5CYII=\n",
       "latex": "$$\\begin{bmatrix}\\pi r^{2}, & 2 \\pi h r\\end{bmatrix}$$",
       "metadata": {}
      }
     ],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "[r.evalf() for r in res['partial_term_values']]",
     "prompt_number": 327,
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 327,
       "text": "[113.097335529233, 376.991118430775]",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAWcAAAAZBAMAAAAfyfK7AAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAu90izRCZq3bvZolU\nRDJfgOGnAAAFDUlEQVRYCe1YT2gcdRT+sjOTzf7J7lpQEdSssdiDgitKoXgwl1YRpWORQvGwW5QG\nCTYrSEpF7EJrm15sJBZR2jqgYFuFTEVi2oVmL71UcBdbRU9ZPHjJIbHNmlpT1++930zWTSqENj0E\nfIRv3+/N733vmzdvfrsE63pzWFP2Q6+P7JpSLGKt4hoWvVnugGCVz1TEdQYO1gVySG0YKJfxdPl+\nGLDfS70BBWzc+SCwfeJVGE9zQwIhCUxDvw6Vyz4SE+clqMQKgO0ysAjO8HiN61QhqLZ9YnegQAlM\n3aDTTzS4U2Ajktfpwq44D2OHZ40i3mw2560Cqp4CYs1myYBTxGTN+QSbcuqZ3JCAHKFpaIw0RedF\nvCZRJVaA1UPRLegE3uGGSAaWlMQL2JZTBVACLR6Mx0v91KtwpIIrTMIp4BvsBd5EJ1/TYtcC0lkF\n2N/u5z0JRF3ERqIjiGbUg+aGBEISmIb4OKKwfdwlQSVWcO7Z56IFuAg8xQ2HM9Bq8REk+lQBDIHU\nDWe6Szot8EvNmaOLe4GqS2/STQCJevdx5PsUWJqmEMsi+ne6D4mGeiY3JJBdgWnIB85jjA0QU2IF\n4Ag73YIzwDaq+jwDrdZRRLKhCqAEWneZaATjMQ9MXfidbDVS/iyF+KwUWqLtBkXn2ZMb6vGyjpaZ\nL65apqHuEu4zIUeIfwzY20VP3Y1z7JKdkZ1VL11E9wI9VUCC/xJtF7jJ+YOiSx+y0x5XyrCFDgjR\n03xFFbiONGLstLwFEXlWmqvAxb9MQx3s8PRBaQKUOGBvFx1tfukBnxnRW5AvoPsaE1QBCUzd4EVc\nHA9MPyKs1l/AjD8JbC3xvj2+zl8XAuhwrRtQ4L4pPzELW2infJhcQ8BIy0zoHJw5D8ckrMQBe7to\nDM5xxksiWko+4yNFKaIA8ghM3WWikTjB60a0XUpKp/OSgbP1ELBO1goPAc/hO+k0PZJLroKsWiYh\nZxZO08XrLsNKHLC3i47v3nocCZhOn60HolUBCcRYd7lofExah7c3VcLJ9ycrwAHdHDvKDwUM5ugK\ndGaA+PA0p048muQa0OUiMN7VBzmYxmoSVGLD3i76C8T/dL8PRMeO5n0dD1UgBDTWXSr6SaCnwkuc\n6ap88kCSJlo5dCwo4C0+3ZoCcIg7OGgcZ/E0NyTQKwZMKMLb4hiPeSYoxMreLpqdmbngi2itli4g\nyZboYxQCU3epaD7AHqHltjGhfZZd50mSbqDjuoIUHnQVkCyIUNg8l8TT3JCA8dBMKF+EnB6m00ps\noE20M8sefLVnz76PSlqNQuMNVcApJYGpu1T0B8DjOVY7DH53ncqlrvKFYNdjGUSuKqAArDeA34D9\n8VHM5NSD5oYEoWJ+mtAME6ucabmgxIa9dUSrN8qXjfVZSqt1jaCzTxXwCZCAf+uXz/RpJK/xewSR\ninMCB9yXPSBJ0XEP1ZICv6ypU8F6tNyf7drsPA/1oLkhAaWFpiFUC9TjO8eEXYkN+xLRl12cZF6a\nHZaSeAw76qqAd0wCrRvMdOfeK29DoXvXzhoiRTgTh+pI7HqABJacUj+Jq2ANDLlQsPlzIotxJhhP\ncw2BHNqLpiFsqjAw3l8XdiU27P2fbvDQgtSQ/GBK9Mx7ptorw+8GCpRA64Yv4mKBVXIurRLPzWmC\nmb75xVuPereeuoLMOyPaWYui4yto121suTOdvg1BK0n9X/RKurQae9jp3rX3f4/iP7D0IBsnGcmh\nAAAAAElFTkSuQmCC\n",
       "latex": "$$\\begin{bmatrix}113.097335529233, & 376.991118430775\\end{bmatrix}$$",
       "metadata": {}
      }
     ],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "markdown",
     "source": "Dataframe with partial terms and other info in "
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "res['df']",
     "prompt_number": 328,
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 328,
       "text": "                              h            r\nvariable_values              10            6\nuncertainties               0.1          0.1\npartial_term_values  113.097336   376.991118\nproduct_term_values  127.910073  1421.223034",
       "html": "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>h</th>\n      <th>r</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>variable_values</th>\n      <td>         10</td>\n      <td>           6</td>\n    </tr>\n    <tr>\n      <th>uncertainties</th>\n      <td>        0.1</td>\n      <td>         0.1</td>\n    </tr>\n    <tr>\n      <th>partial_term_values</th>\n      <td> 113.097336</td>\n      <td>  376.991118</td>\n    </tr>\n    <tr>\n      <th>product_term_values</th>\n      <td> 127.910073</td>\n      <td> 1421.223034</td>\n    </tr>\n  </tbody>\n</table>\n</div>",
       "metadata": {}
      }
     ],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "markdown",
     "source": "Final value ( + /- sqrt of sum of squared error terms )\n"
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "pprint('%f +/- %f ' %(res['fval'], res['uval']) )",
     "prompt_number": 329,
     "outputs": [
      {
       "output_type": "stream",
       "text": "1130.973355 +/- 39.359028 \n",
       "stream": "stdout"
      }
     ],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "markdown",
     "source": "Comparison with uncertainties package:"
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "from uncertainties import ufloat # .umath import \n\nuh  = ufloat(10, 0.1, 'h')\nur = ufloat(6,0.1, 'r')\nupi = ufloat(3.142, 0, 'pi')\n\nuV = uh*upi*ur**2;",
     "prompt_number": 330,
     "outputs": [],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": true
      }
     },
     "cell_type": "code",
     "input": "pprint('Sympy final value: %f +/- %f ' %(res['fval'], res['uval']) )\nprint 'Uncertainties final value: %s +/- %s' %(uV.n, uV.s)",
     "prompt_number": 331,
     "outputs": [
      {
       "output_type": "stream",
       "text": "Sympy final value: 1130.973355 +/- 39.359028 \nUncertainties final value: 1131.12 +/- 39.3641316612\n",
       "stream": "stdout"
      }
     ],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "markdown",
     "source": "Partial derivatives"
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "res['partial_term_values']",
     "prompt_number": 332,
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 332,
       "text": "[36⋅π, 120⋅π]",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAHEAAAAZBAMAAAD099zUAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAu90iiVSZZnZE7xDN\nMqtC8pX5AAAB0klEQVQ4EcWUP0jDQBTGv5pWU66COjpYreBqVHCuUEoHC5lcrSAKIpJNcWlB6CCo\nxdGlKHYpFuLkILqLgx101SKCoINQxD9YjC+5O7GJGdTBt7zeu++X7/r1UnTFNPy8NmNx9P4cswml\n9jdSSYzpQDCV9rE/JYuZzirYdL4qJdxzFKFnsBPMy3FzH3oDjhF8RURn/XKLkxs66ojE0SnHTX1/\ngch7oIxlYFxucXKrwBrI+kasErkDPFR7gElToDIhOm23fJqn2+SZSeQ7UKRAnJJkJIPGeb4A1aLK\niE3ZbJJq8PCJyJyQCPJ8AKxhoITExcrSiCRk52T4UXkEruJCIj2DFWaZmDvSg9qaBD47J1s6HDIt\nJJJE2awD2QICdoqu4uQUGHkWc0LCSTpfVN8j0sDlN3fKIdUaQN9zUhcSTtJBowZlS56z4HF8tXXI\nCShaH0lMIeHkNgWn0S81B/aG16+Q89km2zNQtVXQNRMSTi4i9IJAnJWg9mIXlHMTbZPDieQ1WnRW\nERLxroRnpgpAcqGKgIEbogblTbGf0HpXv8WZZdXBUuufEu7Z5GAvwppn5B74kG1unXftQx54le6J\nD2m4dd61D+kVeib/RMZ++X9b+wDo+JtEGTbXbgAAAABJRU5ErkJggg==\n",
       "latex": "$$\\begin{bmatrix}36 \\pi, & 120 \\pi\\end{bmatrix}$$",
       "metadata": {}
      }
     ],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "uV.derivatives",
     "prompt_number": 333,
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 333,
       "text": "{< pi = 3.142+/-0 >: 360.0,\n < r = 6.0+/-0.1 >: 377.03999999999996,\n < h = 10.0+/-0.1 >: 113.112}",
       "metadata": {}
      }
     ],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "for v_it, v in enumerate(variables): pprint(\"%s : %f:\" %(v,res['partial_term_values'][v_it]) ) ",
     "prompt_number": 334,
     "outputs": [
      {
       "output_type": "stream",
       "text": "h : 113.097336:\nr : 376.991118:\n",
       "stream": "stdout"
      }
     ],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "markdown",
     "source": "***Perfect match  :)*** \n\n\n(Pi is a constant without uncertainty so doesn't matter that it doesn't appear here)"
    },
    {
     "metadata": {},
     "cell_type": "heading",
     "source": "My functions",
     "level": 2
    },
    {
     "metadata": {},
     "cell_type": "markdown",
     "source": "Mode of a gamma distribution of conduction delays, computed from axon diamters"
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "k_n = 5.3316; k_s = 0.1\nth_n = 2.05E-1; th_s = 0.01   \nL_n = 100.0; L_s = 10.0\ng_n = 0.6; g_s = 0.1\nc_n = 5.5; c_s = 0.1",
     "prompt_number": 335,
     "outputs": [],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "# define yout independent variables here\nk, th, c, g, L = symbols('k th c g L', positive=True)\nvariables = [k, th, c, g, L]; # make sure this matches the variables above\n\n# these are the initial guesses for the independent variables\n# in the order they are listed above\nvariables_values = [k_n, th_n, c_n, g_n, L_n]\n\n# this is the uncertainty in the variables above\nuncertainties = [k_s, th_s, c_s, g_s, L_s]\n\n# this is your function definition\nfunction_definition = (k-1) * (c*th / g*L)  # = (k-1)*theta   - i.e. the mode of a gamma distribution\n\n#h*pi*r**2;\nfunction_name = 'DDMode'; # name of the function\n\nres = run_error_propagation(variables, variables_values, uncertainties, function_definition, function_name)\n\npprint('Final value: %f +/- %f ' %(res['fval'], res['uval']))",
     "prompt_number": 336,
     "outputs": [
      {
       "output_type": "stream",
       "text": "Final value: 813.979833 +/- 164.860283 \n",
       "stream": "stdout"
      }
     ],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "markdown",
     "source": "Dataframe"
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "res['df']",
     "prompt_number": 337,
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 337,
       "text": "                              k           th           c             g            L\nvariable_values          5.3316        0.205         5.5           0.6        100.0\nuncertainties               0.1         0.01         0.1           0.1         10.0\npartial_term_values  187.916667  3970.633333  147.996333  -1356.633056     8.139798\nproduct_term_values  353.126736  1576.592907  219.029147  18404.532474  6625.631691",
       "html": "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>k</th>\n      <th>th</th>\n      <th>c</th>\n      <th>g</th>\n      <th>L</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>variable_values</th>\n      <td>     5.3316</td>\n      <td>       0.205</td>\n      <td>        5.5</td>\n      <td>          0.6</td>\n      <td>       100.0</td>\n    </tr>\n    <tr>\n      <th>uncertainties</th>\n      <td>        0.1</td>\n      <td>        0.01</td>\n      <td>        0.1</td>\n      <td>          0.1</td>\n      <td>        10.0</td>\n    </tr>\n    <tr>\n      <th>partial_term_values</th>\n      <td> 187.916667</td>\n      <td> 3970.633333</td>\n      <td> 147.996333</td>\n      <td> -1356.633056</td>\n      <td>    8.139798</td>\n    </tr>\n    <tr>\n      <th>product_term_values</th>\n      <td> 353.126736</td>\n      <td> 1576.592907</td>\n      <td> 219.029147</td>\n      <td> 18404.532474</td>\n      <td> 6625.631691</td>\n    </tr>\n  </tbody>\n</table>\n</div>",
       "metadata": {}
      }
     ],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "markdown",
     "source": "Comparison with uncertainties package"
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "from uncertainties import ufloat\n\nuk = ufloat(k_n, k_s, 'k')\nuth = ufloat(th_n, th_s, 'th')\nuc = ufloat(c_n, c_s, 'c')\nug = ufloat(g_n, g_s, 'g')\nuL = ufloat(L_n, L_s, 'L')\n\nu_DDmode = (uk-1) * (uc*uth / ug*uL)  # = (k-1)*theta   - i.e. the mode of a gamma distribution",
     "prompt_number": 338,
     "outputs": [],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "pprint('Sympy final value: %f +/- %f ' %(res['fval'], res['uval'] ) )#  * res['fval']) )\nprint 'Uncertainties final value: %s +/- %s' %(u_DDmode.n, u_DDmode.s)",
     "prompt_number": 339,
     "outputs": [
      {
       "output_type": "stream",
       "text": "Sympy final value: 813.979833 +/- 164.860283 \nUncertainties final value: 813.979833333 +/- 164.860283133\n",
       "stream": "stdout"
      }
     ],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "markdown",
     "source": "Partial derivatives: "
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "for v_it, v in enumerate(variables): pprint(\"%s : %f:\" %(v,res['partial_term_values'][v_it]) )",
     "prompt_number": 340,
     "outputs": [
      {
       "output_type": "stream",
       "text": "k : 187.916667:\nth : 3970.633333:\nc : 147.996333:\ng : -1356.633056:\nL : 8.139798:\n",
       "stream": "stdout"
      }
     ],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "u_DDmode.derivatives",
     "prompt_number": 341,
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 341,
       "text": "{< th = 0.205+/-0.01 >: 3970.6333333333337,\n < g = 0.6+/-0.1 >: -1356.6330555555553,\n < k = 5.3316+/-0.1 >: 187.91666666666666,\n < c = 5.5+/-0.1 >: 147.99633333333333,\n < L = 100.0+/-10.0 >: 8.139798333333333}",
       "metadata": {}
      }
     ],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "markdown",
     "source": "***Another perfect match! :)***"
    },
    {
     "metadata": {},
     "cell_type": "heading",
     "source": "Mixture distribution",
     "level": 2
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "L_n = 100.0; L_s = 10.0\ng_n = 0.6; g_s = 0.1\n\nkm_n = 5.3316; km_s = 0.1\nthm_n = 2.05E-1; thm_s = 0.01   \ncm_n = 5.5; cm_s = 0.1\n\nku_n = 3; ku_s = 0.1\nthu_n = 1.5E-1; thu_s = 0.01   \ncu_n = 3.2; cu_s = 0.1\n\nfm_n = 0.7; fm_s = 0\nfu_n = 1-fm_n; fu_s = 0\n\n# define yout independent variables here\nL, g, km, thm, cm, fm, ku, thu, cu, fu = symbols('L g km thm cm fm ku thu cu fu', positive=True)\n\n\nD_m = Gamma(\"DelsM\", km, (cm*thm)/L)\nD_u = Gamma(\"DelsU\", ku, (cu*thu)/g*L)  # i.e. the mode of a mixture distribution\nD_mu = fm * D_m + fu * D_u\n\n\n\n# Uncertainty propagation\n\n\n# this is your function definition\nD_m_function_definition = E(D_m)\nD_m_function_name = 'D_m_expectation'\nD_m_variables = [L, g, km, thm, cm]\nD_m_variables_values = [L_n, g_n, km_n, thm_n, cm_n]\nD_m_uncertainties = [L_s, g_s, km_s, thm_s, cm_s]\nD_m_up = run_error_propagation(D_m_variables, D_m_variables_values, D_m_uncertainties,\n                               D_m_function_definition, D_m_function_name)\n\n\nD_u_function_definition = E(D_u)\nD_u_function_name = 'D_u_expectation'\nD_u_variables = [L, g, ku, thu, cu]; # make sure this matches the variables above\nD_u_variables_values = [L_n, g_n, ku_n, thu_n, cu_n]\nD_u_uncertainties = [L_s, g_s, ku_s, thu_s, cu_s]\n\nD_u_up = run_error_propagation(D_u_variables, D_u_variables_values, D_u_uncertainties, \n                               D_u_function_definition, D_u_function_name)\n\n\nD_mu_function_definition = E(D_mu)#E(fm * D_m + fu * D_u) # i.e. the mode of a mixture distribution\nD_mu_function_name = 'D_mu_expectation'\nD_mu_variables = [L, g, km, thm, cm, fm, ku, thu, cu, fu]; # make sure this matches the variables above\nD_mu_variables_values = [L_n, g_n, km_n, thm_n, cm_n, fm_n, ku_n, thu_n, cu_n, fu_n]\nD_mu_uncertainties = [L_s, g_s, km_s, thm_s, cm_s, fm_s, ku_s, thu_s, cu_s, fu_s]\n\nD_mu_up = run_error_propagation(D_mu_variables, D_mu_variables_values, D_mu_uncertainties,\n                                D_mu_function_definition, D_mu_function_name)\n\n\n#pprint('Final value: %f +/- %f ' %(res['fval'], res['uval'] * res['fval']) )\n\n\n\n# Get expectations and samples of the myl, unmyl, and mixture\nD_m_eval = D_m\nfor i in range(0, len(D_m_variables)): D_m_eval = D_m_eval.subs(D_m_variables[i], D_m_variables_values[i]);\n\nD_u_eval = D_u\nfor i in range(0, len(D_u_variables)): D_u_eval = D_u_eval.subs(D_u_variables[i], D_u_variables_values[i]);\n\nD_mu_eval = D_mu\nfor i in range(0, len(D_mu_variables)): D_mu_eval = D_mu_eval.subs(D_mu_variables[i], D_mu_variables_values[i]);\n    \nD_m_eval_E_ana = E(D_m_eval)\nD_m_eval_E_samp = E(D_m_eval, numsamples=10000)\n\nD_u_eval_E_ana = E(D_u_eval)\nD_u_eval_E_samp = E(D_u_eval, numsamples=10000)\n\nD_mu_eval_E_ana = E(D_mu_eval)\nD_mu_eval_E_samp = E(D_mu_eval, numsamples=10000)\n\ndf_samples = pd.DataFrame(data=[list(sample_iter(D_m_eval, numsamples=10000)),\n                                list(sample_iter(D_u_eval, numsamples=10000)),\n                                list(sample_iter(D_mu_eval, numsamples=10000))],\n                          index=['D_m', 'D_u', 'D_mu']).T\n \n",
     "prompt_number": 342,
     "outputs": [],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "markdown",
     "source": "Sympy uncetainty propagation result:"
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "pprint('D_m - sympy final value: %f +/- %f ' %(D_m_up['fval'], D_m_up['uval']))\npprint('D_u - sympy final value: %f +/- %f ' %(D_u_up['fval'], D_u_up['uval']))\npprint('D_mu - sympy final value: %f +/- %f ' %(D_mu_up['fval'], D_mu_up['uval']))",
     "prompt_number": 343,
     "outputs": [
      {
       "output_type": "stream",
       "text": "D_m - sympy final value: 0.060114 +/- 0.006870 \nD_u - sympy final value: 240.000000 +/- 50.519798 \nD_mu - sympy final value: 72.042080 +/- 15.153941 \n",
       "stream": "stdout"
      }
     ],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "markdown",
     "source": "Analytic vs. sampling comparisons:"
    },
    {
     "metadata": {},
     "cell_type": "markdown",
     "source": "Value(expectation) of the mode\n"
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "fig = df_samples.hist(layout=[1,3], figsize=(12, 3), bins=100)",
     "prompt_number": 344,
     "outputs": [
      {
       "output_type": "display_data",
       "text": "<matplotlib.figure.Figure at 0x64a8ad0>",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAs0AAADSCAYAAACxUP3mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX2YJVV95z8N3SNM6LanH7Izw4s2KixOVm2yAq4vocco\nOyQGcN0VjMaZgHl4HvMIg0aZIZsd1MhbFCfRhefZBcxohJXHFxw2AUHSlagJkBcaUGRhDL2KwmBk\ncFqBlYa7f5xTXXWr695b996qU6eqvp/nuU+9narfqbr3d8+pc77nd0AIIYQQQgghhBBCCCGEEEII\nIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIUStWACeAvYD+4Bv\nAecAIyXmSQjRnQXkt0LUkQXk20J4y8PAG+36OPBbwL8A15aWIyFEL+S3QtQT+bYQHhN30JDjgeeA\nX+ly3p8DVwJ/BSwC3wDWAX+KeTv+LjCTc16FEAZXfvs88JLE+R8dPNtCiB4M6tsBcHZsewvGv4WH\nHFB2BkSu/APwCPCGHun+C/CHwKHAL4A77LlTwBeBKwrMoxCiHRd+27IfIYQ7svi2fLNCqNJcP34E\nrOlyvAV8Gbgb+H/AV4CfA39hj90AHFdwHoUQ7bjwW2krhXBPL98WFUKV5vpxOPBEjzSPx9afSWw/\nDRySd6aEEF2R3wpRT7L4tqgIqjTXi+MxDvrNsjMihMhMEX77FLA6tr0edQEL4Zosvv1z4Jdi2+sK\nzZEYClWaq03Y3ToBvAW4Hvgc8J0M5wghysGF384D7wQOBDYBv9bn+UKI/hnEt+eB/wQcDLwMMyhQ\nL7iekrXSfCBGS3eT3Z4CbgMeBG4FJmNptwMPAQ8AJ+eTTdGBmzAxIb+Pee6fAH63xznJQQdpgxDk\nsPUk6ccXYQap3G0/p8TSyo+Lw4XfnocJebUP+G2MBlrUA5XH/jKIb38SM7B3L/AZzDgFUXHeD3we\n2G23Lwc+ZNcvAC616xswb01jwDSwB7VmC+ELST/eYfclkR8L4S8qj4UoiSwOdATwG8DVRF0PpwK7\n7Pou4HS7fhqmO+JZzOw4e4ATcsqrEGJw0vx4hPRuf/mxEH6i8liIEslSaf4k8EFMsPyQtZiuBOxy\nrV0/DNPdG/IIRgQv3PIdzAQIyc87ysyUKJU0P24B7wPuAa4h6taVH5eD/Fb0QuVxNZFv14Relea3\nYMIa3U3ngSi9AnOvOHbYYYeF5+hTzGcDJvxU8nOdB3lr0mcPftDJj68CjsLMJPcoRn/XiVZyh/w4\n94/81s+P734cJ8xzJ1Yckx87+ci3y//k4se9Ks2vxXT9PIzp5nkjZiToXqKwKOuJ4oX+EDgydv4R\ndl8bP/rRj2i1Wl59duzYUXoeqpIv5SnbB3jpoI6ZM2l+/FmM34Z/KFcTdd1658euvl+XvyPdk/92\nWi3v/bjS5XGZ/9uy3Szb5OTHvSrNF1qnOwo4E/hr4HcwAxA22zSbgRvt+m6bbpU952jgrjwyKoQY\nmDQ/fjemgA15K3CfXZcfC+EfKo+FKJnRPtO37PJSzLStZ2MGGLzd7r/f7r8fWALeGzvHaxYWFsrO\nQio+5kt5qjQjRD55OfAqu/0wcI7d750fu/p+Xf6OdE/+2/GcypfHZX6Pst0s23nRT6X5b+wHzJSQ\nb+qQ7mL7qRQzMzNlZyEVH/OlPFWawH7AtFJ1wis/dvX9uvwd6Z78t+MxXpfHExNTLC7uY3x8Dfv3\nd55BuszvUbabZTsvypodrmU1JkLUlpGREaj3DIzyY1F75Mf9Y55ZCxhB/xHCB/LyYwU6F0IIIYQQ\nogeqNFuCICg7C6n4mC/lSbjG1ffr8neke/LfjiiWMr9H2W6W7bxQpVkIIYQQQogeSNMsREFICylE\n9ZEf9480zcI3pGkekomJKUZGRpiYmCo7K0IIIURtMeXtKpW5ovI0ttK8uLgPaNmlv1obH/OlPAnX\n1FErq3vy347IB1POPku8zIXm6mtlu7o0ttKchlqfhRBCCCFEGo3VNEeaqzFgifHxNcutz9JhiTyQ\nFlKI6iM/7p+kptk+Q5JlbreJT4TIE2mac2OJZJeREEIIIbKTvadWZa6oLqo0e46PGiDlqbIcCNwN\n3GS3p4DbgAeBW4HJWNrtwEPAA8DJDvOYSh21sron/+2I7CTHCWWhqfpa2a4uqjQL0RzOA+7H9JEC\nbMNUmo8BbrfbABuAM+xyE3AlJf5XTExMsXHjRo01EEIIUSq99B0HAX8DvABYBXwV0wJ1EfAe4Mc2\n3YXAzXZ9O3AW8BxwLqYFK0npWsi45ipagjTNIi8800IeAfw58DHg/cBvYVqRTwL2AuuAADgW48PP\nA5fZc2/B+PwdiWs68WPFfBVl4pEfe10et5epENcvx0wRL3Nbrdbyy7D0zaJI8vLj0R7HnwE2Ak/Z\ntN8EXo/5xV9hP3HiLVSHA1/HtGI9P2xGi2eUducWolZ8EvggMBHbtxZTYcYu19r1w2ivID+C8Wch\nRHlUpDyOl6VLtFekVyJts6gSvSrNYBwUzJvtgUD4C0/zgtOA6zEBGReAPcAJrGyh8hA/K8xBEDA7\nO1t2NtpQnirHW4DHMXrm2Q5pWkSyjU7HV7Blyxamp6cBmJycZGZmZvl7CPVrw24bguV9eV8/vj0/\nP8/WrVsLu358e+fOnYU8r7TnNzs7W/j9uHx+yXvL8/rh+sLCAh5SgfI4e1latr62rDJDtuvNAcA8\nsAhcbvftwDjhPcA1RAOIPgW8M3bu1cDbUq7ZKhugBclltO4Lc3NzZWdhBcpTNuheCXXJxcAPgIeB\nR4GfA5/DyDPW2TTr7TYYbfO22Pm3ACemXNfhc5xz4pcuf0eubOmehgN//Bg8Lo/pUJZ22zc3N7d8\nzDVllhmy7R5y8uMsLc3PAzPAC4GvYVqqrgI+Yo9/FPgEcHaH80trocrWggVhK1bbnoJbtLJuu2oh\n6mc73OdLfpLfZ5n2PW2hutB+wGiY/wD4HUyhuxmjXd4M3GjT7Aauw3T3Hg4cDdzlML9MTEyxuLiP\n8fE1ds+sE7suW0Fc2dI91Qpvy+POpB1vv0bbEU/Kk7qWV+G+utdnwvW8y+N+RdF/BDwNfDy2bxoT\nwuoVRK1Tl9rlLZi34DsT17EV//LQQEBRNB4NIIpzEvAB4FRMyLkbgBdhWqreDjxp012IGUC0hIm6\n8bWUaxXmx+mDiuSXwj2e+jF4Vh538tlu+1qxiU/k26JIXE1ucihRV8/BwJsxush1sTRvBe6z67uB\nMzF6q6MooYUqC1lCV/kypXbvt3j3KE+V5m8wFWaAJ4A3YQYHnUxUYQYj6XgZJppGWoXZMYEbKw5/\nR65s6Z5qQ+3K4+h7HHVe3pb5G5Lt6tJLnrEe2IWpXB+A0UHeDnwW00XUwugkz7Hp78e0XN2PaaF6\nL37pwYBso3WjQO0+NjAIUS/icgyFnhIilVqWx4ZwlkCVt8JvyvqFlirPCLuDuncfRevqNhKD4HG3\nbl7k5sfJWMySZwhfkB9nY1h5hnxcFIkreYYQQniC+y5cIYQQIqRRleZQp1wlfNQAKU/CPQFRF25x\nkyFI/1sNW/L3OjBaanncVG1vU23nRWMqzaFm0ltJlxAihbRhF2pxFqL6LAFzZWdCiL5ojKZ5pY45\nvi5Ns8gfaSF7E73MQhZNc5omUogikR9nYxBNs8pb4QppmoUQlUe9P0IIIaqCKs2e46MGSHkS7gnc\nWJH+txK25O91ISjPckO1vU21nReqNAshhBBCCNEDaZqlaRYF4ZEW8iDMTIAvwMwO9lVgO3AR8B7g\nxzbdhcDNdn07Zhrt54BzgVtTrju0H6dNZy9Ns/AJj/y4KDzSNI8BS5rkSOROXn7ca0ZAYUfqy4lF\nhXkG2Ag8hfH5bwKvx5RSV9hPnA3AGXZ5OPB1zFTbzzvKrxCiIrQP5h0WzQwo/EbyjJ4UHxu2Gz5q\ngJSnSvKUXa4CDgTCH3Ra6XQacD3wLLAA7AFOKDh/PQh6ppiYmBo6DJ30v9WwJX/3h+EG8wY55qRP\nyw3V9jbVdl6o0ixEMzgAmAf2YoKjfsfufx9wD3ANMGn3HQY8Ejv3EUyLs9csLu4r7eVWCCFE/ZGm\nOYOmWfpJMQieaiFfCHwN2AbcT6Rn/iiwHjgb+BRwB/B5e+xq4K+ALyeu5ZWmOfRx+anIE0/9OE+G\n8uM0H7aXzbBP5a1wgytNc6cBRFPAF4AXY7pv3w48ac/JMoBICFEOPwX+Eng17X2jVwM32fUfAkfG\njh1h961gy5YtTE9PAzA5OcnMzAyzs7NA1BXXazsiAA6IveAmj6/cFwTB8vVWpMxoX9vajm+H6wsL\nC3hGg8rjUSYmpjSOSFSS1XY5iml9ej1wOfAhu/8C4FK7vgHTBTwGTGO0kGkSkJZrMK+vLUhbz3J8\ntDU+vsZ5vufm5pzb7IXylA38mbXjUCLpxcHA3wK/DqyLpTkfuM6uh368CjgK+B7pb+gDPZfx8TXL\nvkSqH6b55VzqsZTnPdR35vJ35MqW7mk48MePwaPyeHx8TcaytNu+ua7pi6TMMkO23UNOfpwlekba\nAKJTgZPs/l2Y5p9tdB5AdEcemS2XJeklRVVZj/HTA+znc8DtwGeBGcyfycPAOTb9/cANdrkEvJcc\nKw7yIyEGxpvyOBoAWGflihDtZPm1HwD8M/BS4CrMG+0+YE3sGk/Y7TQt5M3AlxLXtBV/d+ShaQbp\nJUV2pIVMJ649TtdDZtdCtlqtFSGv5KMiTzzzY2/K4+6+S8Z9Km+FG1zGaX4e0xoVDiDamDjeq9k7\n9VgeWsjBtZNJ+jvui/ZO235th+seaiE9ZDT2IjscavESDcKb8ridtH3djmVP78v/u7artR2ul10e\n/xHwB8ADRHrI9XYbTJfQtlj6W4ATU65Tlp5lCE1z8RqrNHzU6ipP2cAvLWQRDPNcMmiZs2mak2mH\nQfrfathqsKY5TqnlMT19Vppm2fbHNjn5ca84zckBRG8G7gZ2A5vt/s3AjXZ9N3Am0QCio4G78sio\nEEII0WBUHgtRMr36M1/BygFEf4IJcXMD8CJWhri5EBPiZgk4D9OFlMRW/N0hTbNwjWdayCIYyI9X\n+uJwmuZkTGf5qMgTj/zYq/K4eE3zKLDE+PgahZ4TQ5OXH2tyE1WaRUF4VNgWhSrNovbIj9NxMRAw\n7udCDENeftxLniFKJn3QRbkoT8I9gRsrDn9HrmzpnoS/BOVZLvE3JNvVpfaV5omJqQFG6WcJKiKE\nKJf8InAIIbIzMTFVdhaEKIXayzOSXbd5yDPCP4y4ziptn2g26tZNJ095hnSQomjkxyvp7cNpxwZP\nL3mGGBaXcZrFMp1btjTLmRA+sAS0WFyscx1HCCFEGdRenpEvpkB2iY8aIOVJuCdwY0X630rYkr/X\nhaA8yw3V9jbVdl7UstI8MTElzZUQEQcBdwLzwP3AJXb/FHAb8CBwK1EMWIDtwEOYiRJOdpZTIYQQ\nwlNqqWk2EgqjbbTmyEvTnKZvjkszpL0SIZ5pIVcDT2Ec45uYmcROBf4VuBy4AFiDmUFsA3AdcDxw\nOPB14BjMFL5xPNQ0R2nliyIPPPPjIpCmWdQehZzriRsphakwy6GF9zxll6uAA4F9mErzLrt/F3C6\nXT8NuB54FjNZwh7gBFcZFUIIIXykxpXmeuCjBkh5qiQHYOQZe4E54DvAWruNXa6164cBj8TOfQTT\n4lwiQZ/pzaDdfmVa0v9Ww5b8vS4E5VluqLa3qbbzQtEzBkYxYkWleB6YAV6ImUp3Y+J4i+5dJqnH\ntmzZwvT0NACTk5PMzMwwOzsLRH+Qye2IfrezpgHT0zTH4mJ0m53yE9+en5/vmf+8tufn5wu9fvJ5\nF30/rp9fkc8rCAIWFhYQQog4NdY0w7CaZU2zLYbBYy3kHwFPA+8BZoHHgPWYFuhjMbpmgEvt8hZg\nB2YwYRyvNc3SQ4o88NiP88JzTfMYir0uhkWaZiFEVg4lioxxMPBm4G5gN7DZ7t8M3GjXdwNnYvTP\nRwFHA3e5yqwQQkSEsdc1F4IonyyV5iOJNJDfBs61+y/CaB3vtp9TYucoXFVO+KgBUp4qx3rgrzGa\n5juBm4DbMS3Jb8aEnHsjUcvy/cANdnkz8F5KH+0auLEi/W8lbDXU32tYFgd9pB1snEJHyw3V9jbV\ndl5k0TQ/C5yPKXAPAf4JE9u1BVxhP3E2AGfYZbdwVUIIN9wH/GrK/ieAN3U452L7qTimoFXXrqgB\nDS+LNdunKJ9Bfn03Ap8GXgf8DPhE4vh2jFNeZrdvwbwJ3xFLU4imORkzWZpmUSbSQqbjWtMsbbMY\nBo/9OI+yGPrw495lbDFxmuXLYljK0jRPA8cROd37gHuAa4g0k6WFq1LMZCGEEA1gmhLKYpWxoun0\nE3LuEOCLwHmYt9qrgI/YYx/FvOWe3eHcFV42SKiq7KGskrg9nmcopPi9lR2KKdzeuXNnLt9X3qGu\ntm7dWmp+wnWFqiqCwI2VIFj+XutiS/dUO3Iti6G/8jhbKMhex4IMxzrvi3//w/xfx8vYppRXZZbf\nLusz4XpZ5fEYJrbr1g7HpzG6STDhqrbFjt0CnJhI3yoCoAXhMr6etq+Y43kzNzeX+zWHRXnKBvVv\nkhnmuaT4T6clLZjrI236sSy4/B25sqV7Gg788uO8y2Low4/J6GvZ9831mT67L/eizDJDtt1DTn6c\nRd8xgpli9yeYQQgh64FH7fr5wPHAb2MGHVyHmXY3HHzwskSG7T3ki9GsZNVIFXO8iPsS1cRjLWRe\nDOTH0jSLKuGRHxdRFkMffty7jM1afg6XXr4s+iUvP84iz3gd8C7gXkw4G4ALgXdgZhhrAQ8D59hj\n8XBVS3gRrkoIIYSoNCqLhSiZLAMBv2nTzWAGHhyHid36buCVwKuA04G9sXMuxrzRHovpShIDEtfn\n+ILyJNwTuLGimMaVsNVQf69hWRwMcM5oLrGay/wNyXZ10YyAQggnTExMxaQZQggxCEuaHVCURlkl\nmDTNovZ4pIUsir78OLseUppm4Q/y4whfNM020/3co2g4ZcVpFkIIIYQQonGo0uw5PmqAlKfKcSQw\nB3wH+DZwrt1/EWbCg7vt55TYOduBh4AHgJNdZbQzgRsr0v9Wwpb8vS4E5VluqLa3qbbzop/JTbwm\nj4EBQtSUZzGhqOYxEyP8E3Abpq/zCvuJswE4wy7DUFXHYKbkFUIIIRpJbTTNK2O/xtddaZpHgSXG\nx9ewf/8Ted2aqCgeayFvBD6NCWH1M8wMYnG2YyrIl9ntWzCt0nck0knTLGqPx36cF9I0i9ojTbOX\nLAEtjewVPjONCVUVVoDfB9wDXANM2n2HYWQbIY9gWpyFEMIbwog86mkWrqiNPKOuBEGwPKe6LyhP\nleUQ4IvAeZgW5quAj9hjH8W0OJ/d4dzUZp0tW7YwPT0NwOTkJDMzM8vfQ6hfa/9egg7rvbaDDGmS\ntB9L5idte35+nq1bt2ZOP8z2zp07uz6vvLbDfUXfj8vnl7y3PK8fri8sLCCKJhjqbNNA1WJxsf8G\nxDLLDNkW/VLUvOItus5b7+54XpQ5V3snlKds0KGiWRJjmMkNtnY4Pg3cZ9e32U/ILcCJKecM8Dy6\n+U+nJS2Y6yNt2rHRFtAaH1/TNY8uf0eubOmehgO//LgI+nwWWXw46765PtNHx5L56ZcyywzZdg85\n+bE0zQUdz/v+RPXwSAs5AuwCfoIZEBiyHnjUrp8PHA/8NmYA4HXACUQDAV/Gyj+dvvy4bE2zWY6h\ncQeiHzzy46LI7Me+aZrj+VGZK7qRlx9LniFE/Xkd8C7gXkxoOYALgXdgpuRtAQ8D59hj9wM32OUS\n8F5q09oWjjuocx1IiLozmjq7aKht1guxKAoNBCyE0WXnHXaggo9xDZWnyvFNjK/PYAYBHgfcDLwb\neCXwKuB0YG/snIsxrcvHYmQdJRO4saKYxpWwJX93S1iO5U8w4Hnm5TfJ4uK+zAPxmxqvuKm28yJL\npbnTxAhTmFivDwK3Eo28B+8mRnDN0rLjRgMVFFFDCCHEwJRWFoflmBBNJ8ur4zr7iU+McDrwu8C/\nApcDFwBrMIOHQj3k8XSeGKH2mmaQ5qrpSAvZjh+aZvmj6A+P/LiIshgy+HG77/qhaU7z57AeIN8W\nSVzGaX4M46RgwlR9F+OAp2IGF2GXp9v104DrMbOQLQB7MAOKhBBCCDEYKouFKJl+Nc3TGD3kncBa\nIg3kXrsNmhghV3zUAClPwj2BGyvS/1bClvy9LmVxUJ7lhmp7m2o7L/qJnnEI8CXMxAiLiWO9YuCt\nONbPpAj9BKVPx6/jLiYZKHJ7fn7eq/wEdlKFsvMTrmtShHYmJqak6RciP3ItiyFbeRyR3O60r7z0\n8Uk0etUXmlZe+Vh+F7EdruddHmfVd4wB/xsz4n6n3fcAMIvpMlqPGaBwLNGkCJfa5S3ADswbcYg0\nzaL2eKSFLIpMfty/HtJN2vHxNSwu7lPMZtEVz/w477IYpGkWDcClpnkEuAYTs3VnbP9uYLNd3wzc\nGNt/JrAKOAo4Grhr2IwKIUSeKLKNqBgqi4UomSyV5nBihI2YiRHuBjZh3l7fjAlz80ait9n4xAg3\nU6uJEdyzsmusfJQn4Z4g5+ulK9Nc/o5c2dI91YYalsVBeZYbqu1tqu28yKJpDidGSONNHfZfbD9i\nGTPhibqBhfCBpbIzIES/qCwWomTK0mk1UtMc7hPNwCMt5JHAZ4F/g/kx/g/gzzCTInwBeDEmJNXb\ngSftOduBs4DnMJMo3Jpy3UprmqO0Y8CStM0iFY/8uCikaRa1x6WmWQhRbZ4Fzgd+BXgN8PvAyzED\nhW7DTHhwO9HAoQ3AGXa5CbiSWv9XmCl5pW0WQgjRjRoXhPXARw2Q8lQ5ajApQuDGivS/lbAlf68L\nQXmWG6rtbartvKh8pXliYiomzfCJUU/zJRrONLWYFEEIIYRwSz+Tm3hJGDbKP8mZ6fIdNl9hwG6f\nUJ4qS8mTIkB7y1LyeLft2YznZL1e7/OLDsIf7it7EoAiJhUo2t7s7GzlJkUQaczmdJ3+G6nKLDNk\nu7pUfiBgcqCdTwMB045rgEJz8GwAUSmTIkAVBgK2DyYSIo5nflwEtRgIqHJWdEMDARuCjxog5aly\njFD5SRECN1ak/62ELfl7XQjKs9xQbW9TbedF5eUZQoiehJMi3IuZEAFMSLlLMZMfnE0Ucg7aJ0VY\nwstJEYQQQgi3SJ4heYYoiKZ3605MTMXCuEmeIapJ0/0YqibPGEVx10WSxssz/I2aIYSA+CDdqmAG\nE01MTJWdESHEwCjuuiiOylaaq1cgD4aPGiDlSbgncGBjCZhzVtjWUf9bx3sSRRIUeO3uL8FN1fY2\n1XZeSNMshBBCiJoRtjirR1rkR5Zf07XAbwKPA6+w+y4C3gP82G5fiAllBWaA0VnAc8C5wK0p1xxa\n09xbY4WXx6WZbA5N10Km+6jfmmZpm0USj/y4iLIYaqdpXnlM/ixcapo/A2xK7GsBV2BmFjuOyEk3\nAGfY5Sbgyow2GoI0k0IIIQZCZbEQJZPFib4BpIn80mrspwHXA89iQljtAU4YNHP1o/8BCj5qgJQn\n4Z7AoR3zcjsysqrQl9w66n/reE8eUcOyOCjPckO1vU21nRfDvHm+D7gHM2nCpN13GPBILM0jwOFD\n2BBCCMeYl1tT39AofOE9KouFcMSgAwGvAj5i1z8KfAIzQUIaqWKiLVu2MD09DcDk5CQzMzPL85KH\nbyO9tiOS20n8Ox4EQab7nZ2dzfw8XG33k3+X2/G8lWU/CAIWFhYQeTNbMzvt/lQHOy5tubwnzxm6\nLIZs5XFEcrvTvizpZ/tM3+/17ZbKKy/Kb5f1mXA97/I4qyh6GriJaPBBp2Pb7L5L7fIWYAdwZ+Kc\nxg4E1MCE5uDRAKKiqO1AQA0kEiGe+fE0+ZbFoIGAogGUPbnJ+tj6W4H77Ppu4ExgFXAUcDRw18C5\nEylv+eWjPAn3BDWzU0/9bx3vyXMKL4uLHbgeFHjtHpYbqu1tqu28yCLPuB44CTgU+AHmbXUWmMG8\nxj0MnGPT3g/cYJdLwHvp0iXUXMxAI03zKRxSVLgqIYQbSimLpekXIqKsLqfGyzPUbVR/POvWfQPw\nM+CzRJXmHcAiJmRVnA3AdcDxmMFDXweOAZ5PpJM8Q9Qez/y4CDL4MUieIapM2fIMkQujitksXFHD\ncFVCCNGL0cLDR4rmoEpzqSz17PryUQOkPNWKioSrCmpmp5763zrekyiSwIGNJaLwkYvLleemanub\najsvBg05J4SoPoWGjjQEiTOCjMfStrOmGeZ68x2P5x0aaX5+Ptfr+RDaan5+vvRQXr6GqhI+EE4w\nVme1jSgSaZpL1jQD0lvVFA+1kNM4DB1ZH03zGLCkgbsNxUM/zpvaa5rT9qncbRbSNAshhkWhIzMR\ntk4pioBoDhMTU7EKsxACVGn2gNGuAxR81AApT5XkeuDvgH+LCVd1FnAZcC9G03wScL5NGw9XdTN9\nhqsqprANcr5e2Xbqqf+t4z01FfOSWHRrbFDw9btYbqi2t6m280Ka5tKRxko44R0p+67tkv5i++mb\nqLDVb1oIIUR9qJSmeWJiisXFfYyPr0kUzD5qlnsdX7lPGqt60VQtZHctsy865X7TykebivxYmmZR\nffLy40q1NIcV5cXFsbKzIoQQQgghGkRFNc1LZWfAGT5qgJQn4Z7AAzv5TkZUR/1vHe9JFElQnuWG\nanubajsvKlppFkII1/SejEgIIUR9qZSmuT+Nle/HpWmuO9JC+q5T7idtdI78tFnIj6VpFtXHZZzm\na4G9RDFcAaaA24AHgVuJpt8F2A48BDwAnDxsBptDvl2/QojiCcPryXeFA1QWC1EyWSrNnwE2JfZt\nwzjqMcDtRDOIbQDOsMtNwJUZbYgOXb8+aoCUJ+GewBM77XHVo8HJ/cs26qj/reM9eUQNy+KgJLvd\n50comqbqiuvgt1mc6BtAskQ4Fdhl13cBp9v10zCTKDwLLAB7gBOGzqUQQniBZgcUpaGyODeWgDn5\nseibrPoQq8TiAAAS/klEQVSOaeAm4BV2ex+wJnaNJ+z2p4A7gM/bY1djZhT7UuJ60jR3OEc6q/og\nLaTvOuV+0qaPQYjfq3y3nnjmx9PkWxZDQzXNZjnG+Pg4+/c/kfJYRJ3wKU5zy366HV/Bli1bmJ6e\nBmBycpKZmRlmZ2eBqAk/ud2ZKh4PUtZH7Rd7IPAc4+Nr2L//iY7PQ9t+bYfrCwsLeMi1wG8CjxMV\nuFPAF4AXY1qj3g48aY9tx0y1/RxwLkYvKTpifDf0WSFKYKCyGNLL484EGfdVIb2RRQZBUHr5oe18\nt8P1ssrjadoHHzwArLPr6+02GD3Vtli6W4ATU67XGgSgBeEyvp62z/fj2c6Zm5sb6FkVifKUDboX\nYK55A3Ac7X58OfAhu34BcKld3wDMA2MY399DupSry333WmZJE087V9B1k/t62QnXRzteIyuufrMu\nfaOO94RffjxNvmUxZPLjfn2tX3/L4/qD23ZNmeVVU22Tkx8POjBgN7DZrm8GboztPxNYBRwFHA3c\nNUwGhRC5ID1krhhtsxAlo7J4aModFCiqRRZ9x/XAScChmHA3/w34KnAD8CJWduteiOnWXQLOA76W\nck1b8e8zsw3QNCePD/KchB94poUER2MTmqJpTl+OYf76xhgfP0RSjRrgkR8XURZDJj+G/nyhjPTD\nX0vlbX1xqWl+R4f9b+qw/2L7EUJUh17dVypNMhG2QI9oZL7IG5XFhaJxCaI3eQwEFAUSH6DgC8pT\nbdiL0UM+htFDPm73/xA4MpbuCLtvBb0HEAVdtrsd65S21zlZr9dp385cr99tgNHOnTszDYDOY0DM\n7OyskwE48/PzbN26tdD7id9LEdcP1z0d0FszAo9sh+Eki+9UKLO8aqrtqjOEkNuXgXzDHtdAwDzx\nMU/41zo7zcqBgBfY9W2sHAgY6iG/R3q3Vpf7rvtAwGzX70YdB83V8Z7wz4/zpst9N2MgYPJY0TR1\nMF4dBgKWpdOy99AfzdM0K4ZklfFICwkOxyY0W9PcnnaQ/znhF575cRFk8GPozxfKSJ/fteS39SMv\nP1al2etKswreKqPC1veKcD9pVWluKvLjplWa1VhVR/LyYw/nohe+E9f++YKPeWoaExNT4R9TAQQF\nXbcsO+5+sy59o473JIok8ND2EouLi4WGoSvz99tU23mhSrMQIhdMtAi1rApRZYp9+a0K4aBARcAR\n7VRGnjExMRUrlH2QVwx7XPKMutO0bt1ssgyfJBf9pJU8o6nIj5smz2jfJx+uB42TZzT3jU+zFQkh\nhBBClE1lKs3NJewmWlyuOIfdZ2VVpH3UJfmYJ5EnQc3s1FP/W8d7EkUSNNJ2U3XFdfBbVZorQzQ4\nIZSpNLf1XYgqoF4iIaqNfFi0UxlNczQwwRdN8rDH87mm9Fb+Ii2k7zrlftIOd335aXWRHzdb0ywf\nrgeN0zQLIUSVKVtWJYQQYjiGrTQvAPcCdwN32X1TwG3Ag8CtwOSQNkRHRkspgH3UJfmYJ5EnQeXt\nJGVVddT/1vGeKsQClSuPg0babqquuA5+O2yluQXMAscBJ9h92zBOegxwu90WhbAkXbMQQghQeVw4\nExNT6ilqOMPqOx4GXg38JLbvAeAkYC+wDvM6d2zivMya5ig+8/Kp+KFJHvZ4fteU1spPKqSFXAD2\nA88Bz2IK3CngC8CL7fG3A08mzpOmuQ9Nc/z5yGerQ4X8OJfyWJrmtH1jmEhWyw8MUS180TS3gK8D\n/wj8nt23FuOg2OXaYQxoljEhCkctVEJUn8LL4+ZiQr8aFFGjyYwOef7rgEeBX8YUsA8kjrfoUOPd\nsmUL09PTAExOTjIzM8Ps7CwQ6V7C7d7aoyoeD3I6PmrfoA5kfHyC/fufWPH88t7euXNn1++rjO35\n+Xm2bt1aan7C9YWFBSpI8g38VEwLFcAuzA+uxIpzUDM75rcS/cdV345LWy7vqULkUh4bgi5m0o4N\nmj7rtQa9frf0Qcq+LNdfAuZYXNxojlasvCqz/I6XlSqPYQfwAYyjrrP71rPSccF0B2UCaEG47LRe\nxePF2HTB3NycEzv94GOeqE4Xyb9gBg/FW6jimqiRxHZIyv1mXfabdq6g6yb39bIz2PXbn89oC2gd\nfPAhTn6HLn3DlS2X90R1/DjODgYsjxmo7Bo0/VyO13dve1DKLK+aapuc/HiYlubVwIHAIvBLwMnA\nh4HdwGbgMru8ccg8Cs/wsYXHxzxViKFbqC655DK7J0ik6LbdT9rZjOdkvV6W8/O6/mgsznxA2NX7\n9NMjbS2mPvTY5LG9fOcF2ou3WjWxhSqFipbHsxW2Hfr1GOPjh7B//xPZLZdYXjXVdl4MI4o+CviK\nXR8FPg9cghlAdAPwIjIOIOqawZFO4vz4ehWPF2Mz63MVxVOhAURxdgA/w7Q4zwKPYVqo5ugygGjl\n5EO9llVLW8z15a/+UxE/zq081kDAQdKPYl6GxzBjqc1yfHxNX5VpURw+DAR8GJixn3+HcVCAJ4A3\nYQYQncxKBxUVJ9mi5AM+5qkirAbG7XrYQnUfUQsVeNFCFdTMDsABjIyMMDKyqtBBRS59w5Ut+fsK\nhi6Pw8l33BI4tleU7XCg4LNty04hYcv8/TbVdl4MOxBQeEgYpk9vuSIDa1nZQnUrRt98A3A2UQuV\nyJXnCVuqFG9dlE0Uqcr3RnUhyqMs75A8o5Brhl1E0T51/5ZHRbp1h0HyjBzTtlotvfB6SFP8uHuc\ndTLua6o8o3N6lcF+4IM8o1DK6SqqOvFYkkKIqpGcalsIIYQ/eFtp1qQm/uKjLsnHPIk8CWpmJ8lo\naiNB2HgwjOZZmmbhL0EjbTdVV1wHv5WmudaMMjExpW5eIbwn7CUKK86jbZKXxUX1ugkhRNl4q2nu\nHfaGih93ZTMKhdNvLEkxHE3RQoI0zS7SShtZDk3xY2mai0g/ZtcVfq5saq9pFnkRhcIJdZJ5dPkK\nIVwyKp8VonIsEYWfW5QP1wCvKs2qzBWNKXiHHWzkoy7JxzyJPAlqZqdfzMvv4uK+5f/JkZFVmf4v\npWkW/hI0yHbkw03VFdfBb73SNEeVuTr3hJVJUjcphKga7fF09X8phmd6+lVlZ6FBjLJx40bGx9cA\noT9rBsGq4JWmuXNMZh81ycMe9yNP0kkWR1O0kCBNs7u0Y7S//MqPi6YJfgz/Hfh9qqET9tF2Xnkd\nY3x8XBXnAsjLj71qaRauiaJrhF284bomWBBZ+MAHLuCgg15QdjYaxFLKvjDSRntrVdynhejOS8rO\nQAOJT0YWsrSsfVb56ydeaZqFa5aWdc2Li/uWnTWueY7rzMP11avHy812CnXQSlWRK654jI9//NMO\nLAUObLi0kyfRYN/QbyH0abf6SWmaRX8EDbad9gIMkfZ5sePYhWHGf0nTPBxFVZo3AQ8ADwEX9Er8\nwQ/+V972ts2xPWoAd0c8HuzKGQXjFehw/emnf+Y4j72Zn58vOwt1JIMf/2dWrZp0kBVX328dfkft\nE6XEfSNZ2OY9+NqVH8rf+6Kv8tgtZX6PvttOvgwvLleghxnMX6bv1MFvi6g0Hwh8GuOoG4B3AC/v\ndsKVV17Jl78cl5p0egMT+TPI1NsjqW++ZUY9efLJJ0uzXVP69uNicfX91uF31O7TTz755HLlOFnY\nRtuLmaNxdMOVH8rfM+OZHycp83usmu0ofF2cfl98y/SdOvhtEZXmE4A9wALmG/5fwGm9T3ttAVkR\nwzHaYb1F2puv4lDWigH9WPjFKB/+8IdjETei/e3TdneOJ5sMcTcysqpjS/Ull1ym/wC/kB/XnPYX\n33QZR7cX4nga+W1viqg0Hw78ILb9iN3XORMHHMBBB7nQRYr+WOqwHt8Xf/NN02KtLGhD4q3TyfWk\njrqTM4fHP/axi1OP9foTUGzwjmTy49WrP8IzzzzmIDsLDmy4tOOKJWBzh/2depja/TiqcD+7vFw5\n/sFs/+IXz9DtPyBt2e736cdCevl7L1/OK02FyOjHO5xlqJ2FkuzWxXa8MWtl+Zv03cXFRT784Y+0\n+Vg8TVJHnfaC3O6jaT7cuQK+sLCQeq2svlYz31zmbcD/jG2/C/hUIs0ewuZKffSp72cP1UV+rI8+\n5iM/1kef6n9y8eMiRtz9EDgytn0k5u02zssKsCuEyA/5sRDVR34shOeMAt8DpoFVmGGiHg08EEJk\nQH4sRPWRHwtRAU4B/g+mOXx7yXkRQgyG/FiI6iM/FkIIIYQQQghRXbIEUv8ze/we4Lg+z3WdpwXg\nXuBu4C6HeToW+HvgGeADfZ5bRp4WKOc5vRPznd0LfAt4ZR/nlpWvBYp5Vq4oerKEBVY+nyngNuBB\n4FZgkBlVrgX2AvfF9nW77nbMPT4AnJyDrYswetK77eeUHGwdCcwB3wG+DZxr9xdxX51sXUT+93UQ\ncCdGTnA/cIndn/d9dbJzEfnfk28U7ccuf5udOBDz/d3k2PYk8EXgu5jf1YkObW/HPPP7gOuAFxRo\nO6//1H9vr/EQ8KdD2P4TzDO/B/gy8MKCbOfGgZguoGlgjHT91G8Af2XXTwTu6ONc13kCeBjzI8iT\nLHn6ZeDVwB/TXkEt8zl1yhOU95z+A5FjbKL439Ow+YJinpUrinyuIWnP53LgQ3b9AuDSAa77BswL\ncfxPttN1N2DubQxzr3voL0Rnmq0dwPtT0g5jax0wY9cPwXTDv5xi7quTrSLuC2C1XY5i/Of1FHNf\naXaKuidfcOHHLn+bnXg/8Hlgt912ZXsXcJZdH8WUBS5sTwP/gqkoA3wBE3PS5f9cP7bCYPF3YWKK\ng6mLbRrQ9ptj+b+0CNt5O3uWQOqnYn5QYN7wJzHOVVQQ9kHztDZ2fIR8yZKnHwP/aI/3e67rPIWU\n8Zz+HvipXb8TOKKPc8vIV0jez8oVriZLSD6fuI/uAk4f4JrfAJLzzna67mnA9Zh7XMDc8wlkJ80W\npH/vw9h6jGhO3p9hWlkOp5j76mQL8r8vgKfschWmkrePYu4rzQ4Uc0++4MKPXf420zgC0yB2NdF3\n6cL2CzEVumvt9hKmLHBhe7+9zmpMZX018KMCbQ/7n3oisB4YJ+pV/CzZ/t/TbN8GPG/X4+Vubrbz\nrjRnCaTeKc1hGc51nScw8f2+jqks/l4O+cmapyLOLfK6Pjyns4l6DIp6TsPmC4p5Vq4o8rmGpD2f\ntZiuOOxybcp5g9DpuofRHporr/t8H6br8Bqibsu8bE1jWl7upPj7Cm2FPShF3NcBmErXXqJu/iLu\nK80OFPtdlY0LP44zjbvfZsgngQ8SVaJwZPsoTAPTZ4B/xsTJ/iVHtp8APgF8H1NZfhJTkXT5P9ev\nreT+H+aQBzAt/WG5m5vtvCvNrYzpXLayDZun12Oc/RTg9zFvkK7ylPe5RV73dZT7nDZinCTU5hX1\nnPq9djJfUMyzckWRzzWk1/MJg9XnTa/rDmvzKkyBOgM8iinc8rJ1CPAl4DxgMeVaed7XIRi95nmY\nFsSi7ut5e80jgF/D+FLyWnncV9LOLMV+Vz7gMs8uf5shbwEex+iZO5XtRdkeBX4VuNIufw5sc2T7\npcBWzEvKYZhn/y5Htjtdqwz/+EPgFxhNd67kXWnOEkg9meYImybLuS7z9EO7/iO7/DHwFfLpihvm\nXst8Tt141C7LeE6vxLzNn0rUXVPUcxo2X1DMs3JFkc81JO357MXIuMB0qT2ek61O1+32nzAojxMV\nIlcTfe/D2hrDVEo+B9xo9xV1X6Gtv4jZKuq+Qn4K/CVmwE6R31do59UUf09l48KPwe1vM85rMf+7\nD2O65d9o8+DC9iP28w92+4uYyvNjDmy/Gvg74CcYWciXMeNrXNgO6ecZh3W/IxL7h8nDFows552x\nfa5s902WQOrxQXevIereKyoI+zB5Wo3Ru4DpXvkW+Yzo7edeL6J90F2Zz6lTnsp8Ti/C6JNeM8C5\nZeSrqGfliqInS+j0fC4naq3fxmADAcHkOzloJe264cCRVZgWx+/Rfw9Z0tb62Pr5RK0gw9gawejw\nPpnYX8R9dbJVxH0dSiSJOBj4W+DXyf++OtlZF0uT1z35hItJT1z+NrtxElH0DFe2/xY4xq5fZO26\nsP0qTKSSg+01dmF661z+zw1i606MxniE7AMB02xvwsirDk2kK8J2bqQFUj/HfkI+bY/fg3kD63Zu\nmXl6CeZBz2N+iC7ztA6jOfspppXy+5iulk7nlpmnMp/T1Zi36jA01F09zi07X0U+K1cU+VyPIv35\nTGF0zsOEnLse03P0C8zv+Hd7XPdCzD0+APzHIW2dhalA3Iv5j7mRdl32oLZej5EXzBP91jZRzH2l\n2TqloPt6BUYPOm+v/UG7P+/76mSniHvyjaInPXH52+zGSUTRM1zZfhWmpTke+syV7Q8RhZzbhWnt\nd/U/N+h/ahj2bQ8m/O8gts/ChI37v0S/tysLsi2EEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBC\nCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgjROP4/imEdW3J766EAAAAASUVORK5CYII=\n",
       "metadata": {}
      }
     ],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "df_Es = pd.DataFrame(data=[[D_m_eval_E_ana, D_u_eval_E_ana, D_mu_eval_E_ana],\n                           [D_m_eval_E_samp,D_u_eval_E_samp, D_mu_eval_E_samp]],\n                     columns=['D_m', 'D_u', 'D_mu'],\n                     index=['Analytic', 'Sampling'])\ndf_Es",
     "prompt_number": 345,
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 345,
       "text": "                         D_m               D_u              D_mu\nAnalytic  0.0601137900000000  240.000000000000  72.0420796530000\nSampling  0.0598427529438614  239.854544301056  72.0536470558759",
       "html": "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>D_m</th>\n      <th>D_u</th>\n      <th>D_mu</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>Analytic</th>\n      <td> 0.0601137900000000</td>\n      <td> 240.000000000000</td>\n      <td> 72.0420796530000</td>\n    </tr>\n    <tr>\n      <th>Sampling</th>\n      <td> 0.0598427529438614</td>\n      <td> 239.854544301056</td>\n      <td> 72.0536470558759</td>\n    </tr>\n  </tbody>\n</table>\n</div>",
       "metadata": {}
      }
     ],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "markdown",
     "source": "Uncertainty of the mode"
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "tmp1 = pd.DataFrame(['%f' %f.evalf() for f in [D_m_up['uval'], D_u_up['uval'], D_mu_up['uval']] ],\n                    index=['D_m', 'D_u', 'D_mu'], \n                     columns=['Analytic'])\ntmp2 = pd.DataFrame(df_samples.std(), columns=['Sampling'])\ndf_Us = pd.concat([tmp1, tmp2], axis=1).T\ndf_Us",
     "prompt_number": 346,
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 346,
       "text": "                 D_m        D_u       D_mu\nAnalytic    0.006870  50.519798  15.153941\nSampling  0.02629003   138.9449   41.68986",
       "html": "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>D_m</th>\n      <th>D_u</th>\n      <th>D_mu</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>Analytic</th>\n      <td>   0.006870</td>\n      <td> 50.519798</td>\n      <td> 15.153941</td>\n    </tr>\n    <tr>\n      <th>Sampling</th>\n      <td> 0.02629003</td>\n      <td>  138.9449</td>\n      <td>  41.68986</td>\n    </tr>\n  </tbody>\n</table>\n</div>",
       "metadata": {}
      }
     ],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "markdown",
     "source": "Partial derivatives"
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "D_m_up['df'].ix[['partial_term_values', 'product_term_values']]",
     "prompt_number": 347,
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 347,
       "text": "                             L         g        km       thm        cm\npartial_term_values  -0.000601  0.000000  0.011275  0.293238  0.010930\nproduct_term_values   0.000036  0.000000  0.000001  0.000009  0.000001",
       "html": "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>L</th>\n      <th>g</th>\n      <th>km</th>\n      <th>thm</th>\n      <th>cm</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>partial_term_values</th>\n      <td> -0.000601</td>\n      <td> 0.000000</td>\n      <td> 0.011275</td>\n      <td> 0.293238</td>\n      <td> 0.010930</td>\n    </tr>\n    <tr>\n      <th>product_term_values</th>\n      <td>  0.000036</td>\n      <td> 0.000000</td>\n      <td> 0.000001</td>\n      <td> 0.000009</td>\n      <td> 0.000001</td>\n    </tr>\n  </tbody>\n</table>\n</div>",
       "metadata": {}
      }
     ],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "D_u_up['df'].ix[['partial_term_values', 'product_term_values']]",
     "prompt_number": 348,
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 348,
       "text": "                              L            g         ku          thu         cu\npartial_term_values    2.400000  -400.000000  80.000000  1600.000000  75.000000\nproduct_term_values  576.000000  1600.000000  64.000000   256.000000  56.250000",
       "html": "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>L</th>\n      <th>g</th>\n      <th>ku</th>\n      <th>thu</th>\n      <th>cu</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>partial_term_values</th>\n      <td>   2.400000</td>\n      <td> -400.000000</td>\n      <td> 80.000000</td>\n      <td> 1600.000000</td>\n      <td> 75.000000</td>\n    </tr>\n    <tr>\n      <th>product_term_values</th>\n      <td> 576.000000</td>\n      <td> 1600.000000</td>\n      <td> 64.000000</td>\n      <td>  256.000000</td>\n      <td> 56.250000</td>\n    </tr>\n  </tbody>\n</table>\n</div>",
       "metadata": {}
      }
     ],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "D_mu_up['df'].ix[['partial_term_values', 'product_term_values']]",
     "prompt_number": 349,
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 349,
       "text": "                             L            g        km       thm        cm  \\\npartial_term_values   0.719579  -120.000000  0.007893  0.205267  0.007651   \nproduct_term_values  51.779423   144.000000  0.000001  0.000004  0.000001   \n\n                           fm         ku         thu         cu          fu  \npartial_term_values  0.060114  24.000000  480.000000  22.500000  240.000000  \nproduct_term_values  0.000000   5.760000   23.040000   5.062500    0.000000  ",
       "html": "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>L</th>\n      <th>g</th>\n      <th>km</th>\n      <th>thm</th>\n      <th>cm</th>\n      <th>fm</th>\n      <th>ku</th>\n      <th>thu</th>\n      <th>cu</th>\n      <th>fu</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>partial_term_values</th>\n      <td>  0.719579</td>\n      <td> -120.000000</td>\n      <td> 0.007893</td>\n      <td> 0.205267</td>\n      <td> 0.007651</td>\n      <td> 0.060114</td>\n      <td> 24.000000</td>\n      <td> 480.000000</td>\n      <td> 22.500000</td>\n      <td> 240.000000</td>\n    </tr>\n    <tr>\n      <th>product_term_values</th>\n      <td> 51.779423</td>\n      <td>  144.000000</td>\n      <td> 0.000001</td>\n      <td> 0.000004</td>\n      <td> 0.000001</td>\n      <td> 0.000000</td>\n      <td>  5.760000</td>\n      <td>  23.040000</td>\n      <td>  5.062500</td>\n      <td>   0.000000</td>\n    </tr>\n  </tbody>\n</table>\n</div>",
       "metadata": {}
      }
     ],
     "language": "python",
     "collapsed": false
    },
    {
     "metadata": {
      "run_control": {
       "breakpoint": false
      }
     },
     "cell_type": "code",
     "input": "#for v_it, v in enumerate(D_m_variables): pprint(\"%s : %f:\" %(v,D_m_up['partial_term_values'][v_it]) )",
     "prompt_number": 350,
     "outputs": [],
     "language": "python",
     "collapsed": false
    }
   ],
   "metadata": {}
  }
 ],
 "metadata": {
  "gist_id": "8442130",
  "name": "",
  "css": [
   ""
  ]
 },
 "nbformat": 3
}