Source code for burnman.optimize.composition_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.

from __future__ import absolute_import

import numpy as np
from .linear_fitting import weighted_constrained_least_squares

[docs]def fit_composition_to_solution(solution, fitted_variables, variable_values, variable_covariances, variable_conversions=None, normalize=True): """ Takes a Solution object and a set of variable names and associates values and covariances and finds the molar fractions of the solution which provide the best fit (in a least-squares sense) to the variable values. The fitting applies appropriate non-negativity constraints (i.e. no species can have a negative occupancy on a site). Parameters ---------- solution : burnman.Solution object The solution to use in the fitting procedure. fitted_variables : list of strings A list of the variables used to find the best-fit molar fractions of the solution. These should either be elements such as "Fe", site_species such as "Fef_B" which would correspond to a species labelled Fef on the second site, or user-defined variables which are arithmetic sums of elements and/or site_species defined in "variable_conversions". variable_values : numpy array Numerical values of the fitted variables. These should be given as amounts; they do not need to be normalized. variable_covariances : 2D numpy array Covariance matrix of the variables. variable_conversions : dictionary of dictionaries or None A dictionary converting any user-defined variables into an arithmetic sum of element and site-species amounts. For example, {'Mg_equal': {'Mg_A': 1., 'Mg_B': -1.}}, coupled with Mg_equal = 0 would impose a constraint that the amount of Mg would be equal on the first and second site in the solution. normalize : boolean (default: True) If True, normalizes the optimized molar fractions to sum to unity. Returns ------- popt : numpy array Optimized molar fractions. pcov : 2D numpy array Covariance matrix corresponding to the optimized molar fractions. res : float The weighted residual of the fitting procedure. """ n_vars = len(fitted_variables) n_mbrs = len(solution.endmembers) solution_variables = solution.elements solution_variables.extend(solution.solution_model.site_names) solution_matrix = np.hstack((solution.stoichiometric_matrix, solution.solution_model.endmember_noccupancies)) n_sol_vars = solution_matrix.shape[1] if variable_conversions is not None: solution_matrix = np.hstack((solution_matrix, np.zeros((solution_matrix.shape[0], len(variable_conversions))))) for i, (new_var, conversion_dict) in enumerate(variable_conversions.items()): assert (new_var not in solution_variables) solution_variables.append(new_var) for var in conversion_dict.keys(): solution_matrix[:, n_sol_vars+i] += solution_matrix[:, solution_variables.index(var)] # Now, construct A using the fitted variables A = np.zeros((n_vars, solution_matrix.shape[0])) for i, var in enumerate(fitted_variables): A[i, :] = solution_matrix[:, solution_variables.index(var)] b = variable_values Cov_b = variable_covariances # Define the constraints # Ensure that element abundances / site occupancies # are exactly equal to zero if the user specifies that # they are equal to zero. S, S_index = np.unique(A, axis=0, return_index=True) S = np.array([s for i, s in enumerate(S) if np.abs(b[S_index[i]]) < 1.e-10 and any(np.abs(s) > 1.e-10)]) equality_constraints = [S, np.zeros(len(S))] # Ensure all site occupancies are non-negative T = np.array([-t for t in np.unique(solution.solution_model.endmember_occupancies.T, axis=0) if any(np.abs(t) > 1.e-10)]) inequality_constraints = [T, np.zeros(len(T))] popt, pcov, res = weighted_constrained_least_squares(A, b, Cov_b, equality_constraints, inequality_constraints) if normalize: sump = sum(popt) popt /= sump pcov /= sump * sump res /= sump # Convert the variance-covariance matrix from endmember amounts to # endmember proportions dpdx = (np.eye(n_mbrs) - popt).T # = (1. - p[i] if i == j else -p[i]) pcov = return (popt, pcov, res)
[docs]def fit_phase_proportions_to_bulk_composition(phase_compositions, bulk_composition): """ Performs weighted constrained least squares on a set of phase compositions to find the amount of those phases that best-fits a given bulk composition. The fitting applies appropriate non-negativity constraints (i.e. no phase can have a negative abundance in the bulk). Parameters ---------- phase_compositions : 2D numpy array The composition of each phase. Can be in weight or mole amounts. bulk_composition : numpy array The bulk composition of the composite. Must be in the same units as the phase compositions. Returns ------- popt : numpy array Optimized phase amounts. pcov : 2D numpy array Covariance matrix corresponding to the optimized phase amounts. res : float The weighted residual of the fitting procedure. """ n_phases = len(phase_compositions[0]) inequality_constraints = [-np.eye(n_phases), np.zeros(n_phases)] popt, pcov, res = weighted_constrained_least_squares(phase_compositions, bulk_composition, None, None, inequality_constraints) return (popt, pcov, res)