<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 <h1>The BurnMan Tutorial</h1>
Part 4
Step1:  Fitting parameters for an equation of state to experimental data
BurnMan contains least-squares optimization functions that fit model parameters to data.
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. 
An extended example of fitting various kinds of data, outlier removal and detailed analysis can be found in examples/example_fit_eos.py.
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)].
Step2:  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.
The fitting function also takes covariance matrices as input, so we have to build those matrices.
Step3:  The next code block creates a Mineral object (stv), providing starting guesses for the parameters.
The user selects which parameters they wish to fit, and which they wish to keep fixed.
The parameters of the Mineral object are automatically updated during fitting.
Finally, the optimized parameter values and their variances are printed to screen.
Step4:  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.
Step5:  We now plot our optimized equation of state against the original data.
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.
Step6:  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.
Step7:  Finding the best fit endmember proportions of a solution given a bulk composition
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.
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)
Step8:  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).
Step9:  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.
Step10:  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.
Step11:  Finally, the optimized parameters can be used to set the composition of the solution model, and the optimized parameters printed to stdout.
Step12:  As in the equation of state fitting, a corner plot of the covariances can also be plotted.
Step13:  Fitting phase proportions to a bulk composition
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.
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.
First, some tedious data preparation...
Step14:  The following code block loops over each of the compositions, and finds the best weight proportions and uncertainties on those proportions.
Step15:  Finally, we plot the data.
<ASSISTANT_TASK:>
Python Code:
import burnman
import numpy as np
import matplotlib.pyplot as plt
Explanation: <h1>The BurnMan Tutorial</h1>
Part 4: Fitting
This file is part of BurnMan - a thermoelastic and thermodynamic toolkit
for the Earth and Planetary Sciences
Copyright (C) 2012 - 2021 by the BurnMan team,
released under the GNU GPL v2 or later.
Introduction
This ipython notebook is the fourth in a series designed to introduce new users to the code structure and functionalities present in BurnMan.
<b>Demonstrates</b>
burnman.optimize.eos_fitting.fit_PTV_data
burnman.optimize.composition_fitting.fit_composition_to_solution
burnman.optimize.composition_fitting.fit_phase_proportions_to_bulk_composition
Everything in BurnMan and in this tutorial is defined in SI units.
End of explanation
PV = np.array([[0.0001, 46.5126, 0.0061],
               [1.168, 46.3429, 0.0053],
               [2.299, 46.1756, 0.0043],
               [3.137, 46.0550, 0.0051],
               [4.252, 45.8969, 0.0045],
               [5.037, 45.7902, 0.0053],
               [5.851, 45.6721, 0.0038],
               [6.613, 45.5715, 0.0050],
               [7.504, 45.4536, 0.0041],
               [8.264, 45.3609, 0.0056],
               [9.635, 45.1885, 0.0042],
               [11.69, 44.947, 0.002],                
               [17.67, 44.264, 0.002],
               [22.38, 43.776, 0.003],
               [29.38, 43.073, 0.009],
               [37.71, 42.278, 0.008],
               [46.03, 41.544, 0.017],
               [52.73, 40.999, 0.009],
               [26.32, 43.164, 0.006],
               [30.98, 42.772, 0.005],
               [34.21, 42.407, 0.003],
               [38.45, 42.093, 0.004],
               [43.37, 41.610, 0.004],
               [47.49, 41.280, 0.007]])
print(f'{len(PV)} data points loaded successfully.')
Explanation: Fitting parameters for an equation of state to experimental data
BurnMan contains least-squares optimization functions that fit model parameters to data.
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. 
An extended example of fitting various kinds of data, outlier removal and detailed analysis can be found in examples/example_fit_eos.py.
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)].
End of explanation
from burnman.tools.unitcell import molar_volume_from_unit_cell_volume
Z = 2. # number of formula units per unit cell in stishovite
PTV = np.array([PV[:,0]*1.e9,
                298.15 * np.ones_like(PV[:,0]),
                molar_volume_from_unit_cell_volume(PV[:,1], Z)]).T
# Here, we assume that the pressure uncertainties are equal to 3% of the total pressure, 
# that the temperature uncertainties are negligible, and take the unit cell volume
# uncertainties from the paper.
# We also assume that the uncertainties in pressure and volume are uncorrelated.
nul = np.zeros_like(PTV[:,0])
PTV_covariances = np.array([[0.03*PTV[:,0], nul, nul],
                            [nul, nul, nul],
                            [nul, nul, molar_volume_from_unit_cell_volume(PV[:,2], Z)]]).T
