PEC Monte Carlo simulation

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.

[2]:

import MagmaPEC as mpc
import MagmaPandas as mp

import pandas as pd

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 PEC model example

Import melt inclusion and olivine data:

[3]:

melt_file = "./data/melt.csv"
olivine_file = "./data/olivine.csv"

melt = mp.read_melt(melt_file, index_col=["name"])
olivine = mp.read_olivine(olivine_file, index_col=["name"])

Import inclusion internal pressures or calculate them if you have measured melt CO2 (and H2O). See the PEC model example for details on how to do the calculation. Here we import them from a file.

[4]:

pressure_file ="./data/pressure.csv"
pressure = pd.read_csv(pressure_file, index_col = ["name"]).squeeze()

Set up the melt initial FeO prediction model:

[5]:

wholerock_file = "./data/wholerock.csv"
wholerock = mp.read_melt(wholerock_file, index_col=["name"])

x = wholerock.drop(columns=["FeO"])
FeOi_predict = mpc.FeOi_prediction(x=x, FeO=wholerock["FeO"])

do_not_use = ["MnO", "P2O5", "Cr2O3", "total"]

model_fits = FeOi_predict.calculate_model_fits(exclude=do_not_use)
FeOi_predict.select_predictors(idx=3)

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:

  • melt_errors

    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).

  • olivine_errors

    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).

  • FeOi_errors

    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 can be provided to propagate errors on predictions models.

  • Fe3Fe2

    Propagate errors on modelled melt Fe2+/Fe3+ ratios. Errors are automatically calculated by MagmaPandas based on the selected model. Pass True to this parameter to activate it.

  • Kd

    Propagate errors on modelled olivine-melt Fe-Mg partition coefficients. Errors are automatically calculated by MagmaPandas based on the selected model. Pass True to this parameter to activate it.

By default errors are not propagated - you explicitely need to tell MagmaPEC to do so when initialising the PEC_MC_parameters object

In this example we will use all error propagation options, which means we need to provide melt and olivine composition errors. 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.

[6]:

melt_errors_file = "./data/melt_errors.csv"
olivine_errors_file = "./data/olivine_errors.csv"

melt_errors = pd.read_csv(melt_errors_file, index_col=[0])
olivine_errors = pd.read_csv(olivine_errors_file, index_col=[0])

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. We can force this by sorting the columns of the error dataframes (or series) with the elements attributes of the melt and olivine MagmaFrames:

[7]:

melt_errors = melt_errors[melt.elements]
olivine_errors = olivine_errors[olivine.elements]

Here’s what they look like:

[8]:

melt_errors.head()

[8]:
SiO2 Al2O3 MgO CaO FeO Na2O K2O MnO TiO2 P2O5 Cr2O3 CO2 H2O F S Cl
PI032-04-01 1.02 0.64 0.14 0.47 0.69 0.10 0.08 0.05 0.18 0.05 0.05 0.17 0.21 0.22 0.10 0.12
PI032-04-02 1.06 0.84 0.29 0.46 0.54 0.05 0.03 0.01 0.06 0.18 0.04 0.24 0.12 0.13 0.12 0.04
PI041-02-02 1.04 0.90 0.09 0.55 0.49 0.22 0.21 0.01 0.16 0.04 0.00 0.15 0.23 0.04 0.17 0.06
PI041-03-01 0.98 0.64 0.10 0.40 0.53 0.12 0.12 0.08 0.01 0.03 0.00 0.05 0.36 0.06 0.14 0.04
PI041-03-03 1.02 0.54 0.31 0.42 0.68 0.18 0.03 0.02 0.11 0.06 0.10 0.18 0.08 0.13 0.06 0.15 
[9]:

olivine_errors.head()

[9]:
SiO2 FeO MgO NiO MnO Al2O3 CaO
PI032-04-01 1.06 0.60 1.22 0.01 0.09 0.08 0.02
PI032-04-02 0.85 0.57 1.28 0.14 0.07 0.24 0.06
PI041-02-02 1.11 0.44 1.30 0.08 0.10 0.01 0.19
PI041-03-01 0.99 0.50 1.31 0.14 0.16 0.15 0.01
PI041-03-03 0.95 0.64 1.21 0.13 0.11 0.10 0.15

Together with the FeOi_prediction object, we pass these as arguments to the PEC_MC_parameters object. We also set Fe3Fe2 and Kd to True in order to propagate their model errors.

[10]:

mc_parameters = mpc.PEC_MC_parameters(melt_errors=melt_errors, olivine_errors=olivine_errors, FeOi_errors=FeOi_predict, Fe3Fe2=True, Kd=True, temperature=True)

Now we can create the Monte Carlo model with the PEC_MC object

[11]:

pec_mc_model = mpc.PEC_MC(inclusions=melt, olivines=olivine, P_bar=pressure, FeO_target=FeOi_predict, MC_parameters=mc_parameters)

