Merge pull request #554 from bobmyhill/fitsolution2

make fit solution more accessible

Author: bobmyhill

Diff for: burnman/classes/solution.py

@@ -8,6 +8,8 @@
 
 import numpy as np
 from sympy import Matrix, nsimplify
+from collections import OrderedDict
+
 from .material import material_property, cached_property
 from .mineral import Mineral
 from .solutionmodel import MechanicalSolution
@@ -116,9 +118,49 @@ def set_state(self, pressure, temperature):
     def formula(self):
         """
         Returns molar chemical formula of the solution.
+        :rtype: Counter
         """
         return sum_formulae(self.endmember_formulae, self.molar_fractions)
 
+    @material_property
+    def site_occupancies(self):
+        """
+        :returns: The fractional occupancies of species on each site.
+        :rtype: list of OrderedDicts
+        """
+        occs = np.einsum(
+            "ij, i", self.solution_model.endmember_occupancies, self.molar_fractions
+        )
+        site_occs = []
+        k = 0
+        for i in range(self.solution_model.n_sites):
+            site_occs.append(OrderedDict())
+            for j in range(len(self.solution_model.sites[i])):
+                site_occs[-1][self.solution_model.sites[i][j]] = occs[k]
+                k += 1
+
+        return site_occs
+
+    def site_formula(self, precision=2):
+        """
+        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
+
+        :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):
+            formula += "["
+            for species, occ in site_occs.items():
+                formula += f"{species}{occ:0.{precision}f}"
+            formula += split_empty[i + 1]
+        return formula
+
     @material_property
     def activities(self):
         """

Diff for: burnman/optimize/composition_fitting.py

@@ -6,7 +6,40 @@
 from __future__ import absolute_import
 
 import numpy as np
+from copy import deepcopy
+from collections import namedtuple
 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
+
+
+class DummyCompositionSolution(Solution):
+    """
+    This is a dummy base class for a solution object to
+    facilitate composition fitting when no solution model has
+    been prepared. 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,), {"formulas": endmember_site_formulae}
+        )()
+        self.solution_model.endmembers = [None for f in endmember_site_formulae]
+        process_solution_chemistry(self.solution_model)
 
 
 def fit_composition_to_solution(
@@ -59,17 +92,15 @@ def fit_composition_to_solution(
     :rtype: tuple of 1D numpy.array, 2D numpy.array and float
     """
 
-    n_vars = len(fitted_variables)
-    n_mbrs = len(solution.endmembers)
-
-    solution_variables = solution.elements
+    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_sol_vars = solution_matrix.shape[1]
+    n_vars = len(fitted_variables)
+    n_mbrs, n_sol_vars = solution_matrix.shape
 
     if variable_conversions is not None:
         solution_matrix = np.hstack(

Diff for: burnman/utils/chemistry.py

@@ -308,6 +308,12 @@ def process_solution_chemistry(solution_model):
 
         * solution_formulae [list of dictionaries]
             List of endmember formulae in dictionary form.
+        * empty_formula [string]
+            Abbreviated chemical formula with sites denoted by empty
+            square brackets.
+        * general_formula [string]
+            General chemical formula with sites denoted by
+            square brackets filled with a comma-separated list of species
         * n_sites [integer]
             Number of sites in the solution.
             Should be the same for all endmembers.
@@ -442,6 +448,14 @@ def process_solution_chemistry(solution_model):
         "ij, ij->ij", endmember_occupancies, site_multiplicities
     )
 
+    solution_model.empty_formula = re.sub(
+        "([\[]).*?([\]])", "\g<1>\g<2>", solution_model.formulas[0]
+    )
+    split_empty = solution_model.empty_formula.split("[")
+    solution_model.general_formula = split_empty[0]
+    for i in range(n_sites):
+        solution_model.general_formula += f"[{','.join(sites[i])}{split_empty[i+1]}"
+
 
 def site_occupancies_to_strings(
     site_species_names, site_multiplicities, endmember_occupancies

Diff for: examples/example_fit_composition.py

@@ -34,7 +34,9 @@
 import itertools
 
 from burnman import minerals
+from burnman.utils.chemistry import formula_to_string
 from burnman.optimize.composition_fitting import fit_composition_to_solution
+from burnman.optimize.composition_fitting import DummyCompositionSolution
 from burnman.optimize.composition_fitting import (
     fit_phase_proportions_to_bulk_composition,
 )
@@ -96,7 +98,6 @@
     popt, pcov, res = fit_composition_to_solution(
         gt, fitted_species, species_amounts, species_covariances, species_conversions
     )
-
     # We can set the composition of gt using the optimized parameters
     gt.set_composition(popt)
 
@@ -108,9 +109,36 @@
             f"{gt.molar_fractions[i]:.3f} +/- "
             f"{np.sqrt(pcov[i][i]):.3f}"
         )
-    print()
-    print(f"Weighted residual: {res:.3f}")
 
+    # If you don't have a full solution model prepared,
+    # you can instead use the helper class
+    # DummyCompositionSolution
+    element_formulae = [
+        "Mg3Al2Si3O12",
+        "Fe3Al2Si3O12",
+        "Ca3Al2Si3O12",
+        "Ca3Fe2Si3O12",
+        "Mg3Cr2Si3O12",
+    ]
+    site_formulae = [
+        "[Mg]3[Al]2Si3O12",
+        "[Fe]3[Al]2Si3O12",
+        "[Ca]3[Al]2Si3O12",
+        "[Ca]3[Fef]2Si3O12",
+        "[Mg]3[Cr]2Si3O12",
+    ]
+    gt = DummyCompositionSolution(element_formulae, site_formulae)
+    popt1, pcov1, res1 = fit_composition_to_solution(
+        gt, fitted_species, species_amounts, species_covariances, species_conversions
+    )
+
+    assert res == res1
+    gt.set_composition(popt)
+
+    print("\nSite formula:")
+    print(gt.site_formula(2))
+
+    print(f"\nWeighted residual: {res:.3f}")
     """
     Example 2
     ---------

Diff for: misc/ref/example_fit_composition.py.out

@@ -23,4 +23,7 @@ gr: 0.619 +/- 0.004
 andr: 0.048 +/- 0.004
 knor: 0.000 +/- 0.004
 
+Site formula:
+[Mg0.00Fe0.33Ca0.67]3[Al0.95Fef0.05Cr0.00]2Si3O12
+
 Weighted residual: 0.519