{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "613f0b35",
   "metadata": {},
   "source": [
    "# Model comparison"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0141b07",
   "metadata": {},
   "source": [
    "In this notebook we compare PEC corrections calculated with different ol-melt Fe-Mg Kd models. These codes can be use to quickly identify how different models change PEC correction results."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42c08d54",
   "metadata": {},
   "source": [
    "import the relevant packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4c281c68",
   "metadata": {},
   "outputs": [],
   "source": [
    "import MagmaPEC as mpc\n",
    "import MagmaPandas as mp\n",
    "\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import geoplot as gp\n",
    "\n",
    "# this is needed to clean up printed results during looped calculations\n",
    "from IPython.display import clear_output"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c94d92da",
   "metadata": {},
   "source": [
    "load data and set up FeOi according to the [FeOi](https://magmapec.readthedocs.io/en/latest/notebooks/FeOi.html) and [PEC correction](https://magmapec.readthedocs.io/en/latest/notebooks/pec_corr.html) examples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e1a90f4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "melt_file = \"./data/melt.csv\"\n",
    "olivine_file = \"./data/olivine.csv\"\n",
    "pressure_file =\"./data/pressure.csv\"\n",
    "wholerock_file = \"./data/wholerock.csv\"\n",
    "\n",
    "melt = mp.read_melt(melt_file, index_col=[\"name\"], units=\"wt. %\")\n",
    "olivine = mp.read_olivine(olivine_file, index_col=[\"name\"], units=\"wt. %\")\n",
    "pressure = pd.read_csv(pressure_file, index_col = [\"name\"]).squeeze()\n",
    "wholerock = mp.read_melt(wholerock_file, index_col=[\"name\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "26079a0b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAAFaCAYAAADVfgw3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABIp0lEQVR4nO3deVhUZf8/8PcABq5g4IIijLK6goCFiiCgJriEiru4YfmISVq5PJqplZn1M79q7oqmloklLiiioiiglhsZ+bgzIGouWAMiwzb37w+aiWFmYAaYOefA53VdXJczc5bPPYPvOdznnPsWMcYYCCGECIIJ1wUQQgjRHYU2IYQICIU2IYQICIU2IYQICIU2IYQICIU2IYQICIU2IYQICIU2IYQIiBnXBRiKXC7Ho0eP0LRpU4hEIq7LIYQQJcYY8vLy0KZNG5iY6HfsXGdD+9GjR2jXrh3XZRBCiFYPHjyAnZ2dXuvU2dBu2rQpgLI3pVmzZhxXoz/GGAoKCtCwYUPB/qVAbeAHobdB6PUD6m3Izc1Fu3btlDmljzob2ooPt1mzZoIN7QYNGgj+F5XawD2ht0Ho9QPa21Cd9tTZ0NakuLgYpaWlXJehE8YYCgsLIRKJBP2LWrENpqamaNCgAceVESJc9Sa0i4uLce/ePQhpUEO5XK73SQq+qdgGkUgER0dHCm5CqqnehHZpaSkYY2jTpg3Mzc25LqdKjDFl4An5SLt8GwoLC/Ho0SOUlpZSaBNSTfUmtBXMzc1hYWHBdRlVqouhTQipOWH/7U20mj59OtclEEIMgELbSORyuVHWUdi8eXO11yWE6E4ikWBYcB9IJBKj7I9C24CSkpIwePBgDB06FLt370afPn3Qq1cv7N27FwBw7do1eHt7Y+jQoXj77beRlJSkXOftt9/Gnj17cPz4cbX1Pv74Y/Tq1Qt+fn64ePEiLly4gDfffBP+/v745JNPAADe3t4Ayq5TDwwMRJ8+fTBjxgwAwM6dOzF8+HAMGjQIPXr0wKNHjzh4dwgRPolEgqkjAxDlnoqpIwOMEtz1rk+7vMePH+Px48cqzzVv3hzt27eHTCbDjRs31Nbx9PQEANy6dQvNmjWDra1tpfvIzc3FmTNn0KdPHyQlJcHMzAwBAQEYNWoUPv74Y/zwww9wdnaGv7+/yjpJSUkoLS2Fv7+/2noJCQm4cOECzMzMIJfLsWTJEixevBiDBw9WOzr/8ssvMW/ePAwcOBARERE4e/YsAMDKygrR0dHYunUrfvrpJ0RFRVXrPSSkvlIE9rZRmRDbiOBgnYmpIwMQvf8MxGKxwfZbr4+0N2/eDC8vL5WfxYsXAwCys7PVXvPy8lKuO3nyZJ26ILy9vZGTk4M7d+5gwIABCAwMxPPnz/Hs2TM8ffoULi4uEIlEKtv29vaGSCTC8+fPNa73+eefY/r06Zg+fTqePn2KmTNn4uTJk5g4cSKOHz+usv979+6hR48eAIAePXrg7t27AIDu3bsDANq1a4e//vqrZm8kIfVMxcAGALGNCNtGZRr8iLteH2lPnz4dQ4cOVXmuefPmAAA7OztcuXJF67o7d+7U6U5LExMT2NjYoGPHjjh58iQaNGiA4uJiNGjQAK1atcKdO3fg5OSEq1ev4u2331auA0DrepaWlnjrrbfwww8/YMuWLZg7dy7WrFmDoqIieHl5ISQkRLl/JycnXLp0CQMHDsSlS5cwadIk3L9/X+VqDiFdu04I1zQFtkL54DbUEXe9Dm1bW1ut3RsWFhbKrhBNXF1ddd6PiYkJFi1ahH79+sHExAQtWrRATEwMPv30U4wdOxatW7dGkyZNlMFcfr2FCxeqrRcaGoqCggIUFhZi27Zt2Lx5Mw4cOID8/HxMnjxZZd/z58/HpEmTsHz5cnTp0gV+fn64f/++zrUTQlTNmRGOxUHqga0gthFhcVAm5swIR2x8cq3vX8Tq6GFWbm4uLC0tIZVK0axZM8hkMmRkZKB9+/a8uU5bceQsl8sRGBiIvXv3Kr9E6sI1zhXbwMfPoCp1cbAioeFb/ZUdaQOA5DnDtBgHlSNtTQNGlc8nfdTrPm2u/fLLL/D398ebb76Jfv36VXlSkxDCPbFYjOj9ZzAtxgGS56rHvJoCu7bV6+4Rrvn6+iqv5iCECIciuMsfcRsjsAE60iaEkGopf8R99pZxAhug0CaEkGpTBPfa33obJbABCm2tjH1rKiFEmMRiMWLjk40S2ACFtkZc3JpKCCG6oNCuoPzlPP6uhr/DaenSpYiLi8Off/6JJUuWAAD69u2Lly9f6r2tnTt3oqioSPnvCxcu1EqNPj4+Wl9LSkrC7du3a2U/hJCqUWiXw+Wtqa1bt8ayZct0Wlbb6H/lQ3vy5Mno2bNnrdWnDYU2IcZFof0PXW9N1TW4GWOYOXMm+vTpA39/f7x48QK7d+9GYGAgPD09sXv3brX9h4WFKR8vW7YMAQEBmDVrFoCyQB49ejQGDRqEU6dO4cMPP0Tfvn3xxhtvIC0tDRcuXEBaWhqCg4OxZs0a5RE8AMyZMwe+vr4ICAhARkYGAKBjx44YP348unfvrlYLAHzxxRfo2bMn3nvvPeW8mhXrLygowM6dO/Hf//4XU6ZMwZMnTxAUFAQ/Pz+EhYUJZj5OQgSF1VFSqZQBYFKplDHGWEFBAbtx4wYrKCjQuHzoQF92er6IlUabaP05PV/EQgf66rT/gwcPsqioKOXj0tJSlp+fzxhjTCaTMW9vb8YYY0uWLGFHjhxhGRkZbMSIEYwxxvz9/dnx48dZSUkJGzduHLt8+TLbsWMHmzBhgnJ7im399ttvbNy4ccr18vLyVLb766+/stGjRzPGGDt37hybPHkyY4wxKysrJpVKWV5eHvPx8VGp/fHjx8zPz4/J5XJ269Yt1r59e5V9aqqfMcYKCwtZcXExY4yxOXPmsISEBFZSUsLkcrlOnwEfyeVylp+fr2yDEAm9DUKvnzH1NlTMJ33QzTX/WL1xN6aODICDtfZbUz9LdED0fvWjUk1u3ryJPn36KB+bmJjg5MmT+OabbwCgyi4Fxah/5UfmU4zWBwCrVq1CQkICTExMYGpqqnU7FUf5W7hwIQCgQ4cOyttnWYWRDCQSCbp16waRSAQXFxdYWloCQJX1v3jxAv/5z3/w119/4fHjx+jWrVulbSSE6I+6R/5R27emduzYEampqcrHjDEsW7YMhw8fRkJCApo0aVLp+teuXQMAXLlyBU5OTgD+Hf0vJycHcXFxSE5OxrfffqsM3QYNGqh1SShG+QOAS5cuwdnZGQAqHcNBLBYjPT0djDHcvXsXUqkUADTWX36f33//PQYMGICzZ89i8ODBNHogIQZAR9rl1OatqUOGDEF8fDx8fX3x2muvYf/+/RgxYgQCAgLg4eGhHAJWm/j4eHz66adwd3eHl5cXfv/9d+VrzZs3R6tWrRAQEKByZcfQoUMxatQojBo1Svmct7c3bG1t4evrCzMzM+zYsaPK2lu3bo3+/fujZ8+e8PT0hLW1NQBorD8wMBDz58/H2bNnMXHiRISHhyMhIQGNGjVC165ddX6/CCE6qsVumyrNmjWLOTg4MADs999/Vz4/ZcoU5uLiwtzd3VmfPn3YtWvXtG7jyJEjzNXVlTk6OrLhw4cr+3Ar0rdPu7yMjAwW4C1mp+eLWIC3mGVkZOjVztogl8tV+oOFqGIbqE+bG0Jvg9DrZ6x2+7SN2j0SFhaGlJQUODg4qDwfGhqKP/74A2lpaZg3b57KkWJ5L1++REREBA4ePIi7d+/C1tYWy5cvr/U6ubg1lRBCdGHU7hE/Pz+Nz5efPcbHxweZmZnKcZjLi4+Ph7e3N9zc3AAAkZGRCAkJwYoVK7TuUyqVgjGGoqIiyOVyyOVynfpaHRwccODYOQDcz+zC9f5rA2NM2Y7y/+a78jULldDbIPT6AfU21KQtvOvTXrNmDUJCQtQCGwCysrJUjtLFYjEePnyoMeAV7O3tAZSF8KZNm9CkSRPY2dkZpngD0HYjjZAo2qD40pTJZIL7DyiTybguocaE3gah1w/824aCgoJqb4NXob1nzx7ExMQgOVn7FD36zlyRlZWFZs2aoaioCE+fPlX+Wwgq+zISivJtKC4uhomJCSwsLAQ1c41MJoOFhQUvZk2pDqG3gc/1S6VS5SWxlanYhvLTCuqLN6G9b98+LFu2DImJiWjZsqXGZezt7XH69GnlY4lEgrZt21YabJaWlmjWrBmKi4uRk5ODP//8s9ZrN5S6FtpA2ZeumZkZ7/7zVUUkEgmu5oqE3ga+1b9161YsXLgQly9fVjtPp42iDTVpBy9COyYmBh9//DFOnTql7M7QZODAgZg5cyZu3rwJNzc3bNiwAWPGjNFpHw0aNICjo6Ngbq3m89GFrjS1wdTUFA0aNOC4MkJqZuXKlViwYAFmzZqFdu3aGXfn1b2EpToiIyNZ27ZtmampKWvVqhVzdHRkjDFmZmbG7OzsmLu7u/Ln+fPnjDHGFi9ezDZu3KjcxqFDh5SX/IWGhmq9ZKYml9TwQV28zEmIqA3c41P9crmczZ07lwFgS5Ys0bmm2rzkr97Mxi40jGczUFcHtYEfhN4GPtUvkUjg4eGBZcuW4f3339d5vYptqEk+8aJ7hBBC+KywsBCMMYjFYty9exc2Njac1SLss1yEEGJgL1++xODBgzF16lQA4DSwATrSJoQQrXJycjBo0CDcuHEDhw8f5rocABTahBCi0aNHjzBgwAA8efIEZ86cUQ6XzDUKbUII0WD37t2QSqVITk5WDp3BB9SnTQgh5eTl5QEA5s2bh6tXr/IqsAEKbUIIUUpNTUWHDh1w6tQpiEQitGjRguuS1FBoE0IIgOPHj6N///7o1KkT3njjDa7L0YpCmxBS7+3btw9Dhw5FUFAQjh8/zusb8ii0CSH1WlFREZYuXYrRo0fjwIEDaNiwIdclVYquHiGE1EuMMeTn56NJkyY4d+4crK2tBTGqJv8rJISQWsYYw7x589CzZ0/IZDK0aNFCEIEN0JE2IaSeKSkpwfTp0xEdHY01a9YIZkIOBQptQki9UVhYiHHjxuHQoUPYtWsXwsPDuS5JbxTahJB649y5c4iPj0dsbCyGDBnCdTnVQqFNCKnz8vPz0ahRI/Tv3x/3799H69atuS6p2oTR804IIdX08OFDvPHGG/jmm28AQNCBDdCRNiGkDrtz5w769+8PuVyOwYMHc11OraAjbUJInZSWlgZfX180bNgQqampcHV15bqkWkGhTQipk1asWIF27drh3Llzxp8x3YCoe4QQUqfk5+ejcePGiI6ORmlpKa/HEakOOtImhNQZe/fuhZOTEzIyMtC4ceM6F9gAhTYhpI7YuHEjxo8fjwEDBtSp7pCKKLQJIYLGGMPy5csRGRmJqKgo7NixA2Zmdbfnl0KbECJo2dnZWLlyJT799FOsXr1aMAM/VVfd/ToihNRpJSUlKC0tRbt27XDz5k20adOG65KMom5/JRFC6iSZTIaRI0di8uTJAFBvAhugI21CiMDk5eUhNDQU58+fx08//cR1OUZHoU0IEYycnBwEBwfj1q1bOHHiBPr06cN1SUZHoU0IEYwdO3ZAIpEgKSkJ3bt357ocTlCfNiGE9woKCgAAH374IdLS0uptYAMU2oQQnktLS0O3bt1w4sQJiESienXSURMKbUIIbyUnJ6Nv375o3bo1PD09uS6HFyi0CSG8dOzYMQwYMADe3t44duwYbGxsuC6JFyi0CSG8U1JSgo8++ggDBw5EXFwcmjZtynVJvEFXjxBCeEUmk8HCwgKJiYlo0aIFTE1NlSciiZGPtKOioiAWiyESiZCenq58/osvvoCrqytMTEwQFxendX2JRAIzMzN4eHgof+7du2eM0gkhBsYYw+effw4fHx+8evUKtra2dXrgp+oyamiHhYUhJSUFDg4OKs8HBQXh2LFj8PPzq3IbVlZWSEtLU/44OjoaqlxCiJHI5XJ88MEHWLx4MUaOHImGDRtyXRJvGfVrTFsov/nmm8YsgxDCIyUlJZg2bRp27dqFDRs2YMaMGVyXxGuC+9sjNzcXPXr0QGlpKUJDQ7Fo0SKYmppqXV4qlYIxpnxsbm4Oc3NzY5RaI4qay9cuNNQGfuB7G86ePYu9e/diz549GDt2rFqdfK9fFxXbUJO2CCq0bW1tkZ2djZYtW+LFixcYPXo0Vq1ahXnz5mldx97eXuXxwoULsWjRIkOXWmtkMhnXJdQYtYEf+NYGmUwGc3Nz9OzZE9evX0e7du0qPeHIt/qrQ9GGmpxYFVRom5ubo2XLlgCA119/HVOnTsUPP/xQaWhnZWWpzBMnpCNtxVl0kUjEdTnVQm3gBz624dmzZwgJCcHo0aPx0UcfwcXFReuyfKxfXxXbUFxcXO1tCSq0nz59iubNm6NBgwYoLCzEgQMHqhyDwNLSUtCTe4pEIsH+oipQG/iBL2148OABBgwYgBcvXmDAgAE618SX+mtC0YaatMOoV4/MnDkTdnZ2yM7ORr9+/eDk5AQAWLFiBezs7HDhwgVMnjwZdnZ2ePbsGQDgk08+waZNmwAAKSkp6N69O9zd3eHp6YnWrVsLqquDkPru9u3b8PX1RUFBAVJSUuDu7s51SYIjYkLu3a9Ebm4uLC0tIZVKBXmkzRhDQUEBGjZsKNijC2oDP/CpDePGjcNvv/2GEydOoG3btjqtw6f6q6tiG2qST4LqHiGECFNhYSHMzc2xZcsWFBYWwtramuuSBIvGHiGEGFRcXBycnZ1x7949NGnShAK7hii0CSEG8/333yM0NBTe3t46d4eQylFoE0KUJBIJhgX3gUQiqfG21q1bhwkTJmDixImIiYmBhYVFzQskFNqEkDISiQRTRwYgyj0VU0cG1Ci4Hz9+jIULF+LDDz/E9u3baeCnWkTvJCFEGdjbRmVCbCOCg3Umpo4MQPT+MxCLxTpvRy6Xo7S0FLa2tvj999/h4OAg2Cs++IqOtAmp5yoGNgCIbUTYNipTryPu4uJiTJ48GeHh4WCMKYdhJrWLQpuQekxTYCvoE9wFBQUYMWIE9u7di9DQUAprA6LQJqQemzMjHIuD1ANbQWwjwuKgTMyZEa51G7m5uQgODsapU6dw+PBhjBkzxlDlElBoE1Kvrd64G58lOkDyXPON0ZLnDJ8lOmD1xt1at7F9+3akpaXh5MmTCA4ONlSp5B8U2oTUY2KxGNH7z2BajHpwS54zTItx0HoysqioCADw/vvvIy0tDb179zZGyfUehTYh9Zym4K4qsG/duoWOHTvixIkTMDEx0esKE1IzFNqEEJXgPnur8sC+cuUKfH190bBhQ3Tu3Nn4xdZzFNqEEAD/Bvfa33prDeykpCQEBATA0dER586do1vTOUA31xBClMRiMWLjkzW+VlpaipkzZ+LNN99EbGwsmjRpYuTqCEBH2oTwRm2O+1HbioqKYGpqioSEBMTFxVFgc4hCmxAeqM1xP2rbmjVr0LNnT7x8+RJ2dnaCmGO1LqPQJoRj5e9K9HfV//ZxQ2GMYenSpZg9ezaCgoLQuHFjTushZSi0CeFQbY37Udvkcjnef/99LFu2DF9++SW++uorujWdJyi0CeFIbY37YQgpKSnYsGEDNm/ejPnz5xt9/0Q7unqEECOTSCSYMyMcskKZzuN+aLuio7YVFxejQYMG8PPzw//+9z84OzsbZb9Ed3SkTYgRlT/h+Pfzx1iS0LZG437UtJbyV6tIpVL069cPq1atAgAKbJ6i0CbESCqecPw+/BGKSxkm/6Ae3FXdRl5btUS5p2KQnxsOHz6MgIAAXL9+HT179qz1/ZHaQ6FNiBFoO+GoKbiNFdiKL48js4rw0btvQyKR4OzZs+jVq1et75PUHgptQgysqhOO5YO7qnE/DFGL2EaE4x+K0KphLnJzc2t9n6R2UWgTYmC6TDTwRcgjWDRrXem4HzVV1ZfH0dlyvDM6ACkpKbW+b1J7KLQJMTBdJxrY8t1+xMYnG2yY0zmRVX95bAgvxeSR/ZQnJ/l8a319RaFNiIHVZKKB2rR6w26898NrlX95HGL4LLQQc2aE8/rW+vqMQpsQI6jORAOGqKHvkAgMXMU0f3lEMywNBbZeFuPDRSt4eWs9odAmxGj0mWjAEK5cuYL4/Rvx2TBg0Gqm+uXxT2BP32WGGfO+widzwnl3az0pQ6FNiBHpMtFAbWOMoaSkBJ8vno3tU4GRb5hg86Sy4D5769/AXnoQmNG3BEs/DOflrfWkDIU2IUammGjAGIEtl8sxa9YsTJs2Dd+s36U8IerrUhbc70QzvNO3LLCXhgJrT5vj23FFOt1aT7hBoU1IHVVcXIzw8HBs2rQJ/v7+aN++vUq/uq+LCU7MFWHfL2WBvfSUGDv3n9LpShdD3VpPqkahTUgd9OrVKwwbNgz79+/Hjz/+iClTpgBQPyEqthHhm7EiLD1V9ryvry8vrnQh2lFoE1IHbdu2DUlJSTh69CjCwsJUXqvqhCgfrnQh2lFoE1KHlJSUAADee+89XLlyBf3799e4XFUnRLm+0oVoZ9TQjoqKglgshkgkQnp6uvL5L774Aq6urjAxMUFcXFyl24iLi4ObmxucnJwwYsQIvHz50tBlEyIImZmZ8PDwwIkTJ2BiYgJXV9dKl6/shKhizO9PV+826pUupGrVCu3bt2/j9OnTuHDhAvLy8nReLywsDCkpKXBwcFB5PigoCMeOHYOfn1+l6798+RIRERE4ePAg7t69C1tbWyxfvrw6TSCkTrlx4wZ69+6NgoICODk51Whb5e+E/GROOFZv3E2BzSM6z1yTl5eHb775Btu2bYO5uTlatWoFmUyG+/fvw8fHB3PnzkVgYGCl29AWym+++aZONcTHx8Pb2xtubm4AgMjISISEhGDFihW6NoOQOufSpUsIDg5GmzZtkJCQAFtb22pvq+KgUg7WZddl05E2f+h8pB0QEABLS0tcunQJd+/eRWpqKq5cuYKcnBwsWLAAW7duxZYtWwxZK7KyslSO0sViMR4+fAi5XK51HalUqvIjk8nAGBPEDwDOa6A2GL4NGRkZCA32RUZGht7Ll5aWIiIiAi4uLkhKSkLr1q31akP5bWVkZFQ6ybCu9QnxM+CiDdWl85F2amoqzM3N1Z43MTGBv78//P39UVhYWO1CdKXvjND29vYqjxcuXIhFixbVZkkGJZPJuC6hxqgN2mVmZiIyPBif9MvC1JEB2LA7Xq37UNvyU0YGYOPueMTExMDa2hoWFhYoKCjQum7FNpTfVviwPigpBX6Y+KjSOyGrqs+Q6tLvUWWfU1V0Dm1NgS2RSJCfn4/OnTtrXaY22dvb4/Tp0yr7b9u2LUxMtP/BkJWVhWbNmikfm5ubG7zO2sAYg0wmg4WFhd5fVHxBbaicRCLBzIkh2D46C2IbEbZbZ2HaxBBsjzmt9eRg+eWjrbMwNTwYO/afgY2NjV5tqLitXdaPMH4zA6D9TshP+mXhvx+8g9hjxplkWKEu/h4VFxdXe1vVno193bp1iI2NhUgkgouLCzZu3FjtInQ1cOBAzJw5Ezdv3oSbmxs2bNiAMWPGVLqOpaWlSmgLjUgkEuwvqgK1QV1KSgomj+yH6MmFENuUHXQojmgjRgWq9SFLJBJEjApU67qIHp2lcfnK2qBtW99PByZtY/huGtSOthV3Qkbv38PZZ1mXfo9q0g6d+7QvXLig8jglJQWnT59GYmIiUlNTddrGzJkzYWdnh+zsbPTr1095lnvFihWws7PDhQsXMHnyZNjZ2eHZs2cAgE8++QSbNm0CADRt2hTbtm1DaGgonJyc8PDhQyxcuFDXJhDCCykpKXhndAC2TirC0oNQufNQ06BMVc04o2kQJ22TF1S1re+miTBpm+rQrXRjDb+ImI494lOmTEGTJk3w5ZdfonHjxnjvvffQsmVLiEQiJCYmIikpycCl6ic3NxeWlpaQSqWCPNJmjKGgoAANGzYU7NEFtUGdRCLBID83HJlVNiiTYljUbVNFKiF69hbD2t96IzY+GcOC+yDKPRX+rtr3X355RTAvDsrEZ4kOWPbNLnz9+QKs2bQHH0RO1GlbXx1jODrHhBeBXRd/j2qSTzofae/YsQNDhgxBcHAwjhw5gtWrV8POzg5WVlaIjY3VuxGE1DeKMFUENvDPkfJUEaZFq45vXX5Qpv/37U6dB3GqONP6tlGZeHdMIMZ2OI+IUYH4cNGKKre1JKEt/obhJxkm1aPXzTUDBgxAQkICzp07h0mTJmHw4MGYNWsWmjdvbqj6CNELX+c0rLKL45/gTrktVwnKV69eISoqCr3eGl/lIE4ANF6yd3R2KbYmAUv7SfDJnLK7HCvb1q7YZOw9nEJ3QvIV01FWVhabO3cuW7RoEXv69Cm7fPkyCwwMZNu3b9d1E0YllUoZACaVSrkupVrkcjnLz89ncrmc61KqzdhtyMjIYAHeYnZ6vogFeItZRkZGjbdZW20IHejLTs8XsdJoE60/p+eLmGNrc2Xdf/31F/P19WWNGzdmJ0+eVLbv3ldl27n31b/trPhaxZ97X4lYgBvY2QVgAd5ilpycrHVbfFMX/y/UJJ90PtIePXo0unbtCjs7O0yaNAleXl44ceIEnj9/jkGDBhnuW4UQHWjqFuDTDCu6zMgeudsUO/efglgsxpMnT9C3b1/88ccfSExMRL9+/SodxGnOjKpnWl/8tgjfJADbRmWqHHFTN4jA6JrunTt3ZowxJpPJmIeHh8pr9+7d0/vbwtDoSJt7xmqDtqPM2jh6rI02ZGRksNCBvmpHt+Xr7NTOnCUnJyvXmTJlCmvTpg1LT0/Xur3y7dL1SFvx+un5IuU2Km6Lb+ri/4Wa5JPOof3pp58yV1dX5uzszDZu3Kj3joyNQpt7xmiDTmFVg+CuaRsqdtlU1S1RWlrKGCv7/dW35kq/vMoFNp+7QjSpi/8XapJPOl/yB5QNGmVqaopGjRoZ6sC/1tAlf9wzRhv0vRxOXzVpQ8WTj4oTfZ+u3o1P5oQrL8lTdEv8+uuvmDJlCg4fPgxHR0e9a9W6z3KXFPLhEj591cX/C0a55A8ou7lFCIFN6g9d+oq5mNNQ09UiihthFP3J5a/OSExMRGBgIJo3bw5ra2ut26zqypiK/d6D/s8US0Mh2MAm6mplEgRPT8/a2AwhetM0NZYCVyGlyx2M5cepjo2NRUhICPr06YMTJ07AyspK6zaj3FOrPMFaflaarfvOYOkpMZ1srEP06h7R5vHjxzUaw9cQqHuEe8Zsg7auiJqGVHXaoE+XzbY9B9G+fXuEhIRg165deO2119SWrWnbMjIy8P5/JmDNpj1o3769Tm3gk7r4f8Fo3SMAEB0drfZcfHy8vpshpFbxaU5DXbtsVq3/DtbW1khNTcX333+vU2AD2scb0UYsFuPH2JN0hF1H6B3a3377rU7PEWJsVU1Wa+w6Kuuy6fhGCJYsWQLGGLp27QpTU1O17VRnoChS9+k8NOvly5fxyy+/4Pnz59iwYYPyealUiqKiIoMUR4i+FJPVck0R3BW7NSJiHNDGpRc2bNiAr7/+utI/93W6YSYoE3NmhPOizcQ4dA7thw8f4vLly8jPz8elS5eUzzdr1gw7d+40RG2ECFr54F4clIlPEx3QpHUX/Pjjj4iOjsaUKVMqXX/1xt2YOjIADtaag/vfMa6Ne2UM4ZbeJyLj4+MRHBxsqHpqDZ2I5B61oYxEIsGcGeHw6v0WPvvsM/z4448YNmyYzutq6iLR52Sk0D8HodcP1O6JSL1nrlm1ahWuX7+Ofv36wcvLS9/VCal37O3tERufDLlcjiFDhsDd3V3ndbV1s9Dle/WX3iciFy5cCKlUihkzZqBly5YYOXKkwWdhJ0So/vzzT/j4+ODkyZMwMTHRK7AV+HRlDOGe3qEdGBiIL774AqdOncJXX32Fy5cvIyoqyhC1ESJoGRkZ8PX1xcOHD9GmTZsabYsvV8YQ7undPbJ48WKcPn0aBQUFCAgIwPr16+Hv72+I2ggRrPT0dAwYMACNGzdGampqrYQsX66MIdzSO7Q3bdoEFxcXvPPOO+jfv79ycl5CSBnGGCZOnIiWLVsiISEBrVq14rokUofoHdrPnj1DWloaTp06haioKGRmZqJXr17YunWrIeojRFDkcjlMTEywf/9+WFtbaxxHhJCaqNaAUTY2NrC2tkbz5s2Rk5Ojct02IfXVxo0b4dTOCunp6XB0dKTAJgahd2i7ubnB19cXKSkpGDRoEK5fv460tDQDlEaIcKxcuRJffRKJ7eEvMWvKELq1nBiM3t0jhw8fhouLiyFqIYQTiptfFEOl6uu///0vftz2JRLnlU00sN26bEwQusqDGILOR9oPHjwAgEoD+9GjRzWviBAj0mecak0OHTqkEtgADeZEDEvn0B47diymTZuGM2fOQCaTKZ+XSCTYsGED3nzzTVy8eNEgRRJiCDWZwZ0xBolEgjWfz1YJbAUKbmIoOod2SkoK3nrrLSxfvhyvv/46mjdvjoYNGyIwMBAZGRnYv38/hg8fbshaCak1NRmnuqioCKNHj8aotwN1HoWPkNqiV5/2yJEjMXLkSJSUlOD58+do1KiRIAdjIvWbruNUR+8/AwcHB5XX8/PzMXz4cCQlJWHdunX4bOsKGoWPGFW1LvkzMzND69atKbCJIOkzTnV5L168QP/+/XH+/HnEx8fj3Xff5d38lKTuq5WJfQkRkurO4L5gwQLcvn0bp0+fRmBgIADNs9RQYBNDotAm9Y6+M7grhpz/+uuvkZqaih49emjdHo3CRwyNQpvUS7oeIaenp8PX1xf37t2DpaUlXF1dK90ejcJHDE2vE5FPnjzBhg0bcPXqVQCAp6cnZsyYgdatWxukOEIMqeJ0YGUnDf8N3AsXLiAkJAT29vZo0qSJTtujUfiIoel8pH3z5k1069YNN2/eRFBQEAIDA3Hz5k24u7vj5s2bhqyREIPRdoR84sQJ9OvXD127dsXx48dppD7CGzrPETlkyBCMHz8eY8aMUXl+79692LNnD44ePWqQAquL5ojknlDb8Pfff0MsFqN3796IiYmBSCQSXBvKE+rnoCD0+oHanSNSryPtioENlN0pefv2bb12SghfMcZgZWWFxMREHDx4EI0aNeK6JEJU6BzalR2Q6zmhOyG8tHLlSkyZMgVyuRxeXl5o0KAB1yURokbn0HZzc8O+ffvUnv/xxx91HvUvKioKYrEYIpEI6enpyuefPn2KgQMHwtnZGV26dEFKSorG9SUSCczMzODh4aH8uXfvnq5NIEQjxhjmz5+PBQsWwN7eXrB/gpP6QeerR77++mv07dsXsbGx6NmzJ0QiEVJTU5GUlISkpCSdthEWFoZ58+bB19dX5fkFCxbAx8cHx48fx6VLlxAWFoZ79+7BzEy9PCsrKxq/m9Sa0tJSzJgxA1u3bsXq1asxe/ZsrksipFI6h3bHjh2RlpaGDRs24MSJEwDKLvlbvXq1zjNN+/n5aXw+JiYGGRkZAIAePXqgVatWSElJQd++fXUtTyupVKrSfWNubg5zc/Mab9fQFDULuetJCG2Ijo5GdHQ0duzYgUmTJqnVKoQ2VEXobRB6/YB6G2rSFr2u07a1tcVnn30GACgpKdF4JKyvnJwcyOVytGjRQvmcWCxGVlaWxuVzc3PRo0cPlJaWIjQ0FIsWLYKpqanW7dvb26s8XrhwIRYtWlTjuo2l/DC4QsXHNjDGIBKJMGbMGLi4uKBHjx4oKChQWy4zMxPzZ0/Dyv/bpjZ4lNDw8XPQh9DrB/5tg6bfNV3pnbo3btzAuHHjkJOTgwcPHuDKlSuIiYnBypUrq11ExT5Ebd9Ctra2yM7ORsuWLfHixQuMHj0aq1atwrx587RuOysrS+WSGiEdactkMlhYWAi2j5WvbXjx4gVGjBiBjz/+GEFBQVr/ApRIJJg5MQSLgzIxc2IItsecFuSdjnz9HHQl9PoB9TYUFxdXe1t638Y+c+ZMfPvtt7CxsQFQ1kVSk2u0ra2tAZTN8q6QmZmpdoQMlAVuy5YtAQCvv/46pk6diuTkyu9As7S0VPlRvGlC+AHAeQ11rQ2PHz9G3759kZ6eDisrK63LZWZmImJUoMoECRGjApGZmcl5G+rC51Df6tfUhurSO7Tz8vJUTiSKRKIaXxo1cuRIrF+/HgBw6dIl/Pnnn2onK4Gyq0wU31CFhYU4cOAAunfvXqN9k/rj3r178PX1xd9//43k5GR4e3trXK4mEyQQYmh6h7aZmRmKi4uV3xTZ2dkwMdFtMzNnzoSdnR2ys7PRr18/ODk5ASi7Pvb8+fNwdnbG5MmTsXv3bmV/+SeffIJNmzYBKJs9p3v37nB3d4enpydat24tqP5pwh3GGMaNG4cGDRogNTUVHTt21LicrhMkUHATruh8G7vCnj178OOPP+L69euYOnUqdu3ahS+++ELj3ZJcotvYuceXNihOOt6+fRtWVlZ49eqV1tnXhwX3QZR7Kvxdtdd79hbD2t96C2ZwKL58DtUl9PoBjm5jV5gwYQIWLlyIsWPH4tWrV/juu+94F9iEKCQkJCAoKAh5eXlwcXHBq1evKp19vboTJBBiLDqHdmRkpPLfz549w8qVK/HVV1+hT58+BimMkJqKiYnBkCFD0LhxY5iZmek0+7q+EyQQYmw6h/bFixeV/162bJlBiiGktmzZsgVjxozB6NGjceDAATx58kTnk4s0hRjhs2oNGCXkO5NI3ZeWlobp06fjvffew3fffYeHDx/qfXKRphAjfKXzzTWFhYX43//+B8aYyr8VOnXqZJACCdGV4oSjh4cHzp07B19fX4hEIr1mXy9/clEsFmN7zGm8/58J2B6zhwKb8ILOV48oRufTuBGRCPfv36/VwmqKrh7hnjHbUFJSgunTp8Pd3R1RUVEqr1V2GR9QedcHfQ7cE3r9QO1ePaLzkTZdl0r4qrCwEOPGjcOhQ4ewc+dOtdfLzwVZMbipr5oIDc3GTgTt5cuXGDx4MI4ePYrY2FhMmDBB43J0cpHUFRTaRNAWLFiAX375BQkJCRgyZEily9LJRVIXUGgTQVKcivn888+RnJwMf39/ndbTNvs6IUJBoU0E5+7du+jTpw8yMjJgZWUFd3d3vdYXi8WIjU+mwCaCRKFNBOX69evw9fXF8+fPK538gpC6ikKbCEZqair8/PzQtm1bJCcnaxxznZC6jkKbCEJeXh6GDh0KDw8PnDlzRmV6OkLqk5pP8kiIgTHG0LRpUxw9ehTu7u5o2LAh1yURwhk60ia8tmnTJrzzzjuQy+Xw8fGhwCb1HoU24SXGGFasWIEZM2agcePGXJdDCG9QaBPeYYxh7ty5WLhwIZYtW4b/+7//03lKO0LqOurTJryze/durFq1CmvXrsWsWbO4LocQXqHQJrwzfvx42Nvbo2/fvlyXQgjv0N+chBfy8vIwePBgnDlzBqamphTYhGhBoU04l5OTg6CgICQnJ9NdjoRUgbpHCKeys7MxYMAAPH/+HElJSejevTvXJRHCa3SkTTjDGMPYsWORn5+PlJQUgwa2RCLBsOA+NJkHETwKbcIZkUiErVu3IiUlBS4uLgbbj2K6sSj3VLUJfAkRGgptYnTJyckIDg5Gfn4+3Nzc0K5dO4Ptq/z8kP6ummdeJ0RIKLSJUR09ehQDBgxAYWEhSktLDbovTRP6im0ouImwUWgTo/nhhx8QGhqKt956C8eOHdN7Fmp9VDYDOwU3ETIKbWIUf/zxByZMmIDx48fjp59+goWFhUH3N2dGOBYHqQe2gthGhMVBmZgzI9ygdRBS2+iSP2IUnTt3xvHjx9GvXz+jjCOyeuNuTB0ZAAdrzcEtec7wWaIDovfvNngthNQmOtImBiOXyzF//nysX78eADBgwACjDfxUfuZ1yXOm8prkOc3EToSLQpsYRElJCSIiIrB+/XqIRJq7KAxNU3BTYBOho9AmtU4mkyEsLAzff/89tm/fjsjISM5qKR/cZ29RYBPho9AmNaLpTsNFixYhISEBsbGxGD16NHfF/UMR3Gt/602BTQSPQptUm7Y7DT/++GOcOXMGgwYN4rbAcsRiMWLjkymwieBRaJNq0XSnYYhfR6SkpKB58+bw8fHhukRC6iSjhnZUVBTEYjFEIhHS09OVzz99+hQDBw6Es7MzunTpgpSUFK3biIuLg5ubG5ycnDBixAi8fPnSGKWTcrTdaRg3qxD/nTWWblghxICMGtphYWFISUmBg4ODyvMLFiyAj48P7ty5gx07dmD8+PEoKSlRW//ly5eIiIjAwYMHcffuXdja2mL58uXGKp+g6jsNvxv3kO40JMSAjHpzjZ+fn8bnY2JikJGRAQDo0aMHWrVqhZSUFLXZS+Lj4+Ht7Q03NzcAQGRkJEJCQrBixQqt+5RKpWDs3+t0zc3NYW5uXsOWGJ6i5vK188HsGRN0utNw9owJOHD0HAD+tUEffP0c9CH0Ngi9fkC9DTVpC+d3RObk5EAul6NFixbK58RiMbKystSWzcrKUjlKF4vFePjwIeRyudabNuzt7VUeL1y4EIsWLaql6g1PJpNxXYKKFd9sRWR4MLZbZ2m90/DTU/bYsHursna+taE6qA3cE3r9wL9tKCgoqPY2OA9tAGo3X1T2LaTvjRpZWVkqAxMJ6UhbJpPBwsKCs5tTNHFzc8PQ8e8jdN0SHJzxUiW4/71x5TTEYjFv26APagP3hF4/oN6G4uLiam+L89C2trYGADx79kx5tJ2Zmal2hAyUHTWfPn1a+VgikaBt27aV3hptaWlp0NHkDE0kEvHqF/Xbb7/FBx98gBEjRiAi5gq2/9O3XdmdhnxrQ3VQG7gn9PqBf9tQk3bw4pK/kSNHKsenuHTpEv7880/4+vqqLTdw4EBcunQJN2/eBABs2LABY8aMMWqt9RVjDJ9++ilmzZqFDz74ADExMdhBdxoSYnzMiCIjI1nbtm2Zqakpa9WqFXN0dGSMMfbnn3+y/v37MycnJ9apUyeWlJSkXGfx4sVs48aNyseHDh1irq6uzNHRkYWGhjKpVKpxX1KplAHQ+jrfyeVylp+fz+RyOdelMMYY27t3LwPAli9frlJTRkYGCx3oyzIyMtTW4VsbqoPawD2h18+Yehtqkk8ixgR8SrYSubm5sLS0hFQqFWT3CGMMBQUFaNiwIS/+JCwpKcGJEycQEhKi8zp8a0N1UBu4J/T6AfU21CSfeNE9QvhJJpNhzJgxOHfuHMzMzPQKbEKIYVBoE41yc3MRHByMw4cPIz8/n+tyCCH/4PzqEcI/z549Q3BwMO7evYsTJ05oPClMCOEGhTZRM2bMGGRnZ+Ps2bNwd3fnuhxCSDkU2kTN2rVr8dprr8HZ2ZnrUgghFVCfNgEAXL16FcOGDcOrV6/QuXNnCmxCeIpCm+Ds2bPo27cvHj16VCfGdyCkLqPQrkM0Tf1VlSNHjmDgwIF44403cOrUKbz++uuGK5AQUmMU2nWEtqm/KnP79m0MGzYMISEhOHr0KJo2bWr4QgkhNUKhXQdomvpLl+B2cXHBgQMHsG/fPkGMfEgIodAWPG1Tf2kLbsYYli1bhk2bNgEAhg4dCjMzuoiIEKGg0Bawqqb+qhjccrkcs2fPxtKlS/H3338bv2BCSI1RaAvYnBnhOk39NWdGOIqLizFp0iSsW7cOmzZtwoIFC4xcLSGkNlBoC9jqjbvxWaIDJM81D9Qoec7wWaIDVm/cjSVLlmDfvn3Yu3cvpk+fbuRKK1edq14Iqa8otAVMLBYj+p+JCCoGd8WZZD766CMkJCRg9OjRHFWrWXWueiGkPqPQFjhNwa0I7K83/YTZs2cjKysLr7/+OgICAjiuVlV1r3ohpD6j0K4Dyge3Yuqvz/5vD8aPH4+LFy9CKpVyXaIafa96IYSUodCuIxTBvfa33lj05TaMGTMGRUVFSE1NRdeuXbkuT4W+V70QQv5FoV2HiMVi7I09icmTJ8PS0hIpKSlwdHTkuiw1+lz1QghRRaFdx1hYWGDnzp04d+4c2rRpw3U5Gulz1QshRBWFdh1x+PBhvP/++2CMISgoiNcDP+lz1QshRBWFdh2wa9cuDB8+HI8ePUJJSQnX5eiksqteKLAJ0Y5CW+DWrFmDSZMmYcqUKfjxxx/RoEEDrkvSmaarXiiwCakchbaAHTx4ELNnz8a8efOwZcsWmJqacl2S3spf9UKBTUjVaHg3ARs8eDD279+PsLAwrkupEbFYjNj4ZK7LIEQQ6EhbYIqLi/HOO+8gJSUFZmZmgg9sQoh+KLQFpKCgAMOGDcN3332HJ0+ecF0OIYQD1D0iEFKpFEOHDsXly5dx5MgRvPXWW1yXRAjhAIW2QIwbNw7Xr1/HqVOn0LNnT67LIYRwhEJbIL788ksA4N04IoQQ46I+bR67efMmxo0bh1evXqFr164U2IQQOtLmq8uXLyM4OBi2trbIy8tDo0aNuC6JEMIDdKTNQ2fOnEFgYCAcHR2RlJSEVq1acV0SIYQnKLR55v79+wgODkavXr0QFxfH64GfCCHGR90jPNOhQwfs2rULQ4cOhVwu57ocQgjP0JE2T6xevRpbtmwBAIwaNQrm5uYcV0QI4SPehPbx48fh7e2Nbt26wcfHB7/99pvaMklJSWjUqBE8PDyUPwUFBRxUW3sYY/j444/xwQcf0PRahJAq8aJ75K+//sKECROQnJyMjh074uzZsxg/fjzS09PVlu3UqRMuX77MQZW1Ty6X47333sPGjRuxcuVKzJs3j+uSCCE8x4sj7Xv37qFly5bo2LEjAMDf3x+ZmZm4evUqx5UZ1meffYbNmzdj69atFNiEEJ3w4kjb2dkZz549w8WLF+Hj44PY2Fi8fPkSEokEnp6eKsveunULnp6eMDU1xZQpUxAZGVnptqVSKRj7d0orc3Nz3vQXR0ZGwsvLC4MGDVKpEYDyccXnhYTawA9Cb4PQ6wfU21CTtvAitC0tLfHzzz9jwYIFyMvLg6+vLzp16qQ2C4unpyeys7NhaWmJ7OxshISEwMbGBqNGjdK6bXt7e5XHCxcuxKJFiwzSDl38/fffeP/997F8+XLY2dkhMDCw0n55mUxmxOoMg9rAD0Jvg9DrB/5tQ03OxYkYD7++CgsL0bp1a1y6dAlOTk5al1uxYgUePXqEdevWqb2Wm5sLS0tLZGVloVmzZsrnuTzSfvLkCQYOHIisrCwkJCTA29tb67KMMchkMlhYWEAkEhmxytpDbeAHobdB6PUD6m3Izc2FlZUVpFKpSj7pghdH2gDw+PFj2NraAijr6w0MDFQL7MePH6NVq1YwMTFBXl4e4uLiEBERUel2LS0t9X5TDEEikaB///7Iz8/HuXPn0KVLF53WE4lEgv1FVaA28IPQ2yD0+oF/21CTdvDiRCQALF68GG5ubnByckJmZia2b98OAJg2bRoOHz4MAPj555/RtWtXuLu7w8fHB/3798eUKVO4LFsnhYWFCAwMhFwuR2pqqs6BTQghFfGye6Q2KLpHqvPnhyHEx8fDw8ND+ddEVRhjKCgoQMOGDQV7dEFt4Aeht0Ho9QPqbahJPvHmSLsuSkxMxEcffQTGmHLEPkIIqQkKbQOJjY1FSEgI0tPTUVRUxHU5hJA6gkLbAKKjoxEWFoZhw4bh8OHDvLkunBAifBTatSw+Ph4RERF455138P333+O1117juiRCSB3Cm0v+6op+/frhu+++Q3h4uGBPmhBC+IuOtGtBaWkp5syZg/Pnz6NBgwaYOHEiBTYhxCAotGuoqKgI48ePx9q1a3Hnzh2uyyGE1HHUPVIDr169QlhYGBITExETE4MRI0ZwXRIhpI6j0K6BiRMn4ty5czh69Cj69evHdTmEkHqAukdqYPHixTh16hQFNiHEaCi09ZSRkYEpU6agoKBAOQYKIYQYC4V2BRKJBMOC+2icr/GPP/6Ar68vkpOTkZOTY/ziCCH1HoV2ORKJBFNHBiDKPRVTRwaoBPcvv/wCPz8/2NjYICUlBXZ2dtwVSgipt+hE5D8Ugb1tVCbENiI4WGdi6sgARO8/AxMTEwQFBcHDwwNxcXGwsrLiulxCSD1FR9pQD2wAENuIsG1UWXDL5XJs2LABJ06coMAmhHCq3oe2psBWKB/cfn5+aNSoEUdVEkJImXof2nNmhGNxkHpgK4htRFgclIk5M8KNXBkhhKir96G9euNufJboAMlzzRP4SJ4zfJbogNUbdxu5MkIIUVfvQ1ssFiN6/xlMi1EPbslzhmkxDojefwZisZibAgkhpJx6H9qA5uCmwCaE8BGF9j/KB/fZWxTYhBB+otAuRxHca3/rTYFNCOElurmmArFYjNj4ZK7LIIQQjehImxBCBIRCmxBCBIRCmxBCBIRCmxBCBIRCmxBCBIRCm6cKCwuxfPlyFBYWcl1KtVEb+EHobRB6/UDttkHEGNM86IbA5ebmwtLSElKpFM2aNeO6HL1JpVJYWVnh77//hqWlJdflVAu1gR+E3gah1w+ot6Em+URH2oQQIiAU2oQQIiB19o5IRa9Pbm4ux5VUj6Lu3NxciESax/rmO2oDPwi9DUKvH1Bvg+JxdXqn62yfdnZ2Ntq1a8d1GYQQotWDBw/0niS8zoa2XC7Ho0eP0LRpU8F+OxNC6ibGGPLy8tCmTRuYmOjXS11nQ5sQQuoiOhFJCCECQqFNCCECQqHNIwsXLkTHjh3h7u6ON954A6dPn9a67C+//AIPDw+4uLggKCgIjx8/NmKl2kVHR6Nr164wMzPDt99+W+myIpEI3bp1g4eHBzw8PJCczI9xzPVpAx8/h1evXmHs2LFwcnKCi4sLDhw4oHVZPn0Gd+7cQa9eveDi4oI33ngDN27c0Ljc9u3b4ezsDEdHR7z77rsoKSkxcqXa6dKGpKQkNGrUSPmee3h4oKCgQPedMMIbx44dY69evWKMMZaWlsasrKxYQUGB2nJyuZw5OjqyM2fOMMYY+/rrr9mYMWOMWapWaWlp7MaNGyw8PJytW7eu0mUBsLy8PCNVpjtd28DXz2HZsmVs0qRJjDHG7t+/z1q1asVevHihcVk+fQYBAQFsx44djDHG9u/fz3x8fNSWuX//PrO1tWV//vknk8vlbMiQIWzTpk1GrlQ7Xdpw5swZ5uXlVe19UGjzVGlpKWvatCl78OCB2mu//vor69Spk/Jxbm4us7CwYEVFRcYssVKTJk0SbGgrVNUGvn4OnTp1Yr/++qvy8ciRI5VBUhFfPoMnT54wS0tLVlxczBgr+0Js1aoVy8jIUFnuq6++YpGRkcrHR48eZf7+/kasVDtd21DT0KbuEZ7asWMHHB0dNV7DmZWVBQcHB+Xjpk2bomnTprz401xfffv2hbu7Oz744APk5+dzXY5e+Po5VKxLLBYjKytL6/J8+AwePHiANm3awMys7H4/kUgEe3t7tbr1bZsx6doGALh16xY8PT3Ro0cPbNiwQa/91Nk7IvmoT58++N///qfxtWvXrilvBkpMTMSyZctw8uRJrduqeO05M9KVm7q2QReZmZmwt7dHfn4+/vOf/2Du3Ll6/wJXR222gYvPoar6K9ZVWU1cfQaa6Ppe6to2LujSBk9PT2RnZ8PS0hLZ2dkICQmBjY0NRo0apdM+KLSNSJeTPGfPnsWUKVNw5MgRuLq6alzG3t4eEolE+TgvLw95eXmwtbWtrVK1qs0TVfb29gCAxo0bIzIyEu+++26tbbsytdUGrj6HqupX1NWiRQsAZcEcEhKidVnA+J9BRe3atUN2djZKSkpgZmYGxhgePHigrE+h4nuu+NLhA13bUH5UPzs7O4wdOxbJyck6hzZ1j/DIuXPnEB4ejkOHDsHd3V3rcl5eXpDJZEhKSgIAbN68GaGhoWjQoIGRKq25v/76C69evQJQdvfqvn370L17d46r0g9fP4eRI0di/fr1AICMjAycPXsWQ4cOVVuOT59By5Yt0b17d+zZswcA8PPPP0MsFkMsFqssN2LECMTGxuLJkydgjGHTpk0YM2YMBxWr07UNjx8/hlwuB1D2RR8XF6ff+17t3nBS65ycnFjLli2Zu7u78uf69euMMcY2btzIFi9erFz2/PnzrFu3bszZ2Zn17duXZWdnc1W2it27d7O2bduyRo0aMSsrK9a2bVt29epVxphqG86fP8+6du3KunXrxjp16sQmTJjAcnJyuCxdSdc2MMbPz+Hly5ds1KhRzNHRkTk7O7P9+/crX+PzZ3Dz5k3m4+PDnJ2dmZeXF0tPT2eMMRYREcEOHTqkXG7Lli3M0dGRtW/fnkVERHB+4rc8Xdqwbt061qlTJ+X7vmTJEiaXy3XeB93GTgghAkLdI4QQIiAU2oQQIiAU2oQQIiAU2oQQIiAU2oQQIiAU2oQQIiAU2gIQFRUFsVgMkUiE9PT0KpeXyWTo1KkTvL29VZ7PysrCkCFD4OrqCjc3N6xbtw6PHj1SDg/p5OSkMmTknDlzMG3aNJU78OLi4tCjRw+4urqiQ4cOePfddyGVSpWvT506Fa6urvDw8ICfnx/S0tJq7X3QRNdhOisb9nb9+vXo2rUrPDw80LVrV6xdu1Zt/U8++QSmpqbIzMxUeV4sFis/k6qG1q3qvZsyZYpymNQePXogMTGx2u+LLnR97yobqjYqKkpliFELCwuV9+/KlSsYOHAgOnTogC5duqBnz544ePBglbUVFRVh/vz5cHJyQseOHdGlSxfs2LGjRu2tMwxwfTmpZWfPnmUPHjxgDg4O7Pfff69y+Q8++IBNnTpVZSQxuVzOPD09WUxMjPLx48ePVdaravSx+Ph41qZNG+WNJsXFxey9995jvXv3Vt4ccOjQIeUoZ0eOHGHOzs46tdHBwUGn5crTZ5jOyoa9/fvvv5XLSaVS1q5dO/bbb78pnystLWXt2rVjffv2ZUuXLlWrW/GZVLYPXd67v/76S7nda9euMWtra51uujD0e6frULWPHz9mFhYWyt+r9PR0Zm1tzQ4fPqxcJjs7m+3cubPK+saOHcuGDx/OXr58yRhjLCMjg7m5ubHNmzfr08w6iUJbQHQJ7XPnzrEhQ4aoBfDJkydZ7969K11XU2j7+/uzI0eOMMYY6927N1u/fr3K64WFhaxt27bs1KlTatt79uwZe+2111hpaWml+2WsesFT3WE6Kxv29vHjx8zW1lZ5JypjZYHr7e3Nrly5whwcHFTao+0zqbgPfd+7M2fOMBsbG4OFdnXeu6qGql25ciV7++23lY8nTJjA5s6dq3X5VatWMW9vb+bh4cF69OjBLl68yBhj7M6dO6xhw4bs+fPnKssfPXqU2dnZVVpjfUDdI3VIfn4+Zs+ejY0bN6q9duPGDbRo0QJjxoxB9+7dMWzYMNy/f1+v7V+9ehU9e/ZUee61116Dl5cXrl69qrb8mjVrEBISovds07qq7jCdmoa9/emnn9C5c2c4ODhg7ty56Nq1q/K17du3Y+rUqfD09ETz5s116raouA9d37sFCxbA0dERw4cPx/79+9VGjasthhjiNDo6GhEREcrHV65cUWtzeeHh4bh06RKuXbuGtWvXKte9evUqnJ2dYW1trbJ8z549kZ2djWfPntWoTqGjUf7qkLlz52LmzJlo27Yt7ty5o/JacXExTp06hYsXL6Jz587YsmULxowZg19//VWvfWgKEaZhJIQ9e/YgJiZG64h0paWl8PLyUj5W9K0DQKtWrZCQkKB3PZrqqEjbsLdhYWEICwuDRCLBsGHDEBISAldXVzx//hwnT57E1q1bAQARERHYvn07+vfvr/c+dHnvvvzyS3z55Zc4deoU5s6di9TUVLz22msqy3D13lUmNTUVubm5WkcT1OTatWtYvnw5cnJyYGZmhhs3bqCoqEitNqKKQrsOSUlJwbFjx/Dpp59CJpPhr7/+QufOnfHHH3/AwcEB3bt3R+fOnQEAEyZMwIwZM1BaWgpTU1Odtu/p6Ynz588rAwIoO2F09epVvP/++8rn9u3bh2XLliExMREtW7bUuC1TU1OVk5RisbjKk5ZhYWG4e/cugLJg1HeYTl2GvRWLxXjzzTcRFxcHV1dX7N69GyUlJco2l5aWIicnBzk5OWpHgpXtQ9f3TqFfv35477338Pvvv6sENMDNe1eV7du3Y9KkSSq/S15eXrhw4QKGDRumtnxRURFGjBiBpKQkeHl5ITc3F5aWligqKkL37t1x+/Zttff4woULsLOzUw45W29x2ztD9KHriUjG1PunX758yTp06KAche7nn39m3bp1q3QdxlT7tI8ePcpsbW3ZtWvXGGP/nkzr2bOnsu913759zMnJiUkkEr3bpq979+6pnUzbuHGjxmXPnj3L2rVrpzwRWN6NGzeU/3769ClzcnJiJ06cYIwx1qVLFxYfH6+y/IgRI9iaNWuUdSs+k8r2UdV7V1xczG7fvq1c/pdffmHNmzfXOrdjeYZ+7xS09Wnn5eWxpk2bqtTPGGPXr19n1tbWLC4uTvncgwcP2ObNm5lUKmXm5ubsyZMnjDHGPv/8c5Wpz0aNGsWGDx/O8vPzGWP/noisqsb6gEJbACIjI1nbtm2Zqakpa9WqFXN0dFS+VnHYSgVNAXz8+HHm7u7OunXrxvz8/JTDRla2TvnQZoyxgwcPMk9PT+bi4sLEYjGLiIhQuerBzMyM2dnZqQwvW/GEkibVCR7GtA/T+fDhQ+bu7q5crrJhb2fMmME6deqkfG8UJwwvXrzIWrZsqbwapvx7oPjCKx/ale1DsZ62904mk7FevXqxzp07s27durGePXuyxMREnd4DQ793lQ1Vyxhj27ZtY35+fhr38euvv7L+/fuz9u3bsy5durBevXopf19XrlzJHBwcWJ8+fdjXX3+tEtoymYx99NFHrEOHDszNzY117tyZbd26tVrtrGtoaFZCCBEQunqEEEIEhEKbEEIEhEKbEEIEhEKbEEIEhEKbEEIEhEKbEEIEhEKbEEIEhEKbEEIEhEKbEEIEhEKbEEIEhEKbEEIE5P8DnHpWR6udNW4AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 354.331x334.646 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = wholerock.drop(columns=[\"FeO\"])\n",
    "FeOi_predict = mpc.FeOi_prediction(x=x, FeO=wholerock[\"FeO\"])\n",
    "\n",
    "do_not_use = [\"MnO\", \"P2O5\", \"Cr2O3\", \"total\"]\n",
    "\n",
    "model_fits = FeOi_predict.calculate_model_fits(exclude=do_not_use)\n",
    "FeOi_predict.select_predictors(idx=3)\n",
    "FeO_model = FeOi_predict.model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9dde75ee",
   "metadata": {},
   "source": [
    "Here are all the Kd models we're going to test. You can also do this for Fe3Fe2 or melt thermometers instead with `mpc.Fe3Fe2_models` or `mpc.melt_thermometers`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "67874ae8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['blundy2020',\n",
       " 'fixed',\n",
       " 'putirka2016_8a',\n",
       " 'putirka2016_8b',\n",
       " 'putirka2016_8c',\n",
       " 'putirka2016_8d',\n",
       " 'saper2022',\n",
       " 'sun2020',\n",
       " 'toplis2005']"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Kd_models = mpc.Kd_ol_FeMg_models\n",
    "Kd_models"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "974a0d3e",
   "metadata": {},
   "source": [
    "Now we calculate PEC corrections iteratively for each Kd model. \n",
    "\n",
    "For a fixed Kd we also need to define the value and its error and pass them as a list to the configuration\n",
    "\n",
    "Crystallization extents are stored in the `pec_results` dataframe and corrected melt compositions in the `corrected_inclusions` dictionary. This dictionary uses Kd model names as keys and stores compositions in dataframes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b51c84f6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "model: toplis2005\n",
      "009/009\n",
      "Equilibrating ... |██████████████████████████████| 100% [10/10] in 3.6s \n",
      "Correcting    ... |██████████████████████████████| 100% [10/10] in 13.0s \n"
     ]
    }
   ],
   "source": [
    "pec_results = pd.DataFrame(index=melt.index, dtype=float)\n",
    "corrected_inclusions = {}\n",
    "\n",
    "Kd_value = 0.35\n",
    "Kd_error = 0.02\n",
    "\n",
    "for i, model in enumerate(Kd_models):\n",
    "\n",
    "    # add Kd value and error if the model is 'fixed'\n",
    "    set_model = model if model != \"fixed\" else (model, Kd_value, Kd_error)\n",
    "    # set the Kd model\n",
    "    mpc.model_configuration.Kd_model = set_model\n",
    "\n",
    "    clear_output()\n",
    "    print(f\"model: {model}\\n{i+1:03}/{len(Kd_models):03}\")\n",
    "\n",
    "    pec_model = mpc.PEC(inclusions=melt, olivines=olivine, P_bar=pressure, FeO_target=FeO_model)\n",
    "    melts_corrected, pec, checks = pec_model.correct()\n",
    "\n",
    "    pec_results[set_model] = pec[\"total_crystallisation\"]\n",
    "    corrected_inclusions[set_model] = melts_corrected.copy()\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd67ffd4",
   "metadata": {},
   "source": [
    "These are the crystallization extents calculated with the different models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "9141cded",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>blundy2020</th>\n",
       "      <th>(fixed, 0.35, 0.02)</th>\n",
       "      <th>putirka2016_8a</th>\n",
       "      <th>putirka2016_8b</th>\n",
       "      <th>putirka2016_8c</th>\n",
       "      <th>putirka2016_8d</th>\n",
       "      <th>saper2022</th>\n",
       "      <th>sun2020</th>\n",
       "      <th>toplis2005</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>name</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>PI032-04-01</th>\n",
       "      <td>9.669983</td>\n",
       "      <td>12.221924</td>\n",
       "      <td>10.917725</td>\n",
       "      <td>10.917725</td>\n",
       "      <td>10.910095</td>\n",
       "      <td>10.841650</td>\n",
       "      <td>6.590302</td>\n",
       "      <td>12.917725</td>\n",
       "      <td>9.796851</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PI032-04-02</th>\n",
       "      <td>10.807312</td>\n",
       "      <td>13.931787</td>\n",
       "      <td>12.135925</td>\n",
       "      <td>12.056689</td>\n",
       "      <td>12.250586</td>\n",
       "      <td>12.050586</td>\n",
       "      <td>7.784198</td>\n",
       "      <td>14.326660</td>\n",
       "      <td>10.807312</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PI041-02-02</th>\n",
       "      <td>1.159460</td>\n",
       "      <td>2.836523</td>\n",
       "      <td>1.929163</td>\n",
       "      <td>1.754883</td>\n",
       "      <td>1.927637</td>\n",
       "      <td>1.756409</td>\n",
       "      <td>-1.006763</td>\n",
       "      <td>2.918481</td>\n",
       "      <td>-0.005502</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PI041-03-01</th>\n",
       "      <td>14.807861</td>\n",
       "      <td>17.614844</td>\n",
       "      <td>16.084595</td>\n",
       "      <td>15.304321</td>\n",
       "      <td>15.492114</td>\n",
       "      <td>14.007861</td>\n",
       "      <td>11.345508</td>\n",
       "      <td>19.764160</td>\n",
       "      <td>12.645239</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PI041-03-03</th>\n",
       "      <td>14.341699</td>\n",
       "      <td>17.221924</td>\n",
       "      <td>15.762085</td>\n",
       "      <td>14.461047</td>\n",
       "      <td>15.250366</td>\n",
       "      <td>13.350854</td>\n",
       "      <td>11.000000</td>\n",
       "      <td>19.539307</td>\n",
       "      <td>12.131775</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PI041-05-04</th>\n",
       "      <td>-2.556934</td>\n",
       "      <td>-0.966846</td>\n",
       "      <td>-1.805542</td>\n",
       "      <td>-2.505273</td>\n",
       "      <td>-2.147778</td>\n",
       "      <td>-2.496118</td>\n",
       "      <td>-4.984302</td>\n",
       "      <td>-1.109082</td>\n",
       "      <td>-4.403699</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PI041-05-06</th>\n",
       "      <td>3.488623</td>\n",
       "      <td>5.703125</td>\n",
       "      <td>4.443066</td>\n",
       "      <td>3.703882</td>\n",
       "      <td>4.091675</td>\n",
       "      <td>3.344861</td>\n",
       "      <td>0.800000</td>\n",
       "      <td>5.579199</td>\n",
       "      <td>1.994775</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PI041-07-01</th>\n",
       "      <td>13.533203</td>\n",
       "      <td>16.312500</td>\n",
       "      <td>14.638818</td>\n",
       "      <td>13.720996</td>\n",
       "      <td>14.251025</td>\n",
       "      <td>12.945410</td>\n",
       "      <td>9.790521</td>\n",
       "      <td>17.044922</td>\n",
       "      <td>12.212329</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PI041-07-02</th>\n",
       "      <td>12.960449</td>\n",
       "      <td>15.674243</td>\n",
       "      <td>14.072168</td>\n",
       "      <td>13.151294</td>\n",
       "      <td>13.681323</td>\n",
       "      <td>12.419739</td>\n",
       "      <td>9.193353</td>\n",
       "      <td>17.053857</td>\n",
       "      <td>11.118481</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PI052-01-02</th>\n",
       "      <td>-7.058838</td>\n",
       "      <td>-5.447388</td>\n",
       "      <td>-6.301343</td>\n",
       "      <td>-6.634424</td>\n",
       "      <td>-6.479980</td>\n",
       "      <td>-6.276929</td>\n",
       "      <td>-9.327405</td>\n",
       "      <td>-5.450439</td>\n",
       "      <td>-8.164673</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             blundy2020  (fixed, 0.35, 0.02)  putirka2016_8a  putirka2016_8b  \\\n",
       "name                                                                           \n",
       "PI032-04-01    9.669983            12.221924       10.917725       10.917725   \n",
       "PI032-04-02   10.807312            13.931787       12.135925       12.056689   \n",
       "PI041-02-02    1.159460             2.836523        1.929163        1.754883   \n",
       "PI041-03-01   14.807861            17.614844       16.084595       15.304321   \n",
       "PI041-03-03   14.341699            17.221924       15.762085       14.461047   \n",
       "PI041-05-04   -2.556934            -0.966846       -1.805542       -2.505273   \n",
       "PI041-05-06    3.488623             5.703125        4.443066        3.703882   \n",
       "PI041-07-01   13.533203            16.312500       14.638818       13.720996   \n",
       "PI041-07-02   12.960449            15.674243       14.072168       13.151294   \n",
       "PI052-01-02   -7.058838            -5.447388       -6.301343       -6.634424   \n",
       "\n",
       "             putirka2016_8c  putirka2016_8d  saper2022    sun2020  toplis2005  \n",
       "name                                                                           \n",
       "PI032-04-01       10.910095       10.841650   6.590302  12.917725    9.796851  \n",
       "PI032-04-02       12.250586       12.050586   7.784198  14.326660   10.807312  \n",
       "PI041-02-02        1.927637        1.756409  -1.006763   2.918481   -0.005502  \n",
       "PI041-03-01       15.492114       14.007861  11.345508  19.764160   12.645239  \n",
       "PI041-03-03       15.250366       13.350854  11.000000  19.539307   12.131775  \n",
       "PI041-05-04       -2.147778       -2.496118  -4.984302  -1.109082   -4.403699  \n",
       "PI041-05-06        4.091675        3.344861   0.800000   5.579199    1.994775  \n",
       "PI041-07-01       14.251025       12.945410   9.790521  17.044922   12.212329  \n",
       "PI041-07-02       13.681323       12.419739   9.193353  17.053857   11.118481  \n",
       "PI052-01-02       -6.479980       -6.276929  -9.327405  -5.450439   -8.164673  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pec_results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6eacab79",
   "metadata": {},
   "source": [
    "And here are two examples of the associated corrected melt compositions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "626e439b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>SiO2</th>\n",
       "      <th>Al2O3</th>\n",
       "      <th>MgO</th>\n",
       "      <th>CaO</th>\n",
       "      <th>FeO</th>\n",
       "      <th>Na2O</th>\n",
       "      <th>K2O</th>\n",
       "      <th>MnO</th>\n",
       "      <th>TiO2</th>\n",
       "      <th>P2O5</th>\n",
       "      <th>Cr2O3</th>\n",
       "      <th>CO2</th>\n",
       "      <th>H2O</th>\n",
       "      <th>F</th>\n",
       "      <th>S</th>\n",
       "      <th>Cl</th>\n",
       "      <th>total</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>name</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>PI032-04-01</th>\n",
       "      <td>48.943719</td>\n",
       "      <td>14.013614</td>\n",
       "      <td>7.630094</td>\n",
       "      <td>9.705938</td>\n",
       "      <td>10.314919</td>\n",
       "      <td>3.602243</td>\n",
       "      <td>0.672540</td>\n",
       "      <td>0.138315</td>\n",
       "      <td>2.456433</td>\n",
       "      <td>0.279344</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.605778</td>\n",
       "      <td>1.388135</td>\n",
       "      <td>0.076006</td>\n",
       "      <td>0.140193</td>\n",
       "      <td>0.032730</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PI032-04-02</th>\n",
       "      <td>48.479654</td>\n",
       "      <td>14.629497</td>\n",
       "      <td>7.547376</td>\n",
       "      <td>9.386726</td>\n",
       "      <td>10.340466</td>\n",
       "      <td>3.449630</td>\n",
       "      <td>0.910037</td>\n",
       "      <td>0.140879</td>\n",
       "      <td>2.567158</td>\n",
       "      <td>0.325500</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.650701</td>\n",
       "      <td>1.288021</td>\n",
       "      <td>0.081179</td>\n",
       "      <td>0.160070</td>\n",
       "      <td>0.043106</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PI041-02-02</th>\n",
       "      <td>49.035956</td>\n",
       "      <td>16.830785</td>\n",
       "      <td>5.174158</td>\n",
       "      <td>9.123172</td>\n",
       "      <td>10.257030</td>\n",
       "      <td>3.761868</td>\n",
       "      <td>1.070182</td>\n",
       "      <td>0.157860</td>\n",
       "      <td>2.784544</td>\n",
       "      <td>0.557070</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.460645</td>\n",
       "      <td>0.650866</td>\n",
       "      <td>0.047317</td>\n",
       "      <td>0.067540</td>\n",
       "      <td>0.021008</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PI041-03-01</th>\n",
       "      <td>45.805062</td>\n",
       "      <td>15.410386</td>\n",
       "      <td>7.866142</td>\n",
       "      <td>10.765738</td>\n",
       "      <td>10.734843</td>\n",
       "      <td>3.273348</td>\n",
       "      <td>1.119837</td>\n",
       "      <td>0.117426</td>\n",
       "      <td>3.084403</td>\n",
       "      <td>0.539931</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.780285</td>\n",
       "      <td>0.286932</td>\n",
       "      <td>0.077762</td>\n",
       "      <td>0.084965</td>\n",
       "      <td>0.052940</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PI041-03-03</th>\n",
       "      <td>45.225725</td>\n",
       "      <td>15.624567</td>\n",
       "      <td>7.725401</td>\n",
       "      <td>10.983367</td>\n",
       "      <td>10.806821</td>\n",
       "      <td>3.356069</td>\n",
       "      <td>1.146547</td>\n",
       "      <td>0.094428</td>\n",
       "      <td>3.203774</td>\n",
       "      <td>0.508575</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.806345</td>\n",
       "      <td>0.305534</td>\n",
       "      <td>0.080243</td>\n",
       "      <td>0.078065</td>\n",
       "      <td>0.054541</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PI041-05-04</th>\n",
       "      <td>47.739154</td>\n",
       "      <td>18.499127</td>\n",
       "      <td>4.168579</td>\n",
       "      <td>9.356278</td>\n",
       "      <td>9.387706</td>\n",
       "      <td>4.577251</td>\n",
       "      <td>1.598948</td>\n",
       "      <td>0.142979</td>\n",
       "      <td>2.474740</td>\n",
       "      <td>0.821254</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.507617</td>\n",
       "      <td>0.459434</td>\n",
       "      <td>0.087337</td>\n",
       "      <td>0.120858</td>\n",
       "      <td>0.058737</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PI041-05-06</th>\n",
       "      <td>46.350575</td>\n",
       "      <td>16.864911</td>\n",
       "      <td>5.158561</td>\n",
       "      <td>8.871758</td>\n",
       "      <td>11.454093</td>\n",
       "      <td>3.976689</td>\n",
       "      <td>1.407932</td>\n",
       "      <td>0.173516</td>\n",
       "      <td>3.611212</td>\n",
       "      <td>0.615599</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.636053</td>\n",
       "      <td>0.582065</td>\n",
       "      <td>0.106742</td>\n",
       "      <td>0.127315</td>\n",
       "      <td>0.062978</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PI041-07-01</th>\n",
       "      <td>45.738018</td>\n",
       "      <td>15.077100</td>\n",
       "      <td>7.642742</td>\n",
       "      <td>9.589114</td>\n",
       "      <td>11.633964</td>\n",
       "      <td>3.098567</td>\n",
       "      <td>1.259341</td>\n",
       "      <td>0.150583</td>\n",
       "      <td>3.496301</td>\n",
       "      <td>0.549176</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.462927</td>\n",
       "      <td>1.010345</td>\n",
       "      <td>0.074330</td>\n",
       "      <td>0.161746</td>\n",
       "      <td>0.055747</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PI041-07-02</th>\n",
       "      <td>45.806314</td>\n",
       "      <td>15.375532</td>\n",
       "      <td>7.401513</td>\n",
       "      <td>9.952743</td>\n",
       "      <td>11.490914</td>\n",
       "      <td>3.149061</td>\n",
       "      <td>1.342889</td>\n",
       "      <td>0.139523</td>\n",
       "      <td>3.489543</td>\n",
       "      <td>0.608231</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.367453</td>\n",
       "      <td>0.608398</td>\n",
       "      <td>0.072640</td>\n",
       "      <td>0.138611</td>\n",
       "      <td>0.056636</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PI052-01-02</th>\n",
       "      <td>49.177883</td>\n",
       "      <td>16.933432</td>\n",
       "      <td>4.257739</td>\n",
       "      <td>10.328566</td>\n",
       "      <td>8.486900</td>\n",
       "      <td>4.901528</td>\n",
       "      <td>1.524944</td>\n",
       "      <td>0.240910</td>\n",
       "      <td>1.754921</td>\n",
       "      <td>0.681335</td>\n",
       "      <td>0.047285</td>\n",
       "      <td>0.286681</td>\n",
       "      <td>1.116948</td>\n",
       "      <td>0.081782</td>\n",
       "      <td>0.125735</td>\n",
       "      <td>0.053411</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  SiO2      Al2O3       MgO        CaO        FeO      Na2O  \\\n",
       "name                                                                          \n",
       "PI032-04-01  48.943719  14.013614  7.630094   9.705938  10.314919  3.602243   \n",
       "PI032-04-02  48.479654  14.629497  7.547376   9.386726  10.340466  3.449630   \n",
       "PI041-02-02  49.035956  16.830785  5.174158   9.123172  10.257030  3.761868   \n",
       "PI041-03-01  45.805062  15.410386  7.866142  10.765738  10.734843  3.273348   \n",
       "PI041-03-03  45.225725  15.624567  7.725401  10.983367  10.806821  3.356069   \n",
       "PI041-05-04  47.739154  18.499127  4.168579   9.356278   9.387706  4.577251   \n",
       "PI041-05-06  46.350575  16.864911  5.158561   8.871758  11.454093  3.976689   \n",
       "PI041-07-01  45.738018  15.077100  7.642742   9.589114  11.633964  3.098567   \n",
       "PI041-07-02  45.806314  15.375532  7.401513   9.952743  11.490914  3.149061   \n",
       "PI052-01-02  49.177883  16.933432  4.257739  10.328566   8.486900  4.901528   \n",
       "\n",
       "                  K2O       MnO      TiO2      P2O5     Cr2O3       CO2  \\\n",
       "name                                                                      \n",
       "PI032-04-01  0.672540  0.138315  2.456433  0.279344  0.000000  0.605778   \n",
       "PI032-04-02  0.910037  0.140879  2.567158  0.325500  0.000000  0.650701   \n",
       "PI041-02-02  1.070182  0.157860  2.784544  0.557070  0.000000  0.460645   \n",
       "PI041-03-01  1.119837  0.117426  3.084403  0.539931  0.000000  0.780285   \n",
       "PI041-03-03  1.146547  0.094428  3.203774  0.508575  0.000000  0.806345   \n",
       "PI041-05-04  1.598948  0.142979  2.474740  0.821254  0.000000  0.507617   \n",
       "PI041-05-06  1.407932  0.173516  3.611212  0.615599  0.000000  0.636053   \n",
       "PI041-07-01  1.259341  0.150583  3.496301  0.549176  0.000000  0.462927   \n",
       "PI041-07-02  1.342889  0.139523  3.489543  0.608231  0.000000  0.367453   \n",
       "PI052-01-02  1.524944  0.240910  1.754921  0.681335  0.047285  0.286681   \n",
       "\n",
       "                  H2O         F         S        Cl  total  \n",
       "name                                                        \n",
       "PI032-04-01  1.388135  0.076006  0.140193  0.032730  100.0  \n",
       "PI032-04-02  1.288021  0.081179  0.160070  0.043106  100.0  \n",
       "PI041-02-02  0.650866  0.047317  0.067540  0.021008  100.0  \n",
       "PI041-03-01  0.286932  0.077762  0.084965  0.052940  100.0  \n",
       "PI041-03-03  0.305534  0.080243  0.078065  0.054541  100.0  \n",
       "PI041-05-04  0.459434  0.087337  0.120858  0.058737  100.0  \n",
       "PI041-05-06  0.582065  0.106742  0.127315  0.062978  100.0  \n",
       "PI041-07-01  1.010345  0.074330  0.161746  0.055747  100.0  \n",
       "PI041-07-02  0.608398  0.072640  0.138611  0.056636  100.0  \n",
       "PI052-01-02  1.116948  0.081782  0.125735  0.053411  100.0  "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "corrected_inclusions[\"blundy2020\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "17769b67",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>SiO2</th>\n",
       "      <th>Al2O3</th>\n",
       "      <th>MgO</th>\n",
       "      <th>CaO</th>\n",
       "      <th>FeO</th>\n",
       "      <th>Na2O</th>\n",
       "      <th>K2O</th>\n",
       "      <th>MnO</th>\n",
       "      <th>TiO2</th>\n",
       "      <th>P2O5</th>\n",
       "      <th>Cr2O3</th>\n",
       "      <th>CO2</th>\n",
       "      <th>H2O</th>\n",
       "      <th>F</th>\n",
       "      <th>S</th>\n",
       "      <th>Cl</th>\n",
       "      <th>total</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>name</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>PI032-04-01</th>\n",
       "      <td>49.206233</td>\n",
       "      <td>14.365838</td>\n",
       "      <td>6.636341</td>\n",
       "      <td>9.940199</td>\n",
       "      <td>10.233057</td>\n",
       "      <td>3.692838</td>\n",
       "      <td>0.689454</td>\n",
       "      <td>0.132211</td>\n",
       "      <td>2.518211</td>\n",
       "      <td>0.286370</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.621013</td>\n",
       "      <td>1.423046</td>\n",
       "      <td>0.077918</td>\n",
       "      <td>0.143719</td>\n",
       "      <td>0.033553</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PI032-04-02</th>\n",
       "      <td>48.691045</td>\n",
       "      <td>14.986061</td>\n",
       "      <td>6.571461</td>\n",
       "      <td>9.606162</td>\n",
       "      <td>10.303377</td>\n",
       "      <td>3.533640</td>\n",
       "      <td>0.932200</td>\n",
       "      <td>0.135733</td>\n",
       "      <td>2.629678</td>\n",
       "      <td>0.333427</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.666547</td>\n",
       "      <td>1.319389</td>\n",
       "      <td>0.083156</td>\n",
       "      <td>0.163968</td>\n",
       "      <td>0.044156</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PI041-02-02</th>\n",
       "      <td>49.228087</td>\n",
       "      <td>17.112982</td>\n",
       "      <td>4.413466</td>\n",
       "      <td>9.268368</td>\n",
       "      <td>10.248313</td>\n",
       "      <td>3.825082</td>\n",
       "      <td>1.088165</td>\n",
       "      <td>0.149435</td>\n",
       "      <td>2.831335</td>\n",
       "      <td>0.566431</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.468385</td>\n",
       "      <td>0.661803</td>\n",
       "      <td>0.048112</td>\n",
       "      <td>0.068675</td>\n",
       "      <td>0.021361</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PI041-03-01</th>\n",
       "      <td>45.991137</td>\n",
       "      <td>15.828645</td>\n",
       "      <td>6.792065</td>\n",
       "      <td>11.047336</td>\n",
       "      <td>10.675971</td>\n",
       "      <td>3.362487</td>\n",
       "      <td>1.150332</td>\n",
       "      <td>0.111175</td>\n",
       "      <td>3.168397</td>\n",
       "      <td>0.554635</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.801534</td>\n",
       "      <td>0.294746</td>\n",
       "      <td>0.079879</td>\n",
       "      <td>0.087279</td>\n",
       "      <td>0.054382</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PI041-03-03</th>\n",
       "      <td>45.388617</td>\n",
       "      <td>16.031367</td>\n",
       "      <td>6.639883</td>\n",
       "      <td>11.258630</td>\n",
       "      <td>10.804800</td>\n",
       "      <td>3.443792</td>\n",
       "      <td>1.176516</td>\n",
       "      <td>0.087660</td>\n",
       "      <td>3.287515</td>\n",
       "      <td>0.521868</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.827421</td>\n",
       "      <td>0.313520</td>\n",
       "      <td>0.082340</td>\n",
       "      <td>0.080105</td>\n",
       "      <td>0.055966</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PI041-05-04</th>\n",
       "      <td>47.880455</td>\n",
       "      <td>18.844295</td>\n",
       "      <td>3.375884</td>\n",
       "      <td>9.522744</td>\n",
       "      <td>9.336220</td>\n",
       "      <td>4.662865</td>\n",
       "      <td>1.628855</td>\n",
       "      <td>0.133975</td>\n",
       "      <td>2.521028</td>\n",
       "      <td>0.836615</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.517112</td>\n",
       "      <td>0.468028</td>\n",
       "      <td>0.088971</td>\n",
       "      <td>0.123118</td>\n",
       "      <td>0.059836</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PI041-05-06</th>\n",
       "      <td>46.481016</td>\n",
       "      <td>17.212464</td>\n",
       "      <td>4.258297</td>\n",
       "      <td>9.046185</td>\n",
       "      <td>11.480980</td>\n",
       "      <td>4.059090</td>\n",
       "      <td>1.437106</td>\n",
       "      <td>0.163918</td>\n",
       "      <td>3.686041</td>\n",
       "      <td>0.628355</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.649232</td>\n",
       "      <td>0.594126</td>\n",
       "      <td>0.108954</td>\n",
       "      <td>0.129953</td>\n",
       "      <td>0.064283</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PI041-07-01</th>\n",
       "      <td>45.867609</td>\n",
       "      <td>15.532570</td>\n",
       "      <td>6.515228</td>\n",
       "      <td>9.866277</td>\n",
       "      <td>11.598683</td>\n",
       "      <td>3.192366</td>\n",
       "      <td>1.297463</td>\n",
       "      <td>0.143335</td>\n",
       "      <td>3.602140</td>\n",
       "      <td>0.565800</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.476941</td>\n",
       "      <td>1.040930</td>\n",
       "      <td>0.076580</td>\n",
       "      <td>0.166642</td>\n",
       "      <td>0.057435</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PI041-07-02</th>\n",
       "      <td>45.927100</td>\n",
       "      <td>15.826762</td>\n",
       "      <td>6.307860</td>\n",
       "      <td>10.232779</td>\n",
       "      <td>11.449861</td>\n",
       "      <td>3.241603</td>\n",
       "      <td>1.382352</td>\n",
       "      <td>0.133202</td>\n",
       "      <td>3.592089</td>\n",
       "      <td>0.626105</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.378252</td>\n",
       "      <td>0.626277</td>\n",
       "      <td>0.074774</td>\n",
       "      <td>0.142684</td>\n",
       "      <td>0.058300</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PI052-01-02</th>\n",
       "      <td>49.345297</td>\n",
       "      <td>17.251422</td>\n",
       "      <td>3.447559</td>\n",
       "      <td>10.522525</td>\n",
       "      <td>8.426774</td>\n",
       "      <td>4.993573</td>\n",
       "      <td>1.553581</td>\n",
       "      <td>0.233274</td>\n",
       "      <td>1.787877</td>\n",
       "      <td>0.694130</td>\n",
       "      <td>0.048173</td>\n",
       "      <td>0.292065</td>\n",
       "      <td>1.137923</td>\n",
       "      <td>0.083317</td>\n",
       "      <td>0.128096</td>\n",
       "      <td>0.054414</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  SiO2      Al2O3       MgO        CaO        FeO      Na2O  \\\n",
       "name                                                                          \n",
       "PI032-04-01  49.206233  14.365838  6.636341   9.940199  10.233057  3.692838   \n",
       "PI032-04-02  48.691045  14.986061  6.571461   9.606162  10.303377  3.533640   \n",
       "PI041-02-02  49.228087  17.112982  4.413466   9.268368  10.248313  3.825082   \n",
       "PI041-03-01  45.991137  15.828645  6.792065  11.047336  10.675971  3.362487   \n",
       "PI041-03-03  45.388617  16.031367  6.639883  11.258630  10.804800  3.443792   \n",
       "PI041-05-04  47.880455  18.844295  3.375884   9.522744   9.336220  4.662865   \n",
       "PI041-05-06  46.481016  17.212464  4.258297   9.046185  11.480980  4.059090   \n",
       "PI041-07-01  45.867609  15.532570  6.515228   9.866277  11.598683  3.192366   \n",
       "PI041-07-02  45.927100  15.826762  6.307860  10.232779  11.449861  3.241603   \n",
       "PI052-01-02  49.345297  17.251422  3.447559  10.522525   8.426774  4.993573   \n",
       "\n",
       "                  K2O       MnO      TiO2      P2O5     Cr2O3       CO2  \\\n",
       "name                                                                      \n",
       "PI032-04-01  0.689454  0.132211  2.518211  0.286370  0.000000  0.621013   \n",
       "PI032-04-02  0.932200  0.135733  2.629678  0.333427  0.000000  0.666547   \n",
       "PI041-02-02  1.088165  0.149435  2.831335  0.566431  0.000000  0.468385   \n",
       "PI041-03-01  1.150332  0.111175  3.168397  0.554635  0.000000  0.801534   \n",
       "PI041-03-03  1.176516  0.087660  3.287515  0.521868  0.000000  0.827421   \n",
       "PI041-05-04  1.628855  0.133975  2.521028  0.836615  0.000000  0.517112   \n",
       "PI041-05-06  1.437106  0.163918  3.686041  0.628355  0.000000  0.649232   \n",
       "PI041-07-01  1.297463  0.143335  3.602140  0.565800  0.000000  0.476941   \n",
       "PI041-07-02  1.382352  0.133202  3.592089  0.626105  0.000000  0.378252   \n",
       "PI052-01-02  1.553581  0.233274  1.787877  0.694130  0.048173  0.292065   \n",
       "\n",
       "                  H2O         F         S        Cl  total  \n",
       "name                                                        \n",
       "PI032-04-01  1.423046  0.077918  0.143719  0.033553  100.0  \n",
       "PI032-04-02  1.319389  0.083156  0.163968  0.044156  100.0  \n",
       "PI041-02-02  0.661803  0.048112  0.068675  0.021361  100.0  \n",
       "PI041-03-01  0.294746  0.079879  0.087279  0.054382  100.0  \n",
       "PI041-03-03  0.313520  0.082340  0.080105  0.055966  100.0  \n",
       "PI041-05-04  0.468028  0.088971  0.123118  0.059836  100.0  \n",
       "PI041-05-06  0.594126  0.108954  0.129953  0.064283  100.0  \n",
       "PI041-07-01  1.040930  0.076580  0.166642  0.057435  100.0  \n",
       "PI041-07-02  0.626277  0.074774  0.142684  0.058300  100.0  \n",
       "PI052-01-02  1.137923  0.083317  0.128096  0.054414  100.0  "
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "corrected_inclusions[\"saper2022\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79fa84e5",
   "metadata": {},
   "source": [
    "Here we plot PEC results against uncorrected melt FeO contents. It shows that using the Sun et al. (2020) Kd model results in the largest PEC corrections, and the Saper et al. (2022) model the smallest. PEC differences between models for single inclusions are up to ca. 9%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "76cbfc8a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAAFaCAYAAADVfgw3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABblUlEQVR4nO3deXxM1/vA8c/MZCcbktiC2vfYiQgRqqq0FC3VBS3V6peuainVH0rpRqtVtEHRVltLLbW2sS+xK0ooIiUSZF8kM3N/f0xnZCRhkkwymeR5v155MXfunHvmGM+cnHvOc1SKoigIIYSwC2pbV0AIIYTlJGgLIYQdkaAthBB2RIK2EELYEQnaQghhRyRoCyGEHZGgLYQQdkSCthBC2BEHW1egKOj1eq5du4a7uzsqlcrW1RFCiAdSFIXk5GSqVq2KWp13f7pUBu1r167h7+9v62oIIUS+Xb16lerVq+f5fKkM2u7u7oDhzXt4eFi9fEVRSE9Px9XVtUz35KUd7pK2MJB2MChIOyQlJeHv72+KX3kplUHb2EgeHh5FFrQdHR3lgyntYCJtYSDtYFCYdnjQ+XIjUggh7IgEbSGEsCMStIUQwo5I0BZCCDtS4oJ2RkYGffv2pX79+rRo0YKePXty+fJlAGJjY+nZsyf16tWjadOm7Nmzx7aVFUKIYlbigjbAyJEjOXfuHMePH6d3796MHDkSgPHjx9OhQwciIyMJCwtjyJAhaLVaG9dWCCGKT4kL2i4uLvTq1cs07aVDhw78888/AKxatYrRo0cD0LZtW/z8/KS3LYQoU0r8PO158+bRp08fbt26hV6vx8fHx/RcrVq1iIqKyvO1iYmJZN8C09nZGWdn50LXyVhmWd9eU9rhLmkLA2kHg4K0g6Xnluig/eGHHxIZGcmCBQtIT0/PMen8QW+yRo0aZo8nTpzIpEmT8l2P5IxEfj05mzv6OJzVPgxo/i7lXTzIyMjId1mlkbTDXdIWBtIOBvlph/T0dIvOK7FB++OPP2b16tVs374dNzc33NzcAIiLizP1tq9cuZIjMGcXFRVltiKyoD3tZYcn0zBgKRW8k7kd787qU2qebzMNFxeXMr/qKyMjo8y3A0hbGEk7GBSkHbKysiw6r0QG7U8//ZQffviB7du34+XlZTo+cOBA5s+fz9SpU4mIiCAmJoZOnTrlWY6np6fFy9iT0hNYGTGdDF0sLhpfhrSbjLuLJwAZujgqeCcDUME7mQxdHGBYblqWP5hG0g53SVsYSDsY5KcdLD2vxAXt6Oho3nrrLWrXrk3Xrl0BQw/54MGDfPTRRzz33HPUq1cPJycnvv/+exwcrPMWVkZMp36z70y96RWHYFTnjwFw0fhyO97d9JyLxucBpQkhRNEocUG7evXqeY5V+/n5sXXr1iK5boYu9p7edKzpuSHtJrPiEKZe+DNt34OyfZ9FCGEjJS5o20rO3rSv6Tl3F09Trxvupl0UQojiJkH7P/f2poe0m2zrKgkhRA4StP9zb29aCCFKohK3IlIIIUTeJGgLIYQdkaAthBB2RIK2EELYEQnaQghhRyRoCyGEHZGgLYQQdkSCthBC2BEJ2kIIYUckaAshhB2RoC2EEHZEgraFFG0m2p+HoWgzbV0VIUQZJkHbQvqIhSj/rkYfsdDWVRFClGEStC2gaDPRH/4CzWN30B/+QnrbQgibkaBtAX3EQtR1YlE5qlDXiUU5vMjWVRJClFEStLPJbdza2MtWNTDsVKNqkI7+8DyQ3rYQwgYkaGeT27i1qZetMeyUrNKoUNeJQ3UszFbVFEKUYbJzzX+yj1vrNn+Buu1IVA5OKNH7UW5WRX812/b2ej0q7wO2q2wJkpyRyLLDk8nQxZm2aXN38cz13KT0BFZGTDfb0i2vc4UQuZOg/Z97x631EQvRBL6Gw8AVOc5VFIUs2dgXgF9PzqZhwFLThsgrDpHntm0rI6ZTv9l3Fp0rhMidDI+Q17i1zBKxxB19HBW8kwGo4J1Mhi42z3MzdLEWnyuEyJ0EbfIat46VOdkWcFb7cDveHYDb8e64aHzzPNdF42vxuUKI3MnwCOQxbq2gSt0PvGazetmDAc3fZfUptdmYdl6GtJvMikOYjWkLIfJHpSiKYutKWFtSUhKenp4kJibi4eFh9fIVRSE9PR1XV1dUKtWDX1BKSTvcJW1hIO1gUJB2sDRuyfCIKFaSw0WIwpGgLYqV5HARonBkTPsBlPgYSE3McdzZxR1cXW1QI/uV11x4IYTlJGg/SGoijGqS47Dqq1NQsYoNKmS/8poLL4SwnAyPiGIhc+GFsA7paZP3EAjlZIm1tdxvLrz0toWwnARtyHMIhAWni78udkjRZqJbOwpNv2/yHKOWufBCWIcEbVFoyuFF/80IaZ1nrzm3HC5CiPyToP0Ae9SV8VsYRR19stkNAMXFnbK7dCAbbSb6w/NkRogQxUSC9gP8FXuH+Yv1dBnlzb+3/+bKtfNkZmVRzcOLTx4fZOvq2ZzqWJghv7jMCBGiWMjsEWCPviKR3Uahd3DM8Zyj2oGdCzOY0S6Jq6sb8ET7ITzetjtN/KraoKYli6LNRHPiK5kRIkQxkqANxKQ7M+PSJNZ8fIPIhVHoF5w23IQs54m7swsAOi1s/SqZt1pcI+6mD/Xrt7VxrW1PObwIdd04yY4oRDEqcUF7zJgx1KpVC5VKxV9//WU6HhISQu3atWnRogUtWrTgs88+s9o1FTRsXZDK681jmTFbw5qsGkRWqIPi7Wcat9Y4QI/R7nxwqCI+leJo6WK1y9stJXof+stV0G6pY/rRX6mKEr3f1lUTotQqcWPaAwYMYNy4cXTq1CnHc/PmzaN3795Fdm1jbzo8LIV3N/sR38QZlRp6jC5P++EqutT0YvLaxfwweCS6rKwiq4e90AxYQaZkdBOiWJW4oN25c2erlZWYmEj2zLPOzs44Ozvneb7GAbq97E6fMR4093WgjpuGzDt6Jk2vyIG/DzPp18M80SSAJTu/4cDFvXw7omxPYzO2bSnM7ptv0hYG0g4GBWkHS88tcUH7ft555x0mTJhA48aNmTlzJrVr177v+TVq1DB7PHHiRCZNmpTjPEWl0OPV8vQe40HjClDTSUFNFnfSs2hdDkBHz9r1+OXwVuZvfJlKwTEk3axKRkaGFd+d/ZJ2uEvawkDawSA/7ZBu4b6zdhO0v//+e/z9/VEUhfnz59O7d2/OnDlz39dERUWZJRPPq6dd1TuTSTMqUsdNg/qe2deZ2kyW7FzMt3u/pEKnG7QYno7aAc5EqHFxcSnTwwKKopCRkWFZOyTcyDtVgJdf0VSwGOWrLUoxaQeDgrRDloVDrnYTtP39/QFQqVS89tprvP3229y6dYuKFSvm+RpPT0+Ldq5Zf3AdUQm30ajv3pfV6fXU8KpA0vXNhF9bQ8d3M3C4J96rVKoy/cE0sqQdlPukClB5Vy6imhU/+UwYSDsY5KcdLD3PLoK2Vqvl1q1b+PkZemS//vorfn5+9w3Y+dGsSnWaV/Hn2daBpmNhh3bx86EfWDtqAUvCW7P4oy+pEHSD6l0MPW0hhLCFEjflb/To0VSvXp3o6Gi6d+9O3bp1uXPnDo899hjNmjUjICCAr776it9++81q13w6oB2/nDpMlk4HQKZWy8Ijv5GQuBcnBydGdh/NnkkneJJpnPywFlE7XNHr9QW+XlJ6Agt2vc3nfz7Pgl1vk5yRy7BBKXO/BUxCCMuVuD7j/PnzmT9/fo7jhw8fLsKr6qnilMbSwzt5qX0o3x/fSYW2e0necffXFUeNI93TX+GFiSNYsnMRByrtKfDVVkZMp36z76jgnczteHdWHIJRnT+2xhspsWJSHVl8/i36jB9N8/3zqbPzW9RamTYpRH6VuJ62LYwJG8Uf52ez5MxG0jIzCTu1iSoBJ83OSYmGa3sUMmMcGdntVeY/922Br5ehi6WCdzIAFbyTydDFFqr+NpNwA5ebV+Df8yjR5+7+xMfkOFXRK2z9VsfroY7MOP8Wa8YfNfS8bVBtIexZietp28K8YYZx6/kHNjHg+7l4NNuH2sE8nMQd0VG7TQJxR7woX11TqOu5aHy5He9u6mm7aHwLVZ7NpCaierVZzuMLTsO9NxfVhjbTaWHrtzrCVzrx7u8fEu8C7YuhqkKUFhK0wTRu/UzwMJ5d9gn/bkpEn+RKaOrLHJquxcVTj1ql56HgFE7vLs+pL/WkJzhRub2emj3yH8CHtJvM8h3XyMi6ioujP892m1wE76qE+S9omy1g8nOkrlvhvgCFKGtkeCSb8s5urB0xmb3vHedJpnHN/yD+3VSoFD1NO8egUmH4U9FTOVhL9dCCNV95B1de+ucIY5tE8NI/RyjvYJ+7uufn5uKt9BRCXnZl0iEP/PudY93BFUxZvZjxG34uhpoKUXpITzsXxp73SEYDcOuozvwEjRqfdlpj5zHfSsuu5Pm5uejupqXPa2kMalAXNb5AZ74/sq94KyxEKSA97WwUReGf9focOQAUBa5HlufEH1W4HlmewqRVKE27kufn5uJTtX1ZtX8zOp3h2UytliWH/iT86FfFW2kh7JwE7WyMM0RSos2Pu/iooao3TV91gqrehscFdL9dye3OPTcX33rEie3tP+SwS5Ucpzo5ONC/WWt+PH4QgKWHd5HutYW/404Ua5WFsHcyPJJN9hki7v53m6be03f/XiVQQ+UOChbmdsmhVO1Kno+bi5naTJJvHWXJ6av0b9aGpWc30vjZI5z/tG5x11oIu1bmg3bUNh0xB5UcM0ROfqEnI1FN5fYqajxsvRkOpWlXcuPNxc4vO/Lv7b9Zd/A8v2i11PCqwOzeT5mdOyZsFOHX1lDrkYYMXDEXz1ymVQohHqzMB+3qoWocy+u5GaGnaRfDopCmnWM4FV4Z/24a/NrJCFJeyrtp6f1aGoMtuLlonAu/6Pf5POHSgH8Gncz1PCHE/ZX5iKTWqKgSqEGlMW8KlUZNlUANao1kKsvLUw/5sGrv7xbdXDTOyNky9ASd43sRP7MnUTtcTa8VQlimzPe0jYwzRGKvuuPrn1yoGSK2lJSewMqI6WToYnHR+DKk3WTcXTyL5FpODg483qg5Px4/yHOtO969uXgpPs/XJJxQUbdtArMdfuBPFhHhK/tJCpEfErT/Y5oh0kdN7BFvXLT2GbWLMxlVpjaT5NvHWHk+5r43F/O6b9BBGUGLxJeJ2qaz6n0DIUozCdr/uXeGCIHmz499ZwI3bibgV8mLz2d/WMy1s1xxJqMau+QV/vx3DQ/1vP/NRblvIIT1yP8WC924mUBS3b7cuJlg66rclzEZFVDkyajmDv2a/7V/n1u/JxIY0wC/JrnfXJT7BkJYj/S0H8DYw07TGX59T9NpGDz8VSp5e/DFJ7NsXLuchrSbzIpDmI1pFxUnByeGdxnJ4IavcvprNUdm9iS985953lwsLfcNhLAlCdoPYOxhG2kb9CEJ0EeusVmd7sfdxbPYN1RIOK626OZiablvIIQtSdB+AL9KXnBhLWk6DdoGfXA4tx43jY5KFQo2I2N3XCZ+5dTUcdVwZYPCQ72tvwGqos1Et+ZlNP2+QeXgZNWyjYw3F509NGjUWmpbcHOx3tMOht9cFhruDcydM7NI6iZEaSZBO5vsNxuNAcX456Bhr5AEuGl0/PDdV6QXcB37X9HxzP82i0f6u1J7tyfzrq3munMMNTy8mPPk81Z5H/qIhSj/rkYf0brIsgdWD1XjUE7PzcN6mna+AVh2c9H0m8uFtUVSLyFKO7kRmc39bjb6VfLC48JaQ8+7EBxRsXNhBtfnqajXNoE+bk/weJtuNPOyzlxqYxZBzWN3ijR7oOHmohrU5r8l5HVzcew7Exg07BWzewODhr3C2HcmFEn9hCitpKdN7jcbBw17JdceN5AjdaulorbpaLDXmwWvlsPDXaFK/Zuc/9UF/+t10GTUJ0pX+PnKxZ2rW9HDtcjyxD3g5mJe9wakxy1E/khPm7sBRdugD/BfQCmC6X3VQ9UkNb+DeznoOvAmKhWEDLiJq5JKrTY3C7wTjpEtcnU7VVJBlQenrTX+puJwbj0ADufWW+U3FyHKGulpk/fNRmsHFLVGRVojHU4nzJvdOfE6VZs7otLkzEOdH/fL1V1Uve2afXW4ujqhUqlyXZRklNu9gR/Dvi6SOglRmknQpngDSuVyWbg5pHD1jCMX/3GlTu10rPULjz3k6jZ+QUoPW4iCkaCdTXEElE7uWiKru4F7Gr/XWstr7o/gUq8KlLtT6LLtIVe3TPMTonBUSkHvqpVgSUlJeHp6kpiYiIeHh1XK3LorgvGffsusN1/k4eA2pKen4+rqavU51vZEURRph/9IWxhIOxgUpB0sjVvS036AlJQ0Xpr8KTtvXSTpKW9eWPwVXdbUYd7EV3B1dbV19YQQZYwE7QfoNmwcZwIVeNgPNZAyyI+NZ69zcdRkDv0y39bVKxHGvjOR2FuyylGI4iBT/rLZuiuCVn1HsXVXhOmYbyUv9L7mS8H1vk74yI00k9hb9pEBUYjSQHra5D0EsnjamzwZEkT46TXQwdl0vvp8Kk8Ed7VhjXP35fHzXIu/ytWr5wiM6cBe331o0Vt1iXx2Y9+ZSMzNeDL0ho9RbouShBDWJUGbvIdAug0bxy+fT8Fr/c8kdLh7vtc/WkKHBhTqmsZVmD4VPWnW05+Iy/tZ9NLyQpXpqvbmyuryDBwSRMXTjgwZ3JJ1V/+A9ORClZuX2FsJpNTrZ3osqxyFKHoStDEMgfzlm2o2VqT3dcLXpzz+1Xx5SPHk5qKbpucquXlSrYpPoa4ZczOe5Lr9uHLwcyIfXUhybLVClQdQXuPIH98k0iTahaHP6IkIL89fnrf5IKBJocvOjW9FL/SRa8jQOxTpoiQhxF0StCHPIZD+IY8AsGvFZ2bnG6fzFMT/3h5P5NUrpGjBA3D2cidjeSuUZF2B6w+GvCZV95ZjwauOZnlN+pTvzdHjDvgHWn8fxrlzPiQ9PZ3ho9+SVY5CFBO5EQl079Qar3+0Zse8/tHSvVMrq1/rr3/OoW4+FI9WQwFwa/AibnXfgazCTR/MK69JeX2yVfKa3I9vRetkQBRCPJj0tCHPIZDqVa2/v2LT2g2IPLnE0NNuNZS0c9+C5jY4Fq6nnVdeE4ekWKvkNbmfuXM+LNMLKYQoThK0/3PvEEhR+eJjw76STw8bRTJwJyEZzzFHyTpU+DHtyp7g5q3i2vlyXIx0oE49LYpPJSiXUuiyhRAlgwRtG6lcyRvVhbXUqR9AM3rnua9ifgRX9SSyihZVNWc6Pq4m9ogeV52Cyru8FWoshCgJSlzukTFjxvDbb79x5coVTp06RdOmTQGIjY3l+eef5+LFizg7O7NgwQI6deqUaxlFkXskO8mvYCDtcJe0hYG0g0FR5h4pcTciBwwYwJ49e6hZs6bZ8fHjx9OhQwciIyMJCwtjyJAhaLXaPEoRQojSqcQNj3Tu3DnX46tWreLSpUsAtG3bFj8/P/bs2UNISEieZSUmJpptDebs7Iyzs3Oe51vKWGYJ+yXFxJgLxLeiF3PnfFhk1ylIO2zddZgJn33LzDdepEfnNkVVtWJX0j8TxUXawaAg7WDpuSUuaOfm1q1b6PV6fHzuLmipVasWUVFR931djRo1zB5PnDiRSZMmWa1eGRkZVivLmmJuxpNSrx+686uZv3kuR6IOMv+5b4vsepa0Q0pqOqP/70t2J1wi+akKvLB4PsG/PMT8Ka9RvlzpyZZYUj8TxU3awSA/7WDp2g+7CNpAjnEhS76VoqKizMaGrNnTzsjIwMXFpUSN2xl72On/5QK5npLEj0s2otNm4TrS+oExP+3Q9flxnOmoQMPKqIHUwZXZfDaWPq+8z4Gfv7B63YpbSf1MFDdpB4OCtENWVpZF59lF0K5YsSIAcXFxpt72lStXcvSk7+Xp6VkkNyKNVCpVifpg3rhlWBpvZFzAk3BobpHW05J28PXJI1WAb/kS1YaFVdI+E7Yi7WCQn3aw9LwSdyMyLwMHDmT+fEP+6oiICGJiYvKcPVJW3c64xq2DH5N6zjAUknbuW9IuzAHHgi25t6YnQ4JQn081O2ZIFRBkoxoJYZ9KXNAePXo01atXJzo6mu7du1O3bl0APvroI/bt20e9evUYOnQo33//PQ4OdvGLQrHZ+P0vvDimL5kJ/2X109ym/nd7cWh727YVo3hTBQhRmpW4qDd//nxTjzo7Pz8/tm7daoMa2Q8nBydGdh/Nqc1XiTy5hIysLKJ2uKLT6Yv0usYZIbPefJEendvmek5xpgoQojQrcUFbFJ5xqXymNpMl4YusstoyNympaQyf8DF7Ei/n2DyifHm3HOcXV6oAIUqzErci0hpkRWTxaD/gNc501EPDbG18NpnGB1Qc/PlL21XMBuQzYSDtYFCmVkQK++Hj44Xe13wKpWHzCG8b1UiI0k+CtiiwJ0M6oj5nnkFQZoQIUbQkaIsC6x7UGq+LMiNEiOIkNyJFgflX86Wm4kH5RTcxjtrJjBAhipYEbVEoW7+dWeZvOglRnGR4RAgh7IgEbSGEsCMStIUQwo5I0BZCCDsiQVsIIeyIBG0b2borglZ9R7F1V4StqyKEsCMStItZSkoag96YzguLv+LcU268sPgrBr0xnZSUNFtXTQhhByRoF7Nuw8ax0f86KYP8UHs5kTLIj43Vr9Nt2DhbV00IYQckaBcz30pe6H2dzI5JkiUhhKUkaBcz2XZLCFEYErSLmWy7JYQoDMk9Usxk2y0hRGFI0LYB2XZLCFFQMjwihBB2RIK2EELYEQnaQghhRyRol0KyRF6I0kuCdikiS+SFKP0kaJciskReiNKvwEFbq9WydOlS5s+fT2JiojXrJApIlsgLUfoVOGi/+eabXLp0iVu3btGvXz9r1kkUkCyRF6L0szhof/TRR+h0OtPjmJgYpk6dypQpU4iLiyuSyon8kSXyQpR+Fq+I9PX1pUePHsyZM4dWrVoRFBREt27dUKlUdOjQoSjrKCwkS+SFKP1UiqIolp5848YN3nzzTapUqcK0adO4fv06qampNGvWrCjrmG9JSUl4enqSmJiIh4eH1ctXFIX09HRcXV1RqVRWL99eSDvcJW1hIO1gUJB2sDRu5WtM28/PjxUrVtC1a1d69uzJxYsXS1zAFkKUDIo2E+3Pw1C0mbauSqlicdDetm0bbdu2pVOnTlSsWJGNGzeyfv16hg4dyq1bt4qyjkIIO6SPWIjy72r0EQttXZVSxeIx7bfeeouNGzeSkpLC8OHD2b9/P/PmzePAgQMMHDiQP/74oyjrKYSwI4kpcfwQ9SsZjZvhEvUrQ5o9jUd5H1tXq1SwuKetKApqtRq1Wk32YfAOHTqwdevWIqmcEMI+/bDrdRr0OEHXHgdp0OMEK3aOtXWVSg2Le9qzZ8+mb9++ODk58cknn5gX4iBpuYUQ/9FmkpF1jQreyQBU8E4mI+saijYTlYPTA14sHsTiaPvoo4/y6KOPFmVdhBClgOpYGC4qLbfj3angnczteHdcVFr0EQvRBL5m6+rZPYuD9qxZsyhXrhz/+9//zI5/9NFHqFQqxo0rnvwWtWrVwsXFBRcXFwAmTJjA008/XSzXFkI8mOr6AZ5OSOSnbS3IcFZwuaPi6ds3UTL3AxK0C8viedoBAQEcOHAAV1dXs+NpaWl07NiR48ePF0X9cqhVqxYbNmygadOmeZ4j87SLh7TDXdIWBtIOBiVmnva9ARvAzc2NfKzPEUIIUQgWD4+kpKSgKEqObw29Xk9ycrLVK3Y/Q4YMQa/X0759e2bOnImPT+5TiRITE82+UJydnXF2di709Y1llvUvK2mHu6QtDKQdDArSDpaea/HwyPDhw3nooYeYPHmy2fEZM2Zw4cIFwsLCLK5cYURFRVGjRg2ysrJ47733OHXqFJs2bTI7x/hrxr0mTpzIpEmTiqWeQgiRH0lJSVSpUuWBwyMWB+24uDi6dOmCh4cHgYGBABw4cICEhAR27tyJr2/xJyW6fv069evXz9HTNwbtqKgoszdvzZ52RkYGLi4uZX7cTtrBQNrCQNrBoCDtkJSUhJeX1wODtsXDIz4+Phw9epSVK1dy5MgRAEaMGMHgwYNzHesuCqmpqWRlZeHl5QXADz/8QMuWLfM839PTs0huRBqpVKoy/cE0kna4S9oCnNMSUN3OZci0nCcq78rFXyEbys/nwdLz8rUqxsXFheHDhzN8+HCz43fu3LFKD/ZBbty4Qf/+/dHpdCiKQu3atVm2bFmRX1cIYTlVWiK8mksiuQWnoYwF7aJg8eyRTp06mf4+dOhQs+eMwyVFrXbt2hw7doyTJ09y6tQp1q1bR61atYrl2kIIURJYHLRTU+9uY3Xy5Emz58r6nWIhxF17lEpEdhuF3sHR1lUplQq0R+S9Qbqsj+EJIe66kebEjPNvsWb8UQneRcDioJ09MEuQFkLkRdErbP1Wx+uhjubB29YVKyUsvhF58uRJ07S++Ph4098VRSExMbFoaieEsD8qQ19Qp4Wt3+oIX+nEu79/SLwLtLdx1UoDi4P2xYsXSUpKYuvWrfTt21fSsQohchWXng6AxgGCX3Sl88uORN86xpWDsbR/bKCNa2f/LB4e2bhxI126dGHVqlW0bt2ao0ePUrNmTdOPEEIAODs50u0VNz4/XY0PPvLj5aY+NCzvSNPK1WxdtVLB4qD91VdfcfLkSQ4ePMju3btzbIQghBAAtSs5U/WJv+nt74jDDjVZWTqWHPqT8KNf2bpqpYLFYxyOjo5Ur14dgGbNmplNARRCCKMOniqiPVxY98dpqu1pyh7NUdK9tvD3pXhbV61UsDho37lzh7Nnz5qm+937uHHjxkVTQyGEXcnUZpJ8+yi3D7UmuE08R/7KovEbRzj/aV1bV61UsDhop6Wl0atXL7NjxscqlYp//vnHujUTQtidqG16LuxUqOnUHU8Pd6rUT6F1dC2aLjiCLrEiUdt01HhYY+tq2jWLg/bly5eLsBpCiNKgeqganYMjcUd9ad8tCYBO3bLYG+7KepdFdA+dYOMa2r8CrYgUQojcqDVQuYMaN9fyZseT76Sy2+EH1BpZmFdYMtlaiCKQlJ7AykP/R3rSGVw9GjOk/fu4u+TcmKO0UhS4Hlme2Kvu+PgnU14pR2PfAFtXq1SQoC1EEVgZMZ36zZdQwTuZ2/GHWHFIzajOH9u6WsXGxUcNVb1p2kdN7BFv/LQKi55ebutqlQplJmhnZWWh0+msUpaiKNy5c6fMJ7wv6nbQaDQ4OtpnsqEM7Q0qeBs2AqjgnUyG9oaNa1S86j2lMX0mqgRqoHiyN5cJZSJoZ2VlcfHiRaumkNXr9ajVckugKNtBpVJRp04duwzczuk6bse7/9fTdsc53TodBiHKRNA27nRTtWpVq+0RaQxWZb2nXVTtcOfOHa5du4ZOp7O7oK1oM3n60nFWOTYkQ6XGRdHz9KXjKNpMVA5Otq6esHNlImgbOTs74+LiUuhyJGgbSDvkTh+xEM/a13m5/qW7x7Qu6CMWogl8zYY1E6VBmQraQhQHJXo/ys2q6K9m+yLTK6hS9wMStEXhyKBsEQkPD+ftt982O9amTZtClTlo0KA8FzmFhYXRsWNHAgMDmTVrlun422+/TXBwMEOGDCEzM5PU1FR69OhB586d6dq1q6m8M2fO0KlTJwIDA9m+fXuh6lnWOQxcgeMrp3F8+a+7P6+cxmHgCltXrViNfWcCg4a9wth3ZEGNNUnQvo9MbSYLt89nxOJnbV2VBwoJCWHv3r3s27ePDRs2EBsby7Fjx4iJiWH37t00btyYX375BQcHB8LCwti1axcTJkxgzpw5AEycOJGwsDC2bNnClClTbPxu7Nvua4mcv3YNXfQ5lOw/8TG2rlqxunEzgaS6fblxM8HWVSlVJGjnwhisO80IYLVqMmdiTxSonHPnzvH444/Ttm1bs82Qp06dyoYNGwBYsGABS5Ys4fLly3Ts2JH+/fvTvHlzU29369attGzZkgEDBnDjhmHaWL9+/fj3338BQ8rcsLAwHnroIdPUO0dHRzQaDfv376dHjx4A9OzZk3379uHs7Ey1aoa8xo6OjqbNLK5fv069evXw8PCgYsWK3Lx5s0DvWUBMIsyYdoc1V/Uc+y4a3WstYFQTSC0bOzyNe+8DBg9/lTSdIcdImk4jPW4rkjHtbDK1mSwJX8TivV9SodMNmk9MR+0Afx8u2HdbfHw8u3fvJjIyknfeeeeB59+6dYtdu3Zx6dIl3n33Xbp3786UKVPYsWMH5cqVo1GjRgA8++yz/PDDD7z99tusXr2a1atXm8r47bffqFevHhUrViQhIYGqVasC4Onpye3bt03nZWVl8X//938sXrwYMN+s2XhupUqVCvS+y7rY1DS2fqvjn03lmf5UCN//7wqRqTvRHtjH7AENbF29Ihd7O5GUev1Mj7UN+pAEcGGtrapUqkjQzmZM2CjCr62h47sZOBR+ZiAtW7ZEpVJRv359YmNjTcezz7TIHiybNm2Kg4MD/v7+xMcbcg/rdDoqVKgAQECAYRlwnz596NmzJ/3798fHxwcPDw/AsI/n559/zsaNGwHw9vYmKcmQtCchIcFUDsDIkSMZNWoUderUATCba33vuSJ/nDWGHuYjLcpRr108e+c7UuWtR6jkUzZ+e/Gt4In6wlrSdBq0DfrgcG49bhodfpW8bF21UkGGR7KZN2wBbwZN56+PahG1wxW9tnDlHT9+HEVRuHDhgmkjZDAE06tXrwJw5MgR0/HcgrlGoyE+Pp7MzEzTEIuTkxMNGzZk4sSJPPusYbz9+vXrvPzyyyxfvhxXV1cAOnTowNatWwHYsmULQUFBAEyfPp2HHnqIp59+2nS9ypUrExkZSVJSkvSyCyFqm54G67xZ8GoVgpo7UaV+Cj26uNB0t5baP1QjalvpX2Qze/r7/PDdV7hpDO/VTaPjx7CvmTtnpo1rVjpI0M7GycGJkd1Hs2fSCZ5kGic/NARvnU5foPI8PT3p06cPgwcPZtq0aabjAwYM4Ntvv+Wxxx4jISHhvmX83//9H926dWPQoEH4+/ubjj/33HPs2LGDRx55BDCMk8fGxvLMM88QEhLCuXPnaNmyJZUrVyY4OJgzZ87Qv39/rl27xgcffMAff/xBSEgIEyYYxhk//PBDhg0bxiOPPMIHH3xQoPcrDKlJkwMyqVgJug68iUoFbbreoJw2Hf+2t6keWnb+y/lV8sLjwlrpYVuZSrHm2u4SIikpCU9PTxITE/Hw8CAjI4NLly7x0EMP5WtxjXGMO+Lyfha9dDfZTUlYVHLw4EFWrlzJ3LlzbXJ9KNp2KOi/ma0oikJ6ejqurq7suZ6E288OtO4aZ3r+5J+VaPZ8CirvyjasZdHL3g5lecFVQdrh3riVl7LztV8Axp539oBdEqxdu5YxY8bw+uuv27oqIhfBVT1x0jhy9Zwb4Rs8uHrOjaQ76Yz89e0Hv1iIB5AbkXaob9++9O3b19bVEPfh4qNG5+fBl7c/o8ZJLY19GhV46qgQ2UnQFsLKMrWZ/FlxEYv/NEwdTRqeziEH0Myub+uqiVJAhkeEsLIxYaP4dO97NH33MjVCDXP9y6KtuyJo1XcUW3dF2LoqpYoEbSGszNpTR+1NSmo6g9+cwQuLv+LcU268sPgrBr0xnZSUNFtXrVSQoF1ENmzYwLJlywBDXhAwzJVu2bIln376Ka+//jrp6ekFKjv70vehQ4fmed7169fp3r07HTt2NNUlu7Fjx9KlSxfatm3LTz/9BMCSJUuoV68eISEhDBky5L71+OmnnwgKCqJbt26meefZffbZZwQFBdG7d28SEw1LuPv3709wcDBBQUGmOeofffQRERGlpzdm7amj9uaxUZPZ6H+dlEF+qL2cSBnkx8bq1+k2bJytq1Y6KKVQYmKiAiiJiYmKoihKenq6cubMGSU9Pd2i1+uz7ihZq4Yq+qw7uT+v1ytarVbR6/V5lvHEE08od+4YXt+lSxdFURRl5MiRytGjR/PxTnL39ddfK2FhYcqlS5eUF154Ic/zxowZo2zbtk3JyspSOnTooKSlpZk9n5mZqSiKoiQnJysBAQGKoihKWFiY8sUXXzywDpmZmUq7du2UtLQ0Zffu3cqIESPMno+NjVW6du2q6PV6Zfny5cqMGTMURVGUixcvKoqiKOfOnVN69eqlKIqi3Lx5U3n22WfNXp/ffzNb0+v1Smpqaq6fiTtZd5Rvtn2pvLRoiA1qVrz0er3y6Ijxitv3ryvlN75r+nFb9rrS+5X3bF29YnO/z0Ne7o1beZGedi70EQtR/l2NPmJhgV6fkJBAeno6Tk6GXUqWLVvGH3/8wbp16xg5ciQbN24kJCSElJQUFi1axLRp09Dr9fTo0YPz589z+PBhunbtSnBwMB9/bNgMNioqik6dOtGrVy927doFQLVq1Zg5M+9VZhEREYSGhuLg4ECbNm04ffq02fPGHWFSU1Np3Lix6fhXX31FcHAwP/74Y55lR0ZG0rhxY5ycnAgKCuLUqVM5rh0SEoJKpTIlqwKoXbu26drGZFUVK1bk33//tdoeniVNSZ06WlQeD+6A+nyq2TH1+VT6hwTZqEaliwTteyjaTPSHv0Dz2B30h79A0Wbmu4zz589To0YN0+MaNWoQGhpKz549CQsL47HHHjM9N2LECA4fPszLL7/MY489Rv369Xn33XdZvXo1u3fvZu/evdy4cYPZs2czZcoUNm3aZJqs7+joSJUqVfKsh1arNeUUuTdhlNGQIUNo3rw5Dz/8MGCYTvjXX3+xadMmPv30U65fv55r2QkJCWYLAO4NuNmfz+3a77zzjlkSLV9fX6KiovJ8L8J+hAa2wOsf84F8r3+0dO/UykY1Kl3K6H3tvOkjFqKuE4vKUYW6TmyBt4jKzyq+UaNGMXjwYL788ksATp06Rb9+hixp8fHxXL16lQsXLtC6dWsA2rVrZ1G5jo6O6HQ6NBpNnkmgVqxYQXx8PO3bt+eFF17Ay8sLAHd3d7p27crZs2dz/WLInowKDDlS7n3+woULQM4EVO+//z7t27enU6dOpmNK6VuYW2ZVr+rDQ4onNxfdTZBVyc2T6lV97/MqYSm762lHRkbSsWNH6tevT7t27Thz5ozVyjb2slUNDDcIVQ3SC9Tbrl+/PpcuXXrwiRiWa8+aNYspU6Ywffp0wJDNb926dYSHh3P06FFat25N3bp1OXbsGACHDx/OUU5sbCyZmeb1bNOmDX/++SdarZYjR47QtGlTs+fv3LkDgJubG+7u7qjValMg1ul0HDx40DScER0dbfbaunXrcvbsWTIzM9m7dy/NmzfPce3w8HDAPFnV8uXLiY6OzpGqNi4uziy3irBvO5d/ypk1i00/u1Z8ZusqlRp2F7RffvllRo4cyfnz5xk3bhwvvvii1co29bI1huEHleZubzs/vLy8UKvVZGRkPPDcSZMmMXr0aN58803OnTtHREQEs2bN4sknn6Rr16706tWLjIwMxo0bx9SpU+nZsydZWVk5ynnzzTdzzOAYP348M2fOpHPnzrzyyiu4uLhw/Phxvv76awBTcqnQ0FAmT54MGGZ8dOjQgY4dO9K3b19q1aqFVqtl8ODBZmU7OjoyduxYunbtyuTJk3nvvfcAmDVrFpcuXcLHx4c+ffoQFBTEDz/8wKuvvgrAiy++yN9//01ISAjDhg0DDHnEq1atahrjLi1kuy1RFOwqYVRsbCz169fn5s2bODg4oCgKVapU4cCBA9SqVct0XkETRml/HoJy8zio79mQtVILs/39FAsSJW3cuJG4uLj7TsmzphEjRrBo0aIiKfvQoUOcOHGCESNGmB23pB0sMXv2bEJCQsyGfew5YZSxLQYNe4Wkun3xuLCWH8O+tnENi0du7VAWFaQdLE0YZVddm6tXr5r1yFQqFTVq1CAqKsosaBslJiaiKAqZmZno9Xr0ev19x041A/K+u5/X6/I63qtXr/s+b20LFy4ssmu1bduWtm3b3rf8wlzbOFSSvQzj3xVFsYvx7uz1HfvORGJvJeTYbsu3ohdz53xoy2oWueztUJYVpB0sPdeugjaQ41vrfm/UOIOjZs2aLFiwgPLly1O9enWr1UWvLxuLJR6kKNrB+CWbkZFhVwEgIyODmJvxuW63pY9cU+AFVfbGkqHBsiA/7WDpZ8Ougra/vz/R0dFotVrT8MjVq1fNptdlFxUVhYeHB5mZmcTGxuLj42O2rVZhGIcFyrqiage1Wo1arcbFxcVuhkcyMjJwcXGhciVvYnPZbsu3krdpV6HSKns7lPXhkfy2Q273qnJjV0Hb19eXli1bsnz5coYOHcqvv/5KrVq1ch0aAcP8YOOY9s2bN62WrD97z6+sfzCNrN0OxvKMO8zbC5VKxbyPDQueBg17hSTubrdVltjbv1tRyU87WHqeXQVtgG+++YahQ4fy4Ycf4uHhwdKlS21dJSFy5VfJC2S7LWFldvf7fYMGDdi/f79puXeTJk1sXaVCCQ8P5/z58wBm0/GyW7JkiWnhzf2EhYXRsWNHAgMDmTVrlun422+/TXBwMEOGDDHN5R42bBg+Pj5m5d65c4eXX36Z0NDQPDdZSE9Pp3fv3nTp0oUePXrkuspSGMydM1M2tBVWZ3dBuygp8TEo0edy/sTHFNk1swftFi1a8Morr5g9n5+bfCEhIezdu5d9+/axYcMGYmNjOXbsGDExMezevZvGjRvzyy+/AIaNfOfMmWP2+i+++ILHHnuMP/74g7Vr1+Z6jd9//52mTZuyc+dOBg4cyPLlZSOfhhAlhQTt7FITYVSTnD+pifkuKjw8nD59+vD444/Ttm1bTp48SZs2bUzPd+jQgfT0dJYsWcKECRMYNmwY4eHhvP22YR/BVq1a8dprr/HCCy+YXhMfH0/v3r05deoU33//PaGhobRq1Yrvv/8egIceesg0hubo6IhGo2H//v306NEDwCxxU25L0zdv3szu3bsJCQnhm2++yfV91atXj7Q0Q17khIQEfHx8TMmuunTpwsMPP2y2vL0sk00ADKuMtT8PK1AOH5E7CdpFKD4+nnXr1rFixQrTisPsXF1dGTp0KDNnziQsLCzHa19//XVTQI6Li2PIkCHMmTOHZs2a0b9/f/744w/279/PvHnzzF7722+/Ua9ePSpWrPjAxE3ZXb16lcDAQLZv387KlStzLF0HqFOnDn/99RdNmzbl+++/54knnkCtVrNu3Tp27txJnz59TLm5y6qUlDQGvTFdNgGg8BkzRU4StItQy5YtUalU1K9fn9jYWLPnHjT32Nvbm7p165oeL168mG7dutGoUSMAtm3bZhpXNg6vAJw8eZLPP/+cuXPnmsox9nzzShpl5OXlRbdu3XBwcKBjx45m5RotXbqUkJAQ/vrrL6ZOncq0adNITU1lxIgRdO7cmcWLF3Pt2rUHtEzp1m3YONkEAEhMiWNh1K982bgZC6N+JSklztZVKhUkaBeh48ePoygKFy5cwNfXF41GQ1JSEklJSURGRgJ3M/Hd6965z+PGjSMiIoINGzYA8MEHH/Dbb7+xZcsWypcvDxh2qnn55ZdZvny5aT5whw4d2Lp1K2CeuCk3QUFBHD9+HIATJ07kOZXSGPi9vLxITExk8+bNVK1alV27dvHSSy/Z1WKYouBbyQu9r5PZMb2vE74+3jaqkW2s3DmWBj1O0LXHQRr0OMGKnWNtXaVSwe6m/NkTT09P+vTpw40bN/j22285fvw4nTt3pkmTJlStWhWA0NBQ3n33XXbu3Mnjjz+eZ1kajYZly5bx9NNP4+bmRv/+/enatSstWrTA29sQDKZOnUpsbCzPPPMMYJge2bJlSypXrkxwcDA1atQwLRmfMGECv/32GzqdjosXL/LZZ5/x7rvvMnToUCZOnEjPnj1NGf6yGzJkCIMHD+aXX35Bq9WyaNEi3N3d+fDDD+nVqxdVqlQp89n6ngwJIvz0GujgbDpm2ATgERvWqngp2kwysq5RwTsZgAreyWRkXUPRZqJycHrAq8X92FXCKEsVNGGUEh+T+03Hcp6ovCvfPc+CREnh4eFs2LDBtPNMaWSthFG5seeEUdHX4uj47ngSnrmbP9prRSz758wq9Tmlje3gfOI7Fsb+SIPOZ6jgnczteHfO7WrMSN9BBcpPb28kYVQxUXlXhmzBuazbuXMn77//vtkxY45skTf/ar5lfhMAJXofT99K4qdtLchwVnC5o+Lp2zdRMvcDpT9oFyXpaRdAUfYw7Yn0tO+SlKQG0g4G0tMWQtiXhBsoFgw1ivyToC2EsD7jQrV7LTgtQ5CFJFP+hBDCjkjQFkJY3R59RSK7jULv4GjrqpQ6ErSzGbdhFYOWL2DIyoWmn0HLFzBuw6oiu2ZpyPJ3+fJlBgwYkJ+3LUq5mFRHZpx/izXjj0rwtjIJ2tk0q1Kd3o0CWPHMSNPPY42a06yK9bYou5dk+ROlkaJX2PqtjtdDHfkjahpr3v0veNu6YqWABO1sng5oxy+nDpP137LyTK2WX08d4emAdg94ZU5lKcsfQGxsLAMHDqRVq1bs2LEj3+0lSpe4//ZGrFnBmYZqD74YUp5NbaezMP8JM8U9JGhn4+TgQP9mrfnp+CEAfjpxiP7NWuPkULBJNmUly5/xtcuXL2f79u1MmTIlfw0lSh3n/5aqP9a2HPU7xDPlA18qVoxh484vbFwz+ydB+x7G3nZaZmaBe9lGZSXLH0DTpk1xdnamQoUKskt9GXftTw0N1nqwbHxVeoS4UqV+Cr6Jeir+6sn/rr5P1LacCdKE5SRo38PY235m5TeF6mVD2cnyB3D69GkyMzOJj4+XXerLuMqddSiBWtw9s2jT9QYqFTTrfIPybnr+dF1E9VD5fBSGLK7JxdMB7Th5PbpQvWwoO1n+AKpXr87gwYO5dOkSs2fPLlS7Cfum1sDKK6/yWOrbgJvpeKZzAvuyfkKtmWi7ypUCknukACTLn4HkHrlLcm4YGNtB46gh4rMkqtbM4spFFx6ql86FG4n8X+pgdr17wtbVLHKSe0TYhGT5EwXl5OCEX01v1FW0zE/4jBontTT2aYROJ/c7CkuCdhEJCQkhJCTE1tUolC5dukiQFgVW72kHwIFfOk8gU5vJkvBFNL4cYOtq2T0J2kKIIufk4MTI7qMZyWhbV8XuyW1cIYSwIxK0hRDCjkjQtrHSkDBKCFF8JGhnszsuk/NpWvQU3yxISRglhMgPCdrZxCRomTHxFqsvZRQ6eJe1hFEzZsygY8eOdOnShTNnzhS43YQQ9yezR7JR9LD1q2R2LEym20h3+oz1oLmfI3XdNKjJ/4KJ+Ph4du/eTWRkpGklYnbGhFFt2rShd+/eZtPrjAmj6taty5IlS0wJoz755BMaNWpEnTp1eO6557hz5w6dOnXiueeeM7323oRRxtWXliSMevXVV5k5cybdunXjscceo3p187S02RNGqVQq9u3bx4kTJzh8+DB79+5FpVJJ7hEhipD0tHOh0xqC91strrH9RAYRSdoClVNWEkb9/fffBAUFmVZ+Se4RIYqO/O/KhcYBeox255PjVeke4EJbj4L9QlJWEkY1bNjQNOwC+RuHF0LkjwyPZKNSG4J1nzGFGxYxKisJoxo2bEjLli0JDAzE1dWVL7/8ksaNGxe43YQQeZOEUdnsjsvEr5z6gcFaEkYZSMKouyRhlMG97TD2nQncuJmAXyUv5s6ZaevqFRtJGFVMgn2cbF2FEkUSRonCunEzgaS6feHCWltXpdSQoF1EJGGUKMvGvjOR2FsJpOk0AKTpNAwa9kqZ63EXBQnaQgiri731Xw/7P9oGfUgC6XFbgd0E7aFDh7J9+3YqVaoEwMMPP5xjRZ8QomTwregFF9aSptOgbdAHh3PrcdPo8KvkZeuq2T27CdoA48eP57XXXrN1NYQQDzB3zoeoVCr6DBoGgJKZyo8/htm4VqWDzNPOhaIo/LNe/8AFMNZQ2hJGZWRk2P1Yvii8lNQ0Br0xnauJt7hxbBFXE28x6I3ppKSk2bpqds+ugvann35K8+bN6d27t2kRyP0kJiaSmJhIcnIyer0evd4QiB/0k3wVru0x/Jnb89lZUt79fv7880/OnTuHoigEBAQwatQos+d1Op3FZXXp0oU9e/awd+9eNmzYwI0bNzh69CgxMTHs2rWLRo0a8fPPP6MoCjNmzGD27Nlmr583bx69evVix44drFmzJtdrbNq0iSZNmhAeHm6WMOre9ilsu1iznOL6sbf6FmU7dB86jo3+17n4v8acH16Xi/9rzMbq1+k2dJzN61eSPw+WKDHDI8HBwZw9ezbX544dO8aMGTOoUqUKarWaNWvW8OijjxIZGWlaDZibGjVqAFCzZk0WLFhA+fLlc+TSyE3sYT212yQSd8STctXy/l6738q/8PBwPvvsM1QqFTExMSxcuJARI0Zw8OBBADp27MiOHTtYunQpq1ev5tdff+W5555j48aNzJkzh7Zt2xIYGEhiYiKhoaEoisKtW7d44YUXmDFjBidOnGDJkiUkJiYyduxYnn32WWrWrGn6h3dwcDDlBunevTt6vZ4ePXqwdOlSBg0ahJ+fn+mDYnwfmzdvJiYmhk8//ZRBgwYxcuTIHO+rTp06hIeHo9friY+Pp1q1auj1ev73v//x119/0b59e7MyC8r4JZuRkWHxh7kkyMjIsHUVSoSKFTzQ+2aY9Qr1vk5UrOhKenq6zepV3PLzebC0XUpM0N69e7fF5/br14/x48dz7tw5Wrduned5UVFReHh4kJmZSWxsLD4+PnnmxYjapufGIQUXTz1qlZ46wSmc3l2e018pZCSq8WunosbDd19rXFSSF7VaTUJCArt27SIyMpJx48aZjgOoVCrKlSvHCy+8YJYwSqVSoVarcySMunnzJs8//zwff/wxjRo1ol69ejz//PPcuXOH4OBgnn/+edO1jQmjfHx8SExMpFq1aqjVary9vYmPjzerg/F6ANHR0bz66qvMmjWL7t2707t37xxfcvXq1eP06dMEBASYvhSOHj1KQkICO3fuZNu2bURERBQ6/4harUatVuPi4mI3i2syMjJwcXEp84trMjIyGBDaid1n1kEHZ9Nz6vOpDOzWw5RiobRStJno1r5M1iOf41Lew+LPQ1ZWlkXn2c3wSHR0tOnvBw4c4NatW2YJlXLj6emJp6cn7u7upiBgDFT3/vh3U+PfTYVK0dO0cwwqFYY/FT3+3QzPG8/NLq/yVCoVLVu2RK1W06BBA1PCKONziqLc97Xe3t7Uq1fP9Pjbb7+lW7duNG7cGJVKxfbt2wkJCeGRRx7h/PnzpvNOnTrF3LlzmTdvHiqVigoVKpCcnIxKpSIxMZEKFSrkeU0vLy+6d++Oo6MjHTt2JDIyMsc5y5Yty5Ew6sKFC7Ru3RqVSkX79u0f2C6W/lirnOL6sbf6FmU7dO/UBq9/zBOtef2j5eFOrW1ev6L8Sc5I5Jutw/jS8SJLdr1Cyp2kfLfdg9hN0B46dCjNmjWjRYsWvPHGG/z88894enparXy1RkWVQA0qjXmTqDRqqgRqUGvy33sqKwmj6taty7FjxwA4fPiwJU0jSjn/ar48pHjiv+im6echPKle1dfWVStSKw/9Hw06bqZr3wM0DNrCykP/Z/VrlJjhkQfZvn17sVxHUeB6ZHlir7rj659MYYZTy1LCKA8PDzp37mzqaQuxa8Vntq5CsUtPu0gF72QAKngnk5520erXkIRR94j8SUv5Gip8W6uJPaInJUqh3tPm323GG22SMEoSRhkpiiSMgrLdDoo2kwW/PkKDHieo4J3M7Xh3zm0NYFT/LagcHpzXSBJGFVD2AF0lUAOBNqyMjUnCKCEsp49YyOByp/nxcEMyUOOCnsHlTqOPWIgm0HqLAiVoF5EQSRglRJmiRO+n3E0fXlQbdopS9AoqKqGk7wckaAshRIniMHAFYBgiurRBj1/oHdzcrD9MJEFbCCGsQImPgdREkmOcubbHn4qVY6BKFko5T1Tela12HbuZ8ieEECVaaiKMakLcd9uo3SaeW2HbYVQTw3ErkqBtZ2bMmEHHjh1p3749S5cuBUCr1TJ06FCCg4MZO3YsYFgNGhISQpcuXXj00UdJSEgAoH///gQHBxMUFMSRI0ds9TaEKFWitunY921tjra/QnrVHlSpn0JalR6can6KiO9qEbUt51qMgpKgnYutuyJo1XcUW3dF2LoqZvR6PYMHD2bfvn3s3r2bOXPmoCgK69evp3r16uzevZu0tDT27duHh4cHq1evZufOnfTr149FixYBMGfOHHbv3k1YWBhTpkyx8TsSonSoHqrm30appKm8aPq4gkoFTZ+ARPdKHKl0mOqh1gu1MqadTUpKGi9N/pSdty6S9JQ3Lyz+ii5r6rB42puUL+9mcTn79+/n9ddfx8XFhS5dupCamsqRI0dIS0tj4cKFtGjRgpCQENq2bcuhQ4do3rw5X3zxBRkZGbz00ktcu3aN8uXLs3z5cm7fvs1zzz1H5cqVadu2rSmHiaOjIw4ODqbr9e7dG4CePXuyb98+OnbsaKpP9nONC2ayHxNCFI5aoyKpXiZO+8z/T2lcXCjX8GaBVlTnRf7XZtNt2DjOBCrwsB9qIGWQHxvPXqfbsHEc/PnB+ayNNm3axOTJk+ndu7cpU52bmxsnT57ko48+YsUKw11m4+47Q4YM4ciRI+zfv5/Q0FCGDx/Or7/+ysKFCxkwYADXrl1jx44dODndnaA/b948nnrqKVQqFQkJCabJ+J6enmY5rhMTE/nmm2/4/fffzer4zjvvmFZHCiEKz93JCRUq/jlZjgunylG3WSo6rSN1G3V88IvzQYJ2Nr6VvPjLNzVHOklfn7zTv+Zm9OjRzJw5k1WrVjFo0CCOHDnCli1bUKvVaDQa03nGDIXt2rXjwoULnDlzhoiICJYtW0ZWVhbBwcEABAQEmAXsbdu2sXPnTn755RcAvL29SUoyzA1NSEgw5QbJysrimWee4eOPPzYtdQd4//33ad++PZ06dcrX+xJC5E2l1hAZncX6g3r0HZJ5vFw5yrsn09XFeuPZIEHbzJMhQYSfXpMjnWT/kEfyVY6npydz584lMzOT1q1b4+bmxoEDBzh16hRjxowxnXfs2DG6d+/O4cOHef3117lx4waBgYE899xzgCHo/vvvv2bJo06fPs20adPYtGmT6bgxKVTnzp3ZsmULw4cPB+DVV1/lqaeeMgvOy5cvJzo6mm+//Tb/DSSEyJPKQcMp/1Tav6+iS00vJq1eyMO1Vfz06wEWvbTcateRG5HZdO/UOtd0kt07tcpXOd988w2dO3cmMDCQoUOH4ufnR9euXVm5cqXZeb///judO3fGy8uL1q1bM3LkSLZt20ZoaCihoaGm7HzZvfnmm9y6dYvevXsTEhJCYmIiffr04erVqwQHB+Pq6kpgYCD79+9n5cqVhIWFERISwty5cwF48cUX+fvvvwkJCWHYsGH5bCEhRF4qezkwaUYlXO9EMv7nr7gUvZn1mimciT1h1etIwqh7dB7yBjfTkk2PK7m558hWZo1ESSEhIWzYsOG+O++UdJIw6q6ynCgpu7LcDpnaTJaEL2LRnvnoqjWixXObcXDW8/fs+ux698GBWxJGFVBZTCcphCi8MWGjCL+2ho7jM3BwvlRk15HhERsJDw+36162EMLcvGELeDNoOn99VIuoHa7otQ9+TUFI0BZCCCtwcnBiZPfR7Jl0gn5M48SMmkTtcEWnK9wm1/eSoC2EEFbk5ODEyG6vsu3NA/RjGo19A6xavoxpCyFEETAG75dVo61arvS0S5C8kjzt3LmTwMBAgoKCOHHCcBfa0sRRQojSRYJ2CZJXkqf33nuPTZs28cMPPzB+/HgAixNHCSFKFxkeuUdSegIrI6aToYvFRePLkHaTcXfxzFcZ9yaM2rRpE4cPHwYMqxcPHDjA1KlTiYyM5NatW6SlpbF582a8vLxMZRgTOqWlpeHo6Ii3tzfe3t7Ex8cDuSd+elDiKCFE0Rv7zgRu3EygkrcHX3wyy+rlS0/7HisjplO/2XeEdP2B+s2+Y8Whafkuw5gwaufOnUydOjXP8xo0aMDmzZsJDg5m+/btpuPGJE9Dhw41SwYF4ODgQGZmpumxpYmjhBDF48bNBJLq9iX2tnU3PzCSnvY9MnSxVPA2rIis4J1Mhi4232XcmzAqu+wLUFu2bAmAv7+/qQd9b5InFxcXUzIoMIxbG5NHWZo4SghR9Iw97DSdISlcut6BwcNfxa+SF3PnzLTadSRo38NF48vteHcqeCdzO94dF41vvsvILWGUMZhGRkaazsu+zNcYzO9N8uTq6kpWVhbx8fEkJydTsWJFIH+Jo4QQRc/YwzbSNehDEsCFtVa9jgTtewxpN5kVhzAb086vb775htWrV5OamsrQoUPx8fGhc+fONGnShKpVq+b5OmOSp8jISMLCwujXrx9jx45l+vTp9OrVC7VazVdffQWYJ44CWLduHX369GHt2rUEBwfTsmVLAgMDC9YIQoh886vkBRfWkqbToG3QB8259ZTT6AzHrUgSRhVAUSZKsieSMOquspwoKTtpBxg07BWS6valfOQafgr72uJ2kIRRQghhA8Yed6UK+Zt1ZikJ2kIIYUVz58w0/cZRFGTKnxBC2BEJ2kIIYUckaAshhB2RoF0ELl++nOv+jvczdepUNmzYQExMDO+//36u52zatImOHTvSqVMnXnvtNdPxzz77jKCgIHr37k1iomEVVm5JpqZOnUqzZs0ICQnhrbfeKuC7E0Lcz9ZdEbTu9wo79hwtkvIlaOdi7DsTGDTsFca+M6FAry9I0DaqXLkyH3zwQa7PNW3alF27drFnzx5u375NREQEcXFxrF+/nj179jB48GDmz58P5J5kCmDmzJmEh4fzySefFKh+QojcpaSkMeiN6byw+CvOPeXGiO8XM/jNGaSkpFn1OhK0c2Fc2XTjZkKBXv/111/z008/ERISwgcffEBgYCAdO3bkyJEjALRq1YoxY8YQGBjIRx99ZPbay5cvM2DAAACGDRtGcHAwnTt35vLly9SoUcOUHMqYKCoiIoKQkBBUKpUpSVT2JFM1atQwLZEHmDJlCl26dOGPP/4o0HsTQuSu27BxbPS/TsogP9ReTqQOrszG6tfpNmycVa8jU/6yuTd3QJpOw6Bhr+Q7d8Arr7yCv78/b7/9Nk899RR79+4lKiqKl156ie3btxMfH8/o0aOpX78+oaGhDB06NEcZWVlZnD17lv3796NSqdDr725ZdOTIEW7evEnLli05e/ZsjiRReSWZGjNmDFOnTuXGjRs8/PDDHD582JTHRAhROL6VvPjLN9WsJ6z3dcLXx7p7wUpPOxtjD1vboA8A2gZ9CtXjvnz5MgEBAajVamrVqmUaby5fvjwNGjRApVLRokULLl++nOO1jo6OjBkzhuHDh/P666+Tlmb4FSs6OpqxY8eyZMkSIPckUdmPwd0kU8YEUn5+fjRq1Ijo6OgCvS8hRE5PhgShPp9qdkx9PpX+IUFWvY4E7Wz8KnnhcWEtDufWA+Bwbj0eF9bmO3eAo6MjOp2OWrVqcfz4cfR6PZcvXzbly05JSSEyMhJFUTh58iS1atXKUYZOp2PgwIGEhYXh6+vL6tWrSUlJ4ZlnnmHBggX4+PgA0KZNG8LDwwHYsmULQUFBZkmmoqKiTEmmjIE8LS2Nv//+mypVquS/kYQQuereqTVe/5hvwe71j5bunVpZ9Tolanjku+++47PPPuPs2bN8/vnnZjMk9Ho9Y8eOZdOmTahUKt58801effVVq17fOAQyaNgrJAFuGh0/hn2d73KaNWvGhAkTeO211+jRowdBQUGoVCrmzZsHGHrHn3/+OUeOHKFfv374+fnlKCM5OZknnngCvV6PSqXixx9/5IsvvuDixYumdvnggw/o0qULffr0ISgoCG9vb1asWAGQa5Kpd955h1OnTqHT6ZgwYQKurq4FaSYhRC78q/nykOLJzUU3UTDkYfEt50n1qvnPFHo/JSph1IkTJ3BycmLmzJm0a9fOLGgvW7aMpUuXsnXrVhITE2nVqhWbN2+mYcOGOcopbMIo49h2XmPZhU2U1KZNG9NONvZMEkbdJYmSDKQdDArSDnaZMCogwLDVvDE/dHY//fQTo0aNQqPRUKFCBZ566il+/PHH++4Mk5iYiKIoZGZmotfr0ev1WPId9fnsD01/f9D5Bf3OK0HflVZh7fdjLE9RFLtoq+z1LcukHQwK0g6Wnluigvb9REVFUbNmTdPjWrVqPbC3WqNGDQBq1qzJggULKF++PNWrV7danbLP6MiPgwcPFvi1JVFRvBfjl2xGRoZdBYCMjAxbV6FEkHYwyE87WJpgqliDdnBwMGfPns31uWPHjuHv73/f1+e208v9REVF4eHhQWZmJrGxsfj4+OTaiy8I47BAWVdU7aBWq1Gr1bi4uNjN8EhGRgYuLi5lflhA2qFg7ZCVlWXRecUatHfv3l3g19aoUYPLly/Ttm1bAK5cuWLqSefF09PTNKZ98+ZNsrKyrBJgZBMEg6JsB+PmxSqVyq7a2N7qW1SkHQzy0w6Wnmc3wyMDBw7km2++4cknnyQxMZGffvqJzZs3W/RajUaDSqXi2rVrVqmLoigkJibi6elZpj+YRd0OKpUKjUZj9XKLwp07d5gxYwaTJ0+2i98Mioq0g0FRtkOJmj2yfPlyxo8fT3x8PE5OTpQrV47169fTsmVLdDodY8aMMQXqN954w2x2SXa53YXNyspCp9NZpZ7Jycm0bduWiIgI3N3drVKmPSrqdtBoNDg6Olq93KKQmJiIl5cXCQkJeHoWzY4l9kDawaAg7WCXs0eeffZZnn322Vyf02g0pmRIBeHo6Gi1AHDnzh2uXLmCk5NTme9NSDsIUbzkTpoQQtiREtXTthbjiE/2/BvWZCw3KSmpTI9pSzvcJW1hIO1gUJB2ML7mQSPWJWpM21qio6MfOH1QCCFKoqtXr953PUmpDNp6vZ5r167h7u5epr/thRD2Q1EUkpOTqVq16n2nJpfKoC2EEKWV3IgUQgg7IkFbCCHsiATtfNqyZQutW7emZcuWNG3alKVLl9q6SjazefNm2rRpQ/PmzenQoYNp1/fSbsyYMdSqVQuVSsVff/1lOh4bG0vPnj2pV68eTZs2Zc+ePTasZdHLqx0+/PBDGjRogFqtZsOGDTasYfHIqx2GDx9OgwYNaNGiBZ07d+b48ePWuaAiLKbX65UKFSooJ06cUBRFUS5duqQ4OzsrSUlJNq5Z8bt9+7ZSsWJF5cyZM4qiKEp4eLjSpEkTG9eqeOzcuVO5evWqUrNmTeXUqVOm48OGDVPef/99RVEU5dChQ0qNGjWUrKwsG9Wy6OXVDgcOHFAuXLigdOnSRVm/fr0Na1g88mqHdevWmf79169fr9SrV88q1yuV87SLWkJCAmCYV1mxYkWcnZ1tWyEbuHjxIr6+vjRq1AiALl26cOXKFY4ePUqrVtbdXqmk6dy5c67HV61axaVLlwBo27Ytfn5+7Nmzh5CQkGKsXfHJqx3at29fzDWxrbza4fHHHzf9vUOHDly5csUqWTFleCQfVCoVq1at4sknn6RmzZp06tSJpUuXlskdzevVq0dcXBwHDhwAYM2aNaSkpOS6SXFZcOvWLfR6vWnvTjDkfI+KirJhrURJMXfuXNP2f4UlPe180Gq1zJw5k3Xr1hEUFERERAR9+/bl1KlTpp3OywpPT09+/fVXxo8fT3JyMp06daJx48Z2k+CpKNy7JkCR2bQCQyK8VatWFSo1dXYStPPh+PHjXLt2jaCgIMDwK3DVqlU5ceIEXbt2tXHtil/nzp1NO8HfuXOHypUrm4ZLyhrjjvdxcXGm3rYlOd9F6fbTTz/xwQcfsGPHDnx9rbPBrwyP5IO/vz/R0dGcO3cOgAsXLnDx4kXq169v45rZxvXr101/nzZtGqGhodStW9eGNbKtgQMHmjJRRkREEBMTQ6dOnWxcK2Erq1at4r333mP79u3W/fK2yu3MMmTlypVK06ZNlebNmyvNmjVTfvjhB1tXyWZefPFFpUGDBkqdOnWUZ599VomPj7d1lYrFq6++qlSrVk3RaDSKn5+fUqdOHUVRFCUmJkZ5+OGHlbp16yqNGzdWwsPDbVzTopVXO3z44YdKtWrVFCcnJ6VixYpKtWrVlNjYWBvXtujk1Q4ODg5K9erVlYCAANPPzZs3C309WcYuhBB2RIZHhBDCjkjQFkIIOyJBWwgh7IgEbSGEsCMStIUQwo5I0BZCCDsiQVsIIeyIBG1hlxISEpg9e3aBX1+rVi2z3MfWFhISYsolvXbtWg4dOpTnuUOHDqV69eq0aNGCFi1a0KZNmweWn5iYyMsvv0zt2rVp2LAhbdq0eWDu6jfffJMff/zR4vdw/PhxVq1alefzv//+O02bNqVp06Zs2bLFdPy7775j5syZZucOGDCAffv2WXxtcR+FXp4jRD7cm19ap9MpOp0u3+VcunRJqVixYoHrcW/uY2vLnkv6hRdeUL744os8z33Q8/fS6/VKx44dldGjR5va89ixY0rVqlWVzZs35/qa6OhopWHDhoper7f4OmFhYUr//v3zfL5169bKlStXlCtXriitW7dWFMWwKjQkJETJzMw0O/fIkSNKly5dLL62yJv0tMsolUpFSkqK6XGlSpVMaVVr1arFBx98QMeOHXnooYeYPn266bx///2XAQMG0Lx5c5o3b87kyZMBuHHjBv369aNZs2Y0bdqUhQsXml5Tq1YtZsyYQdeuXXnhhReYOnUqzz33HE8++SQtWrTg+vXrbNmyhU6dOtG6dWvat2/Prl27TK8PCwujRYsWBAQE0KZNGy5fvsyoUaNISEgw65nGxMTw1FNP0a5dO5o3b86UKVNMZezevZtmzZrRrl07XnvttTwz8C1ZsoQePXowaNAgGjZsSGhoKKdPn+axxx6jfv36DBo0CL1eD0BycjIjRowwXW/UqFFkZWWZlbdp0yZ+++03Zs2aRYsWLVi8eLHF/0YRERGEhobSpk0bWrVqxa+//grAjh07uHLlCp9++ikODoacby1atGDSpElm/1bZfffddwwYMMCUibBatWpcu3YNgCeffNKUBC09PZ0KFSpw48YNpkyZwvbt22nRogWjRo3KUaajoyNpaWmkpqaa0hO/8cYbzJ49O0e2x1atWhETE0NkZKTF71/kwdbfGsI2ACU5Odn0uGLFisqlS5cURTH0Ql9//XVFURQlNjZW8fDwUKKjoxVFUZSQkBBl9uzZptcZc0o89dRTyvjx4xVFUZQbN24o1atXVw4ePGgqb+TIkaZe3vvvv69Uq1ZNuXHjhqIoinLx4kUlMDBQSUxMVBRFUSIjI5WqVasqmZmZyp9//qnUqVNHuXbtmqIoipKamqqkpqbm2tPu0aOHsnPnTkVRDD36Rx55RFm9erWSkZGhVK1aVfnzzz8VRVGUn376SQFy7WmHhYUpXl5eytWrVxVFUZTHHntMad68uZKQkKBkZWUpzZo1U7Zs2aIoiqKMGDFCWbZsmaIoht7viy++qHz66aeKouS/p12tWjVTfoqJEycq8fHxSsuWLU3vOy4uTqlRo4Zy/fp15aOPPlIef/zxHOUcPXpUcXNzy/UaoaGhyqZNm0yPn332WWXp0qWKTqdT6tatqzRp0kRJSkpSNm/erPTo0cPUFvfrae/atUtp166d0q5dO2Xv3r3K+vXrlTfeeCPP84cNG6Z8/fXXeT4vLCOpWUWuhgwZAoCPjw+1a9fm0qVLeHp6sm/fPrZt22Y6z5iGdPv27aY9In19fXnyySfZsWMH7dq1A2DYsGFm+aZ79+5tSlW5efNmLly4kGMHkKtXr7Jx40aef/55qlSpAoCbm1uu9U1NTeWPP/7gxo0bpmMpKSn8/fff1KlTBzc3N9MOMk899RQjR47M870HBQVRvXp1AFq2bEmtWrXw9PQEICAggH/++QcwjFUfOHCATz75BDD0Ugu6Icb48eN57bXXTI83bdrEP//8w6OPPmo6piiKKcPkvbm7jc/nJTo6msqVK5sed+/ene3bt9OoUSNatmyJn58fO3fuZOfOnXTv3t2iOgcHB3Pw4EHA8FtHr1692Lx5M59//jl79+7Fz8+PTz/91NQmlStXJjo62qKyRd4kaJdRGo0GnU5nepyRkWH2vIuLi9m5Wq32gWXeG0iyPy5fvrzZc9kfK4pCz549WbZsmWWVz4Ver0elUhEREZHjV/P8bjh873vPqy0URWHt2rXUrl27wPXOi6IoNG/e3GyYyCgrK4t58+aRmZlp9iVx4MCBPLd6c3NzIz093fT44YcfZuLEiTRq1Iju3bvj5+fHjh07CA8P57vvvst3fSdOnMikSZOIiYlh3bp1/Pnnn7z//vusWLGCYcOGAYbPmDHvuCg4GdMuo+rUqWPqJa1evZrU1NQHvqZ8+fJ06tSJzz77zHQsLi4OMPTcjOPYcXFxrFmzhtDQUIvq0qNHDzZv3mw2m8M426JPnz4sW7aMmJgYANLS0khLS8PDw4O0tDRTAHV3dyc4OJhZs2aZyrh27RrR0dE0bNiQ9PR0UwD85ZdfSExMtKhu9/P4448za9YsUx3i4+O5cOFCjvM8PDzyfb2OHTsSGRnJH3/8YTp2/PhxMjMz6datG/7+/rz11lumax8/fpzp06czceLEXMtr3rw5f//9t+lx1apV8fDw4JtvvqF79+507dqV3377jX///ZcWLVrkq94HDhwgMTGRnj17kpqaavqyVqvVZvdNzp49S0BAQL7aQeQkQbuM+vzzzxk9ejRBQUEcPXrU4h7Q999/z4EDB2jSpAkBAQF8+eWXAMybN4+TJ0/SvHlzunbtyqRJk0xDIw9Sr149li9fzksvvURAQACNGjVi7ty5gGF3nPfee48ePXoQEBBAly5diIuLo0KFCgwZMoRmzZqZbkSuWLGCs2fP0qxZM5o1a0b//v25desWzs7O/PDDD4wePZp27dpx6NAhqySl//zzz3FwcKBFixY0b96c7t2757pH5nPPPcfKlSvzdSPS29ub9evXM23aNAICAmjcuDHjx483/Ubx+++/k5GRQb169WjQoAEvvvgiX3/9tdlwSnYDBgzg999/Nzv28MMPo1arqV27Nh4eHvj5+dG1a1dT0O3WrRupqakEBASYbkT26tWLw4cPm8rIyspi/PjxpiGi5s2bU7t2bZo2bcqff/7Js88+CxiGr06fPm3xF7nIm+TTFqIM0Ov1tG3blnXr1pnG64vTggUL+Pfff5k2bVqxX7u0kZ62EGWAWq3mm2++yfU3geK6/vjx421y7dJGetpCCGFHpKcthBB2RIK2EELYEQnaQghhRyRoCyGEHZGgLYQQdkSCthBC2BEJ2kIIYUckaAshhB35f685KTIkZ3+tAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 354.331x334.646 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mm = 1/25.4\n",
    "gp.layout(colors=gp.colors.bright)\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(90*mm, 85*mm))\n",
    "\n",
    "for model, pec in pec_results.items():\n",
    "\n",
    "    ax.plot(melt.loc[pec.index, \"FeO\"], pec, lw=0., label=model)\n",
    "\n",
    "\n",
    "ax.legend(frameon=True, fancybox=False)\n",
    "ax.set_xlabel(\"uncorrected melt FeO (wt.%)\")\n",
    "ax.set_ylabel(\"PEC%\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e75f0782",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}