PTV_covariances = np.power(PTV_covariances, 2.)
Explanation: 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.
The fitting function also takes covariance matrices as input, so we have to build those matrices.
End of explanation
stv = burnman.minerals.HP_2011_ds62.stv()
params = ['V_0', 'K_0', 'Kprime_0']
fitted_eos = burnman.optimize.eos_fitting.fit_PTV_data(stv, params, PTV, PTV_covariances, verbose=False)
print('Optimized equation of state for stishovite:')
burnman.tools.misc.pretty_print_values(fitted_eos.popt, fitted_eos.pcov, fitted_eos.fit_params)
print('')
Explanation: The next code block creates a Mineral object (stv), providing starting guesses for the parameters.
The user selects which parameters they wish to fit, and which they wish to keep fixed.
The parameters of the Mineral object are automatically updated during fitting.
Finally, the optimized parameter values and their variances are printed to screen.
End of explanation
import matplotlib
matplotlib.rc('axes.formatter', useoffset=False) # turns offset off, makes for a more readable plot
fig = burnman.nonlinear_fitting.corner_plot(fitted_eos.popt, fitted_eos.pcov,
                                            params)
axes = fig[1]
axes[1][0].set_xlabel('$V_0$ ($\\times 10^{-5}$ m$^3$)')
axes[1][1].set_xlabel('$K_0$ ($\\times 10^{11}$ Pa)')
axes[0][0].set_ylabel('$K_0$ ($\\times 10^{11}$ Pa)')
axes[1][0].set_ylabel('$K\'_0$')
plt.show()
Explanation: 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.
End of explanation
from burnman.tools.misc import attribute_function
from burnman.optimize.nonlinear_fitting import confidence_prediction_bands
T = 298.15
pressures = np.linspace(1.e5, 60.e9, 101)
temperatures = T*np.ones_like(pressures)
volumes = stv.evaluate(['V'], pressures, temperatures)[0]
PTVs = np.array([pressures, temperatures, volumes]).T
# Calculate the 95% confidence and prediction bands
cp_bands = confidence_prediction_bands(model=fitted_eos,
                                       x_array=PTVs,
                                       confidence_interval=0.95,
                                       f=attribute_function(stv, 'V'),
                                       flag='V')
plt.fill_between(pressures/1.e9, cp_bands[2] * 1.e6, cp_bands[3] * 1.e6,
                 color=[0.75, 0.25, 0.55], label='95% prediction bands')
plt.fill_between(pressures/1.e9, cp_bands[0] * 1.e6, cp_bands[1] * 1.e6,
                 color=[0.75, 0.95, 0.95], label='95% confidence bands')
plt.plot(PTVs[:,0] / 1.e9, PTVs[:,2] * 1.e6, label='Optimized fit for stishovite')
plt.errorbar(PTV[:,0] / 1.e9, PTV[:,2] * 1.e6,
             xerr=np.sqrt(PTV_covariances[:,0,0]) / 1.e9,
             yerr=np.sqrt(PTV_covariances[:,2,2]) * 1.e6,
             linestyle='None', marker='o', label='Andrault et al. (2003)')
plt.ylabel("Volume (cm$^3$/mol)")
plt.xlabel("Pressure (GPa)")
plt.legend(loc="upper right")
plt.title("Stishovite EoS (room temperature)")
plt.show()
Explanation: We now plot our optimized equation of state against the original data.
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.
End of explanation
cp_bands = confidence_prediction_bands(model=fitted_eos,
                                       x_array=PTVs,
                                       confidence_interval=0.95,
                                       f=attribute_function(stv, 'K_T'),
                                       flag='V')
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')
plt.plot(pressures/1.e9, (cp_bands[0] + cp_bands[1])/2.e9, color='b', label='Best fit')
plt.ylabel("Bulk modulus (GPa)")
plt.xlabel("Pressure (GPa)")
plt.legend(loc="upper right")
plt.title("Stishovite EoS; uncertainty in bulk modulus (room temperature)")
plt.show()
Explanation: 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.
End of explanation
from burnman import minerals
gt = minerals.JH_2015.garnet()
print(f'Endmembers: {gt.endmember_names}')
print(f'Elements: {gt.elements}')
print('Stoichiometric matrix:')
print(gt.stoichiometric_matrix)
Explanation: Finding the best fit endmember proportions of a solution given a bulk composition
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.
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):
End of explanation
fitted_variables = ['Fe', 'Ca', 'Mg', 'Cr', 'Al', 'Si', 'Fe3+']
variable_values = np.array([1.1, 2., 0., 0, 1.9, 3., 0.1])
variable_covariances = np.eye(7)*0.01*0.01
# Add some noise.
v_err = np.random.rand(7)
np.random.seed(100)
variable_values = np.random.multivariate_normal(variable_values,
                                                variable_covariances)
