Merge pull request #695 from bobmyhill/avPT_w_free_c_vectors

Allow free compositional vectors in optimal thermobarometry problems

Author: bobmyhill

burnman/classes/solution.py

@@ -131,41 +131,62 @@ def formula(self):
     @material_property
     def site_occupancies(self):
         """
-        :returns: The fractional occupancies of species on each site.
-        :rtype: list of OrderedDicts
+        :returns: The fractional occupancies of species on each site
+            as a flat vector.
+        :rtype: numpy.ndarray
         """
         f = self.molar_fractions
-        occs = np.einsum("ij, i", self.solution_model.endmember_occupancies, f)
-        site_occs = []
-        k = 0
-        for site in self.solution_model.sites:
-            site_occs.append(OrderedDict())
-            for species in site:
-                site_occs[-1][species] = occs[k]
-                k += 1
-
-        return site_occs
+        occs = f.dot(self.solution_model.endmember_occupancies)
+        return occs
 
-    def site_formula(self, precision=2):
+    def site_formula(
+        self,
+        precision=2,
+        print_absent_species=True,
+        site_occupancies=None,
+    ):
         """
         Returns the molar chemical formula of the solution with
             site occupancies. For example, [Mg0.4Fe0.6]2SiO4.
 
         :param precision: Precision with which to print the site occupancies
         :type precision: int
 
+        :param print_absent_species: If True, prints site occupancies that are absent (i.e., zero).
+        :type print_absent_species: bool
+
+        :param site_occupancies: Optional site occupancies to use instead of
+            the current ones.
+        :type site_occupancies: numpy.ndarray
+
         :returns: Molar chemical formula of the solution with site occupancies
         :rtype: str
         """
-        split_empty = self.solution_model.empty_formula.split("[")
-        formula = split_empty[0]
-        for i, site_occs in enumerate(self.site_occupancies):
+        if site_occupancies is not None:
+            occs = site_occupancies
+        else:
+            occs = self.site_occupancies
+
+        i = 0
+        # Split the empty formula into site sections.
+        # This allows us to reconstruct the formula including
+        # occupancies and anything outside the sites.
+        split_formula = self.solution_model.empty_formula.split("[")
+        formula = split_formula[0]
+        for j, site in enumerate(self.solution_model.sites):
             formula += "["
-            for species, occ in site_occs.items():
-                if np.abs(occ) < 1.0e-12:
-                    occ = np.abs(occ)
+            for species in site:
+                # avoid printing -0.00
+                occ = occs[i] if np.abs(occs[i]) >= 1.0e-12 else 0.0
+
+                # skip absent species if requested
+                if np.abs(occ) < 1.0e-12 and not print_absent_species:
+                    i += 1
+                    continue
                 formula += f"{species}{occ:0.{precision}f}"
-            formula += split_empty[i + 1]
+                i += 1
+            formula += split_formula[j + 1]
+
         return formula
 
     @material_property

burnman/tools/thermobarometry.py