And finally run the Monte Carlo simulation with n iterations. To get representative errors, a minimum of 50 iterations is recommended.

[12]:

pec_mc_model.run(n=50)

Monte Carlo loop
050/050
Equilibrating ... |██████████████████████████████| 100% [10/10] in 2.2s
Correcting    ... |██████████████████████████████| 100% [10/10] in 13.4s

Results are stored internally in the following attributes:

  • pec: pandas DataFrame

    Average PEC extents (%) of the MC model and their one standard deviation errors.

  • inclusions_corr: MagmaPandas Melt frame

    Averages of corrected melt inclusion compositions (wt. %)

  • inclusions_stddev: pandas DataFrame

    One standard deviation errors on inclusions_corr (wt. %)

  • pec_MC: pandas DataFrame

    PEC extents for individual iterations

  • inclusions_MC: dictionary of MagmaPandas Melt frames

    corrected melt inclusion compositions for individual iterations.

[13]:

pec = pec_mc_model.pec
inclusions_corrected = pec_mc_model.inclusions_corr
inclusions_errors = pec_mc_model.inclusions_stddev

pec_mc = pec_mc_model.pec_MC
inclusions_MC = pec_mc_model.inclusions_MC

[14]:

pec

[14]:
pec stddev
name
PI032-04-01 9.596699 4.571030
PI032-04-02 10.571102 4.999259
PI041-02-02 0.029676 3.332411
PI041-03-01 12.122017 5.592136
PI041-03-03 12.350934 5.439491
PI041-05-04 -4.324584 2.935303
PI041-05-06 1.761258 3.100318
PI041-07-01 11.755184 4.999059
PI041-07-02 10.548042 4.990837
PI052-01-02 -8.153084 2.932983 
[15]:

inclusions_corrected

[15]:
SiO2 Al2O3 MgO CaO FeO Na2O K2O MnO TiO2 P2O5 Cr2O3 CO2 H2O F S Cl
name
PI032-04-01 49.038185 13.965831 7.548374 9.712050 10.270587 3.587690 0.679658 0.131731 2.450270 0.282972 0.00000 0.601138 1.385841 0.138570 0.147149 0.059956
PI032-04-02 48.554467 14.659704 7.477362 9.371954 10.352409 3.458443 0.909402 0.135773 2.574248 0.318282 0.00000 0.591418 1.303731 0.086511 0.166056 0.040240
PI041-02-02 49.194844 16.869946 4.745780 9.219521 10.209481 3.815880 1.081526 0.150459 2.795936 0.565836 0.00000 0.469872 0.668820 0.059603 0.116337 0.036158
PI041-03-01 45.887674 15.778030 6.993274 11.021453 10.686971 3.348819 1.140917 0.110142 3.161059 0.555012 0.00000 0.794686 0.281567 0.086593 0.098787 0.055019
PI041-03-03 45.368272 15.883119 7.006869 11.180753 10.760509 3.410114 1.166775 0.085819 3.245916 0.500428 0.00000 0.807792 0.302871 0.112784 0.093646 0.074334
PI041-05-04 47.905434 18.758468 3.554904 9.522532 9.319479 4.639602 1.619916 0.130669 2.506853 0.830447 0.00000 0.488603 0.448569 0.090882 0.124218 0.059425
PI041-05-06 46.439516 16.934260 4.558820 9.115684 11.523588 3.992881 1.468084 0.171794 3.682809 0.630996 0.00000 0.645502 0.529145 0.113611 0.127881 0.065430
PI041-07-01 45.782869 15.466512 7.094071 9.641606 11.520774 3.130550 1.273328 0.146889 3.518325 0.558488 0.00000 0.465974 1.035865 0.099808 0.170116 0.094825
PI041-07-02 45.851318 15.577426 6.675896 10.248314 11.446601 3.240330 1.380481 0.139995 3.559769 0.619597 0.00000 0.393779 0.600101 0.074141 0.135176 0.057076
PI052-01-02 49.173977 17.128900 3.790439 10.484961 8.417258 4.937496 1.554145 0.235262 1.747572 0.716371 0.07623 0.288786 1.088326 0.169099 0.132997 0.058181 
[16]:

inclusions_errors

