{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Initial FeO prediction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Start by importing MagmaPEC and MagmaPandas (or regular Pandas)" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import MagmaPEC as mpc\n", "import MagmaPandas as mp" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Import your melt or whole-rock data" ] }, { "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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
SiO2TiO2Al2O3MnOMgOCaONa2OK2OP2O5Cr2O3FeOtotal
name
PI03246.9628183.44738116.6501240.1746187.0045589.4007413.1312531.2905360.6353460.00963011.350806100.057812
PI04146.2927093.55397317.1918740.1706716.2216999.4434162.7011101.0760790.6103870.00726210.98633698.255517
PI05347.5157363.08409615.1822940.1699488.2505779.6403572.2982231.0133320.4480630.00922510.95543898.567288
PI05447.5243573.36959416.4565870.1746716.5696229.2163092.6539941.3673370.6485340.02954010.95366098.964204
PI05546.4521683.28701816.0398060.1690847.3960609.4698992.0352041.0722110.5397920.01167610.89614797.369065
\n", "
" ], "text/plain": [ " SiO2 TiO2 Al2O3 MnO MgO CaO Na2O \\\n", "name \n", "PI032 46.962818 3.447381 16.650124 0.174618 7.004558 9.400741 3.131253 \n", "PI041 46.292709 3.553973 17.191874 0.170671 6.221699 9.443416 2.701110 \n", "PI053 47.515736 3.084096 15.182294 0.169948 8.250577 9.640357 2.298223 \n", "PI054 47.524357 3.369594 16.456587 0.174671 6.569622 9.216309 2.653994 \n", "PI055 46.452168 3.287018 16.039806 0.169084 7.396060 9.469899 2.035204 \n", "\n", " K2O P2O5 Cr2O3 FeO total \n", "name \n", "PI032 1.290536 0.635346 0.009630 11.350806 100.057812 \n", "PI041 1.076079 0.610387 0.007262 10.986336 98.255517 \n", "PI053 1.013332 0.448063 0.009225 10.955438 98.567288 \n", "PI054 1.367337 0.648534 0.029540 10.953660 98.964204 \n", "PI055 1.072211 0.539792 0.011676 10.896147 97.369065 " ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "wholerock_file = \"./data/wholerock.csv\"\n", "\n", "wholerock = mp.read_melt(wholerock_file, index_col=[\"name\"])\n", "wholerock.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create the FeOi prediction object and initialise it with your data. Make sure to remove the FeO column from the melt compositions used to predict FeOi. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "x = wholerock.drop(columns=[\"FeO\"])\n", "FeOi_predict = mpc.FeOi_prediction(x=x, FeO=wholerock[\"FeO\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Select the columns in *x* that you do not want to use to predict initial FeO contents. Here we do not want to use minor elements and the totals column:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "do_not_use = [\"MnO\", \"P2O5\", \"Cr2O3\", \"total\"]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The *calculate_model_fits* method calculates best-fit multiple linear regressions for all possible combinations of elements in *x*. For each new regression, the element whose removal results in the lowest regression F-test p-value is removed from the dataset.\n", "\n", "It returns a dataframe with fitted coefficients and misfit statistics (RMSE, cross-validated RMSE and R2). *NaN* values for coefficients indicate that this element has been removed from the calibration dataset and is not used in the regression." ] }, { "cell_type": "code", "execution_count": 5, "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", "
interceptSiO2TiO2Al2O3MgOCaONa2OK2ORMSECV-RMSEdeltaRMSEr2
624.162353-0.1844501.085716-0.337465NaN-0.2922310.344686-0.5795780.1514160.2781020.1266860.930966
521.878478-0.1447301.179186-0.293183NaN-0.296750NaN-0.3007630.1653870.2735530.1081660.917640
422.442147-0.1681071.114827-0.289476NaN-0.260280NaNNaN0.1710320.2323490.0613170.911921
311.585125NaN1.463365-0.230722NaN-0.169262NaNNaN0.2064340.2375330.0311000.871685
29.213242NaN1.640962-0.217524NaNNaNNaNNaN0.2537690.2947480.0409780.806092
15.357205NaN1.725906NaNNaNNaNNaNNaN0.3194300.3556700.0362400.692766
\n", "
" ], "text/plain": [ " intercept SiO2 TiO2 Al2O3 MgO CaO Na2O K2O \\\n", "6 24.162353 -0.184450 1.085716 -0.337465 NaN -0.292231 0.344686 -0.579578 \n", "5 21.878478 -0.144730 1.179186 -0.293183 NaN -0.296750 NaN -0.300763 \n", "4 22.442147 -0.168107 1.114827 -0.289476 NaN -0.260280 NaN NaN \n", "3 11.585125 NaN 1.463365 -0.230722 NaN -0.169262 NaN NaN \n", "2 9.213242 NaN 1.640962 -0.217524 NaN NaN NaN NaN \n", "1 5.357205 NaN 1.725906 NaN NaN NaN NaN NaN \n", "\n", " RMSE CV-RMSE deltaRMSE r2 \n", "6 0.151416 0.278102 0.126686 0.930966 \n", "5 0.165387 0.273553 0.108166 0.917640 \n", "4 0.171032 0.232349 0.061317 0.911921 \n", "3 0.206434 0.237533 0.031100 0.871685 \n", "2 0.253769 0.294748 0.040978 0.806092 \n", "1 0.319430 0.355670 0.036240 0.692766 " ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model_fits = FeOi_predict.calculate_model_fits(exclude=do_not_use)\n", "model_fits" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "RMSE and R2 are both calculated on the entire calibration dataset and indicate how good the model is at predicting FeO. For RMSE lower values are better, while for R2 values close to 1 are best. To check for overfitting, cross-validated RMSE's (CV-RMSE) are also calculated, where large differences between RMSE and CV-RMSE (deltaRMSE) can indicate overfitting issues. As long as RMSE and R2 values are acceptable, the model with the smallest deltaRMSE should be selected. Here, models 1, 2 and 3 have similar RMSE, deltaRMSE and R2 values and any of these models would work well for predicting melt FeO contents. In models 4, 5 and 6, deltaRMSE has increased and overfitting might be an issue here." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next we use the results for the previous step to select our model. Here we use model 3, where TiO2, Al2O3 and CaO are used as predictors. In the *select_predictor* method, you pass the index of your preferred model to the *idx* parameter. It shows a plot of the regression and calibration data. Be mindful of the FeO range, and make sure that you do not extrapolate (too much) beyond this range during PEC corrections." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAAFaCAYAAADVfgw3AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAARmpJREFUeJzt3XlYE+f6N/BvWATcQFEgskhldUEQsIIiIC5V3EBxx13rEa3SVi0/rbX0VKvt61Fr3etua10qoiCuFRS3ulFLPVZREdG4YRsQ2fO8f3BICRBIgGRmwv25Lq/LhJnMfSf4zfjMzDMixhgDIYQQQdDjugBCCCGqo9AmhBABodAmhBABodAmhBABodAmhBABodAmhBABodAmhBABodAmhBABMeC6AE2RyWR4+vQpmjVrBpFIxHU5hBAixxhDTk4O2rRpAz099faddTa0nz59CltbW67LIIQQpR4/fgwbGxu11tHZ0G7WrBmA0jelefPmHFejPsYY8vLyYGJiItj/KVAP/CD0HoReP1C5h+zsbNja2spzSh06G9plH27z5s0FG9qGhoaC/0WlHrgn9B6EXj+gvIfa9KOzoV2VoqIilJSUcF2GShhjKCgogEgkEvQvasUe9PX1YWhoyHFlhAhXgwntoqIi3L9/H0Ka1FAmk6l9kIJvKvYgEong4OBAwU1ILTWY0C4pKQFjDG3atIGRkRHX5dSIMSYPPCHvaZfvoaCgAE+fPkVJSQmFNiG11GBCu4yRkRGMjY25LqNGuhjahJC6E/b/vYlSM2bM4LoEQogGUGhriUwm08o6ZTZt2lTrdQkh/EWhrUGJiYkYNGgQhgwZgt27d6Nnz57o3r079u7dCwC4efMmvL29MWTIEAwdOhSJiYnydYYOHYo9e/bg+PHjldb79NNP0b17d/j7++Py5cu4dOkSunXrhoCAAHz22WcAAG9vbwCl56kHBQWhZ8+emDlzJgBgx44dGDZsGAYOHIiuXbvi6dOnHLw7hJDaaHBj2uVJJBJIJBKF51q0aIF33nkH+fn5uH37dqV1PD09AQB//vknmjdvDrFYXO02srOzcfbsWfTs2ROJiYkwMDBAr169MHLkSHz66af48ccf4eTkhICAAIV1EhMTUVJSgoCAgErrnThxApcuXYKBgQFkMhmWLFmCxYsXY9CgQZX2zpcvX44FCxagf//+mDp1KpKSkgAAZmZm2LZtG7Zs2YKDBw9izpw5tXoPCWnoJBIJoiKnY/nqLTXmQX1o0HvamzZtgpeXl8KfxYsXAwAyMzMr/czLy0u+7qRJk1QagvD29kZWVhbu3buHfv36ISgoCK9evcLLly/x4sULODs7QyQSKby2t7c3RCIRXr16VeV6X375JWbMmIEZM2bgxYsXmDVrFk6dOoUJEybg+PHjCtu/f/8+unbtCgDo2rUr0tLSAABdunQBANja2uKvv/6q2xtJSAMlkUgQHhqAMPExhIcGVNoJ1IQGvac9Y8YMDBkyROG5Fi1aAABsbGxw/fp1pevu2LFDpSst9fT00KpVK7Rv3x6nTp2CoaEhioqKYGhoCEtLS9y7dw+Ojo64ceMGhg4dKl8HgNL1TE1N8d577+HHH3/E5s2bMX/+fKxZswaFhYXw8vJCcHCwfPuOjo64evUq+vfvj6tXr2LixIl48OCBwtkcQjp3nRC+KAvsdUPT4GwlgpNlGsJDA7AnJkmje9wNOrTFYrHSN9fY2Fg+FFIVFxcXlbejp6eHRYsWoU+fPtDT00Pr1q2xf/9+fPHFFxgzZgysrKzQtGlTeTCXX2/hwoWV1gsJCUFeXh4KCgrw/fffY9OmTTh06BByc3MxadIkhW1/8sknmDhxIpYuXYpOnTrB398fDx48ULl2QkjVoiKnI9KnNLABwNlKhEifNERFTsfOfXEa266I6ehuVnZ2NkxNTSGVStG8eXPk5+fj4cOHeOedd3hznnbZnrNMJkNQUBD27t0r/xLRhXOcK/bAx8+gJro4WZHQ8LX+invad58xzIp1rHJPu6oJo8rnkzoa9Jg2165cuYKAgAB069YNffr00cpBDEJI/RCLxdgTk4RZsY6I/015YNe3Bj08wjU/Pz/52RyEEOEpC+6oyOnYE6Ods0cotAkhpA7EYrFGx7ArouERQggREAptQggREAptQggREAptJSQSCSaOGqTxK5w+//xzxMXF4dmzZ1iyZAkAIDAwEG/evFH7tXbs2IHCwkL53y9dulQvNfr4+Cj9WWJiIu7evVsv2yGE1IxCuwpcXJpqZWWF6OholZZVNvtf+dCeNGkSfH19660+ZSi0CdEuCu0Kyp8wP9BdhHVD02oV3IwxzJo1Cz179kRAQABev36N3bt3IygoCJ6enti9e7fC8unp6QgLC5M/jo6ORq9evfDBBx8AKA3kUaNGYeDAgTh9+jQ+/vhjBAYG4t1330VKSgouXbqElJQUDBgwAGvWrJHvwQPAhx9+CD8/P/Tq1QsPHz4EALRv3x7jxo1Dly5dKtUCAMuWLYOvry9mz54tv69mxfrz8vKwY8cO/N///R8mT56M58+fo3fv3vD390dYWJhg7sdJiKAwHSWVShkAJpVKGWOM5eXlsdu3b7O8vLxq15swciA7MlfESrbpyf8cmStiE0YOVGv7hw8fZnPmzJE/LikpYbm5uYwxxvLz85m3tzdjjLElS5awo0ePsocPH7Lhw4czxhgLCAhgx48fZ8XFxWzs2LHs2rVrbPv27Sw8PFz+emWv9dtvv7GxY8fK18vJyVF43V9//ZWNGjWKMcbYuXPn2KRJkxhjjJmZmTGpVMpycnKYj4+PQu0SiYT5+/szmUzG/vzzT/bOO+8obLOq+hljrKCggBUVFTHGGPvwww/ZiRMnWHFxMZPJZIwx1T8DPpHJZCw3N1fegxAJvQeh189Y5R4q5pM66DztCpav3oLw0AA4Wf5zaerqy47YE7NFrde5c+cOevbsKX+sp6eHU6dO4T//+Q8A1DikUDbrX/mZ+cpm6wOAlStX4sSJE9DT04O+vr7S16k4y9/ChQsBAO3atZNfPssqzGSQnp6Ozp07QyQSwdnZGaampgBQY/2vX7/Gv/71L/z111+QSCTo3LlztT0SQtRHwyMV1Nelqe3bt8eFCxfkjxljiI6OxpEjR3DixAk0bdq02vVv3rwJALh+/TocHR0B/DP7X1ZWFuLi4nD+/Hl899138tA1NDSsNCRRNssfAFy9ehVOTk4AUO0cDvb29khNTQVjDGlpaZBKpQBQZf3lt/nDDz+gX79+SEpKwqBBg2j2QEI0gPa0q1Afl6YOHjwYCQkJ8PPzQ6NGjXDgwAEMHz4cvXr1goeHh3wKWGUSEhLwxRdfwN3dHV5eXvj999/lP2vRogUsLS3Rq1cvhTM7hgwZgpEjR2LkyJHy57y9vSEWi+Hn5wcDAwNs3769xtqtrKzQt29f+Pr6wtPTE+bm5gBQZf1BQUH45JNPkJSUhAkTJmD8+PE4ceIEGjduDDc3N7XeM0JIzbQ6y9+cOXNw5MgRPHr0CL///js6deoEAJgyZQouXLgAExMTNG/eHN9++y08PDyqfI24uDjMmzcPxcXFcHd3x86dO6vcaxXCLH/VYTTLHy8wns4wpw6h9yD0+gEBz/IXFhaG5ORktG3bVuH5kJAQ/PHHH0hJScGCBQsU9hTLe/PmDaZOnYrDhw8jLS0NYrEYS5cu1UbphBDCC1odHvH396/y+fJ3j/Hx8cGjR4/ke2jlJSQkwNvbG66urgCAiIgIBAcH46uvvlK6TalUCsYYCgsLIZPJIJPJBDfWKrR6q8IYk/dR/u98V75moRJ6D0KvH6jcQ1164d2Y9po1axAcHFwpsAEgIyNDYS/d3t4eT548qTLgy9jZ2QEA2rZti40bN6Jp06awsbHRTPEaoOxCGiEp66HsSzM/P19w/wDz8/O5LqHOhN6D0OsH/ukhLy+v1q/Bq9Des2cP9u/fj/PnzytdRt0xrYyMDDRv3hyFhYV48eKF/O9CUN2XkVCU76GoqAh6enowNjYW1Jh2fn4+jI2NBT2eKuQe+Fy/VCqVnxJbnYo9lL+toLp4E9r79u1DdHQ0zpw5AwsLiyqXsbOzwy+//CJ/nJ6eDmtr62qDzdTUFM2bN0dRURGysrLw7Nmzeq9dU3QttIHSL10DAwPe/eOriUgkElzNFQm9B77Vv2XLFixcuBDXrl2rdJxOmbIe6tIHL0J7//79+PTTT3H69Gn5cEZV+vfvj1mzZuHOnTtwdXXF+vXrMXr0aJW2YWhoCAcHB8FcWs3nvQtVVdWDvr4+DA0NOa6MkLpZsWIFoqKi8MEHH8DW1la7G1f7Gso6iIiIYNbW1kxfX59ZWloyBwcHxhhjBgYGzMbGhrm7u8v/vHr1ijHG2OLFi9mGDRvkrxEbG8tcXFyYg4MDCwkJUXoZaF0uE+UDXbx0V4ioB+7xqX6ZTMbmz5/PALAlS5aoXFN9XsbeYO7GLjRMB89NFSLqgXt8qj89PR0eHh6Ijo7G3LlzVV6vYg91ySdeDI8QQgifFRQUgDEGe3t7pKWloVWrVpzVIuyjXIQQomFv3rzBoEGDMGXKFADgNLAB2tMmhBClsrKyMHDgQNy+fRtHjhzhuhwAFNqEEFKlp0+fol+/fnj+/DnOnj0rny6ZaxTahBBShd27d0MqleL8+fPyqTP4gMa0CSGknJycHADAggULcOPGDV4FNkChTQghchcuXEC7du1w+vRpiEQitG7dmuuSKqHQJoQQAMePH0ffvn3RoUMHvPvuu1yXoxSFNiGkwdu3bx+GDBmC3r174/jx47y+II9CmxDSoBUWFuLzzz/HqFGjcOjQIZiYmHBdUrXo7BFCSIPEGENubi6aNm2Kc+fOwdzcXBCzavK/QkIIqWeMMSxYsAC+vr7Iz89H69atBRHYAO1pE0IamOLiYsyYMQPbtm3DmjVrBHNDjjIU2oSQBqOgoABjx45FbGwsdu3ahfHjx3NdktootAkhDca5c+eQkJCAmJgYDB48mOtyaoVCmxCi83Jzc9G4cWP07dsXDx48gJWVFdcl1ZowRt4JIaSWnjx5gnfffRf/+c9/AEDQgQ3QnjYhRIfdu3cPffv2hUwmw6BBg7gup17QnjYhRCelpKTAz88PJiYmuHDhAlxcXLguqV5QaBNCdNJXX30FW1tbnDt3Tvt3TNcgGh4hhOiU3NxcNGnSBNu2bUNJSQmv5xGpDdrTJoTojL1798LR0REPHz5EkyZNdC6wAQptQoiO2LBhA8aNG4d+/frp1HBIRRTahBBBY4xh6dKliIiIwJw5c7B9+3YYGOjuyC+FNiFE0DIzM7FixQp88cUXWLVqlWAmfqot3f06IoTotOLiYpSUlMDW1hZ37txBmzZtuC5JK3T7K4kQopPy8/MxYsQITJo0CQAaTGADtKdNCBGYnJwchISE4OLFizh48CDX5WgdhTYhRDCysrIwYMAA/Pnnnzh58iR69uzJdUlaR6FNCBGM7du3Iz09HYmJiejSpQvX5XCCxrQJIbyXl5cHAPj444+RkpLSYAMboNAmhPBcSkoKOnfujJMnT0IkEjWog45VodAmhPDW+fPnERgYCCsrK3h6enJdDi9QaBNCeOnYsWPo168fvL29cezYMbRq1YrrkniBQpsQwjvFxcWYN28e+vfvj7i4ODRr1ozrkniDzh4hhPBKfn4+jI2NcebMGbRu3Rr6+vryA5FEy3vac+bMgb29PUQiEVJTU+XPL1u2DC4uLtDT00NcXJzS9dPT02FgYAAPDw/5n/v372ujdEKIhjHG8OWXX8LHxwdv376FWCzW6YmfakuroR0WFobk5GS0bdtW4fnevXvj2LFj8Pf3r/E1zMzMkJKSIv/j4OCgqXIJIVoik8nw0UcfYfHixRgxYgRMTEy4Lom3tPo1piyUu3Xrps0yCCE8UlxcjGnTpmHXrl1Yv349Zs6cyXVJvCa4/3tkZ2eja9euKCkpQUhICBYtWgR9fX2ly0ulUjDG5I+NjIxgZGSkjVLrpKzm8rULDfXAD3zvISkpCXv37sWePXswZsyYSnXyvX5VVOyhLr0IKrTFYjEyMzNhYWGB169fY9SoUVi5ciUWLFigdB07OzuFxwsXLsSiRYs0XWq9yc/P57qEOqMe+IFvPeTn58PIyAi+vr64desWbG1tqz3gyLf6a6Osh7ocWBVUaBsZGcHCwgIA0LJlS0yZMgU//vhjtaGdkZGhcJ84Ie1plx1FF4lEXJdTK9QDP/Cxh5cvXyI4OBijRo3CvHnz4OzsrHRZPtavroo9FBUV1fq1BBXaL168QIsWLWBoaIiCggIcOnSoxjkITE1NBX1zT5FIJNhf1DLUAz/wpYfHjx+jX79+eP36Nfr166dyTXypvy7KeqhLH1o9e2TWrFmwsbFBZmYm+vTpA0dHRwDAV199BRsbG1y6dAmTJk2CjY0NXr58CQD47LPPsHHjRgBAcnIyunTpAnd3d3h6esLKykpQQx2ENHR3796Fn58f8vLykJycDHd3d65LEhwRE/LofjWys7NhamoKqVQqyD1txhjy8vJgYmIi2L0L6oEf+NTD2LFj8dtvv+HkyZOwtrZWaR0+1V9bFXuoSz4JaniEECJMBQUFMDIywubNm1FQUABzc3OuSxIsmnuEEKJRcXFxcHJywv3799G0aVMK7Dqi0CaEaMwPP/yAkJAQeHt7qzwcQqpHoU0I0Yi1a9ciPDwcEyZMwP79+2FsbMx1STqBQpsQUu8kEgkWLlyIjz/+GFu3bqWJn+oRhTYhRE4ikWDiqEGQSCS1Wl8mk6GoqAhisRi///47vvnmG8Ge8cFXFNqEEAClgR0eGoAw8TGEhwaoHdxFRUWYNGkSxo8fD8aYfBpmUr8otAkh8sBeNzQNA91FWDc0Ta3gzsvLw/Dhw7F3716EhIRQWGsQhTYhBFGR0xHpkwZnq9KwdbYSIdInDVGR02tcNzs7GwMGDMDp06dx5MgRjB49WtPlNmgU2oQQLF+9BasvO+Lus9ILpO8+Y1h92RHLV2+pcd2tW7ciJSUFp06dwoABAzRdaoNHoU0IgVgsxp6YJMyKdUT8bwyzYh2xJyYJYrFY6TqFhYUAgLlz5yIlJQU9evTQVrkNGoU2IQTAP8F9UBJcY2D/+eefaN++PU6ePAk9PT3Y29trr9AGjkKbECInFouxc19ctYF9/fp1+Pn5wcTEBB07dtRidQSg0CaEqCExMRG9evWCg4MDzp07R5emc4BCmxCikpKSEsyaNQvdunXD6dOn0bJlS65LapDo2lJCSI0KCwvRqFEjnDhxAq1btxbELft0Fe1pE0KqtWbNGvj6+uLNmzewsbGhwOYYhTYhPFHXeT/qG2MMn3/+OSIjI9G7d280adKE65IIKLQJ4YW6zvtR32QyGebOnYvo6GgsX74cX3/9NV2azhM0pk0Ix8rP++FsJYKTZem8HzWdK61JycnJWL9+PTZt2oT333+fkxpI1WhPmxCO1WXej/pWVFQEAPD398d///tfCmweotAmhGN1mfejPkmlUvTp0wcrV64EADg5OWl1+0Q1FNqEaFnFA461mfejvr148QK9evXCrVu34Ovrq7XtEvVRaBOiRcoOOKoz70d91jJx1CBcvXoVPXv2hEQiQVJSErp3767xbZPaowORhGhJTQccy+b90GYtkT5pGDP4OPL0LZGcnAwHBwetbJ/UHu1pE6IlfDngWPEuNXFzZWjKXiAnJ0erdZDaodAmREv4csCxqi+P/zeyBCMH+XF+fjipGYU2IVrChwOOANBrQBjm/MAUvzxOMiwakMvJaYZEPRTahGgRFwccy9u1axemTZsGG1c/DFzFSr88djMsGQrs+t2p0l4/3y6tJxTahGidKjca0ISXL19i1qxZGDVqFAwLn+H/jQI+/olhQg9g+i5DrNy4X6Emvl1aT0pRaBOi4xhjKC4uRuvWrXHt2jXoF0vxoe99DPXUw9lPRDhzG1gYXIRVX30qX6fiwcp1Q9MouHmCQpsQHSaTyfDBBx9g2rRpYIzBxcUFK8odEBWbibBwkKjS0AhfznQhlVFoE6KjioqKMH78eGzcuBEBAQHyWfpUOSDKlzNdSGUU2oTooLdv3yI0NBQHDhzATz/9hMmTJyv8vKYDonw504VURqFNiA76/vvvkZiYiPj4eISFhVW5TE0HRLk+04VUTcQYY1wXoQnZ2dkwNTWFVCpF8+bNuS5HbYwx5OXlwcTERLCTz1MP2ldcXAwDAwPIZDLcu3cPLi4uguuhIqHXD1TuoS75pNU97Tlz5sDe3h4ikQipqany55ctWwYXFxfo6ekhLq76uRfi4uLg6uoKR0dHDB8+HG/evNF02YQIwqNHj+Dh4YGTJ09CT08PLi4uXJdENKBWoX337l388ssvuHTpklrzFYSFhSE5ORlt27ZVeL537944duwY/P39q13/zZs3mDp1Kg4fPoy0tDSIxWIsXbq0Ni0QolNu376NHj16IC8vD46OjlyXQzRI5dDOyclBdHQ0bG1tERwcjMWLF2P27Nmws7PDgAED8Msvv9T4Gv7+/rCxsan0fLdu3VSaXSwhIQHe3t5wdXUFAERERGDv3r2qtkCITrp69Sr8/f3RsmVLJCcno127dnV+TboSkr9Unpq1V69eCA8Px9WrV2FlZSV/XiaT4fz589i4cSPS0tI0enuijIwMhb10e3t7PHnyBDKZDHp6VX//SKVSlB+2NzIygpGRkcZqrC9lNQv5kAP1oHkymQxTp06Fs7Mz4uLi0KJFi0q1qtuDRCLB+GEBiPS5j/BQf+w+xO1BSL5/Bqqo2ENdelE5tC9cuFBl2Onp6SEgIAABAQEoKCiodSGqUvdAhJ2dncLjhQsXYtGiRfVZkkbl5+dzXUKdUQ+aUVJSAn19fezfvx/m5uYwNjZGXl6e0uVV6UEikWDamPewIfTB/+b8vo9xIf74fu8Jzs8e4eNnoK6yHqr7nGqicmhXFdjp6enIzc1Fx44dlS5Tn+zs7BSGYdLT02Ftba10Lxso3Tsvf3RWSHva+fn5MDY2FvQRc+qhZhKJBFEfvo/lqzarFIwSiQTjRgzA32/1kJSUVOMBx+p6qLjt6IVz8VH3BwpXQn7U/QGiF87Fzp+O1r7JOtDF36OyGyjXRq3vXLN27VrExMRAJBLB2dkZGzZsqHURqurfvz9mzZqFO3fuwNXVFevXr8fo0aOrXcfU1FSQp/yVEYlEgv1FLUM9KFc6FBGISJ80jB8WWOP50BKJBIN6uSN6QBYWxjZDdna2yr/fFXuoatsrVm9BeGgAnCxLL2EvuxJyT8wWzj9DXfo9qlMfTEUXL15UeDxy5Ej5393c3FR6jYiICGZtbc309fWZpaUlc3BwYIwxtmzZMmZtbc0aNWrEzM3NmbW1NXvx4gVjjLHFixezDRs2yF8jNjaWubi4MAcHBxYSEsKkUmmV25JKpQyA0p/znUwmY7m5uUwmk3FdSq1RD9V7+vQpC+rmxP67TMRKtumx/y4TsaBuTuzp06dVLv/kyRPm7thS5eWr66G6bZf97Mhc1V5f03Tx96gu+aTyxTWTJ09G06ZNsXz5cjRp0gSzZ8+GhYUFRCIRzpw5g8TExNp/c2gAXVzDPepBOYlEgv4BHvgo4AXGd/9neC/+N4aDkuAq7xU5uF8PvN/hEga6i1RavroeJo4ahDDxMaWvJZFIEBU5HctXb+F8LFsXf4+0cnHN9u3bMXjwYAwYMABHjx7FqlWrYGNjAzMzM8TExKjdBCENVdm0p18Gv8SyY4ZIvisDoHxSppKSEgDA5p0HsSLJVqVJnGo6Za+mCaG4mvObqEDdXfO3b9+yefPmsTFjxrDnz5+rvWuvLTQ8wj3qobKqhiWc2xiyHdNQ5VBEbm4uCw4OZitXrlRYv7qhi4rL3Lhxg40L68+ePHlS7XJcD4Moo4u/R3XJJ5VDOyMjg82fP58tWrSIvXjxgl27do0FBQWxrVu3qr1RbaDQ5h71UNmEkQPZkbmlgV3258hcEevsZFEpNP/66y/m5+fHmjRpwk6dOiV//unTp2zCyIHVBnbVXwqOldap7rX4Qhd/j+qSTyoPj4waNQpubm6wsbHBxIkT4eXlhZMnT+LVq1cYOHCgpv4jQIha+H4ln7JhieNJKQpDEc+fP0dgYCD++OMPnDlzBn369JH/rLqhiyrvtD6iGGduA+uG3q909xkaBhEgVdO9Y8eOjDHG8vPzmYeHh8LP7t+/r/a3habRnjb3tN2DJv67r4keVKlz8uTJrE2bNiw1NbVWr61wVkh7sMz/iOR79RNGDqyvVrRCF/8taG1P29XVFW5ubpgxY4bCz+pjrgNC6kJI9zSsbp5qmaz0oOTq1atx4cIF+YVr6r522c0Lhn5ngCVDAbGZiO4+oyvUSfjs7GyWm5ur9jcDF2hPm3va7EHZWHFd9yrrowdVx42vXLnCOnTowNLS0mq9rYrbvHnzJgvq5vi/vfrKY9pCoIv/FrSypw0AzZo1Q+PGjTXz7UFIHfD1noZl/wMIEx+rds//zJkzCAoKQosWLWBubl7n7ZaNVXt4eGD3oSTse/Ie5xM/kfpRLzdB8PT0rI+XIaTW+HhPQ1WHbGJiYhAcHIyePXvi5MmTMDMzq9c6xGIxNu/8mQJbR9RLaMfHx9fHyxBSJ3y7p2FVZ3JE+qQhKnK6fJmsrCxMnDgRoaGhiI2NrfZ/snw/M4Zoh9qhvW3btkrPJSQk1EsxhNQVn05hq2nIRiaTwdzcHBcuXMAPP/yARo0aKX0tVYdZiO5TO7S/++47lZ4jpKFTNmRjZWWFRYsWYeLEiWCMwc3NDfr6+kpfR0hnxhDNU3lq1mvXruHKlSt49eoV1q9fL39eKpWisLBQI8URInRlwR0VOR17YrbAwsICERER2LhxI7755huVJkCqbpiluomiiG5SObSfPHmCa9euITc3F1evXpU/37x5c+zYsUMTtRGiE8qGbAoLCzFu3DgcOHAA27Ztw+TJk1Vaf3k1c1yThkfl0B46dCiGDh2KhIQEDBgwQJM1EaKTtm7dipiYGBw8eBChoaEqr1e2tx4eGoBIn7T/BTb3B1oJN9Qe0165ciVWrFiB69eva6IeQnRO2VWOM2bMwK+//qpWYJfh25kxhDtqh/bChQshlUoxc+ZMWFhYYMSIEdi8ebMmaiNE8J49ewYfHx+cOnUKenp6cHd3r/Vr8enMGMIdtUM7KCgIy5Ytw+nTp/H111/j2rVrmDNnjiZqI0TQHj58CD8/Pzx58gRt2rThuhyiI9QO7cWLF6NHjx4IDAzE77//jnXr1iErK0sTtREiWKmpqejRowdEIlGtJn4iRBm178a+ceNGODs7Y/r06ejbty8cHR01URchgsUYw4QJE2BhYYETJ07A0tKS65KIDlE7tF++fImUlBScPn0ac+bMwaNHj9C9e3ds2UKnHxEik8mgp6eHAwcOwNzcvN7nESGkVnOPtGrVCubm5mjRogWysrIUztsmpKH6+eef4e/vj5ycHDg4OFBgE41QO7RdXV3h5+eH5ORkDBw4ELdu3UJKSooGSiNEOLZu3YqRI0fC1tYWRkZGXJdDdJjawyNHjhyBs7OzJmohRJC++eYbLFiwABMmTADLz0JWVhadlkc0RuU97cePHwNAtYH99OnTuldEiIBcvXoVCxYswNy5c5H550WMoFn4iIapHNpjxozBtGnTcPbsWeTn58ufT09Px/r169GtWzdcvnxZI0USokm1maeasdLpVrt27YqEhAT8fvkY1g29T7PwEY1TObSTk5Px3nvvYenSpWjZsiVatGgBExMTBAUF4eHDhzhw4ACGDRumyVoJqXe1mae6sLAQo0aNwurVqwEAe7d/V+PNDgipL2qNaY8YMQIjRoxAcXExXr16hcaNG6N58+aaqo0QjSo/T7WzlQhOlqV7yNXN7ZGbm4thw4YhMTERY8eOBUCz8BHtqtUpfwYGBrCysqLAJoKmyu3Aynv9+jX69u2LixcvIiEhASEhIQD4eX9Korvq5R6RhAiRundwj4qKwt27d/HLL78gKChI4Wc0Cx/RFhErO6KiY7Kzs2FqagqpVCrI/xEwxpCXlwcTExOV7m7CR0LooWyIRNk81YwxvH37Fo0bN0Z2djaePXsGFxcXDitWnxA+h+oIvX6gcg91ySfa0yYNWk17yKmpqfDz88P9+/dhamoquMAmuketA5HPnz/H+vXrcePGDQCAp6cnZs6cCSsrK40UR4g2lM1TXdGlS5cQHBwMOzs7NG3alIPKCKlM5T3tO3fuoHPnzrhz5w569+6NoKAg3LlzB+7u7rhz544mayRE606ePIk+ffrAzc0Nx48fp5n6CG+ovKc9f/58rFmzBqNHj1Z4fu/evfj4448RHx9f78URwoW///4bI0eORGBgIPbv3y/YcVSim9Ta064Y2EDplZJ3796t16II4QpjDGZmZjhz5gwOHz6Mxo0bc10SIQpUDu3qTjLR0RNQSAOzYsUKTJ48GTKZDF5eXjA0NOS6JEIqUTm0XV1dsW/fvkrP//TTTyrP+jdnzhzY29tDJBIhNTVV/vyLFy/Qv39/ODk5oVOnTkhOTq5y/fT0dBgYGMDDw0P+5/79+6q2QEiVGGP45JNPEBUVBTs7OxoOIbym8pj2N998g8DAQMTExMDX11d+77vExEQkJiaq9BphYWFYsGAB/Pz8FJ6PioqCj48Pjh8/jqtXryIsLAz379+HgUHl8szMzGj+blJvSkpKMHPmTGzZsgWrVq1CZGQk1yURUi2VQ7t9+/ZISUnB+vXrcfLkSQClp/ytWrVK5TtN+/v7V/n8/v378fDhQwCls6ZZWloiOTkZgYGBqpanlFQqVRi+MTIyEsQk9WU1C3noSQg9bNu2Ddu2bcP27dsxceLESrUKoYeaCL0HodcPVO6hLr2odZ62WCzGv//9bwBAcXFxlXvC6srKyoJMJkPr1q3lz9nb2yMjI6PK5bOzs9G1a1eUlJQgJCQEixYtgr6+vtLXt7OzU3i8cOFCLFq0qM51a0v5aXCFio89MMYgEokwevRoODs7o2vXrsjLy1O6PB97UJfQexB6/cA/PVT3u1YTtVP39u3bGDt2LLKysvD48WNcv34d+/fvx4oVK2pdRMUxRGXfQmKxGJmZmbCwsMDr168xatQorFy5EgsWLFD62hkZGQqXiQppTzs/Px/GxsaCHWPlaw+vX7/G8OHD8emnn6J3795K/wcI8LcHdQi9B6HXD1TuoaioqNavpfZl7LNmzcJ3332HVq1aASgdIqnLOdrm5uYASu/yXubRo0eV9pCB0sC1sLAAALRs2RJTpkzB+fPnq319U1NThT9lb5oQ/gDgvAZd60EikSAwMBCpqakwMzOrcflnz57h/YnD8ezZM85r16XPoaHVX1UPtaV2aOfk5CgcSBSJRHU+NWrEiBFYt24dgNLbNz179qzSwUqg9CyTsm+ogoICHDp0CF26dKnTtknDcf/+ffj5+eHvv//G+fPn4e3tXe3yEokE44cFYJT1CYwfRneiIfygdmgbGBigqKhI/k2RmZkJPT3VXmbWrFmwsbFBZmYm+vTpA0dHRwCl58devHgRTk5OmDRpEnbv3i0fL//ss8+wceNGAKV3z+nSpQvc3d3h6ekJKysrQY1PE+4wxjB27FgYGhriwoULaN++fbXL/3ODhLJbiN2nW4gRXlB7atY9e/bgp59+wq1btzBlyhTs2rULy5Ytq/JqSS7R1Kzc40sPZQcd7969CzMzM/kQW3UmjhqEMPExDHT/p+743xgOSoKrnFyKz/jyOdSW0OsH6ndqVrUPRIaHh6Ndu3aIjY3F27dvsXPnTvTs2VPdlyFEK06cOIEVK1YgNjZW5YvAALqFGOEvlYdHIiIi5H9/+fIlVqxYga+//poCm/DW/v37MXjwYDRp0kTh9FRV7r7+zy3EHP53CzEHuiMN4QWVQ/vy5cvyv0dHR2ukGELqy+bNmzF69GiMGjUKhw4dgomJCQD17r4uFoux+1AS9j15D7sPUWATfqjVhFFCvjKJ6L6UlBTMmDEDs2fPxs6dO+VnN5W/+3rpwcU0lYJ7886fKbAJb6g8pl1QUID//ve/YIwp/L1Mhw4dNFIgIaoqO+Do4eGBc+fOwc/PT+HAVXV3XxfawUXScKkc2m/fvkVwcLD8cfm/i0QiPHjwoH4rI0QNxcXFmDFjBtzd3TFnzpwqj7XQwUWiC1QO7fT0dA2WQUjtFRQUYOzYsYiNjcWOHTuULld2cLG6u68Twnd0N3YiaG/evMGgQYMQHx+PmJgYhIeHV7t8TXdfJ4Tv6j5NHyEcioqKwpUrV3DixAkEBASotI6yu68TIgS0p00Eqewg+Jdffonz58+rHNiECB2FNhGctLQ09OzZEw8fPoSZmRnc3d25LokQraHQJoJy69Yt+Pn54dWrV9Xe/IIQXUWhTQTjwoUL8Pf3h7W1Nc6fP1/lnOuE6DoKbSIIOTk5GDJkCDw8PHD27FmF29MR0pDQ2SOE9xhjaNasGeLj4+Hu7i6fR4SQhoj2tAmvbdy4EdOnT4dMJoOPjw8FNmnwKLQJLzHG8NVXX2HmzJlo0qQJ1+UQwhsU2oR3GGOYP38+Fi5ciOjoaKxevVrlW9oRoutoTJvwzu7du7Fy5Up8++23+OCDD7guhxBeodAmvDNu3DjY2dkhMDCQ61II4R36PyfhhZycHAwaNAhnz56Fvr4+BTYhSlBoE85lZWWhd+/eOH/+PF3lSEgNaHiEcCozMxP9+vXDq1evkJiYiC5dunBdEiG8RqFNOMMYw5gxY5Cbm4vk5GQ4OztzXRIhvEehTTgjEomwZcsWNGnSBLa2tlyXQ4gg0Jg20brz589jwIAByM3Nhaurq1YCWyKRYOKoQdXeeZ0QIaDQJloVHx+Pfv36oaCgACUlJVrZpkQiQXhoAMLExxAeGkDBTQSNQptozY8//oiQkBC89957OHbsGJo3b67xbZYF9rqhaRjoLsK6oWkU3ETQKLSJVvzxxx8IDw/HuHHjcPDgQRgbG2tlu1GR0xHpkwZnKxEAwNlKhEifNERFTtfK9gmpb3QgkmhFx44dcfz4cfTp00er84gsX70F4aEBcLIsDe67zxhWX3bEnpgtWquBkPpEe9pEY2QyGT755BOsW7cOANCvXz+tT/wkFouxJyYJs2IdEf8bw6xYR+yJSYJYLNZqHYTUFwptohHFxcWYOnUq1q1bB5FIxGktZcF9UBJMgU0Ej4ZHSL3Lz8/H6NGjER8fj61bt2LSpElclwSxWIyd++K4LoOQOqPQJvVu0aJFOHHiBGJiYhAUFMR1OYToFBoeIfXu008/xdmzZzFw4ECuSyFE51Bokzopu9Lw+vXr6NevHx49eoQWLVrAx8eH69II0UlaDe05c+bA3t4eIpEIqamp8udfvHiB/v37w8nJCZ06dUJycrLS14iLi4OrqyscHR0xfPhwvHnzRhulkyqUv9JwRLAPbt++jcLCQq7LIkSnaTW0w8LCkJycjLZt2yo8HxUVBR8fH9y7dw/bt2/HuHHjUFxcXGn9N2/eYOrUqTh8+DDS0tIgFouxdOlSbZVPyql4peGxyBI4WBmiadOmXJdGiE7T6oFIf3//Kp/fv38/Hj58CADo2rUrLC0tkZycXOnuJQkJCfD29oarqysAICIiAsHBwfjqq6+UblMqlYIxJn9sZGQEIyOjOnaieWU1l6+dTz6p4krDeX6P8EnkdOz86SgA/vegCuqBe0KvH6jcQ1164fzskaysLMhkMrRu3Vr+nL29PTIyMiotm5GRobCXbm9vjydPnkAmkym9aMPOzk7h8cKFC7Fo0aJ6ql7z8vPzuS6hSkuWrcG0MXfgZPlAfqXhfy62w/d71yAvL09hWb72oA7qgXtCrx/4p4eK/0bUwXloA6h08UV130LqXqiRkZGhMDGRkPa08/PzYWxszPnFKRX98MMPSE5Oxp6YJEwYHohIn/tYfdkBPxxWvHCFzz2oinrgntDrByr3UFRUVOvX4jy0zc3NAQAvX76U720/evSo0h4yULrX/Msvv8gfp6enw9rautpLo01NTbUym5ymiEQiXv2ifvfdd/jggw8wadIkWFlZYU/MOURFTseemC1KrzTkWw+1QT1wT+j1A//0UJc+eHHK34gRI+TzU1y9ehXPnj2Dn59fpeX69++Pq1ev4s6dOwCA9evXY/To0VqttaFijOGLL77ABx98gI8++ghbt26Fvr6+/EpDujScEC1hWhQREcGsra2Zvr4+s7S0ZA4ODowxxp49e8b69u3LHB0dWYcOHVhiYqJ8ncWLF7MNGzbIH8fGxjIXFxfm4ODAQkJCmFQqrXJbUqmUAVD6c76TyWQsNzeXyWQyrkthjDG2d+9eBoAtXbpU5Zr41kNtUA/cE3r9jFXuoS75JGJMwIdkq5GdnQ1TU1NIpVJBDo8wxpCXlwcTExNe/JewuLgYJ0+eRHBwsMrr8K2H2qAeuCf0+oHKPdQln3gxPEL4qWzip3PnzsHAwECtwCaEaAaFNqlSdnY2BgwYgCNHjiA3N5frcggh/8P52SOEf16+fIkBAwYgLS0NJ0+erPKgMCGEGxTapJLRo0cjMzMTSUlJcHd357ocQkg5FNqkkm+//RaNGjWCk5MT16UQQiqgMW0CALhx4wZCQ0Px9u1bdOzYkQKbEJ6i0CZISkpCYGAgnj59qhPzOxCiyyi0G7ijR4+if//+ePfdd3H69Gm0bNmS65IIIdWg0G7A7t69i9DQUAQHByM+Ph7NmjXjuiRCSA0otHVI2a2/JBKJSss7Ozvj0KFD2LdvnyBmPiSEUGjrjPK3/goPDVAa3IwxREdHY+PGjQCAIUOGwMCATiIiRCgotHVAxVt/rRuaVmVwy2QyREZG4vPPP8fff//NTbGEkDqh0NYBUVXc+ivSJw1RkdPlyxQVFWHixIlYu3YtNm7ciKioKK7KJYTUAYW2Dli+egtWX3bE3WelEzbefcaw+rIjlq/eIl9myZIl2LdvH/bu3YsZM2ZwVSohpI4otHWAWCzGnpgkzIp1RPxvDLNiHbEnRvHWX/PmzcOJEycwatQoDislhNQVhbaOKAvug5JgeWC/fPkSISEhyMjIQMuWLdGrVy+uy6ySume9ENKQUWjrkPK3/srIyEDPnj1x+fJlSKVSrktTStWzXgghpSi0ddCdO3fQo0cPFBYW4sKFC3Bzc+O6pCqpetYLIeQfFNo6Jj8/H3379oWpqSmSk5Ph4ODAdUlKqXLWCyFEEYW2jjE2NsaOHTtw7tw5tGnThutyqqXKWS+EEEUU2jriyJEjmDt3Lhhj6N27tyAmflLlrBdCiCIKbR2wa9cuDBs2DE+fPkVxcTHX5ailqrNeCCHKUWgL3Jo1azBx4kRMnjwZP/30EwwNDbkuSW3lz3ohhFSPQlvADh8+jMjISCxYsACbN2+Gvr4+1yURQjSMpncTsEGDBuHAgQMICwvjuhRCiJbQnrbAFBUVYfr06UhOToaBgQEFNiENDIW2gOTl5SE0NBQ7d+7E8+fPuS6HEMIBGh4RCKlUiiFDhuDatWs4evQo3nvvPa5LIoRwgEJbIMaOHYtbt27h9OnT8PX15bocQghHKLQFYvny5QDA23lECCHaQWPaPHbnzh2MHTsWb9++hZubGwU2IYT2tPnq2rVrGDBgAMRiMXJyctC4cWOuSyKE8ADtafPQ2bNnERQUBAcHByQmJsLS0pLrkgghPEGhzTMPHjzAgAED0L17d8TFxQli4idCiPbQ8AjPtGvXDrt27cKQIUMgk8m4LocQwjO0p80Tq1atwubNmwEAI0eOhJGREccVEUL4iDehffz4cXh7e6Nz587w8fHBb7/9VmmZxMRENG7cGB4eHvI/eXl5HFRbfxhj+PTTT/HRRx8hPT2d63IIITzHi+GRv/76C+Hh4Th//jzat2+PpKQkjBs3DqmpqZWW7dChA65du8ZBlfVPJpNh9uzZ2LBhA1asWIEFCxZwXRIhhOd4sad9//59WFhYoH379gCAgIAAPHr0CDdu3OC4Ms3697//jU2bNmHLli0U2IQQlfBiT9vJyQkvX77E5cuX4ePjg5iYGLx58wbp6enw9PRUWPbPP/+Ep6cn9PX1MXnyZERERFT72lKpFIwx+WMjIyPejBdHRETAy8sLAwcOVKgRgPxxxeeFhHrgB6H3IPT6gco91KUXXoS2qakpfv75Z0RFRSEnJwd+fn7o0KFDpbuweHp6IjMzE6ampsjMzERwcDBatWqFkSNHKn1tOzs7hccLFy7EokWLNNKHKv7++2/MnTsXS5cuhY2NDYKCgqodl8/Pz9didZpBPfCD0HsQev3APz3U5ViciPHw66ugoABWVla4evUqHB0dlS731Vdf4enTp1i7dm2ln2VnZ8PU1BQZGRlo3ry5/Hku97SfP3+O/v37IyMjAydOnIC3t7fSZRljyM/Ph7GxMUQikRarrD/UAz8IvQeh1w9U7iE7OxtmZmaQSqUK+aQKXuxpA4BEIpHfI/Df//43goKCKgW2RCKBpaUl9PT0kJOTg7i4OEydOrXa1zU1NVX7TdGE9PR09O3bF7m5uTh37hw6deqk0noikUiwv6hlqAd+EHoPQq8f+KeHuvTBiwORALB48WK4urrC0dERjx49wtatWwEA06ZNw5EjRwAAP//8M9zc3ODu7g4fHx/07dsXkydP5rJslRQUFCAoKAgymQwXLlxQObAJIaQiXg6P1Iey4ZHa/PdDExISEuDh4aHyHccZY8jLy4OJiYlg9y6oB34Qeg9Crx+o3ENd8ok3e9q66MyZM5g3bx4YY/IZ+wghpC4otDUkJiYGwcHBSE1NRWFhIdflEEJ0BIW2Bmzbtg1hYWEIDQ3FkSNHeHNeOCFE+Ci061lCQgKmTp2K6dOn44cffkCjRo24LokQokN4c8qfrujTpw927tyJ8ePHC/agCSGEv2hPux6UlJTgww8/xMWLF2FoaIgJEyZQYBNCNIJCu44KCwsxbtw4fPvtt7h37x7X5RBCdBwNj9TB27dvERYWhjNnzmD//v0YPnw41yURQnQchXYdTJgwAefOnUN8fDz69OnDdTmEkAaAhkfqYPHixTh9+jQFNiFEayi01fTw4UNMnjwZeXl58jlQCCFEW2h4RA1//PEH+vXrBxMTE2RlZcHGxobrkgghDQztaavoypUr8Pf3R6tWrZCcnEyBTQjhBIV2BRKJBBNHDYJEIpE/l5GRgd69e8tvOmxlZcVhhYSQhoxCuxyJRILw0ACEiY8hPDRAHtx2dnZYv349Tp48CTMzM26LJIQ0aBTa/1MW2OuGpmGguwjrhqZhcJAHVq1aBaD09L7GjRtzXCUhpKGj0P6fqMjpiPRJg7NV6eXnzlYiLHnvJfZs/Zbjyggh5B8U2v+zfPUWrL7siLvPSm/kc/cZw+L4Fjh68gLHlRFCyD8otP9HLBZjT0wSxu5oifjfGMbuMEfCuVS0adOG69IIIUSOQrscsViMXT8n4tsbXRCfeItuD0YI4R26uKaCTp064dS561yXQQghVaI9bUIIERAKbUIIERAKbUIIERAKbUIIERAKbUIIERAKbUIIERAKbUIIERAKbZ4qKCjA0qVLUVBQwHUptUY98IPQexB6/UD99iBijLF6qIl3srOzYWpqCqlUiubNm3NdjtqkUinMzMzw999/w9TUlOtyaoV64Aeh9yD0+oHKPdQln2hPmxBCBIRCmxBCBERn5x4pG/XJzs7muJLaKas7OzsbIpGI42pqh3rgB6H3IPT6gco9lD2uzei0zo5pZ2ZmwtbWlusyCCFEqcePH6t9k3CdDW2ZTIanT5+iWbNmgv12JoToJsYYcnJy0KZNG+jpqTdKrbOhTQghuogORBJCiIBQaBNCiIBQaPPIwoUL0b59e7i7u+Pdd9/FL7/8onTZK1euwMPDA87OzujduzckEokWK1Vu27ZtcHNzg4GBAb777rtqlxWJROjcuTM8PDzg4eGB8+fPa6nK6qnTAx8/h7dv32LMmDFwdHSEs7MzDh06pHRZPn0G9+7dQ/fu3eHs7Ix3330Xt2/frnK5rVu3wsnJCQ4ODnj//fdRXFys5UqVU6WHxMRENG7cWP6ee3h4IC8vT/WNMMIbx44dY2/fvmWMMZaSksLMzMxYXl5epeVkMhlzcHBgZ8+eZYwx9s0337DRo0drs1SlUlJS2O3bt9n48ePZ2rVrq10WAMvJydFSZapTtQe+fg7R0dFs4sSJjDHGHjx4wCwtLdnr16+rXJZPn0GvXr3Y9u3bGWOMHThwgPn4+FRa5sGDB0wsFrNnz54xmUzGBg8ezDZu3KjlSpVTpYezZ88yLy+vWm+DQpunSkpKWLNmzdjjx48r/ezXX39lHTp0kD/Ozs5mxsbGrLCwUJslVmvixImCDe0yNfXA18+hQ4cO7Ndff5U/HjFihDxIKuLLZ/D8+XNmamrKioqKGGOlX4iWlpbs4cOHCst9/fXXLCIiQv44Pj6eBQQEaLFS5VTtoa6hTcMjPLV9+3Y4ODhUeQ5nRkYG2rZtK3/crFkzNGvWjBf/NVdXYGAg3N3d8dFHHyE3N5frctTC18+hYl329vbIyMhQujwfPoPHjx+jTZs2MDAovd5PJBLBzs6uUt3q9qZNqvYAAH/++Sc8PT3RtWtXrF+/Xq3t6OwVkXzUs2dP/Pe//63yZzdv3pRfDHTmzBlER0fj1KlTSl+r4rnnTEtnbqragyoePXoEOzs75Obm4l//+hfmz5+v9i9wbdRnD1x8DjXVX7Gu6mri6jOoiqrvpaq9cUGVHjw9PZGZmQlTU1NkZmYiODgYrVq1wsiRI1XaBoW2FqlykCcpKQmTJ0/G0aNH4eLiUuUydnZ2SE9Plz/OyclBTk4OxGJxfZWqVH0eqLKzswMANGnSBBEREXj//ffr7bWrU189cPU51FR/WV2tW7cGUBrMwcHBSpcFtP8ZVGRra4vMzEwUFxfDwMAAjDE8fvxYXl+Ziu952ZcOH6jaQ/lZ/WxsbDBmzBicP39e5dCm4REeOXfuHMaPH4/Y2Fi4u7srXc7Lywv5+flITEwEAGzatAkhISEwNDTUUqV199dff+Ht27cASq9e3bdvH7p06cJxVerh6+cwYsQIrFu3DgDw8OFDJCUlYciQIZWW49NnYGFhgS5dumDPnj0AgJ9//hn29vawt7dXWG748OGIiYnB8+fPwRjDxo0bMXr0aA4qrkzVHiQSCWQyGYDSL/q4uDj13vdaj4aTeufo6MgsLCyYu7u7/M+tW7cYY4xt2LCBLV68WL7sxYsXWefOnZmTkxMLDAxkmZmZXJWtYPfu3cza2po1btyYmZmZMWtra3bjxg3GmGIPFy9eZG5ubqxz586sQ4cOLDw8nGVlZXFZupyqPTDGz8/hzZs3bOTIkczBwYE5OTmxAwcOyH/G58/gzp07zMfHhzk5OTEvLy+WmprKGGNs6tSpLDY2Vr7c5s2bmYODA3vnnXfY1KlTOT/wW54qPaxdu5Z16NBB/r4vWbKEyWQylbdBl7ETQoiA0PAIIYQICIU2IYQICIU2IYQICIU2IYQICIU2IYQICIU2IYQICIW2AMyZMwf29vYQiURITU2tcfn8/Hx06NAB3t7eCs9nZGRg8ODBcHFxgaurK9auXYunT5/Kp4d0dHRUmDLyww8/xLRp0xSuwIuLi0PXrl3h4uKCdu3a4f3334dUKpX/fMqUKXBxcYGHhwf8/f2RkpJSb+9DVVSdprO6aW/XrVsHNzc3eHh4wM3NDd9++22l9T/77DPo6+vj0aNHCs/b29vLP5Oaptat6b2bPHmyfJrUrl274syZM7V+X1Sh6ntX3VS1c+bMUZhi1NjYWOH9u379Ovr374927dqhU6dO8PX1xeHDh2usrbCwEJ988gkcHR3Rvn17dOrUCdu3b69TvzpDA+eXk3qWlJTEHj9+zNq2bct+//33Gpf/6KOP2JQpUxRmEpPJZMzT05Pt379f/lgikSisV9PsYwkJCaxNmzbyC02KiorY7NmzWY8ePeQXB8TGxspnOTt69ChzcnJSqce2bduqtFx56kzTWd20t3///bd8OalUymxtbdlvv/0mf66kpITZ2tqywMBA9vnnn1equ+wzqW4bqrx3f/31l/x1b968yczNzVW66ELT752qU9VKJBJmbGws/71KTU1l5ubm7MiRI/JlMjMz2Y4dO2qsb8yYMWzYsGHszZs3jDHGHj58yFxdXdmmTZvUaVMnUWgLiCqhfe7cOTZ48OBKAXzq1CnWo0ePatetKrQDAgLY0aNHGWOM9ejRg61bt07h5wUFBcza2pqdPn260uu9fPmSNWrUiJWUlFS7XcZqFzy1naazumlvJRIJE4vF8itRGSsNXG9vb3b9+nXWtm1bhX6UfSYVt6Hue3f27FnWqlUrjYV2bd67mqaqXbFiBRs6dKj8cXh4OJs/f77S5VeuXMm8vb2Zh4cH69q1K7t8+TJjjLF79+4xExMT9urVK4Xl4+PjmY2NTbU1NgQ0PKJDcnNzERkZiQ0bNlT62e3bt9G6dWuMHj0aXbp0QWhoKB48eKDW69+4cQO+vr4KzzVq1AheXl64ceNGpeXXrFmD4OBgte82raraTtNZ1bS3Bw8eRMeOHdG2bVvMnz8fbm5u8p9t3boVU6ZMgaenJ1q0aKHSsEXFbaj63kVFRcHBwQHDhg3DgQMHKs0aV180McXptm3bMHXqVPnj69evV+q5vPHjx+Pq1au4efMmvv32W/m6N27cgJOTE8zNzRWW9/X1RWZmJl6+fFmnOoWOZvnTIfPnz8esWbNgbW2Ne/fuKfysqKgIp0+fxuXLl9GxY0ds3rwZo0ePxq+//qrWNqoKEVbFTAh79uzB/v37lc5IV1JSAi8vL/njsrF1ALC0tMSJEyfUrqeqOipSNu1tWFgYwsLCkJ6ejtDQUAQHB8PFxQWvXr3CqVOnsGXLFgDA1KlTsXXrVvTt21ftbajy3i1fvhzLly/H6dOnMX/+fFy4cAGNGjVSWIar9646Fy5cQHZ2ttLZBKty8+ZNLF26FFlZWTAwMMDt27dRWFhYqTaiiEJbhyQnJ+PYsWP44osvkJ+fj7/++gsdO3bEH3/8gbZt26JLly7o2LEjACA8PBwzZ85ESUkJ9PX1VXp9T09PXLx4UR4QQOkBoxs3bmDu3Lny5/bt24fo6GicOXMGFhYWVb6Wvr6+wkFKe3v7Gg9ahoWFIS0tDUBpMKo7Tacq097a29ujW7duiIuLg4uLC3bv3o3i4mJ5zyUlJcjKykJWVlalPcHqtqHqe1emT58+mD17Nn7//XeFgAa4ee9qsnXrVkycOFHhd8nLywuXLl1CaGhopeULCwsxfPhwJCYmwsvLC9nZ2TA1NUVhYSG6dOmCu3fvVnqPL126BBsbG/mUsw0Wt6MzRB2qHohkrPL49Js3b1i7du3ks9D9/PPPrHPnztWuw5jimHZ8fDwTi8Xs5s2bjLF/Dqb5+vrKx1737dvHHB0dWXp6utq9qev+/fuVDqZt2LChymWTkpKYra2t/EBgebdv35b//cWLF8zR0ZGdPHmSMcZYp06dWEJCgsLyw4cPZ2vWrJHXXfaZVLeNmt67oqIidvfuXfnyV65cYS1atFB6b8fyNP3elVE2pp2Tk8OaNWumUD9jjN26dYuZm5uzuLg4+XOPHz9mmzZtYlKplBkZGbHnz58zxhj78ssvFW59NnLkSDZs2DCWm5vLGPvnQGRNNTYEFNoCEBERwaytrZm+vj6ztLRkDg4O8p9VnLayTFUBfPz4cebu7s46d+7M/P395dNGVrdO+dBmjLHDhw8zT09P5uzszOzt7dnUqVMVznowMDBgNjY2CtPLVjygVJXaBA9jyqfpfPLkCXN3d5cvV920tzNnzmQdOnSQvzdlBwwvX77MLCws5GfDlH8Pyr7wyod2ddsoW0/Ze5efn8+6d+/OOnbsyDp37sx8fX3ZmTNnVHoPNP3eVTdVLWOMff/998zf37/Kbfz666+sb9++7J133mGdOnVi3bt3l/++rlixgrVt25b17NmTffPNNwqhnZ+fz+bNm8fatWvHXF1dWceOHdmWLVtq1aeuoalZCSFEQOjsEUIIERAKbUIIERAKbUIIERAKbUIIERAKbUIIERAKbUIIERAKbUIIERAKbUIIERAKbUIIERAKbUIIERAKbUIIEZD/D3Ll3SZe0khDAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "FeOi_predict.select_predictors(idx=3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can check if the right predictors are used with the *predictors* attribute" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array(['TiO2', 'Al2O3', 'CaO'], dtype=object)" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "FeOi_predict.predictors" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Coefficients of the linear regression were automatically calculated by the *select_predictors* method and we can access these with the *coefficients*, *intercept_error* and *errors* attributes." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fitted coefficients:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "intercept 11.585125\n", "TiO2 1.463365\n", "Al2O3 -0.230722\n", "CaO -0.169262\n", "dtype: float64" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "FeOi_predict.coefficients" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and their errors as standard deviations:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "intercept 1.203255\n", "TiO2 0.164396\n", "Al2O3 0.048417\n", "CaO 0.049363\n", "dtype: float64" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "FeOi_predict.errors" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The *model* attribute stores the selected model as a callable function, which we can use to predict melt FeO contents. The function accepts a pandas DataFrame or MagmaFrame with melt major element compositions in oxide wt. % as input. Let's try that out on the whole-rock data:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "name\n", "PI032 11.197161\n", "PI041 11.220927\n", "PI053 10.963646\n", "PI054 11.159202\n", "PI055 11.091599\n", "PI056 10.944143\n", "PI057 10.590489\n", "PI058 10.757390\n", "PI060 11.154730\n", "PI063 11.988362\n", "PI065 11.165823\n", "PI066 11.055243\n", "PI067 11.146931\n", "PI068 11.240555\n", "PI069 10.846550\n", "PI070 10.381338\n", "PI071 10.801455\n", "PI074 10.470945\n", "PI077 10.007513\n", "PI078 9.718330\n", "PI079 9.230219\n", "PI081 10.813904\n", "PI084 10.438283\n", "PI094 10.511939\n", "PI097 10.183464\n", "PI111 10.866402\n", "PI114 10.701429\n", "dtype: float32" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model = FeOi_predict.model\n", "model(wholerock)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The *random_sample_coefficients* method randomly samples fitted coefficients within their errors and calculates the matching x-intercept value. This method is used internally in the Monte Carlo PEC correction model to propagate FeOi prediction errors." ] }, { "cell_type": "code", "execution_count": 10, "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", "
TiO2Al2O3CaOintercept
01.555664-0.186890-0.12927210.195312
11.704102-0.263428-0.22900411.937500
21.424805-0.151489-0.27514611.398438
31.440430-0.228882-0.12719711.226562
41.457031-0.237793-0.17614711.789062
\n", "
" ], "text/plain": [ " TiO2 Al2O3 CaO intercept\n", "0 1.555664 -0.186890 -0.129272 10.195312\n", "1 1.704102 -0.263428 -0.229004 11.937500\n", "2 1.424805 -0.151489 -0.275146 11.398438\n", "3 1.440430 -0.228882 -0.127197 11.226562\n", "4 1.457031 -0.237793 -0.176147 11.789062" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "FeOi_predict.random_sample_coefficients(n=5)" ] } ], "metadata": { "kernelspec": { "display_name": "py310", "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.10.17" } }, "nbformat": 4, "nbformat_minor": 2 }