From 36739ed056e80e8a45b4d0a5be66ecbfa5ba6e36 Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Fri, 31 Oct 2025 00:13:34 +0000 Subject: [PATCH 1/4] added optimal PT calculations --- burnman/classes/composite.py | 34 ++++ burnman/optimize/__init__.py | 2 +- burnman/tools/__init__.py | 3 +- burnman/tools/polytope.py | 72 +++++++++ burnman/tools/thermobarometry.py | 261 +++++++++++++++++++++++++++++++ docs/tools.rst | 6 + tests/test_tools.py | 141 +++++++++++++++++ 7 files changed, 517 insertions(+), 2 deletions(-) create mode 100644 burnman/tools/thermobarometry.py diff --git a/burnman/classes/composite.py b/burnman/classes/composite.py index 54e5be8b9..ffb38b15f 100644 --- a/burnman/classes/composite.py +++ b/burnman/classes/composite.py @@ -473,6 +473,40 @@ def endmember_partial_gibbs(self): j += n_endmembers return partial_gibbs + @material_property + def endmember_partial_entropies(self): + """ + Returns the partial entropies for all + the endmember minerals in the Composite + """ + partial_entropies = np.empty(self.n_endmembers) + j = 0 + for i, n_endmembers in enumerate(self.endmembers_per_phase): + if n_endmembers == 1: + partial_entropies[j] = self.phases[i].molar_entropy + else: + partial_entropies[j : j + n_endmembers] = self.phases[ + i + ].partial_entropies + j += n_endmembers + return partial_entropies + + @material_property + def endmember_partial_volumes(self): + """ + Returns the partial volumes for all + the endmember minerals in the Composite + """ + partial_volumes = np.empty(self.n_endmembers) + j = 0 + for i, n_endmembers in enumerate(self.endmembers_per_phase): + if n_endmembers == 1: + partial_volumes[j] = self.phases[i].molar_volume + else: + partial_volumes[j : j + n_endmembers] = self.phases[i].partial_volumes + j += n_endmembers + return partial_volumes + @material_property def reaction_affinities(self): """ diff --git a/burnman/optimize/__init__.py b/burnman/optimize/__init__.py index 520981767..3ba5a370b 100644 --- a/burnman/optimize/__init__.py +++ b/burnman/optimize/__init__.py @@ -1,6 +1,6 @@ # 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 +# Copyright (C) 2012 - 2025 by the BurnMan team, released under the GNU # GPL v2 or later. """ diff --git a/burnman/tools/__init__.py b/burnman/tools/__init__.py index 9b56bd22b..873d84175 100644 --- a/burnman/tools/__init__.py +++ b/burnman/tools/__init__.py @@ -1,6 +1,6 @@ # 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 +# Copyright (C) 2012 - 2025 by the BurnMan team, released under the GNU # GPL v2 or later. @@ -23,4 +23,5 @@ from . import plot from . import polytope from . import solution +from . import thermobarometry from . import unitcell diff --git a/burnman/tools/polytope.py b/burnman/tools/polytope.py index d53c7dc67..0b86b65af 100644 --- a/burnman/tools/polytope.py +++ b/burnman/tools/polytope.py @@ -312,3 +312,75 @@ def simplify_composite_with_composition(composite, composition): return Composite(new_phases) else: return composite + + +def greedy_independent_endmember_selection( + endmember_site_occupancies, site_occupancies, small_fraction_tol=0.0, norm_tol=1e-12 +): + """ + Greedy algorithm to select independent endmembers from a solution to approximate given + site occupancies through a non-negative linear combination of endmember site occupancies. + + This function starts with the full site occupancies and then iteratively selects endmembers + to approximate those site occupancies. It loops through all possible endmembers in a number of steps. + At each step the algorithm selects the endmember that can be subtracted in the largest amount from the + current residual site occupancies without making any site occupancy negative. + The process continues until either no endmember can be subtracted in an amount greater than fraction_tol, + or the norm of the residual site occupancies is less than tol. + + :param endmember_site_occupancies: A 2D array of shape (m, n), where m is the number of endmembers + and n is the number of sites. Each row corresponds to the site occupancies of an endmember. + :type endmember_site_occupancies: np.ndarray + + :param site_occupancies: A 1D array of length n, representing the target site occupancies to approximate. + :type site_occupancies: np.ndarray + + :param small_fraction_tol: Algorithm stops if no endmember can be added with a molar fraction larger than this value. + :type small_fraction_tol: float, optional, default 0.0 + + :param norm_tol: Algorithm stops if the norm of the residual site occupancies is less than this value. + :type norm_tol: float, optional, default 1e-12 + + :return: indices of selected endmembers, their fractions, and the final residual site occupancies. + :rtype: tuple(list[int], list[float], np.ndarray) + """ + + site_occupancy_residuals = site_occupancies.copy() + indices = [] + fractions = [] + + for _ in range(endmember_site_occupancies.shape[0]): + # compute fraction for each candidate endmember + # fraction_i = min_j r_j / s_ij over s_ij > 0; if no s_ij>0 -> fraction_i = 0 + with np.errstate(divide="ignore", invalid="ignore"): + ratios = np.where( + endmember_site_occupancies > 0, + site_occupancy_residuals[np.newaxis, :] / endmember_site_occupancies, + np.inf, + ) # shape (m,n) + + # For each row, take minimum ratio over columns where S>0; if all entries inf -> set 0 + fractions_all = np.min(ratios, axis=1) + fractions_all[np.isinf(fractions_all)] = 0.0 + + # pick largest alpha + i_max = int(np.argmax(fractions_all)) + fraction_max = float(fractions_all[i_max]) + if fraction_max <= small_fraction_tol: + break # no further progress possible or desirable + + # update residual + site_occupancy_residuals = ( + site_occupancy_residuals - fraction_max * endmember_site_occupancies[i_max] + ) + + # ensure non-negativity in elements of the residual + site_occupancy_residuals[site_occupancy_residuals < 0] = 0.0 + indices.append(i_max) + fractions.append(fraction_max) + + # check if done + if np.linalg.norm(site_occupancy_residuals) <= norm_tol: + break + + return indices, fractions, site_occupancy_residuals diff --git a/burnman/tools/thermobarometry.py b/burnman/tools/thermobarometry.py new file mode 100644 index 000000000..b38cd8f73 --- /dev/null +++ b/burnman/tools/thermobarometry.py @@ -0,0 +1,261 @@ +import numpy as np +from scipy.linalg import sqrtm, block_diag +from scipy.optimize import minimize +import sympy as sp + +from ..classes.solution import Solution +from ..classes.combinedmineral import CombinedMineral + +from .polytope import ( + greedy_independent_endmember_selection, + solution_polytope_from_endmember_occupancies, +) + + +def _build_endmember_covariance_matrix(assemblage, dataset_covariances): + """ + Build the Gibbs free energy covariance matrix for the given assemblage. + + :param assemblage: The mineral assemblage for which to build the covariance matrix. + :type assemblage: Assemblage + + :param dataset_covariances: The covariance data from the thermodynamic dataset. + :type dataset_covariances: dict, with keys 'endmember_names' and 'covariance_matrix' + """ + # For endmembers that are made up of other endmembers (the "made" endmembers in HP-speak) + # we need to build a transformation matrix from the "raw" endmembers (those in the covariance matrix) + # to the "transformed" endmembers (those in our assemblage) + cov = dataset_covariances + + A = [] # transformation matrix to go from raw to transformed endmembers + + # Loop over phases in the assemblage to build the matrix A + for phase in assemblage.phases: + if isinstance(phase, Solution): + for endmember, _ in phase.endmembers: + if isinstance(endmember, CombinedMineral): + indices = [ + cov["endmember_names"].index(name) + for name in endmember.mixture.endmember_names + ] + b = endmember.mixture.molar_fractions + A.append([0] * len(cov["endmember_names"])) + for i, idx in enumerate(indices): + A[-1][idx] = b[i] + else: + A.append([0] * len(cov["endmember_names"])) + A[-1][cov["endmember_names"].index(endmember.name)] = 1.0 + else: # single endmember phase + A.append([0] * len(cov["endmember_names"])) + A[-1][cov["endmember_names"].index(phase.name)] = 1.0 + + # Convert A to a numpy array + A = np.array(A) + + # Transform the covariance matrix using matrix A + cov_transformed = A.dot(cov["covariance_matrix"]).dot(A.T) + eigenvalues, _ = np.linalg.eig(cov_transformed) + assert np.all( + eigenvalues > -1.0e-5 + ), f"Transformed covariance matrix has large negative eigenvalues: {min(eigenvalues)}!" + + return cov_transformed + + +def _build_solution_covariance_matrix(assemblage): + n_mbrs_all = len(assemblage.endmember_names) + cov_RTlna = np.zeros((n_mbrs_all, n_mbrs_all)) + + j = 0 + for phase in assemblage.phases: + if isinstance(phase, Solution): + H_G = phase.gibbs_hessian + n_mbrs_phase = H_G.shape[0] + cov_X_phase = np.array(phase.compositional_covariances) + + assert ( + cov_X_phase.ndim == 2 + ), "Compositional uncertainties must be provided as a covariance matrix." + assert ( + cov_X_phase.shape[0] == n_mbrs_phase + ), "Compositional uncertainties covariance matrix must have the same number of rows as phases." + assert ( + cov_X_phase.shape[1] == n_mbrs_phase + ), "Compositional uncertainties covariance matrix must be square." + + cov_RTlna_phase = H_G.dot(cov_X_phase).dot(H_G.T) + cov_RTlna[j : j + n_mbrs_phase, j : j + n_mbrs_phase] = cov_RTlna_phase + j += n_mbrs_phase + else: # Add formula if it is an endmember + cov_RTlna[j, j] = 1.0e-12 # small uncertainty for pure phases + j += 1 + return cov_RTlna + + +def estimate_conditions( + assemblage, + dataset_covariances, + guessed_conditions=np.array([1.0e9, 873.0]), + pressure_bounds=[1.0e5, 400.0e9], + temperature_bounds=[300.0, 4000.0], + P_scaling=1.0e6, + small_fraction_tol=0.0, + max_it=100, +): + """ + Perform the avPT thermobarometric inversion to find the optimal pressure and temperature + for a given mineral assemblage. Algorithm based on Powell and Holland (1994). + + :param assemblage: The mineral assemblage for which to perform the inversion. + :type assemblage: Assemblage + + :param dataset_covariances: The covariance data from the thermodynamic dataset. + :type dataset_covariances: dict, with keys 'endmember_names' and 'covariance_matrix' + + :param guessed_conditions: Initial guess for pressure (Pa) and temperature (K). + :type guessed_conditions: np.array of shape (2,) + + :param pressure_bounds: Bounds for pressure (Pa) during optimization. + :type pressure_bounds: list of two floats + + :param temperature_bounds: Bounds for temperature (K) during optimization. + :type temperature_bounds: list of two floats + + :param P_scaling: Scaling factor for pressure to improve numerical stability. + :type P_scaling: float + + :param small_fraction_tol: If > 0.0, reduces the number of endmembers in solution phases by + transforming to a smaller set of independent endmembers using a greedy algorithm + and excluding those with molar fractions smaller than this value during the inversion. + :type small_fraction_tol: float, optional, default 0.0 + + :param max_it: Maximum number of iterations for the optimization algorithm. + :type max_it: int, optional, default 100 + + :return: Result object from the optimization containing optimal P and T (x) and other properties + of the solution. These include the covariance matrix of the estimated parameters (xcov), + the correlation coefficient between P and T (xcorr), + the affinities of the independent reactions at the optimal conditions (affinities), + the covariance matrix of the affinities (acov), + the affinities weighted by the inverse square root of their covariance matrix (weighted_affinities), + the partial derivatives of the affinities with respect to P and T (dadx), + the number of independent reactions (n_reactions), the number of fitted parameters (n_params), + the degrees of freedom (degrees_of_freedom), the chi-squared value (fun), + the reduced chi-squared value (reduced_chisqr, + given by the reduced chi-squared divided by the number of degrees of freedom), + and the fit quality (fit, given by the square root of the reduced chi-squared value). + :rtype: OptimizeResult + """ + # TODO: Implement compositional unknown parameters like redox state + # to join P and T in the optimization + + P_fixed = np.isclose(pressure_bounds[0], pressure_bounds[1]) + T_fixed = np.isclose(temperature_bounds[0], temperature_bounds[1]) + if P_fixed and T_fixed: + raise Exception("Both pressure and temperature cannot be fixed!") + + # Build the reaction matrix R, possibly reducing the number of endmembers + # by excluding components with small molar fractions + transformation_matrix = [] + if small_fraction_tol > 0.0: + transformation_matrices = [] + for phase in assemblage.phases: + if isinstance(phase, Solution): + occs = phase.solution_model.endmember_occupancies + poly = solution_polytope_from_endmember_occupancies(occs) + S = poly.endmember_occupancies + T = poly.endmembers_as_independent_endmember_amounts + f = phase.molar_fractions + v = f.dot(occs) + mbrs = greedy_independent_endmember_selection( + S, v, small_fraction_tol=small_fraction_tol, norm_tol=1e-12 + ) + idx, _, _ = mbrs + transformation_matrices.append(T[idx, :]) + else: + transformation_matrices.append(np.array([[1.0]])) + + transformation_matrix = block_diag(*transformation_matrices) + else: + transformation_matrix = np.identity(len(assemblage.endmember_names)) + + # Here we build the reaction matrix R in the reduced endmember basis + # i.e. after applying the transformation matrix to the stoichiometric matrix + S = sp.Matrix(transformation_matrix.dot(assemblage.stoichiometric_array)) + S = S.applyfunc(lambda x: sp.nsimplify(x)) + R = np.array([v[:] for v in S.T.nullspace()], dtype=float) + R = R.dot(transformation_matrix) # transform back to original endmember basis + if len(R) == 0: + raise Exception("No independent reactions found for the given assemblage!") + + def calculate_cov_a(assemblage, theta, dataset_covariances): + assemblage.set_state(*theta) + cov1 = _build_endmember_covariance_matrix(assemblage, dataset_covariances) + cov2 = _build_solution_covariance_matrix(assemblage) + cov_mu = cov1 + cov2 + acov = R.dot(cov_mu).dot(R.T) + return acov + + def chisqr(args): + """ + Compute the objective misfit function (chi-squared) for given P and T. + :param theta: Array of pressure (Pa) and temperature (K). + :type theta: np.array of shape (2,) + + :return: Chi-squared misfit value. + :rtype: float + """ + # These are the chemical potential covariances for all endmembers + cov_a = calculate_cov_a( + assemblage, [args[0] * P_scaling, args[1]], dataset_covariances + ) + a = R.dot(assemblage.endmember_partial_gibbs) + chisqr = a.T.dot(np.linalg.inv(cov_a)).dot(a) + return chisqr + + # Minimize the chi-squared function + x0 = [guessed_conditions[0] / P_scaling, guessed_conditions[1]] + bounds = [ + (pressure_bounds[0] / P_scaling, pressure_bounds[1] / P_scaling), + (temperature_bounds[0], temperature_bounds[1]), + ] + res = minimize( + chisqr, + x0, + method="SLSQP", + bounds=bounds, + options={"maxiter": max_it}, + ) + + # Rescale pressure back to Pa + res.x[0] *= P_scaling + res.jac[0] /= P_scaling + + # Post-process results + # First we reset the assemblage to the optimal P and T + # and recalculate the affinity covariance matrix + res.acov = calculate_cov_a(assemblage, res.x, dataset_covariances) + + # Compute additional diagnostics + res.dadx = R.dot( + np.array( + [ + assemblage.endmember_partial_volumes, + -assemblage.endmember_partial_entropies, + ] + ).T + ) + + res.xcov = np.linalg.pinv(res.dadx.T.dot(np.linalg.pinv(res.acov)).dot(res.dadx)) + res.xcorr = res.xcov[0, 1] / np.sqrt(res.xcov[0, 0] * res.xcov[1, 1]) + + res.affinities = R.dot(assemblage.endmember_partial_gibbs) + res.weighted_affinities = np.linalg.pinv(sqrtm(res.acov)).dot(res.affinities) + + res.n_reactions = len(res.affinities) + res.n_params = len(res.x) + res.degrees_of_freedom = res.n_reactions - res.n_params + res.reduced_chisqr = res.fun / res.degrees_of_freedom + res.fit = np.sqrt(res.reduced_chisqr) + + return res diff --git a/docs/tools.rst b/docs/tools.rst index 53f274a54..d19647f17 100644 --- a/docs/tools.rst +++ b/docs/tools.rst @@ -31,6 +31,12 @@ Burnman has a number of high-level tools to help achieve common goals. Several of these have already been described in previous sections. +Optimal thermobarometry +----------------------- + +.. automodule:: burnman.tools.thermobarometry + + Plotting -------- diff --git a/tests/test_tools.py b/tests/test_tools.py index 3247483fa..53441dd3e 100644 --- a/tests/test_tools.py +++ b/tests/test_tools.py @@ -4,6 +4,7 @@ import matplotlib.pyplot as plt import burnman +from burnman import Composite from burnman.tools.chemistry import equilibrium_temperature from burnman.tools.chemistry import equilibrium_pressure from burnman.tools.chemistry import hugoniot @@ -17,6 +18,70 @@ from burnman.utils.math import _pad_ndarray_inverse_mirror from burnman.utils.math import is_positive_definite +from burnman.minerals import HP_2011_ds62, mp50MnNCKFMASHTO +from burnman.tools.thermobarometry import estimate_conditions + + +def setup_assemblage(scale=1.0): + """ + Set up a mineral assemblage for testing thermobarometry. + :param scale: Scale factor for compositional noise. + """ + np.random.seed(42) + + # Define observed mineral assemblage + mu = mp50MnNCKFMASHTO.mu() + bi = mp50MnNCKFMASHTO.bi() + g = mp50MnNCKFMASHTO.g() + ilmm = mp50MnNCKFMASHTO.ilmm() + st = mp50MnNCKFMASHTO.st() + q = HP_2011_ds62.q() + + # These compositions correspond to an equilibrated assemblage at P = 0.4 GPa and T = 873 K + # In practice, these would be the measured compositions of the minerals + compositions = [ + np.array( + [0.55510672, 0.00393874, 0.00365819, 0.43232719, 0.00339884, 0.00157032] + ), + np.array( + [ + 0.12024528, + 0.55795845, + -0.23022757, + 0.2928378, + 0.19457943, + 0.06260454, + 0.00200206, + ] + ), + np.array([0.11996805, 0.75282947, 0.10051702, 0.02013999, 0.00654547]), + np.array([-0.05461279, 0.84687573, 0.1062428, 0.05404236, 0.04745189]), + np.array([0.06925327, 0.80706976, 0.01346876, 0.03074615, 0.07946205]), + ] + + # Assumed compositional uncertainties (1 sigma) + mu.compositional_covariances = [0.01, 0.001, 0.001, 0.01, 0.001, 0.0005] + bi.compositional_covariances = [0.01, 0.01, 0.05, 0.01, 0.01, 0.01, 0.0005] + g.compositional_covariances = [0.01, 0.01, 0.01, 0.01, 0.001] + ilmm.compositional_covariances = [0.01, 0.01, 0.01, 0.01, 0.01] + st.compositional_covariances = [0.01, 0.01, 0.01, 0.01, 0.01] + + assemblage = Composite([mu, bi, g, ilmm, st, q], [1.0 / 6.0] * 6) + for phase, comp in zip(assemblage.phases[:-1], compositions): + noise = np.random.normal(0, 1.0, len(comp)) * phase.compositional_covariances + comp_with_noise = comp + noise * scale + # Normalize to ensure valid molar fractions + comp_with_noise /= np.sum(comp_with_noise) + phase.set_composition(comp_with_noise) + + # Convert uncertainties to covariances + for phase in assemblage.phases[:-1]: + phase.compositional_covariances = np.diag( + np.array(phase.compositional_covariances) ** 2 + ) + + return assemblage + class test_tools(BurnManTest): def test_eqm_T(self): @@ -325,6 +390,82 @@ def psi_func(f, Pth, params): self.assertTrue(len(ticks) == 11) self.assertTrue(len(lines) == 11) + def test_estimate_conditions(self): + assemblage = setup_assemblage() + theta = np.array([0.5e9, 500.0]) + dataset_covariances = HP_2011_ds62.cov() + res = estimate_conditions( + assemblage, dataset_covariances, guessed_conditions=theta + ) + self.assertTrue(res.success) + self.assertEqual(res.n_reactions, 18) + + def test_estimate_conditions_at_equilibrium(self): + assemblage = setup_assemblage(scale=0.0) + theta = np.array([0.5e9, 500.0]) + dataset_covariances = HP_2011_ds62.cov() + res = estimate_conditions( + assemblage, dataset_covariances, guessed_conditions=theta + ) + self.assertTrue(res.success) + self.assertArraysAlmostEqual( + [res.x[0] / 1.0e6, res.x[1]], [400.0, 873.0], tol=0.1 + ) + + def test_estimate_conditions_exclude_small_fractions(self): + assemblage = setup_assemblage() + theta = np.array([0.5e9, 500.0]) + dataset_covariances = HP_2011_ds62.cov() + res = estimate_conditions( + assemblage, + dataset_covariances, + guessed_conditions=theta, + small_fraction_tol=0.01, + ) + self.assertTrue(res.success) + self.assertEqual(res.n_reactions, 11) + + def test_estimate_conditions_fixed_P(self): + assemblage = setup_assemblage() + theta = np.array([0.5e9, 500.0]) + dataset_covariances = HP_2011_ds62.cov() + res = estimate_conditions( + assemblage, + dataset_covariances, + guessed_conditions=theta, + pressure_bounds=[1.0e9, 1.0e9], + ) + self.assertTrue(res.success) + self.assertEqual(res.x[0], 1.0e9) + + def test_estimate_conditions_fixed_T(self): + assemblage = setup_assemblage() + theta = np.array([0.5e9, 500.0]) + dataset_covariances = HP_2011_ds62.cov() + res = estimate_conditions( + assemblage, + dataset_covariances, + guessed_conditions=theta, + temperature_bounds=[500.0, 500.0], + ) + self.assertTrue(res.success) + self.assertEqual(res.x[1], 500.0) + + def test_estimate_conditions_fail_fixed_P_and_T(self): + assemblage = setup_assemblage() + theta = np.array([0.5e9, 500.0]) + dataset_covariances = HP_2011_ds62.cov() + + # assert that an exception is raised when both P and T are fixed + with self.assertRaises(Exception): + _ = estimate_conditions( + assemblage, + dataset_covariances, + guessed_conditions=theta, + pressure_bounds=[1.0e9, 1.0e9], + temperature_bounds=[500.0, 500.0], + ) + if __name__ == "__main__": unittest.main() From 3e88955cd2f9b8130075b674ce037f04d113211c Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Fri, 31 Oct 2025 16:58:08 +0000 Subject: [PATCH 2/4] improve weighted_constrained_least_squares --- burnman/optimize/linear_fitting.py | 39 ++++++++++---- tests/test_fitting.py | 86 ++++++++++++++++++++++++++++++ 2 files changed, 116 insertions(+), 9 deletions(-) diff --git a/burnman/optimize/linear_fitting.py b/burnman/optimize/linear_fitting.py index d3e00b7ff..4394dd13d 100644 --- a/burnman/optimize/linear_fitting.py +++ b/burnman/optimize/linear_fitting.py @@ -6,7 +6,7 @@ import importlib import numpy as np -from scipy.linalg import inv, sqrtm +from scipy.linalg import sqrtm import warnings try: @@ -18,7 +18,12 @@ def weighted_constrained_least_squares( - A, b, Cov_b=None, equality_constraints=None, inequality_constraints=None + A, + b, + Cov_b=None, + equality_constraints=None, + inequality_constraints=None, + allow_rank_deficient=False, ): """ Solves a weighted, constrained least squares problem using cvxpy. @@ -45,6 +50,10 @@ def weighted_constrained_least_squares( in the objective function above. :type inequality_constraints: list containing a 2D array and 1D array + :param allow_rank_deficient: If True, allows the problem to be solved + even if the design matrix is rank-deficient. + :type allow_rank_deficient: bool + :returns: Tuple containing the optimized phase amounts (1D numpy.array), a covariance matrix corresponding to the optimized phase amounts (2D numpy.array), and the weighted residual of the fitting procedure @@ -58,11 +67,23 @@ def weighted_constrained_least_squares( # Create the standard weighted least squares objective function # (https://stats.stackexchange.com/a/333551) n_vars = A.shape[1] - m = inv(sqrtm(Cov_b)) - mA = m @ A - mb = m @ b + M = np.linalg.pinv(sqrtm(Cov_b)) + MA = M @ A + Mb = M @ b x = cp.Variable(n_vars) - objective = cp.Minimize(cp.sum_squares(mA @ x - mb)) + objective = cp.Minimize(cp.sum_squares(MA @ x - Mb)) + + # Add a check for rank deficiency + rank_MA = np.linalg.matrix_rank(MA) + if not allow_rank_deficient and rank_MA < n_vars: + raise Exception( + f"The weighted design matrix is rank-deficient " + f"(Cov_b^(-1/2).A={rank_MA} < n_vars={n_vars}). " + "This probably means that you haven't supplied sufficient " + "independent constraints to yield a unique solution to the " + "problem. If you wish to proceed anyway, set " + "allow_rank_deficient=True in the function call." + ) constraints = [] if equality_constraints is not None: @@ -83,7 +104,7 @@ def weighted_constrained_least_squares( # Set up the problem and solve it warns = [] - if len(constraints) > 1: + if len(constraints) > 0: prob = cp.Problem(objective, constraints) else: prob = cp.Problem(objective) @@ -104,7 +125,7 @@ def weighted_constrained_least_squares( # Calculate the covariance matrix # (also from https://stats.stackexchange.com/a/333551) - inv_Cov_b = np.linalg.inv(Cov_b) - pcov = np.linalg.inv(A.T.dot(inv_Cov_b.dot(A))) + inv_Cov_b = np.linalg.pinv(Cov_b) + pcov = np.linalg.pinv(A.T.dot(inv_Cov_b.dot(A))) return (popt, pcov, res) diff --git a/tests/test_fitting.py b/tests/test_fitting.py index f27b9df72..3b9c01aab 100644 --- a/tests/test_fitting.py +++ b/tests/test_fitting.py @@ -7,6 +7,7 @@ import burnman from burnman.optimize.eos_fitting import fit_XPTp_data +from burnman.optimize.linear_fitting import weighted_constrained_least_squares from burnman.optimize.nonlinear_fitting import ( NonLinearModel, nonlinear_least_squares_fit, @@ -108,6 +109,91 @@ def normal(self, x, flag): nonlinear_least_squares_fit(model=fitted_curve, param_tolerance=1.0e-5) self.assertArraysAlmostEqual([fitted_curve.WSS], [10.486904577]) + def test_basic_unweighted_least_squares_exact_answer(self): + """Test simple unweighted least squares with exact answer.""" + A = np.array([[1, 0], [1, 1]]) + b = np.array([1, 3]) + popt, pcov, res = weighted_constrained_least_squares(A, b) + self.assertArraysAlmostEqual(popt, [1.0, 2.0]) + self.assertArraysAlmostEqual(pcov[0], [1.0, -1.0]) + self.assertArraysAlmostEqual(pcov[1], [-1.0, 2.0]) + self.assertAlmostEqual(res, 0.0) + + def test_basic_unweighted_least_squares(self): + """Test simple unweighted least squares.""" + A = np.array([[1, 1], [1, 2], [1, 3]]) + b = np.array([1, 2, 2]) + popt, pcov, res = weighted_constrained_least_squares(A, b) + self.assertArraysAlmostEqual(popt, [2.0 / 3.0, 1.0 / 2.0]) + self.assertArraysAlmostEqual(pcov[0], [7.0 / 3.0, -1.0]) + self.assertArraysAlmostEqual(pcov[1], [-1.0, 1.0 / 2.0]) + self.assertAlmostEqual(res, 1.0 / 6.0) + + def test_basic_unit_weighted_least_squares(self): + """Test simple unit weighted least squares.""" + A = np.array([[1, 1], [1, 2], [1, 3]]) + b = np.array([1, 2, 2]) + Cov_b = np.diag([1, 1, 1]) + popt, pcov, res = weighted_constrained_least_squares(A, b, Cov_b) + self.assertArraysAlmostEqual(popt, [2.0 / 3.0, 1.0 / 2.0]) + self.assertArraysAlmostEqual(pcov[0], [7.0 / 3.0, -1.0]) + self.assertArraysAlmostEqual(pcov[1], [-1.0, 1.0 / 2.0]) + self.assertAlmostEqual(res, 1.0 / 6.0) + + def test_unevenly_weighted_least_squares(self): + """Test weighted least squares with non-identity covariance.""" + A = np.array([[1, 1], [1, 2], [1, 3]]) + b = np.array([1, 2, 2]) + Cov_b = np.diag([1, 1, 2]) # higher weight on 3rd observation + popt, pcov, res = weighted_constrained_least_squares(A, b, Cov_b=Cov_b) + self.assertArraysAlmostEqual(popt, [4.0 / 7.0, 4.0 / 7.0], tol_zero=1.0e-12) + self.assertArraysAlmostEqual(pcov[0], [19.0 / 7.0, -9.0 / 7.0]) + self.assertArraysAlmostEqual(pcov[1], [-9.0 / 7.0, 5.0 / 7.0]) + self.assertAlmostEqual(res, 1.0 / 7.0) + + def test_equality_constraints_least_squares(self): + """Test equality constraint (force x1 + x2 = 1).""" + A = np.array([[1, 2], [3, 4], [5, 6]]) + b = np.array([7, 8, 9]) + C = np.array([[1, 1]]) # x1 + x2 = 1 + c = np.array([0.5]) + popt, _, res = weighted_constrained_least_squares( + A, b, equality_constraints=[C, c] + ) + self.assertArraysAlmostEqual(popt, [-6.0, 6.5], tol_zero=1.0e-12) + self.assertAlmostEqual(res, 0.0) + + def test_inequality_constraints_least_squares(self): + """Test inequality constraint (x1 >= 2, x2 >= 0).""" + A = np.array([[1, 1], [1, 2], [1, 3]]) + b = np.array([1, 2, 2]) + D = np.array([[-1, 0], [0, -1]]) + d = np.array([-2, 0]) + popt, pcov, res = weighted_constrained_least_squares( + A, b, inequality_constraints=[D, d] + ) + self.assertArraysAlmostEqual(popt, [2.0, 0.0], tol_zero=1.0e-12) + self.assertArraysAlmostEqual(pcov[0], [7.0 / 3.0, -1.0]) + self.assertArraysAlmostEqual(pcov[1], [-1.0, 1.0 / 2.0]) + self.assertAlmostEqual(res, 1.0) + + def test_rank_deficient_least_squares(self): + """Test rank-deficient A with allow_rank_deficient=False → expect error.""" + A = np.array([[1, 2], [2, 4], [3, 6]]) + b = np.array([1, 2, 3]) + with self.assertRaises(Exception): + _, _, _ = weighted_constrained_least_squares(A, b) + + def test_rank_deficient_allowed_least_squares(self): + """Test rank-deficient A with allow_rank_deficient=True.""" + A = np.array([[1, 2], [2, 4], [3, 6]]) # rank 1 + b = np.array([1, 2, 3]) + popt, pcov, res = weighted_constrained_least_squares( + A, b, allow_rank_deficient=True + ) + self.assertAlmostEqual(res, 0.0) + self.assertArraysAlmostEqual(A @ popt, b) + def test_fit_PTV_data(self): fo = burnman.minerals.HP_2011_ds62.fo() From 9546b993924b5ab82841482864c98ceed8ce5412 Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Fri, 31 Oct 2025 17:01:39 +0000 Subject: [PATCH 3/4] added example for optimal thermobarometry --- docs/examples.rst | 1 + examples/example_optimal_thermobarometry.py | 147 ++++++++++++++++++++ 2 files changed, 148 insertions(+) create mode 100644 examples/example_optimal_thermobarometry.py diff --git a/docs/examples.rst b/docs/examples.rst index 6f1ade968..c3ca71994 100644 --- a/docs/examples.rst +++ b/docs/examples.rst @@ -167,6 +167,7 @@ Advanced examples: - :mod:`~examples.example_fit_data`, - :mod:`~examples.example_fit_eos`, - :mod:`~examples.example_fit_solution`, + - :mod:`~examples.example_optimal_thermobarometry`, - :mod:`~examples.example_equilibrate`, and - :mod:`~examples.example_olivine_binary`. diff --git a/examples/example_optimal_thermobarometry.py b/examples/example_optimal_thermobarometry.py new file mode 100644 index 000000000..c3440e8c3 --- /dev/null +++ b/examples/example_optimal_thermobarometry.py @@ -0,0 +1,147 @@ +# 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. + + +""" +example_optimal_thermobarometry +------------------------------- + +This example script is intended to demonstrate optimal thermobarometry +using BurnMan. The technique is based on the work of Powell and Holland +(1994; American Mineralogist 79 (1-2): 120-133). +We cover importing BurnMan modules, creating a composite +material representing a mineral assemblage, fitting the compositions of +the constituent minerals to their solution models, and estimating the +pressure and temperature conditions of equilibration based on the +mineral compositions and their uncertainties. Finally, we print +the estimated conditions along with their uncertainties and correlations. + +*Uses:* + +* :doc:`mineral_database` +* :func:`burnman.optimize.composition_fitting.fit_composition_to_solution` +* :func:`burnman.optimize.thermobarometry.estimate_conditions` + + +*Demonstrates:* + +* creating mineral assemblages +* fitting mineral compositions to solution models +* estimating pressure and temperature conditions of equilibration +""" + +import numpy as np +from burnman import Composite +from burnman.minerals import mp50MnNCKFMASHTO, HP_2011_ds62 +from burnman.optimize.composition_fitting import fit_composition_to_solution +from burnman.tools.thermobarometry import estimate_conditions + +if __name__ == "__main__": + + # 1) Define observed mineral assemblage + mu = mp50MnNCKFMASHTO.mu() + bi = mp50MnNCKFMASHTO.bi() + g = mp50MnNCKFMASHTO.g() + ilmm = mp50MnNCKFMASHTO.ilmm() + st = mp50MnNCKFMASHTO.st() + q = HP_2011_ds62.q() + + assemblage = Composite([mu, bi, g, ilmm, st, q]) + + # 2) Fit the measured compositions (in molar amounts) to the solution + # models for each mineral in the assemblage. + # The amounts do not need to sum to any particular number, + # as the fitting function will normalize them appropriately. + # The compositions correspond approximately to the expected + # compositions at about 4 kbar and 873 K. + + # a) Muscovite + fitted_species = ["Na", "Ca", "K", "Fe", "Mg", "Al", "Si"] + species_amounts = np.array([0.40, 0.01, 0.55, 0.01, 0.01, 3.00, 3.00]) + species_covariances = np.diag( + np.array([0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01]) ** 2 + ) + popt, pcov, res = fit_composition_to_solution( + mu, fitted_species, species_amounts, species_covariances + ) + mu.set_composition(popt) + mu.compositional_covariances = pcov + + # b) Biotite (requires parameter constraining order state of + # Mg and Fe on the first two sites) + fitted_species = ["Mn", "Fe", "Mg", "Al", "Si", "Ti", "Mg_M3"] + species_amounts = np.array([0.01, 1.50, 1.00, 1.65, 2.65, 0.20, 0.10]) + species_covariances = np.diag( + np.array([0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.1]) ** 2 + ) + variable_conversions = {"Mg_M3": {"Mgmthree_A": 1.0}} + popt, pcov, res = fit_composition_to_solution( + bi, fitted_species, species_amounts, species_covariances, variable_conversions + ) + bi.set_composition(popt) + bi.compositional_covariances = pcov + + # c) Garnet + fitted_species = ["Mn", "Fe", "Mg", "Ca", "Al", "Si"] + species_amounts = np.array([0.25, 2.30, 0.40, 0.05, 2.00, 3.00]) + species_covariances = np.diag(np.array([0.01, 0.01, 0.01, 0.01, 0.01, 0.01]) ** 2) + popt, pcov, res = fit_composition_to_solution( + g, fitted_species, species_amounts, species_covariances + ) + g.set_composition(popt) + g.compositional_covariances = pcov + + # d) Ilmenite (requires parameter constraining order state of Fe and Ti) + fitted_species = ["Mn", "Fe", "Ti", "Mg", "Fe2+_A"] + species_amounts = np.array([0.05, 1.0, 0.90, 0.05, 0.4]) + species_covariances = np.diag(np.array([0.01, 0.01, 0.01, 0.01, 0.2]) ** 2) + species_conversions = {"Fe2+_A": {"Fea_A": 1.0}} + popt, pcov, res = fit_composition_to_solution( + ilmm, fitted_species, species_amounts, species_covariances, species_conversions + ) + ilmm.set_composition(popt) + ilmm.compositional_covariances = pcov + + # e) Staurolite + fitted_species = ["Mn", "Fe", "Mg", "Al", "Si", "Ti"] + species_amounts = np.array([0.05, 3.30, 0.72, 17.78, 7.50, 0.12]) + species_covariances = np.diag(np.array([0.01, 0.01, 0.01, 0.01, 0.01, 0.01]) ** 2) + popt, pcov, res = fit_composition_to_solution( + st, fitted_species, species_amounts, species_covariances + ) + st.set_composition(popt) + st.compositional_covariances = pcov + + # f) Quartz (no fitting needed) + + # 3) Estimate the pressure and temperature conditions of equilibration + # based on the fitted compositions and their uncertainties, and + # also the endmember covariances provided by the underlying dataset. + res = estimate_conditions( + assemblage, + dataset_covariances=HP_2011_ds62.cov(), + guessed_conditions=np.array([0.5e9, 500.0]), + small_fraction_tol=0.01, + ) + + # 4) Print the estimated conditions along with their uncertainties + # and correlations. + print( + f"Estimated Pressure: {assemblage.pressure/1.e9:.2f} " + f"+/- {np.sqrt(res.xcov[0, 0])/1.e9:.2f} GPa" + ) + print( + f"Estimated Temperature: {assemblage.temperature:.2f} " + f"+/- {np.sqrt(res.xcov[1, 1]):.2f} K" + ) + print(f"Correlation between P and T: {res.xcorr:.4f}") + print(f"Number of Reactions: {res.n_reactions}") + print(f"Number of Parameters: {res.n_params}") + print(f"Degrees of Freedom: {res.degrees_of_freedom}") + print(f"Reduced Chi-squared: {res.reduced_chisqr:.4f}") + print(f"Fit (sqrt reduced chi-squared): {res.fit:.4f}") + print("Weighted reaction affinities:") + np.set_printoptions(precision=2) + print(res.weighted_affinities) From e06fa8b3fe297442b42028f562e15a611f3f3cae Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Fri, 31 Oct 2025 17:02:51 +0000 Subject: [PATCH 4/4] added test output --- examples/example_optimal_thermobarometry.py | 18 ++++++++++-------- .../ref/example_optimal_thermobarometry.py.out | 8 ++++++++ 2 files changed, 18 insertions(+), 8 deletions(-) create mode 100644 misc/ref/example_optimal_thermobarometry.py.out diff --git a/examples/example_optimal_thermobarometry.py b/examples/example_optimal_thermobarometry.py index c3440e8c3..449d6e553 100644 --- a/examples/example_optimal_thermobarometry.py +++ b/examples/example_optimal_thermobarometry.py @@ -133,15 +133,17 @@ f"+/- {np.sqrt(res.xcov[0, 0])/1.e9:.2f} GPa" ) print( - f"Estimated Temperature: {assemblage.temperature:.2f} " - f"+/- {np.sqrt(res.xcov[1, 1]):.2f} K" + f"Estimated Temperature: {assemblage.temperature:.0f} " + f"+/- {np.sqrt(res.xcov[1, 1]):.0f} K" ) - print(f"Correlation between P and T: {res.xcorr:.4f}") + print(f"Correlation between P and T: {res.xcorr:.2f}") print(f"Number of Reactions: {res.n_reactions}") print(f"Number of Parameters: {res.n_params}") print(f"Degrees of Freedom: {res.degrees_of_freedom}") - print(f"Reduced Chi-squared: {res.reduced_chisqr:.4f}") - print(f"Fit (sqrt reduced chi-squared): {res.fit:.4f}") - print("Weighted reaction affinities:") - np.set_printoptions(precision=2) - print(res.weighted_affinities) + print(f"Reduced Chi-squared: {res.reduced_chisqr:.2f}") + print(f"Fit (sqrt reduced chi-squared): {res.fit:.2f}") + + # Uncomment to see the weighted reaction affinities + # print("Weighted reaction affinities:") + # np.set_printoptions(precision=2) + # print(res.weighted_affinities) diff --git a/misc/ref/example_optimal_thermobarometry.py.out b/misc/ref/example_optimal_thermobarometry.py.out new file mode 100644 index 000000000..e94427c34 --- /dev/null +++ b/misc/ref/example_optimal_thermobarometry.py.out @@ -0,0 +1,8 @@ +Estimated Pressure: 0.43 +/- 0.04 GPa +Estimated Temperature: 887 +/- 11 K +Correlation between P and T: 0.11 +Number of Reactions: 11 +Number of Parameters: 2 +Degrees of Freedom: 9 +Reduced Chi-squared: 0.20 +Fit (sqrt reduced chi-squared): 0.45