{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##Getting Started:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is an IPython notebook.  IPython notebooks are open-source programming tools that run in a web browser, and offer a lot of flexibility and power for scientific computing tasks.  Because they are portable and self-contained, they are ideal for education and [reproducable research](http://reproducibleresearch.net/). Notebooks are made of a single column of 'cells' that usually consist of markdown or code.  Markdown cells like this one and the one above support basic [markdown formatting](https://en.wikipedia.org/wiki/Markdown), as well as $\\LaTeX$ statements.  The cell below is a code cell, which runs Python code interactively if you use your personal computer or a cloud service like [Domino](https://app.dominodatalab.com/jakealbrecht/ChE_problemSolving).  The Python commands can be run over and over, or the background processes can be stopped or restarted without affecting this notebook.  Did I mention flexibility? Cells can be reordered, and the whole notebook can be saved as a PDF document, or can be hosted as a static webpage on sites like [nbviewer](http://nbviewer.ipython.org/github/chepyle/ChE_ProblemSolving/blob/master/ChE%20Problem%20Solving.ipynb).\n",
    "\n",
    "Scientific applications of the Python progamming language have increased dramatically over the past few years.  Its simplicity, versatility, and power has made it a first-choice language for many scientists and engineers.  Just enter \"`import antigravity`\" in a code cell or at a Python prompt to see the zeal Python inspires.\n",
    "\n",
    "To get started, the main function libraries in Python for working with numbers are `numpy` and `pandas`. Standard numerical methods are often in the `scipy` library.  The `sklearn` library has fitting and regression, and `matplotlib` is used for matlab-style plots.  To begin our Python session, these libraries and some specific functions are imported in the cell below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.optimize import fsolve\n",
    "from scipy.optimize import curve_fit\n",
    "from scipy.integrate import odeint\n",
    "\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "# turn on inline plotting for this notebook:\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reference"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This hastily written notebook attempts to replicate the workflow and results of M. Shacham's \"Using Templates for Demonstrating Good Programming Practices\" submitted to the 2015 Annual AIChE meeting: [ftp site](ftp://ftp.bgu.ac.il/shacham/templates/)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#Basic Numerical Problem Solving with Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##One Nonlinear (Implicit) Algebraic Equation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve the followng system of equations:\n",
    "\n",
    "$$P=\\frac{RT}{\\left(V-b\\right)}-\\frac{\\alpha a}{V\\left(V+b\\right)}$$\n",
    "\n",
    "where:\n",
    "\n",
    "$$a=0.42747\\left(\\frac{R^2T_c^2}{P_c}\\right)$$\n",
    "\n",
    "$$b=0.08664\\left(\\frac{RT_c}{P_c}\\right)$$\n",
    "\n",
    "$$m=0.48508+1.55171\\bar{\\omega}-0.1561\\bar{\\omega}^2$$\n",
    "\n",
    "$$\\alpha=\\left[1+m\\left(1-\\sqrt{T/T_c}\\right)\\right]$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a Python function for the first equation, set equal to zero:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def f(V):\n",
    "    return R*T/(V-b)-alph*a/(V*(V+b))-P"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Enter some constants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "R=0.08206 #Gas constant (L-atm/gmol-K)\n",
    "P_c=72.9\n",
    "T_c=304.2\n",
    "w_bar=0.225\n",
    "T=300.\n",
    "P=100.\n",
    "a=0.42747*(R**2.0*T_c**2.0/P_c)\n",
    "b=0.08664*(R*T_c/P_c)\n",
    "m=0.48508+1.55171*w_bar-0.1561*w_bar**2\n",
    "alph=(1+m*(1-np.sqrt(T/T_c)))**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print out values, and solve equation with initial guess"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "initial f(V0)= -19.366385\n",
      "a = 3.653924\n",
      "b = 0.029668\n",
      "m = 0.826312\n",
      "alpha = 1.011481\n",
      "V= 0.064113\n"
     ]
    }
   ],
   "source": [
    "print('initial f(V0)= %f' %f(0.07))\n",
    "print('a = %f' %a)\n",
    "print('b = %f' %b)\n",
    "print('m = %f' %m)\n",
    "print('alpha = %f' %(alph))\n",
    "print('V= %f' %(fsolve(f,0.07))) # Print out solution result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final result of the optimization is $V = 0.0641$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A System of Nonlinear Algebraic Equations (NLE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following is a system of ultrafiltration equations from Cutlip and Shackam, 2013\n",
    "\n",
    "$$Q_1c_1=Q_0c_0$$\n",
    "\n",
    "$$Q_1=Q_0-j_1A$$\n",
    "\n",
    "$$j_1=k\\ln\\left(c_g/c_1\\right)$$\n",
    "\n",
    "$$Q_2c_2=Q_1c_1$$\n",
    "\n",
    "$$j_2=k\\ln\\left(c_g/c_2\\right)$$\n",
    "\n",
    "$$Q_3c_3=Q_2c_2$$\n",
    "\n",
    "$$Q_3=Q_2-j_3A$$\n",
    "\n",
    "$$j_3=k\\ln\\left(c_g/c_3\\right)$$\n",
    "\n",
    "$$Q_0=1/(60\\times1000)$$\n",
    "\n",
    "$$k=3.5e-6$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cg=300\n",
    "c0=10\n",
    "A=0.9\n",
    "k=3.5e-6\n",
    "Q0=1/(60.0*1000.0)\n",
    "Const=[cg,c0,A,k,Q0]\n",
    "xguess=[40.0,70.0,150.0]\n",
    "\n",
    "def Evarfun(x,Const):\n",
    "    c1,c2,c3 = x # unpack concentrations\n",
    "    cg,c0,A,k,Q0 = Const # unpack constants\n",
    "    j1=k*np.log(cg/c1)\n",
    "    Q1=c0*Q0/c1\n",
    "    j2=k*np.log(cg/c2)\n",
    "    Q2=c1*Q1/c2\n",
    "    j3=k*np.log(cg/c3)\n",
    "    Q3=c2*Q2/c3\n",
    "    Evar=[j1,j2,j3,Q1,Q2,Q3]\n",
    "    return Evar\n",
    "\n",
    "def MNLEfun(x,Const):\n",
    "    c1,c2,c3=x# unpack concentrations\n",
    "    cg,c0,A,k,Q0 = Const # unpack constants\n",
    "    j1=k*np.log(cg/c1)\n",
    "    Q1=c0*Q0/c1\n",
    "    j2=k*np.log(cg/c2)\n",
    "    Q2=c1*Q1/c2\n",
    "    j3=k*np.log(cg/c3)\n",
    "    Q3=c2*Q2/c3\n",
    "    fx=[Q0-Q1-(j1*A),Q1-Q2-(j2*A),Q2-Q3-(j3*A)]\n",
    "    return fx\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check function values at initial guesses:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "j1: 7.05e-06\n",
      "j2: 5.09e-06\n",
      "j3: 2.43e-06\n",
      "Q1: 4.17e-06\n",
      "Q2: 2.38e-06\n",
      "Q3: 1.11e-06\n"
     ]
    }
   ],
   "source": [
    "evars_initial=Evarfun(xguess,Const)\n",
    "var_names=['j1: ','j2: ','j3: ','Q1: ','Q2: ','Q3: ']\n",
    "for var in zip(var_names,evars_initial):\n",
    "    print(var[0]+'{:.2e}'.format(var[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK, it is time to solve, note that the solver tolerance `xtol` must be reduced due to the very small function values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "c1 = 20.35 \n",
      "c2 = 56.67 \n",
      "c3 = 163.12\n",
      "function values = [0.0, 0.0, -4.2351647362715017e-22]\n",
      "j1: 9.42e-06\n",
      "j2: 5.83e-06\n",
      "j3: 2.13e-06\n",
      "Q1: 8.19e-06\n",
      "Q2: 2.94e-06\n",
      "Q3: 1.02e-06\n"
     ]
    }
   ],
   "source": [
    "c=fsolve(MNLEfun,xguess,args=Const,xtol=1e-13)\n",
    "print('c1 = %.2f \\nc2 = %.2f \\nc3 = %.2f' %(c[0],c[1],c[2]))\n",
    "print('function values = ' + str(MNLEfun(c,Const)))\n",
    "evars_final=Evarfun(c,Const)\n",
    "for var in zip(var_names,evars_final):\n",
    "    print(var[0]+'{:.2e}'.format(var[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are able to solve the equations for the values of $c_1,c_2$ and $c_3$ , note that the objective function values are all *very nearly* zero, indicating that we have accurately fit our model. The $c$ values as well as the $Q$ and $j$ values now fully describe the unit operation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Systems of First-Order Ordinary Differential Equations (ODEs) - Initial Value Problems (IVPs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Multicomponent gas diffusion is described by the following system of equations:\n",
    "\n",
    "$$ \\frac{dC_A}{dz}=\\frac{x_AN_B-x_BN_A}{D_{AB}}+\\frac{x_AN_C-x_CN_A}{D_{AC}}$$\n",
    "\n",
    "$$ \\frac{dC_B}{dz}=\\frac{x_BN_A-x_AN_B}{D_{AB}}+\\frac{x_BN_C-x_CN_B}{D_{BC}}$$\n",
    "\n",
    "$$ \\frac{dC_C}{dz}=\\frac{x_CN_A-x_AN_C}{D_{AC}}+\\frac{x_CN_B-x_BN_C}{D_{BC}}$$\n",
    "\n",
    "$$x_A=C_A/C_T$$\n",
    "\n",
    "$$x_B=C_B/C_T$$\n",
    "\n",
    "$$x_C=C_C/C_T$$\n",
    "\n",
    "$$N_A=2.115\\times10^{-4}$$\n",
    "\n",
    "$$N_B=-4.143\\times10^{-4}$$\n",
    "\n",
    "$$N_C=0$$\n",
    "\n",
    "$$D_{AB}=1.47\\times10^{-4}$$\n",
    "\n",
    "$$D_{AC}=1.075\\times10^{-4}$$\n",
    "\n",
    "$$D_{BC}=1.245\\times10^{-4}$$\n",
    "\n",
    "$$C_T=0.2/(0.082057\\times328)$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the function that calculates the derivatives at any point in time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ODEfun(Yfunvec,z,NA,NB,NC,DAB,DBC,DAC,CT):\n",
    "    CA,CB,CC=Yfunvec\n",
    "    xA=CA/CT\n",
    "    xB=CB/CT\n",
    "    xC=CC/CT\n",
    "    dCAdz=(xA*NB-(xB*NA))/DAB+(xA*NC-(xC*NA))/DAC\n",
    "    dCBdz=(xB*NA-(xA*NB))/DAB+(xB*NC-(xC*NB))/DBC\n",
    "    dCCdz=(xC*NA-(xA*NC))/DAC+(xC*NB-(xB*NC))/DBC\n",
    "    return np.array([dCAdz,dCBdz,dCCdz])\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, define the constants and run the integration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "NA=0.00002115\n",
    "NB=-0.0004143\n",
    "NC=0\n",
    "DAB=0.000147\n",
    "DBC=0.0001245\n",
    "DAC=0.0001075\n",
    "CT=0.2/(0.082057*328)\n",
    "Const=NA,NB,NC,DAB,DBC,DAC,CT\n",
    "\n",
    "y0=[0.00022290,0,0.007208]\n",
    "tspan=np.linspace(0,0.001,40)\n",
    "ty=odeint(ODEfun,y0,tspan,args=Const)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x7fe928e66ed0>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAEPCAYAAAANl7AYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4HNWd7//317KNN8A2tmVLlmyMN3nfLdsxaEImMZ4B\nMskDhMkvCdzMmAmXzHMnkAEmN8EkTwbIMEkG+N2EO4b8uDNMnOXOOHBDAgwXBWO874u8Ad5kS17A\neLe1fH9/VHWrJaRWy1J3taTP63nq6a7qc7pOV4w+OVWnTpm7IyIiErUuUTdAREQEFEgiIpIlFEgi\nIpIVFEgiIpIVFEgiIpIVFEgiIpIVIg0kM1tgZjvNbI+ZPdhEmafCzzeb2dTm6prZLDNbY2YbzWyt\nmc3MxG8REZHWiSyQzCwHeAZYAIwD7jSzogZlFgIj3X0UsAj4SQp1fwB8292nAt8J10VEJMtF2UOa\nBex1933uXgUsBW5tUOYW4AUAd18N9DWzwc3UPQJcHb7vC5Sn92eIiEhb6BrhvvOBgwnrh4DZKZTJ\nB/KS1H0IeNvMniQI3Dlt2GYREUmTKHtIqc5ZZC383ueAv3b3QuBvgOdbWF9ERCIQZQ+pHChIWC8g\n6OkkKzM0LNMtSd1Z7v6p8P2vgSUNd2xmmsBPROQyuHtLOwkpi7KHtA4YZWbDzaw7cAfwUoMyLwFf\nBjCzYuCku1c2U3evmd0Qvv8ksLuxnbu7FnceeeSRyNuQLYuOhY6FjkXyJd0i6yG5e7WZ3Qe8CuQA\nz7l7mZndE37+rLu/YmYLzWwvcBa4O1nd8KsXAf+vmV0BnA/XRUQky0V5yg53/x3wuwbbnm2wfl+q\ndcPt6/j44AgREclymqmhkyspKYm6CVlDx6KOjkUdHYvMsUycF8w2Zuad8XeLiLSGmeFpHNQQ6Sk7\nEZH2wCxtf4OzVhT/p12BJCKSgs50ViWqANY1JBERyQoKJBERyQoKJBERyQoKJBERyQoa1CAi0kG8\n/PLLrFmzhry8PHr27EnPnj1ZuXIljz/+OD169Ii6ec1SIImItHO1tbUsWrSIsWPH8r3vfS++fdmy\nZWzdurVdhBEokERE2r1HH30UgAceeKDe9jlz5rB7d6PzS2clzdQgItKMcIaCqJvRqBMnTlBQUMCu\nXbsoKCio95m7c/78eXr16tWi72zq92qmBhGRLNdW95FeTuYtX76cwsLCj4URBAHS0jCKkkbZiYi0\nknvbLJcjJyeH/v37N/rZiy++WG+9qqqKu++++/J2lAGdN5AuXoy6BSIirXbjjTdy4sQJ9u/fH99W\nW1vLkiVLWLhwYb2yO3bsoLy8PNNNTFnnPWXXvz9Mngxz5tQt+flRt0pEpEV69erFyy+/zPe//32K\niorivaXPfe5z9OvXL17u3LlzDB8+nK5ds/fPfqSDGsxsAfBjgqe+LnH3Jxop8xRwE3AOuMvdNyar\na2ZLgTFh9b4Ejz2f2uA73c+cgbVr4Z13YOXKYOndu35ATZkC3bun58eLSLuRzYMaUvXmm29y6tQp\nfvSjH7Fs2TL69u3bZNlON6jBzHKAZ4BPAeXAWjN7KeFR5JjZQmCku48ys9nAT4DiZHXd/QsJ9Z8E\nTjbagN69oaQkWCA4gbtnTxBM77wDzz8Pe/cGoRQLqOJi9aJEpN3ZtWsX8+fPp2vXrrz55ptUVFQk\nDaSoRHkNaRaw1933uXsVsBS4tUGZW4AXANx9NdDXzAanUteC+dNvB36eUmvMYPRo+MpX4NlnYfNm\nqKiA730P+vWDn/0sOMVXWAh33AE/+hGsWqVrUSKS1V599dX4fUrHjh1jz549/PrXv464VY2L8mRi\nPnAwYf0QMDuFMvlAXgp15wOV7v7uZbfwyivhk58MFgh6UXv3BkG0ciX8y7/Arl0wcWLQe4otw4a1\n3ThQEZFW+MxnPsNnPvMZAAYOHMhvf/vbiFvUtCgDKdUTspf7l/1O4N+a+nDx4sXx9yUlJZTETt0l\nbYnBqFHB8qUvBdvOnoX164OQ+uUv4W/+Jgiu4mKYPTt4nTEjCDcRkXaktLSU0tLSjO0vskENZlYM\nLHb3BeH6w0Bt4sAGM/spUOruS8P1ncANwLXJ6ppZV4Je0zR3P9zIvtM3U4M7HDwY9KBWrw6WTZvg\nuuvqAmr2bCgqgpyc9LRBRNpURxjU0BJRDWqIMpC6AruAG4HDwBrgzkYGNdzn7gvDAPuxuxc3Vzcc\ngfegu/9RE/vO7NRBly7Bli1BOK1aFbxWVgY9p9mz65bBgzPXJhFJmQKp3vaOF0gAZnYTdUO3n3P3\nx8zsHgB3fzYs8wywADgL3O3uG5qqm/C9PwNWuvv/bGK/0c9ld+IErFlT14taswb69IFZs+oCavp0\naEfTfoh0VAqkets7ZiBFJSsCqaHYgIlYOK1eDdu2BSP/YiE1a5ZO9YlEQIFUb7sCqS1lZSA15uLF\n4PrTmjV1y5EjMG1aXUDNmgVDh2pUn0gaKZDqbVcgtaV2E0iN+eADWLeuLqBWr4YuXerCaebMYEmY\nMkREWkeBVG+7AqkttetAasgdDhwIpkGKhdSGDcEAiZkz64JqyhTo2TPq1oq0SwqketsVSG2pQwVS\nY2pqYOfOIJxiQVVWFlyPivWgZs6ECRMgiydaFMkWCqR62xVIbanDB1JjLlwIpkNau7ZuOXAgmA4p\nFlAzZgQ3/XbpvE8lEWlMewmkl19+mTVr1pCXl0fPnj3p2bMnK1eu5PHHH6dHjx4pf48CKYM6ZSA1\n5tSp4PReLKDWrQuuUU2bFoRTLKSGD9egCenUsj2QamtrWbRoEWPHjuWBBx6Ib1+2bBlPP/00b7zx\nRou+r9PN9i1Z4Kqr6s94DnD8eDAV0tq18OKL8N/+WzDab8aM+kt+vkJKJEvEJk9NDCOAOXPmsHv3\n7iiadFnUQ5LmHT4c9J5iy9q1wbWnhiGVmxt1S0XSIpt7SCdOnKCgoIBdu3ZRUFBQ7zN35/z58/Rq\n4Q32OmWXQQqkVorN15cYUuvWBbNKzJgRzDARex00KOrWirRac4Fkj7bN32h/pOV/l5YtW8ZDDz3E\nzp0726QNoFN20p6YBc+FKiyEz30u2OYO770XnO5bvx7+8R+D1z596sIptiikpIO5nCBpKzk5OfHH\nljf04osv8sUvfpEPP/yQRYsWMWfOHEaMGMErr7zCN7/5TUaNGpXh1ianQJK2YRbMaH7ddXD77cG2\nxJBatw6efDIYRNGnT/2Amj5dp/tELtONN97IAw88wP79+xk2bBgQDHJ4/vnn+fznPw9Av379uOqq\nq/jGN74BwIoVKzh9+nRkbW6KTtlJZjXsSa1fH4RUr1514TRtWrDk5UXdWhEgu68hAezevZsnn3yS\noqKieG/p5ptvjr93d2677Ta+/vWvs3LlSoYPH84XvvCFJr9P15AySIGUZdzh/ffrAmrjxuC1W7e6\ncJo2LQirggKN7pOMy/ZAas62bdt47bXX+MY3vkFtbS1Tpkxhy5YtTZbXNSTpvMxgxIhgue22YFts\n4MSGDUE4LVkCX/saVFfXD6mpU4PThLqZV6RJy5cvp7i4GIDKysqsDVcFkmSnxIETn/1s3fbDh4Me\n1IYN8POfwze/CSdPBnP1JYbU2LGaFkkE2LRpE7/4xS/o06cPBw8e5O233+Y3v/lN1M1qlE7ZSft3\n4kRdSMWW8vJgrr6pU+uWiRM1waxclvZ+yq6lOuU1pPBR47Gnvi5x9ycaKfMUcBNwDrjL3Tc2V9fM\nvg7cC9QAv3X3Bxt8pwKpozt9Opi7LxZUGzfCrl3B6b1YL2rKlGDRozqkGQqkets7XiCZWQ6wC/gU\nUA6sBe5097KEMguB+9x9oZnNBv7J3YuT1TWzPwL+Dljo7lVmNtDdjzXYtwKpM7p4EbZvrwuoTZtg\nyxYYMKAuoGK9KU2NJAkUSPW2d8hBDbOAve6+D8DMlgK3AmUJZW4BXgBw99Vm1tfMBgPXJqn7NeAx\nd68K69ULI+nErrii7jpTTE1N8Oj4TZuCkHrmmeDVva4HFQuq0aN1XUokjaL8rysfOJiwfgiYnUKZ\nfCAvSd1RwPVm9vfABeABd1/Xhu2WjiQnB8aMCZY77gi2uQePit+0KVh+8xtYvDgYUDF+fP2Qmjgx\nuNFXRFotykBKtf/b0u5hV6BfeGpvJvBLYETDQosXL46/LykpoSRxxmvp3MyCm3Lz8mDhwrrtp08H\np/hivannn4cdO2Do0CCgJk+uC6u8PJ3yk3avtLSU0tLSjO0vymtIxcBid18Qrj8M1DYYnPBToNTd\nl4brO4EbCE7ZNVrXzH4HPO7ufwg/2wvMdvcTCd+ra0jSNqqrg8ESmzYFgyhivara2rqQii1FRdC9\ne9Qtlsuga0j1tnfIQQ1dCQYm3AgcBtaQfFBDMfDjsOfTZF0zuwfIc/dHzGw08J/uXthg3wokSR93\nqKioC6lYUO3bF5waTAypyZODQRWS1RRI9bZ3vEACMLObqBu6/Zy7PxYGCu7+bFjmGWABcBa42903\nNFU33N4NeB6YAlwC7nf30gb7VSBJ5p07F4zyi4VUbOnTJwimSZPqQkoDKLKKAqne9o4ZSFFRIEnW\ncIf9+4Ng2rKl7vXQoeAUXyygJk0KlmuuibrFnZICqd52BVJbUiBJ1jtzpq43lRhUffrUhVMsqMaM\nCSailbRpL4H08ssvs2bNGvLy8ujZsyc9e/Zk5cqVPP744/To0SPl71EgZZACSdoldzhwIAimxJDa\nvz84xRcLqokTg9chQzTSr41keyDV1tayaNEixo4dywMPPBDfvmzZMp5++mneeOONFn1fZ7wxVkRa\nwgyGDQuWm2+u237uXDD8fOvWIKB+//vg1b0unCZODJYJE6B37+h+g6TFo48+ClAvjADmzJnD7t27\no2jSZVEPSaQjcofKyrre1NatwbJzZ3CPVCygYsvIkRpEkUQ295BOnDhBQUEBu3btoqCgoN5n7s75\n8+fp1atXi75TPSQRaTtmMHhwsHz603Xbq6uDqZJiIfXii8HrkSPBIztivahYUOkG39S01TG6jNBb\nvnw5hYWFHwsjCAKkpWEUJQWSSGfStWsQPGPHwu23122PDaLYti0IqN//Pnitrq4fUBMmBEvfvtH9\nhmwUYe8pJycn/qjyhl588UW++MUvUlVVxZIlS+jbty9XXHEFu3fv5qabbmLy5MkZbm1yOmUnIk07\nerTudN/WrUFobd8eBFIsnGJLURG0o/833hLZfMru3LlzTJ06lddee41hw4YBwSCH559/ns9//vP0\n69ePr371q3z7299m+PDhAHznO9/h0UcfxZro2WmUXQYpkERaobY2GO0X601t2xYsu3cH8/rFAmr8\n+OB19Oh2P2VSNgcSwO7du3nyyScpKiqK95Zuvvlm+vfvz69+9SveffddHnrooXj58vJy8vPzm/w+\nBVIGKZBE0qCqKrg+tW1b3em/bduCYekjRtQPqfHjg4cltpOBFNkeSMnce++93H333cycOTO+7eTJ\nk/RNcto1qwc1mNlwYKS7/6eZ9QK6uvupdDVKRNqhbt2C03ZFRXDbbXXbL1wIRvdt3x70qH72s+B9\nRUXQe0oMqfHj4dproUuX6H5HB1NQUEBVVVV8fdWqVQwbNixpIEWl2R6SmS0C/hLo7+7XhROW/sTd\nb8xEA9NBPSSRLHD2LJSV1fWoYsvx48Ggi1hAjRsXvA4fHllQtece0sWLF3n88ccZMWIE/fv3Z+TI\nkYwZMyZpnaw9ZWdmmwme7rrK3aeG27a6+8R0NSrdFEgiWezUqSCoEkNq+3b44IOg9xULqAwGVXsO\npMuRzafsLrr7xdhojPDRD53nfxkRyayrroLZs4Ml0alTwYwU27cHr6WlwfsTJ4L5/BKDaty44LpV\nTk4kP0EuTyo9pH8ATgJfBu4D7gV2uPu30t+89FAPSaQDOXWq7hrVjh11oVVZGVyjigVUbBk5ssWT\n0aqHVG97pKfscoCvArHbvV8FlrTnv+gKJJFO4OzZIKhiIVVWFrwePBj0nmKn/2IDMcaMgZ49G/0q\nBVK97Rr23ZYUSCKd2IULwWPny8rqQqqsDN59N5gqqWFQFRVhffsqkIgwkMxsa5J67u6TWr1zswXU\nPfV1ibs/0UiZp4CbgHPAXe6+MVldM1sM/AVwLPyKh9399w2+U4EkIvVVVcF779UPqbIy2LkTO3tW\ngUS0gTQ8WUV339eqHQenAncBnwLKgbXAne5ellBmIXCfuy80s9nAP7l7cbK6ZvYIcNrdf5hk3wok\nEUlNbS2Wk6NAIv2B1OQou8TAMbPBwGygFljr7hVtsO9ZwN7YfsxsKXArUJZQ5hbghbA9q82sb9iW\na5upq+mJRaRthEPKm5r3TdpOs8O+zewvgO8Ab4abnjGz77r7c63cdz5wMGH9EEHoNVcmH8hrpu7X\nzezLwDrgfnc/2cq2ikgn1uLe0blzwdx+ZWXB9aqdO4Nl92645pq6GdcTFz3qI6X7kP4WmOruJwDM\n7BpgJdDaQEr1f+GW/i/0E+C74fvvAf9IMEqwnsWLF8ffl5SUUFJS0sLdiIg0oVcvmDIlWBLFJqaN\nBdTWrfDLXwahdfZsMNIvFlCx96NGQY8ekfyM0tJSSktLM7a/VIZ9vwP8kbtfDNevAN5097mt2rFZ\nMbDY3ReE6w8DtYkDG8zsp0Cpuy8N13cCNxCcsktaN9w+HHi54awSuoYkIlnnww+DYErsUe3cCe+/\nH/SeYgE1ZkzdMmRIRntVkV1DMrP7w7d7gdVmtixcvxXY0gb7XgeMCkPjMHAHcGeDMi8R3Iy7NAyw\nk+5eaWYnmqprZkPc/UhY/8+AZKMFRUSyQ79+UFwcLImqqoJQioXV+vXBk3537QqGsCcG1Jgxwc3A\no0e3y2dTJRtlt5i602rW8L27P9rqnZvdRN3Q7efc/TEzu4dgB8+GZZ4BFgBngbvdfUNTdcPt/wuY\nErb3feAed69ssF/1kESk/Yv1qmLXp2Kh9e67MHBg/ZCKvS8ouOwplXRjbBookESkQ6upCZ5DlRhS\nsffHjwfPokoMqtjrNdck/drIA8nMBhEMbBgHxObVcHf/ZLoalW4KJBHptM6ehT176gJq9+669zk5\ndaf8EpeRI6F376wIpNeBXwAPAPcAdwHH3P1v09WodFMgiYg04A7HjtUFVOLy7rswYAB26FDkgbTB\n3aeZ2ZbYdEFmts7dZ6SrUemmQBIRaYGaGjh4ELv22sifh3QpfK0wsz8lGNXWL10NEhGRLJOTEzwI\nMc1SCaTvm1lf4H7gaeAq4G/S2ioREel0NMpORERSEuWNsQ+6+xNm9nQjH7u7/3W6GiUiIp1PslN2\nO8LXdQ22J94kKyIi0iaSPX7i5fC5Q5Pc/f6myomIiLSFLsk+dPcaYJ7pQSAiIpJmqYyy2wT8xsx+\nRfAYcQiuIf17+polIiKdTSqB1AP4AGg4VZACSURE2oyGfYuISEoiG/ad0ICeBE9cjU2u6gDu/l/S\n1SgREel8kg5qCP0LkEvwTKJSoAA4k8Y2iYhIJ5RKII10928DZ9z9BWAhMLstdm5mC8xsp5ntMbMH\nmyjzVPj5ZjObmmpdM7vfzGrNrH9btFVERNIrlUCKTa76kZlNBPoCA1u74/Aep9jTYMcBd5pZUYMy\nCwkCcRSwCPhJKnXNrAD4Y2B/a9spIiKZkUog/XPYy/jvwEsEMzj8oA32PQvY6+773L0KWArc2qDM\nLcALAO6+GuhrZoNTqPtDgocKiohIO5HKsO+fuXs18Afg2jbcdz5wMGH9EB8/FdhYmXwgr6m6ZnYr\ncMjdt+h+XhGR9iOVQHrPzH5P8NTY/9uG46VT/Z6UUyUcEfh3BKfrWlxfRESik0ogFQF/CtwHPG9m\nLwO/cPflrdx3OcGIvZgCgp5OsjJDwzLdmqh7HTAc2Bz2joYC681slrsfTfzixYsXx9+XlJRQUlJy\n2T9ERKQjKi0tpbS0NGP7a9GNsWbWD3gK+HN3z2nVjs26AruAGwmeQrsGuNPdyxLKLATuc/eFZlYM\n/Njdi1OpG9Z/H5ju7h802K4bY0VEWijyG2PDRpQAdxCMalsL3N7aHbt7tZndB7wK5ADPuXuZmd0T\nfv6su79iZgvNbC9wFrg7Wd3GdtPadoqISGY020Mys30EE6z+AnjZ3dv9TbHqIYmItFy6e0ipBNLV\n7v5RuhoQBQWSiEjLRR5IHZECSUSk5dIdSKncGCsiIpJ2zQaSmY1IZZuIiEhrpNJD+t+NbPtVWzdE\nREQ6tyaHfYeTlY4DrjazzxHMeODAVQRPkRUREWkzye5DGg3cDFwdvsacBv4ynY0SEZHOJ5Vh33Pd\n/Z0MtScjNMpORKTlIh/2bWaDCHpEw6nrUXl7foS5AklEpOWyYeqg3wBvAa8DteE2/TUXEZE2lUoP\naZO7T8lQezJCPSQRkZbLhhtj/4+Z/Um6GiAiIgKp9ZDOAL2AS0BVuNnd/ao0ty1t1EMSEWm5yK8h\nuXufdO1cREQkJpWpg7qY2ZfM7DvheqGZzUp/00REpDNJ5RrS/wDmAH8erp8Jt4mIiLSZVAJptrvf\nC5wHCB8H3q0tdm5mC8xsp5ntMbMHmyjzVPj5ZjOb2lxdM/teWHaTmb1hZgVt0VYREUmvVALpkpnl\nxFbMbCB19yNdtvA7nyF4LPo44M5w/rzEMguBke4+ClgE/CSFuj9w98nhUPVlwCOtbauIiKRfKoH0\nNPAfwCAz+3tgBfBYG+x7FrDX3fe5exWwFLi1QZlbgBcA3H010NfMBier6+6nE+r3AY63QVtFRCTN\nUhll969mth64Mdx0q7uXtcG+84GDCeuHgNkplMkH8pLVNbPvA18CzgHFbdBWERFJsyZ7SGbWP7YA\nlcDPw6Uy3NZaqd4I1OIx7+7+LXcvBP4/4EctrS8iIpmXrIe0gSA0DCgEPgy39wP2A9e2ct/lQOKA\ngwKCnk6yMkPDMt1SqAvwb8Arje188eLF8fclJSWUlJSk1moRkU6itLSU0tLSjO0vlZka/hn4D3d/\nJVy/Cfgzd1/Uqh2bdQV2EZwKPAysAe5MPB0YDmq4z90Xmlkx8GN3L05W18xGufuesP7XgVnu/qUG\n+9ZMDSIiLRT5TA3AHHePP5DP3X9nZv/Q2h27e7WZ3Qe8CuQAz4WBck/4+bPu/oqZLTSzvcBZ4O5k\ndcOvfszMxgA1wLvA11rbVhERSb9UekivETx+4l8JTt/9OXC9u38m/c1LD/WQRERaLhtm+74TGEQw\n9Pvfw/d3pqtBIiLSOTXbQ+qI1EMSEWm5yK8hhddjHuDjjzD/ZLoaJSIi2eHY2WPsOLaDsuNtcftp\ncqkMavgVwZQ9SwgGCoAeYS4i0mG4O4dOHYoHT9mxMnYc30HZsTJqvIZxA8dRNKCo+S9qpVQGNax3\n9+lpb0kG6ZSdiHRGNbU17P9oPzuO7ai3lB0vo0/3PvHgKRpQFLwfWERu71zMgrN06T5ll0ogLQaO\nEQxouBjbHs763S4pkESkI6uurea9D9+LB872Y9vZcWwHu47vYkCvAYwbOK7eUjSgiH49+zX7vdkQ\nSPto5BSdu7d2pobIKJBEpCOorq1m7wd72X40CJwdx4MA2n1iN0P6DKkXOuMHjmfsgLFcecWVl72/\nyAOpI1IgiUh70jB4th/bzvZj29n7wV7yrsxj/MDxjB84Ph4+YweMpXf33m3ejsgDycx6A98ACt39\nL81sFDDG3f9PuhqVbgokEclGNbU1vPfhe2w7ui0eOtuPbmfPB3viwRPr7YwfFPR4enXrlbH2ZUMg\n/RJYD3zZ3ceHAfWOu09OV6PSTYEkIlGq9Vr2n9wfD57Y667juxjUexATBk2Ih874geMpGliU0eBp\nSjYE0np3n25mG919arhtswJJRCQ5d+fw6cP1gmfb0W2UHS/j6iuujgfPhEETGD9oPEUDilp1jSfd\nIr8xFrhoZj0TGnQdCaPtREQEPjj/AduObmNr5dYgeI4F4dOtSzcmDJrAhEETmJ0/m69O/SrjB42n\nb4++UTc566TSQ/o08C1gHPA6MA+4y93fTH/z0kM9JBG5XOeqzrHj2I548Gw9GryeuXQmHjyJy6De\ng6JucpuJ/JRd2IgBBI8IN2CVux9PV4MyQYEkIs2prq1mz4k98dCJBU/5qXJGXzOaibkTmTBwQvA6\naAIFVxXEbyDtqCIPJAuO8OeATxDcj7Tc3f8jXQ3KBAWSiMTErvNsPbqVrZVb4+Gz6/gu8q7MY2Lu\nRCYOCkJn4qCJjLpmFF27pHK1o+PJhkD6CXAd8HOCHtLtwHvufm+6GpVuCiSRzunMpTNsO7qNLZVb\n2FK5JR5CXbt0jQfPxEETmZg7kfEDx6flXp72LBsCaScwzt1rw/UuwA53H9vqnZstAH5M8NTXJe7+\nRCNlngJuAs4RXLvamKxu+DTbPwUuETwx9m53/6jBdyqQRDqwmtoa9n6wl61Ht9YLn4ozFRQNKGJi\n7kQmDZoUD6HcPrlRN7ldyIZRdnuBQmBfuF4YbmsVM8sBngE+BZQDa83spYRHkWNmC4GR7j7KzGYT\nzDpe3Ezd14AH3b3WzB4HHgYeam17RSQ7nTh3gq1Ht7K5YnMQPke3sOPYDnJ75zIpdxKTcifxxYlf\nZFLuJEb2H0lOl5yomyxNSCWQrgLKzGwNwTWkWQQB8DLBc5Fuucx9zwL2uvs+ADNbCtwKJD504xbg\nBYIdrTazvmY2GLi2qbru/npC/dXA5y+zfSKSRapqqth9YjebKzfHez1bKrdw+tJpJg6ayKTcSczM\nn8lfTPsLxg8az1VXXBV1k6WFUgmk7yT5rDXnvfKBgwnrhwhG8jVXJh/IS6EuwH8huPYlIu3I8XPH\n4z2eWADtPL6TgqsLmJQ7icm5k/mrGX/FpNxJDLt6WIcf3dZZNBtI7l5qZsMJTp39p5n1Arq6+6lW\n7jvVMLusf2lm9i3gkrv/W2OfL168OP6+pKSEkpKSy9mNiLRCbGj15srNbKrYFA+fs5fOxk+3zSuY\nx9dmfI0JgyZokEGGlZaWUlpamrH9pTKoYRHwl0B/d7/OzEYDP3H3G1u1Y7NiYLG7LwjXHwZqEwc2\nmNlPgVIk8AdSAAAUXElEQVR3Xxqu7wRuIDhl12RdM7srbPON7n6hkX1rUINIhn104aN4jycWPjuO\n7WBInyFMHjyZybnhMniyej1ZKhsGNfxXgus9qwDcfbeZtcWtx+uAUWHv6zBwB3BngzIvAfcBS8MA\nO+nulWZ2oqm64ei7bwI3NBZGIpJe7s6Bjw6wqWJTPHg2VWyi8mwlEwZNYEruFKYNmcbdU+5mUu6k\nrJ67TTIrpbns3P1iwiNsu9K6a0cAuHu1md0HvEowdPs5dy8zs3vCz59191fMbKGZ7QXOAncnqxt+\n9dNAd+D1sM0r2/M9UyLZ7FLNJXYc2xEPn1gA9erWi8m5k5kyeAp3jL+Dx258TCPcpFmpnLL7B+Ak\n8GWC3sq9BPchfSv9zUsPnbITabmTF06yuSLo7Wys2Mimik3sPrGbEf1GMGXwFCbnTmbqkKlMzp3M\nwN4Do26upEE23BibA3wV+HS46VWCG1Hb7V90BZJI09yd8tPlbDyyMR48Gys2cuzsMSblTmLq4KlM\nGTyFKYOnMGHQBHp269n8l0qHEHkghY0YBODuR9PVkExSIIkEampr2PPBHjYc2VCv52MYU4dMZerg\nqfEA0ik3iSyQwklVHyE4TRf7V1hDcI3mu+35L7oCSTqji9UX2X5se7zns+HIBrZUbmFQ70H1wmfq\nkKkM6TNEo9zkY6IMpG8QzCG3yN3fD7eNAH4K/N7df5iuRqWbAkk6urOXzrK5cjMbjmxgw5ENbKzY\nyK7ju7iu/3X1gmfK4Cl6UJykLMpA2gT8sbsfa7B9IPC6u09JV6PSTYEkHcnJCyfZeCTo8WyoCALo\nwEcHGDdwHNMGT2PakGlMHTKViYMm6nqPtEqUgbTN3Se09LP2QIEk7dXxc8fZcGQD6w+vj4fP0bNH\nmTJ4ClMHT2XakCCAigYU0S2nW9TNlQ4mykDa6O5TW/pZe6BAkvag4kxFEDxhz2f94fWcungqHjqx\nZVT/URpsIBkRZSDVEDyDqDE93b3dPjJRgSTZ5vDpw6w/vJ71R4Jlw5ENXKi+wLQh05g+ZDrTh0xn\n2pBpXNvvWrpYl6ibK51UVgz77mgUSBKlw6cPs+7wOtYfXs+6I+vYcGQDVTVVTM+bHg+f6XnTNZ+b\nZB0FUhookCRTjpw+wvoj64MACl+ra6uZPmQ6M/JmxMOn4KoChY9kPQVSGiiQJB2Onj0a9HoOr2Pd\nkXWsO7yOC9UX4sETey28ulDhI+2SAikNFEjSWh+c/+Bj4fPRhY+YkTcjvkwfMp3hfYcrfKTDUCCl\ngQJJWuL0xdNsOLKBtYfXsu7wOtYeXsuxs8eYOmQqM/NmxgNoRL8RGnAgHZoCKQ0USNKUC9UX2FSx\nibXla+MBtP+j/UzKnRQPn5l5MxkzYIzCRzodBVIaKJAEgsdnbz+6nbWH18YDaOfxnYwdMDYePDPz\nZzJ+4HjdZCqCAiktFEidj7vz3ofvsaZ8DWvK17D28Fo2VWyi4OqCIHjC8JmcO1nT64g0oUMHUvi4\n8R8TzCa+xN2faKTMUwSTvJ4D7nL3jcnqmtltwGJgLDDT3Tc08p0KpA6u8kxlPHhir7279WZW/ixm\n5c9iZt5MpudN56orroq6qSLtRocNpPDBf7uATwHlwFrgzoRHkWNmC4H73H2hmc0G/sndi5PVNbOx\nQC3wLHC/AqnjO3PpDOsPrw96P4eDHtDpi6eZmT+TWXlhAOXPZHCfwVE3VaRdS3cgRTn9zyxgr7vv\nAzCzpcCtQFlCmVuAFwDcfbWZ9TWzwcC1TdV1953htgz9DMmk6tpqth3dFj/1tqZ8De9++C6Tcicx\nK28Wnx3zWf7+k3/PyP4j9W9ApJ2JMpDygYMJ64eA2SmUyQfyUqgr7Zy7c+CjA6wpX8Pq8tWsLl/N\nxiMbKby6MH7q7a9m/BWTcifRPad71M0VkVaKMpBSPWeWlv+bu3jx4vj7kpISSkpK0rEbaYFTF0+x\ntnxtPHxWH1qN48zOn83s/Nk8csMjzMybydU9ro66qSKdQmlpKaWlpRnbX5TXkIqBxe6+IFx/GKhN\nHNhgZj8FSt19abi+E7iB4JRdc3XfRNeQslZNbQ3bj21n9aHVrDq0itXlq3n/5PtMGTyF4vxiZuXP\nYvbQ2ZpgVCSLdORrSOuAUWY2HDgM3AHc2aDMS8B9wNIwwE66e6WZnUihLqSpdyUtV3GmglWHVsXD\nZ93hdeRdmcfs/NkUDy3m3pn3Mil3ku73EenEoh72fRN1Q7efc/fHzOweAHd/NizzDLAAOAvcHevx\nNFY33P5nwFPAAOAjYKO739Rgv+ohpdHF6otsqtjEqkOrWHloJasOreLUxVPMyp9F8dBiiocGPaD+\nPftH3VQRaYEOO+w7SgqktuPuHDx1kJUHg+BZVb6KLZVbGH3NaIrzi+MBNOqaUZpqR6SdUyClgQLp\n8p2vOs/6I+vjvZ+VB1dS4zXMGTqH4qHFzBk6h+l50+nTvU/UTRWRNqZASgMFUmpiw67fOfhO/NTb\n9mPbGTdwHMX5xcwpmMOcoXP0iAWRTkKBlAYKpMZdqL7AhiMbWHlwJe8ceoeVB1dS67Xx4In1fnp1\n6xV1U0UkAgqkNFAgBQ6fPsw7B9+J94C2VG5h7ICx8fCZWzBXvR8RiVMgpUFnDKTq2mo2V2xm5aGV\n8RA6fek0cwvmMnfoXOYUzGFm3kx6d+8ddVNFJEspkNKgMwTSh+c/jIfPioMrWHd4HYVXFzJ36Nwg\nhArmMvqa0er9iEjKFEhp0NECyd3Z+8FeVhxcEQ+gAx8dYGbeTOYVzGNuwVyKhxbTr2e/qJsqIu2Y\nAikN2nsgXay+yPoj61lxYEU8hHp07cHcgrnxAJo8eDJdu0Q5EYeIdDQKpDRob4F04tyJeM9nxcEV\nbDyykdHXjGZewTzmFc5jXsE8Cq4uiLqZItLBKZDSIJsDKfao7bcPvM2Kgyt4+8DbHDp1iOKhxfEA\nmp0/myuvuDLqpopIJ6NASoNsCqTY6Le3D7zN2wff5u0Db9PFujC/cD6fKPwE8wrmMTF3ok6/iUjk\nFEhpEGUgnb10llWHVsUDaPWh1RReXcgnCj8RX/TIBRHJRgqkNMhkIB09e5QVB4JTb8sPLGfHsR1M\nGTwlHj5zC+Zq1msRaRcUSGlgZl5T43Rp48mn3Z19J/ex/MBylu9fzvIDy6k4U8Hcgrl8ovATzC+c\nz4y8GfTs1rNtdywikgEKpDQwM4cgkLp1C5bu3eveN7Wt4XrXbrWc672dE32Wc6zXciqvWI5bDUNr\n5lPIfK7tMp/8rhO5onvOx767LZeuXUFn+EQk3RRIaWBmXlvr1NRAVVXdculS/fWG285frKLs5AY2\nn1zO1lNvsfPcCnp36c+o7vMZ0TUIoKtrrqO62j72Paksje2/qgqqq5PXq6mBnJzmQ6u5bS0t05L3\nLd0We23rXqyIXL4OHUhmtoC6p74ucfcnGinzFHATcA64y903JqtrZv2BXwDDgH3A7e5+ssF3pnQN\n6XzVeVaXr+at/W+x/MByVh1axYh+I7i+8HrmD5vP/ML5DLlyyGX//rbi3nxoNfy8sfLNlUlcv9z3\nDV+b+8ys5SGW7LW5bamUSaVeU9uaKp+To16uZL8OG0hmlgPsAj4FlANrgTvdvSyhzELgPndfaGaz\ngX9y9+Jkdc3sB8Bxd/+BmT0I9HP3hxrsu9FAOn3xNCsOruCt/W/x1v632FixkYmDJnL9sOu5ftj1\nzCuYp+l3MqympvFgbCrEkr02ty1xvbFtjQVmrH1NfU9T64nfU1NTv5ebLNCSLU2VS7V+Wy45OS37\nXD3h9iHdgRTlzS2zgL3uvg/AzJYCtwJlCWVuAV4AcPfVZtbXzAYD1yapewtwQ1j/BaAUqBdIMR+e\n/5DlB5bz1v63+MP+P1B2rIwZeTO4ftj1PFryKMVDizX7dcRycoLliiuibkl6xXq5DZeGodaSMg0D\nMDE8Gy7nzyf/PJXvaPh5Y2UbKxPrCbc03BK3NyyTbD0dZRurk+prU9s6Y0hHGUj5wMGE9UPA7BTK\n5AN5Sermuntl+L4SyG1s55N/Opn3P3yf4qHF3DDsBn746R8yM38mPbr2uLxfI9IKsVOT3bpF3ZJo\n1NY2H3yJveXY+2QB2JL12PtLl+oHZaz3mqwdzbWl4ecN6zT2PVVVwXG5nIBL5bWp98m2dc1AWkQZ\nSKmeK0yle2iNfZ+7ezCi7uNm7p7JLX1uIefdHOYVzGP+sPkpNkdE2lqXLsEo1u7do25J9qitbTzM\nGnttLOia+7yp94nbdu0qZc+eUmprg158ukUZSOVA4oygBQQ9nWRlhoZlujWyvTx8X2lmg929wsyG\nAEcb2/mSHy5pRdNFRNKrSxfit6ZEpyRcAmaPpnVvUZ6lXAeMMrPhZtYduAN4qUGZl4AvA5hZMXAy\nPB2XrO5LwFfC918BlqX3Z4iISFuIrIfk7tVmdh/wKsHQ7efCUXL3hJ8/6+6vmNlCM9sLnAXuTlY3\n/OrHgV+a2VcJh31n9IeJiMhl6bQ3xnbG3y0i0hrpHvbdCQcWiohINlIgiYhIVlAgiYhIVlAgiYhI\nVlAgiYhIVlAgiYhIVlAgiYhIVlAgiYhIVlAgiYhIVlAgiYhIVlAgiYhIVlAgiYhIVlAgiYhIVlAg\niYhIVlAgiYhIVogkkMysv5m9bma7zew1M+vbRLkFZrbTzPaY2YPN1Q+3v2lmp83s6Uz9HhERab2o\nekgPAa+7+2jgjXC9HjPLAZ4BFgDjgDvNrKiZ+heA/w48kN7mdxylpaVRNyFr6FjU0bGoo2OROVEF\n0i3AC+H7F4DPNlJmFrDX3fe5exWwFLg1WX13P+fuK4CL6Wp4R6P/2OroWNTRsaijY5E5UQVSrrtX\nhu8rgdxGyuQDBxPWD4XbUqmv55OLiLQzXdP1xWb2OjC4kY++lbji7m5mjQVIw23WyLZk9UVEpB0x\n98z/LTeznUCJu1eY2RDgTXcf26BMMbDY3ReE6w8Dte7+RHP1zewrwAx3/3oT+1eAiYhcBne3dH13\n2npIzXgJ+ArwRPi6rJEy64BRZjYcOAzcAdyZYv2kByydB1RERC5PVD2k/sAvgUJgH3C7u580szzg\nn939T8JyNwE/BnKA59z9sWT1w8/2AVcC3YEPgU+7+86M/TgREbkskQSSiIhIQ+1ypoambphtUOap\n8PPNZja1ubrJbtY1s4fD8jvN7NPp/XUtk8ljYWZ/bGbrzGxL+PpH6f+Fqcv0v4vw80IzO2Nm96fv\nl7VcBP+NTDKzlWa2Lfz3cUV6f2HqMvzfSA8z+3l4DHaY2cfusYxSmo7FbWa23cxqzGxag+9q2d9O\nd29XC8Hpu73AcKAbsAkoalBmIfBK+H42sKq5usAPgL8N3z8IPB6+HxeW6xbW2wt0ifo4RHQspgCD\nw/fjgUNRH4OojkXCd/4a+AVwf9THIMJ/F12BzcDEcL1fJ/5v5C7g5+H7nsD7QGHUxyHNx2IsMBp4\nE5iW8F0t/tvZHntIyW6YjYnfOOvuq4G+Zja4mbpN3ax7K8E/sCp330dwUGel5Ze1XEaPhbtvcveK\ncPsOoKeZdUvPT2uxTP+7wMw+C7xHcCyySaaPxaeBLe6+Nfy+D929Nj0/rcUyfSyOAL0tmGmmN3AJ\nOJWWX9ZyaTkW7r7T3Xc3sr8W/+1sj4GU7IbZ5srkJanb1M22eWG5ZPuLSqaPRaLPA+vDf5zZIKPH\nwsz6AH8LLG6Dtre1TP+7GA24mf3ezNab2Tdb/xPaTEaPhbu/ShBARwgGXP2DhwOuskC6jkVTWvy3\nM6ph362R6iiMVIZ2X+7NttkyEiSSY2Fm44HHgT9Ocf+ZkOljsRj4kbufM7Nsu40g08eiK/AJYAZw\nHnjDzNa7+/9NsR3plNFjYWb/D8GpuiFAf2C5mb3h7u+n2I50astjkZY2tMdAKgcKEtYLqJ/CjZUZ\nGpbp1sj28vB9pZkN9rqbbY8m+a5yskOmjwVmNhT4d+BLWfIfWUymj8Us4PNm9gOgL1BrZufd/X+0\nya9pnUwfi4PAW+7+AYCZvQJMA7IhkDJ9LOYC/+HuNcAxM1tBENTZ8N9KWx6Lxuo2t7/m/3ZGfaGt\npQtBiL5LcJGsO81fmCum7sJck3UJLlI+GL5/iI8PaugOXBvWt6iPQ0THoi/BxevPRv3boz4WDb73\nEeAbUR+DCP9d9APWE/QMugKvAzdFfRwiOhZ/DTwfvu8NbAcmRH0c0nksEuq+CUxPWG/x387ID9Jl\nHtibgF0EF8keDrfdA9yTUOaZ8PPN1B/58bG64fb+wH8Cu4HXgL4Jn/1dWH4n8Jmof39Ux4Lg0R5n\ngI0Jy4Coj0FU/y4SymRVIEVxLIAvAtuArTQS2p3lWABXAP8aHoftZNHoyzQeiz8j6CWfByqA3yV8\n1qK/nboxVkREskJ7HGUnIiIdkAJJRESyggJJRESyggJJRESyggJJRESyggJJRESyggJJJAuY2S/M\n7LoWlJ9kZs+ls00imaZAEomYmY0Eerv7u6nWcfctwHVmNih9LRPJLAWSSBqZ2T1mtjFc3jezxuZ3\n+wLwUkKdM2b2g/Bhd6+bWbGZ/cHM3jWzmxPq/Q64Ld2/QSRTFEgiaeTuz7r7VGAmwfQq/9hIsXnA\nuoT1XsAb7j4BOA18F/gkwRQt300otwa4Ph3tFolCe5ztW6Q9eoogZH7byGfDCJ6fE3PJg+fqQDAn\n2gV3rzGzbQSTW8YcabAu0q4pkETSzMzuAgrc/d5kxRLeJz70sJbgqaO4e62ZdW1QR5NRSoehQBJJ\nIzObDtwPzE9SbD/BA90Ot/Drh4R1RToEXUMSSa//SvC8oDfDgQ3/s5EybxM8xC2mYa/Hm3g/C3ir\nTVopkgX0+AmRiJnZCOBpd/+TFtYrBW5396PNlRVpD9RDEomYu78HnG7pjbHAXoWRdCTqIYmISFZQ\nD0lERLKCAklERLKCAklERLKCAklERLKCAklERLKCAklERLLC/w++rQb2rroGzQAAAABJRU5ErkJg\ngg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe928e70ad0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(tspan,ty)\n",
    "plt.legend(['$C_A$','$C_B$','$C_C$'])\n",
    "plt.xlabel('z (m)')\n",
    "plt.ylabel('Dependent variable')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using our ODEs, we have described the gas concentrations over distance for the diffusion problem.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear and Polynomial Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Regressing model parameters for a linear or polynomial equation is very common and straightforward.  The goal in each case is determining the parameter values $b_0, b_1...b_n$ given a set of observations $y_i$ at different inputs $x_i$. \n",
    "\n",
    "The form of a linear equation is:\n",
    "\n",
    "$$\\hat{y_i}=b_0+b_1x_i$$\n",
    "\n",
    "For a polynomial:\n",
    "\n",
    "$$\\hat{y_i}=b_0+b_1x_i+b_2x_i^2+...+b_nx_i^n$$\n",
    "\n",
    "Where $\\hat{y_i}$ is the calculated observation.  The goal of the regression is to find values of $b$ so that $\\hat{y_i}$ matches $y_i$ for each of the values of $x_i$.  To do regression, we need data, which is typically a table.  For working with tables of data, the python `pandas` library is indespensible"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Temperature</th>\n",
       "      <th>Pressure</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>290.08</td>\n",
       "      <td>8634</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>302.39</td>\n",
       "      <td>15388</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>311.19</td>\n",
       "      <td>22484</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>318.69</td>\n",
       "      <td>30464</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>325.10</td>\n",
       "      <td>38953</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>330.54</td>\n",
       "      <td>47571</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>334.89</td>\n",
       "      <td>55511</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>338.94</td>\n",
       "      <td>63815</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>342.95</td>\n",
       "      <td>72985</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>346.24</td>\n",
       "      <td>81275</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>349.91</td>\n",
       "      <td>91346</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>353.47</td>\n",
       "      <td>102040</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>356.19</td>\n",
       "      <td>110850</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>358.87</td>\n",
       "      <td>120140</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>362.29</td>\n",
       "      <td>132780</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>365.23</td>\n",
       "      <td>144530</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>367.90</td>\n",
       "      <td>155800</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>370.53</td>\n",
       "      <td>167600</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>373.15</td>\n",
       "      <td>180060</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>375.84</td>\n",
       "      <td>193530</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>378.52</td>\n",
       "      <td>207940</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>381.32</td>\n",
       "      <td>223440</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    Temperature  Pressure\n",
       "0        290.08      8634\n",
       "1        302.39     15388\n",
       "2        311.19     22484\n",
       "3        318.69     30464\n",
       "4        325.10     38953\n",
       "5        330.54     47571\n",
       "6        334.89     55511\n",
       "7        338.94     63815\n",
       "8        342.95     72985\n",
       "9        346.24     81275\n",
       "10       349.91     91346\n",
       "11       353.47    102040\n",
       "12       356.19    110850\n",
       "13       358.87    120140\n",
       "14       362.29    132780\n",
       "15       365.23    144530\n",
       "16       367.90    155800\n",
       "17       370.53    167600\n",
       "18       373.15    180060\n",
       "19       375.84    193530\n",
       "20       378.52    207940\n",
       "21       381.32    223440"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# input data points\n",
    "T=[290.08, 302.39, 311.19, 318.69, 325.1, 330.54, 334.89, 338.94,\n",
    "   342.95, 346.24, 349.91, 353.47, 356.19, 358.87, 362.29, 365.23,\n",
    "   367.9, 370.53, 373.15, 375.84, 378.52, 381.32,]\n",
    "P=[8634, 15388, 22484, 30464, 38953, 47571, 55511, 63815, 72985,\n",
    "   81275, 91346.0, 102040, 110850, 120140, 132780, 144530, 155800,\n",
    "   167600, 180060, 193530, 207940, 223440]\n",
    "\n",
    "# convert to data table- first stick T to P, then transpose and add labels\n",
    "data=pd.DataFrame(np.transpose(np.vstack((T,P))),columns=['Temperature','Pressure'])\n",
    "data # print out data table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Temperature and pressure are changed to column vectors for sklearn\n",
    "X=data['Temperature'].values[:,np.newaxis] \n",
    "y=data['Pressure'].values[:,np.newaxis]\n",
    "model=make_pipeline(PolynomialFeatures(degree=3), LinearRegression())\n",
    "model.fit(X,y)\n",
    "yhat=model.predict(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make plot of predicted vs. measured pressures, get your magnifying glasses..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEPCAYAAACZcRnqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4lNX5//H3jSiixoQAskNihV+FurC6tTBoZam1alGh\n1agVxYItWO23ArYC7uBSlxal1eJWtCpa1CqLSNC2AoqgIrWACoWAESgglJ3cvz/mSZjESUhCZv+8\nrmuunDnPdg4T5s5ztsfcHRERkXiol+gCiIhI5lDQERGRuFHQERGRuFHQERGRuFHQERGRuFHQERGR\nuIlZ0DGzNmY2x8w+NrMlZjY8yB9rZmvMbFHw6h9xzCgzW25mn5hZn4j8rmb2UbDtgYj8Bmb2lyB/\nnpm1i9h2uZktC16XxaqeIiJSfRareTpm1hxo7u6LzewoYCFwPnAxsNXd76uwf0dgCtAdaAW8AbR3\ndzezBcDP3H2Bmb0GPOju081sGPAtdx9mZgOBC9x9kJnlAu8CXYPTLwS6uvvmmFRWRESqJWZ3Ou7+\nhbsvDtLbgH8RDiYAFuWQ84Bn3H2Pu68EVgCnmFkLIMvdFwT7PUk4eAH8AHgiSE8FzgrSfYGZ7r45\nCDSzgH51VjkREamVuPTpmFke0BmYF2T93Mw+MLPHzCwnyGsJrIk4bA3hIFUxv4j9wasVsBrA3fcC\nW8yscRXnEhGRBIp50Ama1l4ARgR3PA8D+cDJwDrg3liXQUREkkP9WJ7czA4l3Oz1tLv/FcDdv4zY\n/ijwSvC2CGgTcXhrwncoRUG6Yn7pMW2BtWZWH8h2941mVgSEIo5pA7wZpXxaeE5EpBbcPVo3yQHF\ncvSaAY8BS939/oj8FhG7XQB8FKRfBgaZ2WFmlg+0Bxa4+xfAV2Z2SnDOAmBaxDGXB+kLgdlBeibQ\nx8xyzKwRcDYwI1o53T1tX2PGjEl4GVQ/1U31S7/XwYjlnc4ZwKXAh2a2KMgbDfzIzE4GHPgcuAbA\n3Zea2XPAUmAvMMz3124Y8DjQEHjN3acH+Y8BT5nZcmAjMCg413/N7FbCI9gAxrlGromIJFzMgo67\n/53od1KvV3HMHcAdUfIXAidEyd9FeAh2tHNNBiZXt7wiIhJ7WpEgjYVCoUQXIabSuX7pXDdQ/TJZ\nzCaHpgIz80yuv4hIbZgZXsuBBDEdvZaqwuMVJN3pDw6R+FPQqYS+kNKb/rAQSQz16YiISNwo6IiI\nSNwo6IiISNwo6EiNvP3223zzm99MdDGqJS8vj9mzZx94RxGJGwWdFJKXl8cRRxxBVlZW2Wv48OFx\nLcN3vvMdPvnkk7heE2DlypXUq1ePkpKSah9jZhowIJJkFHRqyN0ZOXLCQY1uq+05zIxXX32VrVu3\nlr0efPDBWpejpvbu3Ru3a1VGowpFUpuCTg1NnTqDiRPX8eKLMxN6jkhDhw7lwgsvLHt/44038t3v\nfheAwsJCWrduzZ133knTpk3Jz89nypQpZfvu2rWLX/7yl7Rr147mzZszdOhQdu7cWe7YCRMm0KJF\nCwYPHkxhYSFt2uxfDDwvL4977rmHE088kaysLAYPHkxxcTH9+/cnOzubs88+m82b9y97N2/ePE4/\n/XQaNWrEySefzNy5c8u2hUIhbr75Zr797W9z9NFH07dvXzZu3AhAz549AcjJySErK4v58+fz6aef\ncuaZZ9KkSROaNm3KpZdeypYtW+rk31REYiTRq5UmeKVUjyZa/iOPPOUdO57j7duPdijx9u1He8eO\n5/gjjzwV9RzRHOw58vLy/I033vha/vbt271Dhw7++OOP+1tvveVNmjTxoqIid3efM2eO169f32+4\n4QbfvXu3z50714888kj/97//7e7u1113nZ933nm+adMm37p1q5977rk+atSocseOHDnSd+/e7Tt2\n7PA5c+Z469aty5XptNNO8y+//NKLior8mGOO8c6dO/vixYt9586dfuaZZ/q4cePc3X3NmjXeuHFj\nf/31193dfdasWd64cWPfsGGDu7v36tXLjzvuOF++fLnv2LHDQ6GQjxw50t3dV65c6Wbm+/btK7v2\nihUr/I033vDdu3f7+vXrvWfPnn7dddeVK9vs2bOj/ltW9tmLyIEF/39q971b2wPT4VWToFNSUuLP\nPfeat2kz0sG9TZuR/vzzr3tJSUnUc0RzsOdo166dH3XUUZ6Tk1P2evTRR93dff78+d6oUSNv166d\nP/vss2XHlAaO7du3l+VdfPHFfuutt3pJSYkfeeSR/umnn5Zt++c//+n5+fllxx522GG+a9eucuer\nGHSmTJlS9n7AgAE+bNiwsvcPPfSQn3/++e7uftddd3lBQUG5OvXt29efeOIJd3cPhUJ+++23l22b\nOHGi9+vXz93dP//8868FnYpeeukl79y5c7myKeiI1L2DCTpakaCaSjulN2/eSceO17N6dUmNO6oP\n9hxmxrRp0zjzzDO/tq1Hjx4ce+yxbNiwgYsuuqjctkaNGtGwYcOy9+3atWPdunVs2LCB7du307Vr\n17Jt7l6us75p06YcdthhVZarWbNmZemGDRuWe3/44Yezbds2AFatWsXzzz/PK6+8UrZ979695erT\nvHnzcucqPTaa4uJiRowYwd///ne2bt1KSUkJubm5VZZVRBJLfTo1sHz5aiZP7seSJfcyeXJ/li9f\nnZBzRPP73/+e3bt307JlSyZMmFBu26ZNm9i+fXvZ+1WrVtGyZUuaNGlCw4YNWbp0KZs2bWLTpk1s\n3ryZr776qmzf2oz+Cv8h9HVt27aloKCg7FqbNm1i69at/OpXvzrgOaOVY/To0RxyyCEsWbKELVu2\n8NRTT9VodJuIxJ/udGpg1Kiry9IDBvRNyDmifaEvW7aM3/zmN8ydO5eGDRvSo0cP+vfvz0knnVS2\nz5gxY7jjjjuYN28ef/vb37j11lsxM66++mquu+46fve739G0aVOKior4+OOP6dOnT63qV5VLL72U\n7t27M3PmTM466yz27NnDvHnzaN++Pa1ataq0fhC+46pXrx6ffvop7du3B2Dbtm1kZ2dz9NFHU1RU\nxN13313nZRaRuqU7nRRz7rnnlpun88Mf/pCCggJGjhzJCSecwHHHHccdd9xBQUEBe/bsAcJNVo0a\nNaJly5YUFBQwadIkOnToAMD48eM57rjjOPXUU8tGmy1btqzsetHuMA509xO5PbL5sHXr1kybNo07\n7riDY445hrZt23LvvfeWCzSVHXvEEUdw0003ccYZZ5Cbm8uCBQsYM2YM77//PtnZ2Zx77rkMGDBA\n83Ik6bg7Px59TaV/UGUaPU8nSv2DZ0UkoER1r7CwkIKCAlavrptmvHSRTp+xJLcXXpjOJX8cy5Qh\n42rdQpJsDuZ5OrrTERGJgUmTnqZTp+8zevTb7N7dl1Gj3qJTp+8zadLTiS5aQqlPJwOoyUkk/jr0\naUXHvdnMnDkXuvyD4vfPoG/fPDr0aZXooiWUgk6aC4VC/Oc//0l0MUQyTu/83mxstosZc2bQpME+\nds3pxsDL+tM7v3eii5ZQCjoiIjFSOkXiw9yjOXHg6XU2RSKVaSBBmg8kkOj0GUs8Fa4sJJQXSnQx\n6szBDCRQ0FHQyUj6jEVqT6PXREQkJSjoiIhI3CjoZJB69erx2WefHdQ5xo4dS0FBQR2VqGoVn90j\nIqlPQScFTZkyhW7dupGVlUXLli353ve+xz/+8Y+4XLsmc36uuOIKfvOb38SwNCKSahR0aqFwZWHC\nznHffffxi1/8gl//+td8+eWXrF69mmuvvZaXX375oMtUHep8F5GDoaBTC4kKOlu2bGHMmDFMnDiR\n888/n4YNG3LIIYdwzjnnMH78eBYsWMBpp51Wtrjnz3/+87JFPyvasWMHN9xwA3l5eeTk5PCd73yH\nnTt3Rm3SysvL480334x6nosuuogWLVqQk5NDr169WLp0KQB/+MMfmDJlChMmTCArK4vzzjsPgLVr\n1zJgwACOOeYYjj32WB566KFyZbriiivIzc2lU6dOvPvuuzX+NxKR5Kagk0Leeecddu7cyQUXXBB1\ne/369XnggQfYuHEj77zzDrNnz2bixIlR9/3lL3/JokWLeOedd/jvf//L3XffTb160X8dqmpSO+ec\nc1ixYgXr16+nS5cuXHLJJQAMGTKESy65hBtvvJGtW7cybdo0SkpKOPfcc+ncuTNr165l9uzZ3H//\n/cycOROAcePG8fnnn/PZZ58xY8YMnnjiCS3hI0lPq0jXjFYkqKbClYVldyfj5o4ryw/lhao96etg\nz7Fx40aaNGlSaXDo0qVLWbpdu3YMGTKEuXPnMmLEiHL7lZSUMHnyZObPn0+LFi0AOPXUU6tVh4qu\nuOKKsvSYMWN44IEH2Lp1K1lZWUD55rh3332XDRs28Otf/xqA/Px8rrrqKp599ln69OnD888/z8MP\nP0xOTg45OTmMGDGCW265pVblEomXqVNnMHXhBwx4cWbarCIdSwo61VQxMIwNjY37ORo3bsyGDRso\nKSmJGniWLVvG9ddfz8KFC9m+fTt79+6lW7duX9tvw4YN7Ny5k2984xs1rUI5+/bt46abbuKFF15g\n/fr1ZWXasGFDWdCJtGrVKtauXUujRo3KnaNnz55AuOktsmmvbdu2B1U+kViaNOlpHnzwWfbsOYnd\nrcKrSN9880MMHz6Ia665NNHFS1pqXkshp512Gg0aNOCll16Kun3o0KF07NiRFStWsGXLFm6//fao\nj29u0qQJhx9+OCtWrPjatiOPPLLco6337dvH+vXro15vypQpvPzyy8yePZstW7bw+eefA/vvbio2\njbVt25b8/Pxyj6v+6quvePXVVwFo0aJFucVJtVCpJLMOfVrRcVg2xcfPhdAtFB8/l07X5mT8KtIH\noqBTC3WxhlJtzpGdnc0tt9zCtddey7Rp09i+fTt79uzh9ddf58Ybb2Tbtm1kZWVxxBFH8Mknn/Dw\nww9HPU+9evW48soruf7661m3bh379u3jnXfeYffu3XTo0IGdO3fy2muvsWfPHm677TZ27doV9Tzb\ntm2jQYMG5Obm8r///Y/Ro0eX296sWbNy84J69OhBVlYWEyZMYMeOHezbt48lS5bw3nvvAXDxxRdz\n5513snnzZtasWVNukIFIsumd35uBzQrwOd1p8vGp+JxuDGxWkPGrSB+Igk4tJCroAFx//fXcd999\n3HbbbWWPfJ44cSIXXHAB99xzD1OmTOHoo49myJAhDBo06GuPfy51zz33cMIJJ9C9e3caN27MqFGj\ncHeys7OZOHEiV111Fa1bt+aoo44q1+QV+Qjpyy67jHbt2tGqVSu+9a1vcdppp5W7xuDBg1m6dCmN\nGjXihz/8IfXq1ePVV19l8eLFHHvssTRt2pQhQ4bw1VdfAeE+oXbt2pGfn0+/fv247LLLNJBAklrp\nKtLDhvZh8uT+WkW6GmK24KeZtQGeBI4BHPiDuz9oZrnAX4B2wErgYnffHBwzCrgS2AcMd/eZQX5X\n4HHgcOA1dx8R5DcIrtEF2AgMdPdVwbbLgZuC4tzm7k9GKaMW/MxQ+oylLqXbKtIHkpSrTJtZc6C5\nuy82s6OAhcD5wE+ADe4+wcxuBBq5+0gz6whMAboDrYA3gPbu7ma2APiZuy8ws9eAB919upkNA77l\n7sPMbCBwgbsPCgLbu0DXoDgLga6lwS2ijAo6GUqfsUjtJeUq0+7+hbsvDtLbgH8RDiY/AJ4IdnuC\ncCACOA94xt33uPtKYAVwipm1ALLcfUGw35MRx0SeaypwVpDuC8x0981BoJkF9Kv7WoqISE3EpU/H\nzPKAzsB8oJm7FwebioFmQbolsCbisDWEg1TF/KIgn+DnagB33wtsMbPGVZxLREQSKObzdIKmtanA\nCHffGtkxHDSdJbSNY+zYsWXpUChEKBRKWFlERJJRYWEhhYWFdXKumAYdMzuUcMB5yt3/GmQXm1lz\nd/8iaDr7MsgvAiIX/WpN+A6lKEhXzC89pi2w1szqA9nuvtHMioBQxDFtgKiLh0UGHRER+bqKf5CP\nGzeu8p0PIGbNaxa+pXkMWOru90dsehm4PEhfDvw1In+QmR1mZvlAe2CBu38BfGVmpwTnLACmRTnX\nhcDsID0T6GNmOWbWCDgbmFHnlRQRkRqJ5ei1bwNvAR8SHjINMApYADxH+A5lJeWHTI8mPGR6L+Hm\nuBlBfumQ6YaEh0wPD/IbAE8R7i/aCAwKBiFgZj8BSmcr3ubupQMOIstY6eg1SX8avSbV4e5cctNP\n+fPtj+i7IZCUQ6ZTQWVBR0Sk1AsvTOeSP45lypBxWtAzoKBTSwo6IlKZyAU9l7eqT/uivRx66Ada\n0JODCzpaZVpEJIoOfVrRcW82M2fOhS7/oPj9M+jbN08Leh4kBR0RkSh65/dmY7NdzJgzgyYN9rFr\nTjcGXtZfC3oeJAUdEZFKlC7o+WHu0Zw48HQt6FkH1KeTwfUXkerJtAU9D0QDCWpJQUdEpOaScsFP\nERGRihR0REQkbhR0REQkbhR0REQkbhR0RCTjuTs/Hn2N1uOLAwUdEcl4U6fOYOrCD3jxxZmJLkra\n05DpDK6/SKbT+mq1o7XXRERqQeurxZ+CjohkLK2vFn8KOiKS0bS+WnypTyeD6y8i+2l9terT2mu1\npKAjIlJzWntNRERSgoKOiIjEjYKOiIjEjYKOiGQELXWTHBR0RCQjaKmb5KDRaxlcf5FMoKVu6p6W\nwRERqYSWukkuCjoikta01E1yUdARkbSnpW6Sh/p0Mrj+IplGS93UDS2DU0sKOiIiNadlcEREJCUo\n6IiISNwo6IiISNwo6IhIWtAyN6lBQUdE0oKWuUkNGr2WwfUXSQda5ib+knb0mpn9ycyKzeyjiLyx\nZrbGzBYFr/4R20aZ2XIz+8TM+kTkdzWzj4JtD0TkNzCzvwT588ysXcS2y81sWfC6LJb1FJHE6dCn\nFR2HZVN8/FwI3ULx8XPpdG2OlrlJUrFuXpsM9KuQ58B97t45eL0OYGYdgYFAx+CYiWZWGkkfBga7\ne3ugvZmVnnMwsDHI/y0wPjhXLnAz0CN4jTGznFhVUkQSp3d+bwY2K8DndKfJx6fic7oxsFmBlrlJ\nUtVaBsfMjgfygBJglbt/Up3j3P1tM8uLdsooeecBz7j7HmClma0ATjGzVUCWuy8I9nsSOB+YDvwA\nGBPkTwV+F6T7AjPdfXNQ/lmEA9mz1Sm3iKQWLXOTOioNOmaWD/wC+B5QBKwlHCxamFlr4FXgt+6+\nshbX/XnQ5PUecEMQHFoC8yL2WQO0AvYE6VJFQT7Bz9UA7r7XzLaYWePgXGuinEtE0tCoUVcD0Hhl\nA0K9Q4ktjFSpqua18cArwPHu3svdf+Tug9y9F/BN4G/AhFpc82EgHzgZWAfcW4tziIh8jdZVS36V\n3um4+8VVbNsDzAxeNeLuX5amzexRwoENwncwbSJ2bU34DqUoSFfMLz2mLbDWzOoD2e6+0cyKgFDE\nMW2AN6OVZ+zYsWXpUChEKBSKtpuISMYqLCyksLCwTs5VrSHTZnYC4Q7+wwkPBMDdn6zWBcJ9Oq+4\n+wnB+xbuvi5I/wLo7u4/DgYSTCHc8d8KeAM4zt3dzOYDw4EFhO+wHnT36WY2DDjB3Yea2SDgfHcf\nFAwkeA/oQrhJcCHQpbSPJ6JsGjItIlJDMX1yqJmNBXoBnQh/4fcH/k64Q/9Axz4THNvEzFYT7vQP\nmdnJhIPX58A1AO6+1MyeA5YCe4FhERFhGPA40BB4zd2nB/mPAU+Z2XJgIzAoONd/zexW4N1gv3EV\nA46IpB5355Kbfsqfb3+E/YNbJZUc8E7HzJYAJwHvu/tJZtYM+LO7fzceBYwl3emIpJYXXpjOJX8c\ny5Qh4xgwoG+ii5OxYj05dIe77wP2mlk28CXl+15ERGJq0qSn6dTp+4we/Ta7d/dl1Ki36NTp+0ya\n9HSiiyY1VJ15Ou+aWSPgj4T7Sf4H/DOmpRIRidChTys67s1m5sy50OUfFL9/Bn375mnVgRRUZdAx\ns6aEVxVwd3/EzGYAR7v7B3EpnYgI4VUHNjbbxYw5M2jSYB+75nRj4GX9tepACqpqcuhVwB3Ap8Cx\nZjbE3afFrWQiIhG06kB6qHQggZl9DITcfb2ZHQtMcfdT41q6GNNAApHUU7iyUJNAE+xgBhJUFXQW\nuXvnyt6nAwUdEZGai9U8ndZm9iD7F+dsFfHe3X14bS4oIiKZq6qg838Eqw8EFgbvrUK+iIhItejJ\noRlcf5FkoxUHUkNMJocGT/3sXsX2U8xscm0uKiISzdSpM5i68ANefLHGawlLiqhqIMEJhJvYTgX+\nTfgxBAY0B/4f4Qmi97j7kvgUte7pTkckOUya9DQPPvgse/acxPJW9WlftJdDD/2A4cMHcc01lya6\neFJBTAYSuPtHwGVm1gDoDLQj3JezCvjA3XfW5oIiIhVpxYHMccBlcNx9F+Enes470L4iIrWhFQcy\nR3XWXhMRiTmtOJAZNHotg+svkoy04kDyi8mKBJlAQUdEpOZi/TydaBe8pjbHiYhIZqtV0BEREakN\nNa9lcP1FRGojps1rZnadmWVb2GNmtsjM9HByERGpseo0r13p7luAPkAuUADcFdNSiUhacXd+PPoa\n1LIg1Qk6pbdQ5wBPpfKyNyKSGFpTTUodsE/HzB4HWgLHAicBhwBz3L1rzEsXY+rTEYktramWnmL1\nELdSVxJee+1Td/+fmTUGflKbi4lIZtGaalJRdYLOaYQX+NxmZgVAF+D+2BZLRNKB1lSTiqoTdB4B\nTjSzk4DrgUeBJ4FesSyYiKQHrakmkarTp7PI3Tub2RigyN0fNbP33b1LfIoYO+rTEYkframWPmK6\n9pqZvQVMJ9yP8x1gPbDY3U+ozQWTiYKOiEjNxXrttYHATsLzdb4AWgF31+ZiIpK+NBdHquOAQcfd\n1wEvAg2CrA3AX2NZKBFJPZqLI9VRnea1IcDVQK67f8PMOgAPu/tZ8ShgLKl5TeTgaS5O5on1PJ1r\ngR4Ej6t292VmdkxtLiYi6UdzcaQmqhN0drn7LrNwUDOz+oBuD0QE0FwcqZnqBJ25ZnYTcISZnQ0M\nA16JbbFEJJVoLo5UV3X6dOoBVxFeZRpgBvBoOnSGqE9HpG5pLk5miNk8naApbYm7f7O2hUtmCjoi\nIjUXs3k67r4X+LeZtatlwf5kZsVm9lFEXq6ZzTKzZWY208xyIraNMrPlZvaJmfWJyO9qZh8F2x6I\nyG9gZn8J8udFltPMLg+usczMLqtN+UVkP83DkbpQncmhucDHZvammb0SvF6u5vknA/0q5I0EZrl7\nB2B28B4z60h4ImrH4JiJVjp6AR4GBrt7e6C9mZWeczCwMcj/LTA+OFcucDPhUXc9gDGRwU1Eak7z\ncKQuVKdPp3Rhz8hbKXf3udW6gFke8Erpsjlm9gnQy92Lzaw5UOju3zSzUUCJu5cGjunAWGAV8Ka7\nHx/kDwJC7v7TYJ8x7j4/aApc5+5NzexHQE93Hxoc80hwnWcrlE3NayIHoHk4UlFM5umYWUPgp8Bx\nwIfAn9x9T+2KWE4zdy8O0sVAsyDdkmAuUGAN4SV39gTpUkVBPsHP1RBuCjSzLcHzflpWOGZNxDEi\nUgOahyN1qaoh008Au4G3ge8RbvYaUZcXd3c3M91qiCQxzcORulRV0Dk+oknsMeDdOrpmsZk1d/cv\nzKwF8GWQXwS0idivNeE7lKIgXTG/9Ji2wNqgeS3b3TeaWREQijimDfBmtMKMHTu2LB0KhQiFQtF2\nE8lomoeT2QoLCyksLKyTc1Xap1P6HJ3K3lf7Al/v05lAuPN/vJmNBHLcfWQwkGAK4Y7/VsAbwHHB\n3dB8YDiwAPgb8KC7TzezYcAJ7j406Os5390HBQMJ3iP8lFMDFgJd3H1zhbKpT0ekBjQPRyBG83TM\nbB+wPSKrIbAjSLu7H12Ngj1D+AmjTQj339wMTAOeI3yHshK4uDQYmNlo4EpgLzDC3WcE+V2Bx4My\nvObuw4P8BsBTQGdgIzDI3VcG234CjA6Kcpu7PxGlfAo6IiI1FNOHuKUzBR0RkZqL9UPcRCSNadKn\nxJOCjkiG06RPiSc1r2Vw/SWzadKn1FasH+ImImlIkz4lERR0RDKUJn1KIijoiGQwTfqUeFOfTgbX\nX6SUJn1KTWieTi0p6IiI1Jzm6YiISEpQ0BFJU5r0KclIQUckTWnSpyQj9elkcP0lPWnSp8SaJoeK\nSBlN+pRkpqAjkmY06VOSmYKOSBrSpE9JVurTyeD6S/rTpE+JBU0OrSUFHRGRmtPkUJEMozk4kqoU\ndERSkObgSKpS81oG119Sj+bgSDLQPB2RDKE5OJLqFHREUojm4EiqU9ARSTGagyOpTH06GVx/SW2a\ngyOJonk6taSgIyJSc5qnI5LiNO9GMoWCjkgS0LwbyRRqXsvg+kviad6NpCLN0xFJUZp3I5lGQUck\ngTTvRjKNgo5IgmnejWQS9elkcP0luWjejaQKDZkWSUI1HQatgCOZQEFHJEY0DFrk69S8lsH1l9jQ\nMGhJdxoyLZJENAxapHIJa14zs5Vm9qGZLTKzBUFerpnNMrNlZjbTzHIi9h9lZsvN7BMz6xOR39XM\nPgq2PRCR38DM/hLkzzOzdvGtoWSq3vm9GdisAJ/TnSYfn4rP6cbAZgUaBi1CYvt0HAi5e2d37xHk\njQRmuXsHYHbwHjPrCAwEOgL9gIlmVnpr9zAw2N3bA+3NrF+QPxjYGOT/Fhgfj0qJwP5h0MOG9mHy\n5P4aBi0SSFifjpl9DnRz940ReZ8Avdy92MyaA4Xu/k0zGwWUuPv4YL/pwFhgFfCmux8f5A8iHMh+\nGuwzxt3nm1l9YJ27N61QBvXpSExpGLSko1QdMu3AG2b2npldHeQ1c/fiIF0MNAvSLYE1EceuAVpF\nyS8K8gl+rgZw973AFjPLrfNaSMaozUrQCjgi5SVyIMEZ7r7OzJoCs4K7nDLu7mYW89uQsWPHlqVD\noRChUCjWl5QUVToEesCLMxkwoG+iiyMSN4WFhRQWFtbJuZJiyLSZjQG2AVcTbh77wsxaAHOC5rWR\nAO5+V7BVQrjqAAAJoUlEQVT/dGAM4ea1ORHNaz8Cerr70NImOHefp+Y1ORgaAi1SXso1r5nZEWaW\nFaSPBPoAHwEvA5cHu10O/DVIvwwMMrPDzCwfaA8scPcvgK/M7JRgYEEBMC3imNJzXUh4YIJIjXXo\n04qOw7IpPn4uhG6h+Pi5dLo2R0OgRWohUc1rzYCXggFo9YE/u/tMM3sPeM7MBgMrgYsB3H2pmT0H\nLAX2AsMiblGGAY8DDYHX3H16kP8Y8JSZLQc2AoPiUTFJP1oJWqTuJEXzWqKoeU2q6847/0iHDm35\nMPefnPjf8ErQI0delehiiSTEwTSvKehkcP2l5jQEWkRBp9YUdEREai7lBhKIJEpt5tqISN1R0JGM\noscNiCSWmtcyuP6ZRHNtROqOHm0gcgB63IBIclDQkYyguTYiyUFBRzJG6eMGPsw9mhMHnq7HDYgk\ngPp0Mrj+mUpzbUQOjubp1JKCjohIzWmejqQ9za8RSQ8KOpISNL9GJD2oeS2D658KNL9GJPlono6k\nLc2vEUkvCjqS1DS/RiS9KOhI0tP8GpH0oT6dDK5/qtH8GpHkoCHTklRiNbxZAUck9SnoSJ3T8GYR\nqYya1zK4/nVNw5tFMoOGTEtS0PBmETkQBR2pMxreLCIHoqAjdUrDm0WkKurTybD6uzuX3PRT/nz7\nI5jVqkm2WjS8WSR9aci0VFu8RpYp4IhINLrTyZD6a2SZiNQVjV6TA9LIMhFJBgo6GUIjy0QkGSjo\nZBCNLBORRFOfTpzqH69RY9WhkWUicjA0ei0FJNN6ZAo4IpIoutOJcf01akxE0o1GryUxjRoTEdlP\nQSfGNGpMRGQ/BZ040KgxEZGwtO7TMbN+wP3AIcCj7j6+wva4rkigUWMikg40ei0KMzsE+B3QD+gI\n/MjMjk9kmeIdcAoLC+N6vXhL5/qlc91A9ctkaRt0gB7ACndf6e57gGeB8xJcprhK91/8dK5fOtcN\nVL9Mls5BpxUQ2XmyJsgTEZEESeegk76dVSIiKSptBxKY2anAWHfvF7wfBZREDiYws/SsvIhIjNV2\nIEE6B536wL+Bs4C1wALgR+7+r4QWTEQkg6XtPB1332tmPwNmEB4y/ZgCjohIYqXtnY6IiCSftB1I\nYGZtzGyOmX1sZkvMbHiQ38PMFpjZIjN718y6RxwzysyWm9knZtYncaU/MDM73Mzmm9liM1tqZncG\n+blmNsvMlpnZTDPLiTgmHep3t5n9y8w+MLMXzSw74piUr1/E9hvMrMTMciPyUqJ+VdXNzH4efH5L\nzCyyfzUl6gZV/m6mxXdLKTM7JKjLK8H7uvlucfe0fAHNgZOD9FGE+3eOBwqBvkF+f2BOkO4ILAYO\nBfKAFUC9RNfjAHU8IvhZH5gHfBuYAPwqyL8RuCvN6nd2abmBu9KtfsH7NsB04HMgNxXrV8ln1xuY\nBRwabGuainWron5z0uW7JSj39cCfgZeD93Xy3ZK2dzru/oW7Lw7S24B/EZ6nsw4o/es4BygK0ucB\nz7j7HndfSfgfrkdcC11D7r49SB5GuN9qE/AD4Ikg/wng/CCdDvX7r7vPcveSIH8+0DpIp0X9gvf3\nAb+qsHtK1a+S382fAnd6eLI27r4+2Cel6gaV1u8L0uS7xcxaA98DHgVKR6nVyXdL2gadSGaWB3Qm\n/BfJSOBeM/sPcDcwKtitJeEJpKWSfjKpmdUzs8VAMeG/qj4Gmrl7cbBLMdAsSKdD/ZZW2OVK4LUg\nnRb1M7PzgDXu/mGF3VOqfpX8bnYAeprZPDMrNLNuwe4pVTeotH5p890C/Bb4P6AkIq9OvlvSPuiY\n2VHAC8CI4I7nMWC4u7cFfgH8qYrDk3qUhbuXuPvJhP/a72lmvStsd6quQ6rVL1S6zcxuAna7+5Sq\nThHjIh6UKPX7HuEvqjERu1U1FyJp61fJZ1cfaOTupxL+QnuuqlPEvpS1V0n90uK7xcy+D3zp7ouo\n5PfvYL5b0jromNmhwFTgaXf/a5Ddw91fCtIvsP82sIhwW3qp1uy/PU5q7r4F+BvQFSg2s+YAZtYC\n+DLYLR3q1w3AzK4gfOt/ScRu6VC/LkA+8IGZfU64DgvNrBkpWr8Kn90a4MUg/12gxMyakKJ1g6/V\nL12+W04HfhD8Dj4DnGlmT1FX3y2J7qyK1YtwhH4S+G2F/PeBXkH6LOBdL98Zdhjh//ifEgwpT8YX\n0ATICdINgbeC+kwAbgzyR/L1zr5Ur18/4GOgSYX906J+FfaJNpAg6etXxWd3DTAuyO8A/CfV6lZF\n/b6bLt8tFeraC3glSNfJd0vaTg4FzgAuBT40s0VB3mhgCPB7M2sA7Aje4+H29OeApcBeYJgH/6JJ\nqgXwhJnVI3zH+pS7zw7q+pyZDQZWAhdDWtVvOeFf7llmBvCOuw9Ll/pV2Kes/ClWv8o+u7eAP5nZ\nR8Bu4DJIubpB9Pq9YWbp8t1SUWlZ76IOvls0OVREROImrft0REQkuSjoiIhI3CjoiIhI3CjoiIhI\n3CjoiIhI3CjoiIhI3CjoiBwkM2scLAG/yMzWmdmaIP2+hZ9gG8tr32NmvYJ0oZl1DdL5wRL0Z5vZ\niWb2WCzLIVJd6Tw5VCQu3H0j4QVlMbMxwFZ3vy/W1zWzLKCnu/+ytCiABysEvw5c7+6zgn2/YWbH\nuPuXlZxOJC50pyNS98zMugZ3Hu+Z2fSINasKzey+4CFf/zKz7mb2UnBXcmuwT17wMKyng4eEPW9m\nDaNc5zzgjQp5rQg/on20u78akf86cFEM6ipSIwo6InXPgAeBC929GzAZuD3Y5sAud+8OPAxMI/yc\nmW8BV5hZo2C/DsDv3b0j8BUwLMp1zgDeq3Ddx4GH3P3FCvsuAHoeZL1EDpqCjkjda0A4iMwK1sK7\nifLPF3k5+LkEWOLuxe6+G/iM/av1rnb3d4L004SfTFlRO8IPJSzlhO98CqLcGa0j/FRHkYRSn45I\n3TPgY3c/vZLtu4KfJRHp0vel/ycjF0U0Kn8+ScU/HCcABcDzZnaeu++rxjlE4kZ3OiJ1bxfQ1MxO\nhfBzncysYw3P0bb0eODHwNtR9lkFNK+Q5+5+HeEmucgRay2C/UUSSkFHpO7tAy4ExgePNF4EnBZl\nv6qevvhv4FozWwpkE+7/qejvBA+2i+JyoIWZjQ/e9yD83BeRhNKjDUSSjJnlEX5w1gkH2O8oYE4w\nKOFA5ywELtaQaUk03emIJKcD/jXo7tuAOWbWu6r9zOxEYIUCjiQD3emIiEjc6E5HRETiRkFHRETi\nRkFHRETiRkFHRETiRkFHRETiRkFHRETi5v8DHTvNaTD2Mp8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe928d80c90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(X,y,'b*',label='Experimental')\n",
    "plt.plot(X,y,'g+',label='Calculated')\n",
    "plt.legend(loc='best')\n",
    "plt.xlabel('Temp (K)')\n",
    "plt.ylabel('Press. (Pa)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The experimental and calculated data points nearly overlap!\n",
    "\n",
    "The plot of the residuals shows more diagnostics- the fact that these residual points are not randomly distributed indicates the equation we are fitting is not ideal:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x7fe9260b1250>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAEZCAYAAACjPJNSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XucV3W97/HXW0Wz7RVrey85hSWE10R32XY6FeKWoybl\nJS1TE43c2I5tMmSbsW1iuqWis7XZpXjhqPlQMjVD0AezsR5bKEUFER084QEV0hLFlC3K5/yxvgM/\nf8wMM8xv1vpd3s/H4/dgre+6fdfiN/OZ73d9L4oIzMzM+ttWRWfAzMwagwOOmZnlwgHHzMxy4YBj\nZma5cMAxM7NcOOCYmVkuHHDMzCwXDjhWcyQtk/SGpDWSVkq6WdJOJdtvkLRe0vFlx/0wpZ+Z1reV\ndLWk5elcf5T0wy6u0/GZugX5HSNpqaRXJf1e0idLtm0n6fq07UVJ/1R27NaSLpP0vKTXJD0qaecu\nrrOTpOmSXkqf6ZJ27GLfpvQs1qTzLpH01d7em1lvOOBYLQpgVETsCBwEDAMuKdv+DPCVjgRJ2wAn\nA0vTdoBm4FDg8HSuJuDRzq5T8hnXm4xKOhi4GvhiROwMXAf8UpLSLi3Ah4APAJ8Gvi3pmJJTXAoc\nCRwZETsBZwBru7hcC/A+YFA65+4prSvPp3vaCbgY+JmkA3pzf2a94YBjNS0iVgGzgKFlm+4BjpK0\nS1ofCTwOrCrZ5+PAXRGxMp3ruYi4ucJZHAIsjogFaf1msqDwt2n9K8C/RsSrEbEE+A/gqwCSdgUu\nBM6NiOUpj4sj4r+7uNbQdD+vR8RrwF1s+lw6FRG/Al4BDpB0nKQFqdT1/yRN6uU9m3XKAcdqlQAk\n7UMWTOaVbV8L/Ao4Na1/BbipbJ+HgW9J+rqkYSWljk2us0midJSkV7r5fCLt+hAwSNJwSVsDZwML\nImJVCih7kgXCDk+wMUgMA94Gvpiq256WNLabZ3I/MFrSLunco4H7utm/4162kvR5YBdgIfA6cEYq\nkR0HfF3SCZs7j9nmbFN0Bsy2gIC7JAWwA1lguayT/W4CrpJ0K/D3ZEHnGyXbJ5P9VX868EPgz5Ka\nI6IjMHVc5+2SY/45Iq6LiN8Cu24uoxGxXNIlwO9S0ivAP6TlHdK/r5Yc8hrQ8d5lH2BnYDCwH7A/\n8KCkZyLigU4u9+/AscCf0/oDwLXdZG8vSa8A64HnyIJMO9Bekv+Fkm4DjiZ7zmZbzCUcq0UBnJDe\nPTQB/5Oseuxd+0TE74D3k73fuSci1pbtsD4iromIo8h+sX8fuF7SR8qus2vJ57reZDQ1XBgPHBAR\nA4AvA/dK2oOsJAGwU8khOwNr0vKb6d/vRcR/R8RC4DY2Bqxy/wd4miyQ7QT8X2B6N9l7Id3TbhFx\naETcnvJ8hKQ5kv4kaTVwHrBbL27brFMOOFbTImIu8BPgB13sMh34FptWp5Wf578j4hqyEsiQzV1X\n0qfKWq+Vfzpaoh0D/Doilqbr3A+8CHwiIl5JyweXnPogYFFafqKr7HaRPhJojYg3I+KvQCtdB6fu\n3EL2/mefiNgF+Cn+XWEV4C+R1YMfAcMlHZHWxcZ3L1OBz0bEQ+UHSfqmpKMlbS9pm9RcegdgQelu\nnV0wIh4qa71W/umoQnscOE7SIGU+R1Y11hFUbgIuSe9dDgC+BtyQrvEs2Tug76Qm3AcApwD3dvEc\nngDOlfQeSdsDY3j3+6Ge2gF4JSLekjQc+BJdBzmzHnPAsZoXES8DN5I17YXsl2Okba9ExJwuDv0r\nWZPlF4GXgK8DoyNiWck+95SVXO7sZfZ+TvbuYy7Zu5ofAWMi4pm0fRLwLNk7lDnADyJiVsnxpwEf\nJHsvcy9wScf9SDpd0qKSfb9KFsyeB1aQvfc5s5u8dRVExgLfk/Qa8F3gFz25UbPNUVETsEl6D/Cf\nwHbAtsCvIqJZ0kCyL/gHgWXAyRGxOh3TTNbK5x1gXNkPppmZVbHCAg6ApPdGxBupU95vgX8Gjgde\njogrJV0M7BoREyQNIatbPhzYm6wFzv4Rsb6o/JuZWc8VWqUWEW+kxW2Brcle2B5PVj1C+vfEtHwC\ncGtErEtVHkuB4fnl1szM+qLQgJM6nD1G1vt7TkQ8Ceyeeo+T0ndPy3uR1Ut3WEFW0jEzsxpQaMfP\nVB12cBqM8H5Jny7bHqlzX5en6NcMmplZxVTFSAMR8aqkXwOHAask7RERKyXtCfwp7fY8sG/JYfuk\ntHfZTIAyM7MuRESn3QAqpbAqNUnv6xhYMfUZ+BxZ/4e72diU80yyDmik9FNTf4RBZMN9zO/s3BHh\nTwSTJk0qPA/V8vGz8LPws+j+k4ciSzh7AjdK2oos8N0cEQ9KWgDcLukcUrNoyEbJlXQ7sJhsQMOx\nkddTMjOzPiss4EQ2LtShnaT/BfhsF8dcDlzez1kzM7N+4JEG6lhTU1PRWagafhYb+Vls5GeRr0I7\nfvYHSa5pMzPrJUlEvTYaMDOzxtKwAScimDDhytxaZ5iZNbqGDTh33nk/11zzIjNmePxPM7M8NFzA\naW2dztCho5g48SHWrJlCc/Nchg4dRWtrdxMjmplZX1XFSAN5GjPmdAYO3I3x4+cCYu3a9Vx++QWM\nHn1M0VkzM6trDVfCkYQkVq9ey5Ah32L16jc3pJmZWf9puBIOQHv7cqZNG8lJJ41gxoxZtLcvLzpL\nZmZ1z/1wzMy6ERE0N1/F5MkX1XVNiPvhmJkVzC1aK8cBx8ysE27RWnkN+Q7HzGxz3KK18lzCMTPr\nRH+2aG3UkU4ccMzMutDRonXRoquZNu3YirVobdT3Qm6lZmaWk9bW6Uydehvr1h1Ee/tlDB58CQMG\nPM64cady3nlnFJq3PFqp+R2OmVlOGv29kKvUzMxy0ugjnbiEY2aWo0Ye6cTvcMzMzCMNmJlZ/XDA\nMTOzXDjgmJlZLgoLOJL2lTRH0pOSFkkal9IHSpot6RlJsyTtUnJMs6R2SUskjSgq72Zm1nuFNRqQ\ntAewR0Q8JmkH4BHgROAs4OWIuFLSxcCuETFB0hDgFuBwYG/gAWD/iFhfdl43GjAz66W6bjQQESsj\n4rG0/DrwFFkgOR64Me12I1kQAjgBuDUi1kXEMmApMDzXTJuZ2Rarinc4kvYDDgHmAbtHxKq0aRWw\ne1reC1hRctgKsgBlZmY1oPCOn6k67U7gwohYU9rjNiJCUnf1Y51ua2lp2bDc1NREU1NTRfJqZlYv\n2traaGtry/WahXb8lDQAuBf4TUT8KKUtAZoiYqWkPYE5EfFRSRMAIuKKtN9MYFJEzCs7p9/hmJn1\nUl2/w1FWlLkOWNwRbJK7gTPT8pnAXSXpp0raVtIgYDAwP6/8mplZ3xTZSu0oYC7wBBurxprJgsjt\nwAeAZcDJEbE6HTMROBt4m6wK7v5OzusSjplZL+VRwvFYatatiKC5+SomT76oYUa0NWtEdV2lZrWh\nUWcmNLPKc8BpQD2ZT721dTpDh45i4sSHWLNmCs3Ncxk6dBStrdNzzKmZ1RMHnAbUk1LLmDGn09Ly\nDdauXU/HzISXXnoBY8acnl9GzayuOOA0kN6UWhp9ZkIzq7zCO35afno7n3ojz0xoZpXngNNAykst\ny5ev77bU0tx87oblroJSZ9yyzcw64yq1BtNRalm06GqmTTu2X0otbtlmZp1xPxyrmNbW6Uydehvr\n1h1Ee/tlDB58CQMGPM64cady3nlnFJ09M+tGHv1wXKVmFdPbd0Rm1lhcpWYV45ZtZtYdl3Csotyy\nzcy64nc4Vii3aDOrDh5LzeqeW7SZNQ4HHCuEx2ozazwOOHWkJ4NyVguP1WbWeBxw6kgtVU+5RZtZ\n43HAqQO1Wj2Vx6gHZlY93EqtDkQEd9wxk/Hj57J8+WT23beZKVOOZvToY1xiMLMecSs16xFXT5lZ\nLXDHzzrhDpdmVu1cpWZmZq5SMzOz+uGAY2ZmuSg04Ei6XtIqSQtL0gZKmi3pGUmzJO1Ssq1ZUruk\nJZJGFJNrMzPbEkWXcKYBI8vSJgCzI2J/4MG0jqQhwCnAkHTMNZKKzr+ZmfVQob+wI+Ih4JWy5OOB\nG9PyjcCJafkE4NaIWBcRy4ClwPA88mlmZn1XjSWE3SNiVVpeBeyelvcCVpTstwLYO8+MmZnZlqvq\nfjgREZK6a+Pc6baWlpYNy01NTTQ1NVU2Y2ZmNa6trY22trZcr1l4PxxJ+wH3RMSwtL4EaIqIlZL2\nBOZExEclTQCIiCvSfjOBSRExr+x87odjZtZLjdoP527gzLR8JnBXSfqpkraVNAgYDMwvIH9mZrYF\nCq1Sk3QrcDTwPknLgX8BrgBul3QOsAw4GSAiFku6HVgMvA2MdVHGzKx2FF6lVmmuUjMz671GrVIz\nM7M65IBjNaeWptI2s40ccKzm1NJU2ma2kQOO1YxanUrbzDJV3fHTrNSYMaczcOBujB8/FxBr167n\n8ssvYPToY4rOmpn1gEs4VjM8lbZZbXMJx2qKp9I2q13uh2NmZu6HY2Zm9cMBp4q5v4mZ1RMHnCrm\n/iZmVk8ccKqQ+5uYWT1yK7Uq5P4mZlaPXMKpQu5vYmb1yCWcKuX+JmZWb9wPx8zM3A/HzMzqhwOO\nmZnlwgHHzMxy4YBjZma5cMAxM7NcOOCYmVkuHHCsYXgwVLNi1VzAkTRS0hJJ7ZIuLjo/Vjs8GKpZ\nsbrs+ClpfDfHRURM6Z8sdU3S1sDTwGeB54HfA6dFxFMl+7jjp71La+t0pk69jXXrDqK9/TIGD76E\nAQMeZ9y4UznvvDOKzp5ZVcij42d3Q9vsCHT2m1tdpOdhOLA0IpYBSLoNOAF4qruDrLF5MFSz6tBl\nwImIlhzz0VN7A6WDiq0AjigoL1YjygdDXb58vQdDNSvAZgfvlLQ9cA4wBNieVLqJiLP7N2ud6lHJ\nqqWlZcNyU1MTTU1N/ZQdqxUeDNXs3dra2mhra8v1mpsdvFPSHWRVVqcDlwJnAE9FxLj+z94meTkS\naImIkWm9GVgfET8o2cfvcMzMeimPdzg9CTiPRcTBkp6IiAMlDQB+GxG5V2VJ2oas0cBngBeA+bjR\ngJlZnxXdaKDDW+nfVyUNA1YC7++/LHUtIt6WdAFwP7A1cF1psDEzs+rVkxLOucCdwDDgBmAH4LsR\n8dN+z90WcAnHzKz3qqJKrdY44JiZ9V5VVKlJmlSyuuE3eUR8r19yZGZmdakn73D+ysZAsz0wCljc\nbzkyM7O61OsqNUnbAbMi4uj+yVLfuErNzKz38qhS25LBO/+GrMe/mZlZj/XkHc7CktWtgL8F/P7G\nzMx6pSfNovcrWX0bWBUR6/oxT33iKjUzs94rtFm0pIHdHRgRf+mXHPWRA46ZWe8V/Q7nUeCR9O/L\nQHv6vJzSzeqaZwg1q6wuA05E7BcRg4DZwKiI2C0idgOOS2lmdc0zhJpVVk/e4SyKiI9tLq1auErN\n+sozhFojqoqRBoAXJF0CTCeb7fNLZNM7m9UlzxBq1j960g/nNLKm0L8EZqTl0/ozU2ZFKp8hdPXq\nNz1DqFkFbDbgRMSfI2JcRBySPhdWaws1s0rpmCF00aKrmTbt2B7PEOqGBmZd665Z9I8j4kJJ93Sy\nOSLi+P7N2pbxOxwr0h13zOTss+9n2rSRroKzmlJ0P5zDIuIRSU2dbI6I+M/+zNiWcsCxIrihQXWJ\nCJqbr2Ly5ItcFdpDhTYaiIhH0r9tJRkaCOwTEU/0Z6bMao0bGlSXjibthx8+y/8HVWSz73AktUna\nKQWbR4CfS/ph/2fNrHa4oUF1aG2dztCho5g48SHWrJlCc/Nchg4dRWvr9KKzZvSsWfQuEfGapK8B\nN0XEpLIBPc2MjQ0NTjppBDNmzOpxQwOrHJc0q1tPAs7WkvYETgYuSWl+SWJWprn53A3L/gVXjPKS\n5vLl613SrCI96YfzPeB+4NmImC/pQ2RjqpmZVZ0tbdJu/a/XM35Wu2pvpebWM2ZWjYoeLbojEx+R\n9KCkJ9P6gWmoG9sCHhDSzBpVT6rUfgZMBN5K6wvp49A2kr4o6UlJ70g6tGxbs6R2SUskjShJP0zS\nwrTtx325fhHcesbMGl1PAs57I2Jex0qqr+rrjJ8Lgc8Dc0sTJQ0BTgGGACOBa7Sx3ula4JyIGAwM\nljSyj3nI1Zgxp9PS8g3Wrl1PR+uZSy+9gDFjTi86a2ZmuehJwHlJ0oc7ViR9AXixLxeNiCUR8Uwn\nm04Abo2IdRGxDFgKHJFaye0YEfPTfjcBJ/YlD3lzPw0za3Q9aRZ9AfAfwEckvQD8EeivP8v3Ah4u\nWV8B7E1WolpRkv58Sq8p7qdhZo1sswEnIp4FPiNpB7L5cF4n65OzrLvjJM0G9uhk08SI6GxA0Ipp\naWnZsNzU1ERTU1N/Xq7H3E/DzKpFW1sbbW1tuV6zu8E7dwDOAz4ELAJ+Slbl9X1gaSVGi5Y0Bxgf\nEY+m9QkAEXFFWp8JTAKeA+ZExAEp/TTg6Ig4v5NzVnWzaDOzalR0s+ibgGHA48BnyKq6/gn4UoWn\nJii9wbuBUyVtK2kQMBiYHxErgdckHZEaEXwZuKuCeTAzs37WXQnniYg4MC1vTdZQ4IMR8WafLyp9\nHpgKvA94FVgQEcembROBs4G3gQsj4v6UfhhwA7A9cF9EjOvi3C7hmJn1UtHz4SyIiEO6Wq9WDjhm\ntcsjcRSn6Cq1AyWt6fgAw0rWX+vPTJlZY/JIHPXNY6mZWeE8Y2rxCp3x08wsL57HpjH0ZKQBM7N+\n5ZE4GoNLOGZWFTwSR/3zOxwzMyu8lZqZ1ZiIYMKEK/EfXVaNHHDM6oibFVs1c8AxqwOe4M9qgQOO\nWZXZkmoxT/BntcABx6zKbEm1mJsVWy1wwDGrEn2tFutoVrxo0dVMm3asmxVb1XGzaLMqERHcccdM\nxo+fy/Llk9l332amTDma0aOPcUnF+p2bRZs1EFeLWb3zSANmVcS97a2euUrNzMxcpWZmZvXDAcfM\nzHLhgGNmZrlwwDEzs1w44JhZRXikatscBxwzqwiPVG2b44BjZn3ikaqtpwoJOJKukvSUpMclzZC0\nc8m2ZkntkpZIGlGSfpikhWnbj4vIt1k929IqMY9UbT1VVAlnFjA0Ig4CngGaASQNAU4BhgAjgWu0\ncVyPa4FzImIwMFjSyPyzbVa/trRKzEPyWE8VEnAiYnZErE+r84B90vIJwK0RsS4ilgFLgSMk7Qns\nGBHz0343ASfmmWezelWJKjGPVG09UQ1jqZ0N3JqW9wIeLtm2AtgbWJeWOzyf0s2sj8aMOZ2BA3dj\n/Pi5dFSJXX75BYwefUyPz9HcfO6G5d4cZ42l3wKOpNnAHp1smhgR96R9vgO8FRG3VPLaLS0tG5ab\nmppoamqq5OnN6kp5ldjy5etdJdYA2traaGtry/WahQ3eKemrwLnAZyJibUqbABARV6T1mcAk4Dlg\nTkQckNJPA46OiPM7Oa8H7zTrpcmTf8b++3/gXaNUT5jwtaKzZTnKY/DOQgJOeuF/NVnQeLkkfQhw\nCzCcrMrsAeDDERGS5gHjgPnAr4GpETGzk3M74JiZ9VIeAaeodzg/AbYFZqdi+39FxNiIWCzpdmAx\n8DYwtiR6jAVuALYH7uss2JiZWfXyfDhmZub5cMzMrH444JiZWS4ccMzMLBcOOGZmlgsHHDMzy4UD\njpmZ5cIBx8zMcuGAY2ZmuXDAMTOzXDjgmJlZLhxwzMwsFw44ZmaWCwccMzPLhQNOBUQEEyZciUep\nNjPrmgNOBdx55/1cc82LzJgxq+ismJlVLQecPmhtnc7QoaOYOPEh1qyZQnPzXIYOHUVr6/Sis2Zm\nVnWKmvGzLowZczoDB+7G+PFzAbF27Xouv/wCRo8+puismZlVHZdw+kASkli9ei1DhnyL1avf3JBm\nZmbv5hJOH7W3L2fatJGcdNIIZsyYRXv78qKzZGZWlVRvLaskRb3dk5lZf5NERPRr9Yyr1MzMLBcO\nOGZmlgsHHDMzy0UhAUfSv0p6XNJjkh6UtG/JtmZJ7ZKWSBpRkn6YpIVp24+LyLeZmW25QhoNSNox\nItak5X8EDoqIr0kaAtwCHA7sDTwADI6IkDQfuCAi5ku6D5gaETM7ObcbDZiZ9VLdNhroCDbJDsDL\nafkE4NaIWBcRy4ClwBGS9gR2jIj5ab+bgBPzyq+ZmfVdYf1wJH0f+DLwJjA8Je8FPFyy2wqyks66\ntNzh+ZRuZmY1ot8CjqTZwB6dbJoYEfdExHeA70iaAPwIOKtS125padmw3NTURFNTU6VObWZWF9ra\n2mhra8v1moV3/JT0AeC+iPhYCj5ExBVp20xgEvAcMCciDkjppwFHR8T5nZzP73DMzHqpbt/hSBpc\nsnoCsCAt3w2cKmlbSYOAwcD8iFgJvCbpCGUDlX0ZuCvXTJuZWZ8U9Q5nsqSPAO8AzwJfB4iIxZJu\nBxYDbwNjS4orY4EbgO3JSkSbtFAzM7PqVXiVWqW5Ss3MrPfqtkqtVnjqaDOzynHA6YanjjYzqxwH\nnE546mgzs8rzBGyd8NTRZmaV5xJOJzx1tJlZ5bmE0wVPHW1mVlluFm1mZm4WbWZm9cMBx8zMcuGA\nY2ZmuXDAMTOzXDjgmJlZLhxwzMwsFw44ZmaWCwccMzPLhQOOmZnlwgHHzMxy4YBjZma5cMAxM7Nc\nOOCYmVkuHHDMzCwXDjhmZpaLQgOOpPGS1ksaWJLWLKld0hJJI0rSD5O0MG37cTE5NjOzLVVYwJG0\nL/A54LmStCHAKcAQYCRwjTbO63wtcE5EDAYGSxqZc5ZrTltbW9FZqBp+Fhv5WWzkZ5GvIks4U4Bv\nl6WdANwaEesiYhmwFDhC0p7AjhExP+13E3BibjmtUf5h2sjPYiM/i438LPJVSMCRdAKwIiKeKNu0\nF7CiZH0FsHcn6c+ndDMzqxHb9NeJJc0G9uhk03eAZmBE6e79lQ8zM6sOioh8Lyh9DHgQeCMl7UNW\nYjkCOAsgIq5I+84EJpG955kTEQek9NOAoyPi/E7On+8NmZnViYjo1z/+cw84m2RA+iNwWET8JTUa\nuAUYTlZl9gDw4YgISfOAccB84NfA1IiYWVS+zcysd/qtSq0XNkS8iFgs6XZgMfA2MDY2RsSxwA3A\n9sB9DjZmZrWl8BKOmZk1hroZaUDSyNRZtF3SxUXnp5IkLZP0hKQFkuantIGSZkt6RtIsSbuU7N+r\nzrOStpP0i5T+sKQP5nuHXZN0vaRVkhaWpOVy75LOTNd4RtJX8rjf7nTxLFokrUjfjQWSji3ZVpfP\nQtK+kuZIelLSIknjUnrDfS+6eRbV+b2IiJr/AFuT9dnZDxgAPAYcUHS+Knh/fwQGlqVdCXw7LV8M\nXJGWh6T7H5Cex1I2lmTnA8PT8n3AyLQ8FrgmLZ8C3Fb0PZfc56eAQ4CFed47MBB4FtglfZ4FdqnC\nZzEJ+FYn+9btsyBr/XpwWt4BeBo4oBG/F908i6r8XtRLCWc4sDQilkXEOuA2sk6k9aS89cjxwI1p\n+UY2doTdks6zpee6E/hM5bO/ZSLiIeCVsuQ87v0YYFZErI6I1cBsstEvCtPFs4DOuxXU7bOIiJUR\n8Vhafh14iqyRUcN9L7p5FlCF34t6CTh7A8tL1js6jNaLAB6Q9AdJ56a03SNiVVpeBeyelrek8+yG\n5xcRbwOvqmR8uyrU3/e+Wzfnqkb/KOlxSdeVVCM1xLOQtB9ZqW8eDf69KHkWD6ekqvte1EvAqfeW\nD5+MiEOAY4FvSPpU6cbIyrf1/gw61cj3nlwLDAIOBl4Eri42O/mRtAPZX9wXRsSa0m2N9r1Iz+IO\nsmfxOlX6vaiXgPM8sG/J+r68O/LWtIh4Mf37EvBLsirEVZL2AEjF4T+l3cufxT5kz+L5tFye3nHM\nB9K5tgF2joi/9MvNVEZ/3/ufOzlXVX6nIuJPkQA/J/tuQJ0/C0kDyILNzRFxV0puyO9FybOY3vEs\nqvV7US8B5w9kI0jvJ2lbshdbdxecp4qQ9F5JO6blvyEbEmgh2f2dmXY7E+j4obsbOFXStpIGAYOB\n+RGxEnhN0hGSBHwZ+FXJMR3n+gLZSBDVLI97nwWMkLSLpF3JRja/vz9vakukX6wdPk/23YA6fhYp\n39cBiyPiRyWbGu570dWzqNrvRVGtKyr9IatueprsJVhz0fmp4H0NImtV8hiwqOPeyFqIPAA8k/7j\ndyk5ZmJ6DkuAY0rSD0tfvKVkIzV0pG8H3A60k9X/7lf0fZfk7VbgBeAtsnrks/K693St9vQ5swqf\nxdlkL3efAB4n+wW7e70/C+AoYH36mViQPiMb8XvRxbM4tlq/F+74aWZmuaiXKjUzM6tyDjhmZpYL\nBxwzM8uFA46ZmeXCAcfMzHLhgGNmZrlwwLGqJemdNLT6Qkm3S9q+gDwMk3R93tetNElflfSTzezT\nJOnV9MwXS/qXzew/pXyYJbPuOOBYNXsjIg6JiGFknR3PL92YhtnobxeRjUtVlSRV+md4bmTj9n0c\nOEPSId3sey3Z8zHrEQccqxUPAR+WdLSkhyT9ClgkaStJV0man0bGHQPZ0B6S5paUkD6Z9r0hrT8h\n6ZvdXVDSdsCREfH7tN4i6cZ03mWSTpL0b+lcv+kIgMomsmpLo3vPLBnf69yUz8ck3dFRYpP0xZSn\nxyS1pbR3lUgk3Svp79Py6+m6jwF/J+kMSfPSvf60IwhJOkvS05LmAZ/ozcOOiDeAR9Iz/27K90JJ\nrSX7tAP7qWSiM7PuOOBY1Uu/yP+BbKgOyIZgHxcRHwW+BqyOiOFkAxSeq2yY9tOAmemv9QPJhvg4\nBNgrIoZFxIHAtM1c+hCy4ZJKDQI+TTZHyHRgdjrXm8BxygZS/AkwOiI+nq7x/XTsnRExPCIOJpu3\n5JyU/l1gREo/PqWVDwFSuv5e4OG0/1+Ak4FPpHtdD5yextJqIQs0R5FNvNXjYUWUDT9/JNlwSv87\n5XsYsL2kUSW7LgD+rqfntcaWR5WE2ZbaXtKCtDwXuB74JNlgg8+l9BHAMElfSOs7AR8Gfg9cnwLA\nXRHxuKRngf8haSrwa7LxtrrzQbKh3TsE8JuIeEfSImCriOgYrHAh2QyK+wNDyeYvgmw22hfSPsMk\nXQbsTDY/qWujAAACeklEQVQ748yU/jvgRkm3AzN68FzeIRsdGLLJsA4D/pCu9x5gJVnwbYtsVF8k\n/SLlbXM+JelRssA1OSKekjRa0kVkgW4g8CRwb9r/hXTfZpvlgGPV7M30V/sG6ZfqX8v2uyAiZpcf\nnF5ojwJukDQlIm6WdBDZTIXnk5UMzik/rkSw6ayJbwFExHpJ60rS15P9PAl4MiI6q8K6ATg+IhZK\nOhNoSuf6uqThwHHAI5IOA97m3TUQ7ylZXhvvHgTxxoiYWHbv5TPedjb7Y2ceioj/VXKe9wD/DhwW\nEc9LmlSWF9FA885Y37hKzWrd/cDYkvcn+yub0uEDwEsR8XOy+UAOTdVEW0fEDLJqrEM3c+7nyOaM\n742ngfdLOjLlZ4CkIWnbDsDKVOo6o+MASR+KiPkRMQl4iWwukmXAwcrsy8b5TMo9CHxB0vvTuQam\ne58HHJ3WBwBfJAUGSZ+XdHkP76cjuPxZ2SRfG86T7JnyarZZLuFYNevsL+fymRx/Tlal86iy4s+f\nyOb/aAIuSqWQNcBXyKa/nVbSsmsCgKTzACKilXd7HPhIN3na5D1LRKxL1XtTJe1M9jP2Q2AxWZCb\nRxZU5pEFIIArJQ0mKy08EBFPpHz9MR33FNkL/E2um6q8LgFmpftaB4yNiPmSWoD/AlaTvWvp8CHg\nVTa1ySyZEbFa0s/I3uWsTPkudQgwrpNzmW3C0xOYdUPSDcC1EVH+i7ZmSboZ+GbH+50+nGd/4N8i\n4vjN7myGA45ZtyR9DBgfEWcVnZdqI2kKMCMiflt0Xqw2OOCYmVku3GjAzMxy4YBjZma5cMAxM7Nc\nOOCYmVkuHHDMzCwXDjhmZpaL/w9zR+IShe9wTwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe92618c4d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(y,y-yhat,'b*')\n",
    "plt.xlabel('Press. (measured, Pa)')\n",
    "plt.ylabel('Residual')\n",
    "plt.title('RMSE=%.1f Pa' %(np.sqrt(sum((y-yhat)**2))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results from the polynomial fitting show a non-random trend in the residuals- this model is not ideal to describe the data...\n",
    "\n",
    "To use the Riedel equation, which is a linear combination of some (nonlinear) treaments of the Temperature data:\n",
    "\n",
    "$$ \\log(P)=b_0+\\frac{b_1}{T}+b_2\\log(T)+b_3T^2$$\n",
    "\n",
    "We calculate the $T$-containing terms, then stack them together and perform multiple linear regression:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "T=data['Temperature'].values[:,np.newaxis] # Temperature and pressure are changed to column vectors\n",
    "X=np.hstack((np.ones(T.shape),1/T,np.log10(T),T**2))\n",
    "y=np.log10(data['Pressure'].values[:,np.newaxis])\n",
    "model=LinearRegression(fit_intercept=False)\n",
    "model.fit(X,y)\n",
    "yhat=model.predict(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Spit out $b$ parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b_0 = 4.04e+01\n",
      "b_1 = -2.96e+03\n",
      "b_2 = -1.08e+01\n",
      "b_3 = 4.08e-06\n"
     ]
    }
   ],
   "source": [
    "for ix,b in enumerate(model.coef_.tolist()[0]):  # print out the coefficients\n",
    "    print('b_%d = %.2e' %(ix,b))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot residuals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x7fe925fea950>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAacAAAEZCAYAAAAzL+qdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuYXFWZ7/HvzxAUBYyRmYRcAB9JHJIzOAyS4JUWBULE\nBImHq9xnwpGJcY7oIUEcEo9jAI+IEcP0AIYAg9zFMEBI5NiCjiYid0Ik8RANgTQXCfdokPf8sVc1\nlUpVdVWnu/fu7t/neerpfVlr73dXUvXW3nvttRQRmJmZFclb8g7AzMyskpOTmZkVjpOTmZkVjpOT\nmZkVjpOTmZkVjpOTmZkVjpOTmZkVjpOTWSJpraRXJb0kaYOkKyXtXLb+cklvSJpSUe87afmJaX57\nSd+WtC5t63FJ36mxn9JrfpOxDpe0WNL6tO/dKtY/UrH9zZIWp3VjJf1Y0tOSnpO0RNLYBvY5VNIz\nku6uU6YlxfOSpBclrZJ0UjPHZgZOTmblAjgsInYC3g/8LXB2xfrHgBNKCyRtBxwJrEnrAWYDfw/s\nl7bVAtxbbT9lr5lNxvoGcBswreqBRIwv3z6wDrgurX4ncDMwFhgGrAB+3MA+zwNW8uZx1rI+7Xdn\n4EzgEkl7NbB9sw5OTmZVREQ7sBQYX7HqFuAjkoak+UnAA0B7WZkPADdHxIa0rd9HxJXdHN/TEfFv\nwD2dlZV0ALALcGOq++uIWBgRGyPideBC4H2S3lVnGx8iey8WAmoizh8DzwN7SfqUpPskvSDpD5LO\naXQ7NvA4OZltSQCSRpElnuUV6zeRnWUcneZPAK6oKPMr4EuSPi/pbyVV+zKv+gUv6SOSnq/z+lAX\njulE4IaIeK3G+o8BT0XE8zViGgR8D/inZnYq6S2SPgMMAR4CXgY+FxHvBD4FfF7S1Ga2aQOHk5PZ\nmwTcLOlF4A/A74BvVCl3BXCCpHeSfbHfXLF+HtklsOOAXwNPSDqhbH1pP+VJ51SAiPh5RLyrzuu/\nmjog6e1kl/4ur7F+FHAR8KU6m5kJ/Coi7mtwtyMkPQ88A3yNLCGtjoifRcQjABHxEHANcECD27QB\nZru8AzArkACmRsT/lfQxskt4HyC7J9NRJiJ+IemvyO5H3RIRm8pPjiLiDWABsEDSW4FTgR9IWh4R\nvy3fTy8c0xHAcxFxV+WKdAxLge9HxLXVKksaAXwB2LeJfT4ZEaOrbGsicC7Z5cHtgbfy5n0wsy34\nzMmsivRl/j2yM6BqriI726i8pFe5nT9FxAKy+y7jOtuvpI9WtLKrfH24yUM5sVqM6f7SUrJ7Y/Pq\n1J8A7AqslPQU2f2pCZKerHG5sp6ryc4yR0XEEODf8HeQ1eAzJ7PaLgT+p6SJEbGc7HJc6Qt5PnBX\nRGzVrFrSPwP3kZ1xbSa7vLdjWtZRrNoO0/Z2aiQ4SW/jzc/w2yS9LSI2la0fRdZScHpFvZ2BO4Cf\nR8RZnezmNmD3svmjgWOBKdH8eDs7As9HxJ8lTUjbuaPJbdgA4V8tZjVExLPAIrLm0JBdjou07vmI\n+GmNqq8A3waeIrvv8nlgWkSsLStzS8UZ0Y1dCPFV4MUU06q033LHA/8VEY9XLP8M2eXKk8v2/2JK\nZkg6TtLD6Tj/nFoGPh0RTwMvAH9O07XUSlqnA19P9/S+BlS9lGgGoDwHG5Q0iezX6SDg0ojY6hJK\nejjxULIP4kmlm7K16koaSvaffndgLXBkRGxMv9Ra02YHAf9aus4uaV+yG8ZvA26LiC/2yAGbmVlD\ncjtzSs1TLyJrrjsOOKbyQT1Jk4E9I2IM2aWJixuoOwtYFhFjgTvTPGRNWfeNiH2Ag4Hvp+2Qtntq\n2s+YlPjMzCwneV7WmwCsiYi1EbGZrFlp5TMPU8guq5Cu+Q+RNLyTuh110t/DU/3XUisqgB2AFyLi\nL5J2BXaKiFKLrCtKdczMLB95JqeRZF2qlDyRljVSZkSdusPS0/2QPbU/rFRI0gRJjwCP8OZzHSNT\n/ZL1VeIwM7NelGdyavRmVyPNVVVte6k1UZTNr4iI8WT9nn03PURpZmYFk2dT8vVA+YN6o9nyDKZa\nmVGpzOAqy9en6XZJwyNiQ7pkt1WroohYJel3wJ5pe6NqbKuDpPxajpiZ9WER0ewzcbmeOd1D1vhg\nD0nbA0cBiyvKLCb1AC1pf2BjumRXr+5isgcPSX9vTvX3SD1II2l3YAywOnXO+aKkiemhwuPZujsa\nACKi8K9zzjkn9xgcp+PsqzE6zu5/dVVuySmy3pBnkD2EtxK4NiIelXSapNNSmduA/ydpDVkz8NPr\n1U2bPhc4SNJjwIFpHuAjwP2S7gOuB6ZHxItp3enApcBqsoYWS3rw0M3MCi8imDXr/G1KMNsi1x4i\nIuJ24PaKZa0V8zMarZuW/xH4ZJXlV5F1OVNtW78hG7vHzMyAG2+8gwULnmK//ZYybdohvb5/9xDR\nz7S0tOQdQkMcZ/fqC3H2hRjBcba2XsX48Ydx1ll389JLFzB79l2MH38Yra1Vf9v3mFx7iOhLJIXf\nKzPr7yKCG25Ywhln3MW6dfMYPXo2F1xwANOmHULzff2CJKKPNYgwM7OCkYQkNm7cxLhxX2Ljxtc6\nlvUm90puZmZbWL16HQsXTuKIIw7mppuWsnr1us4rdTNf1muQL+uZmTXPl/XMzKzfcHIyM7PCcXIy\nM7PCcXIyM7PCcXIyM7PCcXIyM7PCcXIyM7PCcXIyM7PCcXIyM7PCcXIyM7PCcXIyM7PCyTU5SZok\naZWk1ZLOrFFmflr/gKR9OqsraaikZZIek7RU0pC0/CBJ90h6MP39eFmdtrSt+9Jrl548bjMzqy+3\n5CRpEHARMAkYBxwjaa+KMpOBPSNiDDAduLiBurOAZRExFrgzzQM8AxwWEXsDJwJXlu0qgGMjYp/0\nerbbD9jMzBqW55nTBGBNRKyNiM3ANcDUijJTgEUAEbEcGCJpeCd1O+qkv4en+vdHxIa0fCWwg6TB\nZfvq3cFKzMyspjyT00igfJCQJ9KyRsqMqFN3WES0p+l2YFiVfU8DfpMSW8midEnv7KaOwszMul2e\ngw02OjhSI2c0qra9iAhJWyyXNB44FziobPFxEfGkpB2BGyUdHxHll/0AmDNnTsd0S0sLLS0tDR2A\nmdlA0dbWRltb2zZvJ8/ktB4YXTY/muwMqF6ZUanM4CrL16fpdknDI2KDpF2Bp0uFJI0CbgKOj4jH\nS8sj4sn092VJV5NdNqybnMzMbGuVP9znzp3bpe3keVnvHmCMpD0kbQ8cBSyuKLMYOAFA0v7AxnTJ\nrl7dxWQNHkh/b071hwC3AmdGxC9LO5A0qNQ6L92D+jTwUHcfrJmZNS7XYdolHQpcCAwCLouIeZJO\nA4iI1lSm1CrvFeDkiLi3Vt20fChwHbAbsBY4MiI2pntJs4DVZSEcBLwG/IzsbGwQsAz4UuWY7B6m\n3cyseV0dpj3X5NSXODmZmTWvq8nJPUSYmVnhODmZmVnhODmZmVnhODmZmVnhODmZmVnhODmZmVnh\nODmZmVnhODmZmVnhODmZmVnhODmZmVnhODmZmVnhODmZmVnhODmZmVnhODmZmVnhODmZmVnhODmZ\nmVnh5JqcJE2StErSakln1igzP61/QNI+ndWVNFTSMkmPSVqahmdH0kGS7pH0YPr78bI6+0p6KG3r\nuz15zGZm1rnckpOkQUBpCPZxwDGS9qooMxnYMyLGANOBixuoOwtYFhFjgTvTPMAzwGERsTdwInBl\n2a4uBk5N+xkjaVJ3H6+ZmTUuzzOnCcCaiFgbEZuBa4CpFWWmAIsAImI5METS8E7qdtRJfw9P9e+P\niA1p+UpgB0mDJe0K7BQRK9K6K0p1zMwsH3kmp5HAurL5J9KyRsqMqFN3WES0p+l2YFiVfU8DfpMS\n28hUv2R9lTjMzKwXbZfjvqPBcmqwzFbbi4iQtMVySeOBc4GDGtx/hzlz5nRMt7S00NLS0uwmzMz6\ntba2Ntra2rZ5O3kmp/XA6LL50Wx5BlOtzKhUZnCV5evTdLuk4RGxIV2ye7pUSNIo4Cbg+Ih4vGwf\no2psawvlycnMzLZW+cN97ty5XdpOnpf17iFrfLCHpO2Bo4DFFWUWAycASNof2Jgu2dWru5iswQPp\n782p/hDgVuDMiPhlaQcR8RTwoqSJkgQcX6pjZmb5UESjV9d6YOfSocCFwCDgsoiYJ+k0gIhoTWVK\nrfJeAU6OiHtr1U3LhwLXAbsBa4EjI2KjpLPJWu6tLgvhoIh4VtK+wOXADsBtETGzSqyR53tlZtYX\nSSIiGrk9s2U9f+E2xsnJzKx5XU1O7iHCzMwKx8nJzMwKx8nJzMwKx8nJzMwKx8nJzMwKx8nJzMwK\nx8nJzMwKx8nJzMwKx8nJzMwKx8nJzMwKx8nJzMwKx8nJzMwKx8nJzMwKx8nJzMwKx8nJzMwKx8nJ\nzMwKJ9fkJGmSpFWSVks6s0aZ+Wn9A5L26ayupKGSlkl6TNLSNDx7aflPJb0k6XsV+2hL27ovvXbp\nqWM2M7PO5ZacJA0CSkOwjwOOkbRXRZnJwJ4RMQaYDlzcQN1ZwLKIGAvcmeYBNgFnA1+uEk4Ax0bE\nPun1bPcdqZmZNSvPM6cJwJqIWBsRm4FrgKkVZaYAiwAiYjkwRNLwTup21El/D0/1X42IXwB/qhFP\n08MIm5lZz8gzOY0E1pXNP5GWNVJmRJ26wyKiPU23A8Mqthk14lmULumd3Vj4ZmbWU7bLcd+1kkSl\nRs5oVG17ERGSGtnPcRHxpKQdgRslHR8RV1YWmjNnTsd0S0sLLS0tDWzazGzgaGtro62tbZu3k2dy\nWg+MLpsfTXYGVK/MqFRmcJXl69N0u6ThEbFB0q7A050FEhFPpr8vS7qa7LJh3eRkZmZbq/zhPnfu\n3C5tJ8/LevcAYyTtIWl74ChgcUWZxcAJAJL2BzamS3b16i4GTkzTJwI3V2xzizMxSYNKrfMkDQY+\nDTzUDcdnZmZdlNuZU0S8LmkGcAcwCLgsIh6VdFpa3xoRt0maLGkN8Apwcr26adPnAtdJOhVYCxxZ\n2qektcBOwPaSDgcOAv4ALEmJaRCwDLikZ4/ezMzqUUSjt34GNknh98rMrDmSiIimW0O7hwgzMysc\nJyczMyscJyczMyscJyczMyscJyczMyscJyczMyscJyczMyscJyczMyscJyczMyscJyczMyscJycz\nMyscJyczMyscJyczMyucmkNmSDqjTr2IiAt6IB4zM7O64zntRPWh1KsOiW5mZtZdch3PSdIk4EKy\nQf4ujYjzqpSZDxwKvAqcFBH31asraShwLbA7abDBiNiYlt8IfAC4PCK+ULaPfYHLgbcBt0XEF6vE\n4fGczMya1GPjOUnaQdIMSQskLZT0A0k/6FqYW2x3EHARMAkYBxwjaa+KMpOBPSNiDDAduLiBurOA\nZRExFrgzzQNsAs4GvlwlnIuBU9N+xqTEZ2bWtIhg1qzz8Y/ZbdNIg4grgWFkiaANGA283A37ngCs\niYi1EbEZuAaYWlFmCrAIICKWA0MkDe+kbked9PfwVP/ViPgF8KfyHUjaFdgpIlakRVeU6piZNevG\nG+9gwYKnuOmmpXmH0qc1kpz2jIivAS9HxCJgMjCxG/Y9ElhXNv9EWtZImRF16g6LiPY03U6WWMtV\n/pwZmeqXrK8Sh5lZXa2tVzF+/GGcddbdvPTSBcyefRfjxx9Ga+tVeYfWJ9VrEFHy5/T3BUl/C2wA\n/qob9t3oOW8j1yqrNtKIiJDkc2sz63HTpx/H0KHv5owz7gLEpk1v8M1vzmDatEPyDq1PaiQ5XZIa\nE5wNLAZ2BL7WDfteT3aJsGQ0W57BVCszKpUZXGX5+jTdLml4RGxIl+yebiCOUTW2tYU5c+Z0TLe0\ntNDS0tLJps1soJCEJDZu3MS4cV9i3bo3OpYNJG1tbbS1tW3zdnJrrSdpO+C3wCeAJ4EVwDER8WhZ\nmcnAjIiYLGl/4MKI2L9eXUnnA89FxHmSZgFDImJW2TZPAvataK23HJiZtnMrMD8illTE69Z6ZlbX\nvHmXMHbsbhxxxMHcdNNSVq9ex6xZ/5B3WLnqamu9TpOTpHPKZjsKR8TXm91ZlW0fypvNwS+LiHmS\nTkvbb01lSq3yXgFOjoh7a9VNy4cC1wG7UdaUPK1bS/b81vbARuCgiFhV1pR8B7Km5DOrxOrkZGbW\npJ5MTl/mzaS0A3AYsDIiTmk6yj7Mycm6IiKYPftbzJv3lQF3eccMejA5VdnRW4GlEXFAszvry5yc\nrCtuuGEJp5xyBwsXTvKNcRuQeuwh3CregZtam9XlZsW2LfwgbwOt9SQ9VDb7FuCvgW2+32TWn7lZ\nsW2L0oO8++23dMD+n2nkzOnTZa9DgBER8b0ejcqsj6tsVrxx42sDslmxNcdn3G+qN2TG0DT5YsWq\nndI1xD/2XFhmfd/q1etYuHDSFs2KzerxGfeb6l3Wu5eslZ7ImmU/n5a/C/g98J6eDc2sb5s9+x87\npgfil0tv6U8tIv0g75tqXtaLiD0i4j3AMuCwiHh3RLwb+FRaZmaWu/7W0WrpjPvhh7/NwoWHDtgz\n7kaec3o4Iv5bZ8v6OzclNyuW1tarmD//GjZvfj+rV3+DMWPOZvDgB5g582hOO+1zeYdnSVebkjfS\nt96Tks4GriK7xHcsNfqeMzPrLb4/07810lrvGLLm4z8CbkrTx/RkUGZmnXGLyP6t0zOniHiOrFNU\nM7NCcYvI/qvmPSdJ342IL0q6pcrqiIgpPRtasfiek5lZ83rintMV6e+3q6zzt7SZmfWYpjp+TQ/m\njoqIB3supGLymZOZWfN6rONXSW2Sdk6J6TfApZK+05UgzczMGtFIa70hEfEicARwRURMAD7Zs2GZ\nmdlA1khyGiRpV+BIsiHMoZvuOUmaJGmVpNWSzqxRZn5a/4CkfTqrK2mopGWSHpO0VNKQsnWzU/lV\nkg4uW96Wlt2XXrt0x/GZmVnXNJKcvg7cAfwuIlZIei+welt3LGkQUBqCfRxwjKS9KspMBvaMiDHA\ndODiBurOApZFxFjgzjSPpHHAUan8JGCB3nwgIoBjI2Kf9Hp2W4/PzMy6rtPkFBHXR8TeEfH5NP+7\niJjWDfueAKyJiLURsRm4BphaUWYKsCjtdzkwRNLwTup21El/D0/TU4EfRsTmiFgLrAEmlu3LT+6Z\nmRVEIw0i3ifpTkmPpPm9U3dG22okUP7E3BNsPcJurTIj6tQdFhHtabodGJamR6Ry5XVGlM0vSpf0\nuuPYzMxsGzRyWe8S4Czgz2n+Ibqn+6JG71s1ckajattLbb8b2c9xqSPbjwIflXR8g7GZmVkPaKTj\n17dHxPLS7ZmICEmbu2Hf64HRZfOj2fLMplqZUanM4CrLS53RtksaHhEbUkOOp+tsaz1ARDyZ/r4s\n6Wqyy4ZXVgY8Z86cjumWlhZaWlo6O0ZL+tOYO2ZWW1tbG21tbdu+oYio+wJuB/YE7kvznwVu76xe\nA9vdDvgdsAewPXA/sFdFmcnAbWl6f+BXndUFzgfOTNOzgHPT9LhUbnuygRJ/R3bGNQjYJZUZDNwA\nTK8Sb1jXXX/97bHTTv8cN9ywJO9QzKwXpe/OpnNEI+M5vRf4d+CDwEbgcbLLYGubS4NVt30ocGFK\nEJdFxDxJp6VM0JrKlFrlvQKcHBH31qqblg8FriMbvXctcGREbEzrzgJOAV4HvhgRd0h6B/AzssQ0\niGwgxS9FxRvjHiK6xmPumA1sXe0houHuiyTtSHam8TLZF/61ze6sL3Ny6pqI4IYblnDGGXexbt08\nRo+ezQUXHMC0aYf48p7ZANDt3RdJ2lHSGZIWSDodeJWsZ4hHgOO6HqoNJB5zx8y6orNeyV8Efgkc\nDJwEbCJ7WPX+ng/N+guPuWNmzao3ntODEbF3mh4EPAXsHhGv9WJ8heHLemZmzeuJXsn/UpqIiL8A\n6wdqYjIzs95V78zpL2T3mUp2AErJKSJi5x6OrVB85mRm1rxuHwk3IgZtW0hmZmZd00j3RWZmZr3K\nycnMzArHycnMzArHycnMzArHycnMzArHycnMBoyIYNas8/FjIcXn5GRmA8aNN97BggVPcdNNS/MO\nxTrh5GRm/V5r61WMH38YZ511Ny+9dAGzZ9/F+PGH0dp6Vd6hWQ2NjIRrZtanTZ9+HEOHvpszzrgL\nEJs2vcE3vzmDadMOyTs0q8FnTmbW73nolr4n1+QkaZKkVZJWSzqzRpn5af0DkvbprK6koZKWSXpM\n0lJJQ8rWzU7lV0k6uGz5vpIeSuu+21PHa2b5KQ3d8vDD32bhwkPrDt3ihhP5a3gk3G7fcTYMx2/J\nBjBcD/waOCYiHi0rMxmYERGTJU0EvhsR+9erK+l84NmIOD8lrXdFxCxJ44Crgf2AkcBPgDEREZJW\npP2skHQbMD8illTE645fzQaIG25Ywimn3MHChZN86W8b9cSQGT1tArAmItZGxGbgGmBqRZkpwCKA\niFgODJE0vJO6HXXS38PT9FTghxGxOSLWAmuAiZJ2BXaKiBWp3BVldcxsAHHDieLIs0HESKD8vPoJ\nYGIDZUYCI+rUHRYR7Wm6HRiWpkcAv6qyrc1pumR9Wm5mA4wbThRHnsmp0WtkjZwOqtr20iW7brsW\nN2fOnI7plpYWWlpaumvTZlYAlQ0n1q17ww0nmtTW1kZbW9s2byfP5LQeGF02P5otz2CqlRmVygyu\nsnx9mm6XNDwiNqRLdk93sq31abratrZQnpzMrH8qNZw44oiDuemmpXUbTtjWKn+4z507t0vbybNB\nxHZkjRo+ATwJrKB+g4j9gQtTg4iadVODiOci4jxJs4AhFQ0iJvBmg4g909nVcmBm2s6tuEGEmVm3\n6PaRcHtaRLwuaQZwBzAIuCwll9PS+taIuE3SZElrgFeAk+vVTZs+F7hO0qnAWuDIVGelpOuAlcDr\nwOll2eZ04HKyoehvq0xMZmbWu3I7c+prfOZkZta8vtiU3MzMrConJzMrrKL21FDUuPoTJyczK6yi\nDnFR1Lj6EycnG3D8q7f4itpTQ1Hj6o+cnArKX6A9x796i2/69OOYM+ef2LTpDUo9NcydO4Pp049z\nXAOEk1NB+Qu0+/lXb99R1CEuihpXf+TkVDD+Au05/tXbtzQzxEVvKmpc/Y2fc2pQbz3nFBHccMMS\nzjjjLtatm8fo0bO54IIDmDbtEP866waloRBGjxbr1r3BwoWHulNPsx7k55z6CV826Fn+1WsDUV+8\nh+0zpwb1Zg8R8+Zdwtixu23R8eSsWf/QK/s2s54XEcye/S3mzftKr/zwzHPwxK6eOTk5NcjdF5lt\nqbe/YPuT3koWra1XMX/+NWze/H5Wr/4GY8aczeDBDzBz5tGcdtrnemy/5XxZz2yAyftSjVuUNq+3\nGzz15UZATk5mfVReycEtSruut5NFX76H7eRk1sfknRz68q/xvOWRLPpqI6A8R8I1sy6YPv04hg59\nN2eccRel5PDNb87otRvdXR3K3PeoMr090u7s2f/YMd2XHptwcjLrY7qaHLpTV75gS5ch99tvaZ/6\nkuxufTVZ9LZcWutJGgpcC+xOGq02IjZWKTcJuJBstNtLI+K8zupLmg2cAvwFmBkRS9PyfclGu30b\n2Wi3X0zLTwK+BTyRdvu9iPhBlVjcWs8Koy89blCEFmOWnz7VlFzS+cCzEXG+pDOBd0XErIoyg4Df\nAp8E1gO/Bo5JQ7lXrS9pHHA1sB8wEvgJMCYiQtIKYEZErJB0GzA/IpZIOhHYNyJmdhKzk5NZF7jX\nk4GtrzUlnwIsStOLgMOrlJkArImItRGxGbgGmNpJ/anADyNic0SsBdYAEyXtCuwUEStSuSvK6ii9\nzKwH9OUWY5afvJLTsIhoT9PtwLAqZUYC5Reyn0jL6tUfwZuX58rrVC5fX7atAKZJelDS9ZJGdeF4\nzKyOvtpizPLTYw0iJC0DhldZ9dXymXTJrdr1ssplqrKsXv1G3QJcHRGbJU0nOxP7RLWCc+bM6Zhu\naWmhpaVlG3ZrNnC4EcDA0dbWRltb2zZvJ697TquAlojYkC65/TQi/qaizP7AnIiYlOZnA29ExHm1\n6kuaBRAR56Y6S4BzgN+nMnul5ccAB0TE/6jY5yDguYgYUiVm33Mys9z1tSb5fe2e02LgxDR9InBz\nlTL3AGMk7SFpe+CoVK9e/cXA0ZK2l/QeYAywIiI2AC9KmqjsX/P4Uh1J5Wd3U4CV3XGAZmY9YaB0\nG5VXcjoXOEjSY8CBaR5JIyTdChARrwMzgDvIEsa1EfFovfoRsRK4LpW/HTi97HTndOBSYDVZQ4sl\naflMSQ9Luj/t76QeO2qzAsq7jz5rTN49g/Q290reIF/Ws/4qz+EUrHF9tUl+X7usZ2Y5G2i/xPu6\ngdYk390XmQ1QeffRZ83r7X758uTkZDZAFaGPPmvOQGqS7+RkNoANpF/i1re4QUSD3CDCzKx5bhBh\nZmb9hpOTmVk/0N+eV3NyMjPrB/pbzxFOTmZmfVh/fV7NrfXMzPqw/vq8ms+czMz6sP7ac4TPnMzM\n+rj++Lyan3NqkJ9zMjNrnp9zMjOzfsPJyczMCsfJyczMCieX5CRpqKRlkh6TtFTSkBrlJklaJWm1\npDMbqS9pdiq/StLBZcv/VdIfJL1UsY+3Sro21fmVpN174pjNzKxxeZ05zQKWRcRY4M40vwVJg4CL\ngEnAOOAYSXvVqy9pHHBUKj8JWKA321P+GJhQJZZTgeciYgzwHeC8bjlCMzPrsryS0xRgUZpeBBxe\npcwEYE1ErI2IzcA1wNRO6k8FfhgRmyNiLbAGmAgQESsiYkMnsdwIfKKrB2VmZt0jr+Q0LCLa03Q7\nMKxKmZFAeWP9J9KyevVHpHLV6tTSsZ+IeB14QdLQRg7CzMx6Ro89hCtpGTC8yqqvls9EREiq9gBR\n5TJVWVavfq3tdNmcOXM6pltaWmhpaemuTZuZ9QttbW20tbVt83Z6LDlFxEG11klqlzQ8IjZI2hV4\nukqx9cDosvlRaRlArfr16tSyHtgNeFLSdsA7I+KP1QqWJyczM9ta5Q/3uXPndmk7eV3WWwycmKZP\nBG6uUua7m5VHAAALQUlEQVQeYIykPSRtT9bQYXEn9RcDR0vaXtJ7gDHAiiZi+SxZAwszs36vyGNA\n5ZWczgUOkvQYcGCaR9IISbdCx/2fGcAdwErg2oh4tF79iFgJXJfK3w6cXupzSNL5ktYBO0haJ+lf\n0rYuA94taTXwz1RpOWhm1h8VeQwo963XIPetZ2b9RWvrVcyffw2bN7+f1au/wZgxZzN48APMnHk0\np532uW7dV1f71nOv5GZmA0xfGAPK3ReZmQ0wfWEMKJ85mZkNQEUfA8r3nBrke05mZs3zeE5mZtZv\nODmZmVnhODmZmVnhODmZmVnDeqtXCScnMzNrWG/1KuHkZGZmnWptvYrx4w/jrLPu5qWXLmD27LsY\nP/4wWluv6pH9+TknMzPrVG/3KuEzJzMz61Rv9yrhMyczM2tIb/Yq4R4iGuQeIszMmuceIszMrN/I\nJTlJGippmaTHJC2VNKRGuUmSVklaLenMRupLmp3Kr5J0cNnyf5X0B0kvVezjJEnPSLovvU7piWM2\nM7PG5XXmNAtYFhFjyYZF32r0WUmDgIuAScA44BhJe9WrL2kc2XDu41K9BXrzbt2PgQlVYgnghxGx\nT3r9oJuOMRdtbW15h9AQx9m9+kKcfSFGcJxFkVdymgIsStOLgMOrlJkArImItRGxGbgGmNpJ/alk\niWZzRKwF1gATASJiRURsqLIfpVe/0Ff+wzrO7tUX4uwLMYLjLIq8ktOwiGhP0+3AsCplRgLlTUGe\nSMvq1R+RylWrU0sA0yQ9KOl6SaMaPAYzM+shPdaUXNIyYHiVVV8tn4mIkFStGVzlMlVZVq9+re1U\nugW4OiI2S5pOdib2iU7qmJlZT4qIXn8Bq4DhaXpXYFWVMvsDS8rmZwNn1qtPdu9pVlmdJcDEiu2+\nVCeuQcDGGuvCL7/88suv5l9dyRN5PYS7GDgROC/9vblKmXuAMZL2AJ4ka+hwTCf1FwNXS7qA7HLe\nGGBFvUAkDS+7FzUFWFmtXFfa6ZuZWdfkdc/pXOAgSY8BB6Z5JI2QdCtARLwOzADuIEsY10bEo/Xq\nR8RK4LpU/nbg9NKTs5LOl7QO2EHSOkn/krY1U9LDku5P+zupZw/dzMw64x4izMyscNxDRBlJb5O0\nXNL9klZKmlelzC6SlqQyD0s6KYdQS7EMSg8O31Jj/fz0QPIDkvbp7fjK4qgZp6TjUnwPSvqFpL2L\nFmNZmf0kvS7piN6MrSKGzv7NW9L6hyW19XJ45XHU+zcvxGdI0tr0/+4+SVUv/xfhM9RZnAX6DHX6\nfqZyDX2O3PFrmYjYJOnjEfGqpO2An0v6SET8vKzYDOC+iJgtaRfgt5KuSpche9sXyS5h7lS5QtJk\nYM+IGCNpInAxWSOTPNSME/h/wMci4gVJk4B/J58468VYeij8PLJGNnnef6z3bz4E+D5wSEQ8kf5/\n5qXe+1mUz1AALRHxx2orC/QZqhsnxfkMdRZnU58jnzlViIhX0+T2ZK33Kt/op4Cd0/TOwHN5JKb0\nPNZk4FKq/yN3PKgcEcuBIZKqPU/WozqLMyJ+GREvpNnlQK8/Z9bAewnwBeAG4JneiqtSA3EeC9wY\nEU8ARMSzvRhehwbiLMRnKKn3BVmIz1BSM84ifIbKdPbDreHPkZNTBUlvSY0j2oGfpkYW5S4Bxkt6\nEniA7BdiHr4DfAV4o8b6ag8x5/GftrM4y50K3Naz4VRVN0ZJI8l6H7k4LcrrRm1n7+UYYKikn0q6\nR9LxvRfaFjqLsyifoQB+kt6rf6yyviifoc7iLJfXZwg6ibPZz5GTU4WIeCMi/o7sP+HHJLVUFDkL\nuD8iRgB/B3xfUtVLQT1F0mHA0xFxH/V/qVSu69Uv1SbiRNLHgVOAM+uV624Nxngh2fNzQU7dXTUY\n52Dg78nOWg4BviZpTC+FCDQcZ+6foeTDEbEPcCjwT5I+WqVMrp+hpJE4c/sMlekszqY+R05ONaTT\n5FuBD1Ss+hBwfSrzO+Bx4H29Gx0fAqZIehz4IXCgpCsqyqwHRpfNj0rLelMjcZJu4F4CTImI5wsY\n477ANanMNLIOhacUMM51wNKIeC0ingPuAt5fwDiL8BkiIp5Kf58BfsTWHUMX4TPUSJx5f4aAhuJs\n7nPUkz1B9LUXsAswJE3vQPbh/kRFmQuAc9L0MLJT/aE5xnwAcEuV5ZOB29L0/sCvcn5va8W5G1kH\nvfsX4N+/aowVZRYCRxQxTuBvgJ+Q3St9O/AQMK6Aceb+GUrvz05p+h3AL4CDK8rk/hlqMM7cP0ON\nxFlRvtPPkVvrbWlXYJGkt5CdVV4ZEXdKOg0gIlqBbwILJT2QyvyvqNM6pZeUHjTuiDMibpM0WdIa\n4BXg5DwDTLaKE/gX4F3AxcpGN9kcEdWGNukt1WIsomr/5qskLQEeJLvfc0lsfc+0t1V7P4vwGRoG\n/Cj9n9sO+I+IWFrAz1CncVKMz1AjcTbFD+GamVnh+J6TmZkVjpOTmZkVjpOTmZkVjpOTmZkVjpOT\nmZkVjpOTmZkVjpOTWRlJL29D3TZJq9JQED+XNLaT8j/JqduebtXseyZpb0mX9VQ81j84OZltaVse\n/Avg2Mj6ZlwEfKtWQUkHAr+NiJe2YX89Jg1t0Kim3rOIeBB4r6S/bi4qG0icnMyqUOZbkh5KA6gd\nmZa/RdICSY9KWirpVknTqmzibmDPOrs4Fvhx2uYe6YxroaTfSvoPSQcrGzjuMUn7pXLvkPQDZQNi\n3lvqlyzVv0vSb9Lrg2n5rmn5fek4PpyWd5zpSPqspIVp+nJJ/ybpV8B5kt4r6fbUy/Rdkt6Xyr1H\n0i/T+/KNLr7FtwP/vYt1bQBwcjKr7giyTlP3Bj4JfEvS8LR894jYCzge+CBbnjmUelr+NFk3QrV8\nGLinbP69wP8h6x/vfcBREfFh4MtkvXgDfBW4MyImAgemmN5ONrzLQRGxL3A0MD+VPxZYEllP0e8n\nG56Cingrz3pGAB+MiC+TDVr3hYj4ANkQGAtSme8C34+IvYEn6xxjPSuAj3Wxrg0A7lvPrLqPAFdH\n1r/X05J+BuxHllSuA4iIdkk/Lasj4D8kvUbW0/YX6mx/REV/co9HxCMAkh4h68AV4GFgjzR9MPBp\nSV9O828l6zV7A3CRpPcDfyEb1wmyBPADSYOBmyOilJxqCeD6iAhJO5Il3utTf2mQDcAJWa/in0nT\nV5GNbNqsp8qOy2wrTk5m1ZXGnKmm1vLSPad7u7C/P5VNvwH8uWy6/HN6RESs3iIYaQ7wVEQcn+4V\nbQKIiLvTmDqHAZdLuiAirmTLs6UdKuIojQT9FmBjOuvaJpIOB85Js6em90ds2/096+d8Wc+suruB\no9I9pr8iuwS1nGwogGnpntQwoKWiXqMDET4p6d1NxnQHMLNjR1IpcexMdvYEcALZkBlI2g14JiIu\nBS4DSuXbJf1N6n3/M1RJEhHxIvC4pM+mbSmNGQTZe3B0mj6uvJ6kVVW2dXNE7JNepcS9K/D7ho/c\nBhwnJ7MtBUBE/IjsntEDwJ3AVyLiaeBGsvGHVgJXAvcCL1TWB5A0QtKtNfbzc7YcyLIyQVS7L/S/\ngcGpIcLDwNy0fAFwoqT7ye5XlRo8fBy4X9K9ZI0PvpuWzwL+kyzJVN4zKt/vccCpabsPA6WB4b5I\nNtLpg2T3qErDYuxS41irmUA2XppZVR4yw6xJkt4REa+kM5/lwIdS4mpmGy1kjR4+3xMx5kHSp4D3\nRMRFDZRtA45s9n2zgcP3nMya95+ShpA1EPh6V75gI6JN0tmSdirqs07NiohaZ4lbSJcH1zgxWT0+\nczIzs8LxPSczMyscJyczMyscJyczMyscJyczMyscJyczMyscJyczMyuc/w/yltMpvZs7TgAAAABJ\nRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe926030290>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(y,y-yhat,'b*')\n",
    "plt.xlabel('logP. (measured, -)')\n",
    "plt.ylabel('Residual')\n",
    "plt.title('RMSE=%.1f Pa' %(np.sqrt(sum((10**y-10**yhat)**2))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These residuals look much better- no clear trends and much lower residual error.  The Riedel equation is a better fit than the polynomial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Nonlinear Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Regress P vs T data to the Antoine equation:\n",
    "    \n",
    " $$ \\log\\left(P\\right)=A+\\frac{B}{T+C}$$\n",
    " \n",
    " We use `scipy`'s `curve_fit` with a custom function coding in the Antoine equation.  *Note: Doing curve fitting with log-transformed observations violates the assumtion of identically-distributed errors*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define Antoine function, regress and print out best fit values of $A$, $B$, and $C$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A: 9.02\n",
      "B: -1207.37\n",
      "C: -52.77\n"
     ]
    }
   ],
   "source": [
    "def NLRfun(X,A,B,C):\n",
    "    T=X\n",
    "    logP=A+B/(T+C)\n",
    "    yhat=logP\n",
    "    return yhat\n",
    "\n",
    "X=data['Temperature'].values\n",
    "y=np.log10(data['Pressure'].values)\n",
    "Beta0=[10, -2000, 0]\n",
    "Beta_names=['A: ','B: ','C: ']\n",
    "\n",
    "beta,covar = curve_fit(NLRfun,X,y,p0=Beta0)\n",
    "for b in zip(Beta_names, beta):\n",
    "    print(b[0]+'{:.2f}'.format(b[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make plot of residuals, and also calculate the RMSE of the predicted pressures:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x7fe925f3cfd0>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAEZCAYAAADR8/HkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuYXFWZ7/HvjwQGFDAEnISQIBwISHIE8RJQUVpFCAED\nEoeLKChoOI4R54Az6TAowaNE8IiAGImIAUFFCIhhgCTIYxtwBETCJUIgYYgmAcJFAghEYfqdP/bq\nUCmqqququ2pXd/0+z1NP78tae79V3bveXnuvvbYiAjMzs2bbJO8AzMysPTkBmZlZLpyAzMwsF05A\nZmaWCycgMzPLhROQmZnlwgnIzMxy4QRkbUfSSkkvSXpB0hOSLpe0dcH6SyV1S5pcVO87afnxaX4z\nSd+WtCpt61FJ3ymzn57XBTXGOlLSfElr0r53LFq/taQrJD2VXldI2qpg/UclLU37/q2kPSrs61JJ\nfyuI9XlJKlO2I8XTU26ZpE/X8t7MnICsHQVwaERsBewFvA04vWj9w8BxPQskDQWOBFak9QAzgHcA\n707b6gDuLrWfgtfJNcbaDdwITCmzfiawHbAzsAswIi1D0ljgCmAq8CbgemC+pCFlthXA2QWxbh2V\n71Rf01MOmA5cXCnBmRVzArK2FhFrgUXA+KJV1wP7SRqW5icC9wJrC8q8C7guIp5I2/pTRFzez/E9\nGREXAXeVKTI+xfDXiHgeuI7X3stBwK0R8Z8R0Q2cDewA7F9hlyVbPFXE+UvgWWAPSYdIWiLpOUl/\nlnRGPdu0wc8JyNqVACSNJksudxStXw/8Ejg6zR8H/LiozO3AKZI+L+ltZU5XlTuFtZ+kZyu83lvl\n+1gITJE0TNI2ZC2lG9O6KNr/Jmm+ONkW+mdJz0i6S9IR1QQgaRNJHwOGAfcDfwU+GRFvAg4BPi/p\nsCrfj7URJyBrRwKuk/Q88GfgEeDrJcr9GDhO0puAD5C1LgrNImtVHAv8Hlgt6biC9T37KUwsJwJE\nxG0RsU2F139W+V6+l34+AzwNvAJ8Py37FbC/pP0lbQacBmwGvKHMti4AdgXeDHwFuLSXRDhK0rPA\nU6n8JyNieUT8JiL+mN7n/cCVVG51WZtyArJ2FMBh6dpFB/AhstNpG5WJiN+SfRmfDlwfEeuLCnRH\nxOyI2I/sGss3gB9J2r1oP4WJ5ZJ+fi8/AR4CtgS2Bv6L7LoPEfEQcDxwIfAYsC3wALC61IYiYklE\nPJve101p25VaQY+l97RtRLwjIq4CkLSPpF9LelLSOuCktG+zjTgBWVuLiMXAd8laMqVcAZzC60+/\nFW/nbxExm+w6yLje9ivp/UW944pf76vyLUwE5kTEyxHxIjAHmFQQ1zUR8baI2I6sc8JOZK21Rvop\nWWtxdEQMAy7C3zVWgv8ozOA8YIKkfdK8eO3ayQXAARFxa3ElSf+STm9tIWlo6p69JbCksFipHUbE\nrUW944pfvy3Yz+bA5ml28zTf4z7gc5I2l7QFWY+3ewvqvlPSEElvBn4A/DIiHi4Vk6SPS9oyXdM5\nkOzU4vzSH1lFWwLPRsTfJU0APsFrPQfNNnACsrYXEU8Dl5F1JYbsyzLSumcj4tdlqr4IfBt4nOw6\nyOeBKRGxsqDM9UUtm2vqCPEl4PkU07K03x6fBnYD1pCdWtuJ7LRbj/PIWmXLyK4Tfa5nhaRjJS0t\nKHty2sazZC3Cz6YWYjnlkso/A19L19i+Avy84ruztqU8H0gnaSLZATIE+GFEvO40SLpx72Cyg/DT\nEbGkUl1Jw8n+4N8CrASOjIh1Bdvbkew8+BkR8e3GvTszM6sktxZQuhnuQrJz2OOAY4pvYpM0Cdg1\nIsaSnVr4fhV1O4GbI2I34JY0X+hc4IaGvCkzM6tanqfgJgArImJlRLxC1lWz+F6ByWSnRoiIO4Bh\nkkb2UndDnfTz8J6NSTqcrJfQA415S2ZmVq08E9AOwKqC+dVpWTVlRlWoOyLd3Q7ZXesjACRtCfwb\naZgSMzPLV54JqNqLT9UMDaJS20vjWPUsnwl8JyJeqnKbZmbWQENz3PcaYEzB/Bhef4NccZnRqcym\nJZavSdNrJY2MiCckbQ88mZZPIBuy5ByyIUO6Jb2c7t3YQJK7i5qZ1Sgiav7HPs8W0F3AWEk7pWFC\njuL19xzMJ41ILGlfYF06vVap7nxe64Z6PGn4lIj4QETsHBE7k/We+0Zx8ukRES3/OuOMM3KPYTDE\n6DgdZ6u/BkKc9cqtBRQRr0qaRjaY4hDgkoh4UNJJaf2ciLhR0iRJK8juffhMpbpp098Erkpjbq0k\nG0LfzMxaTJ6n4IhsvKmbipbNKZqfVm3dtPwvwAG97PfMmoM1M7N+5ZEQBqiOjo68Q+jVQIgRHGd/\nc5z9a6DEWY9cR0JoRZLCn4mZWfUkEQOsE4KZmbUxJyAzM8uFE5CZmeXCCcjMzHLhBGRmZrlwAjIz\ns1w4AZmZtamIoLPznD4Np9MXTkBmZm3qmmsWMnv241x77aJc9u8EZGbWZubMuYLx4w/ltNNu5YUX\nzmXGjMWMH38oc+Zc0dQ4ch0LzszMmm/q1GMZPnxbTj11MSDWr+/mrLOmMWXKQU2Nwy0gM7M2IwlJ\nrFu3nnHjTmHdupc3LGsmt4DMzNrQ8uWrmDt3IkcccSDXXruI5ctXNT0GD0ZaxIORmpnVxoORmpnZ\ngOIEZGZmuXACMjOzXDgBmZlZLpyAzMwsF05AZmaWCycgMzPLhROQmZnlwgnIzMxy4QRkZma5cAIy\nM7NcOAGZmVkunIDMzCwXuSYgSRMlLZO0XNL0MmUuSOvvlbR3b3UlDZd0s6SHJS2SNCwtnyBpSXrd\nJ+moxr9DMzMrJ7fHMUgaAjwEHACsAX4PHBMRDxaUmQRMi4hJkvYBzo+IfSvVlXQO8HREnJMS0zYR\n0SlpC+BvEdEtaSSwFBgREf9dFJcfx2BmVoOB+DiGCcCKiFgZEa8AVwKHFZWZDFwGEBF3AMNS8qhU\nd0Od9PPwVP/liOhOy7cAnitOPmZm1jx5JqAdgMJH8K1Oy6opM6pC3RERsTZNrwVG9BRKp+H+CPwR\nOKWvb8DMzOqXZwKq9jxXNc06ldpeOpcWBfN3RsR44B3A+ZLeVGUMZmbWz4bmuO81wJiC+TFkLZlK\nZUanMpuWWL4mTa+VNDIinpC0PfBk8Y4jYpmkR4BdgT8Ur585c+aG6Y6ODjo6Oqp7R2ZmbaCrq4uu\nrq4+byfPTghDyToSfBh4DLiTyp0Q9gXOS50QytZNnRCeiYizJXUCw1InhJ2A1RHxqqS3ALcC/zsi\nni+Ky50QzMxqUG8nhNxaQCkRTAMWAkOAS1ICOSmtnxMRN0qaJGkF8CLwmUp106a/CVwl6URgJXBk\nWr4f0CnpFeAVYGpx8jEzs+bJrQXUqtwCMjOrzUDshm1mZm3MCcjMzHLhBGRmZrlwAjIzs1w4AZmZ\nWS6cgMzMLBdOQGZmlgsnIDMzy4UTkJmZ5cIJyMzMcuEEZGZmuXACMjOzXDgBmZlZLpyAzMwsF05A\nZmaWCycgMzPLhROQmZnlwgnIzMxy4QRkZma5cAIyM7NcOAGZmVkunIDMzCwXTkBmZpYLJyAzM8uF\nE5CZmeXCCcjMzHLhBGRmZrlwAjIzs1w4AZmZWS5yTUCSJkpaJmm5pOllylyQ1t8rae/e6koaLulm\nSQ9LWiRpWFr+EUl3Sbov/fxg49+hmZmVo4ioXCD7An8PsBMQwErgdxHxXJ92LA0BHgIOANYAvweO\niYgHC8pMAqZFxCRJ+wDnR8S+lepKOgd4OiLOSYlpm4jolPR24ImIeELSeGBhRIwuEVf09pmYmdlr\nJBERqrVe2RaQpPdLmg8sBo4GdiRLQscAt0qaL2m/OuMFmACsiIiVEfEKcCVwWFGZycBlABFxBzBM\n0she6m6ok34enurfExFPpOUPAFtI2rQP8ZuZWR8MrbDuY8CpEbG81EpJuwH/B7itzn3vAKwqmF8N\n7FNFmR2AURXqjoiItWl6LTCixL6nAH9IycvMzHJQNgFFxCmVKkbEw0DFMr2o9jxXNc06ldpeRISk\njZan02/fBD5SbmMzZ87cMN3R0UFHR0eVoZqZDX5dXV10dXX1eTuVWkAbSDoUGA9sTvqij4iv9XHf\na4AxBfNjyFoylcqMTmU2LbF8TZpeK2lkutazPfBkwfsYDVwLfCoiHi0XWGECMjOzjRX/Y37mmWfW\ntZ1ee8FJmgMcCXwxLToSeEtde9vYXcBYSTtJ2gw4CphfVGY+cFyKY19gXTq9VqnufOD4NH08cF2q\nPwy4AZgeEb/rh/jNzKwPqukFd39EvE3SfRGxp6QtgQUR0ZcOCD3bPhg4DxgCXBIRsySdBBARc1KZ\nC4GJwIvAZyLi7nJ10/LhwFVknSZWAkdGxDpJpwOdQOE1rY9ExNNFMbkXnJlZDertBVdNArozIiZI\nup3s4v0zwNKI2LW+UFubE5CZWW3qTUDVXAO6XtI2wLeAP6RlF9e6IzMzs0IVW0Bp5IFdyVo8D0ra\nHNg8ItY1K8BmcwvIzKw2jbgR9avAz4EjgBslTY2I9YM5+ZiZWfOUbQFJegB4V0S8JGlbsqFr3tXU\n6HLgFpCZWW36vQUE/C0iXgKIiGd6KWtmZlaTSi2g58jGgevxfuDWNB0RMbnBseXCLSAzs9r0ezds\nSR0V6kVE/KbWnQ0ETkBmZrVp2H1A7cYJyMysNo3oBXeDpH+S9IYS694g6ShJN9a6QzMzM6h8Cu4f\ngWnAx4H/Bh4nG3V6JNkNrD8HvhcRTzUn1OZwC8jMrDYNPQWXHgLXMwDpnwoe7DboOAGZmdXG14D6\niROQmVltGnEfUM+GXyjxWi3pF5L+V33hmplZu6tmMNLzyR5//bM0fzSwC7AE+BHQ0ZDIzMxsUKvm\ncQz3RcSeRcvuiYi3S7o3IvZqaIRN5lNwZma1adgpOOCl1OV6k/Q6Elif1vmb2szM6lJNC2gXstNw\n+6ZFtwP/AqwB3hkRtzU0wiZzC8jMrDYNawFFxCMRcWhEbJdeh0bEioh4ebAlH7N6RQSdnefgf17M\nqldNL7gxqcfbU+l1jaTRzQjObKC45pqFzJ79ONdeuyjvUMwGjGquAc0F5gOj0uv6tMys7c2ZcwXj\nxx/KaafdygsvnMuMGYsZP/5Q5sy5Iu/QzFpeNd2w3xwRhQnnUkn/t1EBmQ0kU6cey/Dh23LqqYsB\nsX59N2edNY0pUw7KOzSzlldNC+gZSZ+SNETSUEmfBJ5udGBmA4EkJLFu3XrGjTuFdete3rDMzCqr\nJgGdABwJPEE2IOk/AZ9pZFBmA8ny5auYO3ciS5d+m7lzD2b58lV5h2Q2IHgsuCLuhm1mVpt6u2GX\nvQYk6bsV6kVEnFzrzszMzHpU6oTwB0qPdKAyy83MzKrmU3BFfArOzKw2jRwLzszMrN/lmoAkTZS0\nTNJySdPLlLkgrb9X0t691ZU0XNLNkh6WtEjSsILlv07PM6p0fcvMzJogtwQkaQhwITARGAccI2mP\nojKTgF0jYiwwFfh+FXU7gZsjYjfgljQP2QjepwNfbuT7MjOz6tScgCR9IT2eoZpRFCqZAKyIiJUR\n8QpwJXBYUZnJwGUAEXEHMEzSyF7qbqiTfh6e6r8UEb8F/tbHuM3MrB/U0wIS8H7gF33c9w5kT1rt\nsTotq6bMqAp1R0TE2jS9FhhRtE33MDAzawE1t2Ii4sJ+2ne1iaCanhUlu4ZHREiqOeHMnDlzw3RH\nRwcdHR21bsLMbNDq6uqiq6urz9updCPqqRXqRUSc28d9rwHGFMyPIWvJVCozOpXZtMTyNWl6raSR\nEfGEpO2BJ2sNrDABmZnZxor/MT/zzDPr2k6lU3BbAVuWeG2VXn11FzBW0k6SNgOOInvsQ6H5wHEA\nkvYF1qXTa5XqzgeOT9PHA9cVbdOjRJqZtYBcb0SVdDBwHjAEuCQiZkk6CSAi5qQyPb3dXgQ+ExF3\nl6ublg8HrgJ2BFYCR0bEurRuJVny3Ax4FjgwIpYVxeQbUc3MalDvjai9JiBJWwAnknV33oJ0rSUi\nTqgjzpbnBGRmVptGjoRwOVlPsolAF9m1l7/WuiMzM7NC1bSA7omIt0u6LyL2lLQpcFtE7NOcEJvL\nLSAzs9o0sgX09/TzOUlvA4YBb651R2ZmZoWquQ/o4nRh/3SyHmZbAl9paFRmZjbo+XEMRXwKzsys\nNv3+RNSCDZ9RMLvhmzkivlbrzszM6hURzJjxLWbN+lck3843GFRzDehFsl5vfwW6gUnATg2Myczs\nda65ZiGzZz/OtdcuyjsU6yc1n4KT9A/AoojYvzEh5cun4Mxay5w5V3DBBVfyyit7sXz51xk79nQ2\n3fReTj75aE466ZN5h2c08BRcCW/k9aNWm5k1xNSpxzJ8+LaceupiQKxf381ZZ01jypSD8g7N+qia\na0D3F8xuAvwj4Os/ZtYUkpDEunXrGTfuFFat6t6wzAa2alpAHy2YfhVYmx4CZ2bWFMuXr2Lu3Ikc\nccSBXHvtIpYvX9V7JWt5Za8BpXt/yoqIvzQkopz5GpCZWW0acQ3obrJu1yIbWfrZtHwb4E/AzrXu\nzMzMrEfZbtgRsVNE7AzcDBwaEdtGxLbAIWmZmZlZ3aq5D+g9EXFjz0xE3AS8t3EhmZm1toigs/Mc\nfLq+b6pJQI9JOj09fXRnSf/Oa4+/NjNrO74ptn9Uk4COIet6/Qvg2jR9TCODMjNrRXPmXMH48Ydy\n2mm38sIL5zJjxmLGjz+UOXOuyDu0AanXbtgR8QxwchNiMTNrab4ptn+VTUCSzo+IL0m6vsTqiIjJ\nDYzLzKzl+KbY/lWpBfTj9PPbJdb5ypuZtSXfFNt/ahqMNN2cOjoi7mtcSPnyjahm1gyD6fESDXsk\nt6QuSVun5PMH4IeSvlNPkGZmlnFPuup6wQ2LiOeBI4AfR8QE4IDGhmVmNji5J91rqklAQyRtDxwJ\n3JCW+RyVWRP4hsfBZ+rUY5k58wusX99NT0+6M8+cxtSpx+YdWtNVk4C+BiwEHomIOyXtAixvbFg2\nUPkLs3/5NM3gU9yTbt26l9u2J12vCSgiro6IPSPi82n+kYiY0vjQbCDyF2b/8Gmawa2nJ93Spd9m\n7tyD27YnXa+94CTtDswGRkbEeEl7ApMj4uvNCLDZ3AuuPn5scv+KCObNW8Cppy5m1apZjBkzg3PP\n3Z8pUw5qy/+UrbU1rBcccDFwGvD3NH8/HorHivi8dv9qldM0PqVqjVRNAnpDRNzRM5OaB/3yRFRJ\nEyUtk7Rc0vQyZS5I6++VtHdvdSUNl3SzpIclLZI0rGDdjFR+maQD++M9WKZVvjAHk1Y4TeNTqtZQ\nEVHxBdwE7AosSfMfB27qrV4V2x0CrAB2AjYF7gH2KCozCbgxTe8D3N5bXeAc4N/S9HTgm2l6XCq3\naaq3AtikRFxh9TnrrB/EvHkLoru7O+bNWxCzZl2cd0hWp4suujzGjTskxo49LaA7xo49LcaNOyQu\nuujyvEOzFpS+N2vOA70ORgpMA34A7C7pMeBRoD/Oq0wAVkTESgBJVwKHAQ8WlJkMXJaywh2Shkka\nSfY01nJ1JwP7p/qXAV1AZ1r/s4h4BVgpaUWK4fZ+eC8GzJjxuQ3THpxxYPOgm9YM1fSCeyQiPkz2\nGIbdgf3IWiN9tQNQeE5hdVpWTZlRFeqOiIi1aXotMCJNj0rlKu3PzPApVWuOSqNhbwmcBOwCLAUu\nImtFfIPs9NXP+7jvaq9qVvMXr1Lbi4iQVGk/JdfNnDlzw3RHRwcdHR1VhGCtJAbROFt58aCb7a3S\nMdTV1UVXV1f/7KTUi+zhc5eSJaFrgDuBxcDb6znXV2L7+wILCuZnANOLylwEHF0wv4ysRVO2bioz\nMk1vDyxL051AZ0GdBcA+JeLq6+lQawFXX31TbLXVv8S8eQvyDqUldHd3x/TpZ0d3d3feodgAUcsx\nRJ3XgColiPsKpocATwJb1LOTMtsfCjxC1iFgM3rvhLAvr3VCKFuXrBNCTzLq5PWdEDYju4b0COk+\nqKJ91vhrslbii+eltXJCdnJsLfUcQ41IQEsqzffHCzgYeIjslN6MtOwk4KSCMhem9fcC76hUNy0f\nDvwKeBhYRDaYas+601L5ZcBBZWKq/TdmLaO7uzuuuurGGDOmMyBizJjOuPrqm9r2y20gJORGJEcn\ntfrVcwzVm4AqdULYU9ILPS/gbQXzz1eoV7WIuCkido+IXSNiVlo2JyLmFJSZltbvFRF3V6qblv8l\nIg6IiN0i4sCIWFew7qxU/q0RsbA/3oO1Fl8831geNwhHVHfzaiOHG/L9S/Vr5jFUNgFFxJCI2Krg\nNbRgeut+j8Ssn7TCDZytIo+EXO2XfyOSo8fQ6x/NOoZqeiJqO/BYcDbYzJp1MbvttuNGvdk6Oz/b\n7/upZzzAefMWcMIJCxkzRqxa1c3cuQf36V6jCI+hl4d6x4Kr5kZUMxvAmnWDcD03r/Z3V+/iFt+q\nVd1tfQq21TkBmVm/qOfLvxHJ0fcvDRw+BVfEp+DM6tes033WWuo9BecEVMQJyMysNo18HpCZmbW4\naru/txInIDOzQWAg3vvkBGRmbWsgthqKDeR7n5yAzKxtDcRWQ7E8RrvoL05AZtZ2BnKrodhAHn7K\n9wGZWdsZbE98Haj3PjkB5SjCD00zy8NgGzGhWaNd9DefgsvRYDj/bDZQedDa/PlG1CLNuBG1nkEb\nzcxalQcjHUAG2/lnM7N6+BRcDgZyrxUzs/7iFlBOBmqvFTOz/uJrQEU8GKmZgXup1sKDkZqZ9SP3\nUm08JyAzswKDaZSEVudrQGZmBdxLtXncAjIzK+Beqs3jFpCZWRH3Um0O94Ir4l5wZjbYNLpHn3vB\nmZlZSa3ao88JyMxskGr1Hn2+BmRmNki1eo++XFpAkoZLulnSw5IWSRpWptxEScskLZc0vZr6kmak\n8sskHViw/BuS/izphca+OzOzTETQ2XkOeV1XbvUefXmdgusEbo6I3YBb0vxGJA0BLgQmAuOAYyTt\nUam+pHHAUan8RGC2XvukfwlMaNg7MjMr0grXXlr5uUe59IKTtAzYPyLWShoJdEXEW4vKvAc4IyIm\npvlOgIj4Zrn6kmYA3RFxdqqzAJgZEbcXbPeFiNiqQmzuBWdmfdJuz/waaL3gRkTE2jS9FhhRoswO\nQGGqXp2WVao/KpUrVcfMrCmmTj2WmTO/wPr13fRceznzzGlMnXpsVfXzPnXXLA1LQOkazf0lXpML\ny6XmRqlPuXiZSpWrUL/cdszMGqqv115a4dRdMzSsF1xEfKTcOklrJY2MiCckbQ88WaLYGmBMwfzo\ntAygXP1Kdao2c+bMDdMdHR10dHTUugkza3P1jKZQeOou6zZ9Ol/96nc3OnXXCo+J6Orqoqurq8/b\nyesa0DnAMxFxdrq2MywiOovKDAUeAj4MPAbcCRwTEQ+Wq586IfyUrLPBDsCvgF0LL+r4GpCZtaqI\nYN68BZx66mJWrZrFmDEzOPfc/Zky5aANyWbevAWccMJC5s6d2DLdqQfaNaBvAh+R9DDwoTSPpFGS\nbgCIiFeBacBC4AHg5xHxYKX6EfEAcFUqfxPwzz3ZRNI5klYBW0haJemrzXmrZmbVqXTqrtVvKq2H\nx4Ir4haQmeVp1qyL2W23HTc6ddfZ+dmqWkd5qbcF5JEQzMxayIwZn9swXXiKrbh1tGpVd0vdVFoP\nJyAzswFisD0mwqfgivgUnJlZbQZaJwQzM2tzTkBmZpYLJyAzM8uFE5CZmeXCCcjMzHLhBGRmZrlw\nAjIzs1w4AZmZWS6cgMzMLBdOQGZmlgsnIDMzy4UTkJmZ5cIJyMzMcuEEZGZmuXACMjOzXDgBmZlZ\nLpyAzMwsF05AZmaWCycgMzPLhROQmZnlwgnIzMxy4QRkZma5cAIyM7NcOAGZmVkunIDMzCwXuSQg\nScMl3SzpYUmLJA0rU26ipGWSlkuaXk19STNS+WWSDkzLtpB0g6QHJS2VNKvx79LMzCrJqwXUCdwc\nEbsBt6T5jUgaAlwITATGAcdI2qNSfUnjgKNS+YnAbElKdc6JiD2AvYH3SZrYqDfXDF1dXXmH0KuB\nECM4zv7mOPvXQImzHnkloMnAZWn6MuDwEmUmACsiYmVEvAJcCRzWS/3DgJ9FxCsRsRJYAewTES9H\nxG8A0rbuBnbo37fUXAPhj3IgxAiOs785zv41UOKsR14JaERErE3Ta4ERJcrsAKwqmF/Na0mjXP1R\nqVypOgCk03UfJWs5mZlZToY2asOSbgZGllj174UzERGSokS54mUqsaxS/ddtR9JQ4GfA+amFZGZm\neYmIpr+AZcDINL09sKxEmX2BBQXzM4DpleqTXQvqLKizgOwUXM/8j4Dzeokt/PLLL7/8qu1VTy5o\nWAuoF/OB44Gz08/rSpS5CxgraSfgMbLOBcf0Un8+8FNJ55KdehsL3Akg6evA1sCJlQKLCFVab2Zm\n/UPpv/7m7lQaDlwF7AisBI6MiHWSRgEXR8QhqdzBwHnAEOCSiJhVqX5adxpwAvAq8KWIWChpNPBn\n4EHg7ymM70bEj5rwds3MrIRcEpCZmVnbjYQgaXNJd0i6R9IDpW5KlbSdpAWpzFJJn84h1J5Yhkha\nIun6MusvSDfe3itp72bHVxBH2TglHZviu0/SbyXtmUeMKZaKn2cq825Jr0o6opmxFcXQ2++9I61f\nKqmryeH1xFDpd95Kx9DK9Le3RNKdZcrkfhz1FmcrHEfVfJapXFXHUF7XgHITEeslfTAiXkq94m6T\ntF9E3FZQbBqwJCJmSNoOeEjSFRHxag4hfwl4ANiqeIWkScCuETFW0j7A98k6b+ShbJzAfwEfiIjn\n0g3AP6A14+y5Afpssg4seV4PrPR7HwZ8DzgoIlanv9E8VPosW+kYCqAjIv5SamULHUcV46Q1jqPe\nYqzpGGq7FhBARLyUJjcju75U/GE+TtZhgfTzmTwOnHTtahLwQ0r/IjfckBsRdwDDJJW6p6qheosz\nIn4XEc+l2TuA0U0Mb4MqPk+ALwLzgKeaFVexKuL8BHBNRKwGiIinmxgeUFWMLXEMFaj0RdgSx1FS\nNs5WOY4d9EQ/AAAGiklEQVTo/R+zqo+htkxAkjaRdA/ZTay/jogHiopcDIyX9BhwL9l/enn4DvCv\nQHeZ9aVu1s3jj7K3OAudCNzY2HDKqhinpB3IRtP4flqU1wXS3j7PscBwSb+WdJekTzUvtA16i7FV\njiHIfo+/Sp/V50qsb5XjqLc4C+V1HFWMsdZjqC0TUER0R8Tbyf7IPiCpo6jIacA9ETEKeDvwPUkl\nT9k0iqRDgScjYgmV/+MoXtfUL80a4kTSB8l6KE6vVK4RqozzPLL7yCKVafopuCrj3BR4B1kL5CDg\nK5LGNinEamPM/Rgq8L6I2Bs4GPiCpPeXKJPrcZRUE2euxxG9x1jTMdSWCahHas7eALyraNV7gatT\nmUeAR4Hdmxsd7wUmS3qUbPSGD0n6cVGZNcCYgvnRaVkzVRMn6YLpxcDkiHi2yTFCdXG+E7gylZlC\nNpjt5BaMcxWwKLIxDp8BFgN7tViMrXAMkfb/ePr5FPALsnEmC7XCcVRNnLkfR1XEWNsx1MwREFrh\nBWwHDEvTW5AdvB8uKnMucEaaHkHWJB+eY8z7A9eXWD4JuDFN7wvcnvNnWy7OHckGht03799/pTiL\nyswFjmjFOIG3Ar8iu375BuB+YFyLxdgSx1D6fLZK028EfgscWFQm9+OoyjhzPY6qibGofK/HUNv1\ngiMbuucySZuQtQAvj4hbJJ0EEBFzgLOAuZLuTWX+LSr0+miSACiMMyJulDRJ0grgReAzeQaYvC5O\n4KvANsD3lT0d45WIeN1/d01WKs5WVOr3vkzSAuA+smswF8frr2PmGiOtcwyNAH6R/u6GAj+JiEUt\neBz1Gif5H0fVxFgT34hqZma5aOtrQGZmlh8nIDMzy4UTkJmZ5cIJyMzMcuEEZGZmuXACMjOzXDgB\nWVuS9Nc+1O2StCw9auA2Sbv1Uv5XOQ5D029q/cwk7SnpkkbFYwOfE5C1q77cABfAJyIbT/Ay4Fvl\nCkr6EPBQRLzQh/01TBo6v1o1fWYRcR+wi6R/rC0qaxdOQNbWlPmWpPvTg7aOTMs3kTRb0oOSFkm6\nQdKUEpu4Fdi1wi4+AfwybXOn1HKaK+khST+RdKCyh4s9LOndqdwbJf1I2YMT7+4ZSyvVXyzpD+n1\nnrR8+7R8SXof70vLN7RYJH1c0tw0famkiyTdDpwtaRdJN6URjhdL2j2V21nS79Ln8vU6P+KbgH+q\ns64Nck5A1u6OIBvIc0/gAOBbkkam5W+JiD2ATwHvYeMWQM8ovx8lGxKnnPcBdxXM7wL8f7Lx3HYH\njoqI9wFfJhtBGuDfgVsiYh/gQymmN5A9PuQjEfFO4GjgglT+E8CCyEYp3ovs8QcUxVvcehkFvCci\nvkz2YLMvRsS7yB6xMDuVOR/4XkTsCTxW4T1WcifwgTrr2iDXjmPBmRXaD/hpZGNSPSnpN8C7yRLH\nVQARsVbSrwvqCPiJpJfJRnn+YoXtjyoaA+3RiPgjgKQ/kg0qCrAU2ClNHwh8VNKX0/w/kI3W/ARw\noaS9gP8mey4QZF/yP5K0KXBdRPQkoHICuDoiQtKWZMn16jTGF2QPaoRsROuPpekryJ5yWavHC96X\n2UacgKzd9Ty3pJRyy3uuAd1dx/7+VjDdDfy9YLrweDwiIpZvFIw0E3g8Ij6Vrt2sB4iIW9NzWQ4F\nLpV0bkRczsatni2K4uh5KvAmwLrUeuoTSYcDZ6TZE9PnI/p2vc0GMZ+Cs3Z3K3BUuubzZrLTRXeQ\nDTU/JV0jGgF0FNWr9mF1j0natsaYFgInb9iR1JMctiZrBQEcR/Y4BiTtCDwVET8ELgF6yq+V9NY0\n8vvHKJEIIuJ54FFJH0/bUnrmDGSfwdFp+tjCepKWldjWdRGxd3r1JOftgT9V/c6trTgBWbsKgIj4\nBdk1nHuBW4B/jYgngWvInmHzAHA5cDfwXHF9AEmjJN1QZj+3sfEDD4uTQKnrNP8P2DRd/F8KnJmW\nzwaOV/Y4+d2Bnk4GHwTukXQ32QX/89PyTuA/yBJJ8TWcwv0eC5yYtrsU6HmA2JfInnp5H9k1o57H\nLmxX5r2WMoHsmVtmr+PHMZiVIemNEfFiasHcAbw3JadattFB1tHg842IMQ+SDgF2jogLqyjbBRxZ\n6+dm7cHXgMzK+w9Jw8guyn+tni/RiOiSdLqkrVr1XqBaRUS51t5G0qm8FU4+Vo5bQGZmlgtfAzIz\ns1w4AZmZWS6cgMzMLBdOQGZmlgsnIDMzy4UTkJmZ5eJ/ANGubgDXMvnQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe925f68fd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "yhat=NLRfun(X,*beta)\n",
    "plt.plot(y,y-yhat,'b*')\n",
    "plt.xlabel('logP. (measured, -)')\n",
    "plt.ylabel('Residual, log(Pa)')\n",
    "plt.title('RMSE=%.1f Pa' %(np.sqrt(sum((10**y-10**yhat)**2))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try the regression again without the log-transformed pressure observations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A: 20.81\n",
      "B: -2797.54\n",
      "C: -51.85\n"
     ]
    }
   ],
   "source": [
    "def NLRfun(X,A,B,C):\n",
    "    T=X\n",
    "    P=np.exp(A+B/(T+C))\n",
    "    yhat=P\n",
    "    return yhat\n",
    "\n",
    "X=data['Temperature'].values # Temperature and pressure are changed to column vectors\n",
    "y=data['Pressure'].values\n",
    "Beta0=[10, -2000, 0]\n",
    "Beta_names=['A: ','B: ','C: ']\n",
    "\n",
    "beta,covar = curve_fit(NLRfun,X,y,p0=Beta0)\n",
    "for b in zip(Beta_names, beta):\n",
    "    print(b[0]+'{:.2f}'.format(b[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x7fe925e61390>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAEZCAYAAACjPJNSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHoNJREFUeJzt3XuUXGWZ7/HvT0BBQWIYDwJGgxocOhMEL8FRmfSc8UAQ\nBoQoF2FAQaIi4oXlIclB0zgK6gxxBs8AkUsIcAQhhpuamOiiBXUkoiEkBCRwjJMECKgEgpJjoJ/z\nx347KSrd1VWdrrduv89atbJvtfezd6rrqfey362IwMzMrN5e0ugAzMysMzjhmJlZFk44ZmaWhROO\nmZll4YRjZmZZOOGYmVkWTjhmZpaFE461HEmrJf1Z0kZJj0u6VtIrS9ZfLalP0lFl7/tGWn5qmn+p\npIskrUn7+q2kbwxynP7XxTXG+hpJt0lal479urL195ftf7Ok20rWv0fSLyU9LekRSWdUOFbFfZVt\n253i2SjpGUkPSvpwLedmVisnHGtFARwZEbsBbwEmAOeVrX8IOKV/gaQdgeOAh9N6gOnAW4F3pH11\nA78e6Dglr7NrjLUP+AEwZcATiRhfun9gDXBjinkH4GbgWxGxO3A8MEvSAbXuaxDr0ravBM4FLpe0\nf43nZ1Y1JxxraRGxHlgEjC9bdTvwHkmj0vxkYBmwvmSbtwO3RMTjaV+/i4hrRzi+JyLiMuCeobaV\nNAn4K+C7adGewB7AtWlf9wAPAEMmhQH2NVSctwJPAftLOkLS0lSq+i9JM6vZh9lQnHCsVQlA0msp\nksndZes3AbcCJ6T5U4Bryrb5BfA5SZ+QNEGSBjvONguLqq6nKrzeNYxzOhWYFxHPAUTEo8B9wGmS\ndkj7fD3w01r3VYmkl0g6BhgFLAeeBU5OpaojgE9IOnoY52P2Ik441ooE3CLpGeC/gEeALw+w3TXA\nKZJ2B/4OuKVs/YXA14CTgF8CayWdUrK+/zilieR0gIj4aUS8qsLr5zWdkPRyimq3q8tWTQXOp0ig\nPwFmRMS6Ye6r3N6SngKeBL5AkWRWRcRPIuL+dJ7LgRuASbWcj9lAnHCsFQVwdGp76Ab+O0X12Iu2\niYifAa+maN+5PSI2lW3QFxGXRMR7gN2BrwBXSXpz2XFKE8mVdTqnY4E/RMSd/Qsk7QN8D/hQROxE\nUW14rqT31bqvQTyazmmPiHhrRPS3HR0s6Q5JT0jaAHyMomrPbLs44VhLS1+q36QoqQzkOuBzbFud\nVr6f/xcRl1C0Y3QNdVxJh5T1CCt/vbvGUzl1gBjfBayNiMUpxoeA7wOHD2Nftfg2RWnwtRExCrgM\nf1fYCPCHyNrBvwETJR2c5sXWtpeLgfdGxF3lb5L0GUmTJO0iacfUXXpXYGnpZgMdMCLuKuu9Vv76\nWclxdgZ2TrM7p/nSOF5LUVKbW3aYFcCbJf29Cm8EjqTo/DCgCvuqxa7AUxHxF0kTgQ+xtWef2bA5\n4VjLi4jfU3zBntu/KL2IiKci4o5B3von4CLgMYp2jE8AUyJidck2t5eVXKrq9VXmz8AzKaYH03FL\n/RPw84j4bdl5PZBi+g/gaaCXoiPAFQCSTpK0opp9DWKwJHIm8KXURvYF4DtV7MtsSGrkA9gkXUXR\nC+aJiJiQlvUAH6X4AoCikXRBWjcdOA14ATg7IhZlD9rMzIal0QnnEIoumNeUJJyZwMaImFW2bRdF\n3fI7gH2AHwH7RURf3qjNzGw4GlqllurVnxpg1UD15kcD10fE5lTl8TAwsY7hmZnZCGrWNpxPSVom\n6cqSO8X3BtaWbLOWoqRjZmYtoBkTzqXAvsCBFI25F1XY1j1nzMxaxI6NDqBcRDzRPy3pCooxsQDW\nAWNKNn1tWvYikpyEzMyGISIGvA1gpDRdCUfSXiWzx1CM7QRwG3CCiiHl9wXGAUsG2kdE+BXBzJkz\nGx5Ds7x8LXwthnstbrppAbvt9hnmzVvY8Hjr+cqhoQlH0vXAzylublsj6TTga5Luk7SMYvymzwJE\nxEqKodZXAguAMyPXVTKzjjN79nWMH38kM2bcxcaNs5g+/U7Gjz+S2bOva3RoLauhVWoRceIAi6+q\nsP0FwAX1i8jMrDB16kmMHr0H55xzJyA2berjggvOYsqUwxodWstquio1Gznd3d2NDqFp+Fps5Wux\nVaVrIQlJbNiwia6uz7Fhw3NbltnwNPTGz3qQ5Jo2MxsRF154Ofvt9zqOPfZQ5s9fxKpVa5g27aON\nDqsuJBF17jTghGNmZlkSjqvUzMwsCyccMzPLwgnHzMyycMIxM7MsnHDMzCwLJxwzM8vCCcfMzLJw\nwjEzsyyccMzMLAsnHDMzy8IJx8zMsnDCMTOzLJxwzMwsCyccMzPLwgnHzMyycMIxM7MsnHDMzCwL\nJxwzM8vCCcfMzLJwwjEzsyyccMzMLAsnHDMzy8IJx8zMsnDCMTOzLJxwzMwsCyccMzPLwgnHzMyy\ncMIxM7MsnHDMzCwLJxwzM8vCCcfMzLJwwjEzsyyccMzMLIuGJhxJV0laL2l5ybLRkhZLekjSIkmj\nStZNl7RK0oOSDm1M1GZmNhyNLuHMASaXLZsGLI6I/YAfp3kkdQHHA13pPZdIanT8ZmZWpYZ+YUfE\nXcBTZYuPAuam6bnA+9P00cD1EbE5IlYDDwMTc8RpZmbbrxlLCHtGxPo0vR7YM03vDawt2W4tsE/O\nwMzMbPh2bHQAlURESIpKmwy0sKenZ8t0d3c33d3dIxuYmVmL6+3tpbe3N+sxFVHp+zxDANJY4PaI\nmJDmHwS6I+JxSXsBd0TEX0uaBhARX03bLQRmRsTdZfuLRp+TmVmrkUREqJ7HaMYqtduAU9P0qcAt\nJctPkPRSSfsC44AlDYjPzMyGoaFVapKuByYBfyVpDfBF4KvAjZJOB1YDxwFExEpJNwIrgeeBM12U\nMTNrHQ2vUhtprlIzM6tdp1apmZlZG3LCMTOzLJxwzMwsCyccMzPLwgnHzMyycMIxM7MsnHDMzCwL\nJxwzM8vCCcfMzLJwwjEzsyyccMzMLAsnHDMzy8IJx8zMsnDCMTOzLJxwzMwsCyccMzPLwgnHzMyy\ncMIxM7MsnHDMzCwLJxwzM8vCCcfMzLJwwjEzsyyccMzMLAsnHDMzy8IJx8zMsnDCMTOzLJxwzMws\nCyccMzPLwgnHzMyycMIxM7MsnHDMzCwLJxwzM8vCCcfMzLLYsZqNJE0AuoCdgQCIiGvqGJeZmbWZ\nIROOpB5gEjAe+D5wOPBTwAnHzMyqVk2V2geA9wKPRcRHgLcAo+oaFSBptaT7JC2VtCQtGy1psaSH\nJC2SVPc4zMxsZFSTcJ6LiBeA5yXtDjwBjKlvWEBRddcdEQdFxMS0bBqwOCL2A36c5s3MrAVUk3B+\nKelVwOXAPcBS4Od1jWorlc0fBcxN03OB92eKw8zMtpMiYvCV0quBscCqiNggaV/glRGxrO6BSf8X\neBp4AZgdEZdLeioiXpXWC/hj/3zJ+6LSOZmZ2bYkERHlP/JH1KCdBiR9FLgAeAR4g6SpEXFrPYMp\n8+6IeCwlvcWSHixdGREhyZnFzKxFVOql9llgfEQ8KekNwLeBbAknIh5L/z4p6WZgIrBe0msi4nFJ\ne1G0J22jp6dny3R3dzfd3d31D9jMrIX09vbS29ub9ZiDVqlJWhoRBw02X9egpJcDO0TERkmvABYB\n51P0lvtDRHxN0jRgVERMK3uvq9TMzGqUo0qtUsJ5ErierQ33xwM3pPmIiLPrFlTRVnRzmt0R+D8R\ncaGk0cCNwOuA1cBxEbGh7L1OOGZmNWp0wvkwaVSB/kVpvj/hzB3ofY3mhGNmVruGJpxW5YRjZla7\nHAnHg3eamVkWTjhmZpaFE46ZmWVRc8KR9ElJx0uq6tEGZmZmMLwSjoBD2Npt2czMbEjupWZmZg0f\nS+2cCu+LiJhVh3jMzKxNVWqH2Y0X3/jZT4MsNzMzG5Sr1MzMrLFVaiVB7AKcDnQBu5BKNxFxWj0D\nMzOz9lJNL7VrgT2ByUAvxeOln61jTGZm1oaGrFKTdG9EHCjpvog4QNJOwE8j4uA8IdbGVWpmZrVr\nlrHU/pL+fVrSBGAU8Or6hWRmZu2omtECLk/PoTkPuA3YFfhCXaMyM7O2415qZmbWHFVqkmaWvL7Y\n/6pnUNY5IoJp075Ou/xIaLfzMRtJ1bTh/ImiV9qzQB/wPmBsHWOyDvLd7/6QSy55jPnzFzU6lBHR\nbudjNpJqrlKT9DJgUURMqk9I28dVaq1h9uzruPjiG9i8+S2sWvVlxo07j512WsbZZ5/Axz52cqPD\nq1m7nY91nqa48XMArwD2GelArLNMnXoSo0fvwTnn3AmITZv6uOCCs5gy5bBGhzYs7XY+ZvVQzUgD\ny0tmXwL8N+BLdYvIOoIkJLFhwya6uj7HmjV9W5a1onY7H7N6qKaE848l088D6yNic53isQ6yatUa\n5syZzLHHHsr8+YtYtWpNo0PaLu12PmYjbdA2nHTvzaAi4o91iWg7uQ3HzKx2je4W/WvgV+nf3wOr\n0uv3abmZNRl3y7ZmNmjCiYixEbEvsBg4MiL2iIg9gCPSMjNrMu6Wbc2smsE7V0TE3wy1rFm4Ss06\nkbtl2/Zqlm7Rj0o6D7iO4mmfHwLW1TMoM6uNu2VbK6hmpIETKbpC3wzMT9Mn1jMoM6tNebfsDRue\nc7dsazpDlnAi4g/A2RliMbPt4G7Z1uwqdYv+94j4tKTbB1gdEXFUfUMbHrfhtL+IYPr0f+HCCz/v\nX/BmI6TRbTjXpH8vGmCdv9GtYfp7Yr3jHYvcRmHWQip1i/5V+re3/wXcB/wxIn6SKT5rQfW6F2T2\n7OsYP/5IZsy4i40bZzF9+p2MH38ks2dfN6LHMbP6qOZ5OL2SXplGHvgVcIWkb9Q/NGtV9boXZOrU\nk+jp+SSbNvXR3xPr/PPPYurUk0b0OGZWH9X0UhsVEc8AxwLXRMRE4L31DcvqqVVLIO6JZdbaqkk4\nO0jaCzgO+H5a5jacFlZLCaSW5JSjBNLfE2vFiouYM+fwmnpiedgXs8aqJuF8Cfgh8EhELJH0Roox\n1WwYGvmlN5wSSC3JKUcJZPr0M5gy5TAkMWXKYUyb9tGq3+thX8waa8iEExE3RcQBEfGJNP9IREyp\nf2gDkzRZ0oOSVkk6t1FxDFcjv/RqKYEMt3pse0og9eLOBmZNIiIqvoA3Az8G7k/zBwDnDfW+eryA\nHYCHgbHATsC9wP5l20Qzuuyya6Or64gYN25GQF+MGzcjurqOiMsuu3bEjtHX1xfnnvu16OvrG3Sb\nm25aELvt9pno6vps7Lbbp2PevIWD7uvGG38QY8ZMC4gYM2Za3HTTgor7blbtdC5m9ZK+O+v6HV5N\nldrlwAzgL2l+OY0b2mYi8HBErI7iIXA3AEc3KJaaDKd9I2qsfqum9FRtCaSdGujb6VzMWlk1Cefl\nEXF3/0zKhI164uc+QOk35Nq0rOkN50uv2uq3WqqMamkDacbqseFqp3Mxa1lDFYGABcCbgKVp/gPA\ngnoXvQaJZQpwecn8ycA3y7YZZoGy/i644Fsxb97C6Ovri3nzFsaFF14+4Ha1Vr+5ysjMthcZqtSq\neTzBWcC3gDdLehT4LdCoO+3WAWNK5sdQlHJepKenZ8t0d3c33d3d9Y6rKtOnn7FlutKQLLUONV9e\nelqzps9VRmZWUW9vL729vVmPWc1o0Y8A/yBpV4rn4TxLcU/O6vqGNqB7gHGSxgKPAsczQHtSacJp\nRcNJIB4peORFeJBQa1/lP8bPP//8uh9z0ISTEszHgDcCK4DLKBrov0LRU+w7dY+uTEQ8L+ksivuC\ndgCujIgHcseRQ60JpNrSk1XPg4SajaxKjyeYDzwD/CdwKEX11Sbg7Ii4N1uENfLjCWx7+XHN1oka\n/XiCN0XEASmQK4DHgNdHxHP1DMis0fy4ZrP6qNQt+oX+iYh4AVjnZGOdwPftmNVHpYRzgKSN/S9g\nQsn8M7kCbAVR4w2a1vx8347ZyBu0DadVNaINZ968hZx22g+ZM2eyq13MrCXlaMOpZqSBjjVUycWD\nQpqZVc8Jp4KhhpbxEyjNzKrXsQmnUuml2pKLG5fNzKrXsQmnUumllpKLG5fNzKrTcZ0Gqr2pr78j\nwJgxYs2aPubMOdwdAsysbTX6xs+2VO1NfR6bzMxsZHVcwql2YEyPTWaN4kFDrV11ZBuO212smVX7\n4D2zVtNxbThmzcqDhlojuQ3HrIN40FBrdx1ZpWbWjHxfl8clbHdOOGZNpNPbF91+1d7chmNmDef2\nq8ZzG46ZdQS3X3UGV6mZWcO5/aozuIRjZk3Bo3u0P7fhmJmZH8BmZmbtwwnHzMyycMIxM8A3XVr9\nOeGYGeCbLq3+nHDMOly1j1Q3217uFm3W4XzTpeXiEo5Zh/NNl5aLSzhm5psuLQvf+GlmllkzPkbc\nN36ambWhTu0R6IRjZiPC9/EMrdN7BDrhmNmI6NRf7bWYOvUkeno+yaZNffT3CDz//LOYOvWkRoeW\nhROOmW2XTv/VXotO7xHoXmpmtl18H09tOrlHoBOOmW2X8l/ta9b0ddSv9lpNn37GlulOS8pNV6Um\nqUfSWklL0+vwknXTJa2S9KCkQxsZp5lt1f+rfcWKi5gz5/CO+tVu1Wu6+3AkzQQ2RsSssuVdwLeB\ndwD7AD8C9ouIvrLtfB+OmVmNOvk+nIFO+mjg+ojYHBGrgYeBiVmjMjOzYWvWhPMpScskXSlpVFq2\nN7C2ZJu1FCUdMzNrAQ3pNCBpMfCaAVb9L+BS4Etp/p+Bi4DTB9nVgHVnPT09W6a7u7vp7u4eZqRm\nZu2pt7eX3t7erMdsujacUpLGArdHxARJ0wAi4qtp3UJgZkTcXfYet+GYdbhmHKus2XVkG46kvUpm\njwGWp+nbgBMkvVTSvsA4YEnu+Mys+XnUg+bUdAkH+Jqk+yQtAyYBnwWIiJXAjcBKYAFwposyZlbK\nox40t6auUhsOV6mZda6IYN68hZxzzp2sWXMhY8ZMZ9asSUyZcpir1obQkVVqZmbD1eljlTU7D21j\nZm2lk8cqa3auUjMzM1epmZlZ+3DCMTOzLJxwzMwsCyccMzPLwgnHzMyycMIxM7MsnHDMzCwLJxwz\nM8vCCcfMzLJwwjEzsyyccMzMLAsnHDMzy8IJx8zMsnDCMTOzLJxwzMwsCyccMzPLwgnHzMyycMIx\nM7MsnHDMzCwLJxwzM8vCCcfMzLJwwjEzsyyccMzMLAsnHDMzy8IJx8zMsnDCMTOzLJxwzMwsCycc\nMzPLwgnHzMyycMIxM7MsnHDMzCwLJxwzM8uiIQlH0gcl3S/pBUlvLVs3XdIqSQ9KOrRk+dskLU/r\n/j1/1GZmtj0aVcJZDhwD3Fm6UFIXcDzQBUwGLpGktPpS4PSIGAeMkzQ5Y7wtqbe3t9EhNA1fi618\nLbbytcirIQknIh6MiIcGWHU0cH1EbI6I1cDDwMGS9gJ2i4glabtrgPfnibZ1+Y9pK1+LrXwttvK1\nyKvZ2nD2BtaWzK8F9hlg+bq03MzMWsSO9dqxpMXAawZYNSMibq/Xcc3MrDkpIhp3cOkO4JyI+HWa\nnwYQEV9N8wuBmcDvgDsiYv+0/ERgUkR8fIB9Nu6EzMxaWERo6K2Gr24lnBqUnuBtwLclzaKoMhsH\nLImIkPSMpIOBJcA/ARcPtLN6XzAzMxueRnWLPkbSGuCdwPclLQCIiJXAjcBKYAFwZmwtgp0JXAGs\nAh6OiIX5Izczs+FqaJWamZl1jmbrpTZskianm0VXSTq30fGMJEmrJd0naamkJWnZaEmLJT0kaZGk\nUSXb13TzrKSXSfpOWv4LSa/Pe4aDk3SVpPWSlpcsy3Lukk5Nx3hI0ik5zreSQa5Fj6S16bOxVNLh\nJeva8lpIGiPpjnTz+ApJZ6flHfe5qHAtmvNzEREt/wJ2oLhnZyywE3AvsH+j4xrB8/stMLps2deB\n/5mmzwW+mqa70vnvlK7Hw2wtyS4BJqbpHwCT0/SZwCVp+njghkafc8l5HgIcBCzPee7AaOARYFR6\nPQKMasJrMRP43ADbtu21oOj9emCa3hX4DbB/J34uKlyLpvxctEsJZyJFu87qiNgM3EBxE2k7Ke8M\ncRQwN03PZeuNsMO5ebZ0X98F/mHkwx+eiLgLeKpscY5zPwxYFBEbImIDsJhi9IuGGeRawLafDWjj\naxERj0fEvWn6WeABik5GHfe5qHAtoAk/F+2ScPYB1pTM998w2i4C+JGkeySdkZbtGRHr0/R6YM80\nPZybZ7dcv4h4Hnha0ugRP4uRU+9z36PCvprRpyQtk3RlSTVSR1wLSWMpSn130+Gfi5Jr8Yu0qOk+\nF+2ScNq958O7I+Ig4HDgk5IOKV0ZRfm23a/BgDr53JNLgX2BA4HHgIsaG04+knal+MX96YjYWLqu\n0z4X6VrMo7gWz9Kkn4t2STjrgDEl82N4ceZtaRHxWPr3SeBmiirE9ZJeA5CKw0+kzcuvxWsprsW6\nNF2+vP89r0v72hHYPSL+WJeTGRn1Pvc/DLCvpvxMRcQTkVDcNjAxrWrrayFpJ4pkc21E3JIWd+Tn\nouRaXNd/LZr1c9EuCeceihGkx0p6KUXD1m0NjmlESHq5pN3S9CuAQylG274NODVtdirQ/0d3G3CC\npJdK2petN88+Djwj6WBJorh59taS9/Tv6wPAj+t8Wtsrx7kvAg6VNErSq4D/Afywnic1HOmLtd8x\nFJ8NaONrkeK+ElgZEf9WsqrjPheDXYum/Vw0qnfFSL8oqpt+Q9EINr3R8Yzgee1L0avkXmBF/7lR\n9BD5EfBQ+o8fVfKeGek6PAgcVrL8bemD9zBwccnyl1HccLuKov53bKPPuyS264FHgb9Q1CN/JNe5\np2OtSq9Tm/BanEbRuHsfsIziC3bPdr8WwHuAvvQ3sTS9Jnfi52KQa3F4s34ufOOnmZll0S5VamZm\n1uSccMzMLAsnHDMzy8IJx8zMsnDCMTOzLJxwzMwsCycca1qSXkhDqy+XdKOkXRoQwwRJV+U+7kiT\n9GFJ3xxim25JT6drvlLSF4fYflb5MEtmlTjhWDP7c0QcFBETKG52/HjpyjTMRr19nmJcqqYkaaT/\nhu+MYty+twMnSzqowraXUlwfs6o44ViruAt4k6RJku6SdCuwQtJLJP2LpCVpZNypUAztIenOkhLS\nu9O2V6f5+yR9ptIBJb0MeGdE/DLN90iam/a7WtKxkv417WtBfwJU8SCr3jS698KS8b3OSHHeK2le\nf4lN0gdTTPdK6k3LXlQikfQ9SX+Xpp9Nx70X+FtJJ0u6O53rZf1JSNJHJP1G0t3Au2q52BHxZ+BX\n6Zp/IcW9XNLskm1WAWNV8qAzs0qccKzppS/y91EM1QHFEOxnR8RfAx8FNkTERIoBCs9QMUz7icDC\n9Gv9AIohPg4C9o6ICRFxADBniEMfRDFcUql9gb+neEbIdcDitK/ngCNUDKT4TWBKRLw9HeMr6b3f\njYiJEXEgxXNLTk/LvwAcmpYflZaVDwFSOv9y4Bdp+z8CxwHvSufaB5yUxtLqoUg076F48FbVw4qo\nGH7+nRTDKf3vFPcEYBdJR5ZsuhT422r3a50tR5WE2XDtImlpmr4TuAp4N8Vgg79Lyw8FJkj6QJp/\nJfAm4JfAVSkB3BIRyyQ9ArxB0sXA9ynG26rk9RRDu/cLYEFEvCBpBfCSiOgfrHA5xRMU9wPGUzy/\nCIqn0T6atpkg6cvA7hRPZ1yYlv8MmCvpRmB+FdflBYrRgaF4GNbbgHvS8XYGHqdIvr1RjOqLpO+k\n2IZyiKRfUySuCyPiAUlTJH2eItGNBu4Hvpe2fzSdt9mQnHCsmT2XfrVvkb5U/1S23VkRsbj8zalB\n+0jgakmzIuJaSW+heFLhxylKBqeXv69EsO1TE/8CEBF9kjaXLO+j+HsScH9EDFSFdTVwVEQsl3Qq\n0J329QlJE4EjgF9JehvwPC+ugdi5ZHpTvHgQxLkRMaPs3MufeDvQ0x8HcldE/GPJfnYG/gN4W0Ss\nkzSzLBbRQc+dse3jKjVrdT8EzixpP9lPxSMdXgc8GRFXUDwP5K2pmmiHiJhPUY311iH2/TuKZ8bX\n4jfAqyW9M8Wzk6SutG5X4PFU6jq5/w2S3hgRSyJiJvAkxbNIVgMHqjCGrc8zKfdj4AOSXp32NTqd\n+93ApDS/E/BBUmKQdIykC6o8n/7k8gcVD/nasp9krxSr2ZBcwrFmNtAv5/InOV5BUaXzaxXFnyco\nnv/RDXw+lUI2AqdQPP52TknPrmkAkj4GEBGzebFlwJsrxLRNO0tEbE7VexdL2p3ib+wbwEqKJHc3\nRVK5myIBAXxd0jiK0sKPIuK+FNdv0/seoGjA3+a4qcrrPGBROq/NwJkRsURSD/CfwAaKtpZ+bwSe\nZlvbPCUzIjZIupyiLefxFHepg4CzB9iX2Tb8eAKzCiRdDVwaEeVftC1L0rXAZ/rbd7ZjP/sB/xoR\nRw25sRlOOGYVSfob4JyI+EijY2k2kmYB8yPip42OxVqDE46ZmWXhTgNmZpaFE46ZmWXhhGNmZlk4\n4ZiZWRZOOGZmloUTjpmZZfH/AXyHjZZF7kXOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe925eb61d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "yhat=NLRfun(X,*beta)\n",
    "plt.plot(y,y-yhat,'b*')\n",
    "plt.xlabel('Press. (measured, Pa)')\n",
    "plt.ylabel('Residual, Pa')\n",
    "plt.title('RMSE=%.1f Pa' %(np.sqrt(sum((y-yhat)**2))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, the residuals look similar to the log-transformed case, but the scale is different, and in terms of obervation error, the RMSE is lower without fitting the log-transformed pressure observations.  The best fit values of $A$ and $B$ have also changed significantly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "#Conclusion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have gone through many calculations common to chemical engineering (optimization, regression, integration) using Python.  Hopefully this has been an educational excersise in presenting Python as an accessable tool to quickly analyze data and models.  Feel free to [download](https://github.com/chepyle/ChE_problemSolving) and edit this workbook as appropriate for your own needs!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}