[16]:
SiO2_stddev Al2O3_stddev MgO_stddev CaO_stddev FeO_stddev Na2O_stddev K2O_stddev MnO_stddev TiO2_stddev P2O5_stddev Cr2O3_stddev CO2_stddev H2O_stddev F_stddev S_stddev Cl_stddev
name
PI032-04-01 0.684140 0.767229 1.429579 0.528847 0.318606 0.195741 0.082248 0.046226 0.169112 0.057255 0.000000 0.150842 0.203130 0.130984 0.077634 0.068188
PI032-04-02 0.675087 1.020384 1.429734 0.518800 0.277615 0.147678 0.044618 0.014049 0.106298 0.135802 0.000000 0.204367 0.137161 0.093134 0.085342 0.033133
PI041-02-02 0.921795 1.078272 1.085552 0.573999 0.378420 0.214054 0.177747 0.015982 0.136563 0.048991 0.000000 0.148134 0.277951 0.034226 0.113298 0.035241
PI041-03-01 0.683158 0.912316 1.736349 0.649194 0.159993 0.196658 0.120308 0.062801 0.140206 0.040101 0.000000 0.061839 0.228635 0.045563 0.098435 0.036972
PI041-03-03 0.621632 0.789395 1.602801 0.546556 0.215348 0.198945 0.064844 0.021529 0.151684 0.068454 0.000000 0.134119 0.070433 0.110846 0.052602 0.088702
PI041-05-04 0.767873 0.978518 0.898687 0.418429 0.343335 0.253467 0.127715 0.126536 0.093800 0.034276 0.000000 0.249720 0.167424 0.096160 0.010195 0.009643
PI041-05-06 0.859292 0.986402 0.975461 0.432861 0.399822 0.375694 0.213302 0.185424 0.197941 0.043942 0.000000 0.217330 0.236147 0.087190 0.059650 0.030361
PI041-07-01 0.577925 1.078946 1.408947 0.504997 0.248594 0.136789 0.045372 0.075074 0.143228 0.026259 0.000000 0.056172 0.185572 0.096999 0.027073 0.091682
PI041-07-02 0.535806 0.957769 1.528439 0.551789 0.227877 0.166253 0.147385 0.059076 0.138005 0.025977 0.000000 0.166085 0.174210 0.002745 0.045242 0.063498
PI052-01-02 0.906649 0.888955 0.877843 0.435415 0.402362 0.134263 0.073166 0.024505 0.097982 0.228803 0.097102 0.008038 0.151310 0.162627 0.060260 0.017870 
[17]:

pec_mc.head()

[17]:
name PI032-04-01 PI032-04-02 PI041-02-02 PI041-03-01 PI041-03-03 PI041-05-04 PI041-05-06 PI041-07-01 PI041-07-02 PI052-01-02
iteration
0 2.814795 1.06189 -7.713281 3.010327 2.955713 -9.099756 -0.522314 2.185376 0.01543 -14.604492
1 12.524048 4.575262 -1.423328 8.784961 8.0 -5.385857 1.880664 10.16709 10.0 -10.747681
2 3.789526 4.0 -0.382031 6.201294 6.136987 -7.401428 -1.350854 7.808716 6.187897 -11.27063
3 11.8 8.137012 -1.426172 7.44588 9.01008 -6.124573 -2.069324 10.4 9.422449 -8.559125
4 2.243164 2.316284 -2.730713 5.905371 6.262866 -8.778662 0.838452 6.921082 4.375781 -11.380933

The dataframes in inclusions_MC also have isothermal_equilibration and Kd_equilibration columns. These columns show if equilibration during stage [1] and [2] respectively was successful. In extreme cases, error propagation of melt and olivine compositions by random sampling can result in melt-olivine pairs that cannot be equilibrated without crystallising more than the entire inclusion or exchanging more Mg or Fe than is present in the inclusion. If that is the case, no corrected compositions are calculated and eiher the isothermal- of Kd-equilibration column is set to False.

[18]:

inclusions_MC["PI032-04-01"].head()

[18]:
SiO2 Al2O3 MgO CaO FeO Na2O K2O MnO TiO2 P2O5 Cr2O3 CO2 H2O F S Cl isothermal_equilibration Kd_equilibration
iteration
0 49.659755 15.346003 4.765847 9.638317 10.297186 3.761251 0.883066 0.117185 2.604499 0.349611 0.0 0.524141 1.289585 0.384034 0.217611 0.161906 NaN NaN
1 49.021621 13.406203 8.494789 8.91942 10.763932 3.410622 0.6347 0.251858 2.50811 0.199926 0.0 0.692912 1.483525 0.0 0.163488 0.048894 NaN NaN
2 48.247397 15.670665 5.687451 10.352983 10.144835 3.799917 0.579199 0.116463 2.703525 0.313039 0.0 0.678262 1.34723 0.029433 0.053967 0.275634 NaN NaN
3 48.691004 13.595075 7.825504 10.155966 10.394342 3.6004 0.633725 0.139596 2.580742 0.192747 0.0 0.491091 1.477336 0.145906 0.064359 0.012207 NaN NaN
4 49.668286 14.139461 5.626875 10.317219 10.487224 3.856336 0.609107 0.118617 2.712715 0.264858 0.0 0.482999 1.367683 0.159937 0.188682 0.0 NaN NaN