{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# PEC Monte Carlo simulation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Start by importing MagmaPEC and MagmaPandas and any other packages you want to use. Here we also import Pandas for importing pressure data. For details on the use of MagmaPandas, please see it's [documentation](https://magmapandas.readthedocs.io/en/latest/)." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import MagmaPEC as mpc\n", "import MagmaPandas as mp\n", "\n", "import pandas as pd" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Confirm the model and PEC configurations. If you want to change models or PEC settings, follow the [configuration example](https://magmapec.readthedocs.io/en/latest/notebooks/config.html)." ] }, { "cell_type": "code", "execution_count": 2, "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+.............................sun2024\n", "Kd Fe-Mg ol-melt........................toplis2005\n", "Melt thermometer....................putirka2008_15\n", "Volatile solubility model.......iaconomarziano2012\n", "Volatile species.............................mixed\n", "##################################################\n", "\n", "\n", "############ Post-entrapment crystallisation ############\n", "################### correction model ####################\n", "Settings_________________________________________________\n", "Fe2+ behaviour...................................buffered\n", "Stepsize equilibration (moles)...................0.002 \n", "Stepsize crystallisation (moles).................0.05 \n", "Decrease factor..................................5 \n", "FeO convergence (wt. %)..........................0.05 \n", "Kd convergence...................................0.005 \n", "#########################################################\n" ] } ], "source": [ "print(mpc.model_configuration)\n", "print(mpc.PEC_configuration)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the next few steps we import all relevant data and set up the melt initial FeO prediction model. These steps are identical 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": "markdown", "metadata": {}, "source": [ "Import melt inclusion and olivine data:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "melt_file = \"./data/melt.csv\"\n", "olivine_file = \"./data/olivine.csv\"\n", "\n", "melt = mp.read_melt(melt_file, index_col=[\"name\"])\n", "olivine = mp.read_olivine(olivine_file, index_col=[\"name\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Import inclusion internal pressures or calculate them if you have measured melt CO2 (and H2O). See the [PEC model example](https://magmapec.readthedocs.io/en/latest/notebooks/pec_corr.html) for details on how to do the calculation. Here we import them from a file." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "pressure_file =\"./data/pressure.csv\"\n", "pressure = pd.read_csv(pressure_file, index_col = [\"name\"]).squeeze()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set up the melt initial FeO prediction model:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAAFaCAYAAADVfgw3AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAARmpJREFUeJzt3XlYE+f6N/BvWATcQFEgskhldUEQsIIiIC5V3EBxx13rEa3SVi0/rbX0VKvt61Fr3etua10qoiCuFRS3ulFLPVZREdG4YRsQ2fO8f3BICRBIgGRmwv25Lq/LhJnMfSf4zfjMzDMixhgDIYQQQdDjugBCCCGqo9AmhBABodAmhBABodAmhBABodAmhBABodAmhBABodAmhBABodAmhBABMeC6AE2RyWR4+vQpmjVrBpFIxHU5hBAixxhDTk4O2rRpAz099faddTa0nz59CltbW67LIIQQpR4/fgwbGxu11tHZ0G7WrBmA0jelefPmHFejPsYY8vLyYGJiItj/KVAP/CD0HoReP1C5h+zsbNja2spzSh06G9plH27z5s0FG9qGhoaC/0WlHrgn9B6EXj+gvIfa9KOzoV2VoqIilJSUcF2GShhjKCgogEgkEvQvasUe9PX1YWhoyHFlhAhXgwntoqIi3L9/H0Ka1FAmk6l9kIJvKvYgEong4OBAwU1ILTWY0C4pKQFjDG3atIGRkRHX5dSIMSYPPCHvaZfvoaCgAE+fPkVJSQmFNiG11GBCu4yRkRGMjY25LqNGuhjahJC6E/b/vYlSM2bM4LoEQogGUGhriUwm08o6ZTZt2lTrdQkh/EWhrUGJiYkYNGgQhgwZgt27d6Nnz57o3r079u7dCwC4efMmvL29MWTIEAwdOhSJiYnydYYOHYo9e/bg+PHjldb79NNP0b17d/j7++Py5cu4dOkSunXrhoCAAHz22WcAAG9vbwCl56kHBQWhZ8+emDlzJgBgx44dGDZsGAYOHIiuXbvi6dOnHLw7hJDaaHBj2uVJJBJIJBKF51q0aIF33nkH+fn5uH37dqV1PD09AQB//vknmjdvDrFYXO02srOzcfbsWfTs2ROJiYkwMDBAr169MHLkSHz66af48ccf4eTkhICAAIV1EhMTUVJSgoCAgErrnThxApcuXYKBgQFkMhmWLFmCxYsXY9CgQZX2zpcvX44FCxagf//+mDp1KpKSkgAAZmZm2LZtG7Zs2YKDBw9izpw5tXoPCWnoJBIJoiKnY/nqLTXmQX1o0HvamzZtgpeXl8KfxYsXAwAyMzMr/czLy0u+7qRJk1QagvD29kZWVhbu3buHfv36ISgoCK9evcLLly/x4sULODs7QyQSKby2t7c3RCIRXr16VeV6X375JWbMmIEZM2bgxYsXmDVrFk6dOoUJEybg+PHjCtu/f/8+unbtCgDo2rUr0tLSAABdunQBANja2uKvv/6q2xtJSAMlkUgQHhqAMPExhIcGVNoJ1IQGvac9Y8YMDBkyROG5Fi1aAABsbGxw/fp1pevu2LFDpSst9fT00KpVK7Rv3x6nTp2CoaEhioqKYGhoCEtLS9y7dw+Ojo64ceMGhg4dKl8HgNL1TE1N8d577+HHH3/E5s2bMX/+fKxZswaFhYXw8vJCcHCwfPuOjo64evUq+vfvj6tXr2LixIl48OCBwtkcQjp3nRC+KAvsdUPT4GwlgpNlGsJDA7AnJkmje9wNOrTFYrHSN9fY2Fg+FFIVFxcXlbejp6eHRYsWoU+fPtDT00Pr1q2xf/9+fPHFFxgzZgysrKzQtGlTeTCXX2/hwoWV1gsJCUFeXh4KCgrw/fffY9OmTTh06BByc3MxadIkhW1/8sknmDhxIpYuXYpOnTrB398fDx48ULl2QkjVoiKnI9KnNLABwNlKhEifNERFTsfOfXEa266I6ehuVnZ2NkxNTSGVStG8eXPk5+fj4cOHeOedd3hznnbZnrNMJkNQUBD27t0r/xLRhXOcK/bAx8+gJro4WZHQ8LX+invad58xzIp1rHJPu6oJo8rnkzoa9Jg2165cuYKAgAB069YNffr00cpBDEJI/RCLxdgTk4RZsY6I/015YNe3Bj08wjU/Pz/52RyEEOEpC+6oyOnYE6Ods0cotAkhpA7EYrFGx7ArouERQggREAptQggREAptQggREAptJSQSCSaOGqTxK5w+//xzxMXF4dmzZ1iyZAkAIDAwEG/evFH7tXbs2IHCwkL53y9dulQvNfr4+Cj9WWJiIu7evVsv2yGE1IxCuwpcXJpqZWWF6OholZZVNvtf+dCeNGkSfH19660+ZSi0CdEuCu0Kyp8wP9BdhHVD02oV3IwxzJo1Cz179kRAQABev36N3bt3IygoCJ6enti9e7fC8unp6QgLC5M/jo6ORq9evfDBBx8AKA3kUaNGYeDAgTh9+jQ+/vhjBAYG4t1330VKSgouXbqElJQUDBgwAGvWrJHvwQPAhx9+CD8/P/Tq1QsPHz4EALRv3x7jxo1Dly5dKtUCAMuWLYOvry9mz54tv69mxfrz8vKwY8cO/N///R8mT56M58+fo3fv3vD390dYWJhg7sdJiKAwHSWVShkAJpVKGWOM5eXlsdu3b7O8vLxq15swciA7MlfESrbpyf8cmStiE0YOVGv7hw8fZnPmzJE/LikpYbm5uYwxxvLz85m3tzdjjLElS5awo0ePsocPH7Lhw4czxhgLCAhgx48fZ8XFxWzs2LHs2rVrbPv27Sw8PFz+emWv9dtvv7GxY8fK18vJyVF43V9//ZWNGjWKMcbYuXPn2KRJkxhjjJmZmTGpVMpycnKYj4+PQu0SiYT5+/szmUzG/vzzT/bOO+8obLOq+hljrKCggBUVFTHGGPvwww/ZiRMnWHFxMZPJZIwx1T8DPpHJZCw3N1fegxAJvQeh189Y5R4q5pM66DztCpav3oLw0AA4Wf5zaerqy47YE7NFrde5c+cOevbsKX+sp6eHU6dO4T//+Q8A1DikUDbrX/mZ+cpm6wOAlStX4sSJE9DT04O+vr7S16k4y9/ChQsBAO3atZNfPssqzGSQnp6Ozp07QyQSwdnZGaampgBQY/2vX7/Gv/71L/z111+QSCTo3LlztT0SQtRHwyMV1Nelqe3bt8eFCxfkjxljiI6OxpEjR3DixAk0bdq02vVv3rwJALh+/TocHR0B/DP7X1ZWFuLi4nD+/Hl899138tA1NDSsNCRRNssfAFy9ehVOTk4AUO0cDvb29khNTQVjDGlpaZBKpQBQZf3lt/nDDz+gX79+SEpKwqBBg2j2QEI0gPa0q1Afl6YOHjwYCQkJ8PPzQ6NGjXDgwAEMHz4cvXr1goeHh3wKWGUSEhLwxRdfwN3dHV5eXvj999/lP2vRogUsLS3Rq1cvhTM7hgwZgpEjR2LkyJHy57y9vSEWi+Hn5wcDAwNs3769xtqtrKzQt29f+Pr6wtPTE+bm5gBQZf1BQUH45JNPkJSUhAkTJmD8+PE4ceIEGjduDDc3N7XeM0JIzbQ6y9+cOXNw5MgRPHr0CL///js6deoEAJgyZQouXLgAExMTNG/eHN9++y08PDyqfI24uDjMmzcPxcXFcHd3x86dO6vcaxXCLH/VYTTLHy8wns4wpw6h9yD0+gEBz/IXFhaG5ORktG3bVuH5kJAQ/PHHH0hJScGCBQsU9hTLe/PmDaZOnYrDhw8jLS0NYrEYS5cu1UbphBDCC1odHvH396/y+fJ3j/Hx8cGjR4/ke2jlJSQkwNvbG66urgCAiIgIBAcH46uvvlK6TalUCsYYCgsLIZPJIJPJBDfWKrR6q8IYk/dR/u98V75moRJ6D0KvH6jcQ1164d2Y9po1axAcHFwpsAEgIyNDYS/d3t4eT548qTLgy9jZ2QEA2rZti40bN6Jp06awsbHRTPEaoOxCGiEp66HsSzM/P19w/wDz8/O5LqHOhN6D0OsH/ukhLy+v1q/Bq9Des2cP9u/fj/PnzytdRt0xrYyMDDRv3hyFhYV48eKF/O9CUN2XkVCU76GoqAh6enowNjYW1Jh2fn4+jI2NBT2eKuQe+Fy/VCqVnxJbnYo9lL+toLp4E9r79u1DdHQ0zpw5AwsLiyqXsbOzwy+//CJ/nJ6eDmtr62qDzdTUFM2bN0dRURGysrLw7Nmzeq9dU3QttIHSL10DAwPe/eOriUgkElzNFQm9B77Vv2XLFixcuBDXrl2rdJxOmbIe6tIHL0J7//79+PTTT3H69Gn5cEZV+vfvj1mzZuHOnTtwdXXF+vXrMXr0aJW2YWhoCAcHB8FcWs3nvQtVVdWDvr4+DA0NOa6MkLpZsWIFoqKi8MEHH8DW1la7G1f7Gso6iIiIYNbW1kxfX59ZWloyBwcHxhhjBgYGzMbGhrm7u8v/vHr1ijHG2OLFi9mGDRvkrxEbG8tcXFyYg4MDCwkJUXoZaF0uE+UDXbx0V4ioB+7xqX6ZTMbmz5/PALAlS5aoXFN9XsbeYO7GLjRMB89NFSLqgXt8qj89PR0eHh6Ijo7G3LlzVV6vYg91ySdeDI8QQgifFRQUgDEGe3t7pKWloVWrVpzVIuyjXIQQomFv3rzBoEGDMGXKFADgNLAB2tMmhBClsrKyMHDgQNy+fRtHjhzhuhwAFNqEEFKlp0+fol+/fnj+/DnOnj0rny6ZaxTahBBShd27d0MqleL8+fPyqTP4gMa0CSGknJycHADAggULcOPGDV4FNkChTQghchcuXEC7du1w+vRpiEQitG7dmuuSKqHQJoQQAMePH0ffvn3RoUMHvPvuu1yXoxSFNiGkwdu3bx+GDBmC3r174/jx47y+II9CmxDSoBUWFuLzzz/HqFGjcOjQIZiYmHBdUrXo7BFCSIPEGENubi6aNm2Kc+fOwdzcXBCzavK/QkIIqWeMMSxYsAC+vr7Iz89H69atBRHYAO1pE0IamOLiYsyYMQPbtm3DmjVrBHNDjjIU2oSQBqOgoABjx45FbGwsdu3ahfHjx3NdktootAkhDca5c+eQkJCAmJgYDB48mOtyaoVCmxCi83Jzc9G4cWP07dsXDx48gJWVFdcl1ZowRt4JIaSWnjx5gnfffRf/+c9/AEDQgQ3QnjYhRIfdu3cPffv2hUwmw6BBg7gup17QnjYhRCelpKTAz88PJiYmuHDhAlxcXLguqV5QaBNCdNJXX30FW1tbnDt3Tvt3TNcgGh4hhOiU3NxcNGnSBNu2bUNJSQmv5xGpDdrTJoTojL1798LR0REPHz5EkyZNdC6wAQptQoiO2LBhA8aNG4d+/frp1HBIRRTahBBBY4xh6dKliIiIwJw5c7B9+3YYGOjuyC+FNiFE0DIzM7FixQp88cUXWLVqlWAmfqot3f06IoTotOLiYpSUlMDW1hZ37txBmzZtuC5JK3T7K4kQopPy8/MxYsQITJo0CQAaTGADtKdNCBGYnJwchISE4OLFizh48CDX5WgdhTYhRDCysrIwYMAA/Pnnnzh58iR69uzJdUlaR6FNCBGM7du3Iz09HYmJiejSpQvX5XCCxrQJIbyXl5cHAPj444+RkpLSYAMboNAmhPBcSkoKOnfujJMnT0IkEjWog45VodAmhPDW+fPnERgYCCsrK3h6enJdDi9QaBNCeOnYsWPo168fvL29cezYMbRq1YrrkniBQpsQwjvFxcWYN28e+vfvj7i4ODRr1ozrkniDzh4hhPBKfn4+jI2NcebMGbRu3Rr6+vryA5FEy3vac+bMgb29PUQiEVJTU+XPL1u2DC4uLtDT00NcXJzS9dPT02FgYAAPDw/5n/v372ujdEKIhjHG8OWXX8LHxwdv376FWCzW6YmfakuroR0WFobk5GS0bdtW4fnevXvj2LFj8Pf3r/E1zMzMkJKSIv/j4OCgqXIJIVoik8nw0UcfYfHixRgxYgRMTEy4Lom3tPo1piyUu3Xrps0yCCE8UlxcjGnTpmHXrl1Yv349Zs6cyXVJvCa4/3tkZ2eja9euKCkpQUhICBYtWgR9fX2ly0ulUjDG5I+NjIxgZGSkjVLrpKzm8rULDfXAD3zvISkpCXv37sWePXswZsyYSnXyvX5VVOyhLr0IKrTFYjEyMzNhYWGB169fY9SoUVi5ciUWLFigdB07OzuFxwsXLsSiRYs0XWq9yc/P57qEOqMe+IFvPeTn58PIyAi+vr64desWbG1tqz3gyLf6a6Osh7ocWBVUaBsZGcHCwgIA0LJlS0yZMgU//vhjtaGdkZGhcJ84Ie1plx1FF4lEXJdTK9QDP/Cxh5cvXyI4OBijRo3CvHnz4OzsrHRZPtavroo9FBUV1fq1BBXaL168QIsWLWBoaIiCggIcOnSoxjkITE1NBX1zT5FIJNhf1DLUAz/wpYfHjx+jX79+eP36Nfr166dyTXypvy7KeqhLH1o9e2TWrFmwsbFBZmYm+vTpA0dHRwDAV199BRsbG1y6dAmTJk2CjY0NXr58CQD47LPPsHHjRgBAcnIyunTpAnd3d3h6esLKykpQQx2ENHR3796Fn58f8vLykJycDHd3d65LEhwRE/LofjWys7NhamoKqVQqyD1txhjy8vJgYmIi2L0L6oEf+NTD2LFj8dtvv+HkyZOwtrZWaR0+1V9bFXuoSz4JaniEECJMBQUFMDIywubNm1FQUABzc3OuSxIsmnuEEKJRcXFxcHJywv3799G0aVMK7Dqi0CaEaMwPP/yAkJAQeHt7qzwcQqpHoU0I0Yi1a9ciPDwcEyZMwP79+2FsbMx1STqBQpsQUu8kEgkWLlyIjz/+GFu3bqWJn+oRhTYhRE4ikWDiqEGQSCS1Wl8mk6GoqAhisRi///47vvnmG8Ge8cFXFNqEEAClgR0eGoAw8TGEhwaoHdxFRUWYNGkSxo8fD8aYfBpmUr8otAkh8sBeNzQNA91FWDc0Ta3gzsvLw/Dhw7F3716EhIRQWGsQhTYhBFGR0xHpkwZnq9KwdbYSIdInDVGR02tcNzs7GwMGDMDp06dx5MgRjB49WtPlNmgU2oQQLF+9BasvO+Lus9ILpO8+Y1h92RHLV2+pcd2tW7ciJSUFp06dwoABAzRdaoNHoU0IgVgsxp6YJMyKdUT8bwyzYh2xJyYJYrFY6TqFhYUAgLlz5yIlJQU9evTQVrkNGoU2IQTAP8F9UBJcY2D/+eefaN++PU6ePAk9PT3Y29trr9AGjkKbECInFouxc19ctYF9/fp1+Pn5wcTEBB07dtRidQSg0CaEqCExMRG9evWCg4MDzp07R5emc4BCmxCikpKSEsyaNQvdunXD6dOn0bJlS65LapDo2lJCSI0KCwvRqFEjnDhxAq1btxbELft0Fe1pE0KqtWbNGvj6+uLNmzewsbGhwOYYhTYhPFHXeT/qG2MMn3/+OSIjI9G7d280adKE65IIKLQJ4YW6zvtR32QyGebOnYvo6GgsX74cX3/9NV2azhM0pk0Ix8rP++FsJYKTZem8HzWdK61JycnJWL9+PTZt2oT333+fkxpI1WhPmxCO1WXej/pWVFQEAPD398d///tfCmweotAmhGN1mfejPkmlUvTp0wcrV64EADg5OWl1+0Q1FNqEaFnFA461mfejvr148QK9evXCrVu34Ovrq7XtEvVRaBOiRcoOOKoz70d91jJx1CBcvXoVPXv2hEQiQVJSErp3767xbZPaowORhGhJTQccy+b90GYtkT5pGDP4OPL0LZGcnAwHBwetbJ/UHu1pE6IlfDngWPEuNXFzZWjKXiAnJ0erdZDaodAmREv4csCxqi+P/zeyBCMH+XF+fjipGYU2IVrChwOOANBrQBjm/MAUvzxOMiwakMvJaYZEPRTahGgRFwccy9u1axemTZsGG1c/DFzFSr88djMsGQrs+t2p0l4/3y6tJxTahGidKjca0ISXL19i1qxZGDVqFAwLn+H/jQI+/olhQg9g+i5DrNy4X6Emvl1aT0pRaBOi4xhjKC4uRuvWrXHt2jXoF0vxoe99DPXUw9lPRDhzG1gYXIRVX30qX6fiwcp1Q9MouHmCQpsQHSaTyfDBBx9g2rRpYIzBxcUFK8odEBWbibBwkKjS0AhfznQhlVFoE6KjioqKMH78eGzcuBEBAQHyWfpUOSDKlzNdSGUU2oTooLdv3yI0NBQHDhzATz/9hMmTJyv8vKYDonw504VURqFNiA76/vvvkZiYiPj4eISFhVW5TE0HRLk+04VUTcQYY1wXoQnZ2dkwNTWFVCpF8+bNuS5HbYwx5OXlwcTERLCTz1MP2ldcXAwDAwPIZDLcu3cPLi4uguuhIqHXD1TuoS75pNU97Tlz5sDe3h4ikQipqany55ctWwYXFxfo6ekhLq76uRfi4uLg6uoKR0dHDB8+HG/evNF02YQIwqNHj+Dh4YGTJ09CT08PLi4uXJdENKBWoX337l388ssvuHTpklrzFYSFhSE5ORlt27ZVeL537944duwY/P39q13/zZs3mDp1Kg4fPoy0tDSIxWIsXbq0Ni0QolNu376NHj16IC8vD46OjlyXQzRI5dDOyclBdHQ0bG1tERwcjMWLF2P27Nmws7PDgAED8Msvv9T4Gv7+/rCxsan0fLdu3VSaXSwhIQHe3t5wdXUFAERERGDv3r2qtkCITrp69Sr8/f3RsmVLJCcno127dnV+TboSkr9Unpq1V69eCA8Px9WrV2FlZSV/XiaT4fz589i4cSPS0tI0enuijIwMhb10e3t7PHnyBDKZDHp6VX//SKVSlB+2NzIygpGRkcZqrC9lNQv5kAP1oHkymQxTp06Fs7Mz4uLi0KJFi0q1qtuDRCLB+GEBiPS5j/BQf+w+xO1BSL5/Bqqo2ENdelE5tC9cuFBl2Onp6SEgIAABAQEoKCiodSGqUvdAhJ2dncLjhQsXYtGiRfVZkkbl5+dzXUKdUQ+aUVJSAn19fezfvx/m5uYwNjZGXl6e0uVV6UEikWDamPewIfTB/+b8vo9xIf74fu8Jzs8e4eNnoK6yHqr7nGqicmhXFdjp6enIzc1Fx44dlS5Tn+zs7BSGYdLT02Ftba10Lxso3Tsvf3RWSHva+fn5MDY2FvQRc+qhZhKJBFEfvo/lqzarFIwSiQTjRgzA32/1kJSUVOMBx+p6qLjt6IVz8VH3BwpXQn7U/QGiF87Fzp+O1r7JOtDF36OyGyjXRq3vXLN27VrExMRAJBLB2dkZGzZsqHURqurfvz9mzZqFO3fuwNXVFevXr8fo0aOrXcfU1FSQp/yVEYlEgv1FLUM9KFc6FBGISJ80jB8WWOP50BKJBIN6uSN6QBYWxjZDdna2yr/fFXuoatsrVm9BeGgAnCxLL2EvuxJyT8wWzj9DXfo9qlMfTEUXL15UeDxy5Ej5393c3FR6jYiICGZtbc309fWZpaUlc3BwYIwxtmzZMmZtbc0aNWrEzM3NmbW1NXvx4gVjjLHFixezDRs2yF8jNjaWubi4MAcHBxYSEsKkUmmV25JKpQyA0p/znUwmY7m5uUwmk3FdSq1RD9V7+vQpC+rmxP67TMRKtumx/y4TsaBuTuzp06dVLv/kyRPm7thS5eWr66G6bZf97Mhc1V5f03Tx96gu+aTyxTWTJ09G06ZNsXz5cjRp0gSzZ8+GhYUFRCIRzpw5g8TExNp/c2gAXVzDPepBOYlEgv4BHvgo4AXGd/9neC/+N4aDkuAq7xU5uF8PvN/hEga6i1RavroeJo4ahDDxMaWvJZFIEBU5HctXb+F8LFsXf4+0cnHN9u3bMXjwYAwYMABHjx7FqlWrYGNjAzMzM8TExKjdBCENVdm0p18Gv8SyY4ZIvisDoHxSppKSEgDA5p0HsSLJVqVJnGo6Za+mCaG4mvObqEDdXfO3b9+yefPmsTFjxrDnz5+rvWuvLTQ8wj3qobKqhiWc2xiyHdNQ5VBEbm4uCw4OZitXrlRYv7qhi4rL3Lhxg40L68+ePHlS7XJcD4Moo4u/R3XJJ5VDOyMjg82fP58tWrSIvXjxgl27do0FBQWxrVu3qr1RbaDQ5h71UNmEkQPZkbmlgV3258hcEevsZFEpNP/66y/m5+fHmjRpwk6dOiV//unTp2zCyIHVBnbVXwqOldap7rX4Qhd/j+qSTyoPj4waNQpubm6wsbHBxIkT4eXlhZMnT+LVq1cYOHCgpv4jQIha+H4ln7JhieNJKQpDEc+fP0dgYCD++OMPnDlzBn369JH/rLqhiyrvtD6iGGduA+uG3q909xkaBhEgVdO9Y8eOjDHG8vPzmYeHh8LP7t+/r/a3habRnjb3tN2DJv67r4keVKlz8uTJrE2bNiw1NbVWr61wVkh7sMz/iOR79RNGDqyvVrRCF/8taG1P29XVFW5ubpgxY4bCz+pjrgNC6kJI9zSsbp5qmaz0oOTq1atx4cIF+YVr6r522c0Lhn5ngCVDAbGZiO4+oyvUSfjs7GyWm5ur9jcDF2hPm3va7EHZWHFd9yrrowdVx42vXLnCOnTowNLS0mq9rYrbvHnzJgvq5vi/vfrKY9pCoIv/FrSypw0AzZo1Q+PGjTXz7UFIHfD1noZl/wMIEx+rds//zJkzCAoKQosWLWBubl7n7ZaNVXt4eGD3oSTse/Ie5xM/kfpRLzdB8PT0rI+XIaTW+HhPQ1WHbGJiYhAcHIyePXvi5MmTMDMzq9c6xGIxNu/8mQJbR9RLaMfHx9fHyxBSJ3y7p2FVZ3JE+qQhKnK6fJmsrCxMnDgRoaGhiI2NrfZ/snw/M4Zoh9qhvW3btkrPJSQk1EsxhNQVn05hq2nIRiaTwdzcHBcuXMAPP/yARo0aKX0tVYdZiO5TO7S/++47lZ4jpKFTNmRjZWWFRYsWYeLEiWCMwc3NDfr6+kpfR0hnxhDNU3lq1mvXruHKlSt49eoV1q9fL39eKpWisLBQI8URInRlwR0VOR17YrbAwsICERER2LhxI7755huVJkCqbpiluomiiG5SObSfPHmCa9euITc3F1evXpU/37x5c+zYsUMTtRGiE8qGbAoLCzFu3DgcOHAA27Ztw+TJk1Vaf3k1c1yThkfl0B46dCiGDh2KhIQEDBgwQJM1EaKTtm7dipiYGBw8eBChoaEqr1e2tx4eGoBIn7T/BTb3B1oJN9Qe0165ciVWrFiB69eva6IeQnRO2VWOM2bMwK+//qpWYJfh25kxhDtqh/bChQshlUoxc+ZMWFhYYMSIEdi8ebMmaiNE8J49ewYfHx+cOnUKenp6cHd3r/Vr8enMGMIdtUM7KCgIy5Ytw+nTp/H111/j2rVrmDNnjiZqI0TQHj58CD8/Pzx58gRt2rThuhyiI9QO7cWLF6NHjx4IDAzE77//jnXr1iErK0sTtREiWKmpqejRowdEIlGtJn4iRBm178a+ceNGODs7Y/r06ejbty8cHR01URchgsUYw4QJE2BhYYETJ07A0tKS65KIDlE7tF++fImUlBScPn0ac+bMwaNHj9C9e3ds2UKnHxEik8mgp6eHAwcOwNzcvN7nESGkVnOPtGrVCubm5mjRogWysrIUztsmpKH6+eef4e/vj5ycHDg4OFBgE41QO7RdXV3h5+eH5ORkDBw4ELdu3UJKSooGSiNEOLZu3YqRI0fC1tYWRkZGXJdDdJjawyNHjhyBs7OzJmohRJC++eYbLFiwABMmTADLz0JWVhadlkc0RuU97cePHwNAtYH99OnTuldEiIBcvXoVCxYswNy5c5H550WMoFn4iIapHNpjxozBtGnTcPbsWeTn58ufT09Px/r169GtWzdcvnxZI0USokm1maeasdLpVrt27YqEhAT8fvkY1g29T7PwEY1TObSTk5Px3nvvYenSpWjZsiVatGgBExMTBAUF4eHDhzhw4ACGDRumyVoJqXe1mae6sLAQo0aNwurVqwEAe7d/V+PNDgipL2qNaY8YMQIjRoxAcXExXr16hcaNG6N58+aaqo0QjSo/T7WzlQhOlqV7yNXN7ZGbm4thw4YhMTERY8eOBUCz8BHtqtUpfwYGBrCysqLAJoKmyu3Aynv9+jX69u2LixcvIiEhASEhIQD4eX9Korvq5R6RhAiRundwj4qKwt27d/HLL78gKChI4Wc0Cx/RFhErO6KiY7Kzs2FqagqpVCrI/xEwxpCXlwcTExOV7m7CR0LooWyIRNk81YwxvH37Fo0bN0Z2djaePXsGFxcXDitWnxA+h+oIvX6gcg91ySfa0yYNWk17yKmpqfDz88P9+/dhamoquMAmuketA5HPnz/H+vXrcePGDQCAp6cnZs6cCSsrK40UR4g2lM1TXdGlS5cQHBwMOzs7NG3alIPKCKlM5T3tO3fuoHPnzrhz5w569+6NoKAg3LlzB+7u7rhz544mayRE606ePIk+ffrAzc0Nx48fp5n6CG+ovKc9f/58rFmzBqNHj1Z4fu/evfj4448RHx9f78URwoW///4bI0eORGBgIPbv3y/YcVSim9Ta064Y2EDplZJ3796t16II4QpjDGZmZjhz5gwOHz6Mxo0bc10SIQpUDu3qTjLR0RNQSAOzYsUKTJ48GTKZDF5eXjA0NOS6JEIqUTm0XV1dsW/fvkrP//TTTyrP+jdnzhzY29tDJBIhNTVV/vyLFy/Qv39/ODk5oVOnTkhOTq5y/fT0dBgYGMDDw0P+5/79+6q2QEiVGGP45JNPEBUVBTs7OxoOIbym8pj2N998g8DAQMTExMDX11d+77vExEQkJiaq9BphYWFYsGAB/Pz8FJ6PioqCj48Pjh8/jqtXryIsLAz379+HgUHl8szMzGj+blJvSkpKMHPmTGzZsgWrVq1CZGQk1yURUi2VQ7t9+/ZISUnB+vXrcfLkSQClp/ytWrVK5TtN+/v7V/n8/v378fDhQwCls6ZZWloiOTkZgYGBqpanlFQqVRi+MTIyEsQk9WU1C3noSQg9bNu2Ddu2bcP27dsxceLESrUKoYeaCL0HodcPVO6hLr2odZ62WCzGv//9bwBAcXFxlXvC6srKyoJMJkPr1q3lz9nb2yMjI6PK5bOzs9G1a1eUlJQgJCQEixYtgr6+vtLXt7OzU3i8cOFCLFq0qM51a0v5aXCFio89MMYgEokwevRoODs7o2vXrsjLy1O6PB97UJfQexB6/cA/PVT3u1YTtVP39u3bGDt2LLKysvD48WNcv34d+/fvx4oVK2pdRMUxRGXfQmKxGJmZmbCwsMDr168xatQorFy5EgsWLFD62hkZGQqXiQppTzs/Px/GxsaCHWPlaw+vX7/G8OHD8emnn6J3795K/wcI8LcHdQi9B6HXD1TuoaioqNavpfZl7LNmzcJ3332HVq1aASgdIqnLOdrm5uYASu/yXubRo0eV9pCB0sC1sLAAALRs2RJTpkzB+fPnq319U1NThT9lb5oQ/gDgvAZd60EikSAwMBCpqakwMzOrcflnz57h/YnD8ezZM85r16XPoaHVX1UPtaV2aOfk5CgcSBSJRHU+NWrEiBFYt24dgNLbNz179qzSwUqg9CyTsm+ogoICHDp0CF26dKnTtknDcf/+ffj5+eHvv//G+fPn4e3tXe3yEokE44cFYJT1CYwfRneiIfygdmgbGBigqKhI/k2RmZkJPT3VXmbWrFmwsbFBZmYm+vTpA0dHRwCl58devHgRTk5OmDRpEnbv3i0fL//ss8+wceNGAKV3z+nSpQvc3d3h6ekJKysrQY1PE+4wxjB27FgYGhriwoULaN++fbXL/3ODhLJbiN2nW4gRXlB7atY9e/bgp59+wq1btzBlyhTs2rULy5Ytq/JqSS7R1Kzc40sPZQcd7969CzMzM/kQW3UmjhqEMPExDHT/p+743xgOSoKrnFyKz/jyOdSW0OsH6ndqVrUPRIaHh6Ndu3aIjY3F27dvsXPnTvTs2VPdlyFEK06cOIEVK1YgNjZW5YvAALqFGOEvlYdHIiIi5H9/+fIlVqxYga+//poCm/DW/v37MXjwYDRp0kTh9FRV7r7+zy3EHP53CzEHuiMN4QWVQ/vy5cvyv0dHR2ukGELqy+bNmzF69GiMGjUKhw4dgomJCQD17r4uFoux+1AS9j15D7sPUWATfqjVhFFCvjKJ6L6UlBTMmDEDs2fPxs6dO+VnN5W/+3rpwcU0lYJ7886fKbAJb6g8pl1QUID//ve/YIwp/L1Mhw4dNFIgIaoqO+Do4eGBc+fOwc/PT+HAVXV3XxfawUXScKkc2m/fvkVwcLD8cfm/i0QiPHjwoH4rI0QNxcXFmDFjBtzd3TFnzpwqj7XQwUWiC1QO7fT0dA2WQUjtFRQUYOzYsYiNjcWOHTuULld2cLG6u68Twnd0N3YiaG/evMGgQYMQHx+PmJgYhIeHV7t8TXdfJ4Tv6j5NHyEcioqKwpUrV3DixAkEBASotI6yu68TIgS0p00Eqewg+Jdffonz58+rHNiECB2FNhGctLQ09OzZEw8fPoSZmRnc3d25LokQraHQJoJy69Yt+Pn54dWrV9Xe/IIQXUWhTQTjwoUL8Pf3h7W1Nc6fP1/lnOuE6DoKbSIIOTk5GDJkCDw8PHD27FmF29MR0pDQ2SOE9xhjaNasGeLj4+Hu7i6fR4SQhoj2tAmvbdy4EdOnT4dMJoOPjw8FNmnwKLQJLzHG8NVXX2HmzJlo0qQJ1+UQwhsU2oR3GGOYP38+Fi5ciOjoaKxevVrlW9oRoutoTJvwzu7du7Fy5Up8++23+OCDD7guhxBeodAmvDNu3DjY2dkhMDCQ61II4R36PyfhhZycHAwaNAhnz56Fvr4+BTYhSlBoE85lZWWhd+/eOH/+PF3lSEgNaHiEcCozMxP9+vXDq1evkJiYiC5dunBdEiG8RqFNOMMYw5gxY5Cbm4vk5GQ4OztzXRIhvEehTTgjEomwZcsWNGnSBLa2tlyXQ4gg0Jg20brz589jwIAByM3Nhaurq1YCWyKRYOKoQdXeeZ0QIaDQJloVHx+Pfv36oaCgACUlJVrZpkQiQXhoAMLExxAeGkDBTQSNQptozY8//oiQkBC89957OHbsGJo3b67xbZYF9rqhaRjoLsK6oWkU3ETQKLSJVvzxxx8IDw/HuHHjcPDgQRgbG2tlu1GR0xHpkwZnKxEAwNlKhEifNERFTtfK9gmpb3QgkmhFx44dcfz4cfTp00er84gsX70F4aEBcLIsDe67zxhWX3bEnpgtWquBkPpEe9pEY2QyGT755BOsW7cOANCvXz+tT/wkFouxJyYJs2IdEf8bw6xYR+yJSYJYLNZqHYTUFwptohHFxcWYOnUq1q1bB5FIxGktZcF9UBJMgU0Ej4ZHSL3Lz8/H6NGjER8fj61bt2LSpElclwSxWIyd++K4LoOQOqPQJvVu0aJFOHHiBGJiYhAUFMR1OYToFBoeIfXu008/xdmzZzFw4ECuSyFE51Bokzopu9Lw+vXr6NevHx49eoQWLVrAx8eH69II0UlaDe05c+bA3t4eIpEIqamp8udfvHiB/v37w8nJCZ06dUJycrLS14iLi4OrqyscHR0xfPhwvHnzRhulkyqUv9JwRLAPbt++jcLCQq7LIkSnaTW0w8LCkJycjLZt2yo8HxUVBR8fH9y7dw/bt2/HuHHjUFxcXGn9N2/eYOrUqTh8+DDS0tIgFouxdOlSbZVPyql4peGxyBI4WBmiadOmXJdGiE7T6oFIf3//Kp/fv38/Hj58CADo2rUrLC0tkZycXOnuJQkJCfD29oarqysAICIiAsHBwfjqq6+UblMqlYIxJn9sZGQEIyOjOnaieWU1l6+dTz6p4krDeX6P8EnkdOz86SgA/vegCuqBe0KvH6jcQ1164fzskaysLMhkMrRu3Vr+nL29PTIyMiotm5GRobCXbm9vjydPnkAmkym9aMPOzk7h8cKFC7Fo0aJ6ql7z8vPzuS6hSkuWrcG0MXfgZPlAfqXhfy62w/d71yAvL09hWb72oA7qgXtCrx/4p4eK/0bUwXloA6h08UV130LqXqiRkZGhMDGRkPa08/PzYWxszPnFKRX98MMPSE5Oxp6YJEwYHohIn/tYfdkBPxxWvHCFzz2oinrgntDrByr3UFRUVOvX4jy0zc3NAQAvX76U720/evSo0h4yULrX/Msvv8gfp6enw9rautpLo01NTbUym5ymiEQiXv2ifvfdd/jggw8wadIkWFlZYU/MOURFTseemC1KrzTkWw+1QT1wT+j1A//0UJc+eHHK34gRI+TzU1y9ehXPnj2Dn59fpeX69++Pq1ev4s6dOwCA9evXY/To0VqttaFijOGLL77ABx98gI8++ghbt26Fvr6+/EpDujScEC1hWhQREcGsra2Zvr4+s7S0ZA4ODowxxp49e8b69u3LHB0dWYcOHVhiYqJ8ncWLF7MNGzbIH8fGxjIXFxfm4ODAQkJCmFQqrXJbUqmUAVD6c76TyWQsNzeXyWQyrkthjDG2d+9eBoAtXbpU5Zr41kNtUA/cE3r9jFXuoS75JGJMwIdkq5GdnQ1TU1NIpVJBDo8wxpCXlwcTExNe/JewuLgYJ0+eRHBwsMrr8K2H2qAeuCf0+oHKPdQln3gxPEL4qWzip3PnzsHAwECtwCaEaAaFNqlSdnY2BgwYgCNHjiA3N5frcggh/8P52SOEf16+fIkBAwYgLS0NJ0+erPKgMCGEGxTapJLRo0cjMzMTSUlJcHd357ocQkg5FNqkkm+//RaNGjWCk5MT16UQQiqgMW0CALhx4wZCQ0Px9u1bdOzYkQKbEJ6i0CZISkpCYGAgnj59qhPzOxCiyyi0G7ijR4+if//+ePfdd3H69Gm0bNmS65IIIdWg0G7A7t69i9DQUAQHByM+Ph7NmjXjuiRCSA0otHVI2a2/JBKJSss7Ozvj0KFD2LdvnyBmPiSEUGjrjPK3/goPDVAa3IwxREdHY+PGjQCAIUOGwMCATiIiRCgotHVAxVt/rRuaVmVwy2QyREZG4vPPP8fff//NTbGEkDqh0NYBUVXc+ivSJw1RkdPlyxQVFWHixIlYu3YtNm7ciKioKK7KJYTUAYW2Dli+egtWX3bE3WelEzbefcaw+rIjlq/eIl9myZIl2LdvH/bu3YsZM2ZwVSohpI4otHWAWCzGnpgkzIp1RPxvDLNiHbEnRvHWX/PmzcOJEycwatQoDislhNQVhbaOKAvug5JgeWC/fPkSISEhyMjIQMuWLdGrVy+uy6ySume9ENKQUWjrkPK3/srIyEDPnj1x+fJlSKVSrktTStWzXgghpSi0ddCdO3fQo0cPFBYW4sKFC3Bzc+O6pCqpetYLIeQfFNo6Jj8/H3379oWpqSmSk5Ph4ODAdUlKqXLWCyFEEYW2jjE2NsaOHTtw7tw5tGnThutyqqXKWS+EEEUU2jriyJEjmDt3Lhhj6N27tyAmflLlrBdCiCIKbR2wa9cuDBs2DE+fPkVxcTHX5ailqrNeCCHKUWgL3Jo1azBx4kRMnjwZP/30EwwNDbkuSW3lz3ohhFSPQlvADh8+jMjISCxYsACbN2+Gvr4+1yURQjSMpncTsEGDBuHAgQMICwvjuhRCiJbQnrbAFBUVYfr06UhOToaBgQEFNiENDIW2gOTl5SE0NBQ7d+7E8+fPuS6HEMIBGh4RCKlUiiFDhuDatWs4evQo3nvvPa5LIoRwgEJbIMaOHYtbt27h9OnT8PX15bocQghHKLQFYvny5QDA23lECCHaQWPaPHbnzh2MHTsWb9++hZubGwU2IYT2tPnq2rVrGDBgAMRiMXJyctC4cWOuSyKE8ADtafPQ2bNnERQUBAcHByQmJsLS0pLrkgghPEGhzTMPHjzAgAED0L17d8TFxQli4idCiPbQ8AjPtGvXDrt27cKQIUMgk8m4LocQwjO0p80Tq1atwubNmwEAI0eOhJGREccVEUL4iDehffz4cXh7e6Nz587w8fHBb7/9VmmZxMRENG7cGB4eHvI/eXl5HFRbfxhj+PTTT/HRRx8hPT2d63IIITzHi+GRv/76C+Hh4Th//jzat2+PpKQkjBs3DqmpqZWW7dChA65du8ZBlfVPJpNh9uzZ2LBhA1asWIEFCxZwXRIhhOd4sad9//59WFhYoH379gCAgIAAPHr0CDdu3OC4Ms3697//jU2bNmHLli0U2IQQlfBiT9vJyQkvX77E5cuX4ePjg5iYGLx58wbp6enw9PRUWPbPP/+Ep6cn9PX1MXnyZERERFT72lKpFIwx+WMjIyPejBdHRETAy8sLAwcOVKgRgPxxxeeFhHrgB6H3IPT6gco91KUXXoS2qakpfv75Z0RFRSEnJwd+fn7o0KFDpbuweHp6IjMzE6ampsjMzERwcDBatWqFkSNHKn1tOzs7hccLFy7EokWLNNKHKv7++2/MnTsXS5cuhY2NDYKCgqodl8/Pz9didZpBPfCD0HsQev3APz3U5ViciPHw66ugoABWVla4evUqHB0dlS731Vdf4enTp1i7dm2ln2VnZ8PU1BQZGRlo3ry5/Hku97SfP3+O/v37IyMjAydOnIC3t7fSZRljyM/Ph7GxMUQikRarrD/UAz8IvQeh1w9U7iE7OxtmZmaQSqUK+aQKXuxpA4BEIpHfI/Df//43goKCKgW2RCKBpaUl9PT0kJOTg7i4OEydOrXa1zU1NVX7TdGE9PR09O3bF7m5uTh37hw6deqk0noikUiwv6hlqAd+EHoPQq8f+KeHuvTBiwORALB48WK4urrC0dERjx49wtatWwEA06ZNw5EjRwAAP//8M9zc3ODu7g4fHx/07dsXkydP5rJslRQUFCAoKAgymQwXLlxQObAJIaQiXg6P1Iey4ZHa/PdDExISEuDh4aHyHccZY8jLy4OJiYlg9y6oB34Qeg9Crx+o3ENd8ok3e9q66MyZM5g3bx4YY/IZ+wghpC4otDUkJiYGwcHBSE1NRWFhIdflEEJ0BIW2Bmzbtg1hYWEIDQ3FkSNHeHNeOCFE+Ci061lCQgKmTp2K6dOn44cffkCjRo24LokQokN4c8qfrujTpw927tyJ8ePHC/agCSGEv2hPux6UlJTgww8/xMWLF2FoaIgJEyZQYBNCNIJCu44KCwsxbtw4fPvtt7h37x7X5RBCdBwNj9TB27dvERYWhjNnzmD//v0YPnw41yURQnQchXYdTJgwAefOnUN8fDz69OnDdTmEkAaAhkfqYPHixTh9+jQFNiFEayi01fTw4UNMnjwZeXl58jlQCCFEW2h4RA1//PEH+vXrBxMTE2RlZcHGxobrkgghDQztaavoypUr8Pf3R6tWrZCcnEyBTQjhBIV2BRKJBBNHDYJEIpE/l5GRgd69e8tvOmxlZcVhhYSQhoxCuxyJRILw0ACEiY8hPDRAHtx2dnZYv349Tp48CTMzM26LJIQ0aBTa/1MW2OuGpmGguwjrhqZhcJAHVq1aBaD09L7GjRtzXCUhpKGj0P6fqMjpiPRJg7NV6eXnzlYiLHnvJfZs/Zbjyggh5B8U2v+zfPUWrL7siLvPSm/kc/cZw+L4Fjh68gLHlRFCyD8otP9HLBZjT0wSxu5oifjfGMbuMEfCuVS0adOG69IIIUSOQrscsViMXT8n4tsbXRCfeItuD0YI4R26uKaCTp064dS561yXQQghVaI9bUIIERAKbUIIERAKbUIIERAKbUIIERAKbUIIERAKbUIIERAKbUIIERAKbZ4qKCjA0qVLUVBQwHUptUY98IPQexB6/UD99iBijLF6qIl3srOzYWpqCqlUiubNm3NdjtqkUinMzMzw999/w9TUlOtyaoV64Aeh9yD0+oHKPdQln2hPmxBCBIRCmxBCBERn5x4pG/XJzs7muJLaKas7OzsbIpGI42pqh3rgB6H3IPT6gco9lD2uzei0zo5pZ2ZmwtbWlusyCCFEqcePH6t9k3CdDW2ZTIanT5+iWbNmgv12JoToJsYYcnJy0KZNG+jpqTdKrbOhTQghuogORBJCiIBQaBNCiIBQaPPIwoUL0b59e7i7u+Pdd9/FL7/8onTZK1euwMPDA87OzujduzckEokWK1Vu27ZtcHNzg4GBAb777rtqlxWJROjcuTM8PDzg4eGB8+fPa6nK6qnTAx8/h7dv32LMmDFwdHSEs7MzDh06pHRZPn0G9+7dQ/fu3eHs7Ix3330Xt2/frnK5rVu3wsnJCQ4ODnj//fdRXFys5UqVU6WHxMRENG7cWP6ee3h4IC8vT/WNMMIbx44dY2/fvmWMMZaSksLMzMxYXl5epeVkMhlzcHBgZ8+eZYwx9s0337DRo0drs1SlUlJS2O3bt9n48ePZ2rVrq10WAMvJydFSZapTtQe+fg7R0dFs4sSJjDHGHjx4wCwtLdnr16+rXJZPn0GvXr3Y9u3bGWOMHThwgPn4+FRa5sGDB0wsFrNnz54xmUzGBg8ezDZu3KjlSpVTpYezZ88yLy+vWm+DQpunSkpKWLNmzdjjx48r/ezXX39lHTp0kD/Ozs5mxsbGrLCwUJslVmvixImCDe0yNfXA18+hQ4cO7Ndff5U/HjFihDxIKuLLZ/D8+XNmamrKioqKGGOlX4iWlpbs4cOHCst9/fXXLCIiQv44Pj6eBQQEaLFS5VTtoa6hTcMjPLV9+3Y4ODhUeQ5nRkYG2rZtK3/crFkzNGvWjBf/NVdXYGAg3N3d8dFHHyE3N5frctTC18+hYl329vbIyMhQujwfPoPHjx+jTZs2MDAovd5PJBLBzs6uUt3q9qZNqvYAAH/++Sc8PT3RtWtXrF+/Xq3t6OwVkXzUs2dP/Pe//63yZzdv3pRfDHTmzBlER0fj1KlTSl+r4rnnTEtnbqragyoePXoEOzs75Obm4l//+hfmz5+v9i9wbdRnD1x8DjXVX7Gu6mri6jOoiqrvpaq9cUGVHjw9PZGZmQlTU1NkZmYiODgYrVq1wsiRI1XaBoW2FqlykCcpKQmTJ0/G0aNH4eLiUuUydnZ2SE9Plz/OyclBTk4OxGJxfZWqVH0eqLKzswMANGnSBBEREXj//ffr7bWrU189cPU51FR/WV2tW7cGUBrMwcHBSpcFtP8ZVGRra4vMzEwUFxfDwMAAjDE8fvxYXl+Ziu952ZcOH6jaQ/lZ/WxsbDBmzBicP39e5dCm4REeOXfuHMaPH4/Y2Fi4u7srXc7Lywv5+flITEwEAGzatAkhISEwNDTUUqV199dff+Ht27cASq9e3bdvH7p06cJxVerh6+cwYsQIrFu3DgDw8OFDJCUlYciQIZWW49NnYGFhgS5dumDPnj0AgJ9//hn29vawt7dXWG748OGIiYnB8+fPwRjDxo0bMXr0aA4qrkzVHiQSCWQyGYDSL/q4uDj13vdaj4aTeufo6MgsLCyYu7u7/M+tW7cYY4xt2LCBLV68WL7sxYsXWefOnZmTkxMLDAxkmZmZXJWtYPfu3cza2po1btyYmZmZMWtra3bjxg3GmGIPFy9eZG5ubqxz586sQ4cOLDw8nGVlZXFZupyqPTDGz8/hzZs3bOTIkczBwYE5OTmxAwcOyH/G58/gzp07zMfHhzk5OTEvLy+WmprKGGNs6tSpLDY2Vr7c5s2bmYODA3vnnXfY1KlTOT/wW54qPaxdu5Z16NBB/r4vWbKEyWQylbdBl7ETQoiA0PAIIYQICIU2IYQICIU2IYQICIU2IYQICIU2IYQICIU2IYQICIW2AMyZMwf29vYQiURITU2tcfn8/Hx06NAB3t7eCs9nZGRg8ODBcHFxgaurK9auXYunT5/Kp4d0dHRUmDLyww8/xLRp0xSuwIuLi0PXrl3h4uKCdu3a4f3334dUKpX/fMqUKXBxcYGHhwf8/f2RkpJSb+9DVVSdprO6aW/XrVsHNzc3eHh4wM3NDd9++22l9T/77DPo6+vj0aNHCs/b29vLP5Oaptat6b2bPHmyfJrUrl274syZM7V+X1Sh6ntX3VS1c+bMUZhi1NjYWOH9u379Ovr374927dqhU6dO8PX1xeHDh2usrbCwEJ988gkcHR3Rvn17dOrUCdu3b69TvzpDA+eXk3qWlJTEHj9+zNq2bct+//33Gpf/6KOP2JQpUxRmEpPJZMzT05Pt379f/lgikSisV9PsYwkJCaxNmzbyC02KiorY7NmzWY8ePeQXB8TGxspnOTt69ChzcnJSqce2bduqtFx56kzTWd20t3///bd8OalUymxtbdlvv/0mf66kpITZ2tqywMBA9vnnn1equ+wzqW4bqrx3f/31l/x1b968yczNzVW66ELT752qU9VKJBJmbGws/71KTU1l5ubm7MiRI/JlMjMz2Y4dO2qsb8yYMWzYsGHszZs3jDHGHj58yFxdXdmmTZvUaVMnUWgLiCqhfe7cOTZ48OBKAXzq1CnWo0ePatetKrQDAgLY0aNHGWOM9ejRg61bt07h5wUFBcza2pqdPn260uu9fPmSNWrUiJWUlFS7XcZqFzy1naazumlvJRIJE4vF8itRGSsNXG9vb3b9+nXWtm1bhX6UfSYVt6Hue3f27FnWqlUrjYV2bd67mqaqXbFiBRs6dKj8cXh4OJs/f77S5VeuXMm8vb2Zh4cH69q1K7t8+TJjjLF79+4xExMT9urVK4Xl4+PjmY2NTbU1NgQ0PKJDcnNzERkZiQ0bNlT62e3bt9G6dWuMHj0aXbp0QWhoKB48eKDW69+4cQO+vr4KzzVq1AheXl64ceNGpeXXrFmD4OBgte82raraTtNZ1bS3Bw8eRMeOHdG2bVvMnz8fbm5u8p9t3boVU6ZMgaenJ1q0aKHSsEXFbaj63kVFRcHBwQHDhg3DgQMHKs0aV180McXptm3bMHXqVPnj69evV+q5vPHjx+Pq1au4efMmvv32W/m6N27cgJOTE8zNzRWW9/X1RWZmJl6+fFmnOoWOZvnTIfPnz8esWbNgbW2Ne/fuKfysqKgIp0+fxuXLl9GxY0ds3rwZo0ePxq+//qrWNqoKEVbFTAh79uzB/v37lc5IV1JSAi8vL/njsrF1ALC0tMSJEyfUrqeqOipSNu1tWFgYwsLCkJ6ejtDQUAQHB8PFxQWvXr3CqVOnsGXLFgDA1KlTsXXrVvTt21ftbajy3i1fvhzLly/H6dOnMX/+fFy4cAGNGjVSWIar9646Fy5cQHZ2ttLZBKty8+ZNLF26FFlZWTAwMMDt27dRWFhYqTaiiEJbhyQnJ+PYsWP44osvkJ+fj7/++gsdO3bEH3/8gbZt26JLly7o2LEjACA8PBwzZ85ESUkJ9PX1VXp9T09PXLx4UR4QQOkBoxs3bmDu3Lny5/bt24fo6GicOXMGFhYWVb6Wvr6+wkFKe3v7Gg9ahoWFIS0tDUBpMKo7Tacq097a29ujW7duiIuLg4uLC3bv3o3i4mJ5zyUlJcjKykJWVlalPcHqtqHqe1emT58+mD17Nn7//XeFgAa4ee9qsnXrVkycOFHhd8nLywuXLl1CaGhopeULCwsxfPhwJCYmwsvLC9nZ2TA1NUVhYSG6dOmCu3fvVnqPL126BBsbG/mUsw0Wt6MzRB2qHohkrPL49Js3b1i7du3ks9D9/PPPrHPnztWuw5jimHZ8fDwTi8Xs5s2bjLF/Dqb5+vrKx1737dvHHB0dWXp6utq9qev+/fuVDqZt2LChymWTkpKYra2t/EBgebdv35b//cWLF8zR0ZGdPHmSMcZYp06dWEJCgsLyw4cPZ2vWrJHXXfaZVLeNmt67oqIidvfuXfnyV65cYS1atFB6b8fyNP3elVE2pp2Tk8OaNWumUD9jjN26dYuZm5uzuLg4+XOPHz9mmzZtYlKplBkZGbHnz58zxhj78ssvFW59NnLkSDZs2DCWm5vLGPvnQGRNNTYEFNoCEBERwaytrZm+vj6ztLRkDg4O8p9VnLayTFUBfPz4cebu7s46d+7M/P395dNGVrdO+dBmjLHDhw8zT09P5uzszOzt7dnUqVMVznowMDBgNjY2CtPLVjygVJXaBA9jyqfpfPLkCXN3d5cvV920tzNnzmQdOnSQvzdlBwwvX77MLCws5GfDlH8Pyr7wyod2ddsoW0/Ze5efn8+6d+/OOnbsyDp37sx8fX3ZmTNnVHoPNP3eVTdVLWOMff/998zf37/Kbfz666+sb9++7J133mGdOnVi3bt3l/++rlixgrVt25b17NmTffPNNwqhnZ+fz+bNm8fatWvHXF1dWceOHdmWLVtq1aeuoalZCSFEQOjsEUIIERAKbUIIERAKbUIIERAKbUIIERAKbUIIERAKbUIIERAKbUIIERAKbUIIERAKbUIIERAKbUIIERAKbUIIEZD/D3Ll3SZe0khDAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "wholerock_file = \"./data/wholerock.csv\"\n", "wholerock = mp.read_melt(wholerock_file, index_col=[\"name\"])\n", "\n", "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)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we need to set up the object that handles the random sampling of errors in the Monte Carlo simulation. This is done with the `PEC_MC_parameters` class and it includes the following parameters for error propagation:\n", "\n", "- `melt_errors`\n", " \n", " propagate errors on melt composition by providing one standard deviation errors per element as a pandas Series (fixed errors for all inclusions) or DataFrame (errors per inclusion).\n", "\n", "- `olivine_errors`\n", "\n", " propagate errors on olivine composition by providing one standard deviation errors per element as a pandas Series (fixed errors for all inclusions) or DataFrame (errors per inclusion).\n", "\n", "- `pressure_errors`\n", "\n", " propagate errors on inclusion internal pressures by providing one standard deviations arrors as float or ints (fixed errors for all inclusions), or pandas series or numpy arrays (errors per inclusion).\n", "\n", "- `FeOi_errors`\n", "\n", " propagate errors on estimate melt initial FeO contents. Fixed errors can be provided either for the whole dataset, or per inclusion. Alternatively, an [FeOi_prediction object](https://magmapec.readthedocs.io/en/latest/notebooks/FeOi.html) can be provided to propagate errors on predictions models.\n", "\n", "- `Fe3Fe2`\n", "\n", " Propagate errors on modelled melt Fe2+/Fe3+ ratios. Pass *True* to this parameter to activate it. Errors are automatically calculated by MagmaPandas based on the the procedure outlined [here](https://magmapandas.readthedocs.io/en/latest/notebooks/Fe3Fe2_errors.html)\n", "\n", "- `Kd`\n", "\n", " Propagate errors on modelled olivine-melt Fe-Mg partition coefficients. Pass *True* to this parameter to activate it. Errors are calibration errors stated in the original publications - see [here](https://magmapandas.readthedocs.io/en/latest/models.html) for their values \n", "\n", "- `temperature`\n", "\n", " Propagate errors on modelled olivine liquidus temperatures. Pass *True* to this parameter to activate it. Errors are calibration errors stated in the original publications - see [here](https://magmapandas.readthedocs.io/en/latest/models.html) for their values\n", "\n", "\n", "By default errors are not propagated - you explicitely need to tell MagmaPEC to do so when initialising the `PEC_MC_parameters` object" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this example we will use all error propagation options, which means we need to provide analytical errors for all elements measured in melt and olivine. We import these from .csv files containing error data for individual inclusions and olivines. This is just an example with randomly generated errors, normally you should use analytical errors." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "melt_errors_file = \"./data/melt_errors.csv\"\n", "olivine_errors_file = \"./data/olivine_errors.csv\"\n", "\n", "melt_errors = pd.read_csv(melt_errors_file, index_col=[0])\n", "olivine_errors = pd.read_csv(olivine_errors_file, index_col=[0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Make very sure that the elements in the error data have identical sorting to the melt and olivine dataframes, otherwise errors will be applied to the wrong elements. \n", "We can force this by sorting the columns of the error dataframes (or series) with the *elements* attributes of the melt and olivine MagmaFrames:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "melt_errors = melt_errors[melt.elements]\n", "olivine_errors = olivine_errors[olivine.elements]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's what they look like:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
SiO2Al2O3MgOCaOFeONa2OK2OMnOTiO2P2O5Cr2O3CO2H2OFSCl
PI032-04-011.020.640.140.470.690.100.080.050.180.050.050.170.210.220.100.12
PI032-04-021.060.840.290.460.540.050.030.010.060.180.040.240.120.130.120.04
PI041-02-021.040.900.090.550.490.220.210.010.160.040.000.150.230.040.170.06
PI041-03-010.980.640.100.400.530.120.120.080.010.030.000.050.360.060.140.04
PI041-03-031.020.540.310.420.680.180.030.020.110.060.100.180.080.130.060.15
\n", "
" ], "text/plain": [ " SiO2 Al2O3 MgO CaO FeO Na2O K2O MnO TiO2 P2O5 \\\n", "PI032-04-01 1.02 0.64 0.14 0.47 0.69 0.10 0.08 0.05 0.18 0.05 \n", "PI032-04-02 1.06 0.84 0.29 0.46 0.54 0.05 0.03 0.01 0.06 0.18 \n", "PI041-02-02 1.04 0.90 0.09 0.55 0.49 0.22 0.21 0.01 0.16 0.04 \n", "PI041-03-01 0.98 0.64 0.10 0.40 0.53 0.12 0.12 0.08 0.01 0.03 \n", "PI041-03-03 1.02 0.54 0.31 0.42 0.68 0.18 0.03 0.02 0.11 0.06 \n", "\n", " Cr2O3 CO2 H2O F S Cl \n", "PI032-04-01 0.05 0.17 0.21 0.22 0.10 0.12 \n", "PI032-04-02 0.04 0.24 0.12 0.13 0.12 0.04 \n", "PI041-02-02 0.00 0.15 0.23 0.04 0.17 0.06 \n", "PI041-03-01 0.00 0.05 0.36 0.06 0.14 0.04 \n", "PI041-03-03 0.10 0.18 0.08 0.13 0.06 0.15 " ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "melt_errors.head()" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
SiO2FeOMgONiOMnOAl2O3CaO
PI032-04-011.060.601.220.010.090.080.02
PI032-04-020.850.571.280.140.070.240.06
PI041-02-021.110.441.300.080.100.010.19
PI041-03-010.990.501.310.140.160.150.01
PI041-03-030.950.641.210.130.110.100.15
\n", "
" ], "text/plain": [ " SiO2 FeO MgO NiO MnO Al2O3 CaO\n", "PI032-04-01 1.06 0.60 1.22 0.01 0.09 0.08 0.02\n", "PI032-04-02 0.85 0.57 1.28 0.14 0.07 0.24 0.06\n", "PI041-02-02 1.11 0.44 1.30 0.08 0.10 0.01 0.19\n", "PI041-03-01 0.99 0.50 1.31 0.14 0.16 0.15 0.01\n", "PI041-03-03 0.95 0.64 1.21 0.13 0.11 0.10 0.15" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "olivine_errors.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Together with the `FeOi_prediction` object, we pass these as arguments to the `PEC_MC_parameters` object. Errors on pressure are fixed at 2 kbar for al inclusions via the *pressure_errors* argument. We also set *Fe3Fe2* and *Kd* to *True* in order to propagate their model errors." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "mc_parameters = mpc.PEC_MC_parameters(melt_errors=melt_errors, olivine_errors=olivine_errors, pressure_errors=2e3, FeOi_errors=FeOi_predict, Fe3Fe2=True, Kd=True, temperature=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can create the Monte Carlo model with the `PEC_MC` object" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "pec_mc_model = mpc.PEC_MC(inclusions=melt, olivines=olivine, P_bar=pressure, FeO_target=FeOi_predict, MC_parameters=mc_parameters)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and run it with 20 iterations:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Monte Carlo iterations... |█████████████████████████| 100% [20/20] in 1:11.8 \n" ] } ], "source": [ "pec_mc_model.run(n=20)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Results are stored internally in the following attributes:\n", "\n", "- `pec`: pandas DataFrame\n", "\n", " Average PEC extents (%) of the MC model and their one standard deviation errors. Positive values indicate post-entrapment crystallisation and negative melting.\n", "\n", "- `inclusions_corr`: MagmaPandas Melt frame\n", "\n", " Averages of corrected melt inclusion compositions (wt. %)\n", "\n", "- `inclusions_stddev`: pandas DataFrame\n", "\n", " One standard deviation errors on inclusions_corr (wt. %)\n", "\n", "- `pec_MC`: pandas DataFrame\n", "\n", " PEC extents for individual iterations. Positive values indicate post-entrapment crystallisation and negative melting.\n", "\n", "- `inclusions_MC`: dictionary of MagmaPandas Melt frames\n", "\n", " corrected melt inclusion compositions for individual iterations.\n", "- `MC_inputs`: numpy array of dictionaries\n", "\n", " input parameters for each MC iteration. Each dictionary contains \n", " - 'inclusions'\n", " - 'olivines'\n", " - 'P_bar'\n", " - 'FeO_target' \n", " - 'Fe3Fe2_offset_parameters'\n", " - 'Kd_offset_parameters'\n", " - 'temperature_offset_parameters'\n", "\n", " if you want to view the parameters used for iteration _n_, use `pec_mc_model.MC_inputs[n]`, and if you want to run this iteration individually use `mpc.PEC(**pec_mc_model.MC_inputs[n])`" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "pec = pec_mc_model.pec\n", "inclusions_corrected = pec_mc_model.inclusions_corr\n", "inclusions_errors = pec_mc_model.inclusions_stddev\n", "\n", "pec_mc = pec_mc_model.pec_MC\n", "inclusions_MC = pec_mc_model.inclusions_MC" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
pecstddev
name
PI032-04-019.6852775.443236
PI032-04-0210.8687204.281761
PI041-02-020.3684952.965795
PI041-03-0113.5246235.176464
PI041-03-0313.0875505.469071
PI041-05-04-3.8418952.874289
PI041-05-062.5734863.331781
PI041-07-0112.8854254.738129
PI041-07-0211.9960614.764667
PI052-01-02-7.5497423.441786
\n", "
" ], "text/plain": [ " pec stddev\n", "name \n", "PI032-04-01 9.685277 5.443236\n", "PI032-04-02 10.868720 4.281761\n", "PI041-02-02 0.368495 2.965795\n", "PI041-03-01 13.524623 5.176464\n", "PI041-03-03 13.087550 5.469071\n", "PI041-05-04 -3.841895 2.874289\n", "PI041-05-06 2.573486 3.331781\n", "PI041-07-01 12.885425 4.738129\n", "PI041-07-02 11.996061 4.764667\n", "PI052-01-02 -7.549742 3.441786" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pec" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
SiO2Al2O3MgOCaOFeONa2OK2OMnOTiO2P2O5Cr2O3CO2H2OFSCl
name
PI032-04-0148.85473613.9899397.7332519.75635610.2724263.5528500.6642880.1405312.3845490.2700050.0000000.6207481.3732650.1372890.1653660.084401
PI032-04-0247.93847114.9850397.5100949.54591210.3225343.4521540.9140130.1397132.5679650.3525080.0000000.6777131.3097150.0810840.1544270.048657
PI041-02-0249.10267716.9513974.8830789.42506410.1277673.7368711.0230890.1552052.7792840.5523440.0000000.4416300.6142320.0537830.0996360.053943
PI041-03-0145.91713815.6758277.42410810.80675510.6721453.2458291.1207390.1181333.1046180.5394310.0000000.7996630.3666410.0945530.0562600.058159
PI041-03-0345.44474415.6062867.19650311.09544810.8621313.3575121.1639750.0869783.2452600.5122730.0000000.8221700.3013810.1272180.0675150.110604
PI041-05-0447.87757918.6716613.7798609.4575229.2990514.6303101.6365400.1627592.4995040.8246010.0000000.4038480.4843870.0939560.1235630.054859
PI041-05-0646.38538816.8106624.8213818.86049011.5240274.1974721.3683070.1289713.6582010.6292240.0000000.6160490.6741350.1280500.1240940.073549
PI041-07-0145.96181514.8608807.4255459.57862111.6918913.1131031.2658790.1675453.5212560.5479850.0000000.4671221.0054100.1222660.1754410.095240
PI041-07-0245.94133815.5890717.0360139.95262211.4276143.1660681.3552620.1395803.5010120.6113350.0000000.3608880.6570410.0729690.1266360.062551
PI052-01-0249.20734817.0565204.11808710.2760988.4939124.9071111.5412900.2413601.7467240.5984960.1185320.2869581.1246040.1107950.1087280.063437
\n", "
" ], "text/plain": [ " SiO2 Al2O3 MgO CaO FeO Na2O \\\n", "name \n", "PI032-04-01 48.854736 13.989939 7.733251 9.756356 10.272426 3.552850 \n", "PI032-04-02 47.938471 14.985039 7.510094 9.545912 10.322534 3.452154 \n", "PI041-02-02 49.102677 16.951397 4.883078 9.425064 10.127767 3.736871 \n", "PI041-03-01 45.917138 15.675827 7.424108 10.806755 10.672145 3.245829 \n", "PI041-03-03 45.444744 15.606286 7.196503 11.095448 10.862131 3.357512 \n", "PI041-05-04 47.877579 18.671661 3.779860 9.457522 9.299051 4.630310 \n", "PI041-05-06 46.385388 16.810662 4.821381 8.860490 11.524027 4.197472 \n", "PI041-07-01 45.961815 14.860880 7.425545 9.578621 11.691891 3.113103 \n", "PI041-07-02 45.941338 15.589071 7.036013 9.952622 11.427614 3.166068 \n", "PI052-01-02 49.207348 17.056520 4.118087 10.276098 8.493912 4.907111 \n", "\n", " K2O MnO TiO2 P2O5 Cr2O3 CO2 \\\n", "name \n", "PI032-04-01 0.664288 0.140531 2.384549 0.270005 0.000000 0.620748 \n", "PI032-04-02 0.914013 0.139713 2.567965 0.352508 0.000000 0.677713 \n", "PI041-02-02 1.023089 0.155205 2.779284 0.552344 0.000000 0.441630 \n", "PI041-03-01 1.120739 0.118133 3.104618 0.539431 0.000000 0.799663 \n", "PI041-03-03 1.163975 0.086978 3.245260 0.512273 0.000000 0.822170 \n", "PI041-05-04 1.636540 0.162759 2.499504 0.824601 0.000000 0.403848 \n", "PI041-05-06 1.368307 0.128971 3.658201 0.629224 0.000000 0.616049 \n", "PI041-07-01 1.265879 0.167545 3.521256 0.547985 0.000000 0.467122 \n", "PI041-07-02 1.355262 0.139580 3.501012 0.611335 0.000000 0.360888 \n", "PI052-01-02 1.541290 0.241360 1.746724 0.598496 0.118532 0.286958 \n", "\n", " H2O F S Cl \n", "name \n", "PI032-04-01 1.373265 0.137289 0.165366 0.084401 \n", "PI032-04-02 1.309715 0.081084 0.154427 0.048657 \n", "PI041-02-02 0.614232 0.053783 0.099636 0.053943 \n", "PI041-03-01 0.366641 0.094553 0.056260 0.058159 \n", "PI041-03-03 0.301381 0.127218 0.067515 0.110604 \n", "PI041-05-04 0.484387 0.093956 0.123563 0.054859 \n", "PI041-05-06 0.674135 0.128050 0.124094 0.073549 \n", "PI041-07-01 1.005410 0.122266 0.175441 0.095240 \n", "PI041-07-02 0.657041 0.072969 0.126636 0.062551 \n", "PI052-01-02 1.124604 0.110795 0.108728 0.063437 " ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "inclusions_corrected" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
SiO2_stddevAl2O3_stddevMgO_stddevCaO_stddevFeO_stddevNa2O_stddevK2O_stddevMnO_stddevTiO2_stddevP2O5_stddevCr2O3_stddevCO2_stddevH2O_stddevF_stddevS_stddevCl_stddev
name
PI032-04-010.6860300.9220661.4929550.5736100.3316830.1884750.0837190.0539390.1409330.0384650.0000000.1615330.2491220.1402330.1015470.069214
PI032-04-020.5568640.7219401.2734840.4220990.2611760.1506920.0501060.0132960.1350440.1448590.0000000.1986390.0950090.0765900.1140920.036706
PI041-02-020.6578860.9046550.9265790.5892700.3483150.2062510.2407450.0126010.1761040.0351350.0000000.1759570.1843040.0372340.1202620.059369
PI041-03-010.6074370.7670891.5506400.4915810.2168980.1392600.1119880.0627560.1171920.0322720.0000000.0596900.2563210.0553740.0624440.030225
PI041-03-030.7028750.7635891.6250720.6409120.2607890.1933490.0439440.0218750.1574450.0728490.0000000.1853210.0561510.1235190.0437920.114647
PI041-05-040.6728301.1680310.8486820.4398150.4007440.2422450.0916770.1099280.0771350.0255040.0000000.2060500.1078200.1005050.0087680.008581
PI041-05-060.7323371.2123650.9938450.4186590.4149830.3389460.2537960.1508570.1836790.0493700.0000000.2636630.2880440.1154290.0645890.026683
PI041-07-010.7373421.2720261.3941570.4870660.2696090.1414610.0522960.0614350.1535850.0279100.0000000.0362350.1573210.0910540.0273960.131493
PI041-07-020.5073480.8088091.3538310.4756400.2109960.1492890.1429770.0530070.1384950.0243980.0000000.2089350.2084260.0031210.0365000.061395
PI052-01-020.6782650.9778730.8919430.6144130.4789230.1615910.0779470.0226070.0862270.1622370.1392490.0090410.1552340.1279270.0431710.020256
\n", "
" ], "text/plain": [ " SiO2_stddev Al2O3_stddev MgO_stddev CaO_stddev FeO_stddev \\\n", "name \n", "PI032-04-01 0.686030 0.922066 1.492955 0.573610 0.331683 \n", "PI032-04-02 0.556864 0.721940 1.273484 0.422099 0.261176 \n", "PI041-02-02 0.657886 0.904655 0.926579 0.589270 0.348315 \n", "PI041-03-01 0.607437 0.767089 1.550640 0.491581 0.216898 \n", "PI041-03-03 0.702875 0.763589 1.625072 0.640912 0.260789 \n", "PI041-05-04 0.672830 1.168031 0.848682 0.439815 0.400744 \n", "PI041-05-06 0.732337 1.212365 0.993845 0.418659 0.414983 \n", "PI041-07-01 0.737342 1.272026 1.394157 0.487066 0.269609 \n", "PI041-07-02 0.507348 0.808809 1.353831 0.475640 0.210996 \n", "PI052-01-02 0.678265 0.977873 0.891943 0.614413 0.478923 \n", "\n", " Na2O_stddev K2O_stddev MnO_stddev TiO2_stddev P2O5_stddev \\\n", "name \n", "PI032-04-01 0.188475 0.083719 0.053939 0.140933 0.038465 \n", "PI032-04-02 0.150692 0.050106 0.013296 0.135044 0.144859 \n", "PI041-02-02 0.206251 0.240745 0.012601 0.176104 0.035135 \n", "PI041-03-01 0.139260 0.111988 0.062756 0.117192 0.032272 \n", "PI041-03-03 0.193349 0.043944 0.021875 0.157445 0.072849 \n", "PI041-05-04 0.242245 0.091677 0.109928 0.077135 0.025504 \n", "PI041-05-06 0.338946 0.253796 0.150857 0.183679 0.049370 \n", "PI041-07-01 0.141461 0.052296 0.061435 0.153585 0.027910 \n", "PI041-07-02 0.149289 0.142977 0.053007 0.138495 0.024398 \n", "PI052-01-02 0.161591 0.077947 0.022607 0.086227 0.162237 \n", "\n", " Cr2O3_stddev CO2_stddev H2O_stddev F_stddev S_stddev \\\n", "name \n", "PI032-04-01 0.000000 0.161533 0.249122 0.140233 0.101547 \n", "PI032-04-02 0.000000 0.198639 0.095009 0.076590 0.114092 \n", "PI041-02-02 0.000000 0.175957 0.184304 0.037234 0.120262 \n", "PI041-03-01 0.000000 0.059690 0.256321 0.055374 0.062444 \n", "PI041-03-03 0.000000 0.185321 0.056151 0.123519 0.043792 \n", "PI041-05-04 0.000000 0.206050 0.107820 0.100505 0.008768 \n", "PI041-05-06 0.000000 0.263663 0.288044 0.115429 0.064589 \n", "PI041-07-01 0.000000 0.036235 0.157321 0.091054 0.027396 \n", "PI041-07-02 0.000000 0.208935 0.208426 0.003121 0.036500 \n", "PI052-01-02 0.139249 0.009041 0.155234 0.127927 0.043171 \n", "\n", " Cl_stddev \n", "name \n", "PI032-04-01 0.069214 \n", "PI032-04-02 0.036706 \n", "PI041-02-02 0.059369 \n", "PI041-03-01 0.030225 \n", "PI041-03-03 0.114647 \n", "PI041-05-04 0.008581 \n", "PI041-05-06 0.026683 \n", "PI041-07-01 0.131493 \n", "PI041-07-02 0.061395 \n", "PI052-01-02 0.020256 " ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "inclusions_errors" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
namePI032-04-01PI032-04-02PI041-02-02PI041-03-01PI041-03-03PI041-05-04PI041-05-06PI041-07-01PI041-07-02PI052-01-02
iteration
011.46301310.8688962.15205113.9806159.322839-3.4156860.6189716.52055710.301697-6.992139
111.03521713.3327640.0002211.84698515.287744-4.3888790.83902611.16305511.717395-7.155634
22.8195134.801294-3.2656257.81112111.339624-6.0-1.07.61090111.0-7.909094
317.41918911.3571780.616.87378514.450586-3.6479491.81766413.93809812.245459-6.078774
40.05.996069-4.7143557.6357427.634424-8.244373-3.0029545.9393927.183643-12.570471
\n", "
" ], "text/plain": [ "name PI032-04-01 PI032-04-02 PI041-02-02 PI041-03-01 PI041-03-03 \\\n", "iteration \n", "0 11.463013 10.868896 2.152051 13.980615 9.322839 \n", "1 11.035217 13.332764 0.00022 11.846985 15.287744 \n", "2 2.819513 4.801294 -3.265625 7.811121 11.339624 \n", "3 17.419189 11.357178 0.6 16.873785 14.450586 \n", "4 0.0 5.996069 -4.714355 7.635742 7.634424 \n", "\n", "name PI041-05-04 PI041-05-06 PI041-07-01 PI041-07-02 PI052-01-02 \n", "iteration \n", "0 -3.415686 0.61897 16.520557 10.301697 -6.992139 \n", "1 -4.388879 0.839026 11.163055 11.717395 -7.155634 \n", "2 -6.0 -1.0 7.610901 11.0 -7.909094 \n", "3 -3.647949 1.817664 13.938098 12.245459 -6.078774 \n", "4 -8.244373 -3.002954 5.939392 7.183643 -12.570471 " ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pec_mc.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The dataframes in `inclusions_MC` also have *isothermal_equilibration*, *Kd_equilibration*, and *FeO_converge* columns. These columns show if equilibrations during stage [1] and [2] and FeO convergence in stage [2] were successful. Extreme cases of random error sampling may yield melt-olivine pairs that cannot be equilibrated without requiring crystallising exceeding the mass of the inclusion inclusion or exchange of more Mg or Fe than the inclusion contains. If that is the case, no corrected compositions are calculated and the isothermal- of Kd-equilibration column is set to *False*." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
SiO2Al2O3MgOCaOFeONa2OK2OMnOTiO2P2O5Cr2O3CO2H2OFSClisothermal_equilibrationKd_equilibrationFeO_converge
iteration
048.32068713.7310667.96183110.60829110.2978613.4354630.562640.1859082.4802710.2743580.00.5516621.0138080.4496790.1264750.0TrueTrueTrue
148.90236713.7232248.3755559.46586510.4528483.6619730.5933710.1630552.4379560.3062670.00.416021.201490.00.1236610.176348TrueTrueTrue
249.12249514.4684625.96606210.05964910.0307133.7905030.721590.1880322.2820570.3513950.01.0557211.56710.00.3736560.022564TrueTrueTrue
349.08646112.7778849.1793939.03334510.433093.2813050.6545250.2171692.2521730.2520450.00.7948181.8013230.00.1507810.085688TrueTrueTrue
449.57099915.3278175.5037059.9583539.948874.0419570.8941940.0717482.3742530.2467880.00.4637121.1516320.1253790.1666310.153963True<NA>True
\n", "
" ], "text/plain": [ " SiO2 Al2O3 MgO CaO FeO Na2O \\\n", "iteration \n", "0 48.320687 13.731066 7.961831 10.608291 10.297861 3.435463 \n", "1 48.902367 13.723224 8.375555 9.465865 10.452848 3.661973 \n", "2 49.122495 14.468462 5.966062 10.059649 10.030713 3.790503 \n", "3 49.086461 12.777884 9.179393 9.033345 10.43309 3.281305 \n", "4 49.570999 15.327817 5.503705 9.958353 9.94887 4.041957 \n", "\n", " K2O MnO TiO2 P2O5 Cr2O3 CO2 H2O \\\n", "iteration \n", "0 0.56264 0.185908 2.480271 0.274358 0.0 0.551662 1.013808 \n", "1 0.593371 0.163055 2.437956 0.306267 0.0 0.41602 1.20149 \n", "2 0.72159 0.188032 2.282057 0.351395 0.0 1.055721 1.5671 \n", "3 0.654525 0.217169 2.252173 0.252045 0.0 0.794818 1.801323 \n", "4 0.894194 0.071748 2.374253 0.246788 0.0 0.463712 1.151632 \n", "\n", " F S Cl isothermal_equilibration \\\n", "iteration \n", "0 0.449679 0.126475 0.0 True \n", "1 0.0 0.123661 0.176348 True \n", "2 0.0 0.373656 0.022564 True \n", "3 0.0 0.150781 0.085688 True \n", "4 0.125379 0.166631 0.153963 True \n", "\n", " Kd_equilibration FeO_converge \n", "iteration \n", "0 True True \n", "1 True True \n", "2 True True \n", "3 True True \n", "4 True " ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "inclusions_MC[\"PI032-04-01\"].head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Lets view the pressures used for different iterations:" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(name\n", " PI032-04-01 2820.134038\n", " PI032-04-02 2947.689858\n", " PI041-02-02 671.967541\n", " PI041-03-01 2264.566801\n", " PI041-03-03 2216.237876\n", " PI041-05-04 101.998027\n", " PI041-05-06 1227.736645\n", " PI041-07-01 2384.084584\n", " PI041-07-02 1164.839501\n", " PI052-01-02 119.941865\n", " Name: IaconoMarziano, dtype: float64,\n", " name\n", " PI032-04-01 4659.725542\n", " PI032-04-02 4787.281361\n", " PI041-02-02 2511.559044\n", " PI041-03-01 4104.158304\n", " PI041-03-03 4055.829380\n", " PI041-05-04 1941.589530\n", " PI041-05-06 3067.328149\n", " PI041-07-01 4223.676087\n", " PI041-07-02 3004.431004\n", " PI052-01-02 1959.533368\n", " Name: IaconoMarziano, dtype: float64)" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pec_mc_model.MC_inputs[1][\"P_bar\"], pec_mc_model.MC_inputs[4][\"P_bar\"]" ] }, { "cell_type": "code", "execution_count": null, "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": 2 }