You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
LLEPE/docs/_Examples/1_getting_started.ipynb

600 lines
42 KiB

{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"LLEPE: Liquid-Liquid Equilibrium Parameter Estimator\n",
"\n",
"Copyright © 2020, UChicago Argonne, LLC. All rights reserved.\n",
"\n",
"Released under the modified BSD license. See LICENSE for more details."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# LLEPE Tutorial - Getting started\n",
"## Introduction\n",
"In this notebook, you will learn how to use LLEPE to fit thermodynamic parameters to experimental data and explore how well the parameters fit.\n",
"## Installation\n",
"Create a conda environment with the following command. The environment name in this example is \"thermo_env\". <br/> \n",
"```$ conda create --name thermo_env python=3.7```<br/>\n",
"Then run the following line to activate the environment"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```$ conda activate thermo_env```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In your terminal run<br/>\n",
"```$ git clone https://xgitlab.cels.anl.gov/summer-2020/parameter-estimation.git```<br/>\n",
"Navigate into the folder with <br/>\n",
"```$ cd parameter-estimation```<br/>\n",
"And run <br/>\n",
"```pip install -e.```<br/>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Moving neccessary files for example"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Copy and paste the xml file named 'elementz.xml' from the 'data/xmls' folder to the Cantera data file located in 'Anaconda3/envs/thermo_env/Lib/site-packages/cantera' "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Import and instantiate LLEPE"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First, you will need to import the package and instantiate LLEPE with a few parameters."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from llepe import LLEPE\n",
"opt_dict = {'Nd(H(A)2)3(org)_h0': {'upper_element_name': 'species',\n",
" 'upper_attrib_name': 'name',\n",
" 'upper_attrib_value': 'Nd(H(A)2)3(org)',\n",
" 'lower_element_name': 'h0',\n",
" 'lower_attrib_name': None,\n",
" 'lower_attrib_value': None,\n",
" 'input_format': '{0}',\n",
" 'input_value': -4.7e6}}\n",
"llepe_parameters = {'exp_data': '../../data/csvs/Nd_exp_data.csv',\n",
" 'phases_xml_filename': '../../data/xmls/twophase.xml',\n",
" 'opt_dict': opt_dict,\n",
" 'phase_names': ['HCl_electrolyte', 'PC88A_liquid'],\n",
" 'aq_solvent_name': 'H2O(L)',\n",
" 'extractant_name': '(HA)2(org)',\n",
" 'diluant_name': 'dodecane',\n",
" 'complex_names': ['Nd(H(A)2)3(org)'],\n",
" 'extracted_species_ion_names': ['Nd+++'],\n",
" 'aq_solvent_rho': 1000.0,\n",
" 'extractant_rho': 960.0,\n",
" 'diluant_rho': 750.0}\n",
"estimator = LLEPE(**llepe_parameters)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Parameters explanation "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### exp_csv_filename"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"exp_csv_filename is the file name for the csv containing experimental data. <br/>\n",
"Let us get the pandas dataframe created by LLEPE."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>h_i</th>\n",
" <th>h_eq</th>\n",
" <th>z_i</th>\n",
" <th>z_eq</th>\n",
" <th>Nd_aq_i</th>\n",
" <th>Nd_aq_eq</th>\n",
" <th>Nd_d_eq</th>\n",
" <th>Nd_org_eq</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>0.01</td>\n",
" <td>0.088304</td>\n",
" <td>1</td>\n",
" <td>0.921696</td>\n",
" <td>0.050001</td>\n",
" <td>0.0239</td>\n",
" <td>1.0921</td>\n",
" <td>0.026101</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>0.01</td>\n",
" <td>0.105094</td>\n",
" <td>1</td>\n",
" <td>0.904906</td>\n",
" <td>0.099998</td>\n",
" <td>0.0683</td>\n",
" <td>0.4641</td>\n",
" <td>0.031698</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>0.01</td>\n",
" <td>0.109017</td>\n",
" <td>1</td>\n",
" <td>0.900983</td>\n",
" <td>0.150006</td>\n",
" <td>0.1170</td>\n",
" <td>0.2821</td>\n",
" <td>0.033006</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>0.01</td>\n",
" <td>0.106012</td>\n",
" <td>1</td>\n",
" <td>0.903988</td>\n",
" <td>0.200004</td>\n",
" <td>0.1680</td>\n",
" <td>0.1905</td>\n",
" <td>0.032004</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>0.01</td>\n",
" <td>0.118934</td>\n",
" <td>1</td>\n",
" <td>0.891066</td>\n",
" <td>0.300011</td>\n",
" <td>0.2637</td>\n",
" <td>0.1377</td>\n",
" <td>0.036311</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" h_i h_eq z_i z_eq Nd_aq_i Nd_aq_eq Nd_d_eq Nd_org_eq\n",
"0 0.01 0.088304 1 0.921696 0.050001 0.0239 1.0921 0.026101\n",
"1 0.01 0.105094 1 0.904906 0.099998 0.0683 0.4641 0.031698\n",
"2 0.01 0.109017 1 0.900983 0.150006 0.1170 0.2821 0.033006\n",
"3 0.01 0.106012 1 0.903988 0.200004 0.1680 0.1905 0.032004\n",
"4 0.01 0.118934 1 0.891066 0.300011 0.2637 0.1377 0.036311"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"estimator.get_exp_df()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The rows are for experiments, and the columns are for the measured quantaties. <br/>\n",
"LLEPE is looking for the ordering of these columns so it is important your experimental file has this ordering. Column names do not matter.<br/>\n",
"Below is a table explaining the meaning of the column headers and the needed column order.<br/>\n",
"If you have more than one rare earth element, append the data to the end in the same order (aq_i, aq_eq, d_eq)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"| Order | Column | Meaning |\n",
"|-------|------------|--------------------------------------------------------------------|\n",
"| 0 | h_i | Initial Concentration of H+ ions (mol/L) |\n",
"| 1 | h_eq | Equilibrium concentration of H+ ions (mol/L) |\n",
"| 2 | z_i | Initial concentration of extractant (mol/L) |\n",
"| 3 | z_eq | Equilibrium concentration of extractant (mol/L) |\n",
"| 4 | \\{RE\\}\\_aq_i | Initial concentration of RE ions (mol/L) |\n",
"| 5 | \\{RE\\}\\_aq_eq | Equilibrium concentration of RE ions in aqueous phase (mol/L) |\n",
"| 6 | \\{RE\\}\\_d_eq | Equilibrium Ratio between amount of RE atoms in organic to aqueous |"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### phases_xml_filename"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is the xml file containing information to be loaded into Cantera, the thermodynamic modeling package. <br/>\n",
"Please see parameter-estimation/data/xmls for file examples. <br/>\n",
"We can explore what has been loaded."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[<cantera.composite.Solution object at 0x00000278E83A15F8>, <cantera.composite.Solution object at 0x00000278E83A1CF8>]\n"
]
}
],
"source": [
"print(estimator.get_phases())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It is a list of two Cantera solutions so we will dig in a little further and see what species these solutions contain."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"HCl_electrolyte\n",
"['H2O(L)', 'H+', 'OH-', 'Cl-', 'Nd+++']\n",
"PC88A_liquid\n",
"['(HA)2(org)', 'dodecane', 'Nd(H(A)2)3(org)']\n"
]
}
],
"source": [
"for phase in estimator.get_phases():\n",
" print(phase.name)\n",
" print(phase.species_names)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can explore Cantera solutions further by visiting https://cantera.org/ and seeing Cantera's documentation."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### opt_dict"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is a dictionary that contains the information about what species and what thermodynamic properties are to be modified.<br/> \n",
"The number after the thermodynamic property is the initial guess for the optimizer. <br/> \n",
"In this example, we chose to optimize the standard enthalpy (h0) of the neodymium-PC88A complex ('Nd(H(A)2)3(org)') and give it an initial guess of -4.7e6. Thus, <br/> \n",
"```python \n",
"opt_dict={'Nd(H(A)2)3(org)': {'h0': -4.7e6}}```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Say we wanted to also modify the extractant ('(HA)2(org)'), but this time change both the standard enthalpy (h0) and the molar volume (molarVolume), then the dictionary would be\n",
"```python \n",
"opt_dict={'Nd(H(A)2)3(org)': {'h0': -4.7e6, 'molarVolume':1.01},\n",
" '(HA)2(org)': {'h0': -4.7e6, 'molarVolume':1.01}}```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### phase_names"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This a list of the phase names in the xml file and can be found in the field phase id."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Names and rhos"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"| Parameter | Meaning | Example value |\n",
"|---------------------|----------------------------------------------|-------------------|\n",
"| aq_solvent_name | Name of solvent in aqueous phase | 'H2O(L)' |\n",
"| extractant_name | Name of extractant in organic phase | '(HA)2(org)' |\n",
"| diluant_name | Name of diluant in organic phase | 'dodecane' |\n",
"| complex_name | Name of rare earth complex in organic phase | 'Nd(H(A)2)3(org)' |\n",
"| rare_earth_ion_name | Name of rare earth ion name in aqueous phase | 'Nd+++' |\n",
"| rhos | Density of species (g/L) | 1000 for 'H2O(L)' |"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For the variables containing \"rho\", these parameters can be left \"None\", and molecular weight and molar volume will be used to calculate density.<br/> However, molar volume values may be wrong and mess up calculations so it is recommended to find density values and replace the default values."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Fitting thermodynamic properties to data"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now that the thermodynamic properties have been set, we now need to set up the optimizer. <br/> The default optimizer is from scipy.optimize.minimize with the arguments below. The optimizer optimizes a value multiplied by the initial guess. <br/> Say $x$ is the variable controlled by the minimizer, the value that is entering the objective function is $x\\times\\mathrm{Guess\\,value}$. So for our case, the values tested are $(4.6\\times 10^6)x$. This is more important for bounds and constraints."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"minimizer_kwargs = {\"method\": 'SLSQP',\n",
" \"bounds\": [(1e-1, 1e1)],\n",
" \"constraints\": (),\n",
" \"options\": {'disp': True, \n",
" 'maxiter': 1000, \n",
" 'ftol': 1e-6}}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With the minimizer arguments defined, we can perform our fit.<br/>\n",
"This minimizes the log mean squared error between the predicted and experimental Distribution ratio (D)."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimization terminated successfully. (Exit mode 0)\n",
" Current function value: 0.025193288852542232\n",
" Iterations: 4\n",
" Function evaluations: 16\n",
" Gradient evaluations: 4\n",
"{'Nd(H(A)2)3(org)_h0': {'upper_element_name': 'species', 'upper_attrib_name': 'name', 'upper_attrib_value': 'Nd(H(A)2)3(org)', 'lower_element_name': 'h0', 'lower_attrib_name': None, 'lower_attrib_value': None, 'input_format': '{0}', 'input_value': -4704699.156668724}}\n",
"0.025193288852542232\n"
]
}
],
"source": [
"est_enthalpy, obj_value = estimator.fit()\n",
"print(est_enthalpy)\n",
"print(obj_value)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can see that the fit function returns an identical structure to opt_dict"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Updating the xml"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we have our new values, let us write them to our original xml to replace the old values"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"estimator.update_xml(est_enthalpy)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Visualization and analysis"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can also see how well this new xml data fits to the experimental data with a parity plot."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAuAAAAI0CAYAAABccZwwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd5hV1fm38XuJFSuWaBAF2883xgRQFDWKXaOJvQDSUbGgKHajUWJJDBp7A8WOYi8oIEhXkA4aWzSoiCQRY0NA6nr/WEMyEsoZmHP2mXPuz3XNNXP22XvPg5fid9Y8ez0hxogkSZKkwlgj6wIkSZKkcmIAlyRJkgrIAC5JkiQVkAFckiRJKiADuCRJklRABnBJkiSpgAzgkiRJUgGtmXUBklTThRC2BJ4HFgCLgFYxxn9kW5UkqVgFB/FI0uoJIdQCYoxxcQihPVAvxnjdMs4bFmM8oND1SZKKiy0okrQKQghXhRDuBIgxLooxLq54a0Pgnewqq3lCCH8KIZyfp3t/EkI4ZBWuGxtC+Hk+apIkA7ikshdCqBNCiCGE0Usd7xFCuGU5l+0CvF3p3EYhhDHAOcDE/FVbWkIIWwBtgR4F+n51QwjTK75eUTi/CbimEDVJKj8GcEmCRsA/gV1CCD9d6vjk5Vzzc+CtJS9ijJNjjE2B3wOXLzkeQtg2hDAshDAMaLTk6xDCttX9h6ih2gP9YoxzC/T9jgQG5HDeS8CBS/37IEnVwgAuSSlojwcGAUfDf/q6fwFMCiGsEUK4PIQwLYQwI4TQAtgR+GvFuetUute3wJwlL2KM02KMB1T0fk9e8nWMcdqKCqpYqX02hDAzhPBxCKFLpfcahxAmhhBmhRCeDCH0CSH8T8/5Cu69TQjhuYp7/3tJK00I4WcVPxx8E0J4J4RwdKVrPgkhXBRCeCuE8G3F9113ZffMwRHA8KXq+ySEcHHF95odQugVQtgyhNC/4s/8WgihTqXzl1v3MhwJ9FtZUTHGH4AJwGE5/jkkKWcGcEmCxqSV7heAYyuO/T+gFvAecBXwW2A/4GfAucA/YoyzKs7dLYQwIoQwFDgfuHF1igkhrAH0BaYAWwMHA+eHEA4PIaxdUeejwKbA08AJVbh3LeBl4FOgQcX9+4QQ1qr4ngOBn1T8GXuHEHaudPnJwK+B7YBfklavl3vPHEv6BfDBMo6fABwK/B9wFNAf+B2wOen/XV0qvncudVPp3GakH7Ry8R7QMMdzJSlnbkMoSWkF/EVgCHBvCGHDimN/BTYBLgIaxhg/BQghvALsueTiGONoUrBboSrsgLIHsEWMcUkP8tQQwn1AC2AusBZwa0zbWD0TQrggx/tSUXdd4OIY48KKY6+HEPYDNgBuqHigdEgI4WWgJdCt4rzbY4wzAEIIfUn/jJZ7zxzr2QSYtYzjd8QY/1XxvUYCX8QYJ1W8fp70QwnAXjnUvUQzYEqlH5xWZhZgC4qkamcAl1TWKtpHfkZqD/k6hDCW1BaxZFX8YOC9GOPfK122JZUewMyD+kDdEMI3lY7VAkaSgu7n8cd7yH5ahXtvA3xaKSgvURf4rNJuLkvuu3Wl1/+s9PWcimtWdM9cfE3aOWZp/6r09dxlvN6gCnUvkVP7SSUbAt+s9CxJqiJbUCSVu11JgW5qxeslbSiNgUmklocvlpxc0cZwLJUewMyDz4CPY4ybVPrYMMZ4JPAPYOsQQqh0flUe6PwM2DaEsPQCzAxgm4r2l8r3/Xw17pmLt0htJquqKnUfCbxShXv/jNQGJEnVygAuqdw1JrUlLFlRfpEU1JasgH8A7BtC+L8QwsbAPaSAl88V8LHAdyGES0MI64UQaoUQdg0h7AGMBhYCXUIIa4YQjqdSO0yO9/4HcEMIYf0QwrohhF8BY4DZwCUhhLVCCAeQeq9z6eVe3j1z0Q/Yvwr1Ly2nukMI2wHrxBjfX+r6tSrqXfKxZsX56wC7k3u/uCTlzAAuqdz9aKvBij7vT0i9yVNijINIYW48MA6YCfwAfJivgmKMi0ghshHwMfAlcD+wcYxxPnA86QHIr4HmwHOVr6/YLeR3K7n3jsA0YDrQvOK+R5Pab74E7gbaLiOw5nzPXOoBHgGODCGst7Lvs5zvnWvdv2HZ7Sf9SL8BWfLRreL40cCwJT3vklSdHEUvSTVcCOEhYHqM8cqsa1kVIYQ/kh6yvDWP36MfcGeMMace8JCGKp0aY/xrvmqSVL58CFOSlKkY4/JWx6vTMGBoridXDFWSpLwwgEuSSl6MsXvWNUjSEragSJIkSQXkQ5iSJElSARnAJUmSpAIqmR7wzTffPDZo0CDrMiRJkiQmTJjwZYxxi2W9VzIBvEGDBowfPz7rMiRJkiRCCJ8u7z1bUCRJkqQCMoBLkiRJBWQAlyRJkgrIAC5JkiQVkAFckiRJKiADuCRJklRABnBJkiSpgAzgkiRJUgEZwCVJkqQCMoBLkiRJBWQAlyRJkgrIAC5JkiQVkAFckiRJKiADuCRJklRABnBJkiSpgAzgkiRJUgEZwCVJkqQCMoBLkiRJBWQAlyRJkgrIAC5JkiQVkAFckiRJKiADuCRJklRABnBJkiSpgAzgkiRJUgEZwCVJkqQCMoBLkiRJBWQAlyRJkgrIAC5JkiQVkAFckiRJKiADuCRJklRABnBJkiSpgAzgkiRJysy330KXLvDNN1lXUjgGcEmSJGVizBho3BjuvhuGDcu6msIxgEuSJKmgFi+GG2+EffeFRYtgxAg49tisqyqcNbMuQJIkSeXjiy+gXTsYMACOPx7uvx/q1Mm6qsJyBVySJEkFMXgwNGwIQ4emtpNnnim/8A0GcEmSJOXZwoVw5ZVw6KGwySYwdiycdRaEkHVl2bAFRZIkSXkzbRqccgq88QZ07Ai33w7rr591VdkygEuSJCkvnn8+he6FC6F37xTEZQuKJEmSqtkPP8A556SHLHfYASZNMnxXZgCXJElStXn/fdhrL7jrLujaFUaNgh13zLqq4mILiiRJklZbjPDww9C5M6y3HvTtC7/9bdZVFSdXwCVJkrRaZs2Ctm2hQwfYYw+YMqUK4bt3b2jQANZYI33u3TuPlRYHV8AlSZK0yiZOhObNYepU+MMf4IoroFatHC/u3Rs6dYI5c9LrTz9NrwFatcpLvcXAFXBJkiRVWYxw222p33vu3DRc56qrqhC+IaX1JeF7iTlz0vES5gq4JEmSquTLL9P2gn37wlFHwYMPwmabrcKNpk2r2vES4Qq4JEmScjZiBDRqBAMGwK23wosvrmL4Bth226odLxEGcEmSJK3UokWpx/vAA9MuJ6NHw3nnreY4+euvh9q1f3ysdu10vIQZwCVJkrRCn38OBx8M3bpBy5bpwcvdd6+GG7dqBT17Qv36KcnXr59el/ADmGAPuCRJklbg5Zehffv0oOVDD6XtBldr1XtprVqVfOBemivgkiRJ+h/z5qVJlkcdBfXqpVXvdu2qOXyXKVfAJUmS9CMffQQtWsCECXDOOXDjjbDuullXVToM4JIkSfqPxx+HM86AtdaC55+HY4/NuqLSk9cWlBDCr0MIH4QQPgohXLaM9y8IIbwbQngrhDA4hFC/0nuLQgiTKz5eymedkiRJ5W727LS3d6tW0LAhTJ5s+M6XvAXwEEIt4C7gCGAXoGUIYZelTpsENIkx/hJ4Buhe6b25McZGFR9H56tOSZKkcjdlCjRpkh6yvOIKGDas5LfizlQ+V8D3BD6KMU6NMc4H+gDHVD4hxjg0xrhk/uibQL081iNJkqRKYoS774amTeGbb2DQILjuOljTJuW8ymcA3xr4rNLr6RXHludUoH+l1+uGEMaHEN4MISzzFyAhhE4V54yfOXPm6lcsSZJUJr7+Gk48ETp3TsN1pkxJe30r//L5882yNqmJyzwxhNZAE2D/Soe3jTHOCCFsDwwJIbwdY/z7j24WY0+gJ0CTJk2WeW9JkiT92KhRaaDOjBlph5MLLoA13Jy6YPL5j3o6sE2l1/WAGUufFEI4BLgCODrGOG/J8RjjjIrPU4FhQOM81ipJklTyFi+GP/0JmjWDWrXgjTfgoosM34WWz3/c44CdQgjbhRDWBloAP9rNJITQGOhBCt9fVDpeJ4SwTsXXmwO/At7NY62SJEkl7Z//hMMPh9/9Dk44ASZNgj33zLqq8pS3FpQY48IQwjnAq0At4IEY4zshhGuA8THGl4AbgQ2Ap0MaqzStYseTnwE9QgiLST8k3BBjNIBLkiStgldfTSPkZ82C++6DU091omWW8vqMa4yxH9BvqWNXVfr6kOVcNwr4RT5rkyRJKnULFsCVV0L37vDzn8OQIemzsuUmM5IkSSXo44/Tg5ZjxkCnTnDLLVC7dtZVCQzgkiRJJefpp+G009LXTz4JJ5+cbT36MZ95lSRJKhFz58KZZ6bA/f/+Xxonb/guPgZwSZKkEvDOO7DHHtCjB1xyCbz+Omy3XdZVaVlsQZEkSarBYoRevaBLF9hgAxgwIG03qOLlCrgkSVIN9e236UHL00+HffZJ4+QN38XPAC5JklQDjR0LjRvDM8/AH/+Y9vr+6U+zrkq5MIBLkiTVIIsXw003wa9+BYsWwYgRcPnlabS8agYDeA32wgsvcPrpp3PMMccwcODArMuRJEl59sUX8JvfwMUXw1FHpXHy++yTdVWqKgN4DdCjRw+22morGjZsyA477MAjjzwCwLHHHst9993HQw89xJNPPpnz/Tp27MhPfvITdt111xWeN2DAAHbeeWd23HFHbrjhhpUelyRJ+TNkCDRqBEOHwl13wbPPwqabZl2VVoUBvAZ466236NatG1OmTOGJJ57gggsu+NH71113HZ07d875fu3bt2fAgAErPGfRokV07tyZ/v378+677/LEE0/w7rvvLve4JEnKj4UL4fe/h0MOgY02SpMtzz4bQsi6Mq0qA3gN8Pbbb7PzzjsDsN1227H22msDEGPk0ksv5YgjjmC33XbL+X7NmjVj05X8yDx27Fh23HFHtt9+e9Zee21atGjBiy++uNzjkiSp+k2bBgccANddB+3bw4QJ0LBh1lVpdbkPeA2wJIDHGLnzzju5/vrrAbjjjjt47bXX+Pbbb/noo48488wzAdhvv/2YNWvW/9znpptu4pBDDsnpe37++edss802/3ldr149xowZs9zjkiSper3wAnTsCAsWwGOPQatWWVek6mIAL3KfffYZs2bN4sgjj+Tzzz/nl7/8Jd26dQOgS5cudOnS5X+uGTly5Gp/3xjj/xwLISz3uCRJqh4//JAesrzzTthtN+jTB3baKeuqVJ0M4EXurbfeolmzZgwZMoSvv/6aXXfdldGjR7PPCh55ro4V8Hr16vHZZ5/95/X06dOpW7fuco9LkqTV98EH0KIFTJ4M558PN9wA66yTdVWqbgbwIvf222/TuHFjAOrUqcMpp5zCK6+8ssIAXh0r4HvssQcffvghH3/8MVtvvTV9+vTh8ccfZ+edd17mcUmStHoeeSQ9XLnuutC3L/z2t1lXpHzxIcwiVzmAAxx11FH069dvte7ZsmVL9t57bz744APq1atHr169/vPekUceyYwZM1hzzTW58847Ofzww/nZz37GySefzM9//vPlHpckSatm1ixo2xbatYPdd0+r34bv0haW1dNbEzVp0iSOHz8+6zIkSZJyNnFiajn5+9/hqqvgyiudaFkqQggTYoxNlvWeK+CSJEkFFiPcfjvsvTfMnp2G7Fx9teG7XNgDLkmSVED//jd06PDfPu8HH4TNN8+6KhWSK+CSJEkFMmJEGqQzYADceiu89JLhuxwZwCVJkvJs0SK45ho48EBYbz0YPRrOO89x8uXKFhRJkqQ8+vxzaN0ahg1L0yzvuQc23DDrqpQlA7gkSVKevPIKtG8Pc+akXu927Vz1li0okiRJ1W7+fLjwwvSQZd26MGFCCuKGb4Er4JIkSdXqo4/S3t4TJkDnznDTTWm6pbSEAVySJKmaPPEEnHFG2s/7uefguOOyrkjFyBYUSZKk1TR7Npx6KpxyCvziF2mcvOFby2MAlyRJWg1vvQVNmqSHLK+4AoYPh/r1s65KxcwWFEmSpFUQI9x7L3TtCnXqwKBBcPDBWVelmsAVcEmSpCr6+ms46SQ4+2w44IDUcmL4Vq4M4JIkSVUwejQ0bgwvvgjdu0O/frDllllXpZrEAC5JkpSDxYvhhhtgv/3Sft6vvw4XXwxrmKZURfaAS5IkrcQ//wlt2sBrr6XWk549YZNNsq5KNZUBXJIkaQUGDkzh+7vvUvA+7TQnWmr1+EsTSZKkZViwAC67DA4/HDbfHMaNg9NPN3xr9bkCLkmStJRPPoGWLeHNN1PovvVWqF0766pUKgzgkiRJlTzzTGoziRH69IHmzbOuSKXGFhRJkiRg7lw466z0kOXOO8OkSYZv5YcBXJIklb1334U990yTLS++GEaOhO23z7oqlSpbUCRJUtmKER54AM49FzbYAPr3h1//OuuqVOpcAZckSWXpu+/glFNSv/fee8OUKYZvFYYBXJIklZ1x49I4+aefhuuuS3t9//SnWVelcmEAlyRJZWPxYvjLX2CffdI+38OHwxVXQK1aWVemcmIPuCRJKgszZ0K7dqnP+9hjoVcv2HTTrKtSOXIFXJIklbyhQ6FhQxgyBO68E557zvCt7BjAJUlSyVq4EK66Cg4+GDbaKE227NzZcfLKli0okiSpJH32Wdrl5PXXoX17uOOOtNWglDUDuCRJKjkvvggdOqQHLR99FFq3zroi6b9sQZEkSSXjhx+gS5f0kOV228HEiYZvFR8DuCRJKgl/+1saqHPHHXD++TBqFOy0U9ZVSf/LFhRJklTjPfIInH02rLMOvPQSHHVU1hVJy+cKuCRJqrG+/x7atk37e++2Wxonb/hWsTOAS5KkGmnSpBS6e/eGq69Oe3zXq5d1VdLKGcAlSVKNEmPq895rL5g9GwYPhm7dYE0ba1VD+K+qJEmqMf79b+jYMfV5/+Y38NBDsPnmWVclVY0r4JIkqUYYORIaNYL+/eGWW6BvX8O3aiYDuCRJKmqLFsG118IBB6RdTkaPTtsMOk5eNZUtKJIkqWjNmJEG6QwdmsbK33MPbLRR1lVJq8cALkmSilK/fml7wTlz4IEHoH17V71VGmxBkSRJRWX+fLjoovSQZd26MH48dOhg+FbpcAVckiQVjb//HVq0SKH77LPhpptgvfWyrkqqXgZwSZJUFPr0gU6doFYtePZZOP74rCuS8sMWFEmSlKnZs+G006BlS9h1V5g82fCt0mYAlyRJmXn7bdhjj/SQ5eWXw/DhUL9+1lVJ+WULiiRJKrgYoUcP6NoVNt4YBg6EQw7JuiqpMFwBlyRJBfXNN3DSSXDWWdCsGUyZYvhWeTGAS5KkgnnzzTRO/sUX4c9/TmPlt9wy66qkwjKAS5KkvFu8OAXuffdN+3mPHAmXXAJrmERUhuwBlyRJefWvf0GbNjBoEJx4Itx3H2yySdZVSdnx505JkpQ3gwZBw4ZpxbtHD3jqKcO3ZACXJEnVbsGCtK3g4YfDZpvBuHFpyI7j5CVbUCRJUjX75JM0VOfNN+H00+HWW6F27ayrkoqHAVySJFWbZ59NUy0XL06j5Zs3z7oiqfjYgiJJklbb3Llw9tnpIcuddoJJkwzf0vIYwCVJ0mp57z1o2hTuuQcuughefx223z7rqqTiZQuKJElaJTHCgw/CueemHu9+/eCII7KuSip+roBLkqQq++47aNUKTj0V9torjZM3fEu5MYBLkqQqGT8edtsNnnwSrrsOBg6EunWzrkqqOfIawEMIvw4hfBBC+CiEcNky3r8ghPBuCOGtEMLgEEL9Su+1CyF8WPHRLp91SpKklVu8GG6+GfbZB+bPh+HD4YoroFatrCuTapa8BfAQQi3gLuAIYBegZQhhl6VOmwQ0iTH+EngG6F5x7abA1UBTYE/g6hBCnXzVKkmSVmzmTDjqKLjwQvjNb2DyZNh336yrkmqmfK6A7wl8FGOcGmOcD/QBjql8QoxxaIxxTsXLN4F6FV8fDgyKMX4VY/waGAT8Oo+1SpKk5Rg2DBo1gtdegzvugOeeg003zboqqebKZwDfGvis0uvpFceW51Sgf1WuDSF0CiGMDyGMnzlz5mqWK0mSKlu4EK6+Gg46CDbYAMaMgXPOcZy8tLryuQ3hsv7zjMs8MYTWQBNg/6pcG2PsCfQEaNKkyTLvLUmSqu6zz9IuJyNHQrt2cOedKYRLWn35XAGfDmxT6XU9YMbSJ4UQDgGuAI6OMc6ryrWSJKn6vfRSajmZOBEeeQQeesjwLVWnfAbwccBOIYTtQghrAy2AlyqfEEJoDPQghe8vKr31KnBYCKFOxcOXh1UckyRJeTJvHpx3HhxzDNSvnwJ4mzZZVyWVnry1oMQYF4YQziEF51rAAzHGd0II1wDjY4wvATcCGwBPh9RQNi3GeHSM8asQwrWkEA9wTYzxq3zVKklSufvb36BFC5g0Cbp0ge7dYZ11sq5KKk0hxtJonW7SpEkcP3581mVIklTjPPoonHVWCtwPPghHH511RVLNF0KYEGNssqz3nIQpSVKZ+v779IBl27ZpsuWUKYZvqRAM4JIklaHJk2H33dPq91VXwZAhUK/eyq+TtPoM4JIklZEY05aCTZumFfAhQ+APf4A187kxsaQf8T83SZLKxFdfQceO8OKLcOSRaXvBLbbIuiqp/LgCLklSGXj99bS3d79+cPPN8PLLhm8pKwZwSZJK2KJFcN11sP/+sPbaMGoUdO3qOHkpS7agSJJUombMgNatYehQaNkS7r0XNtoo66okGcAlSSpB/funLQa//x569YIOHVz1loqFLSiSJJWQ+fPh4ovTQ5ZbbQUTJqQHLw3fUvFwBVySpBIxdWoaJz9uXJps+Ze/wHrrZV2VpKUZwCVJKgFPPgmdOqWV7meegRNOyLoiSctjC4okSTXYnDlw+ulp5fvnP08TLg3fUnEzgEuSVEP99a+wxx7pIcvLL4fhw6FBg6yrkrQytqBIklTDxAg9e8L558PGG8Orr8Khh2ZdlaRcuQIuSVIN8s030Lw5nHkmNGsGU6YYvqWaxgAuSVIN8eab0LgxPPcc3HBD2ut7yy2zrkpSVRnAJUkqcosXQ/fusN9+qf1k5Ei49FJYw/+LSzWSPeCSJBWxf/0L2raFgQPT7ib33w+bbJJ1VZJWhz87S5JUpF57DRo2hBEj4N574emnDd9SKTCAS5JUZBYsgN/9Dg47DDbdFMaOhTPOcJy8VCpsQZEkqYh8+im0bAmjR8Npp8Gtt8L662ddlaTqZACXJKlIPPccnHoqLFoETzyRpltKKj22oEiSlLEffoDOndNDljvuCJMmGb6lUmYAlyQpQ++9B02bwt13w4UXwhtvwA47ZF2VpHyyBUWSpAzECA89BOecA7VrwyuvwJFHZl2VpEJwBVySpAL77jto3Ro6dkyr31OmGL6lcmIAlySpgMaPh912gz594NprYdAgqFs366okFZIBXJKkAogRbrkF9tkH5s2DYcPgyiuhVq2sK5NUaPaAS5KUZ19+Ce3bpz7vY46BXr1gs82yrkpSVlwBlyQpj4YPT+PkBw2C22+H5583fEvlzgAuSVIeLFwI3brBQQfBBhvAm2/Cuec6Tl6SLSiSJFW76dOhVSsYMQLatoW77kohXJLAAC5JUrXq2zf1e8+bBw8/nAK4JFVmC4okSdVg3jw4/3w4+mjYdluYONHwLWnZDOCSJK2mDz9M2wvedht06ZL6vf/v/7KuSlKxsgVFkqTV8NhjcNZZsPba8MILaZtBSVoRV8AlSVoF33+fer3btIFGjWDyZMO3pNwYwCVJqqLJk6FJE3jkEbjqKhg6FLbZJuuqJNUUtqBIkpSjGOHuu+HCC2HTTWHwYDjwwKyrklTTuAIuSVIOvvoKjj8ezjknDdeZMsXwLWnVGMAlSVqJN95Ifd6vvAJ/+Qu8/DJssUXWVUmqqQzgkiQtx6JFcP31sP/+sNZaKYhfcAGs4f89Ja0Ge8AlSVqGf/wDWreGIUOgRQvo0QM22ijrqiSVAgO4JElLGTAgTbH8/nvo1Qs6dIAQsq5KUqnwl2iSJFWYPx8uuQSOOAK23BLGj4eOHQ3fkqqXK+CSJAFTp0LLljB2LJx5Jtx8M6y3XtZVSSpFBnBJUtl76ik4/fS00v3003DiiVlXJKmU2YIiSSpbc+ZAp07QvDnsskuacGn4lpRvBnBJUln6619hjz3gvvvgsstgxAho0CDrqiSVA1tQJEllJUa4/37o0iVtK/jqq3DYYVlXJamcuAIuSSob336b9vTu1An22y+Nkzd8Syo0A7gkqSyMGZPGyT/7LPzpT2mv7622yroqSeXIAC5JKmmLF8ONN8K++6b2k5EjU8+34+QlZcUecElSyfriizTR8tVX4YQTUu/3JptkXZWkcufP/5KkkjR4MDRsCMOGwT33pP29Dd+SioEBXJJUUhYuhCuugEMPhTp1YNy4NNnScfKSioUtKJKkkjFtWhonP2oUnHoq3HYbrL9+1lVJ0o+5Ai5JKgnPP59aTt5+Gx5/PPV7LzN89+6dJu6ssUb63Lt3gSuVVO4M4JKkGu2HH6BzZzj+eNhxR5g0Ka2CL1Pv3mkT8E8/TVuifPppem0Il1RABnBJUo31/vvQtCncfTdceCG88QbssMMKLrjiCpgz58fH5sxJxyWpQOwBlyTVODHCww+nle/ateGVV+DII3O4cNq0qh2XpDxwBVySVKPMmgVt2kCHDrDnnjB5co7hG2Dbbat2XJLywAAuSaoxJkyA3XaDJ56Aa66B116Drbeuwg2uvz4tmVdWu3Y6LkkFYgCXJBW9GOHWW2HvvdNDl8OGwe9/D7VqVfFGrVpBz55Qv37aGLx+/fS6Vat8lC1Jy2QPuCSpqH35ZWo3efllOPpoeOAB2Gyz1bhhq1YGbkmZcgVcklS0hg9Pe3sPHAi33w4vvLCa4VuSioABXJJUdBYtgj/8AQ46KA3TGT0azj3XcfKSSoMtKJKkojJ9OrRunVa/27SBu+6CDTfMuipJqj4GcElS0Xj5ZWjfPj1o+fDD0LZt1hVJUvWzBas//0sAACAASURBVEWSlLl586BrVzjqKNhmm7TdoOFbUqlyBVySlKkPP4QWLWDixNTn3b07rLtu1lVJUv4YwCVJmXn8cTjjDFhrrbTDyTHHZF2RJOWfLSiSpIKbPRs6dkzbcTdqBFOmGL4llQ8DuCSpoKZMgd13h4cegiuvhKFDU9+3JJULW1AkSQURI9xzD1xwAWy6Kbz2WtrnW5LKjSvgkqS8+/prOOEE6Nw5he7Jkw3fksqXAVySlFejRqU+77594aab0l7fP/lJ1lVJUnYM4JKkvFi0CP74R2jWDNZcMwXxCy+ENfw/j6Qyt8Ie8BDCBSt6P8Z4c/WWI0kqBf/8ZxonP3gwNG8OPXrAxhtnXZUkFYeVrUNsWPHRBDgL2Lri40xgl5XdPITw6xDCByGEj0IIly3j/WYhhIkhhIUhhBOXem9RCGFyxcdLuf6BJEnZevVVaNgwrXjffz888YThW5IqW+EKeIzxDwAhhIHAbjHGWRWvuwFPr+jaEEIt4C7gUGA6MC6E8FKM8d1Kp00D2gMXLeMWc2OMjXL7Y0iSsrZgQdpWsHt32HXXtL3gLitdqpGk8pPrNoTbAvMrvZ4PNFjJNXsCH8UYpwKEEPoAxwD/CeAxxk8q3lucYx2SpCL08cfQsiWMGQNnngk33wzrrZd1VZJUnHIN4I8CY0MIzwMROA54ZCXXbA18Vun1dKBpFWpbN4QwHlgI3BBjfGHpE0IInYBOANtuu20Vbi1Jqi5PPw2nnQYhwFNPwUknZV2RJBW3nAJ4jPH6EEJ/YL+KQx1ijJNWcllY1q2qUNu2McYZIYTtgSEhhLdjjH9fqq6eQE+AJk2aVOXekqTVNGcOdO0KPXtC06ap13u77bKuSpKKX1U2g6oNfBdjvA2YHkJY2V+z04HKw4XrATNy/WYxxhkVn6cCw4DGVahVkpRH77wDe+6Zwvell8LIkYZvScpVTgE8hHA1cClwecWhtYDHVnLZOGCnEMJ2IYS1gRZATruZhBDqhBDWqfh6c+BXVOodlyRlI0a47z7YYw+YOTPteHLDDbDWWllXJkk1R64r4McBRwOz4T+r0xuu6IIY40LgHOBV4D3gqRjjOyGEa0IIRwOEEPYIIUwHTgJ6hBDeqbj8Z8D4EMIUYCipB9wALkkZ+vbb9KBlp07wq1/BlClw2GFZVyVJNU+uD2HOjzHGEEIECCGsn8tFMcZ+QL+ljl1V6etxpNaUpa8bBfwix9okSXk2diy0aAHTpsGf/gSXXOJES0laVbn+9flUCKEHsEkI4XTgNeD+/JUlSSoGixfDTTelFe/Fi2HECLjsMsO3JK2OXHdBuSmEcCjwHbAzcFWMcVBeK5MkZeqLL6BdOxgwAI4/Pk21rFMn66okqebLKYCHEP4cY7wUGLSMY5KkEjN4MLRuDV9/DXffnYbrhGVtLitJqrJcf4l46DKOHVGdhUiSsrdwYRonf+ihsMkmqff7rLMM35JUnVa4Ah5COAs4G9ghhPBWpbc2BEblszBJUmFNmwannAJvvAEdO8Ltt8P6OT1yL0mqipW1oDwO9Af+BFxW6fisGONXeatKklRQL7yQQvfChdC7dwrikqT8WGELSozx2xjjJ8BtwFcxxk9jjJ8CC0IITQtRoCQpf374Ac49F447DrbfHiZONHxLUr7l2gN+D/B9pdezK45JkmqoDz6AvfaCO++Erl1h1CjYccesq5Kk0pfrIJ4QY4xLXsQYF4cQcr1WklREYoRHHoHOnWHddeHll+E3v8m6KkkqH7mugE8NIXQJIaxV8XEeMDWfhUmSqt+sWdC2LbRvD02apHHyhm9JKqxcA/iZwD7A58B0oCnQKV9FSZKq38SJsNtu8Pjj8Ic/pL2+t94666okqfzkOgnzC6BFnmuRJOVBjGlLwUsugS22gKFDoVmzrKuSpPK1sn3AL4kxdg8h3AHEpd+PMXbJW2WSpNX2739Dhw7Qty8cdRQ8+CBstlnWVUlSeVvZCvh7FZ/H57sQSVL1GjEibSk4cybcdlvabtCJlpKUvRUG8Bhj34rPDxemHEnS6lq0CK6/PvV5b789jB6der8lScVhZS0ofVlG68kSMcajq70iSdIq+/xzaNUKhg+H1q3h7rthww2zrkqSVNnKWlBuqvh8PLAV8FjF65bAJ3mqSZK0Cl55Bdq1g7lz4aGH0teSpOKzshaU4QAhhGtjjJWfme8bQhiR18okSTmZPx8uuwxuuQUaNoQnn4Sdd866KknS8uS6D/gWIYTtl7wIIWwHbJGfkiRJufroI9hnnxS+zzkH3nzT8C1JxS7XcfJdgWEhhCXTLxsAZ+SlIklSTh5/HM44A9ZaC55/Ho49NuuKJEm5yHUQz4AQwk7A/6s49H6McV7+ypIkLc/s2dClCzzwAPzqVymIb7tt1lVJknKVUwtKCKE2cDFwToxxCrBtCOG3ea1MkvQ/3noLmjRJA3WuvBKGDTN8S1JNk2sP+IPAfGDvitfTgevyUpEk6X/ECPfcA3vuCd98A6+9BtdeC2vm2kgoSSoauQbwHWKM3YEFADHGuYDz1CSpAL7+Gk48Ec4+Gw48EKZMgYMOyroqSdKqyjWAzw8hrEfFUJ4Qwg6APeCSlGejRkGjRvDSS3DjjWmv75/8JOuqJEmrI9cAfjUwANgmhNAbGAxckreqJKnMLV4Mf/oTNGsGtWrBG2/ARRfBGrn+rS1JKlor7R4MIQTgfdI0zL1IrSfnxRi/zHNtklSW/vlPaNMm9XmffDL07Akbb5x1VZKk6rLSAB5jjCGEF2KMuwOvFKAmSSpbAwem8D1rFtx3H5x6KgSfuJGkkpLrLzPfDCHskddKJKmMLViQxskffjhssQWMGwennWb4lqRSlOsGVgcCZ4YQPgFmk9pQYozxl/kqTJLKxccfQ8uWMGYMdOqUxsrXrp11VZKkfMk1gB+R1yokqUw980xa6Y4RnnoKTjop64okSfm2wgAeQlgXOBPYEXgb6BVjXFiIwiSplM2dC127Qo8eabhOnz6w3XZZVyVJKoSV9YA/DDQhhe8jgL/kvSJJKnHvvptCd48ecMkl8Prrhm9JKicra0HZJcb4C4AQQi9gbP5LkqTSFCP06gVdusAGG8CAAemhS0lSeVnZCviCJV/YeiJJq+7bb9ODlqefDvvsk8bJG74lqTytbAW8YQjhu4qvA7Bexeslu6BslNfqJKkEjBsHLVrAp5/CH/8Il17qREtJKmcrDOAxxlqFKkSSSs3ixWlLwcsug7p1YcSItPotSSpvuW5DKEmqgpkzoV076N8fjjsu9X7XqZN1VZKkYuAvQSWpmg0ZAg0bps933QXPPmv4liT9lwFckqrJwoXw+9/DIYfARhulyZZnn+04eUnSj9mCIknVYNo0aNUq7endoQPccQesv37WVUmSipEBXJJW0wsvQMeOsGABPPZYCuKSJC2PLSiStIp++CEN1TnuuDTJcuJEw7ckaeUM4JK0Cj74APbeO7WadO0Ko0bBTjtlXZUkqSawBUWSquiRR9LDleuuC337wm9/m3VFkqSaxBVwScrRrFnQtm3a37tJkzRO3vAtSaoqA7gk5WDSJNh9d+jdG7p1g8GDYeuts65KklQTGcAlaQVihNtvh732gjlz0nCdq6+GWrWyrkySVFPZAy5Jy/Hvf6ftBV96KbWaPPggbL551lVJkmo6V8AlaRlGjoRGjaB/f7j11hTCDd+SpOpgAJekJXr3ZlH97bk2XMUBzRax7oLvGD0azjvPcfKSpOpjAJckgN69mXHaVRwyrRdXcQ0teYKJ3+3E7u/3zroySVKJMYBLEtDvgkE0/OFNxrInD9GOR2nDhnO/gCuuyLo0SVKJMYBLKmvz58OFF8JvvniIusxgArvTjkf4T8fJtGlZlidJKkHugiKpbP3979CiBYwfD503fJibZp3Busz78UnbbptNcZKkkuUKuKSy9MQT0LgxfPQRPPcc3HnPmqxbe6nNvWvXhuuvz6ZASVLJMoBLKiuzZ8Opp8Ipp8AvfgGTJ8NxxwGtWkHPnlC/ftrypH799LpVq6xLliSVGFtQJJWNt9+G5s3h/ffTs5XdusGalf8WbNXKwC1JyjsDuKSSFyP06AFdu8Imm8CgQXDwwVlXJUkqV7agSCpp33wDJ50EZ50F++8PU6YYviVJ2TKASypZo0encfIvvgjdu0O/fvCTn2RdlSSp3BnAJZWcxYvhhhtgv/3S85Svvw4XXwxr+DeeJKkI2AMuqaT861/Qpk3q8z7pJLjvPth446yrkiTpvwzgkkrGoEEpfH/7bdpB8LTT0gq4JEnFxF/ISqrxFiyAyy+Hww6DzTeHcePg9NMN35Kk4uQKuKQa7ZNPoGVLePNN6NQJbrklDbCUJKlYGcAl1VjPPpumWsYITz4JJ5+cdUWSJK2cLSiSapy5c9O+3ieeCDvvnMbJG74lSTWFAVxSjfLee9C0Kdx7b9pacORI2G67rKuSJCl3tqBIqhFihAcfhHPOgQ02gP794de/zroqSZKqzhVwSUXvu++gVavU77333mmcvOFbklRTGcAlFbXx46FxY3jqKbj+ehg4EH7606yrkiRp1RnAJRWlxYvh5pthn33SPt/Dh8Pvfge1amVdmSRJq8cecElFZ+ZMaN8e+vWD446D+++HTTfNuipJkqqHK+CSisrQodCwIQweDHfdlfb6NnxLkkqJAVxSUVi4EK66Cg4+GDbaKE22PPtsx8lLkkqPLSiSMvfZZ2mXk5EjU+vJHXekrQYlSSpFBnBJmXrpJejQAebPh0cfhdats65IkqT8sgVFUibmzYMuXeCYY6BBA5g40fAtSSoPeQ3gIYRfhxA+CCF8FEK4bBnvNwshTAwhLAwhnLjUe+1CCB9WfLTLZ52SCutvf0sDde64A84/H0aNgp12yroqSZIKI28tKCGEWsBdwKHAdGBcCOGlGOO7lU6bBrQHLlrq2k2Bq4EmQAQmVFz7db7qlVQYjz4KZ50F66yT2k+OOirriiRJKqx8roDvCXwUY5waY5wP9AGOqXxCjPGTGONbwOKlrj0cGBRj/KoidA8CHDwt1WDffw/t2kHbtrD77mmcvOFbklSO8hnAtwY+q/R6esWxars2hNAphDA+hDB+5syZq1yopPyaNCmF7sceg6uvhiFDoF69rKuSJCkb+Qzgy9q9N1bntTHGnjHGJjHGJltssUWVipOUfzGmPu+99kor4IMHQ7dujpOXJJW3fAbw6cA2lV7XA2YU4FpJReCrr9IY+S5d4NBDU8vJAQdkXZUkSdnLZwAfB+wUQtguhLA20AJ4KcdrXwUOCyHUCSHUAQ6rOCapBnj9dWjUCPr1g1tugb59YfPNs65KkqTikLcAHmNcCJxDCs7vAU/FGN8JIVwTQjgaIISwRwhhOnAS0COE8E7FtV8B15JC/DjgmopjkorYokVw3XWw//6w9towenTaZtBx8pIk/VeIMde27OLWpEmTOH78+KzLkMrWjBlpkM7QoXDKKXDPPbDRRllXJUlSNkIIE2KMTZb1nqPoJa22/v3T9oJz5sCDD6btBl31liRp2RxFL2mVzZ8PF10ERx4JdevChAnQvr3hW5KkFXEFXNIqmToVWrSAcePg7LPhpptgvfWyrkqSpOJnAJdUZU8+CaefnvbzfvZZOP74rCuSJKnmsAVFUs7mzEnBu0UL2HVXmDzZ8C1JUlUZwCXl5O23oUkT6NULfvc7GD4c6tfPuipJkmoeW1AkrVCM0LNn2s97441h4EA45JCsq5IkqeZyBVzScn3zDZx8Mpx5ZhquM2WK4VuSpNVlAJe0TG++mcbJv/ACdO+exspvuWXWVUmSVPMZwCX9yOLF8Oc/w377pf28R46Eiy+GNfzbQpKkamEPuKT/+Ne/0kTLgQPhpJNS7/cmm2RdlSRJpcUALgmA116D1q3h22+hR4+03aATLSVJqn7+UlkqcwsWpG0FDzsMNtssTbbs1MnwLUlSvrgCLpWxTz+Fli1h9Oi04n3rrVC7dtZVSZJU2gzgUpl69lk47bT00GWfPtC8edYVSZJUHmxBkcrM3Llw9tlw4omw004waZLhW5KkQjKAS2XkvfegaVO45x646CJ4/XXYfvusq5IkqbzYgiKVgRjhwQfh3HNh/fXTUJ0jjsi6KkmSypMr4FKJ++67tL3gqafCXnvB5MmGb0mSsmQAl0rY+PGw227w5JNw3XVpwE7dullXJUlSeTOASyUoRrjlFthnH5g/H4YNgyuugFq1sq5MkiTZAy6VmJkzoUMHeOUVOPZY6NULNt0066okSdISroBLJWTYMGjUCAYNgjvvhOeeM3xLklRsDOBSCVi4EK6+Gg46CDbYAMaMgc6dHScvSVIxsgVFquGmT4dWrWDECGjXLq18b7BB1lVJkqTlMYBLNVjfvtC+PcybB48+mrYblCRJxc0WFKkGmjcPzj8fjj4a6teHiRMN35Ik1RSugEs1zIcfQvPmMGkSnHce/PnPsM46WVclSZJyZQCXapDHHoOzzoK114YXX0wr4JIkqWaxBUWqAb7/PvV6t2kDjRvDlCmGb0mSaioDuFTkJk+G3XeHRx5JWw0OGQL16mVdlSRJWlW2oEhFKka46y646CLYbLMUvA84IOuqJEnS6nIFXCpCX30Fxx8P554LhxySWk4M35IklQYDuFRk3ngjjZN/5RW4+ea01/fmm2ddlSRJqi4GcKlILFoE118P+++fdjkZNQq6dnWcvCRJpcYecKkI/OMfaZDOkCHQsiXcey9stFHWVUmSpHwwgEsZ698f2rWD2bPhgQfSdoOuekuSVLpsQZEyMn8+XHwxHHkkbLUVjB8PHToYviVJKnWugEsZmDo1tZqMHZsmW/7lL7DeellXJUmSCsEALhXYk09Cp06wxhrwzDNwwglZVyRJkgrJFhRpVfXuDQ0apCTdoEF6vQJz5qTg3aIF/PznMGmS4VuSpHLkCri0Knr3Tml6zpz0+tNP02uAVq3+5/S//hWaN4f33oPLL4c//AHWWquA9UqSpKLhCri0Kq644r/he4k5c9LxSmKEnj1hjz3g3/+GV1+FP/7R8C1JUjkzgEurYtq0lR7/5pu06n3GGdCsWRonf+ihBapPkiQVLQO4tCq23XaFx8eMgcaN4fnn4c9/Tnt9b7llAeuTJElFywAurYrrr4fatX98rHZtFl97Pd27w777pkMjR8Ill6TnNCVJksAALq2aVq1Sc3f9+mlyTv36/OvGRziydysuvRSOPTbtcrLXXlkXKkmSio27oEirqlWr/+x48tpr0KZN6vu+9960IYoTLSVJ0rK4Ai6thoUL08Ynhx0GdeqkyZZnnGH4liRJy+cKuLSKPv0UTjkFRo2C006D227737ZwSZKkpRnApVXw3HNw6qmwaBE88USabilJkpQLW1CkKvjhB+jcOY2Q32mn9KCl4VuSJFWFAVzK0fvvQ9OmcPfdcNFF8PrrsMMOWVclSZJqGltQpJWIER56CM45J/V49+sHRxyRdVWSJKmmcgVcWoFZs6B1a+jYMa1+T5li+JYkSavHAC4tx4QJsNtu0KcPXHstDBoEdetmXZUkSarpDODSUmKEW2+FvfdOD10OGwZXXgm1amVdmSRJKgX2gEuVfPkldOgAL78MxxwDDzwAm26adVWSJKmUuAIuVRg+HBo2hIED4Y474PnnDd+SJKn6GcBV9hYtgm7d4KCDYIMNYMyYtOOJ4+QlSVI+2IKisjZ9OrRqBSNGQNu2cNddKYRLkiTliwFcZevll6F9+/Sg5cMPpwAuSZKUb7agqOzMmwddu8JRR8E228DEiYZvSZJUOK6Aq6x8+CG0aJFCd5cu0L07rLNO1lVJkqRyYgBX2ejdG848E9ZeG158EY4+OuuKJElSObIFRSVv9uy0t3fr1tC4MUyebPiWJEnZMYCrpE2ZArvvnh6yvOoqGDIk9X1LkiRlxRYUlaQY4e674cIL0zCdwYPhwAOzrkqSJMkVcJWgr76CE05Iw3QOPjitghu+JUlSsTCAq6S88Ubq8375ZfjLX6BvX9hii6yrkiRJ+i8DuErCokXwxz/C/vvDmmumIH7BBbCG/4ZLkqQiYw+4arx//APatEl93i1aQI8esNFGWVclSZK0bAZw1WgDBqQplt9/D716pe0GQ8i6KkmSpOXzF/SqkebPh0sugSOOgK22gvHjoWNHw7ckSSp+roCrxpk6FVq2hLFj4ayz0sOW662XdVWSJEm5MYCrRnn6aTjttLTS/cwzabtBSZKkmsQWFNUIc+bAGWfAySfDLrukcfKGb0mSVBMZwFX03nkH9twTevaEyy6DESOgQYOsq5IkSVo1tqCoaMUI998P550HG24Ir74Khx2WdVWSJEmrJ68r4CGEX4cQPgghfBRCuGwZ768TQniy4v0xIYQGFccbhBDmhhAmV3zcm886VXy+/Tbt6d2pE+y7bxonb/iWJEmlIG8r4CGEWsBdwKHAdGBcCOGlGOO7lU47Ffg6xrhjCKEF8GegecV7f48xNspXfSpeY8em8D1tGtxwA1x8sRMtJUlS6chnrNkT+CjGODXGOB/oAxyz1DnHAA9XfP0McHAI7uRcrhYvhhtvhF/9Kn09ciRceqnhW5IklZZ8Rputgc8qvZ5ecWyZ58QYFwLfAptVvLddCGFSCGF4CGG/ZX2DEEKn8P/bu/cgrer7juPvb0DFWxANGVsVwYgzXqKiGybWkWqihoYUmoREkCoTjbd4+YPYasfYqWRoamLM6HivMKIjxcuMDGITJghehlRlQcXgxBa1EsQEFMcaEXTh2z/OoV3XBR7Yfc6z++z7NbPDc27P8935cpbP/vidcyJaI6J13bp13Vu9KrV2LYwZUzxcZ9y44i4nJ53U6KokSZK6Xz0DeGcj2VnjPm8BQzJzBDAFmBURn/3Ujpl3ZWZLZrYMHjy4ywWrMR5/HI47Dp54Au64o7jX9377NboqSZKk+qhnAF8NHNJu+WBgzbb2iYj+wEBgfWZuysx3ADJzKfAqcEQda1UDtLXBj34EZ5wBgwYVc78vusjHyUuSpOZWzwC+BBgeEcMiYndgAjC3wz5zgcnl6/HAwszMiBhcXsRJRBwGDAdeq2OtqtiqVXDqqTBtGpx3HixZAl/8YqOrkiRJqr+63QUlM9si4jJgPtAPmJGZKyJiKtCamXOB6cB9EbESWE8R0gFGAVMjog3YDFycmevrVauq9cgjRejevBlmzYKJExtdkSRJUnUis+O07N6ppaUlW1tbG12GtmPjRrjySrj1Vmhpgdmz4QtfaHRVkiRJ3S8ilmZmS2fbvMGbKvG738GXv1yE7x/+EBYvNnxLkqS+yUfRq64yYeZMuPRS2GsveOwx+PrXG12VJElS4zgCrrp5/30491z43vdg5Mji3t6Gb0mS1NcZwFUXy5bBCScUF1lOnQoLFsBBHR/DJEmS1AcZwNWtMuGmm4r53hs3Fg/XufZa6Nev0ZVJkiT1DM4BV7d5++3i9oKPPgpjx8KMGXDAAY2uSpIkqWdxBFzd4qmn4PjjYf58uPlmmDPH8C1JktQZA7i6ZPNmuO46OO204i4nzzwDl1/u4+QlSZK2xSko2mVvvgmTJsGTT8I55xT3+N5330ZXJUmS1LMZwLVLHnsMJk8uLrScObO43aAkSZJ2zCko2imbNsGUKfCNb8AhhxS3GzR8S5Ik1c4RcNVs5UqYMAGWLi3mef/0pzBgQKOrkiRJ6l0M4KrJrFlw0UWw227FHU7GjWt0RZIkSb2TU1C0XR98UNzbe9Kk4jaDL75o+JYkSeoKA7i2aflyaGmBe+4pnma5aFEx71uSJEm7ziko+pRMuP324mLL/feHxx8v7vMtSZKkrnMEXJ/w7rswfjxceil85SvFlBPDtyRJUvcxgOv//OY3xTzvRx+Fn/8c5s2DwYMbXZUkSVJzMYCLLVvgJz+BUaOgf39YvLiYfvIZ/3ZIkiR1O+eA93F/+EPxGPkFC+Css+DOO2HgwEZXJUmS1LwM4H3Y/PnFUyzffx/uvru43WBEo6uSJElqbk4y6IM+/hiuugpGj4bPfx5aW+H88w3fkiRJVXAEvI95/XWYOBGefRYuvhhuvBH23LPRVUmSJPUdBvA+5KGH4PvfL0a6H3wQvvOdRlckSZLU9zgFpQ/48MNitPu734Ujj4Tnnzd8S5IkNYoBvMmtWAEjRxZ3N7nqKnj6aRg2rNFVSZIk9V1OQWlSmTB9OlxxBey7b3HHkzPPbHRVkiRJcgS8Cb33XnGh5QUXwMknF4+TN3xLkiT1DAbwJvPcczBiBDz8cPF0y/nz4cADG12VJEmStjKAN4ktW+CGG4oR7y1birneV1/t4+QlSZJ6GueAN4G1a2HyZPjVr+Bb3yqeajloUKOrkiRJUmccH+3lFi6E44+HRYvgttuKqSeGb0mSpJ7LAN5LtbXBtdfC6afDwIHF3O9LLvFx8pIkST2dU1B6oVWr4OyzYfFiOO88uPlm2HvvRlclSZKkWjgC3hPdfz8MHVpcQTl0aLFcmjOnmHKyfDnMmlXc69vwLUmS1Hs4At7T3H8/XHghbNhQLL/xBlx4IRs/+gx/t2wit9wCJ54Is2fD4Yc3tlRJkiTtPAN4T3PNNf8fvkuvbDiYCRcdywsfw5Qpxf29d9+9QfVJkiSpSwzgPc2qVZ9YvJdz+AG3MeDjjcybB2PGNKguSZIkdQvngPc0Q4YA8D77cC4zmcy9fIklvHjQGMO3JElSEzCA9zTTprFmwGGcyFLuZxJTuZYFe47loOuvaHRlkiRJ6gZOQelpJk3iwC1w0iUvcfcHFzDq0Ddg2h0waVKjK5MkSVI3iMxsdA3doqWlJVtbWxtdhiRJkkRELM3Mls62OQVFkiRJqpABXJIkSaqQAVySJEmqkAFckiRJqpABXJIkSaqQAVySJEmqkAFckiRJqpABXJIkSaqQAVySJEmqkAFcdt92eAAACDpJREFUkiRJqpABXJIkSaqQAVySJEmqkAFckiRJqpABXJIkSaqQAVySJEmqkAFckiRJqpABXJIkSaqQAVySJEmqkAFckiRJqpABXJIkSaqQAVySJEmqkAFckiRJqpABXJIkSaqQAVySJEmqkAFckiRJqpABXJIkSaqQAVySJEmqkAFckiRJqpABXJIkSaqQAVySJEmqkAFckiRJqpABXJIkSaqQAVySJEmqkAFckiRJqpABXJIkSaqQAVySJEmqkAFckiRJqlBdA3hEjI6IVyJiZURc3cn2PSLigXL7sxExtN22fyjXvxIRX6tnnZIkSVJV6hbAI6IfcCvwV8BRwMSIOKrDbucD72bm4cAvgOvLY48CJgBHA6OB28r3kyRJknq1eo6AjwRWZuZrmfkRMBsY12GfccDM8vXDwFcjIsr1szNzU2a+Dqws30+SJEnq1eoZwA8Cft9ueXW5rtN9MrMNeA84oMZjiYgLI6I1IlrXrVvXjaVLkiRJ9VHPAB6drMsa96nlWDLzrsxsycyWwYMH70KJkiRJUrX61/G9VwOHtFs+GFizjX1WR0R/YCCwvsZjP2Hp0qVvR8QbXS1aPcLngLcbXYTqxv42L3vb3Oxvc7O/3e/QbW2oZwBfAgyPiGHAmxQXVZ7dYZ+5wGTgP4DxwMLMzIiYC8yKiBuBPweGA89t78My0yHwJhERrZnZ0ug6VB/2t3nZ2+Zmf5ub/a1W3QJ4ZrZFxGXAfKAfMCMzV0TEVKA1M+cC04H7ImIlxcj3hPLYFRHxIPAy0AZcmpmb61WrJEmSVJXI/NTUaqmh/C28udnf5mVvm5v9bW72t1o+CVM90V2NLkB1ZX+bl71tbva3udnfCjkCLkmSJFXIEXBJkiSpQgZwSZIkqUIGcFUqIkZHxCsRsTIiru5k+x4R8UC5/dmIGFquHxoRH0bEC+XXHVXXru2robejImJZRLRFxPgO2yZHxH+VX5Orq1q16mJ/N7c7d+dWV7VqVUN/p0TEyxGxPCIej4hD223z/O3Buthbz906cQ64KhMR/YD/BM6geNjSEmBiZr7cbp8fAMdm5sURMQH4ZmaeVQbxeZl5TPWVa0dq7O1Q4LPAlcDczHy4XL8/0Aq0UDzxdilwYma+W+G3oO3oSn/LbX/KzH2qrFm1q7G/pwHPZuaGiLgEOLX82ez524N1pbflNs/dOnEEXFUaCazMzNcy8yNgNjCuwz7jgJnl64eBr0ZEVFijds0Oe5uZ/52Zy4EtHY79GvDrzFxf/qP9a2B0FUWrZl3pr3q+Wvq7KDM3lIvPUDyhGjx/e7qu9FZ1ZABXlQ4Cft9ueXW5rtN9MrMNeA84oNw2LCKej4gnI+KUehernVJLb+txrKrR1R4NiIjWiHgmIv6me0tTN9jZ/p4P/HIXj1W1utJb8Nytm3o+il7qqLOR7I5zoLa1z1vAkMx8JyJOBOZExNGZ+T/dXaR2SS29rcexqkZXezQkM9dExGHAwoh4KTNf7aba1HU19zci/pZiuslf7uyxaoiu9BY8d+vGEXBVaTVwSLvlg4E129onIvoDA4H1mbkpM98ByMylwKvAEXWvWLWqpbf1OFbV6FKPMnNN+edrwBPAiO4sTl1WU38j4nTgGmBsZm7amWPVMF3preduHRnAVaUlwPCIGBYRuwMTgI5XVc8Ftl5FPx5YmJkZEYPLi0kofxMfDrxWUd3asVp6uy3zgTMjYlBEDALOLNep59jl/pZ93aN8/TngZODl7R+liu2wvxExAriTIqCtbbfJ87dn2+Xeeu7Wl1NQVJnMbIuIyyh+OPcDZmTmioiYCrRm5lxgOnBfRKwE1lP8sAAYBUyNiDZgM3BxZq6v/rtQZ2rpbUR8CXgEGAT8dURcl5lHZ+b6iPgxxT8UAFPtbc/Slf4CRwJ3RsQWikGff2l/BwY1Xo0/m38G7AM8VF4Xvyozx3r+9mxd6S2eu3XlbQglSZKkCjkFRZIkSaqQAVySJEmqkAFckiRJqpABXJIkSaqQAVySJEmqkAFcknqZiMiIuK/dcv+IWBcR8xpZ145ExBMR0dLoOiSp0QzgktT7fAAcExF7lstnAG82opDyibWSpJ1gAJek3umXwJjy9UTg37ZuiIi9I2JGRCyJiOcjYly5fmhEPB0Ry8qvvyjX/1lEPBURL0TEbyPilHL9n9q95/iIuKd8fU9E3BgRi4Drt/N5e0bE7IhYHhEPAFt/YZCkPs2RC0nqnWYD/1hOOzkWmAGcUm67BliYmedFxH7AcxGxAFgLnJGZGyNiOEVobwHOBuZn5rSI6AfsVcPnHwGcnpmbI+Kft/F5FwEbMvPYiDgWWNZt370k9WIGcEnqhTJzeUQMpRj9/vcOm88ExkbEleXyAGAIsAa4JSKOBzZThGgoHiM+IyJ2A+Zk5gs1lPBQZm7eweeNAm5uV+/ynfsuJak5GcAlqfeaC9wAnAoc0G59AN/OzFfa7xwR/wT8ETiOYgriRoDMfCoiRlFMabkvIn6WmfcC2e7wAR0++4MaPo8O7yFJwjngktSbzQCmZuZLHdbPBy6PMgFHxIhy/UDgrczcApwD9Cu3Hwqszcx/BaYDJ5T7/zEijoyIzwDf3E4d2/q8p4BJ5bpjKKbKSFKfZwCXpF4qM1dn5k2dbPoxsBuwPCJ+Wy4D3AZMjohnKKafbB3FPhV4ISKeB74NbH3Pq4F5wELgre2Usq3Pux3Yp5x68vfAczv9TUpSE4pM/3dQkiRJqooj4JIkSVKFDOCSJElShQzgkiRJUoUM4JIkSVKFDOCSJElShQzgkiRJUoUM4JIkSVKF/heR4T7ERtJ4oQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 743.226x572.185 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"(<Figure size 743.226x572.185 with 1 Axes>,\n",
" <matplotlib.axes._subplots.AxesSubplot at 0x278e4f1df48>)"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"estimator.parity_plot(\"Nd_aq_eq\", print_r_squared=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can also find what the r-squared value is. The closer to 1, the better the prediction model."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.9970803631106648\n"
]
}
],
"source": [
"print(estimator.r_squared())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Yay! Good job! That is an amazing fit."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.4"
}
},
"nbformat": 4,
"nbformat_minor": 4
}