{ "cells": [ { "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\".
\n", "```$ conda create --name thermo_env python=3.7```
\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
\n", "```$ git clone https://xgitlab.cels.anl.gov/summer-2020/parameter-estimation.git```
\n", "Navigate into the folder with
\n", "```$ cd parameter-estimation```
\n", "And run
\n", "```pip install -e.```
" ] }, { "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.
\n", "Let us get the pandas dataframe created by LLEPE." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
h_ih_eqz_iz_eqNd_aq_iNd_aq_eqNd_d_eqNd_org_eq
00.010.08830410.9216960.0500010.02391.09210.026101
10.010.10509410.9049060.0999980.06830.46410.031698
20.010.10901710.9009830.1500060.11700.28210.033006
30.010.10601210.9039880.2000040.16800.19050.032004
40.010.11893410.8910660.3000110.26370.13770.036311
\n", "
" ], "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.
\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.
\n", "Below is a table explaining the meaning of the column headers and the needed column order.
\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.
\n", "Please see parameter-estimation/data/xmls for file examples.
\n", "We can explore what has been loaded." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[, ]\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.
\n", "The number after the thermodynamic property is the initial guess for the optimizer.
\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,
\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.
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.
The default optimizer is from scipy.optimize.minimize with the arguments below. The optimizer optimizes a value multiplied by the initial guess.
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.
\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": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "data": { "text/plain": [ "(
,\n", " )" ] }, "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.7" } }, "nbformat": 4, "nbformat_minor": 4 }