Merge pull request #693 from bobmyhill/optimal_pt

added optimal PT calculations

Author: bobmyhill

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):
         """

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.
 
 """

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)

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

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