{ "cells": [ { "cell_type": "markdown", "id": "175ddf88", "metadata": {}, "source": [ "# Performance benchmark" ] }, { "cell_type": "markdown", "id": "b6016dc5", "metadata": {}, "source": [ "MagmaPEC is a Python implementation of the algorithm also ued by Petrolog3. Compared to Petrolog, MagmaPEC has\n", "\n", " * a more up-to-date selection of melt Fe3/Fe2, ol-melt Fe-Mg Kd and melt thermometers and allows for PEC corrections\n", " * the option to perform PEC corrections at any pressure\n", " * the option to predict melt inclusion inition FeO from melt compositions.\n", "\n", "Additionally, MagmaPEC is open-source and users can extend it with their own models. As a performance baseline, we compare MagmaPandas results with those from Petrolog3. Petrolog and MagmaPandas have not implemented equivalent models for all parameters, but we settings are kept as similar as possible. In Petrolog we use:\n", "\n", "\n", "\n", "| Property | Value |\n", "|---------------|---------------------------|\n", "| Kd | Toplis et al. (2005) |\n", "| fO2 | QFM+1 |\n", "| melt Fe3Fe2 | Kress & Carmichael (1988) |\n", "| pressure | 1 bar. |\n", "| olivine | Putirka (2005, C-D) |\n", "| FeO initial. | 10.9 wt.% |\n", "\n", "Note that 1 bar is the only pressure available in petrolog. MagmaPEC settings are shown in the cells below.\n" ] }, { "cell_type": "code", "execution_count": 1, "id": "bfdd85f8", "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", "mm = 1/25.4" ] }, { "cell_type": "code", "execution_count": 2, "id": "6b39ecc8", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "################## MagmaPandas ###################\n", "##################################################\n", "General settings__________________________________\n", "fO2 buffer.....................................QFM\n", "ΔfO2.............................................1\n", "Melt Fe3+/Fe2+................kress_carmichael1991\n", "Kd Fe-Mg ol-melt........................toplis2005\n", "Melt thermometer....................putirka2008_15\n", "Volatile solubility model.......iaconomarziano2012\n", "Volatile species.............................mixed\n", "##################################################\n", "\n" ] } ], "source": [ "mpc.model_configuration.Kd_model = \"toplis2005\"\n", "mpc.model_configuration.Fe3Fe2_model = \"kress_carmichael1991\"\n", "mpc.model_configuration.melt_thermometer = \"putirka2008_15\"\n", "mpc.model_configuration.fO2buffer = \"QFM\"\n", "mpc.model_configuration.dfO2 = 1\n", "\n", "print(mpc.model_configuration)" ] }, { "cell_type": "markdown", "id": "c3edce08", "metadata": {}, "source": [ "Import all input files" ] }, { "cell_type": "code", "execution_count": 3, "id": "eb469397", "metadata": {}, "outputs": [], "source": [ "petrolog_file = \"./data/petrolog/petrolog_output.csv\"\n", "input_melt_file = \"./data/petrolog/petrolog_melt_input.csv\"\n", "input_olivine_file = \"./data/petrolog/petrolog_olivine_input.csv\"" ] }, { "cell_type": "code", "execution_count": 4, "id": "7f58445b", "metadata": {}, "outputs": [], "source": [ "petrolog = pd.read_csv(petrolog_file, index_col=\"SAMP_NO\").query(\"number == 3\")\n", "petrolog.index.name = \"name\"\n", "\n", "melt = mp.Melt(pd.read_csv(input_melt_file, index_col=[0]), units=\"wt. %\", datatype=\"oxide\")\n", "olivine = mp.Olivine(pd.read_csv(input_olivine_file, index_col=[0]), units=\"wt. %\", datatype=\"oxide\")" ] }, { "cell_type": "markdown", "id": "7d34fe96", "metadata": {}, "source": [ "Run MagmaPEC corrections at 1 bar and a fixed FeO(initial) of 10.9 wt.%" ] }, { "cell_type": "code", "execution_count": 5, "id": "b6d71509", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Equilibrating ... |██████████████████████████████| 100% [172/172] in 24.8s \n", "Correcting ... |██████████████████████████████| 100% [172/172] in 1:11.1 \n" ] } ], "source": [ "P_bar = 1\n", "FeOi = 10.9\n", "\n", "pec_model = mpc.PEC(inclusions=melt, olivines=olivine, P_bar=P_bar, FeO_target=FeOi)\n", "melts_corrected, pec, checks = pec_model.correct()" ] }, { "cell_type": "code", "execution_count": null, "id": "ff000b7b", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAAFaCAYAAADVfgw3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABOJUlEQVR4nO3deXyTVfb48c99krQpdIO2lF0GoTAMKIggsiuoIKBWUJkZF8AR5+fCuMFXwRmdcXdGVEYcRkbFZUZHlrKIouAGijgiiihQaEHK1hYKbVmaNslzf3+kTRuallSaJmnP+/XyZZM8Sc4T4+ntfe49R2mtNUIIISKCEeoAhBBCBE6SthBCRBBJ2kIIEUEkaQshRASRpC2EEBFEkrYQQkQQSdpCCBFBJGkLIUQEsYY6gJ/LNE0OHDhAXFwcSqlQhyOEEGdEa82xY8do27YthlHzeDpik/aBAwfo0KFDqMMQQoh6tXfvXtq3b1/j4xGbtOPi4gDPCcbHxwf8PK01JSUlxMTERMQIXeINLok3eCIpVgh9vMXFxXTo0MGb22oSsUm74kONj4+vc9K22WwR9UWSeINH4g2eSIoVwife0723XIgUQogIIklbCCEiiCRtIYSIIJK0hRAigkjSFkKICCJJWwghIogkbSGEiCCStIUQIoJI0hZCiAgiSVsIISKIJG0hhKgnf/vb37j//vvRWgftPSRpCyHEGdJa8/DDDzN9+nQsFktQ3ytiC0YJIUQ40Fozffp0nnnmGR577DFmzpwZ1PeTpC2EEGdg/vz5PPPMMzz33HP84Q9/CPr7SdIWQogzcOONN9KuXTvGjBnTIO8nc9pCCFFHZWVl3HzzzWzatAm73d5gCRskaQshRJ2UlJSQnp7Om2++yYEDBxr8/WV6RAghAnT8+HGuuOIKNmzYwLvvvssll1zS4DFI0hZCiABNnDiRjRs38sEHHzBkyJCQxCBJWwghAvTHP/4RwzDo169fyGKQOW0hhKjFgQMH+MMf/kBZWRkXXHBBSBM2SNIWQoga7dmzh6FDh7JkyRJyc3NDHQ4gSVsI0YRoRzFm7la0o/i0x+7YsYMhQ4agtWbdunV07NixASI8vZDMaTscDiZOnMjWrVtp1qwZrVu3Zt68eXTq1Inhw4eTk5NDfHw8ADfddBN33313KMIUQjQS2lWGe8VUdN5aVGIZujAKlToUy7iXUNaoasfn5uYydOhQWrZsyZo1a2jbtm0IovYvZBcip06dyujRo1FK8cILLzB16lQ+/PBDAObMmcPYsWNDFZoQIoJpRzG6cB8qsT3K7hn8uVdMRcVnYOlT5j3OzM7AvQKs6QuqvUZqairTp0/nxhtvJCUlpaFCD0hIpkfsdjuXX345SikABgwYwK5du0IRihCikdCuMlwZk3At6I356aW4FvTGlTEJ8/hhdN5ajLPLfI43zi5D5631mSpZv349ixYtQinFvffeG3YJG8Jkyd+cOXMYN26c9/b06dN54IEH6NGjB0888QSdO3eu8blFRUU+tWujo6OJjo6u8fiKY4NZ77Y+SbzBJfEGT0PH6l4+FZXgZzS9rAiVWOb3OSqxDH10L7TuwYcffkh6ejqDBg1i/PjxUHoMCvdBlRF7MAX6OYU8aT/++OPs3LmTefPmAfDGG2/QoUMHtNbMnTuXsWPHsnXr1hqff+rFgZkzZzJr1qzTvq/D4TizwBuYxBtcEm/wNEisjmKseZ9hO6/6aNqduQmtLfircm0etVEWncTKRYu4/vrrGT58OP95YwFli2/COPw5qoUTfdSGmTwY87K/g5/57/pSUlIS0HFKh/BX9t/+9jfefvtt1qxZQ2Jiot9j7HY7+/fvJykpyef+4uJiEhISfC5aQmAjbYfDgd1u907PhDOJN7gk3uBpyFh17lbMzy7DMqCg2mPuDUloVz+MNh/7TJGY2VHoonTes17N1VdfzRVXXMHLL79M7Cd3YSQs9XusJf3VoJ1DcXExiYmJFBUV+eS0U4VspD179mzeeustn4TtcrkoKCggNTUVgMWLF5OamlotYVeVkJBQ6wnWRCkV9l/6qiTe4JJ4g6dBYm3RAV3ofxSsC6OwXD8fc/V9uFatRSU60YU2z+qRK16ib94h7rnnHh577DGcxwogfx2GnxG7a9VaKD0WtKmSQD+jkCTtffv2ce+999K5c2cuuugiwDNC/vjjjxkzZgylpaUYhkFycjLLly8PRYhCiAii7PGo1KGY2RnVRsgqdShGbDJG+gKflSVvLlzGqKNFtGvXjqeffhqtNc7i/bXMfzs9z23do6FOy6+QJO327dvXOOm+cePGBo5GCNEYWMa9hHsF1UfT417yHqPs8ajWPXjqqae4//77efbZZ7nrrrsqXyS+XS0jdhsqsX2Qz+L0Qn4hUggh6oOyRmE9ZTR96lSG1po//elPPProozz00EPV24PZ41Gtah6xN8QqktORpC2EaFQqRtP+VDTgffrpp5k+fbrfY4xx/8R8t/YReyhJ0hZCNBm/+tWvmDt3LrfddluNxwQyYg8lSdpCiEbN5XKxdOlSJkyYwOTJkwN+Xm0j9lCSKn9CiEartLSUa6+9ll//+tds37491OHUCxlpCyEapZMnT3L11Vfz6aefkpGRQffu3UMdUr2QpC2EaHSOHTvG2LFj2bhxIytXrmTEiBGhDqneSNIWQjQ6VquVpKQkPvzwQwYNGhTqcOqVJG0hRKORn59Pfn4+PXv2ZMmSJaEOJygkaQshGoX9+/czYsQI7HY7mzZtwjAa5zoLSdpCiIjib/307t27GTFiBG63m3fffbfRJmyQpC2EiBDm8cO4l92CLvoWo6XL2+cxK+1eRl42mpiYGD755BPOOuusUIcaVJK0hRBhraIpr/nTUlSyA6U1OMFyiUL/lEHCF4dp3749GRkZtGnTJtThBp0kbSFEWKtoymtLLwMUoDCzNOZ6jWVIGck7t7L+k28xYhJCHWqDaLwTP0KIiKcdxf6b8nZR6GLQZRqV6ISi/SGKsOFJ0hZChC3PBccamhK0AI6FT53rhiJJWwgRtlRi+5qbEhwF81D41LluKJK0hRBhq7KNmG/idu/Q6OMxUHJ12NS5bihyIVIIEdYs416i+O0TNNv5DSSWoY9aUQl9sN4+HyM2OdThNThJ2kKIsPb8Cy/yf/+3nG+//IzubePDrilBQ5OkLYQIS1prHn/8cR588EFmzJjBL/tcgFIq1GGFnMxpCyHCjtaamTNn8uCDD/KXv/yFJ598UhJ2ORlpCyHCztGjR3n77bd55plnuOeee0IdTliRpC2ECBtut5vjx4/TsmVLtmzZQmxsbKhDCjsyPSKEOC3tKMbM3Yp2FAftPZxOJ7/97W+57LLLME1TEnYNZKQthKhRRbEmnbcWlVjmraxnGfcSyup/08vP4XA4uPbaa1m1ahVvv/12oy6teqYkaQshalRRrMnSp3IruZmdgXsFWNMX1Mt7nDhxgquuuorPP/+cZcuWMXr06Hp53cZKfp0JIfyqsVjT2WXovLX1NlWyevVqvvrqK95//31J2AGQkbYQwq9aizUlOj2Pt+7xs1+/pKSEmJgYrrrqKrKysmjVqtXPfq2mREbaQgi/ai3WdIaV9XJzc+nfvz8vvPACgCTsOpCkLYTwq6ZiTWb2mVXWy8nJYejQoRw5coQRI0bUR6hNikyPCCFqZBn3Eu4V4Fq1tnxKxOZdPfJzZGVlMWLECAzDYN26dXTu3LmeI278JGkLIWqkrFFY0xf47YBe1eker/Dggw9it9tZs2YNHTp0CGbojVZIkrbD4WDixIls3bqVZs2a0bp1a+bNm0enTp3Iz8/nxhtvJDs7m+joaObNm8fgwYNDEaYQTcbpkq6yx/u96FjbOm4sNu9xbrcbi8XC/PnzOXnyJKmpqUE9n8YsZHPaU6dOJTMzk++++46xY8cydepUAO6//34GDBjAzp07efXVV/ntb3+Ly+UKVZhCNGraVYZz4W9x/vNXuD8aiWtBb1wZk9Au/6tGTlWxjts66iCWAQVYRx1ExWfgXjHVe8yXX35Jz549yc7OJi4uThL2GQpJ0rbb7Vx++eXeql0DBgxg165dALzzzjvcfvvtAPTr14/U1FQ+//zzUIQpRKOmXWU45/aCA4sxUvKhsADVfD/ELvFJujU+P4B13J999hmXXnopKSkppKSkBOtUmpSwmNOeM2cO48aNo6CgANM0ff7jdurUiZycnBqfW1RUhNbaezs6Opro6Ogaj684tupzwpnEG1xNOV5XxmSM7jlYuinAM4AyszT6YCn6xFrMkqJa56f10b21rOMu47Pl/+bqSfcyZMgQMjIyaNasWVh/zqH+LgT6viFP2o8//jg7d+5k3rx5lJSUVKuZe7oT6dixo8/tmTNnMmvWrNO+r8PhqHuwISTxBleTi9dRjDX/M6z9fO82uihcmRqd4MCRlwWtflnza0QnYT1qw+LnIVeBhRsefoCRI0fy+uuvo5SipKTkzGJuIKH6LgT6+YQ0af/tb39jyZIlrFmzhmbNmtGsWTMADh065B1t79mzp1pirionJ4f4+MrRQCAjbYfDgd1uj4ii6hJvcDXVeHXRbtxJ/gdEqgXoQxbsqV1Q9piaXyQmBnfqMMzsDJ8pEnd2FNY2F7Fs1Z106dKFuLi4JvXZ/lxOpzOg40KWtGfPns1bb73FmjVrSExM9N5/zTXXMHfuXB5++GG+/vprcnNza109kpCQ4JO0A6WUiogvUgWJN7iaXLwtOkCh/8GNPqRQbYdixCSc9mUsV/iu4z62r4wtuQkMenwe59mivX89N6nP9gzeNxAhuRC5b98+7r33XgoLC7nooovo3bs3F1xwAQBPPfUU69evp2vXrkyaNIk33ngDqzXkszhCRLyqNbFr2u3ozgTNWVjTXw3oNSvWcVsnfcfruydy1qx8lrgux7DV/NeuODMhyYbt27evca46NTWVDz/8sIEjEqLxqmkttTH6Bcz3y0fJCaXoIxZUqyHYbn+1TrWytdb8+cnZ/PnPf2PmzJk8+uijKKXC+qJjJJMhrBCNXE01sc33qXG3Y6A7HAEyMjL485//zOOPP84DDzwQ1HMRkrSFaNQq1lJXTdjgWUvtWlW9JnZdOtVUJPYrLruI999/n1GjRgX9fIQkbSEatdprYpfhWnwTHNvsTdC61IbR4yDWPpUrGU7tVKNdZTiX3ULhjuXEtrdhOxHLyNShaNfF9dqCTPgnSVuIRqy2mtjmwRMYv/wIy6DKBO3OBHJN6Fq5RqHqqFzZ43Eu/R3u5otI+bUGHMCxem9BJmom9bSFaMRqq4lNmcbS3XdtsKUb6GLQZb4XESs61ZQczaNg+1Lsv/R9vL5bkImaSdIWopGzjHsJXZyOa1Ub3BuSca1qg5k7HFLdfo9XLYBjvvdVdKr58z1TIaHU//PKE7sILknaQjRyVddSG8M/wDrpO4iKRRX534GnDymIq7xdtVPN7+57mHid7P95Z9iCTARGkrYQTYSyx2NU1MQu+BLVUmNm+U5zuDNB6w64P27rHZU78i/n9pWKY8eO0eVXfYjucEm9tyATgZMLkUI0MRUrSox+CnO9xpWpIRE4ArosGsvERRhJv0AX7iOvxMLIsVdz+PBh7r7v/+jevXu9tyATdSNJW4hG7tSNMhUrSpRFYQwE9zqNzgeVDOS70OufhvRX2VvanBGXjKC0tJR169aRlpYGBN6CTASHJG0hGqlaN8qkDsXMzkAfcGC0VxjDK4oVaczsd3EsmcLQ6euwWq2sW7eOTp06VXv9mlqQieCSpC1EI1XT9nX3Cs+KEleGE52fgW1I9eV71p3refTBGYwYk07btm0bOnRRC7kQKUQjdLpWYLgcWIbMwmjXwu/zVaKT3467SBJ2GJKkLUQjVPv2dad3LlrXUFPblOV7YUuSthCNUG3b1yvWU9dYUzsrCkOW74UtSdpCNEK1bV+vWE+tHcUY/e8mb3t/9rxsUrDKivP91nAsXZbvhTG5EClEI1XTempj9Au4MiZ5V5UkG1a20oU25z+EtcsgjMR2oQ5d1EKSthCNVE3rqV0Zk3xWlViAnplgrvsd7o3JmDXUzxbhQaZHhGjkKravV0yJ+FtVYukGKqoUy8UHUPEZuFdMDVG04nQkaQvRhOjCfZREF/p9rKK6n5RZDW+StIVoQv760n8oyD7u9zF9FG91PymzGr4kaQvRhJzbfzDHYvpUX1WSpVHxoKI829mlzGr4kqQtRCPndrt5/fXXMU2TUaNG0ev/1nmbIrg+i8K1TKMPaoyBnoQtZVbDW0BJ+8cff2T16tXV7v/ggw/Ytm1bvQclhKiddhRj5m497byzy+XipptuYvLkyfzvf/8DfJsiWEauhbbj0SfaYX6dgmtVG3SxrNMOZwEt+Zs5cyYPP/xwtfsTEhJ44IEHWLp0aT2HJUTTUNfyptpVhrHyVtwFX1Sv3HfKEr3S0lJ+/etfs2LFCt5++20GDBjg87iyx6PanYtxzb+lzGoECShpZ2dn06dPn2r3DxgwgOzs7HoPSojGrtayqX7WR1ckVfe6x7C2fBdLv+qV+6p2Qnc4HKSnp/PJJ5+wdOlSxowZU2s8UmY1cgSUtN1u/w1AwfPnlxCibmorm1p1QwyxrTBX3+epzBfnQOcWYu1XvZSqa5VniV7FKNlms9GuXTvee+89Lr744gY9NxFcASXt2NhYMjMz6datm8/9mZmZNG/ePCiBCdFYVWxwqZqwoTL5Ohf+1tPDMbEMc89xjD5OrKNMdIHGVBp/l6IqlugdjWpNZmYmF154If/6178a6IxEQwroQuR9993HlVdeyQcffEBRURFFRUWsWrWK9PR07rvvvmDHKESjUlE2VZdpdIFGl1UZOccUoGzLsY46iHHeYVRsCZY00/NYXPlaar+vaeNQWRQXXXQRv/71rykr81+WVUS+gEba1113HSdPnmTq1Kns2+dZcN++fXseeughJk6cGNQAhWh0Ylth7jmOLtCoFp5ErOI9S+70YSeWC01AwbHyXYrlVJRCxXs6qBtdlPd+MzuKktjzGXbpWIqKili9ejVRUVI3pLEKuGDU5MmTmTx5MocOHQIgJSUlaEEJ0ZiZq+/D6OPEklb5h66ZpXF9oMCwoqKcnjv9jKyNgZ4O6s7NBkbbFuiiaI43O48LHv4fDpfJ2rVr6dq1awOejWhoAU2PvPnmm96fs7OzfRL2888/X/9RCdFIeeezK6Y8yhldFJRGQ1Tl0NozsvYkdO99FgVt7agO6RgXrcY66Tvy+j9KSpt2rFu3ThJ2ExBQ0p49e7b359tuu83nsddee61+IxKiEautDZjRJg6VeJ7PFnNjoMLcoXBm2HFvSPbsYjxyBZb0V9h+BE64DLp3787nn3/OWWed1VCnIUIooKSttfb7s7/bgZg2bRqdOnVCKcUPP/zgvX/48OF07tyZ3r1707t3b5599tk6v7YQ4ex0bcAsV873bjF3b0jGvbotqs11WG/ZiTH8Ayw3fYs55p98vek7Bg8ezMyZMz2vq5Tf1xSNT0Bz2lW/EKd+OX7Ol2XChAnMmDGDwYMHV3tszpw5jB07ts6vKUQ4OnWnYWUbsAyfmtYV9T6M2GQMP40LAIhNRmvN56tXM378eM455xweeeSREJ2ZCJWAkvahQ4d48cUXq/0McPjw4Tq/6dChQ+v8HCEiSW07HmtqA1a13kdNOxQ/+OADrr76agYOHMiyZctkn0QTFFDSHjlyJF9//XW1nwFGjBhRrwFNnz6dBx54gB49evDEE0/QuXPnWo8vKirymaKJjo4mOjq6xuMrjv050zqhIPEGV7DidS+fikrws+NxOVjSX8Vy1aueYk+F+6DKaPp0cWzevJnhw4ezaNEiYmJiwvpzlu/Cz3v/01E6hJ9op06dePfdd+nZsycAe/fupUOHDmitmTt3Li+++CJbt271+9zi4mISEhKq3T9z5kxmzZoV1LiFqJWjGOt/L8Q2OrfaQ873WuEa+QYkp0EdCjNlZ2dz9tlnA56yEhaLpd7CFeGhuLiYNm3aUFRURHx8zd+NgEba+/bt4+677yYzM5Pzzz+fv/3tb7Rs2bLegq3QoUMHwDNPfscdd3DfffdRUFBAUlJSjc/JycnxOcFARtoOhwO73R4RF28k3uAKRry6aDdmC6ffx1RMPrYvroXSZqhWQzHG/fO0DXTnz5/PbbfdxkcffcSQIUMi5vOV70LdOJ3+vzOnCihp33rrrXTr1o0pU6awcOFCZsyYUe91DVwuFwUFBaSmpgKwePFiUlNTa03Y4CkPW9tvpZoopSLii1RB4g2ueo23RYcaV4hw0sRySSEqqggzOwPzXd/qfKd67rnnuPvuu7ntttsYMmSIN8ZI+nwjKVYIXbyBvmdASTsnJ4eVK1cCcOmll9KvX7+fHxlw++23s2zZMnJzcxk5ciSxsbFs3ryZMWPGUFpaimEYJCcns3z58jN6HyFCJv5czOwC3xUip7T08ledr6rHHnuMBx98kBkzZvDkk0+ilIqY+WERPAElbZvN5v25PubS5s6dy9y5c6vdv3HjxjN+bSFCpeqKEeJLcW9UuL+3o1Lt6NxCjJTKll4VKqrznbpSpKSkhIULF/KXv/yFBx98MKJGqiK4Akrau3fv5tprr63x9jvvvFP/kQkRYfzVyHZvt2Hm9QX7D1iG5FV7jnnwOEZsq8rbpklBQQEpKSl8+eWXxMTENEjsInIElLSfe+45n9un64IhRFNTU41sS3cn+qetqJRBuHcs9ak5YmZpVLTDU0AqfQFut5tbbrmFdevWsWXLFknYwq+AkvZNN90U7DiEiGi11RRRiU7oexvmwpXonSVVyrGC5TJwr15L2bECbrzldhYuXMhrr72G3W5v4DMQkSLg0qwVlixZwnfffYfD4fDe9/TTT9drUEJEmtPWFLHY4KxYjPMccAyIq7wgSUIZd02+hiXLP2fhwoVcffXVDRe4iDgBFYyqcNddd/Hqq6/yr3/9C7fbzdtvv01BQUGwYhMiYlTWFPFN3BU1RVRqd89W9iiFSlKVCRsozYeln3zNsmXLJGGL06pT0v7oo49YtmwZKSkpPPPMM3z99dfk5+cHKzYhIopl3Es+Ffpcq9qgi9M9HdZrSer2jiP5YeceRo8eHaLIRSSp0/SI3W7HMAyUUjidTlJTU9m/f3+wYhMiLPmrwFdxn2X0HM/tUyv0gU+hKJ1QyqGdRRx0/4K+f3yJlqfZFSlEhTol7bi4OE6ePMngwYO56aabaN26tc8abiEaM3+V+0gZhEKhD31erZpf1e3pVZN6Xl4et143hi17y8h477XTbmMXoqo6Je233noLq9XKX//6V2bPns3Ro0dZuHBhsGITIqz4W4ftem8RdDGw9q2ylC87A/cKz/b0UxO9q8DCF18fY/N+O6s+Wkf37t1DcSoigtUpaVfUBQGkkp5oUvytw9ZlGrTGkuZ7bNXt6e73p/kkegsw6hdw8cUX00IStvgZ6pS0t23bxmOPPcauXbtwuVze+//3v//Ve2BChFrVuWvPv0t9DzgGqoX/56pEJzpvu98NN816QNSqTTXWHBGiNnVK2tdeey033ngjU6ZMkXq+otGqPndtQzus4DqCz7c+zrNJxu9rHLViFuyGOIffx2uqOSLE6dQpadtsNqZPnx6sWIQIC35riOzUmN9rzCyF0aW8PGqUQqNwZ4KlW5XnZ4IuPoracQ/64FHc60yMgQplqVybrQttqMT2DXZOovGoU9IeNWoUq1atYtSoUcGKR4iQqrGGSFeFztSY+zRmpka1APOgBYy26Oz9uLJM7/Z0DajWJVgGlGIB3DsV5nqNZYgnaXs33MjUiPgZ6pS0R4wYwZVXXonFYiE6OhqtNUop2WAjGg1dsBti/DerVi3B6KbQ0RrygdJ41MlSrJeDLlM+29NdK010mUZFKSxdFc7vDfQXLeBYdLUmvkLURZ2S9q233sqCBQs477zzZE5bNDraVYbrnQngLsXfZmF9BNxbNJwovwCZVwRtPPsUVJSCKk2WVAs8Sbz8PqNtC1SveRhdhskIW5yROiXtpKQkJkyYEKxYhAgp15KbMHrshXxP2dSKuWsAd6aJdoCle5U57T4m7pU1JPijQFyV24VRWCRhi3pQp9oj6enpzJs3jyNHjnDy5EnvP0JEOu0oRu97D0s3hTFQoQ9qXCtN3OtNnItMzM0GaOWTyPVGUM08Cb4qd6bp01bMvUNDQm9J2KJe1GmkPXPmTABuu+02b786pRRutzsowQkRLBVrsEloB9jQedtRLTwXH5VFYRmiPJtnjoF2a3ReAqqdAXjW+OkyjS4GyyjPRUZX+cVJfRR0IZ5R9nrTc/t4DNbb54foTEVjU6ekbZrm6Q8SIoxVrME2D36KinGgS+wYKUPQ593qSbZVVMxT6wINyb2haEflg+Uba05N8MSBuVGjOoOyKcxDUaiSqzBikxvwLEVjVqek7W8qpFmzZvUWjBDB5lp2M+QtQpkm2ECdAOPQIsyvnVAWjXtnCZaulbOG7kwTjlmw/u5NzNX3YWZneDqsn7KxpuqFSPOggdKyUkQER52SdmxsbLWu0Dabjf79+zN//ny6detWwzOFCD3tKEZnrcDoY2JJq5KYd5iY365CdRqN3rkUV6YJCaALAIuB6jEeIzYZNe4l3BlOXNvXQUs3pUeOojK1z8YaMzsKo+NYjCGzqpVmFaI+1ClpP/LII8TGxjJ58mS01rz22mucPHmS1q1bc+utt/Lpp58GKUwhzpyZtx0sJT4JG8CSZmB+XwI9r0ev/gaO74ETeEbOhQpMjek4jvn+HeiCL9HxTsq2F+K2mji3WDF/0KhkG5QkeUbW6S9JuVURNHVK2kuWLOGbb77x3p42bRqDBw/m888/55lnnqn34ISoL9pRjM79ERL9P65agl46ERKdWAZUXSWice9YgmvuOoyzj2C52In5lSZmiMLoYgCeTTTmVhPd7EKs6Qsa5oREk1WnJX8nT55k165d3tu7du3y9oi0WuvcI1iIoNOuMlwZk3At6I3OfgAKwfWpiXnIk2y9xx0Bmpeh8F2fDWBJ02Dmoo+V4f5QY+4H9YvKx1WUwtLbBQVfoh3FDXNiosmqU6Z99NFH6d+/P3379gVg06ZNzJs3j+PHj3PNNdcEJUAhzkTV4k/arTHXg1kA5o/ACY2K15AKnATaetZd+6PagOWXnqa87h2mTy0R7zFSuU80gDol7fHjxzNkyBC++uortNYMGDCAVq1aAZVruIUIF+bxw5g/LcWW7ll/ba7XqDYK2xDfnY7mJqAFUAi61O9LQSHeHY6WNAPXzsraIhWkcp9oCHWaHgFo1aoVo0ePZuTIkcTGxsqOSBG23MtuQSV76llXbIapNvXRzYAosIwEXJ4KfdV2OO4wUQn4JGji8KzLLieV+0RDqdNI+3//+x8333wz27ZtQ+vKL7bsiBThRjuK0UXforQGVO1dZlqA+TnQFjgM7q807u81KtGz7E+lgmWQb7LXeQrXITtGu1ifZr5CBFudkva0adP417/+xe9//3vWrl3LnDlziImJCVZsQvxsunAfKqEMTnqmQIxfqBq7zHAcVFdPLRFsQAfgELh2g0oGW1vfBgbuHQqj09VYxv3D245MRtiiodRpesTpdHLBBRfgcrmIi4tj1qxZLF++PFixCfGzaFcZ5rrH0LmF6GIwd4BruUaXlu9wrMKdaUI8WLobYAdrusI21MA23sA22MB0tsP9jR3nEnB9As4MOxybgCX9VZQ9HqN1D0nYokHVaaRdsawvKSmJ7777jvbt27Nnz56gBCZEoKo24FX2eM+KkZbvYj3XxL0arGMMXJ+a6BjQe8C10/Ssyz4E2gXWq8pLrbbBpwa2pRvYdoNl4o/oov0oQKV2lyQtQqpOSXvixIkUFBQwc+ZMhg4disvl4i9/+UuwYhOiVtUb8EZB0oVQsN6zxK/AM1+t3drzN+V+oAXoItC7gLZgHVZl6qMQnxrY4FnGh6MIy1n9G/TchKhJnZL23XffDcCll15KQUEBDoeDuLi40zyrumnTprF8+XL27NnDli1b6NmzJwD5+fnceOONZGdnEx0dzbx58xg8eHCdX180DX4b8H63HN2sPAmXF3Uy12uMtgpjaNV5aRPy8CZsd6b2qYFdQZbxiXATUNLeunVrrY/36FG3zQQTJkxgxowZ1RLy/fffz4ABA1i1ahVff/01EyZMIDs7W3ZbCh+e5rvbMQ9+iu2UBrxGDyeu5ZXd0mmu0UeothHGkmbg/N5ErzPReQZ7D2g6jISqTfTcOwxUymCZDhFhJaBs2LNnT9q3b4/VavVZ6geglPLZ2h6IoUOH+r3/nXfeYffu3QD069eP1NRUPv/8c4YPH16n1xeNU9XpEOJLoLQQ9zoTY6ACN9561mgTd6ZnDbbxq/Ldj36oNqBNMFUbiuPi4adMXLsru6qjgFTt/8lChEhASfv666/nyy+/ZMKECUyZMoWuXbvWeyAFBQWYpklKSor3vk6dOpGTk1Pr84qKinx+kURHRxMdHV3j8RXHnvrLJ1xJvJXcy6eiEk6ZDtkJrqUaFY2nGNRRoAzMjaCzPCtDtP/m6uiDniV9Ro9D/Gr7EWzVuqqDa9UXmCVFYTPajqTvQyTFCqGPN9D3DShpv/766xw7doy3336bG2+8EZvNxpQpU7j22mvrtQnCqbW6AzmJjh07+tyeOXMms2bNOu3zHA5H3YILsSYfr6MYa95n2M7znQ6xdDUwfzSxXKoqezJmmp6krcHSHUzTT6PenSYqBazDDHSBE3euC6jeVZ2EMhx5WdDql/V7Pmcokr4PkRQrhC7ekpKSgI4LeLI4Li6OW265hVtuuYUPP/yQ66+/nvz8fGbMmPGzg6wqKcnzf8qhQ4e8o+09e/ZUS8qnysnJIT6+chQUyEjb4XBgt9ur/ZIIRxJv+esW7cZs4fT7mNGKU5bqeRK56gzuD4F24N5YvssxFXQeqHiwXFTlgmVhDW9cFIU9tQvKHh6byCLp+xBJsULo43U6/X+/TxVw0na5XCxbtoyXX36ZnJwc7r33XiZPnvyzA/TnmmuuYe7cuTz88MN8/fXX5Obmnnb1SEJCgk/SDpRSKiK+SBWafLwtOniW9Pmhj1J9qV4LMFqCOwpUR7AMLI/lGLiOaCxDK5f6qSgFZdGY2drTSqxcRT0RIyah/s6jnkTS9yGSYoXQxRvoewaUtO+55x6WLl3KxRdfzIMPPsjAgQPPKLjbb7+dZcuWkZub6y08lZWVxVNPPcUNN9xA165diYqK4o033pCVIwIAZY9HJV2IO3OJT3svd6bpf6le+ZZ1lQSWhCpTJ4cVGKbP8e7sKNTZ49DFFlyr1paXWLVJPRERlpQOYOLYMAy6devmt0ckeApJNbTi4mISEhIoKiqq00hba01JSQkxMTER8dtf4q3k3r8Z9ztDUFGl3hUe2gFGL3xaiJlZGvdGjfVqhWulHaITKW12nOiSeIxWg9FoOPSFNzm7kwYRddW/MGzR1XZXhptI+j5EUqwQ+ngDzWkBDWM/+eSTegtMiJ/LSPoFpr0llhG5lcv7LJ7NM87vTc829EJAlS/n2xuN6ngl931osHrxazz36j+5dNx4oHLrOwntcGqbt6ejssdLEwMR1gJK2sOGDQt2HELUqOrol4Q+uLPew5JaPh9tUZ611Dme7emcAGzRoJJwFw7llgwnb771X1555VVvwobK5Ky1hgCv2gsRDmTCWIStU2uLmEdtUGLAfnDnac+8dYmGaFCjQX+UgDHxPYyoZrhjW3Pdjb9jxYoVvP3229IOTzQakrRF2KpWW2SdCZ2Vp9tMxTE7TcxvQH8MKrkYc/FI6HgFlivmk5aWxtKlSxkzZkyIzkCI+idJW4QlT32Rtd6Ercs0ugisp9YQ6Wqgt5pYRlesECnFnbkQ812DJ59c0PCBCxFkde4RCZ4126+99hpz586lqKiovmMSwnORMM6BLtDoMu258Jjo/1iVgk+/Rks3MA9+gnYUN0CkQjSsn5W077nnHnbv3k1BQQHp6en1HZNo4rydZw4W4t6mca/SuDZrz8oQf8f721yTWOZJ/EI0MgEl7aeeesqneW9ubi4PP/wwf/rTnzh06FDQghNNk3vZ7yAhA9sEjXWwgfUKA6OdZ3WIv3Zh/jbXUGSXOtiiUQooabdq1YpLL72UTZs2ATBo0CBGjBjByJEjGTBgQFADFE2LdhRj5izHkua758vSzYAoML9phjMj2tOvcYnG/Ap0iu9ruHcYqNbDwnJzjBBnKqALkZMnT+byyy/nnnvuoU2bNjzyyCOMGzeOEydO0KtXr2DHKBqxU3cgmnu/QSWU4G88oZJAF8Zgnfw1FO1nz+5skve8Q/TmT3BuLd8lWWjH6DgOy5Wy/Vw0TgGvHklNTeXf//43K1euZNSoUTz44INccsklwYxNNGLV+zva0KVR4CqCGgor6CNAWQHu9+5g2y8fZOS1f+Diiy/mrQU5mHnbpfGuaBICmh5ZvXo1/fr1Y/DgwSQlJbFy5UpWrFjBpEmTKCgoCHaMopHRjmJc71yHis/AOuoglgEFWEflotJ+AutRcHvqh1RlZmnQYLkYSFnDD08PoUOHDrzwwgsoezyWs/pjnNVfErZo9AIaad97772sXLmS48ePM2XKFL788kvmzJnDhg0buOaaa/j444+DHadoBCpG1+bBT8GRh1FeLFK7NeZ6DcWeOtfaCe7NGvd2jdHSszpExQMGqJYKa2sXw35pYczzy0hITg7lKQnR4AJK2lprDMPAMAyfbjIDBgzgww8/DFpwonFxr5gKsUtQtlKIgYo/9Mz1GtVG+TTfde800QfA6KYgDsyfNEpXrhJp1TURiynrsEXTE1DSfvrpp7nqqquIiorimWee8X0BqXctAlCxw1EdL0V1AJ1dfn+ZRhf76Zbe1cC52cS1SaMcoBLBMqjKMUVRsqRPNEkBZdzRo0czevToYMciGjHvDsdCsA4xcB8yMbM0qoWny4w/qhXo/Snocwux9qjcJ1DRUUbmr0VTFFDSfvLJJ2nevDl33nmnz/1PPfUUSql66xMpGi+V2B592IJRvqZanQ/mFxq9FXCBxc9zdL6FdxLvwvXyA4w+P5aWv2gOhVHSUUY0aQGtHnnrrbf43e9+V+3+O++8k//85z/1HpRofJQ9HtV2KGY+uNeZuNeAdgMOoNTPTsedJq5Sg9vums7GNreS/IfdWIZ/iHXSd1jTF3ibFgjR1AQ8IR0TU70bdbNmzQigW5kQAFjTX8X51HuY1pOeZgVOT7EnXQDmBjB3mBhJlatFTraEh+66mXue/Lun/VOzxFCfghAhF1DSPn78OFrran3TTNPk2LFjNTxLiFO4HGDVKAOM8xVGlyqrRTJNzM2Vq0VUlCJ+VSvufviZiOgvKERDCWh6ZNiwYTz66KPV7n/iiSekFZkImJm3HawO0PgkbCivLaLwVPOzgDsrCiN1KEZMQoiiFSI8BTTSfuqppxg2bBgrV67kwgsvBGDDhg0UFhby2WefBTVA0XgogGagasjDqrXnccfbVqJ6pMvFRiH8CGiknZKSwqZNm5g6dSplZWWUlZVxyy23sGnTJlq1ahXsGEUEMgv34/5hBWbhfu99KrU7nIz21L/2pwgsvRS2lilYRs+Ri41C+BHwhUi73c6UKVOYMmWKz/2lpaVER0fXe2AiMpmO41heOg+3cx/EayhWENUR663fYNjjMTpdiZn7X8ws7TunvcNEJXjmslULl6fyX+seITwTIcJTQCPtwYMHe3+eNGmSz2MV0yVCALj/2Rdly0E1MzHiNKqZCdafcP2jN2buVozLZqNaTcD9lcK52MT1iYlziYnOBWOgJ4nrQpvsdhSiBgGNtE+cOOH9+fvvv/d5TJb8iQpm4X44sQfLBQrVEU/fxjjQOeDesBfXhxejTjZDpQ7l5JQf+fGR/pzX/jhRgyprishuRyFqV+fCIacmaVmOJSqYuz6HaNAHNWamZ3u6t0KfHYzUo1i6FeHOymDtn97n5iUuvu9zCUkf/4hKdHpG2LLbUYhaBZS0qyZmSdKiRtHNQVOtYp+ZpSFXo22e25YuZZzTppSPVn1E636DqnWvEULULKCkvWXLFu8qkaNHj3p/1lpTVFQUvOhEZIlqDlRfg210Ubi/1Shb5X2p3VrQoYOnUpSyx8tFRyECFFDSzsrKCnYcojFwO1FJ/h9SLX1vW47FyMVGIX6GgJL2WWedxQ8//MCOHTs499xzOfvss4MdlwhTtU1lmBuehRr6EujC8s01yMVGIc5EQEn7xRdfZNasWaSlpZGZmcmrr75Kenp6sGMTYaR6I97KEqnKGoV2FKNObEElUX0NdqYJJWBui0f/GCcXG4U4AwEn7S1bttC+fXu2bNnC//t//0+SdhPjXjEVFZ+BpU8ZukzDMTAPLcG9Aiyj52BmfQYtTYz+CnO9xpWpvatHtAPQNiyXrcSQbulCnJGAkrbNZqN9e8/8Y69evXzWbQdDp06dsNvt2O12AB544AGuu+66oL6nqFlFqzDjnFLc6zztwTwJuRRd+Bbm/o89uxhzC1EWz8qRisROHLiWa0gbi+Ws/qE+FSEiXkBJu7S0lG3btnnXaJ96u0eP+r/yv2jRInr27FnvryvqzjOHXea/Ae8OE/O7PKxjDU9zg50KS1fl2SyTBO5MwHYWtvGvh+4EhGhEAkraJ0+e5PLLL/e5r+K2Uopdu3bVf2QBKioq8tnwEx0dXWstlIpjI2UnZ6jj1Y5itOMY5iFQhp8GvGkG5vcm5nETY6BnaqT0W43RviXqqA1aDcF62ytgsYXlZx7qz7euIineSIoVQh9voO+rdBh+op06dSIhIQHTNLngggt44oknSElJ8TmmuLiYhITqNT5nzpzJrFmzGirUxstVhvHe7zEKPke1NNEHjgFlWK9SKItv4nZ9aqLalNfEBo68byXqvPlYug4Hmb8WIiDFxcW0adOGoqIi4uNr/v+mztvYG8LatWvp2LEjTqeTBx98kJtuuon33nvP77E5OTk+JxjISNvhcGC32yNid2co4tWuMlzzL8DonoNlQOX97kww1+tqo21dALoYyPeMtuNIxtJrVEQ0MJDvQ/BEUqwQ+nidTmdAx4Vl0u7YsSPguQB61113kZaWVuOxCQkJtf5WqolSKiK+SBUaMl7X0iko9mDpVr27jHObiVGmKws8ZWmMVmAZYmBmaVwfKMyUwdhiEuTzDaJIijeSYoXQxRvoe4Zd0j5x4gROp5PExETA0wm+T58+oQ2qCdGOYvSBtRgpmsrtMJVUvIFrsYZkE457VpFUjLyNLgr3Fjvm8McbOGohmo6wS9p5eXmMHz8et9uN1prOnTvz+uuy8qCh6MJ9qGR3jd1l9BFAWaDY5WkPVuiZMjEGgrIojDbN4UQ+JMsWdSGCIeySdufOnfn2229DHUaTpRLbwzE7Kr6GnY2A0dfE0rWyf4aZpb1z3bowCuLbNXjcQjQVYZe0RWgpezwqdSjELkEfLPXubDTzARdgxydhg2daxJWpcW+3oVoNlRUjQgRRQO3GRNNiGfcSHL8afaIdOiYacx/gAstQMFrW8KRmBrpgBMa4fzZkqEI0OZK0RTXKGoU1fQGWievgWDyWy0HFgEpUNXdSd7TGOv416aAuRJBJ0hY1cxShWhRjxBqoeE+vx4q57qrM7ChU62FSCEqIBiBz2qJm9gR0vgvAu0XdLAIOaU8nmpRocCRLqVUhGpAk7SaualMDwKfBwdEDu4l2ujG+BeNXCssQA6NMY/6oMbNsGOf/G6OLjLCFaEiStJuoiqYG5sFPUTEO9JEScAIpdiiNQbUawterVzH8bNBF5eVVDQ1RYCQAsSmSsIUIAUnaTZRr2c2QtwilTVQ0YAdtB446UPZCdNY7jLza8BaBAnDvNNEHQLW1YxRLwhYiFCRpN0HaUYzOWoGlr8bo4rtJxr1RYwwH9xqwdPN9nqWrgfN7hdlsLNZ0mcMWIhRk9UgToR3FmLlbPf/O2w5RpT67HcGzSYYo0PvBaOX/dYy2LbAMmSVL+4QIERlpN3L+GvIS0wOV6P941QJoBjqrhtcrivZetBRCNDxJ2o1c1Ya83vu2H8bc4r+Knz4KllSFGa9xZ5o+c9pmtqcDu8xlCxE6krQbMbNwP2bOR1hHl1I1QVu6OzE3adw7NJa0Khcad5ioZDy1slM15lcW9K5UT9PeQpusxxYiDEjSboS8UyIHPkIl5ONeDSpeYwysbBWmOnp2OLp2mNDc030GF6j24FppgjJQPcZjHTPXZ+22ECK0JGk3Qr5TIp6RdNXyqQAUgeUSz8+6SGP+CPqwQpfGgWqO0WY4lnEvoaxRqNY9QnQmQohTSdJuZLSjGJ231mcOGyrLp+oyjblbQyzelmEqRWEMB+d7qVhHvoNK7S6jaiHClCTtRsbM2w7RJ9Fl5UWdjgFx5Qm6ObhWanADrcHcZUIrMGINzOwojDbDMc7qH8rwhRCnIUm7kfDOY+d+io4qxLXEs+VctQaOeua09WFgGPA+kAumC/QmcJco6D4Gm2yYESLsSdJuJKrOY7vXadT5yrdV2E4T9gOrwLjQt/uMe4eJufVr2TAjRASQHZGNQMU8tnF2GbpMo4upttvR0tWA5kB09XZhljQFrhycb12LdvnOhQshwosk7UbAsySvPNkeK9/V6IdKAGq4vmi0Bezv4V4xNRghCiHqiSTtRkAltvdsTweIo+aWYEVAsf+H9FGw/MqFzluLdtRwkBAi5CRpNwLKHo9KGYx7h4GKUn5bgrl3mmgX4PDMYVdlZmlUvGeFiUp0ogv3NWD0Qoi6kAuRjYRGo3eamJtNaAl8Vd4SrCXoI+UHGYANzK/A3GxitPOMsFW8p50Y4NmuLgWhhAhbkrQbAe0ohkNfwEkT1RMsrRTEgXu9hgSwDlHejTRmlsbcodHHQceA5fwqj0lBKCHCniTtRsAs2I0uzgcDrL/0zHjpMg0nwDrcdwbM6KIwM7WnU80BA9dPdow2sVIQSogIIUk7wlRtxFsxIi5Z8yeiWjjBUuXA2laRtABdpFDxo7BcOR+O50tBKCEihCTtCKFdZbjfvdWnmYFucQGqtJioQx9DK2A/uD41PUWhallF4rk/Fev41zyJOja5Ac9ECHEmJGlHCHPFragEz45H7fZU7NPZS+A8MBToYlBngZkPrqUa61WeVSTunabPZhozS4MyMNpcJCNrISKQJO0w42/6A0cxOn8tlvM8G2jM9RpSQBWDOgSqjfKWXLUA7kzTU3tEAzkxmD+ASiz1jLCddlSXcTJ3LUSEkqQdJvz1clSpQzHG/hOK93t3PFZsU7f0ULgPlW9ZvwB0gfZW87N0MzB/MCH1cqKuzyjf5r4dDRhSdlWIiCZJO0z46+VoZmdgrgCGP1m547HiAmMc6EPlz12tPRcXj1bpUJMcjeWih4HyzTdSclWIRkGSdhioqXGB6lCK+4ePoPQYqtVQzOwMVAfPNIeKUmBqjF74Nt8t71BDSTIq6RcNfSpCiCALy23sO3fuZODAgaSlpdG/f3+2bt0a6pCCyqfgE6DdGvc60zOCTjyMdclFHD1SgPvoONwft0WXReP6UaPsvgkbPOuw9REFSRfKNIgQjVBYJu1bb72VqVOnsmPHDmbMmMHNN98c6pCCqqLgky7T6AKNe62G1mAdY2AdBrYxBTRL+ZCczK1YJq7Dcu06ODwSmtfwny8xGsvAGQ17EkKIBhF20yP5+fls2rSJDz/8EIDx48dzxx138NNPP9GpU6dqxxcVFaF1ZXGk6OhooqOja3z9imOrPifUTGcp+oTG9Z7nYiJHwXpR+c7G8uV9UcVwVtttuN8ahEodjnH1q5hv9gPyqr9gSRK07BSScwzHz7c2Em/wRFKsEPp4A33fsEvae/fupW3btlitntCUUnTs2JGcnBy/Sbtjx44+t2fOnMmsWbNO+z4Oh6Ne4j0jrjKMD+7E2LUCo48DS5qBLtCebeblzPXaZ0kf5OHOzsC5yg3JQ1DZy7GcXTm14s6Owp00CKe2QUlJA59QpbD4fOtA4g2eSIoVQhdvSYD/v4Zd0gZPoq6qtt9AOTk5xMdXzt0GMtJ2OBzY7fZq79PQ3Bm3oeOWQWwplrTyqY44MAvAKNDo6PLlfUNO6UJzdhl65xcYv/0KvdqCa1WVZYKthhJ11T9D1josnD7fQEi8wRNJsULo43U6nQEdF3ZJu0OHDuzbtw+Xy4XVakVrzd69e6uNqCskJCT4JO1AKaVC+kUyjx/G3LMU69AyzPLSqdqtMb/S4AT3j55GvKql/+erRCfqxCEsVy/wvyEnxEL9+daVxBs8kRQrhC7eQN8z7C5EtmrVij59+vDmm28CsHjxYjp16uR3aiSSuZfdgkp2+NQIcX+hIRVsVxtYhxpYxyp0of/nV617rezxGK17hE3CFkIET9iNtAH++c9/MmnSJB5//HHi4+N57bXXQh1SvdKOYnTRtyitUVEGKl7j3m6iD4N1aOXvURWlMJI07kyNpVvlb2Gpey1E0xWWSbtbt258+eWXoQ4jaHThPoyWLnB6NsMYAxWu8l2NpzIGKlxLNDorGZWkpO61EE1cWCbtxq5iXbblEoW5XmNmAgboPE9tkYpOMgDKosCiUJe8jhHfOqzmrYUQDS/s5rSbAmWPR6UOxfwp2tObMRaUE1QrcL/v2Q2p3Z4VM2aWBmc0Roe+Mm8thJCRdqiUXvwsn8/8iIGph4kZAJZhlaNrd6aJa5lGRXtqX5udLpdkLYQAZKQdEkVFRYwaewV3LjqCERePJe2UddjdDNAGWqeiWl+HOfrFEEUqhAg3MtJuYAUFBYwdfSl3dN7GVVfGE9WyyO9xqk0ilkGLUB37hXRnoxAivMhIu4EtW7aMe7tu55rrDOyjjsJx0/+BxTGo1O4NG5wQIuzJSLuBnDx5Ert5khtsGeg+JkaJifkZaIefPo5V1mFHSrEdIUTDkKTdALJ3bOPbRy/k8m5uovo5sQ5RVPyR495hYn4PeofpKbV6IgXV9mJZhy2E8EuSdpBt3bqVLY9dwNjrnNj26srCUOUsaQZ6p4kxDNwfJWO96UuMxHYhilYIEe5kTjuIvv32W8aMHMqgNJOYdvjd8Qie+3VWFEbHEZKwhRC1kqQdRPfddx/909rSKi3RpzDUqcyDFrTzCpkSEUKclkyPBEFFWdn//ve/ROlS1OIhqCiFiteeWiNdqmykybJhdByH9Zp/hzBiIUSkkJF2PVu5ciXnnHMOBw8eJDk5mfiUdp4t69lRGAMV+qDGtdLEtVbjzLDDsauxpL8a6rCFEBFCknY9WrRoEenp6aSlpdGyZWX3Asu4l9DF6bhXtwVbK7RqA5bRWG/ZiTV9Qci6zAghIo9Mj9ST119/ncmTJzNx4kQWLFiAzWbzPqasUVjTKzvMWKRSnxDiZ5KRdj3Yt28ft9xyC1OmTOH111/3SdhVSYcZIcSZkpH2GdJa0759ezZs2EDv3r0jqheeECLyyEj7Z9Ja8+c//5lp06ahtaZPnz6SsIUQQSdJ+2fQWjNjxgwefvhh2rZtK8laCNFgZHqkjkzT5I477uAf//gHc+bM4c477wx1SEKIJkSSdh3Nnz+fefPm8fLLLzNlyhSfxypWh0gfRyFEsEjSrqMpU6aQlpbGRRdd5L1Pu8pwr5iKzluLSixDF0Z5O6bLGmwhRH2SOe0AlJSUcN111/Hll19is9l8EjaAe8VUVHwG1lEHsQwowDrqICo+A/eKqSGKWAjRWEnSPo3jx48zZswYVqxYwfHjx6s9rh3F6Ly1GGeX+dxvnF2GzluLdhQ3VKhCiCZAknYtCgsLufTSS9m4cSMffPABl1xySbVjPHPYZX6eDSrRiS7cF+wwhRBNiMxp12LixIls376djz76iH79+vk9RiW2Rxf6n7fWhTZUYvtghiiEaGIkadfiySefxDAMzjnnnBqPUfb48ip+GT5TJFX7PAohRH2R6ZFT7Nmzh5tvvhmHw0Hv3r1rTdgVKqr4uVa1wb0hGdeqNujidGlqIISodzLSrmLHjh2MHDkSm83G4cOHad8+sKmNU6v4yTptIUSwyEi73A8//MDQoUOJjY1l3bp1ASfsqqSKnxAi2CRpA7m5uQwbNow2bdrw2Wef0bZt21qP145izNytspxPCNHgZHoEaN26NU888QTXXnstiYmJNR4nOx+FEKHWpJP2mjVrOHjwIDfccANTp55+92LFzkdLn6qrRDJwrwBr+oIgRiqEEB5hNT0yadIk2rdvT+/evenduzfTp08P2nutWLGCMWPGsHDhQrTWp53ykJ2PQohwEHYj7fvvv5877rgjqO+xaNEibr75Zq688kr+8/oC3Esnn3bKI5Cdj6p1j6DGLYQQYZe066qoqAittfd2dHQ00dHRNR6/ePFiJk+ezG9+8xteeeUV1IpbUAl+pjyWgyX91conJrSrdecjCe184qgvFa8ZjNcOBok3uCIp3kiKFUIfb6DvG3ZJe/bs2bz00kt07NiRRx99lN69e9d6fMeOHX1uz5w5k1mzZtV4/Pnnn8+sWbOYMWMGzmMFWPM+w3Ze9SkP5/ufUXY0D7zL92wYSYNQ2cuxVJkicWdH4U4ahFPboKSkTudaFw6HI2ivHQwSb3BFUryRFCuELt6SAPOH0g34a2XIkCFs27bN72PffvsthmHQpk0bDMMgIyOD2267jZ07dxIbG1vt+OLiYhISEsjJySE+vnJd9OlG2lprHA4Hdrsd8rZhfnYZlgEF1Y5zb0jCGPaBz5SHdpVhrrgVnV9lKqXVUIxx/wza6pGq8UZCWzOJN7giKd5IihVCH29xcTGJiYkUFRX55LRTNehIe926dQEfm56ezv33309mZiZ9+/at8biEhIRaT7AmSilo0aGWKY8oVIsOPv/xlC0a4+rQ7HxUSkXEF7+CxBtckRRvJMUKoYs30PcMq9Uj+/ZVljHdsGEDBQUFdOnSJWjvV1nsyTdxn67Yk+x8FEKESljNaU+aNIm8vDwsFgsxMTEsXLiQhISEoL6nZdxLuFeAa9Xa8lUgNu/qESGECDdhlbTXrFnT4O8pxZ6EEJEkrJJ2KCl7vKyzFkKEvbCa0xZCCFE7SdpCCBFBJGkLIUQEaXJJu7S0lMcee4zS0tJQhxIQiTe4JN7giaRYIXLibdAdkfWpYkfk6XYPnaqoqIjExEQKCwuDvpywPki8wSXxBk8kxQqhjzfQnNbkRtpCCBHJJGkLIUQEidh12hWzOsXFdWs+UHF8cXFxRNRDkHiDS+INnkiKFUIfb8X7n27GOmLntPft20eHDh1CHYYQQtSrvXv30r59+xofj9ikbZomBw4cIC4uLiJ+iwshRG201hw7doy2bdtiGDXPXEds0hZCiKZILkQKIUQEkaQthBARpMkm7UmTJtG+fXt69+5N7969mT59eqhDqmbnzp0MHDiQtLQ0+vfvz9atW0MdUq06depE9+7dvZ/pf//731CH5GPatGl06tQJpRQ//PCD9/78/HxGjRpF165d6dmzJ59//nkIo6xUU7zDhw+nc+fO3s/52WefDWGUHg6Hg6uuuoq0tDR69+7NqFGj+Omnn4Dw/HxrizccP18fuom66aab9N///vdQh1Griy66SL/66qtaa60XLlyoBwwYENqATuOss87SW7ZsCXUYNfrss8/03r17q8U5efJk/dBDD2mttf7f//6nO3bsqJ1OZ4iirFRTvMOGDdMrVqwIYWTVlZSU6JUrV2rTNLXWWv/973/Xl1xyidY6PD/f2uINx8+3qiY70g53+fn5bNq0ieuvvx6A8ePHs3v3bu9oQNTd0KFD/S6leuedd7j99tsB6NevH6mpqWExGqwp3nBkt9u5/PLLvSu5BgwYwK5du4Dw/HxrizfcNemkPXv2bM455xzGjh3Ld999F+pwfOzdu5e2bdtitXr2Pyml6NixIzk5OSGOrHa//e1v6dWrF7/73e84dOhQqMM5rYKCAkzTJCUlxXtfp06dwv5znj59Or169eK6664Ly2QzZ84cxo0bFzGfb0W8FcL58220SXvIkCEkJyf7/Wfv3r089thjZGVl8f3333PzzTczevRojh8/HuqwfZy6/lyH+erMtWvXsnnzZjZt2kRSUhI33XRTqEMKSKR9zm+88Qbbtm3j+++/Z8iQIYwdOzbUIfl4/PHH2blzJ4899hgQ/p/vqfGG++fbZOe0T5WWlqY3btwY6jC88vLydHx8vHfuzzRNnZqaqnfv3h3awAJ04MABHRsbG+ow/Dp1jrhZs2Y6Pz/fe7tfv376k08+CUFk/p3uWkF0dLQ+fPhwA0ZUs7/+9a+6b9+++ujRo977wvnz9RfvqcLp89W6Cc9p79u3z/vzhg0bKCgooEuXLiGMyFerVq3o06cPb775JgCLFy+mU6dOdOrUKbSB1eDEiRMUFhZ6b7/11lv06dMndAHVwTXXXMPcuXMB+Prrr8nNzWXw4MEhjso/l8tFXl6e9/bixYtJTU0lKSkphFF5zJ49m7feeovVq1eTmJjovT9cP19/8Ybz51uhye6IHDlyJHl5eVgsFmJiYnj88ce56KKLQh2Wj8zMTCZNmkRBQQHx8fG89tpr/OpXvwp1WH7t2rWL8ePH43a70VrTuXNnnn/++bD6JXP77bezbNkycnNzSU5OJjY2lqysLPLy8rjhhhvYvXs3UVFRvPjiiwwbNizU4fqNd/PmzQwbNozS0lIMwyA5OZnZs2dz7rnnhjTWilpAnTt3Ji4uDoDo6Gi++uqrsPx8a4r3448/DsvPt6omm7SFECISNdnpESGEiESStIUQIoJI0hZCiAgiSVsIISKIJG0hhIggkrSFECKCSNIWIVG1jGuPHj28my9++uknrFartyxm7969GThwoPd5+fn5TJ48mc6dO9OrVy969erF448/DsD7779Pz5496dmzJx988IH3Oa+88gpPPPFEjbF8+umnNGvWjN69e3POOecwePBgvv/++1rjLyws5Omnn/7Z5z9p0iReeOGFn/18gFmzZtGrV6+wLYUrgiSk+zFFk1V1a3ZOTo5OSEjQmzdv1rt379ZJSUl+n3Py5EndvXt3/dBDD2mXy6W11vr48eP6ueee01pr3bdvX71nzx69Z88e3bdvX6211rm5uXr48OG6rKysxlg++eQT7/Faa/3cc8/p8847r9b4a4tTa33a0qP1URq46tbr/fv367i4OH3kyJEzek0R/mSkLUKuQ4cOpKWlsWPHjlqP+89//kNcXBwPP/wwFosFgObNm/OHP/wBAJvNxsmTJzlx4gRRUVEA3H333Tz99NPYbLaA47nkkkvIzMwEPNuuL774Ys4//3zOO+88Fi9eDMDvf/97CgsL6d27N+effz7gKZ4/a9YsRowYwWWXXYbb7ea+++7zjv7vvPNOysrKqr3f8ePHmTJlive4P//5z97Htm7dygUXXEDPnj35zW9+w4ABA3j33XcBfLaKHzt2DKUUpmkGfJ4iMllDHYAQW7ZsYfv27d6twhXJsMI555zD66+/zjfffMOFF15Y4+s8/fTT3sqCzz77LO+++y6tW7emX79+dYrn7bffpm/fvhQWFnLrrbeycuVK2rRpw+HDh+nbty+DBg1i3rx5nH/++dVK+n733XesWrUKm83GP/7xD7755hu++eYbLBYLV1xxBc8//3y1LkmPPPIIZWVlfP/995SUlDB48GB69OjBNddcww033MDdd9/N9ddfzzfffEP//v19njtnzhzmzp3Lvn37eOWVV8KqRoYIDknaImQmTJiA3W6nWbNmvPLKK3Tt2pWffvqJxMTEn1XffMiQIXz11VeAZ+R5+eWXs2rVKp577jm++OILUlNTmT17tncUXtXWrVu9vyjS0tJ47bXXWL9+Pbt27WL06NHe47TWZGZmctZZZ/mN4YYbbvCO6tesWcPNN99MdHQ0ALfccgvz5s2rlrTXrFnD888/j2EYNG/enBtvvJE1a9Zw2WWX8cMPP/Cb3/wGgL59+3LOOef4PHfatGlMmzaNzZs3c/311zNy5EhJ3I2cJG0RMosWLaJnz54BH9+3b19eeumlgI6dOXMms2bNIjc3l2XLlvHJJ5/w0EMP8e9//5vJkydXO75Hjx5s3LjR574ff/yRc845h7Vr11Y7vqYOQrGxsd6ftdbVakmferu24yru9/ecU5177rm0a9eOTz/9lPHjx5/2eBG5ZE5bRIxf//rXFBYW8sgjj+B2uwE4efIkTz75pM9xGzZsoKioiFGjRnHixAlv0jMMo06NLgYOHMjOnTv5+OOPvfd99913lJWVER8fz8mTJ3G5XDU+/5JLLmHBggWUlZXhcrl4+eWXGTlypN/j5s+fj9aaEydO8OabbzJy5EgSEhLo0aMHb731FgDffvstW7Zs8T5v27Zt3p+zs7P59ttv6dGjR8DnJyKTjLRF2Dl1Thvgyy+/pFmzZnz22Wfcf//9dOnShdjYWJRS3ukDAKfTyf3338/ChQsBz3x4586d6dmzJ8nJyWRkZAQcR4sWLVixYgXTp0/n7rvvxul00rFjR5YuXUrLli29rdWaN29ebZQOMHXqVLKzsznvvPMAz4XKadOmVTvuj3/8I3feeSe9evUCPPWnJ0yYAMDrr7/O5MmTeeaZZ+jTpw/nnnsuCQkJANx///1kZWVhs9mwWq288MIL/PKXvwz4/ERkktKsQoSxEydO0KxZM5RSbN26leHDh5OZmUmLFi1CHZoIERlpCxHGvvjiC6ZPn+7tqzh//nxJ2E2cjLSFECKCyIVIIYSIIJK0hRAigkjSFkKICCJJWwghIogkbSGEiCCStIUQIoJI0hZCiAgiSVsIISKIJG0hhIgg/x+mFwHP/0atOgAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "gp.layout(colors=gp.colors.bright)\n", "limits = (-5, 27)\n", "\n", "fig, ax = plt.subplots(figsize=(90*mm, 85*mm))\n", "\n", "ax.plot(limits, limits, \"--\", c=\"k\")\n", "ax.plot(- petrolog.loc[pec.index, \"OL_PER\"], pec[\"total_crystallisation\"], \"o\")\n", "\n", "ax.set_xlabel(\"PEC% Petrolog3\")\n", "ax.set_ylabel(\"PEC% MagmaPEC\")\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "0082993b", "metadata": {}, "source": [ "With these similar settings results plot very close to the 1:1 line." ] } ], "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 }