{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "C30Xe311ePd3"
},
"source": [
"The BurnMan Tutorial
\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Part 4: Fitting\n",
"This file is part of BurnMan - a thermoelastic and thermodynamic toolkit\n",
"for the Earth and Planetary Sciences\n",
"\n",
"Copyright (C) 2012 - 2021 by the BurnMan team,\n",
"released under the GNU GPL v2 or later."
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "13QqqnL6UzBk"
},
"source": [
"### Introduction\n",
"\n",
"This ipython notebook is the fourth in a series designed to introduce new users to the code structure and functionalities present in BurnMan.\n",
"\n",
"Demonstrates\n",
"\n",
"1. burnman.optimize.eos_fitting.fit_PTV_data\n",
"2. burnman.optimize.composition_fitting.fit_composition_to_solution\n",
"3. burnman.optimize.composition_fitting.fit_phase_proportions_to_bulk_composition\n",
"\n",
"Everything in BurnMan and in this tutorial is defined in SI units. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import burnman\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "XkGTkGDwNppm"
},
"source": [
"### Fitting parameters for an equation of state to experimental data\n",
"\n",
"BurnMan contains least-squares optimization functions that fit model parameters to data.\n",
"There are two helper functions especially for use in fitting Mineral parameters to experimental data; these are `burnman.optimize.eos_fitting.fit_PTp_data` (which can fit multiple kinds of data at the same time), and `burnman.optimize.eos_fitting.fit_PTV_data`, which specifically fits only pressure-temperature-volume data. \n",
"\n",
"An extended example of fitting various kinds of data, outlier removal and detailed analysis can be found in `examples/example_fit_eos.py`.\n",
"In this tutorial, we shall focus solely on fitting room temperature pressure-temperature-volume data. Specifically, the data we will fit is experimental volumes of stishovite, taken from Andrault et al. (2003). This data is provided in the form [P (GPa), V (Angstrom^3) and sigma_V (Angstrom^3)]."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 1000
},
"id": "Fyh8uDFdNr0D",
"outputId": "65478de5-262a-4b95-9036-46817f8331f6"
},
"outputs": [],
"source": [
"PV = np.array([[0.0001, 46.5126, 0.0061],\n",
" [1.168, 46.3429, 0.0053],\n",
" [2.299, 46.1756, 0.0043],\n",
" [3.137, 46.0550, 0.0051],\n",
" [4.252, 45.8969, 0.0045],\n",
" [5.037, 45.7902, 0.0053],\n",
" [5.851, 45.6721, 0.0038],\n",
" [6.613, 45.5715, 0.0050],\n",
" [7.504, 45.4536, 0.0041],\n",
" [8.264, 45.3609, 0.0056],\n",
" [9.635, 45.1885, 0.0042],\n",
" [11.69, 44.947, 0.002], \n",
" [17.67, 44.264, 0.002],\n",
" [22.38, 43.776, 0.003],\n",
" [29.38, 43.073, 0.009],\n",
" [37.71, 42.278, 0.008],\n",
" [46.03, 41.544, 0.017],\n",
" [52.73, 40.999, 0.009],\n",
" [26.32, 43.164, 0.006],\n",
" [30.98, 42.772, 0.005],\n",
" [34.21, 42.407, 0.003],\n",
" [38.45, 42.093, 0.004],\n",
" [43.37, 41.610, 0.004],\n",
" [47.49, 41.280, 0.007]])\n",
"\n",
"print(f'{len(PV)} data points loaded successfully.')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"BurnMan works exclusively in SI units, so we have to convert from GPa to Pa, and Angstrom per cell into molar volume in m^3.\n",
"The fitting function also takes covariance matrices as input, so we have to build those matrices."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from burnman.tools.unitcell import molar_volume_from_unit_cell_volume\n",
"\n",
"Z = 2. # number of formula units per unit cell in stishovite\n",
"PTV = np.array([PV[:,0]*1.e9,\n",
" 298.15 * np.ones_like(PV[:,0]),\n",
" molar_volume_from_unit_cell_volume(PV[:,1], Z)]).T\n",
"\n",
"# Here, we assume that the pressure uncertainties are equal to 3% of the total pressure, \n",
"# that the temperature uncertainties are negligible, and take the unit cell volume\n",
"# uncertainties from the paper.\n",
"# We also assume that the uncertainties in pressure and volume are uncorrelated.\n",
"nul = np.zeros_like(PTV[:,0])\n",
"PTV_covariances = np.array([[0.03*PTV[:,0], nul, nul],\n",
" [nul, nul, nul],\n",
" [nul, nul, molar_volume_from_unit_cell_volume(PV[:,2], Z)]]).T\n",
"PTV_covariances = np.power(PTV_covariances, 2.)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The next code block creates a Mineral object (`stv`), providing starting guesses for the parameters.\n",
"The user selects which parameters they wish to fit, and which they wish to keep fixed.\n",
"The parameters of the Mineral object are automatically updated during fitting.\n",
"\n",
"Finally, the optimized parameter values and their variances are printed to screen."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"stv = burnman.minerals.HP_2011_ds62.stv()\n",
"params = ['V_0', 'K_0', 'Kprime_0']\n",
"fitted_eos = burnman.optimize.eos_fitting.fit_PTV_data(stv, params, PTV, PTV_covariances, verbose=False)\n",
"\n",
"print('Optimized equation of state for stishovite:')\n",
"burnman.utils.misc.pretty_print_values(fitted_eos.popt, fitted_eos.pcov, fitted_eos.fit_params)\n",
"print('')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `fitted_eos` object contains a lot of useful information about the fit. In the next code block, we fit the corner plot of the covariances, showing the tradeoffs in different parameters."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib\n",
"matplotlib.rc('axes.formatter', useoffset=False) # turns offset off, makes for a more readable plot\n",
"\n",
"fig = burnman.nonlinear_fitting.corner_plot(fitted_eos.popt, fitted_eos.pcov,\n",
" params)\n",
"axes = fig[1]\n",
"axes[1][0].set_xlabel('$V_0$ ($\\\\times 10^{-5}$ m$^3$)')\n",
"axes[1][1].set_xlabel('$K_0$ ($\\\\times 10^{11}$ Pa)')\n",
"axes[0][0].set_ylabel('$K_0$ ($\\\\times 10^{11}$ Pa)')\n",
"axes[1][0].set_ylabel('$K\\'_0$')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We now plot our optimized equation of state against the original data.\n",
"BurnMan also includes a useful function `burnman.optimize.nonlinear_fitting.confidence_prediction_bands` that can be used to calculate the nth percentile confidence and prediction bounds on a function given a model using the delta method."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from burnman.utils.misc import attribute_function\n",
"from burnman.optimize.nonlinear_fitting import confidence_prediction_bands\n",
"\n",
"T = 298.15\n",
"pressures = np.linspace(1.e5, 60.e9, 101)\n",
"temperatures = T*np.ones_like(pressures)\n",
"volumes = stv.evaluate(['V'], pressures, temperatures)[0]\n",
"PTVs = np.array([pressures, temperatures, volumes]).T\n",
"\n",
"# Calculate the 95% confidence and prediction bands\n",
"cp_bands = confidence_prediction_bands(model=fitted_eos,\n",
" x_array=PTVs,\n",
" confidence_interval=0.95,\n",
" f=attribute_function(stv, 'V'),\n",
" flag='V')\n",
"\n",
"plt.fill_between(pressures/1.e9, cp_bands[2] * 1.e6, cp_bands[3] * 1.e6,\n",
" color=[0.75, 0.25, 0.55], label='95% prediction bands')\n",
"plt.fill_between(pressures/1.e9, cp_bands[0] * 1.e6, cp_bands[1] * 1.e6,\n",
" color=[0.75, 0.95, 0.95], label='95% confidence bands')\n",
"\n",
"plt.plot(PTVs[:,0] / 1.e9, PTVs[:,2] * 1.e6, label='Optimized fit for stishovite')\n",
"plt.errorbar(PTV[:,0] / 1.e9, PTV[:,2] * 1.e6,\n",
" xerr=np.sqrt(PTV_covariances[:,0,0]) / 1.e9,\n",
" yerr=np.sqrt(PTV_covariances[:,2,2]) * 1.e6,\n",
" linestyle='None', marker='o', label='Andrault et al. (2003)')\n",
"\n",
"plt.ylabel(\"Volume (cm$^3$/mol)\")\n",
"plt.xlabel(\"Pressure (GPa)\")\n",
"plt.legend(loc=\"upper right\")\n",
"plt.title(\"Stishovite EoS (room temperature)\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can also calculate the confidence and prediction bands for any other property of the mineral. In the code block below, we calculate and plot the optimized isothermal bulk modulus and its uncertainties."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"cp_bands = confidence_prediction_bands(model=fitted_eos,\n",
" x_array=PTVs,\n",
" confidence_interval=0.95,\n",
" f=attribute_function(stv, 'K_T'),\n",
" flag='V')\n",
"\n",
"plt.fill_between(pressures/1.e9, (cp_bands[0])/1.e9, (cp_bands[1])/1.e9, color=[0.75, 0.95, 0.95], label='95% confidence band')\n",
"plt.plot(pressures/1.e9, (cp_bands[0] + cp_bands[1])/2.e9, color='b', label='Best fit')\n",
"\n",
"plt.ylabel(\"Bulk modulus (GPa)\")\n",
"plt.xlabel(\"Pressure (GPa)\")\n",
"plt.legend(loc=\"upper right\")\n",
"plt.title(\"Stishovite EoS; uncertainty in bulk modulus (room temperature)\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "fow-vPnOSVvu"
},
"source": [
"### Finding the best fit endmember proportions of a solution given a bulk composition\n",
"\n",
"Let's now turn our focus to a different kind of fitting. It is common in petrology to have a bulk composition of a phase (provided, for example, by electron probe microanalysis), and want to turn this composition into a formula that satisfies stoichiometric constraints. This can be formulated as a constrained, weighted least squares problem, and BurnMan can be used to solve these problems using the function `burnman.optimize.composition_fitting.fit_composition_to_solution`.\n",
"\n",
"In the following example, we shall create a model garnet composition, and then fit that to the Jennings and Holland (2015) garnet solution model. First, let's look at the solution model endmembers (pyrope, almandine, grossular, andradite and knorringite):"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "g7C_nKvthteT",
"outputId": "267e6a2b-fbfd-4ac9-ea6e-4cc428e9ae34"
},
"outputs": [],
"source": [
"from burnman import minerals\n",
"\n",
"gt = minerals.JH_2015.garnet()\n",
"\n",
"print(f'Endmembers: {gt.endmember_names}')\n",
"print(f'Elements: {gt.elements}')\n",
"print('Stoichiometric matrix:')\n",
"print(gt.stoichiometric_matrix)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, let's create a model garnet composition. A unique composition can be determined with the species Fe (total), Ca, Mg, Cr, Al, Si and Fe3+, all given in mole amounts. On top of this, we add some random noise (using a fixed seed so that the composition is reproducible)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fitted_variables = ['Fe', 'Ca', 'Mg', 'Cr', 'Al', 'Si', 'Fe3+']\n",
"variable_values = np.array([1.1, 2., 0., 0, 1.9, 3., 0.1])\n",
"variable_covariances = np.eye(7)*0.01*0.01\n",
"\n",
"# Add some noise.\n",
"v_err = np.random.rand(7)\n",
"np.random.seed(100)\n",
"variable_values = np.random.multivariate_normal(variable_values,\n",
" variable_covariances)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Importantly, Fe3+ isn't an element or a site-species of the solution model, so we need to provide the linear conversion from Fe3+ to elements and/or site species. In this case, Fe3+ resides only on the second site (Site B), and the JH_2015.gt model has labelled Fe3+ on that site as Fef. Therefore, the conversion is simply Fe3+ = Fef_B."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"variable_conversions = {'Fe3+': {'Fef_B': 1.}}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we're ready to do the fitting. The following line is all that is required, and yields as output the optimized parameters, the corresponding covariance matrix and the residual."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from burnman.optimize.composition_fitting import fit_composition_to_solution\n",
"popt, pcov, res = fit_composition_to_solution(gt,\n",
" fitted_variables,\n",
" variable_values,\n",
" variable_covariances,\n",
" variable_conversions)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, the optimized parameters can be used to set the composition of the solution model, and the optimized parameters printed to stdout."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# We can set the composition of gt using the optimized parameters\n",
"gt.set_composition(popt)\n",
"\n",
"# Print the optimized parameters and principal uncertainties\n",
"print('Molar fractions:')\n",
"for i in range(len(popt)):\n",
" print(f'{gt.endmember_names[i]}: '\n",
" f'{gt.molar_fractions[i]:.3f} +/- '\n",
" f'{np.sqrt(pcov[i][i]):.3f}')\n",
"\n",
"print(f'Weighted residual: {res:.3f}')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As in the equation of state fitting, a corner plot of the covariances can also be plotted."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig = burnman.nonlinear_fitting.corner_plot(popt, pcov, gt.endmember_names)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Fitting phase proportions to a bulk composition\n",
"\n",
"Another common constrained weighted least squares problem involves fitting phase proportions, given their individual compositions and the overall bulk composition. This is particularly important in experimental petrology, where the bulk composition is known from a starting composition. In these cases, the residual after fitting is often used to assess whether the sample remained a closed system during the experiment.\n",
"\n",
"In the following example, we take phase compositions and the bulk composition reported from high pressure experiments on a Martian mantle composition by Bertka and Fei (1997), and use these to calculate phase proportions in Mars mantle, and the quality of the experiments.\n",
"\n",
"First, some tedious data preparation..."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import itertools\n",
"\n",
"# Load and transpose input data\n",
"filename = '../burnman/data/input_fitting/Bertka_Fei_1997_mars_mantle.dat'\n",
"with open(filename) as f:\n",
" column_names = f.readline().strip().split()[1:]\n",
"data = np.genfromtxt(filename, dtype=None, encoding='utf8')\n",
"data = list(map(list, itertools.zip_longest(*data, fillvalue=None)))\n",
"\n",
"# The first six columns are compositions given in weight % oxides\n",
"compositions = np.array(data[:6])\n",
"\n",
"# The first row is the bulk composition\n",
"bulk_composition = compositions[:, 0]\n",
"\n",
"# Load all the data into a dictionary\n",
"data = {column_names[i]: np.array(data[i])\n",
" for i in range(len(column_names))}\n",
"\n",
"# Make ordered lists of samples (i.e. experiment ID) and phases\n",
"samples = []\n",
"phases = []\n",
"for i in range(len(data['sample'])):\n",
" if data['sample'][i] not in samples:\n",
" samples.append(data['sample'][i])\n",
" if data['phase'][i] not in phases:\n",
" phases.append(data['phase'][i])\n",
"\n",
"samples.remove(\"bulk_composition\")\n",
"phases.remove(\"bulk\")\n",
"\n",
"# Get the indices of all the phases present in each sample\n",
"sample_indices = [[i for i in range(len(data['sample']))\n",
" if data['sample'][i] == sample]\n",
" for sample in samples]\n",
"\n",
"# Get the run pressures of each experiment\n",
"pressures = np.array([data['pressure'][indices[0]] for indices in sample_indices])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following code block loops over each of the compositions, and finds the best weight proportions and uncertainties on those proportions."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from burnman.optimize.composition_fitting import fit_phase_proportions_to_bulk_composition\n",
"\n",
"# Create empty arrays to store the weight proportions of each phase,\n",
"# and the principal uncertainties (we do not use the covariances here,\n",
"# although they are calculated)\n",
"weight_proportions = np.zeros((len(samples), len(phases)))*np.NaN\n",
"weight_proportion_uncertainties = np.zeros((len(samples),\n",
" len(phases)))*np.NaN\n",
"\n",
"residuals = []\n",
"# Loop over the samples, fitting phase proportions\n",
"# to the provided bulk composition\n",
"for i, sample in enumerate(samples):\n",
" # This line does the heavy lifting\n",
" popt, pcov, res = fit_phase_proportions_to_bulk_composition(compositions[:, sample_indices[i]],\n",
" bulk_composition)\n",
"\n",
" residuals.append(res)\n",
"\n",
" # Fill the correct elements of the weight_proportions\n",
" # and weight_proportion_uncertainties arrays\n",
" sample_phases = [data['phase'][i] for i in sample_indices[i]]\n",
" for j, phase in enumerate(sample_phases):\n",
" weight_proportions[i, phases.index(phase)] = popt[j]\n",
" weight_proportion_uncertainties[i, phases.index(phase)] = np.sqrt(pcov[j][j])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we plot the data."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig = plt.figure(figsize=(6, 5))\n",
"ax = [fig.add_subplot(3, 1, 1)]\n",
"ax.append(fig.add_subplot(3, 1, (2, 4)))\n",
"for i, phase in enumerate(phases):\n",
" ebar = plt.errorbar(pressures, weight_proportions[:, i],\n",
" yerr=weight_proportion_uncertainties[:, i],\n",
" fmt=\"none\", zorder=2)\n",
" ax[1].scatter(pressures, weight_proportions[:, i], label=phase, zorder=3)\n",
"\n",
"ax[0].set_title('Phase proportions in the Martian Mantle (Bertka and Fei, 1997)')\n",
"ax[0].scatter(pressures, residuals)\n",
"for i in range(2):\n",
" ax[i].set_xlim(0., 40.)\n",
"\n",
"ax[1].set_ylim(0., 1.)\n",
"ax[0].set_ylabel('Residual')\n",
"ax[1].set_xlabel('Pressure (GPa)')\n",
"ax[1].set_ylabel('Phase fraction (wt)')\n",
"ax[1].legend()\n",
"fig.set_tight_layout(True)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can see from this plot that most of the residuals are below one, indicating that the probe analyses consistent with the bulk composition. Three analyses have higher residuals, which may indicate a problem with the experiments, or with the analyses around the wadsleyite field. \n",
"\n",
"The phase proportions also show some nice trends; clinopyroxene weight percentage increases with pressure at the expense of orthopyroxene. Garnet / majorite percentage increases sharply as clinopyroxene is exhausted at 14-16 GPa.\n",
"\n",
"And we're done! Next time, we'll look at how to determine equilibrium assemblages in BurnMan."
]
}
],
"metadata": {
"colab": {
"collapsed_sections": [],
"name": "BurnMan_1.0_manuscript.ipynb",
"provenance": []
},
"kernelspec": {
"display_name": "Python 3.10.5 ('base')",
"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.5"
},
"vscode": {
"interpreter": {
"hash": "c6e4e9f98eb68ad3b7c296f83d20e6de614cb42e90992a65aa266555a3137d0d"
}
}
},
"nbformat": 4,
"nbformat_minor": 0
}