{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Machine Learning with Python\n", "\n", "We now have our data. We have sanitized it into a csv format. We have explored it.\n", "\n", "Now lets try to predict some properties." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "## Load the data\n", "# df = pd.read_csv('../data/mpdata.csv')\n", "df = pd.read_csv('https://gitlab.com/costrouc/mse-machinelearning-notebooks/raw/master/data/mpdata.csv')" ] }, { "cell_type": "code", "execution_count": 3, "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", " \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", "
material_idenergyvolumensitesenergy_per_atompretty_formulaspacegroupband_gapdensitytotal_magnetizationpoisson_ratiobulk_modulus_voigtbulk_modulus_reussbulk_modulus_vrhshear_modulus_voigtshear_modulus_vrh
2621mp-559777-372.248991732.10601656-6.647303Na5Ca2Al(PO4)494.64962.730746-8.000000e-07NaNNaNNaNNaNNaNNaN
4063mp-603327-355.351207654.01937464-5.552363H3SNO3615.18851.9721481.312010e-02NaNNaNNaNNaNNaNNaN
5443mp-559382-17.05370733.4821243-5.684569CoO21640.00004.5097549.997688e-01NaNNaNNaNNaNNaNNaN
4616mp-558564-281.494573681.18521136-7.819294SiO2125.51131.757625-1.510000e-05NaNNaNNaNNaNNaNNaN
2883mp-667374-1188.8764402214.927434168-7.076645NaAlSiO41694.54142.5559691.502800e-03NaNNaNNaNNaNNaNNaN
\n", "
" ], "text/plain": [ " material_id energy volume nsites energy_per_atom \\\n", "2621 mp-559777 -372.248991 732.106016 56 -6.647303 \n", "4063 mp-603327 -355.351207 654.019374 64 -5.552363 \n", "5443 mp-559382 -17.053707 33.482124 3 -5.684569 \n", "4616 mp-558564 -281.494573 681.185211 36 -7.819294 \n", "2883 mp-667374 -1188.876440 2214.927434 168 -7.076645 \n", "\n", " pretty_formula spacegroup band_gap density total_magnetization \\\n", "2621 Na5Ca2Al(PO4)4 9 4.6496 2.730746 -8.000000e-07 \n", "4063 H3SNO3 61 5.1885 1.972148 1.312010e-02 \n", "5443 CoO2 164 0.0000 4.509754 9.997688e-01 \n", "4616 SiO2 12 5.5113 1.757625 -1.510000e-05 \n", "2883 NaAlSiO4 169 4.5414 2.555969 1.502800e-03 \n", "\n", " poisson_ratio bulk_modulus_voigt bulk_modulus_reuss bulk_modulus_vrh \\\n", "2621 NaN NaN NaN NaN \n", "4063 NaN NaN NaN NaN \n", "5443 NaN NaN NaN NaN \n", "4616 NaN NaN NaN NaN \n", "2883 NaN NaN NaN NaN \n", "\n", " shear_modulus_voigt shear_modulus_vrh \n", "2621 NaN NaN \n", "4063 NaN NaN \n", "5443 NaN NaN \n", "4616 NaN NaN \n", "2883 NaN NaN " ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df.sample(5)" ] }, { "cell_type": "code", "execution_count": 4, "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", " \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", " \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", " \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", "
energyvolumensitesenergy_per_atomspacegroupband_gapdensitytotal_magnetizationpoisson_ratiobulk_modulus_voigtbulk_modulus_reussbulk_modulus_vrhshear_modulus_voigtshear_modulus_vrh
energy1.000000-0.852020-0.9651980.3451710.176241-0.3787380.185483-0.1219210.045966-0.080688-0.059895-0.071393-0.128289-0.056726
volume-0.8520201.0000000.862495-0.110264-0.1623440.302136-0.3093430.0341270.000059-0.208067-0.187330-0.201073-0.110852-0.052474
nsites-0.9651980.8624951.000000-0.194668-0.2172660.371992-0.2482810.095329-0.032862-0.032160-0.028786-0.0309920.0423270.017493
energy_per_atom0.345171-0.110264-0.1946681.000000-0.038482-0.221015-0.334075-0.1715080.051704-0.565429-0.465262-0.523804-0.434698-0.211633
spacegroup0.176241-0.162344-0.217266-0.0384821.000000-0.0936780.250417-0.0655600.0284440.1956540.2059960.2044830.1138080.058407
band_gap-0.3787380.3021360.371992-0.221015-0.0936781.000000-0.421409-0.220998-0.068310-0.266141-0.249368-0.262229-0.0366940.000743
density0.185483-0.309343-0.248281-0.3340750.250417-0.4214091.0000000.3221210.0655130.5018740.5380170.5294850.1654600.086553
total_magnetization-0.1219210.0341270.095329-0.171508-0.065560-0.2209980.3221211.0000000.0230070.0801650.0895850.086458-0.006023-0.031027
poisson_ratio0.0459660.000059-0.0328620.0517040.028444-0.0683100.0655130.0230071.000000-0.019499-0.022243-0.021264-0.1434810.082642
bulk_modulus_voigt-0.080688-0.208067-0.032160-0.5654290.195654-0.2661410.5018740.080165-0.0194991.0000000.9304990.9819580.6760280.325090
bulk_modulus_reuss-0.059895-0.187330-0.028786-0.4652620.205996-0.2493680.5380170.089585-0.0222430.9304991.0000000.9829770.5982730.312095
bulk_modulus_vrh-0.071393-0.201073-0.030992-0.5238040.204483-0.2622290.5294850.086458-0.0212640.9819580.9829771.0000000.6479450.324180
shear_modulus_voigt-0.128289-0.1108520.042327-0.4346980.113808-0.0366940.165460-0.006023-0.1434810.6760280.5982730.6479451.0000000.460951
shear_modulus_vrh-0.056726-0.0524740.017493-0.2116330.0584070.0007430.086553-0.0310270.0826420.3250900.3120950.3241800.4609511.000000
\n", "
" ], "text/plain": [ " energy volume nsites energy_per_atom \\\n", "energy 1.000000 -0.852020 -0.965198 0.345171 \n", "volume -0.852020 1.000000 0.862495 -0.110264 \n", "nsites -0.965198 0.862495 1.000000 -0.194668 \n", "energy_per_atom 0.345171 -0.110264 -0.194668 1.000000 \n", "spacegroup 0.176241 -0.162344 -0.217266 -0.038482 \n", "band_gap -0.378738 0.302136 0.371992 -0.221015 \n", "density 0.185483 -0.309343 -0.248281 -0.334075 \n", "total_magnetization -0.121921 0.034127 0.095329 -0.171508 \n", "poisson_ratio 0.045966 0.000059 -0.032862 0.051704 \n", "bulk_modulus_voigt -0.080688 -0.208067 -0.032160 -0.565429 \n", "bulk_modulus_reuss -0.059895 -0.187330 -0.028786 -0.465262 \n", "bulk_modulus_vrh -0.071393 -0.201073 -0.030992 -0.523804 \n", "shear_modulus_voigt -0.128289 -0.110852 0.042327 -0.434698 \n", "shear_modulus_vrh -0.056726 -0.052474 0.017493 -0.211633 \n", "\n", " spacegroup band_gap density total_magnetization \\\n", "energy 0.176241 -0.378738 0.185483 -0.121921 \n", "volume -0.162344 0.302136 -0.309343 0.034127 \n", "nsites -0.217266 0.371992 -0.248281 0.095329 \n", "energy_per_atom -0.038482 -0.221015 -0.334075 -0.171508 \n", "spacegroup 1.000000 -0.093678 0.250417 -0.065560 \n", "band_gap -0.093678 1.000000 -0.421409 -0.220998 \n", "density 0.250417 -0.421409 1.000000 0.322121 \n", "total_magnetization -0.065560 -0.220998 0.322121 1.000000 \n", "poisson_ratio 0.028444 -0.068310 0.065513 0.023007 \n", "bulk_modulus_voigt 0.195654 -0.266141 0.501874 0.080165 \n", "bulk_modulus_reuss 0.205996 -0.249368 0.538017 0.089585 \n", "bulk_modulus_vrh 0.204483 -0.262229 0.529485 0.086458 \n", "shear_modulus_voigt 0.113808 -0.036694 0.165460 -0.006023 \n", "shear_modulus_vrh 0.058407 0.000743 0.086553 -0.031027 \n", "\n", " poisson_ratio bulk_modulus_voigt bulk_modulus_reuss \\\n", "energy 0.045966 -0.080688 -0.059895 \n", "volume 0.000059 -0.208067 -0.187330 \n", "nsites -0.032862 -0.032160 -0.028786 \n", "energy_per_atom 0.051704 -0.565429 -0.465262 \n", "spacegroup 0.028444 0.195654 0.205996 \n", "band_gap -0.068310 -0.266141 -0.249368 \n", "density 0.065513 0.501874 0.538017 \n", "total_magnetization 0.023007 0.080165 0.089585 \n", "poisson_ratio 1.000000 -0.019499 -0.022243 \n", "bulk_modulus_voigt -0.019499 1.000000 0.930499 \n", "bulk_modulus_reuss -0.022243 0.930499 1.000000 \n", "bulk_modulus_vrh -0.021264 0.981958 0.982977 \n", "shear_modulus_voigt -0.143481 0.676028 0.598273 \n", "shear_modulus_vrh 0.082642 0.325090 0.312095 \n", "\n", " bulk_modulus_vrh shear_modulus_voigt shear_modulus_vrh \n", "energy -0.071393 -0.128289 -0.056726 \n", "volume -0.201073 -0.110852 -0.052474 \n", "nsites -0.030992 0.042327 0.017493 \n", "energy_per_atom -0.523804 -0.434698 -0.211633 \n", "spacegroup 0.204483 0.113808 0.058407 \n", "band_gap -0.262229 -0.036694 0.000743 \n", "density 0.529485 0.165460 0.086553 \n", "total_magnetization 0.086458 -0.006023 -0.031027 \n", "poisson_ratio -0.021264 -0.143481 0.082642 \n", "bulk_modulus_voigt 0.981958 0.676028 0.325090 \n", "bulk_modulus_reuss 0.982977 0.598273 0.312095 \n", "bulk_modulus_vrh 1.000000 0.647945 0.324180 \n", "shear_modulus_voigt 0.647945 1.000000 0.460951 \n", "shear_modulus_vrh 0.324180 0.460951 1.000000 " ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df.corr()" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAARAAAAD1CAYAAACY9KqQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGdpJREFUeJzt3XuQHeV55/Hvj5EYSVwkgbQsSDKChRhzCzgq28FVtoMwJg6FXbWsLSrEEOOlnHVshzi+aKlatsi6iqxTwXHVLrYWcwsUtlfGCRvbgMIlyZYtYmRY7hHiasmyhSyQjJA0mpln/+ge+zCamdNvd5/WnNbvU9Wlc+k+b/dozjNvv7dHEYGZWRkH7e8TMLP+5QBiZqU5gJhZaQ4gZlaaA4iZleYAYmalOYCY9RFJN0jaIunxSd6XpK9I2iDpUUlv7XjvEknP5NsldZyPA4hZf7kJOG+K938XODHfLgeuA5B0BHAV8HbgbcBVkuZXPRkHELM+EhH/BGybYpcPALdEZi0wT9LRwPuANRGxLSJeAdYwdSAqxAHErF0WAT/peL4xf22y1yuZUfUDzGxq7/udQ+IX20YK7bvu0T1PALs7XloVEat6cmI1cAAx67Gt20Z48O7FhfadefSzuyNiWYXiNgFLOp4vzl/bBLxn3OsPVCgH8C2MWQOCkRgttNXgTuAjeW/MO4DtEbEZuBs4V9L8vPH03Py1SlwDMeuxAEapZ9a7pNvJahILJG0k61mZCRARXwW+B7wf2AC8Dvxh/t42SX8O/Cj/qKsjYqrG2EIcQMwaMEottQsi4qIu7wfwiUneuwG4oZYTyTmAmPVYEIy0dN2d/d4GIuk8Sf+aj5z7Qg8+f4mk+yU9KekJSZ+uu4xx5Q1IeljS3/fo8+dJWi3paUlPSfrtHpVzRf7zelzS7ZJm1fS5+4yklHSEpDX5CMk1VQc4TVLGl/Kf2aOSviNpXpUyUo0ShbZ+s18DiKQB4H+QjZ47GbhI0sk1FzMMfCYiTgbeAXyiB2V0+jTwVA8//6+BuyLiJOA3e1GWpEXAp4BlEXEqMACsqOnjb2LfAUxfAO6NiBOBe/PndZexBjg1Ik4H1gMrK5ZRWAAjRKGt3+zvGsjbgA0R8VxEDAHfIBtJV5uI2BwRP84f/5LsC1d5AM1EJC0Gfg+4vkefPxd4F/B1gIgYiohXe1EW2e3tbEkzgDnAT+v40ElGUn4AuDl/fDPwwbrLiIh7ImI4f7qWrBuzEQHsjdFCW7/Z3wGkJ6PjJiNpKXAm8GCPivgy8DmoqcVsX8cBLwM35rdJ10s6pO5CImIT8JfAS8Bmsq7Ae+oup8NReVcjwM+Ao3pYFsBHge/3uIw3GC249Zv9HUAaI+lQ4NvAn0TEjh58/vnAlohYV/dnd5gBvBW4LiLOBHZSvbq/j7wN4gNkAesY4BBJF9ddzkTyXoSe1eUlXUl2W3tbr8oYLwrevvgWJt1ko+ZqJWkmWfC4LSLuqPvzc+8ELpD0Atmt2NmSbq25jI3AxogYq0GtJgsodTsHeD4iXo6IvcAdwFk9KGfMz/MJX+T/bulFIZIuBc4Hfj+aTEcQMFJw6zf7O4D8CDhR0nGSDiZrqLuzzgIkiazN4KmI+Ks6P7tTRKyMiMURsZTsOu6LiFr/akfEz4CfSHpz/tJy4Mk6y8i9BLxD0pz857ec3jYM3wmMrU9xCfB3dRcg6Tyy28sLIuL1uj9/KtlAsnbewuzXcSARMSzpj8mG1A4AN0TEEzUX807gD4DHJD2Sv/afI+J7NZfTlE8Ct+UB9znykYZ1iogHJa0GfkxW3X8YqGVC1yQjKa8BviXpMuBF4EM9KGMlMAisyWIiayPi41XKSTgjRlAzRTVMTixl1lunnn5wfPu7Cwrte9KbNq+rOJmuUR6JatZjAQzt99aC3nAAMWvAaLTzFsYBxKzHspGoDiBmVkIgRlp6CzMtrkrS5W0px9cyPctp6lomMxoqtPWbaRFAyJafb0s5vpbpWc5+CyBjtzBFtn7jWxiznhMjMV3+Vter0QCy4IiBWLpk5j6vv2nRDJb95qwJB6Q88+RhyeXE8MQrYM9iDofriH3KGTlxMLmM4eGBCV8fOHIeg8cvnvBaZuxI/wszfPjE4xMHjpzL4PGL9ilHu9N/UWOS34KB+fMYfNOSCa9FxRYZH3fQxC/PmDefWYsnLid5esgk+8+YO59Zx0xSRonv9p5NG7dGxMKip7SXiX9f+l2jAWTpkpn8y91Luu/Y4f2nL08uZ2TrL5L23/6VE5LL2LL18ORjFq5JD1Rbz93dfacOg+tnJ5exe2F6NJi5vcS3rswf4cQAMrCnRJCenT6Y8tmVn3mx6L4RroGYWQWjfdi+UUSlsNjr5QjN2iBrRD2o0NZNt++cpGslPZJv6yW92vHeSMd7tUxaLV0D6ViO8L1k08x/JOnOiOjF7FCzPlbPLUyR71xEXNGx/yfJFtAasysizqh8Ih2qXFXPlyM0a4NsOv9BhbYuUr9zFwG313MVE6sSQBpdjtCsXwViKAYKbWRLEDzUsXWOXyn8nZN0LNmKcvd1vDwr/8y1kiqtOzum542o+Q/gcsi6a80ORKPFb2G21jSdfwWwOiI6u9iOjYhNko4H7pP0WEQ8W6WQKjWQQssRRsSqiFgWEcsWHtnOvnCzqdTYiJqyBOgKxt2+5ItlExHPkSXWPnPfw9JUCSA9X47QrA0CMRLFti4KfecknQTMB37Y8dp8SYP54wVkK/VV7vAofU/R0HKEZq1QoIG0q8m+c5KuBh6KiLFgsgL4xriFo98CfE3SKFnF4Zo6ekwrNUrk64r269qiZo2IoLaRqBN95yLiv4x7/l8nOO4HwGm1nEQHt2qa9ZxaOxK10QDyzJOHJc9t+d6j9yaX88TQrqT9L7wxPeXJ7KHkQxi6cHxGx+7mfTctz/Qrp6QnBzhoKP2Xe/DV9GOiTBt64uUM156nr7oAhiabsdjn2nlVZtNI0J+LBRXhAGLWgLYuaegAYtZjQdJAsr7iAGLWc/25XGERDiBmPeYaiJlV4hqImZUSIfaOtvOr1s6rMptGsvVAXAMxs1K8qLKZlZQ1oroGYmYleSCZmZXioew1ieGR5KRPqRPjAE45OC250ozXk4solZntsFl7ko/ZMSvtFy8G009Mu9Jnue1amJ6MqQkjc9InE2qk91/uOtYDmY5cAzHrsQjYO9rOAFL6qiQtkXS/pCclPSHp03WemFlbZLcwBxXa+k2VGsgw8JmI+LGkw4B1ktY4sZTZvjwSdZyI2Axszh//UtJTZDkqHEDMOrS5G7eWOpOkpWRLxD9Yx+eZtUt9tzAFcuNeKunljhy4H+t47xJJz+TbJXVcWeVGVEmHAt8G/iQidkzw/q8SS81iTtXizPpSHUPZE/JRfzMi/njcsUcAVwHLyCpF6/JjX6lyTpVqIJJmkgWP2yLijon26UwsNZPBKsWZ9aVsVfZa8sJUyUf9PmBNRGzLg8Ya4LzSF5Wr0gsj4OvAUxHxV1VPxKytAjE8OlBo66Jobtx/L+lRSasljWWy60ku6yo1kHcCfwCc3XG/9f6qJ2TWRqN5aoduG1Mn1y7i/wBLI+J0slrGzXVfS6cqvTD/F1raN2VWo8RemKmSa3fNjRsRnUO9rwf+e8ex7xl37ANFT2oy/TdyxawP1dQL0zU3rqSjO55eADyVP74bODfPkTsfODd/rZJGh7KPnDjI9q+ckHRMmaRPqXNbHrvifyaX8Ru3/FHyMdvvOrr7TuPM3Jk25+SQF9L/S19flD5/5ri/HU4+Zmhu+rkN7Eqb27L3sPR5Pb9c0uO/o1HPZLqCuXE/JekCsoGe24BL82O3SfpzsiAEcHVEpGc6G8dzYcx6rM4Vybrlxo2IlcDKSY69AbihlhPJOYCYNaCtI1EdQMx6LIDhls7GdQAx6zEvKGRmlXhVdjMrJ9wGYmYltXk6vwOIWQMcQMyslECMuBfGzMpyI6qZlRJuRDWzKsIBpLrh4QG2bD086ZjZQ+nlpCZ9KjMxbv1Hrks+5qT/9Z+Sjxk667Wk/UfXH5pcRpnESi9+ND2BE0pPrJVqJO3HBYCGep0kywPJzKwC10DMrBSPA5lCvlL0Q8CmiDi/+imZtUy+qHIb1VED+TTZqkdpjRtmB4igvbcwVdM6LAZ+j2ztRTObUNaIWmTrN1VrIF8GPgccVsO5mLVW9LqjZz+pkhfmfGBLRKzrst/lY0vUj+zYWbY4s74WoUJbv6maF+YCSS+QZcg6W9Kt43fqzEw3cPghFYoz608R9QWQArlx/1TSk3liqXslHdvx3khHDqc7xx9bRukAEhErI2JxRCwlW17+voi4uI6TMmubOtpAOnLj/i5wMnCRpJPH7fYwsCxPLLWaX+eFAdgVEWfk2wV1XFc7pwiaTTOjoyq0ddE1N25E3B8RY4lN1pIlkOqZWgJIRDzgMSBmEwuK3b4UuIVJzW97GfD9juez8vbItZI+WO5q3sgjUc0akNAJs0DSQx3PV0XEqtTyJF0MLAPe3fHysRGxSdLxwH2SHouIZ1M/u1OjAWTGDrFwzWDSMUMXpifPOmxW2qStMhnjykyMe/o/pmfAO+uKjyftv/OY9Jb83SW6GOdsmJ18TPfMjftS4py9nYvTL2Z0Ro/7WCNpIFml3LgAks4BrgTeHRG/+jJExKb83+ckPQCcCVQKIG4DMWtCFNymViQ37pnA14ALImJLx+vzJQ3mjxeQ9aI+WfWyfAtj1oA6xngUzI37JeBQ4H9LAngp73F5C/A1SaNkFYdrIsIBxKwf1DUStUBu3HMmOe4HwGn1nMWvOYCY9VgEhBdVNrOy2joXxgHErAkOIGZWTn9OlCvCAcSsCa6BmFkpaQPJ+ooDiFkTXAMxs9JcAzGz0lwDqW748FG2nrs76Zh5352fXM6OWWnRfubO9P/d1IxxkD4xDuAH1341af/fWveh5DL2vpT+Mz5t+VPJx8we2Jt8zGjiDLx7n35zchmxq8dfg8A1EDMrr60DyaqmdZgnabWkpyU9Jem36zoxs1apZzbutFO1BvLXwF0RcWE+vXhODedk1j6+hXkjSXOBdwGXAuRrNA7Vc1pmLRLpCyP1iyq3MMcBLwM3SnpY0vWSnLfBbB/KaiBFtj5TJYDMAN4KXBcRZwI7gYnyVDixlFlL20CqBJCNwMaIeDB/vposoLyBE0uZ4QAyXkT8DPiJpLGO9+XUsMaiWSu1NIBU7YX5JHBb3gPzHPCH1U/JrGU8kGxiEfEIWe4JM5uC+rB2UUQ7F2o0m25quoUpkFx7UNI38/cflLS0472V+ev/Kul9NVxVs0PZtfsgBtenJSR65ZT0DvQYHEna/5AX0n8Mo+sPTT6mTNKn1Lkt637rW8ll/LsN6XN0/vmRk5KPoUwtPvEv98BrA8lFKO3XpZQ6aiAdybXfS9aJ8SNJd45Lz3AZ8EpEnCBpBfAXwIfzJNwrgFOAY4B/kPQbEVHp6l0DMWtCPeNAuibXzp/fnD9eDSxXliDmA8A3ImJPRDwPbMg/rxIHELNeK3r70r2WUiS59q/2iYhhYDtwZMFjk3k2rlkTit/C1JJcuykOIGYNSGgDqZpce2yfjZJmAHOBXxQ8NplvYcya0FBy7fz5JfnjC4H7IiLy11fkvTTHAScC/1LxqlwDMes11TQbt2By7a8DfyNpA7CNLMiQ7/ctstHiw8AnqvbAgAOIWTNqGolaILn2buA/THLsF4Ev1nIiOQcQsya0dCSqA4hZA9o6lN0BxKwJDiBmVkq4BmJmVTiAVBczYPfCtJ6jg4bSW6+1K21C1euL0nuzNJJ+XrtL/BKlJn0qMzHu2Q+nJa8COPHWP0o+pokv0cwd6f8vuxYP9+BM3siLKpuZjeNbGLMmtPQWpmpmuiskPSHpcUm3S5pV14mZtUbeiFpk6zelA4ikRcCngGURcSrZ0NoVdZ2YWat4UeVJj58taS9ZWsufVj8lsxbqw+BQRJW0DpuAvwReAjYD2yPinrpOzKwthG9h9iFpPtkyaceRrbF4iKSLJ9jv15npXnut/Jma9at8Nm6Rrd9UaUQ9B3g+Il6OiL3AHcBZ43d6Q2a6Q9MXIjZrhZa2gVQJIC8B75A0J1+0dTnwVD2nZdYyLQ0gpRtRI+JBSauBH5MtUPIwMG3XbjTbn/qxfaOIqpnprgKuqulczNrLAcTMSunT25Mims1MNwIzt6c1uwy+WmJy1MK0/63j/jZ9MtWLH01vMp+zIS0rH8Bpy9OalcpkjCszMe6Zi69LPmbLyM7kY+YflDa4+auvHp9cxi3Pvz35mJcS9+/HHpYiPJnOrAFNjAORdISkNZKeyf/dZyq3pDMk/TCfgvKopA93vHeTpOclPZJvZ3Qr0wHErAnN9MJ8Abg3Ik4E7s2fj/c68JGIOAU4D/iypHkd7382Is7It0e6FegAYtZr9aW27KYzL+7NwAf3OZWI9RHxTP74p8AWYGHZAh1AzHpMCVtFR0XE5vzxz4Cjpjwv6W3AwcCzHS9/Mb+1uVbSYLcC3Qtj1oSacuNK+gfg305w3JVvKC4ipMlbVSQdDfwNcElEjDXxriQLPAeTjen6PHD1VCfrAGLWgJpy4xIR50xahvRzSUdHxOY8QGyZZL/Dge8CV0bE2o7PHqu97JF0I/Bn3U7WtzBmTRgtuFXTmRf3EuDvxu+Q59T9DnBLRKwe997R+b8iaz95vFuBDiBmvdbcimTXAO+V9AzZZNdrACQtk3R9vs+HgHcBl07QXXubpMeAx4AFwH/rVqBvYcya0MBI1Ij4Bdmk1vGvPwR8LH98K3DrJMefnVqmA4hZAzyZzszKcwAxs7JcA6mDSG62jbQkc6UMzS3xY9Ce5EOiRJP17IG9aQeUGY1U4pe7zMS4fzNwSPIxeyMta+DmobnJZYyM1jCEayqejWtmZYkDeDaupBskbZH0eMdrXWf9mVmHli5pWKRSfRPZrL1ORWb9mVlOEYW2ftM1gETEPwHbxr3cddafmeWam43buLJtIEmz/swOdO6FmUSBWX+XA5cDzJjnphI7QLU0gJSdC/Pzjok3k876g3GJpQ5J78YzawOntnyjrrP+zCx3IKe2lHQ78EPgzZI2SrqMSWb9mdkkDtRG1Ii4aJK39pn1Z2b7Ev15e1KER6KaNaEPx3gU4QBi1gDXQOpQ5j6vgYalgV3NtF6VaSQbTZ2B19AvamrGOEifGAcwU2mzKefPTJ/kN8UohHr0aftGEa6BmDWgH3tYinAAMWuAA4iZlRO0thHVq7KbNWC6JNfO9xvpWJH9zo7Xj5P0oKQNkr6Zp4CYkgOIWROmT3JtgF0dCbQv6Hj9L4BrI+IE4BXgsm4FOoCY9djYQLIG5sKUXmYjTyZ1NjCWbKrQ8Q4gZr0WUXzLc+N2bJcnlFR0mY1Z+WevlTQWJI4EXo2I4fz5RmBRtwLdiGrWgIRemClz49aUXPvYiNgk6Xjgvjwb3fbCZ9jBAcSsAXWNVasjuXZEbMr/fU7SA8CZwLeBeZJm5LWQxcCmbufjWxizXgtgNIpt1RRJrj1f0mD+eAHwTuDJiAjgfuDCqY4fzwHErAnN9MIUSa79FuAhSf+PLGBcExFP5u99HvhTSRvI2kS+3q3AxufCDOxJS+IzXGIRs5E5acP+9h6Wnr1q5LXkQ9i5OP035N6n35y0/8Br6dcyc0d6YqWvvnp88jFlkj6lzm357BHPJpfxw23p1/JI4v5NTKYrmFz7B8Bpkxz/HPC2lDLdBmLWhJaORHUAMWtAW6fzl81M9yVJT0t6VNJ3JM3r7Wma9S8FaDQKbf2mbGa6NcCpEXE6sB5YWfN5mbXLaMGtz5TKTBcR93SMWFtL1mdsZpM4YFNbFvBR4Ps1fI5ZOzm15cQkXQkMA7dNsc+vM9PNdWY6OxCFe2HGk3QpcD6wPB/FNqGIWAWsAph1zJJ2/hTNumhrL0ypACLpPOBzwLsj4vV6T8mshQ7UGkieme49ZNOMNwJXkfW6DAJrsmUEWBsRH+/heZr1rwCNHKABZJLMdF3HyJtZh3bGD49ENWtCP3bRFtFsADkIhmf3/gepkbTJYb9ckt6braH06xidkX5M7Er7L1J67iZ2LR7uvtM4tzz/9uRjRkbTJ+2lJn0qMzHujhPWJB+TPGXRAcTMSgn6cpRpEQ4gZj0m+nOUaREOIGZNcAAxs1ICOFC7cc2sOt/CmFl5LQ0gXlTZrOeSEkuVViQ3rqTf6ciL+4ik3WPJpSTdJOn5jvfO6FamA4hZrwWNBBAK5MaNiPvH8uKSpbJ8HbinY5fPduTN7bp2tAOIWROaWZEsNTfuhcD3q0yIdQAxa0BDK5IVzY07ZgVw+7jXvpivdXztWAKqqbgR1azXAhgpXL1YIOmhjuer8jV1gNpy45KnvjwNuLvj5ZVkgedgsjV8Pg9cPdXJOoCY9VxS+8aUybXryI2b+xDwnYjY2/HZY7WXPZJuBP6s28k2GkD2bNq49dmVn3lxgrcWAFsbOIUmyjkgruWlhspJNUWr36RlpOfyA+DYpL2b6cYdy417Dd1z217EuGwKHcFHZO0nj094ZIdGA0hELJzodUkPTRV169JEOb6W6VlOU9cyqWYCyDXAtyRdBrxIVstA0jLg4xHxsfz5UmAJ8I/jjr9N0kJAZLG46yJhvoUx67UAGkgaVSQ3bv78BWDRBPudnVqmA4hZzwVEO+fzT5cAsqr7Ln1Tjq9lepbT1LXsK60Xpq9oiowMZlaDuQcfFWcdtaLQvndt/Mq6/dpWk2i61EDM2q2lf6gdQMx6zpnpzKysAEbb2QbiAGLWBNdAzKw0BxAzKyWCGCmRsKcPOICYNaGBkaj7gwOIWRN8C2NmpUS4F8bMKnANxMzKCtdAzKwcj0Q1s7ICcDeumZURQLgb18xKCS8oZGYVtLUG4gWFzHpM0l1kq8IXsTUizuvl+dTJAcTMSnNqSzMrzQHEzEpzADGz0hxAzKw0BxAzK80BxMxKcwAxs9IcQMysNAcQMyvt/wNLeDG7G2ksYwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.matshow(df.corr())\n", "plt.colorbar()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Lets choose a very simple example to show methodology\n", "\n", "How about we try to predict the `energy_per_atom`. You can see from the correlation plot that there are two very highly correlated values in purple.\n", "\n", "We will simplify our model and only use the first four columns. Obviously `volume` is not usefull in the calculation but we want to see if our algorithm can automatically determine this." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "simplified_df = df[['energy', 'volume', 'nsites', 'energy_per_atom']]" ] }, { "cell_type": "code", "execution_count": 7, "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", "
energyvolumensitesenergy_per_atom
0-4.06460011.8527651-4.064600
1-16.38209647.2641584-4.095524
2-8.18695923.6173882-4.093479
3-4.06414211.8747031-4.064142
4-2.157191603.4752101-2.157191
\n", "
" ], "text/plain": [ " energy volume nsites energy_per_atom\n", "0 -4.064600 11.852765 1 -4.064600\n", "1 -16.382096 47.264158 4 -4.095524\n", "2 -8.186959 23.617388 2 -4.093479\n", "3 -4.064142 11.874703 1 -4.064142\n", "4 -2.157191 603.475210 1 -2.157191" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "simplified_df.head(5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "All (99%) of machine learning algoritms need the data as arrays of floating point numbers. Scikit learn is no different. This is how easy it is to convery a pandas dataframe from a numpy array.\n", "\n", "Not covered here but you most likely will need it at one point [preprocessing data](http://scikit-learn.org/stable/modules/preprocessing.html) and how to handle categorical data." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(6928, 3) (6928,)\n", "[[ -4.0645998 11.85276501 1. ]\n", " [-16.38209642 47.26415795 4. ]\n", " [ -8.18695876 23.61738783 2. ]] [-4.0645998 -4.0955241 -4.09347938]\n" ] } ], "source": [ "# convert from pandas dataframe to numpy array\n", "X = simplified_df[['energy', 'volume', 'nsites']].values\n", "y = simplified_df['energy_per_atom'].values\n", "\n", "print(X.shape, y.shape)\n", "print(X[:3], y[:3])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Scikit Learn\n", "\n", "Very quick overview. [Scikit](http://scikit-learn.org/stable/) learn provides a unified framework for working with machine learning algorithms. It includes classification, regression, clustering, dimensionality reduction, model tuning, pre and post processing of data.\n", "\n", "Is that a lot? **YES** scikit learn is huge and you cannot expect to use and learn everything.\n", "\n", "The flow chart gives some good advice for which algorithms to use for your problem. See their [flow chart](http://scikit-learn.org/stable/tutorial/machine_learning_map/index.html)\n", "\n", "![sklearn flowchart](../images/scklearn-flowchart.png)\n", "\n", "There are a ton of algoritms over 100! This is where sklearn really shines. All algorithms have the exact same api (this is the pseudocode).\n", "\n", "```python\n", "from sklearn import MyImportantModel\n", "\n", "model = MyImportantModel()\n", "model.fit(X, y)\n", "```\n", "\n", "Once you have `fit` your model you can using is to predict future data.\n", "\n", "```python\n", "y_predict = model.predict(X_predict)\n", "```\n", "\n", "We will be using a linear model to fit our data. Always start with the simplest model! That way you know what sort of improvement a complex one can get you.\n", "\n", "[sklearn.LinearRegression](http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "# Lets use a simple linear model\n", "from sklearn.linear_model import LinearRegression\n", "\n", "model = LinearRegression()" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEKCAYAAAAMzhLIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXl4U1X6x7/n3tzs6ZKkLZQuuACKCio4I26AO0uRxRlQcN8AcRRZnN+IigrMCFRk3NDRGVGK2wjIVkSURRlGR1xAK+BGFwpt03TLntx7fn/c5Da3SdoCLW3hfJ6nT9vk3txzT5LznnPe9/2+hFIKBoPBYDASwXV0AxgMBoPRuWGGgsFgMBjNwgwFg8FgMJqFGQoGg8FgNAszFAwGg8FoFmYoGAwGg9EszFAwGAwGo1mYoWAwGAxGszBDwWAwGIxm0XR0A9oCu91Oe/bs2dHNYDAYjC7F7t27HZTStJaO67SGghByEEADABFAiFI6MNGxPXv2xFdffXWimsZgMBgnBYSQ4tYc12kNRZihlFJHRzeCwWAwTmWYj4LBYDAYzdKZDQUFsJkQspsQcm9HN4bBYDBOVTrz1tNllNJDhJB0AB8TQvZRSndEngwbj3sBICcnp6PayGAwGCc9nXZFQSk9FP5dCWA1gN81ef5VSulASunAtLQWnfYMBoPBOEY6paEghJgIIZbI3wCuBfB9x7aKwWB0RgoKVuK0XmeB43mc1ussFBSs7OgmnXR01q2nDACrCSGA3MaVlNJNHdskBoPR2SgoWIkp02fBeNU0ZI/uC29ZEaZMnwUAmDjx5g5u3ckDORlKoQ4cOJCyPAoG49TjtF5nwTvwNuhz+ymP+Yr3wPDVcvz2074ObFnXgBCyu7kctQidcuuJwWAwWkPxrz9Bl9VX9Zguqy+Kf/2pg1p0csIMBYNxCtPV9/dzT+8Ff1mR6jF/WRFyT+/VQS06OWGGgsE4RYns73sH3obsh1fBO/A2TJk+q0sZi3lzH4fnkxfgK94DKobgK94DzycvYN7cxzu6aScVzFAwGF2MtloFzJn7FIxXTYM+tx8Ir4E+tx+MV03DnLlPtXGLE3O89zJx4s14eckiGL5ajtJnx8Lw1XK8vGQRc2S3MZ016onBYMShLaN8in/9Cdmj4+zvv3di9vfb6l4mTryZGYZ2hkU9MRhdiLaM8unoiKGOvj6DRT0xGCclbRnl09H7+yxiqevAtp4YjC5E7um94C0rUs3CjzXKJ7JdM2fuUyh+7yfknt4L+Sdwf78t74XRvrAVBYPRhWjrVcDEiTfjt5/2QRJF/PbTvhO619/RKxpG62ErCgajC9HRq4C25GS6l5Md5sxmMBiMUxTmzGYwGAxGm8AMBYPBYDCahRkKBoPBYDQLMxQMBoPBaBZmKBgMBoPRLMxQMBgMBqNZmKFgMBgnhK5e++JUhhkKBqMDOVUGz5Oh9sWpDEu4YzA6iGiZbV1WX/jLiuD55IWTsp4CU4rtnLQ24Y4ZCgajgziVBk+O55H98CoQvlE1iIohlD47FpIodmDLTm1YZjaD0ck5lWS2WW3rrg0zFAxGB3E8g2dX820wpdiuDVOPZTA6iHlzH5dLfzbxUeQvWdTseW1ZDvVEwZRiuzbMR8FgdCAFBSvlwfNXefCcN/fxFgfPU8m3wWhfmDObwThJYY5hRlvBnNkMxkkKcwwzTjTMUDAYXYxox7Dr+09x6JV7UPHOo3C73R3u1O5qTnZG62DObAajixHxYTw0YxZqGjxIGzVbcYZ3pFO7KzrZGa2D+SgYjC5KZ3Nqd7b2MFqG+SgYjJOYgoKVOPjzgTZL2GuLLaNTKYHwVIMZCgajixHZ4tGkdGsTp3ZbCfa1h5Od+Tw6B53WUBBCrieE7CeE/EwI+XNHt4fB6CzMmfsUjFdNQ8rlk1BduPS4s50jr6fP7QfCa6DP7QfjVdMwZ+5TR/U6bZ19zRRnOw+d0lAQQngALwIYBqAvgJsIIX2bP4vB6Ny01ew4ssVj6jsYKVfcCueWZSjJH4vKD54+JuXZttoymjjxZry8ZBEMXy1H6bNjYfhqedz2tLYf2sqAnWwEg0HU19ef0Gt21qin3wH4mVL6KwAQQt4BcAOAombPYjA6KW0ZEZR7ei94y4qgz+0HU9/BMPUdrDiNjyW6KPr1IhzrltHEiTc324aj6YfiX39C9ug4Buy9U8vnEQqF8M0332Dr1q3YunUrPv/8c0ydOhXPPPPMCWtDp4x6IoTcCOB6Sund4f9vAfB7Sum0qGPuBXAvAOTk5AwoLi7ukLYyGK2hLSOCpt4/Da/+8w2Ifg8EazYMvQeB/vzZMdexOJF1MY6mH07VKCpJkvDdd98phmHHjh0xK4iLLroIX3755XFfq7VRT511RdEilNJXAbwKyOGxHdwcBqNZ2mp2XFCwEiveXw376EeVQb16/WLcc+tNxzyon0jBvqPph2MVTexqSJKEH374QTEM27dvR01NTbPn7N69G3V1dUhOTj4hbeyshuIQgOyo/7PCjzEYXZK22t6J3rcHAH1uP9hGzkTh5uXH1b6WtozaiqPph1NBcfaOO+7A+vXr4XA4juo8u92OX375BRdeeGE7tUxNp3RmA/gfgF6EkNMIIVoAEwCs7eA2MU5B2soB3VYRQV09V+Fo+2HixJvx20/7IIkifvtp30llJACgvLy8VUbCZrNh3LhxeOGFF/DDDz/gyJEjJ8xIAJ10RUEpDRFCpgH4CAAP4J+U0h86uFmMU4y2dEC31ey4LR3PHcGpsEqglOK3335TtpL27t2Lb775BhwXOy8fOnQoNm/eHPN4SkoKBg8ejKFDh2Lo0KE499xz455/ouiUzuyjhUl4MNqDjnamxqtVAQB3T5mGgEhBg37QoA+czoD77roDL734Qru3iRGfkpISxTBs3boVJSUlque//fZb9O/fP+a8L774AhdffDGSkpJwxRVXKIahX79+4Hm+3dt90juzGYz2piPDMxOtZib9YQxESQLRaFVigG+8vRSXXnLJSTUz78yUl5erDMOvv/7a7PFbt26NaygGDBiAL7/8EhdccAE0ms47HHdWHwWD0eG0RpKivSQmEiWbvfbGm6CCEfbh01XPJV334CmfiNbe7N69G5MnT0afPn3Qo0cPTJo0Ca+//nqLRgIAPvvss7iPazQaXHTRRZ3aSADMUDAYCWnJ8doWEhOJDE0ip3XQ40KorqJLO7S7KgcPHsQrr7yCAwcOtHisTqfDkCFD8OSTT2LHjh1YubJry450bjPGYHQgLTle44WqIiwx0ZotoOac5Ymc1pzOAPBa+LuwQ7szUlNTgx07dmDr1q0YOXIkrr766phjBg8enPB8QRDw+9//XvExDBo0CHq9vj2bfEJhzmwG4xg53trVzTnLI8lm0dnSjsKlMPUdAveezZAkEek3/Fl5rv6jpfjH8892WMGipk73zu4rqa+vVwzDtm3b8M033yAyFk6dOhUvvvhi3PP69++PPXv2gOd5XHTRRYphuOSSS2AymU7kLbQJzJnNYLQzxxuq2pyzPDLQ3nHvZAS9Lgi2bKRecStMfQfDkNsfjjXzUfnB06BBH9K69+hQI3G0IcQdYVhcLhc+//xzxfm8e/duSJIU99itW7cmfJ2nn34agiDgsssug8Viaa/mdjrYioLBOAYKClZi+qxHUHX4EDQp3ZB86U3QWOxHpZHUmvDb4121tDeJ7kHa/hIqy0tjjk+kKzXpD2NQuHlLmxkPURRVUUn/+9//EAqFWn3+4cOH0a1bt2O+fleBrSgYjHYierDLCQ92VWsXwmoxHpWQXmu0jBKtWuwZmTit11kdvt2TaFVUcvgQCgpWxrRpztynQM68HM4tyxCsLgWnM0HyufHK8gKk5c1us1rbhBBMmDAB1dXVrT7nvPPOw9ChQzFkyBAkJSUd03VPVtiKgsE4So43ES9668WekQkqSaiuOhJ3wI83A6//aCloKIDkEbPaXe21ubbnnt4Lbrcb3OCpMX1Rvel5ZNmTYvqDcBz4pHTYhz2otL1y9Xykj3n0qPozEAjgyy+/xMGDBzFp0qS4x4wbNw6rVq1KeC9nn3224mMYPHgw0tLSjqYrTgpau6JghoLBOEqOZzuooGAl7p4yDSFej1BdBTTJGdCIPrz28gut3tN3NTSAHzqtXTPGE2WFNzVadRsWISBSVfJfdeFSJF82Ec7C52L6Q2tKgnXU/6naXrwwDzkzVjfbn8FgELt371a2knbu3AmPxwOj0YiamhpotdqYe3jhhRfwwAMPKP/37t0bQ4YMUVYNp8LWUku01lCwPApGp6Gr1Ec+ntrQD82YhSDRwHb9A7ANfwgA4GuoxW1339vq+3VUHm42j+J4+zFRfshDM2bFJAEmj5gFEvKjetPzKMkfC+eWZUi54lZoLPa4/RHyumPaLlizY/rTV/I9umVmY9GiRRg+fDisVisGDRqEv/zlL/j444/h8XgAAB6PJ2FdhmuvvRZ33XUXVqxYgbKyMuzfvx+vvPIKJkyYwIzEUcJ8FIxOQVsK8LUHTbeLGjYsAqK2fuo/WgqPGACnNYIGfbCkpEJvMMJRUa7aUqp21iB93GMQ3TWo+7wAtqgtmHj3G69fNEcWo3bn20i94hbluIihaot+TJQfUvnB08iOY6DEoA8WgcA4fl6LdSNyz4j1uRh6D0Ll2oVIvXwSJL8HngP/QeDwfoBSzJ49u8X2bt26FZdddlnM471798Zrr73WqntmNA/bemJ0CjpagK85EvkJTFoe1VVHYM/IRG1NDSStLK0RanCg9rO3YB8+XVVcKOSuASiQM3M1Dr/xJ1ivntzi/Sbql8pV80CD3pgKd3PmPnXc/Zhoa61k8RhkTJifMO8jeqtq2LVXx41iiteXjsKlILwGIefRlZzJzMzE0KFDMXHiRAwbNuyozo2mK+aBtBVs64nRpejMdRbi6S4lXfcgzBYLJFGEyWQC1Scp+kv1X/w7RovJNnImBGs2eLMV/rIiBKvLYu431OBAcXGJassoUb/QoBc5M1bDes1kePZ+jEl/GIOJE29uk35MtLWW1r1HQkmT6LoR8+Y+jhXvr1ZtXU1+aCYWLlyE4cOH4eUli2D4ajlK88fCufavEOsrYdORFtuVkZGB8ePHY9myZdi/fz/KysqwYsWK4zYSxyvDcirADAWjU3A8+/7tSUHBShz8+UDcwffgzweUwTxaf6mpEXAXbYfzYzkcVAp4UPHvJ6FJzlDdr7toO2o/ewvp4x5TDVi2tG5x+0WwZauMUOHmLQDaph8TaVwtWfRM4yD/7FgYvloeN9JqztynYLjyfvAWG1x7t6Dhu03wuFx45JHZWLNmTaNRkUQE3PWgkoSfD/wYI4xnt9tx4403KsV6Dh8+jHfeeQf33XcfevfuDUJaNi5A8z6bROKLTGBRDfNRMDoFHVUfublth8hsU5PSLa62kialG6ZMnwVLig0eSaMcI9iylL/dRdtRu+NNlS+iat0ihGorUbV2oRItVPvZCmUV4i7ajrpd7yJYfRgenRFSOHw0eqsm9YpblbZES5+3RT+2pHGVaFsmUqzn4M/7wR95FqIrNodh69atuOOOO2IeN5vNuP7668HzPK688koMHToU55xzznEX62nJZ9ORUvJdCkppl/8ZMGAAZXR9VqwooD3P7EMJx9GeZ/ahK1YUtPv1LGmZNGPCApozcw3NmLCAWtIylev2PLMPzZiwgNrzZlFNcobqOM6USi0XjqQZExZQIugoMSRR3mKjGRMWUNuIhymflEYzJiyggi2HZkxYQHMfWa/8ZExYQInORG0jHqaalO4UIBQgNGfmmrjX4i02CsFIiaCnAKG82UbtebPUrycYlD5T+pFwVDBaKAiJ6c+j7et4x5eUlNDly5fT22+/nebm5lIALf5kZWVRSZLa9X2NJvIeNu3/nmf2Sfh80qDxVDBaTtjnsCMB8BVtxRjLnNmMU5Z4juKaHW/B+91GhHxugNfBeu0UmM+9MmqWXwZwGvAGC0SXE3xSGsT6Kuh7ng9/+QHwxiSE6irAGZJAQwHQgBc5M2NzBEoWj4E9byZMfQfDV7wHjtXzYB8zB84ty+I6uStXz0P6mDmNq4qNS5By+S3QWOyqv6tWLwCVQqBBPzidEeYLRyLl0ptUSXlAbD5E9HMt5U94f/kKjg//CiodnYSITqfD0KFD8d57750wnaSWcl4KClbingceRtJ18oqvdufbcO/ZDHveiU9m7AiYhAeD0QJNtx3cRdvhLtoGU//h8B7YhaCzFM4tryBw+ACs10yGqe9g1Ox4C649m2EfOUO1leQrLULGH56Qw17DBkWTnAERNO62FdHqUbvjTQCAsc+lEENBODYugdjgiO+8DnhV4ar24dNR+cFTAMfDfM5QEI6Ho3ApJL9bpT1VXbgUWnsOTH0HKxLoAOKGv06f9QgaXG6EeD0oBUqr6nD3lGmwmE2q4w29fg+iN4N66prtX0IIKKVIsdrx4APT8Oijf4EgCKpj2jviqDXCjTQUQPWm5xGqqwDRGlRZ4kcrHX+ywpzZjFOWpo7ful3vwtR3CDxF22C9ZjJyZqxG+phH4d6/E67vP4WveA8avl6PtLxZKudnWt4sEI7IYbE73oT16snImbEKtusfAA2FZEMS5RiuWrcI4DSwDXsQNTvelNsghpBy+S0gWmNcZzRvUctL6LL6goYCSB/zKFw/bEfN9jdgH/Ygcmauhu36B1D3eQFEdw1swx5E3a53lXOKf/0pYWRU1eEyBEIiBHsONKmZ4AQtgkSDqsOHVMcTQmDoeWFMfzYt1uP1ekEpRU11FebOfSKukWjviKOWik/NmfsUkkfMQo/7/oHc2WtBA95OG33XkbAVBeOk4Whnp00dv8HqUoACtmEPqmaUaXmzlNk79XsShKv6UbfzbRj7DgkL3pVBsGXBcMYAeH/7BpWr54P6PSA6IyBR2K6/X05Wq6tE5ar54JPSoLHYAcKhat0ipEVtfTg2LgENBVXX9JcVgWh0cG5ZBk7QKo7wSJttwx6Eo3ApetzzirxdBvVM2ltWBCHjdPhL9sJXsgfen/8HAJC89fD+/IV8EY5H2pg5cKxbGLMq4s1W5e/kVBvun3IfHnvssaMq1nO8hZ9aQ0uO+YO/HIBQswxBp/x+8UlprChUHNiKgnFScCyz04kTb1aFexJBj6CzVDEE7qLtKH99KirenQOAwHLhSBBd7Iy/dufbIFoDQrVH0PD1ehh6DYJt+EOgoSC8P30hb8H4fSBaPajfA3AE/kM/KpFTIBSmc4ai6sNnQDge5n7XwvnxMpTkj0Hl6vkwnH4RJE+dalbs2LgE1munwHr1ZIguZ1zjJdZVonbn2xCsWfAV74H7479j3A0jcXav01H17ydQtvQmVK2ej4bd6xCqOxLbQZIIgIIGfDGzctf3n8Dyu3GwjXgYbqrFggV/RU44MzzyfrQkI3KsOR+tee3oY+bMfQrz5j4OSRTx20/7VFFtGlNqePW4CtarJ4MGfHBsyE+4AjlVYc5sxknB8WR2FxSsxEMzZqG6pg5UCkFjtiJUewScKUU1s68uXAoh/TT4y/chLW+24vx07dkcuwIIBmA+/3rZ11FdGvNaVesWQfJ5wJuSIdZXARotCM8jfcycJiGypSBaIxAKARoe1O8FEXTg9GakDrkDpr6DceiVe2C7/oG4Cq6it042ToQAR/tdJxws5w+D/sh3WLLoGWW1xmkNoJwAyVsHojXC0sRhPukPY/Cvt1bGCB/eccvNqmztYxE3TFTPItrZ3JpjmvvMVH7wJOz2tBj5lZboihneTD2WcUpxrIquU++fhlf/+QZEnwecMRmE42AfORPOj5fBek1s9FHVmgWQfG4QnRHU7wGnNyFt9F9ijnOszwfhZT9E1ZoFcY+pXPU0eGMKQrVHQLQG0IAXgi0b+tx+8P7yP1XuhWPdIkhhn0S04Uq54lb4Sn+A58B/4iq4Vm9YAnA8ILWmaA8BCIGx9yWy85vwqN70d1gtRjyXLw+y55x3Hoq+/xFEqwUN+OR2iyLsw//UGMG16mlAb1ZJmDjWL4YU9KvaX//RUoQCPljzHokyoAsx9JKLsGXz5rgtbKl87Jy5T+HgzwegSemGlMsnyfeB+AYo4Wcmfyyko4joikw0nA0e1XsQqVES6bvOCDMUjFOK1q4oomd9trRuqHF5lVDI8temKDPz4oWjkDMjjt5R/hjYhk9H3c63EaqVt2sShb8CAKeXC/MkOsbQ6/cIVP6mqs9QtW4RzP2uVYn+RVYIPe77h/KY99fdcKzPBxVD4HQmEI5HqL4CRGuUI6F0JjTsXgca8CTuOF6A5fzroc/pB23mWTj00q0gGj1o0BdTuc+i41FeVg6i1YEGvOD0Juiyz4W/fD+kgA+5D78vD7RLJyB93ONxVzjR7Y9oVoEQ0IBH0a1y79mMe2+7GS+9+EJMc5vTodKYU2EbOVO1sgOlchizxQ7RXQtCQ8psvy10sZQtzyBNuKozCKTThtcyrSfGKUVL0S1ArB+j1k9hj4pgipbhiGRXR+MvKwKflC6rvl7/AHJmrlaytpseRwQdBFs2zBeMANEa4h+jM8Jf+j3sYed5dBSV98Au1fG6rL4I1R5G8TN5KH1+Ispfm4rKVfMgeetBAx6IDVWwXHQDcmevkyO19n0O157NSB87B9nT3wd4OeKIMySB05thv+HP6DF1OSCJsF59H4y9ByFUXQpNcjdoBA0yJsxHj/v+AfO5V8ptO/NylB+uBGdKQvrYOciZuRppo/8Cf/l+6DL7gBCgeOFolOTfCBr0x/U9hOoqYh6jAS/SxzyK3NnrkHn3S0i94hbY82bhtTfejPs+N41UcxdtR+nzt4DoDLCNnKnqR/vw6SC8gJwZq2Af/hA4nRGU06H0iAOTJk2C0+lE3YZFx+SPiPhAbrnjbhivmqb67DS955NBEoQZCsZJQVPHdDwdoqa6Pk2/3NHGQZ/TLyas1VG4FJBEJSqK8BqkXD4Jjo1LYhzN0Ghlwb6ibdBYM+OGyJrOHgzJF1ufQZfVF0FnKagYgv/Qj6jb9R6OvDUz/CyF5KlDsLoEENWRUL6SPcr5NOhTwng5rQEZf3gS9hv+D0TQg4oiHGsXomLln+XEwKh2a0Qfgl5XTJvcez+Rrx7wouKdR3HoH/dBdNcgLW8W/KXfgwZ8AJUAKiY0jJrkjJjHiKCPf/8eV8x7XFCwEk6nExXvPIqSZ29Ecf4fUL35RUAKJgxrDdVVKIbD3P86cDodTOddDcGWjfoaBwKBADwbFzWrXdWUqfdPw+2TH4B34G2gIdkoJppYCLaskyK8loXHMk4aJk68udkvedMEu2hNJkA2DpWr54MGPOB0JpgvGKGEunI6I8wXjED9f99XRUXV7XoXYr0DlaueBg36INiykXL5LajesEQJU3V+vAySu06WBg94QbQG6DL7wHTWZfDs26Fqg+itR+32NwEQlD73R9BQoNX37y/5HpRK8JcVxczq9bn9QMUQHB/+DRkT5qu2uUoWj1F8JJKgAzgNyl66XW7f/v9AdNeAaI3gTUkqv0PV2oWQPPUgWh0AyFtSkgTLgDw4Cpeqt9PWLgQVg/AV71E/RqW4tTUEo1l1b40Z1NOR0mSLrv6/7yvFj5qGtfJJaShbdhfEuko56izgiwk+qF6/GG+9+VarHdav/vMN2Ec/qtL1Sh40HtWFS1V+pYgP6WQIr2U+CsYpw2m9zoIz7Xx4f9qFYHUZ+KQ0SAEv0m/4c0wNiZL8MarynO6i7XB++ho4QQ/b9Q9AdNfEiv19+AyIVg+xvlKeuQd8IIIOVBLB6UwqR6dj4xJIfg/0OeepfBSen76A48O/HsVdEYBwSLpoNMz9r4PY4EB14VLQUBD2vJkt+glqdrwF13cfqdu2IR+ixwWIAWhSuoGGgiAaIe4efFNpkap1i2Dqcyl0Pc5WorY4neynAQCiM4D6vSr/R+Xq+YAkggb94M1WSEEfJt91u8pHkcgH5dyyDABg6DUInqJtqvejcvV8gFJV3Y6Gr9fHrc/t3bwEddWVrfoMHfzlgPLZiBZ9DDU4ULt9OUSXM8a/09V9FGxFwejSJKrtHC9Mcdi1V+OV5SvVYaprF6Li/SdAOF7lgCU6M0qfnygnyQk6JdmO6gyoXDUfRCOA0xpQ8e4cOflMkiB566HRGWEbPl0eAFfNBxG04MPGRZdzLoKOEgSqDsJ01uVo+KYQ/tLvIfncqPxAXpEQQQeAQNbQiwPhobH1gOhyIvniP0KbcQYc6xej/n9rUf/lKnlFc8WtqNn2LzjWL4Y9yrlbtXYhUofeqXo574FdSBs1Wy0PMmIGHOvzwelMCDpL5dWGK/7WTlNpkbS8WahaswC6HmdHGgwqhgDCgzNaZEMbDgEGZAPM6YyqlYpzQz4uveQS1bUSqbwGq8tgHzkDtTvelJMdP16GoLMU4ARwOmNMJFiihMkKpwNT75+mMk4FBSsxfdYjqDp8CETQw2ZNhaPysGr1Eomqqt70PEK1RyDYsqHt3guB0r1wbnwOuWeoE/y6KsxQnEJ0xTjv5oiOl7f2caBs59uYNGkSeFMK7HmzYmSlC955T9m3B8ID26jZ8nZTKKAM1hAM4LQ6pOXNilutrvLDv8mZ0loDAAop4FXlElQXLoXhjItAQUHdNZAAeWsq4FXaLthzQYM+0CAAwQCi1SJ93GOyE/r7T4CQHwDk2XgoAMGWBfvImRDdtXBu+jssF4xA3c63oUlOg7nfteHophBCDQ441i2Wt1lEMaxhdAREawThBTn7O4roBMMIuqy+EF3OJnpWC+NvEVmzY86VfO6Y1Vbl6nlIy5uFinfnqK5Xt+vdmKxy64gZuOWOu5REuYkTb06o2STYsmDqOxj+Qz+i4ev1cs4IrwXR8DEGMBKqnEgy/rU33lQMxdT7p+HV5Sthz5uFnKhVIKc1QWPNVG0z8aZUUCmkiDwCnac6Y1vR6ZzZhJC5hJBDhJBvwz/DO7pNJwMnYyWviHNaqT99/QMQbNmqSKZIIZrpsx5BvTNWcM+973PQQACEF+QZvdYABL3gBD1Ed03canWEcOD0RjnyKawH5SnaBve+z8FbbDD0GoSGbwvBcbzzDqHsAAAgAElEQVRynWgjAQBBRzH4pDQQQQfCAWl58qDm2bcDqUPugPWayeh+10vIevAdZPzhSYRqDuPIWzNRvXEJjH2HwHtgF2jID+vVk+Eu2ib7PnQGJSIpfexjgBhE5t0vQ7BlI33Mo7BedQ+qC5eqnOqJtKU0Kd0gumtw+I0/oeLdOeAEA+r/t6aJQ34hDL0HxZxLtAaVw1+f209xNjd1+sar9BeRRCkpP4JJkyaB0xpxpPwQqtYujAkIMPQaBNf3n8Lz8xdIH/OoHIlmsSZ0bkt+d0xgQXXhUiRfepPiQFf8EE0+R/bh08EZLPAV74EU9KPyg6dRsngMKj94CimX36IYici1uroDO5rOuqJYQild3NGNOJk4Ebo6J5JI5bmc0X1x+I0/KQNT0CkPPNGy4II1C8HqQwCvQflrU+QEt3DCHBH0IFotLBeOhLtoG+zhvea6nW/LM3NBh4p350CwZSF50HiY+g4GDfqV2WqorgKhuiPQpHRH9cYlYdkLGcnX0Ow9EK0R8DWA+r2oWrMApr5DIPncsJx/vSpPICIAmDF+Hhwbl6Dh202g3noIdjk5zz7sQTkfQZJUbY2E7kYGY8/+naAAKt55VPahSCGYz7smxglbtXYhjL0vifXBrF0Ix4ZnIbqcEKxZkNx1cO3ZDENuf5WPgga84S2gMqUtke2apk7fSKW/GHVdnRFpo+REvKoP/yYP7IHGgABNSjdIAT/ce7dA8rtV24ahugrwlviaTURrhC6zjxK0EMk50VjsigN9ztynICbYogrVVQCUIn1slF9m9fyYlVo8B3ZXXtF3VkPBaGPas5LXif4CRFZH8epPC7Ys1O58O8axWbVuISR3HWgoAGJIQvoNj6gGt/pd78lG4f0nQDRqwxE9UFZvel4VUeTc8mqjiF4rIBodKJUAMQjR7UT62MdU7SA6Q9wBTrBlKbPa6k3PQxSDSB40HgAQanAAVAIN+WWnba9BqN3xJvSnDUDV2oXQJGcofWJv0iecwYKUK25Vort4UyqkgBe+g9/EiiOOmg3HusXInb1WcYyHXLXhbTV5NcbpTOBMybBeMzlG+qRq3UKk5c2WM8bDe/pEZ44RQaxauxCWC0dCn9sPzo+XwV++Xz0wr1sEw+kDYOx9CZxblkF0q7WuIg5xx8YlMVFa1O9BoPI32K6dqmSSV296Hq4ftoKAghAORNCBN9ngLyuKkY3nDEngtAZVv5gvHInq9YtVyX5Nqwq2VGmvuc96ZzAunS7qiRAyF8DtAOoBfAVgBqW0Js5x9wK4FwBycnIGFBcXn8BWdj0SRY1I219CZXnpMb9ua3V12pLTep2FCpIK7y+7AY4DKAXhODlqxpQKKehD+tg5MfcakeVwbMgHQCA2VMuZ0363Eg7r/uFTSAEfIAaRPu5xaGxZ8JfshbHPpajd+bayDx6RiBDdtaj59B+JG4uwcQj5lRBUTUp3hBqcSPrd6Jjs64r3ngBnMMdoTKVccau8mhFDKMkfC1AKPjkdYl0lOGNyjNPW2HcIPD/uQKjOASJo5FlwnGzpeFFLvMWG4JFfkDNzNTz7dzZZmZXC0Ov38JXsDavhGmC5ME/xz1SGy7bG6ic9BRBO9u3QkBIJRXRGUEkCAVWinkSXU8lkL106Ia78SdWaBciatgIl+WMhWLNgvWayMqiH6iqRPu7x2Nogfjeoz4Xc2WtVK06i0QK8gPTR/6eSHBF9bvDGJiHB6xaChkTYrp2sbDVFMsN7ntk74YB+LFpkJ+K71akzswkhWwgh38f5uQHAywDOAHA+gMMA8uO9BqX0VUrpQErpwLS0tHiHMKKYN/dx1H+0NCYxrMHlPi4/xbEWp1fUPTkeWlMSCMfFKIFGH6PRGUE0WnBaAw7+fAC+g9/JQnehAIhGgGXgDciZuRr2vJkJ96eD1aWoeHcOCK+F6JUHKcnnlgdyMYT6Xe9B9LlAfQ2gYhBV6xbh0Iu3wrFuEUoWj0bD7vXgDUnydQFUb3kFUtAX/wZ5AeA0IIZkpN/4BHJmrkH62MfAW9JgPPsKcDp5z99dtF3VRkgh8BYbqtYsQMniMaje9LxiJNxF21H+2hTZMOoMMPUdAsGWrWyDRfrfNuxBuL7ZIAsbGkywDBiVMFua+r2oXD1f3mtfPR9SwAfJXQsi6FC7821VfQ3rNZPBGZPhK96r+APSx8yBa89mlCweA+eWZaCBxDLslgF54AxmgNcibfRfwuc/Co0xCdZrpwKEoMfk10EEveLHSJSQKPncykrL0HsQKlfPR832N2QF2FAAuqy+MPUdjMy7XkLu7LXIvPtlUF8DNMkZSkhr5L7Sxz0OTiv7pBR/xMiZIBwf459Ky5sN3mBG7Y43lffOX1YEXm+Mq1Ab4ViUcpv7brVGQbct6XQrimgIIT0BrKeUntvccR2ZRzH1/ml47Y03EfS4wOuNMOp1cNXXthiq2d5ElqwHfzkQE9MPEPDGJBh7XQzvr7sRqj2CtO49QCUJ1VVHYEvrBsJxcdUzo1+X0xoghdVMqUQBMaDM0IlGj8n33KmohVpSbHC7XBAjQnIBrzx79HlAODQKzAW8SmIUAOVvojWAUoqMcY1bNY4N+Ui54ja5BOjahZC8LhCNoISZUoqwOmsliFb2SQi2bGismfCV7AXhNZA8dfKgTsNlnROIwRGtEZzeFJtw5nMDBNB1742AoxjU55JnyQEviEYXdxYfSUaTa1xwoH4vBLu8ZdTw9XpQiYJwRBYdDK8W4kVfVRcuRaiuMqGOFBH0SB/3mDyAh4IJ8iDmq4T6IrUvCKeBFPQoarbR58TTbIqs1hKtKKo3PQ+iEWC9Wj4m56F3Ez5ftXYhACgRaXFXKFEaUURrBOF4pN3wZ+hz+6H89akJysnOBzgevM4Uty+cW5Yh866XVH0Yt2/zxyJj/DxUfvAUOJ0JVAzJSsHfbcS/Xl0W9/t9LCuK5kQLNToDQpQqasKgFJPvvTuuPlZztMmKghDycHM/R9WiVkII6R717xgA37fHdYDmde3jzXjTM7NhNCfLUS9aAwgn4NXlK2Ed9X/ybHb0o/BIPKzDHoIz7XxMuv1O3Hrv/apIo1vvux+EE+TX0Bmh0eqR1q0HCMdBYzCHX9uItG49Eursp3XrIV+fcCAanXJepJ2coFeumzNDjoDhk9NhGXgDeFOKnBTmaYB73+eKZhE3eCqcLi+oRFFdU4Oqw4cAXovi0nI58kRnBCEcbrnzHhz8+QA4YzLMF+ZBsGXL2yoaOcKH01tgGz4d+p79sez1f+HgzwdAtEbU1zgggYDozLAMkM8T3U4QnkBj7SG3b0AeOFNK2EgQcKYUpI99LBzFMwecVg/Pgf80zvpGzED9F/9W9s+JoIVl4Ci5TUE/OJ1ervoWjkzik9OhP30AvL98BUgiOK1B7lQxKKurNqMYSoO+2NnlqNnQJNnB6UwIVPwCffa54EypsFyYB9vw6Qln8WJdpTLLTh8jRylZr54M157N4LRGcFqdMmMHAMf6fDg3vxRzfduwB8HpTQmjliLyEsHqMiRfelNMxFO0LyA6socGfRBd1aD+xLIYTR8LOsugz+0Hy4UjY6KTlKiisC+J+j2xr1l7BIZeg2RDJQYheerg2PgcaMAXR/5kIQAi91H4vZW8DUpbIw7zphFS+pzzAACh2iMJVpxlqj5MJEUSkeWgoQDsI2eA8BoI1h4Iel0JIwlbo0XWlKa6VpHrc1o9RF5QPjvp4x4HZ0zCstffaLeVRUtbT5bwz0AAUwD0CP9MBhBbC7FtWEgI2UsI2QNgKIDp7XGR5sJFo7VcsmesgnXU/4FPSofLnAVfsFFfhwhC3BC62m1vwLVnMzRma8y2QFrebBCtDrYR05E+9jFQrREOZw2IzgL76PByftxjqPEGccud9yoGJfJ70qRJqK6thWXgDTCfPwyE1ygS2Yb+w+W4dkJirmsf9iC8P+2Cffh08DoTCM/HaduscIjlY8oHkEoh1YBtGZAHaGSJ6fpd74UdqQAnGACNFqHaI6je9CJ8B79VZjtSKAhz/2HyuRDlLR5vHWzDpyN9zByI9Q5wBgvce7cgLW+WPIvT6uOWHHUXbVP6P/rLHUn+ipQxFWzZSshp5PykAXlo+N9qQAqBBn2K+mtLaDP7AFRKOGimjZoNTXI60sfOkQX9ftqF+i/+DSLoEggG6pXBvOn9SX636r4lbz16TH5d2U5pen3J744ZFB0blyDl8klKKKpgy4LGYlec1iX5Y1G96XlIfjdSLr0p5jVp0I+MCfOPSrNJsGUBAFIuvUlu06bnUZI/Fs4ty5Byxa1yVFG4PURnjNMnOjR89SFET73ss9DqYR/+EAAKY69BylZc1ZoFoJKE9LGPqvouWpzR1HcwUq64VW7D4jGyeKK7FoHyAyAgCYUcNckZasOS2w+OJn1bXbgUyYPGqwMMRs5A3c63IdiyE267tkaLrCmJjIsYCsHc7zrlvXRuWQbTOVeCNyS1m/hgs1FPlNInAYAQsgPAhZTShvD/cwFsaI8GUUpvafmo4ydRuOhDM2ahpq5O0XKJPGcf9qC8DI6KUql459H4M0a37HuXwv83fZ4GPKj/4t/IvOsluczm6vng9SbocxsL1ogNDoDTyjP3/tfFROBEJCDSx86BY+NzkIKBcNKVV7lO0+tGZnSRED/3vs9VZTuTfn+jsmdN/R4QrR4I+kAJp3qMcBqkj43arihcClPfIXDt2QyJArzRDPuIGaq2ur7/BLzerHKcOtbnI3XonfIWw6p5Kgd0wjj4sBQEoB6g5BmgHqlXy0YiEiYbjemcoaj59LXmPxiEgzbjDOhzzgOnt6Dhu03ofks+ypbdlTAaqanBivzNm6wxkTeOjUtABF3cLOmIsYtud/RgH/f6YWkK5f0RdLCGI3oAKI7tyHZX99v/rmxbRXIoYvSRzFY4CpdCl9knNiJp3UJQSVJpNkWc7cr7oNFDCvqQMX5ezGekat0iQKKxmk9BPwR7NpIHjQdvSkXF+3NlkcLkDJjOugy266YqbSxemBfz3iZfehOq1i5UnPqRJDiiNYKKQSQN+iNSr7gF7qLtqNn2L9Wx8jbmEgAUJfljQQQdLANGKcfL35FSEEEP6zVTwJtSVfccWRHZ82Y2G0nYkhZZvOOB2DKukyZNih0PCpdCrKtEcX1FC696bLQ2PDYDQLQ6WSD8WJclUbhoqbMGNORLOMBHG4/IzCRelmfm3S+j/LUpCb/cQUcpyl+fiqCjFESrR6j2CMqW3QUa8CEtHLpZ/toUGM++Aq5vNkDyu+HY8CwgSWGRNnnvPhJXzplSkHaDPAiXLp2QcFBTZk6hADw/7YqpzAaNVrVnXfHBPBBBG3OcyvE37EE4tywLG715MJ17tdrIDp+OytXzYR85Q/34yBlwblmG7rf/PWaAjMxm48XXUzGk+CiSL70ZdV+sQu3nBUDIj6rV82A658q4A2uwqhhEo20itEcgpJ8GIbU7PL98hbRRs2E4fUDjACaJ8BXvQcplE+HYkK8ygNGib9EGK/K3odcguL7dFM6Oliu9ST4POJ0BVAzFjb1vet/Jg8bDsXEJTOdcGVdoT/LUIfRVpewjEvQAIcrrRoxF7WcrINZVRhVdkiOVBGuPmHuqWrsQkt+NpN+NReoVt8D58TLlPE5vgqnvEOh6nN0oOaKTK9wZ+1wqr2Y25Murr0iuQnhyQQM+NOxeC31uf5jOulyZoGiSM0A4DXIfWaf0AxVDgBhE6uDb4fzkHzGDejwDp7HYQUPBRuFFQS/v3ROA8ALcP3wKQ25/GPtcCt6UisoP/xbOpZCPtQzIUyLQip/JU1Zapr6DG6PNFo9B9YYlcsRbOMAg8r7xyelKuG1bCgDGMy633HmPIk0PqCeyudlZbXbtaFprKN4E8CUhZHX4/9EAlrdLi04Q0ZIA0WUnQQHBFl+JsqlcQdNZTGQQTR18u0qCuqkDkhgsIDoDgo7YEpnRg3Co9gjcRduQNvovsjNzx3KY+l0TM5uoXDVPJU3B6c2xM9nwjC5SppPwQoychX34dDjW56seIxyJe5xzyzLlixK9WqEBL9x7P1GE9xpXKvGjYYLVZcpqILrP489mF4H6vShZPBrQ6IBQAM5Nz6vfWEmEv2SPvPXQJJHMUbgUQvrpCBzeD/ACki/+I8z9r0OoujRc5AZwrF8sO891Rmi79UKg8lc41udDdDkBXqMMkJqUbki+bCJ4U6riVI9sWZj7XQvB2gO1n70F8/nXw/3DVqUMqfWa+2QdqA//puQVqGbW4T35yH3zplRIfg8avvoQNOiPylnQwzJgVGNY6od/Aw345WitJjNrKeAFeI1sJAQdaMCvkruo/OApJbTY3P86uH/YCk/RNhhy+yP1yrth7H2JLMERFaZat+tdxRB6ftyB+v++LxtCvw80GEDg8E9KQiOnM8E0YBQavi1EsPI38ANGKSubqrULYe5/neotjB54jX0uRcniMcqMWZPSDaazB8fNkQAAy4A8ufyssxTgOJjPGQpt995wfvyKymATEGXldeSdOarEwaafxUibIgaPSiF5xRI1YYm8/03zJ9oDGkygu+X3tFtt71ZHPRFCLgRwefjfHZTSb9qlRcfAsUQ9RXwU5MzLVQNvZBYfL2HL3O+6mLj3qrULASpB8jYoy37C8UqMNm+yQvK7lGQoCFqEasrlVYAufhnNSPRFdAx5JJLDuWVZTERH8TN5qugMd9F2VG95FbzOqOj8RLYlKKVyFJS3PmFER+7stQlfO95xkTZbr54sOyBDgRjjKXrqkXHjE3GjYagUguTzgGgE1XkV78+Vy3gG5YEOYmvKecpk3vMqXN9/goav18m+kkgklUYHhPwA0YAImphoMIjBcNayKDu5NVoQwjVGUgWDAM8reRvyuXICXWOEVvhavKAcF/23cj1JAhFkf0+02mjFv5+Ujw0bLNPZg2E66zJUrVsEY69BcH1bCN5ig+RzNb4epRBSukP01kGX2UcRGyRag9xmhOTKdZKIpN+Ngbtomyy93SQPImLotPYcVQ5CqK4KnNEcUyvc3O9aZWAmWiNAAX3uefCX749JojP3vy7mdUWfC5xWH7M9lzr49sbSqmvmwz76USVjPlR7BMRgAccLEF1Oub+DARjOvAi+4u8UufdIJr2veA8q3p8LTtBC8rnjlkmt+OBpEELCW7ccOJNFbcTXLQRvsUPjqYbPI9cgp0E/tCYLksxmVFcdOWFRjW2ZE9Ue6rFGAPWU0n8RQtIIIadRSn87qlZ1IiJv5h33ToZ91P8pnZ5y+STUbH8DpnOuVJQoiVaWRnZ995FKrsCxQS61SMUQMibMl5fTzkNxjQylVJZcNiUr+/Ql+WOajb6Q/I0x5JEZezxtnKbbFaa+gxFwlKBh97qwCClF0qA/QvI2wP3jdlBfPTi9Ka7IW7Sj0l20PeEWUMTxF71akeWtA7FqpMOno3LVvISzQEolgOchBbyoWvNXSD6XrIFECKgY3iY6CiMBQlD+j3tBdEZoUjPRfdJiZeYnuuvDA74fnDYJEuFAA36AyhFPCcua8ppw/oSE2K8NafyT1yrGBxotLOcPk303Pjeg0YNwPHhDkjKzFb314JMzIHrrUL1hCYjWAMJpYBkwEvVfrgFvSILru03wl30P65V3gzelwrNvB+wjHpY/V0E/wHFICg/2jXIXct0LgCBj/JPKZAMaLVx7PlIN+A1fr0f9f98DEQygkihv0wyfrsz6HesWIenicfCX70flanlrR7BmQ5fZJ8oQG8IrFA38h/bB3P+6qO+PETQUVK6rrCbWLYLlguEQrD2UTG1OZ4D5wjxlK8vzyQu4987bUbh5eWOYtVYP0dsAzmjGvVOn4qUXX4BGZ0Sg4hdYBoyCu2gbrFfLmeER577t+mnQWOyo27AIgaBXtSKoLlwK+/XTFKNi+Go5jHoBP0bVD9FYM6EXPXj5pRc6XHZj3tzH5azuOIl47QaltMUfAE8AWAfgQPj/TAA7W3PuifgZMGAAPVYIx9GcmWto7iPrlR/biIcpEQyUcBzteWYfumJFAaWUUoBQojXIv3UmmjRoPLXnzaIAoTkz11B73ixKdCaaMWGB6vWSBo2nRGeiRNCrnhPsOTRjwgJqz5tFBXsOBeGoJqU75YwpNGfmGqpJ6a4cHzk28jvyGva8WRREoJwplWZMWEBzZq6hGRMWUE1yhtwerUG5RtNjOFMqTRo0Xvmft9goMSQrx2hSulNDr4vjnJdCiSFZ3R+CnoIXKHiB2kY8rLr/nJlrKAiJuU/wAoVGR4khmdrzZsltNKbQpEHjqWDLCSc2tPKH8PJvQd/YLsJTIkS1T6OjQrczKac3h48xUnACBTgKEMrrzTHvXcaEBZQIekq0RmrodTHNmbmG8snd5HPDr5GZlU0ppXTFigLa88w+6n7RGqnGniv3b3IGhWBSnSt/lowU4ChnTKEZExZQ8/nDw4+RuO9Z5PMHEMqbbRREkO+vyXVBeEq0JsonpdGMCQuobcTDlE9Ko0K3M6PaIB+f1j2LrlhRQFesKKD2jEyl39K6Z9EpU++nPc/sQwnHUUtSKuX1prjnRpgy9f7GY8KfC0uKlVqSUhvfj/C5vN5MCWn8nkX6sOl3ryWmTL2fEp1J/lxF9Y8hKVW+n6jXu+qaa+T+JY3f40j/WtIylWsea1tOBG3VNgBf0dbYgFYdBHwLedr0TdRje1pz7on4OR5D0fPMPnEHh55n9ok5VjBaKAihIBy1jXiYapIz5IFEa1ReA4SoDI89bxblw8eBcDHPEUOy8kVuOhjwFjvlzbbGL7nFJhsdfVLUF8JIkwaNpxp7rvLl16R0p/a8WXJ7wkasqYGJ3GfjoGmQB+6oL1nkXMuFI5XjIr9zH1lP9T0voJwpRdV2PimNEp1Zub4y2OrkPsqesYpar5sWHtg4SgQdzX7438qxnDFFvjfCUd5sTWgYiEYnt1mrV/qH6MzKuZqU7pToLZRoTRS8VjYGTQwHEfR0ytT7lfd3xYoCaknLVN1PZOBYsaKApnXPCp9roPaMzGa/nD3P7CMbPHsOBWnsN6IzUUtSatwvedMv/5Sp98uDa9SgOmXq/XEHiUQDR9OB35JipWndszrl4He8TJl6v/wdBaGC0aJ6b5s7lgsbjZOtP1pDaw1Fq3wUhJAvKaW/I4R8TSm9kBBiArCLUtqvxZNPAMeTmX00eiqE4+SSkIYkpcpZzbZ/yZpBpmSk5c1WMlQjWy/RWaLxMkbLXrxNFQ0ENGaR5jz0LipXzVO2EsALgCSp9oqbhic61i8BZ7Qo+8OlSycgfdzjqHh3DnJmxGZ5luSPkZP2bFkQGxzQZZ6FYOVvSgnP6HuJtC3iQylZ8se4mkrVYQdz5t0vw1f6A6rWPiNv0fk9yhZPNClD7kTSwFGKVhAn6GEfOQOcORWHX5sS3vIhAC/AfN7V8JXsQeYdcow8CB/e5/fG8QNQWPpfg9Qr70bJ4jGwjZALClWvX4x7br0pbhZrW4mwNZbuVIc0CzSE117u+O0LBgNovY+itYZiJoBeAK4B8FcAdwJ4m1L69+NtaFtwvBIerR0c0rr1kJPjODn5zXTOlXAXbQMhnFLdSvZDNEYyRcsARJdNjM7FSLQvnjTojzH1fZvmGwDRJSE5iN5amPtd1xj5AR5Eq1XKTzZ14lWung/CCwAh4C02hGrKYbkwHDkS9qmonXqyszPl0pviShxIoSBK88ciYYW2ePAaQBTlRCwKWZZBZwT1u8Fb0iAF/bBcMAyeH3cAVELKFbeCN6Wiat0iSO5aVQ3oSHRM0yIyzrV/RcjnPuEyKk0rpD2X3/WrnTFOHtrUUIRf8BoA10LegvqIUvrx8TWx7ThRWk/pmdngBk+F6K5B9cfLACmE9LGPyf9/9ALSxz6mhNvW7HgzXNDdoBrY3UXbUfvZCjlyQ9ADBMp5EWQtm6cBwsXo3DQXhUQEvepazo+Xwb3v85gIpJTLb5H1kdYtAg34Ybv+fjm2fNU8ub6wLVsJbdWkyrpINOCVM1Fz+sH1w1Z5MBca781/5GfUf7kKvt++abEOQwycRo5uEgNKAplrz2ZI7lolCSsSJqmsCqIUUpvqDlWung/ekITMu19uF8VNBuNkoU2jngghz1BKHwHwcZzHThkcFeXIzpILwPA6I0L1smYP4TVwrFusDFCRiBZZ6sKrSpTiTalylA9k/SAAMYlUjsKliihe0winiLZPtPGo3fk2uLASqpKApzeBiiLSxz0WG4H0wVMAiJxAdd6VKrll3mRThNGA+CGz2u69UbP1n+CT7EouAKVUnvG3hvAWGqgICAZYBuSB+lxw/7gdwepShHY7oM/th/Sxc5RT5GpsBlRvXALBmg1j3yFwffdRguxmD0J+D0qfHatkszIjwWAcO60Nj70GQFOjMCzOYyc1kSS9ul3vynv4W5YpgzbRm0E4HrbrH1CFfxK9Baa+Q1RSGaZwKcugsxS8JS3u8w1frwcnxCb+6LLPVSVkReLZddnnxsSuJ5IYiWj5iO6amK0wx4Z8uIu2q7NOLXa49n4CX8leGE67ABqLXV7RhAKQ3HUx9aDjQ6DN7APzuVeh5vMCeSvMkq5S8bRdN1UOZ1yfD3/5/hiZB2333hDrKsLGpCphdrMmuRuy7EknTb1iBqOjadZQEEKmAJgK4IywSF8EC4D/tGfDOiOR+OVg9WGEGhyQ/B4lCxZURFpUPoZSEWx9PjxF22DsOwSg8r5/qMEhS0sbzTCcMxSu7z6K2R6yXTsVVBJjMr8Dh38Cb7Gpqoqlj52DqjULVMl7+tx+iihdTIapoEN14VKAcDFVzOwjZqBi1dPwlX6PUH0VfL99DVBJLvMJIBQeqK1X3RNTuAUcD0ghEEGv+Gto0Aeit8jZuuX7UeMokWP/qYRQXUV8rSyXE+C1cKxbDNFTK2fS8hp0mzAPgLoQTrzMeIGGMG8uq6TLYLQVLa0oVhKe1eAAACAASURBVAIohOzA/nPU4w2UUme7taqTEtm+uO3Ou+XaACMeRqjBIWcXJ5BkFl1OCN3OUCU5+cuKUFv4LLREkktwag2Kf4BPSgcohWN9viyJ4HUrEgucKQXmftfC/cOnsF4zBc4tryiFYuIVeAEQV5SOSrJD2LFusSJi6CvZC1/JHviK9wABL1zfFsbtA3/ZD+FV0uCox2QNJsJpZAnzusqwkfCDT86IFS8LeEF0cvRYc3pOKUNuh8Zih6NwKVKH3KHqVxr0g7fYoO12Jhxr5st1LRSH8XNsq4nBaENaUo+tA1BHCFkKwEkb1WOTCCG/p5S2vljwScLEiTdj+qxHwA2eqgxw5nOvRMmSPyScvYdqylUFYPS5/ZAy7GFUb3oeOTNXh+UqngxLIEiwj3g4SmFVLslItHpI7lp4ftyhyBtUb1ii6FLF813QkB8pl09VbWulXH4Lqjc8C3/5fnB6Ew69ei/E+srWdwCVgLBQXnQkFCQK67X3KtmtlavnKyJyMdtqu9cCFBA99bGGbEM+wGmQdOFI1Gz9JyS/F0m/Gx1jmAR7NjLvekku5PLsWMXvw2Aw2p7W+ihehrr+hCvOY6cMEae2CqKJL2IniUAofgGbUF1Fo56+xSoXQgkXpwEiCqsz4d28BA211fJleAGO9fmo2/UuiM4AQ+9BqC5cGuO7iCiRaix2lXPaV7wHmpTuaPhmo1ysx+dq+Ya5cGGWmnIkDbwBFJClHCIaSqEQkn4/VpFdkMX7PKBAzErKsXEJiEYHTtDBcMZFcO/7XFkxEUEHTm+G7WrZ4Bhy+6Ny1dNw7VFLpzSVtW5LtU4GgxFLaw0FoVFxtJRSiRByNDpRJxXRyrOAHPIKKkFy18orA04WDCNaAyCGEqvR2holgSNVw+IZlMraatjSusHp8sJ6zeRGZ/W/n1QJs0nu2nDUk0cWvxNDMXv41YVLkXzZRFRveDbxDRIOgrUHQi4n7HmzQAA4Cp8DANR/8QE4UwoIr0X6hDkq9d36/74HTmeC+YIRaPh6PXIeelet5x8WjSMaAcmXTYTGYof3l/8h5ap74Fifj+wH31GF/coRTD7YrpsGx4ZnFRE83mxV5L9PhFong3Gq09rB/ldCyJ8gryIA2cH9a/s0qfMTLcoVqWVMNAJAjOB0xhjhO401M0byumrdQpj6XKa8ZkSML55ByT29F6odVeAEPSrenQPBloXkQeORceMTqFz917BktBvg5XrRIJyskAoCcJwiryzYspRkNWi0QHRdBk4DY+9BMPUdAhBOFvjze+D48G+KOmp0wp67aLuybRTR+I/M9I19LkX9rvfgK96jPBcRAJS8DQCoHOZqy1Z0/Ws/W5HYXyGJACFIH/e4ql+tFiPLj2AwTgCtzcxOB/B3AFdCTrn9BMBDlNKj2NxuP05Uwl00SjZ3cQnSxz2Gyg+eBm9KSVDAfl5jtnNYTVOfcx6CVcVKUljlh38DAHCCTmVo6j9aittvuhGvvFGgrAx8B7+Do3AJdN37wPvzF+Foo/j1ni0DR8NzYGeMQ1v0uaG1ZSHp4hvhL/sRDd9slNVag36lGE1EetqxcQnEekdMop/r+0/h/PjluLLOVWufARWDoH4v+KQ0mM4ZCvcPn0J018kquI5idcb56vmyzHj0NtW6RTDygMvrU+WDRPq1ucL0DAajZdo8M7sz0xGGIgLH88h+eFVYd4jE11NaPAZ8UhrEBoeyGohkGhOOwJ6RiQaXG9pzr4V77ydy/YqgH0lWO176+3N49Im5qM8dDMnvhq9kL/zl+1otu20ZkAfPgV0gvKCsKiLlJiPJedG1EBJJhESH36oKPXFacHpDTAEeye8FOA4aU0pjsZigB8lmE/ih0+A58B+4i7aFK67JW1KmvoNlGXS/B7zehHvvvA0vvfiC0sdN+7X02bGQxPgGksFgtEybZGYTQmZTShcSQp5HHPEeSumfjqONJwURf0Vz4Z6alG7ysVHZzb7iPRCMZgTc9QAaVygNHifsGZmgkoTqqiPyquWXn4Ff4tfhTQRnsICKIfAWO0SXM64BowEvMibMV0dKJaieJfndimxGTIW9D/+mqrxGAz7wplSIAS9ETy1AKURPHWjQC5JkQd2GRUgeMQupV96trJpMWh7VezbF1WJq6hOK9CtzYjMYJwauhed/DP/+CsDuOD+nPPPmPg7PJy+A+j0Q/R44Ni6R1V7FkFIBL/nSmxCqPaJ63FG4FEGvHHEUCoVwxhmn4/Ntn+CtN9+CLySBHzoN2Q+vgnfgbeC0uhbbkZycDF5rQNKg8Ui96l5wOjNowIe6Xe+BaLTwlxWpjq/d+TaIVvZ5lL8+VXbIAxCs2THHRsrAGvsOQcPudUq93kjEVvoNf1amEbwxBbYR02HPmwlOq4c+tz/45HSkj52DnBmrwQ2eiqAEVLzzKEqXToBj1dO4/aYbUXXkECRRxG8/7YvxOUT6OLr/PJ+80G5lHxkMhhq29dQGFBSsxG133wvTBSPh/v4TSD43aNAH3mwFCIF9xMOoXDUPnMECsb4KGmsP6DLPgrb8W1z8u4H47LPP4HK5sGjRIrz4ymsxZQ6d25ej4b//hmpRR3gI9hyI9ZWYOP5GfPaf/8Ytj+jcsgz6nH5w798ZI/sR7SOIOKIDjpKY56rWLoTkqZO3z+qrEogSjkHOjNWq1+JNqTEZ40q7whLmrZXebiv5bwaD0Uib+CgIIevQjF40pXTUsTWvbeloQwHIEuQ13qDaaRwuEdrw7SZQvwfm/tdArK+S1VjDgoDRDB8+HIWbNsXdjy9ZPBqE14CKofD2jh89z+ytDJiJ9vEjgn7Oj5fJqq9+T4yiLRBVR4JKMJxxEXwle+SQVkEPKkkg4aijeDW7o2tURP/f/fa/x5Uij25X5LpMm4nBOPG01lC0tPW0GPj/9u49Osr62v/4e2dyDxAIAaRy1yCipVbRFrtUKNbLqTdsPVIv1Vq1cvFYrdjjwePhKNouEPFWa9XDT7sIWj0VFRVFK1aX1WOpWqqo4IUAViGA3HK/7N8fMxkymckkhMlMMvm81spaeZ65PPuJ8uyZ5/v97s0C4DOgCngg9LMH+GR/g0wn28q/pO9xF7L9pfvYcNtZbHniZhp2bmbXm4/j1bugsY497zxH1Sd/jZkkAF577TWGjTw45q2fQG4BA8+5iWHXPsnAs2+k94DBEZ+qiwd9Lebrmnpb548+Fgtk0WfCv+K1scch6nd8Sf7YifT77qUUnXgFgd7FeGMjg374XxSdNJ2tzy0kr2QCW5ffGXUbrXDCuRHvVbdtU8SK8ZZxNa0haTpuWdkGSkuXdOyPn0SlpUsYWTKGjECAkSVjukXMIvurrRIefwYwswUtss4yM0vtR/guZvioEqp6FzP4ojv4/L6f0lDxVfCBdpaWGDhwIBMnTmTSpO9y3Y03RTRO3/bMbRR887SIFdtMnskNc27i/PPPo7R0CTt27aaxZTmMUPmPDQumkFl4AI01FfT9zo/Yveqp2GsWsvOo/OBVdr35eLDOVE0lfY4+K/y82i/WsvvtZXhN1d7B65w8eh95elSJjczCQZQvm0/O0MOjyqi3XFmd1X8oRd+7Irg2BbrsLaXm3RCHnjWWqk1runzMIonQ3nUUHwDfd/dPQ9sjgefc/dBOjq9dkn3ryd3Zvn07/fv3D+8rLV3CxVdcSf/TrmX7ygep2xx/PWJRUVEoMUxi0qRJjB07FjMLv1fz+/HrP1nLsF9E375pmh46smQMVeMvoqHiq3BviczCQdTv2U7/k2cEC+s9u4CGyt0MOmcOW5+5HTIsam1F3qijqV7/DvU7vgzOWqr4KqI731d/fijqNY21NZgZA878ZcSYhnsjuUPGUrvxHzRUV5KV34u6qj1kZOfTq9kajaak0bT+oiuvjWj6O2s9h6SLRLdCPQW4n+BqbAOGAz9z9xf2N9BESEaiWL9+PStXruTll1/mlVdeIS8vj7Vr10Zc1Hv1KqSippbGmoqo11tOPtnFw6ktX08AZ9GDD3DhhRe069htXaBaHZ+4bUr42N7oUFdFRkEhjRU7ycgvxLJyaNhVTiC/Lw01FVBfi2Xnkj34EA6YOrfNft9N4wvujWRk5exNULu3BkuXhJLDiIP2jqU0bw+a1X8ohceeG1GuvCuvjdB6Dkk3Ce1w5+7Pm1kJMCa060N3r9mfALu6jRs3snLlSlauXMkrr7zC+vXro55z1113c8PcX4dvRdRsWkPds/PJzs9l11fbICODnKFfp9/En5A9cCSWEQhXVr3xprntThTNS4Y0fWpvXuOotXUGmX0PCK/83vrcQhoa6+k17mR2vfEYQ6Y/FNHHu8+4s8Irx2u/XMvmP/wnhRPODZfpqNu2qdXChuAcGJrxtPW5hWTkFEBDLZ7dC6oq2LR1F5dOmxl+XUFBAVuzcin63hXdam2E1nNIT9XeVqj5wDXAcHe/zMxKzOwQd3+mc8NLni+++CKcGFauXMknn7Q9Vn/zLbeSP/nfIscOvj+L3L8+xLYtX5KVlc2gc/47utBdTSVln7Z/AV3T/e8b5txE2WProtp7zp1zI5ddeQ2cfFXEbaF+J1wcXuvQ1AJ115uPY9l54TGKnW/8Idh7es0rUbWovLEBr69j6zMLwmsxopJR4SAaqnaxYcHZwYHzulr6ffenbHt2IQOm3BARz7TpMyEnj/zJMyk6ZGtUifGuXuCvrYQtkq7aWxTw/xFcYDchtP058DjQ7RPFokWLmDdvHh999NE+vS43N5et5ZsZ1uxTdvA+/sM07ConO78PlpnNjtcfod/xF4af01TobvjQIbHetlXnn39eqwOmTfuvnvVLNnzxOZaVQ9FJ0yMGmJua/Qy7dinlT/06XJK8btsm8DeiOt0NOP06tvzx5nAdp4KvnxizCVJjTSUWyMKpwjKzKDrucgIF/cjse0B0r+4nbmbgv8wK77eMQLBg4Y4vGXHw6C7f27qthC2SrtqbKA5y93PN7EcA7l5pTSOv3VxdXV27kkROTg4TJkwIDz4fc8wxjDn8G+FbEa0N9u5+O5hLmwZvy5fNJ9NI+KripkQysmQMZRs2xOwlbdm5WCCTgjHHUfP5B2x+dDaWnUvd9o2xe2vXV5ORk0fd9o0M/s5dZBcPC3ao27kltL6inrxRR1G35TMGnDE3YjC736RLot+vtjriOAVjTyD/kO+w8fazu81gcLyELZKu2psoas0sj9DiOzM7COjwGIWZnQPMAQ4FjnH3Vc0eux74KdAA/FtnD5hPnDgx5v6srCy+9a1vhRPDt7/9bfLy8iKe0/xWxM6//CG66dC/XM225+9m99+WhVqe5tIrN5ff/u7+TrvYlH26jkDR0Kg+FOVPz8Nrq9l078V4Qz0Dztg7S2nL0rmxp8tm5tCroIDdjQ3UbFpDwdgTwt9SqstWU7ViIRX//IiCcSex/cX7QpVx84BAzEQViNGFT/f4Rbq+9iaK/wKeB4aaWSnwHeDi/Tjue8DZwO+a7zSzscBU4DDga8BLZjba3TttSsno0aMZPHgwW7Zs4eijjw4nhmOPPZaCgoK4r21+K6JuW+xP5fU7N2NG0lp1Dh9VwoZ/bsa9MdyHIrNwUPD4gUwskEXB4ZMj2pPmDhtH+bJ5UZ3oik6aTmbvYmqW30750lsZMOU/Iu7N//auO3j9L3/hwYd+T13lHrLye5FpRta4U6P6b2xdNp/LL7mIxY/fo3v8It2Nu8f9ITgddijQH/g+cBpQ3Nbr2vMDvAKMb7Z9PXB9s+0XgAltvc9RRx3l++Pdd9/1Xbt27dd7jDj4EB809VYf/stnwj+Dpt7qmX0H+4iDD+nw+y5eXOojDj7ELSPDRxx8iC9eXNrm8y0rN2YslpXrYB4oHOSDpt7qw6590gdNvdUDhYMczEccfIiDeWbfwV58+qyI1w4YPCQijmnTZ3jxoK+5ZeU5mA8YPMSnTZ/hub37eaDPAO8z4VzP6j/MMfOMnAKfNn1Gh85HRDoPsMrbca1u7zqKf7j71xOdpMzsFeBaD916MrN7gDfdfXFo+3+A5e7+vzFeezlwOcCwYcOOKisrS3R4+6S0dAmXXXkNfVrMPGpPwbt479m0Ejjik3wbA6hmGWT1H0rdto1k5BbQWFNBVtHQUDvS2HWetjxxM401le1aK1BauoRLp82kzjIjx2SWzadg3ElkFw+LWPjXN8co//LzfT7/9v6NVCxQpGMSuo4CeNvMjnb3v+5DAC8BB8R4aLa7P9Xe92mNu99PcBEg48ePT3kJ3OiZR7n0L+rHHQvu6PCF64Y5N5E/eWarpTtimT5jJhn5hWQWfY2G6t2RVWCXzaMx1GGuuaaBZmjfWoEb5txEfSCX4mbd/HKHj6P49Flsf/E++h1/YdQius6gkhoiydHeRPEt4AIzWw9UELwd5e4+rrUXuPuJHYjnc4K3uZoMCe3rFhI9I6bs03UMPSv6ol72WOtrMO5f9BC9vnkau99+hoFTZkdPeX0i9sD1gMEHUlq6hIqKCsofnR3R9a7lOELZp+twJ2bCqdu+MWJfZw5WdySRisi+a6t6bJOTgVEEe2afTnCc4vROiOdpYKqZ5YTqSZUAb3XCcbqF4aNKYlZejXfhbaiupGrtG3htZSvfHKoof3peVHOlH045k2lXzyLjhOkMu3Yp/U+5kq9WLqJhZfStruGjSsgsHBS7ym1OftIaDJV9ui7mOe7LYkYRaVvcRGFmuWb2c2AWcArwubuXNf109KBmNsXMNhFcwPesmb0A4O7vA48BawjOsprhnTjjqTMksgx1Rzq7WVYOdds3Eug9IOaF3HLyKJp8WbAc+oKz2fb83RT1zmf5ipfCn86bVnMPOOM6tu/YwQUXXEB2QR+mzwiW4Tj1pBNprNrJ5kdn8/nvLmPPey9TXbaaXS/cyeWXXEzeqofZePvZ5K16uM3xlP3RkUQqIvuurVtPDwN1wGvAqcBY4Kr9Pai7LwWWtvLYLcAt+3uMVEj0PfOOrATuVZDP7j1OQ9Xu8OrrvWMU8/GaKgIF/Rh88V3hwfE7Fsznwh9fGPM2V2NtFcOuDdZxuv/h+axdt5a33n0/ojxH+dPzKMjO4IF77k7qLR+V1BBJjrY63IVnO5lZJvCWux+ZrODaqyt0uINgldftA46gat0b4TUKeSUTKCp/N2krj0tLl3DhTy4lUNCX/EOPj4ql4p1nGTrkwKhZQq1VqG3Zua586VwGTImeNZWqUtua9STScYlqhfp288TQcrur6CqJwjIyCPQZGNGkZ+vyO2nYtQVvTM6COwhOj8WMYb+IXXp88eLF4Ytp04V2/SdrySzoR//Tro2IvV+oV0Tz18dqbapS2yLdT6Kmx37DzHY1vSeQF9pumvXUZz/jTCtZeb0oalFcr/jUq9j+9K+SGodl5xHIL4xd7bXvAfz48hn8/BezOOcHU1j8+FLyJ89k2Flj2fH6I2x98hYaa6qC5UZidK7LyMlTGQ6RHibuYLa7B9y9T+int7tnNvtdSaKF+qqK2GU8qqIbGXUmr62ioWoXW59bGDm7adl8+h53AQPOuI4dNc4Dv38EO/i48AB2v+MvpPis2Qw/qISf/fQnVKxeEdkbe9l8Jh3/nX0eYBeR7q296yikHYYf1MpitYOS+2nbsnLpfeRpVLz3p3Cp8ECvIhqrKygYewLeUE/9zs0MOndueIFck6Z1Gvf+Jjje8OBDvwrXcbr84h9z72/u2TsuoFLbIj1Cu0p4dHVdZYyioyU3Es0sI+Y4wobbpjD8l8uC3w6W30lGZg512zaSVTyUwgnndou+1SKSOIku4SHt0FUa21hWTuyy4Vk5wSTx3EK8rpaiFoPutVs3UPveChp7FZARCGgWkYgA+kaRliwzh0BBYVQTpYaKnWT27o/X11F82jVRU1y3PnkL2Xm9IgobpuIbkYgkR3u/UbS3hId0Jw11NNZUsu35u8OrrxtrKqGhjiHFfWjYsy3moHtDdSV9Tr4qYnV2fqh2koj0XEoU6SiQhddU01C5E9xpqNyJ11RDIIvP1n3IiINHxy7vkZWr2kkiEkWJIg0VFxdjufkE8gvBjEB+IZabj5lhlsHGDRvYsfz2qCmu/Yv6qXaSiERRokhD5/xgCrRsveqN5I78JsOuXUrx2f9JbXUVVSsWRhTvu2PBfK2REJEomvWUROH1B5+sIzOvgLqqPYw4aHTCZxYtX/FSRC8K2Fu3qWnsYeCU2Wx/+lcxy26ketaWiHQtShRJ0rKybNOU1O0Djkh4V7bWGh7VbdsUuV25J+q1iW6+JCLdn249JUnzbmxNn+qLT72KqnVvJHxmUWt9GrL6D4nczu+VsGOKSPpSokiS1rqx1W3blPCZRbEaHpUvm09eyYSIuk2XXvzjhB1TRNKXbj0lyfBRsetAZfUfkpCZRS37MlxwzhSWr3g4PNYw6djxvPr6c+x647GIuk0iIm1RokiSWN3Yti6/k4KxE/e7K1usznqLH9eKahFJDJXwSKLOmvXUWnc6FfcTkXgS0uGuu+guiaKzZAQCDL0mupudus6JSDyq9dSDtDbLSSuqRSQRlCjSQKxZTlpRLSKJokTRiUpLlzCyZAwZgQAjS8ZQWrqkU45z/vnn8duF88lb9XBESQ4NZItIIihR7KOWF//pM2bGTAZNM5Gqxl/E0GueoGr8RUy7elanJovP1n1IY0MDn637UElCRBJGg9n7oGWr0x2vP0LF6hUUnz4rqtHPDXNu0kwkEenSNJjdCVqW4aha9wbFp8+K2eintZXY6u0gIt2NEsU+aHnxbyq/0VxTMtBMJBFJF0oU+6Dlxb+p/EZzTclAM5FEJF0oUeyDlhf/vJIJbF02P2Yy0EwkEUkXGszeRy2L75160oksX/FSeDvRTYhERDqLSniIiEhcXXrWk5mdY2bvm1mjmY1vtn+EmVWZ2buhn/tSEZ+IiOyVqjLj7wFnA7+L8dgn7n5EkuMREZFWpCRRuPsHAGaWisOLiMg+6Iqznkaa2Ttm9mczO661J5nZ5Wa2ysxWlZeXJzM+EZEepdO+UZjZS8ABMR6a7e5PtfKyL4Bh7r7NzI4CnjSzw9x9V8snuvv9wP0QHMxOVNwiIhKp0xKFu5/YgdfUADWh3/9mZp8AowFNaRIRSZEudevJzAaYWSD0+yigBPg0tVGJiPRsqZoeO8XMNgETgGfN7IXQQ8cDq83sXeB/gSvcfXsqYhQRkaBUzXpaCiyNsf+PwB+TH5GIiLSmS916EhGRrkeJQkRE4lKiEBGRuJQoREQkLiUKERGJS4lCRETiUqIQEZG4lChERCQuJQoREYlLiUJEROJSohARkbiUKEREJC4lChERiUuJQkRE4lKiEBGRuJQoEqy0dAkjS8aQEQgwsmQMpaVLUh2SiMh+SUnjonRVWrqEaVfPIn/yTIaeNZaqTWuYdvUsAM4//7wURyci0jHm7qmOYb+NHz/eV61aleowGFkyhqrxF5E7fFx4X3XZavJWPcxn6z5MYWQiItHM7G/uPr6t5+nWUwKVfbqOnCFjI/blDBlL2afrUhSRiMj+U6JIoOGjSqjZtCZiX82mNQwfVZKiiERE9p8SRQLNnXMjlX+6h+qy1XhDPdVlq6n80z3MnXNjqkMTEekwDWYnUNOA9Q1zbqLssXUMH1XCgoXzNZAtIt2aBrNFRHooDWaLiEhCKFGIiEhcShQiIhKXEoWIiMSlRCEiInEpUYiISFxKFCIiEpcShYiIxJWSRGFm883sQzNbbWZLzaxvs8euN7OPzewjMzs5FfGJiMheqfpG8SJwuLuPA9YC1wOY2VhgKnAYcApwr5kFUhSjiIiQokTh7ivcvT60+SYwJPT7mcCj7l7j7p8BHwPHpCJGEREJ6gpjFJcAy0O/HwhsbPbYptC+KGZ2uZmtMrNV5eXlnRxi+6gNqoiko06rHmtmLwEHxHhotrs/FXrObKAeKN3X93f3+4H7IVgUcD9CTQi1QRWRdJWy6rFmdjHwM2Cyu1eG9l0P4O6/Cm2/AMxx9zfivVdXqB6rNqgi0t106eqxZnYKcB1wRlOSCHkamGpmOWY2EigB3kpFjPtKbVBFJF2laoziHqA38KKZvWtm9wG4+/vAY8Aa4Hlghrs3pCjGfaI2qCKSrlLS4c7dD47z2C3ALUkMJyHmzrkxOCYxeSY5Q8ZSs2kNlX+6hwUL56c6NBGR/aJWqAmiNqgikq7UClVEpIfq0oPZIiLSfShRiIhIXEoUIiISlxKFiIjEpUQhIiJxpcWsJzMrB8qa7SoGtqYonK6gJ5+/zr1n0rl3zHB3H9DWk9IiUbRkZqvaM+UrXfXk89e569x7mmScu249iYhIXEoUIiISV7omivtTHUCK9eTz17n3TDr3TpSWYxQiIpI46fqNQkREEiStEoWZnWNm75tZo5mNb/HY9Wb2sZl9ZGYnpyrGZDCzI8zszVCvj1VmdkyqY0omM7vSzD4M/b8wL9XxpIKZ/cLM3MyKUx1LspjZ/NB/99VmttTM+qY6ps5mZqeErmkfm9m/d9Zx0ipRAO8BZwOvNt9pZmOBqcBhwCnAvWYWSH54STMP+G93PwK4MbTdI5jZJOBM4BvufhhwW4pDSjozGwqcBGxIdSxJ9iJwuLuPA9YC16c4nk4Vuob9BjgVGAv8KHStS7i0ShTu/oG7fxTjoTOBR929xt0/Az4G0vlTtgN9Qr8XAv9MYSzJNg34tbvXALj7lhTHkwoLCbYa7lEDkO6+wt3rQ5tvAkNSGU8SHAN87O6funst8CjBa13CpVWiiONAYGOz7U2hfenq58B8M9tI8BN1Wn+yamE0cJyZ/Z+Z/dnMjk51QMlkZmcCn7v731MdS4pdAixPdRCdLGnXtW7X4c7MXgIOiPHQbHd/KtnxpEq8vwMwGbja3f9oZv8K/A9wYjLj60xtnHsmUAR8GzgaeMzMRnkaTe9r4/z/g+Btp7TUnn//ZjYbqAdKz0oRvwAAA4lJREFUkxlbOut2icLdO3LB+xwY2mx7SGhftxXv72BmvweuCm0+DjyYlKCSpI1znwY8EUoMb5lZI8FaOOXJiq+ztXb+ZvZ1YCTwdzOD4P/nb5vZMe7+ZRJD7DRt/fs3s4uB04DJ6fThoBVJu671lFtPTwNTzSzHzEYCJcBbKY6pM/0TOCH0+3eBdSmMJdmeBCYBmNloIJseUizO3f/h7gPdfYS7jyB4K+LIdEkSbTGzUwiOzZzh7pWpjicJ/gqUmNlIM8smOGHn6c44ULf7RhGPmU0B7gYGAM+a2bvufrK7v29mjwFrCH4lneHuDamMtZNdBtxpZplANXB5iuNJpkXAIjN7D6gFLuoBnywl6B4gB3gx9I3qTXe/IrUhdR53rzezmcALQABY5O7vd8axtDJbRETi6im3nkREpIOUKEREJC4lChERiUuJQkRE4lKiEBGRuJQopMcLVVld3Gw708zKzeyZVMbVFjN7pWWVZJHOoEQhAhXA4WaWF9r+HilauR9a+yLSpShRiAQ9B3w/9PuPgEeaHjCzAjNbZGZvmdk7ocJ7mNkIM3vNzN4O/Rwb2j/YzF4N9QN5z8yOC+3f0+w9f2hmD4V+f8jM7jOz/wPmxTlenpk9amYfmNlSoCmxiXQqfXoRCXoUuDF0u2kcwRXex4Uemw287O6XhJrhvBUqTrcF+J67V5tZCcHkMh44D3jB3W8J9QzIb8fxhwDHunuDmd3ayvF+BlS6+6FmNg54O2FnLxKHEoUI4O6rzWwEwW8Tz7V4+CTgDDO7NrSdCwwjWFPrHjM7AmggWOIcgjV4FplZFvCku7/bjhAeb1ZWprXjHQ/c1Sze1ft2liIdo0QhstfTBPt3TAT6N9tvwA9aNsUysznAZuAbBG/jVgO4+6tmdjzBW1kPmdnt7v57IhsJ5bY4dkU7jtexsxLZTxqjENlrEcEWsv9osf8F4EoLXanN7Juh/YXAF+7eCFxIsDAbZjYc2OzuDxAs8X5k6PmbzexQM8sApsSJo7XjvUrwthZmdjjBW2QinU6JQiTE3Te5+10xHroZyAJWm9n7oW2Ae4GLzOzvwBj2fiuYSLAnxDvAucCdof3/DjwD/AX4Ik4orR3vt0AvM/sAuAn42z6fpEgHqHqsiIjEpW8UIiISlxKFiIjEpUQhIiJxKVGIiEhcShQiIhKXEoWIiMSlRCEiInEpUYiISFz/H2pBkqI1wRffAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from sklearn.model_selection import cross_val_predict\n", "\n", "predicted = cross_val_predict(model, X, y, cv=10)\n", "\n", "fig, ax = plt.subplots()\n", "ax.scatter(y, predicted, edgecolors=(0, 0, 0))\n", "ax.plot([y.min(), y.max()], [y.min(), y.max()], 'k--', lw=4)\n", "ax.set_xlabel('Measured')\n", "ax.set_ylabel('Predicted')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "# lest do the cross validation by hand\n", "import sklearn" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "X_train, X_test, y_train, y_test = sklearn.model_selection.train_test_split(X, y, test_size=0.1)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model = LinearRegression()\n", "model.fit(X_train, y_train)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.5827284267819932" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y_predict = model.predict(X_test)\n", "\n", "# calculate mean square error\n", "sklearn.metrics.mean_squared_error(y_test, y_predict)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.6.5" } }, "nbformat": 4, "nbformat_minor": 2 }