# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit for
# the Earth and Planetary Sciences
# Copyright (C) 2012 - 2025 by the BurnMan team, released under the GNU
# GPL v2 or later.
import numpy as np
from copy import deepcopy
from .linear_fitting import weighted_constrained_least_squares
from ..utils.chemistry import dictionarize_formula
from ..utils.chemistry import process_solution_chemistry
from ..classes.solution import Solution
[docs]
class DummyCompositionSolution(Solution):
"""
A mock Solution class to enable compositional fitting when no thermodynamic
model is available. Useful for pure composition analysis.
The model is initialized with appropriate
chemical formulae for each endmember, and can do all basic
compositional processing that doesn't involve any material
properties.
:param endmember_element_formulae: Formulae for each of the independent endmembers.
e.g. ['Mg2SiO4', 'Fe2SiO4'].
:type endmember_element_formulae: list of str
:param endmember_site_formulae: Site formulae for each of the independent endmembers,
in the same order as endmember_element_formulae. e.g. ['[Mg]2SiO4', '[Fe]2SiO4'].
:type endmember_site_formulae: list of str
"""
def __init__(self, endmember_element_formulae, endmember_site_formulae):
self.endmember_formulae = [
dictionarize_formula(f) for f in endmember_element_formulae
]
self.solution_model = type(
"Dimension", (object,), {"site_formulae": endmember_site_formulae}
)()
self.solution_model.endmembers = [None for f in endmember_site_formulae]
process_solution_chemistry(self.solution_model)
[docs]
def fit_composition_to_solution(
solution,
fitted_variables,
variable_values,
variable_covariances,
variable_conversions=None,
normalize=True,
):
"""
Fit a set of molar fractions to a solution model that best matches
observed compositional variables (in a least-squares sense).
The fitting applies appropriate non-negativity constraints
(i.e. no species can have a negative occupancy on a site).
:param solution: The solution to use in the fitting procedure.
:type solution: :class:`burnman.Solution`
:param fitted_variables: 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".
:type fitted_variables: list of str
:param variable_values: Numerical values of the fitted variables.
These should be given as amounts; they do not need to be normalized.
:type variable_values: numpy.array
:param variable_covariances: Covariance matrix of the variables.
:type variable_covariances: 2D numpy.array
:param variable_conversions: 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.
:type variable_conversions: dict of dict, or None
:param normalize: Whether to normalize the molar fractions to sum to one.
:type normalize: bool, default=True
:returns: (molar_fractions, covariance_matrix, weighted_residual)
- molar_fractions: Optimized molar fractions of endmembers.
- covariance_matrix: Covariance matrix of the optimized fractions,
computed from the pseudo-inverse of the weighted design matrix.
- weighted_residual: The weighted residual of the fit, calculated as the
square root of the weighted sum of squared differences
between observed and modeled variables.
:rtype: tuple of 1D numpy.array, 2D numpy.array and float
.. seealso::
:func:`burnman.optimize.linear_fitting.weighted_constrained_least_squares` -
Performs the weighted constrained least squares optimization used internally.
"""
solution_variables = deepcopy(solution.elements)
solution_variables.extend(solution.solution_model.site_names)
solution_matrix = np.hstack(
(solution.stoichiometric_matrix, solution.solution_model.endmember_noccupancies)
)
n_vars = len(fitted_variables)
n_mbrs, n_sol_vars = solution_matrix.shape
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.0e-10 and any(np.abs(s) > 1.0e-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.0e-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 = dpdx.dot(pcov).dot(dpdx.T)
return (popt, pcov, res)
[docs]
def fit_phase_proportions_to_bulk_composition(phase_compositions, bulk_composition):
"""
Determine the proportions of phases that best fit a given bulk composition
(in a weighted constrained least squares sense).
The fitting applies appropriate non-negativity constraints
(i.e. no phase can have a negative abundance in the bulk).
:param phase_compositions: The composition of each phase.
Can be in weight or mole amounts.
:type phase_compositions: 2D numpy.array
:param bulk_composition: The bulk composition of the composite.
Must be in the same units as the phase compositions.
:type bulk_composition: numpy.array
:returns: (phase_proportions, covariance_matrix, weighted_residual)
- phase_proportions: Optimized proportions of each phase.
- covariance_matrix: Covariance matrix of the optimized phase proportions,
computed from the pseudo-inverse of the weighted design matrix.
- weighted_residual: The weighted residual of the fit, representing the square root
of the weighted sum of squared differences between observed
and modeled bulk compositions.
:rtype: tuple of 1D numpy.array, 2D numpy.array and float
.. seealso::
:func:`burnman.optimize.linear_fitting.weighted_constrained_least_squares` -
Performs the weighted constrained least squares optimization used internally.
"""
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)