Explanation: 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).
End of explanation
variable_conversions = {'Fe3+': {'Fef_B': 1.}}
Explanation: 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.
End of explanation
from burnman.optimize.composition_fitting import fit_composition_to_solution
popt, pcov, res = fit_composition_to_solution(gt,
                                              fitted_variables,
                                              variable_values,
                                              variable_covariances,
                                              variable_conversions)
Explanation: 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.
End of explanation
# We can set the composition of gt using the optimized parameters
gt.set_composition(popt)
# Print the optimized parameters and principal uncertainties
print('Molar fractions:')
for i in range(len(popt)):
    print(f'{gt.endmember_names[i]}: '
            f'{gt.molar_fractions[i]:.3f} +/- '
            f'{np.sqrt(pcov[i][i]):.3f}')
print(f'Weighted residual: {res:.3f}')
Explanation: Finally, the optimized parameters can be used to set the composition of the solution model, and the optimized parameters printed to stdout.
End of explanation
fig = burnman.nonlinear_fitting.corner_plot(popt, pcov, gt.endmember_names)
Explanation: As in the equation of state fitting, a corner plot of the covariances can also be plotted.
End of explanation
import itertools
# Load and transpose input data
filename = '../burnman/data/input_fitting/Bertka_Fei_1997_mars_mantle.dat'
with open(filename) as f:
    column_names = f.readline().strip().split()[1:]
data = np.genfromtxt(filename, dtype=None, encoding='utf8')
data = list(map(list, itertools.zip_longest(*data, fillvalue=None)))
# The first six columns are compositions given in weight % oxides
compositions = np.array(data[:6])
# The first row is the bulk composition
bulk_composition = compositions[:, 0]
# Load all the data into a dictionary
data = {column_names[i]: np.array(data[i])
        for i in range(len(column_names))}
# Make ordered lists of samples (i.e. experiment ID) and phases
samples = []
phases = []
for i in range(len(data['sample'])):
    if data['sample'][i] not in samples:
        samples.append(data['sample'][i])
    if data['phase'][i] not in phases:
        phases.append(data['phase'][i])
samples.remove("bulk_composition")
phases.remove("bulk")
# Get the indices of all the phases present in each sample
sample_indices = [[i for i in range(len(data['sample']))
                    if data['sample'][i] == sample]
                    for sample in samples]
# Get the run pressures of each experiment
pressures = np.array([data['pressure'][indices[0]] for indices in sample_indices])
Explanation: Fitting phase proportions to a bulk composition
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.
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.
First, some tedious data preparation...
End of explanation
from burnman.optimize.composition_fitting import fit_phase_proportions_to_bulk_composition
# Create empty arrays to store the weight proportions of each phase,
# and the principal uncertainties (we do not use the covariances here,
# although they are calculated)
weight_proportions = np.zeros((len(samples), len(phases)))*np.NaN
weight_proportion_uncertainties = np.zeros((len(samples),
                                            len(phases)))*np.NaN
residuals = []
# Loop over the samples, fitting phase proportions
# to the provided bulk composition
for i, sample in enumerate(samples):
    # This line does the heavy lifting
    popt, pcov, res = fit_phase_proportions_to_bulk_composition(compositions[:, sample_indices[i]],
                                                                bulk_composition)
    residuals.append(res)
    # Fill the correct elements of the weight_proportions
    # and weight_proportion_uncertainties arrays
    sample_phases = [data['phase'][i] for i in sample_indices[i]]
    for j, phase in enumerate(sample_phases):
        weight_proportions[i, phases.index(phase)] = popt[j]
        weight_proportion_uncertainties[i, phases.index(phase)] = np.sqrt(pcov[j][j])
Explanation: The following code block loops over each of the compositions, and finds the best weight proportions and uncertainties on those proportions.
End of explanation
fig = plt.figure(figsize=(6, 6))
ax = [fig.add_subplot(4, 1, 1)]
ax.append(fig.add_subplot(4, 1, (2, 4)))
for i, phase in enumerate(phases):
    ebar = plt.errorbar(pressures, weight_proportions[:, i],
                        yerr=weight_proportion_uncertainties[:, i],
                        fmt="none", zorder=2)
    ax[1].scatter(pressures, weight_proportions[:, i], label=phase, zorder=3)
ax[0].set_title('Phase proportions in the Martian Mantle (Bertka and Fei, 1997)')
ax[0].scatter(pressures, residuals)
for i in range(2):
    ax[i].set_xlim(0., 40.)
ax[1].set_ylim(0., 1.)
ax[0].set_ylabel('Residual')
ax[1].set_xlabel('Pressure (GPa)')
ax[1].set_ylabel('Phase fraction (wt %)')
ax[1].legend()
fig.set_tight_layout(True)
plt.show()
Explanation: Finally, we plot the data.
End of explanation
<END_TASK>