@@ -107,6 +107,12 @@ def estimate_conditions(
     for a given mineral assemblage. Algorithm based on Powell and Holland (1994).
 
     :param assemblage: The mineral assemblage for which to perform the inversion.
+        Each solution phase in the assemblage must have its composition set along with
+        its compositional covariance matrix (called `compositional_covariances`).
+        If there are compositional degrees of freedom, they can be added by setting
+        the attribute `free_compositional_vectors` on the relevant solution phases, where each
+        row of the array corresponds to a free compositional vector, and the columns correspond
+        to the amounts of endmembers of that phase in each vector.
     :type assemblage: Assemblage
 
     :param dataset_covariances: The covariance data from the thermodynamic dataset.
@@ -146,13 +152,21 @@ def estimate_conditions(
         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
 
+    # Count the number of free compositional vectors across all phases
+    n_free_vectors = 0
+    for phase in assemblage.phases:
+        if hasattr(phase, "free_compositional_vectors"):
+            n_free_vectors += phase.free_compositional_vectors.shape[0]
+            phase.baseline_composition = phase.molar_fractions.copy()
+
+    # Check if P and/or T are fixed
     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!")
+    if P_fixed and T_fixed and n_free_vectors == 0:
+        raise Exception(
+            "Both pressure and temperature cannot be fixed if there are no free compositional vectors!"
+        )
 
     # Build the reaction matrix R, possibly reducing the number of endmembers
     # by excluding components with small molar fractions
@@ -189,7 +203,18 @@ def estimate_conditions(
         raise Exception("No independent reactions found for the given assemblage!")
 
     def calculate_cov_a(assemblage, theta, dataset_covariances):
-        assemblage.set_state(*theta)
+        if n_free_vectors > 0:
+            vi = 2
+            for phase in assemblage.phases:
+                if hasattr(phase, "free_compositional_vectors"):
+                    n_free_vectors_in_phase = phase.free_compositional_vectors.shape[0]
+                    v = np.array(theta[vi : vi + n_free_vectors_in_phase])
+                    dc = phase.free_compositional_vectors.T.dot(v)
+                    new_composition = phase.baseline_composition + dc
+                    phase.set_composition(new_composition)
+                    vi += n_free_vectors_in_phase
+
+        assemblage.set_state(*theta[:2])
         cov1 = _build_endmember_covariance_matrix(assemblage, dataset_covariances)
         cov2 = _build_solution_covariance_matrix(assemblage)
         cov_mu = cov1 + cov2
@@ -207,18 +232,36 @@ def chisqr(args):
         """
         # These are the chemical potential covariances for all endmembers
         cov_a = calculate_cov_a(
-            assemblage, [args[0] * P_scaling, args[1]], dataset_covariances
+            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)
+        chisqr = a.T.dot(np.linalg.pinv(cov_a)).dot(a)
         return chisqr
 
     # Minimize the chi-squared function
-    x0 = [guessed_conditions[0] / P_scaling, guessed_conditions[1]]
+    x0 = [
+        guessed_conditions[0] / P_scaling,
+        guessed_conditions[1],
+        *[0.0] * n_free_vectors,
+    ]
     bounds = [
         (pressure_bounds[0] / P_scaling, pressure_bounds[1] / P_scaling),
         (temperature_bounds[0], temperature_bounds[1]),
+        *[(None, None)] * n_free_vectors,
     ]
+
+    if n_free_vectors > 0:
+        # Nelder-Mead is more robust for problems where the feasible region
+        # may be complex due to compositional degrees of freedom
+        res = minimize(
+            chisqr,
+            x0,
+            method="Nelder-Mead",
+            bounds=bounds,
+            options={"maxiter": max_it},
+        )
+        x0 = res.x
+
     res = minimize(
         chisqr,
         x0,

burnman/utils/chemistry.py

@@ -561,11 +561,20 @@ def ordered_compositional_array(formulae, elements):
     return formula_array
 
 
-def formula_to_string(formula):
+def formula_to_string(formula, use_fractions=True, significant_figures=2):
     """
     :param formula: Chemical formula
     :type formula: dict or Counter
 
+    :param use_fractions: Whether or not to use fractions in the output string.
+        If False, numbers are rounded to the number of significant figures
+        given by the significant_figures parameter.
+    :type use_fractions: bool
+
+    :param significant_figures: Number of significant figures to use if
+        use_fractions=False.
+    :type significant_figures: int
+
     :returns: A formula string, with element order as given in the list
         IUPAC_element_order.
         If one or more keys in the dictionary are not one of the elements
@@ -574,21 +583,40 @@ def formula_to_string(formula):
     """
 
     formula_string = ""
-    for e in IUPAC_element_order:
-        if e in formula and np.abs(formula[e]) > 1.0e-12:
-            if np.abs(formula[e] - 1.0) < 1.0e-12:
-                formula_string += e
-            else:
-                formula_string += e + str(nsimplify(formula[e]))
-
-    for e in formula:
-        if e not in IUPAC_element_order:
+    if use_fractions:
+        for e in IUPAC_element_order:
             if e in formula and np.abs(formula[e]) > 1.0e-12:
                 if np.abs(formula[e] - 1.0) < 1.0e-12:
                     formula_string += e
                 else:
                     formula_string += e + str(nsimplify(formula[e]))
 
+        for e in formula:
+            if e not in IUPAC_element_order:
+                if e in formula and np.abs(formula[e]) > 1.0e-12:
+                    if np.abs(formula[e] - 1.0) < 1.0e-12:
+                        formula_string += e
+                    else:
+                        formula_string += e + str(nsimplify(formula[e]))
+
+    else:
+        for e in IUPAC_element_order:
+            if e in formula and np.abs(formula[e]) > 1.0e-12:
+                if np.abs(formula[e] - 1.0) < 1.0e-12:
+                    formula_string += e
+                else:
+                    formula_string += e + str(round(formula[e], significant_figures))
+
+        for e in formula:
+            if e not in IUPAC_element_order:
+                if e in formula and np.abs(formula[e]) > 1.0e-12:
+                    if np.abs(formula[e] - 1.0) < 1.0e-12:
+                        formula_string += e
+                    else:
+                        formula_string += e + str(
+                            round(formula[e], significant_figures)
+                        )
+
     return formula_string
 
 

examples/example_optimal_thermobarometry.py

@@ -37,6 +37,7 @@
 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
+from burnman.utils.chemistry import formula_to_string
 
 if __name__ == "__main__":
 
@@ -70,7 +71,9 @@
     mu.compositional_covariances = pcov
 
     # b) Biotite (requires parameter constraining order state of
-    # Mg and Fe on the first two sites)
+    # Mg and Fe on the first two sites).
+    # Here, we will allow the Mg-Fe ordering to vary to improve the fit
+    # by adding a free compositional vector.
     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(
@@ -82,6 +85,7 @@
     )
     bi.set_composition(popt)
     bi.compositional_covariances = pcov
+    bi.free_compositional_vectors = Composite([bi]).reaction_basis
 
     # c) Garnet
     fitted_species = ["Mn", "Fe", "Mg", "Ca", "Al", "Si"]
@@ -94,6 +98,8 @@
     g.compositional_covariances = pcov
 
     # d) Ilmenite (requires parameter constraining order state of Fe and Ti)
+    # Here, we will allow the Fe-Ti ordering to vary to improve the fit
+    # by adding a free compositional vector.
     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)
@@ -103,6 +109,7 @@
     )
     ilmm.set_composition(popt)
     ilmm.compositional_covariances = pcov
+    ilmm.free_compositional_vectors = Composite([ilmm]).reaction_basis
 
     # e) Staurolite
     fitted_species = ["Mn", "Fe", "Mg", "Al", "Si", "Ti"]
@@ -128,20 +135,74 @@
 
     # 4) Print the estimated conditions along with their uncertainties
     #    and correlations.
+    print("Results of thermobarometric inversion:")
     print(
-        f"Estimated Pressure: {assemblage.pressure/1.e9:.2f} "
+        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:.0f} "
+        f"    Estimated Temperature: {assemblage.temperature:.0f} "
         f"+/- {np.sqrt(res.xcov[1, 1]):.0f} K"
     )
-    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:.2f}")
-    print(f"Fit (sqrt reduced chi-squared): {res.fit:.2f}")
+    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:.2f}")
+    print(f"    Fit (sqrt reduced chi-squared): {res.fit:.2f}")
+
+    # The names of the sites in these models are a bit cumbersome,
+    # so we will do some string replacements to make them more
+    # readable in the output.
+    replace_strings = [
+        ["Fethree", "Fe3+"],
+        ["monetwo", ""],
+        ["mtwoa", ""],
+        ["mtwob", ""],
+        ["mthree", ""],
+        ["x", ""],
+        ["y", ""],
+        ["tone", ""],
+        ["t", ""],
+        ["Ka", "K"],
+        ["Naa", "Na"],
+        ["Caa", "Ca"],
+        ["Fea", "Fe"],
+        ["Mga", "Mg"],
+        ["Mna", "Mn"],
+        ["Tia", "Ti"],
+        ["Fe3+a", "Fe3+"],
+        ["Oh", "(OH)"],
+        ["b", ""],
+        ["v", ""],
+    ]
+
+    print("\nMineral compositions and site occupancies:")
+    for phase in assemblage.phases:
+        np.set_printoptions(precision=2, formatter={"all": lambda x: f"{x:0.2f}"})
+        print(f"{phase.name}:")
+
+        if hasattr(phase, "free_compositional_vectors"):
+            free_site_occupancy_vectors = phase.free_compositional_vectors.dot(
+                phase.solution_model.endmember_occupancies
+            )
+            for v in free_site_occupancy_vectors:
+                site_formula = phase.site_formula(
+                    precision=2, site_occupancies=v, print_absent_species=False
+                )
+                for old, new in replace_strings:
+                    site_formula = site_formula.replace(old, new)
+                print(f"    free site occupancy vector: {site_formula}")
+
+        print(
+            f"    composition: {formula_to_string(phase.formula, use_fractions=False)}"
+        )
+        if hasattr(phase, "molar_fractions"):
+            site_formula = phase.site_formula(2)
+            for old, new in replace_strings:
+                site_formula = site_formula.replace(old, new)
+            print(f"    site occupancies: {site_formula}")
+    print("")
 
     # Uncomment to see the weighted reaction affinities
     # print("Weighted reaction affinities:")