Merge pull request #32 from ycexiao/multiple-contribution

feat: support multiple-contribution refinement

Author: sbillinge

news/multiple-contribution.rst

@@ -0,0 +1,23 @@
+**Added:**
+
+* No news added: Support multiple contribution refinement.
+
+**Changed:**
+
+* <news item>
+
+**Deprecated:**
+
+* <news item>
+
+**Removed:**
+
+* <news item>
+
+**Fixed:**
+
+* <news item>
+
+**Security:**
+
+* <news item>

src/diffpy/apps/refinebase/parametric_model.py

@@ -1,13 +1,17 @@
 import logging
 import re
+from pathlib import Path
 
 import networkx as nx
+from pyobjcryst import loadCrystal
 
 from diffpy.srfit.fitbase import FitContribution
 from diffpy.srfit.fitbase.parameter import Parameter, ParameterProxy
 from diffpy.srfit.pdf.pdfgenerator import PDFGenerator
 from diffpy.srfit.structure import constrain_as_space_group
-from diffpy.structure import Structure
+from diffpy.srfit.structure.diffpyparset import DiffpyStructureParSet
+from diffpy.srfit.structure.objcrystparset import ObjCrystCrystalParSet
+from diffpy.structure.parsers import get_parser
 
 # NOTE: MCP server prefers logging for output
 logger = logging.getLogger(__name__)
@@ -152,9 +156,13 @@ def residual(self):
 
 
 class ParametricModelEquation(ParametricModel):
-    def __init__(self, name, equation_str=None):
+    def __init__(self, name, equation_str=None, from_model_name=None):
         super().__init__(name=name)
         self.equation_str = None
+        if from_model_name is not None:
+            for name, obj in from_model_name.calc_obj.__dict__.items():
+                if name not in ["name", "profile", "_observers"]:
+                    setattr(self.calc_obj, name, obj)
         if equation_str:
             self.set_equation(equation_str)
 
@@ -170,14 +178,6 @@ def set_equation(self, equation_str):
     def get_equation(self):
         return self.equation_str
 
-    def set_residual_equation(self, residual_equation_str):
-        self.residual_equation = residual_equation_str
-        self.calc_obj.set_residual_equation(residual_equation_str)
-        self._rebuild_graph()
-
-    def get_residual_equation(self):
-        return self.residual_equation
-
     def evaluate(self):
         yc = self.calc_obj._eq()
         if (
@@ -194,60 +194,171 @@ def residual(self):
 class ParametricModelPDF(ParametricModel):
     # NOTE: qmin, qmax, stype(scattering type) are meta handled
     #   throughout the loaded profile in the refinement session
-    def __init__(self, name, structure: Structure):
+    def __init__(
+        self,
+        name,
+        structure_file_path=None,
+        from_model_name=None,
+    ):
         super().__init__(name=name)
         self.calc_obj = PDFGenerator(name)
-        self.calc_obj.setStructure(structure)
+        # NOTE: Certain space groups require dual origin handling.
+        DUAL_ORIGIN_SG_NUMBERS = {
+            48,
+            50,
+            59,
+            68,
+            70,
+            85,
+            86,
+            88,
+            125,
+            126,
+            129,
+            130,
+            133,
+            134,
+            137,
+            138,
+            141,
+            142,
+            201,
+            203,
+            222,
+            224,
+            227,
+            228,
+        }
+        if structure_file_path is not None:
+            stru_parser = get_parser("auto")
+            structure = stru_parser.parse(
+                Path(structure_file_path).read_text()
+            )
+            sg = getattr(stru_parser, "spacegroup", None)
+            self.space_group_symbol = sg.short_name if sg is not None else "P1"
+            if sg.number in DUAL_ORIGIN_SG_NUMBERS:
+                structure = loadCrystal(structure_file_path)
+                self.calc_obj.setStructure(structure)
+            else:
+                self.calc_obj.setStructure(structure)
+        elif from_model_name is not None:
+            self.calc_obj.setPhase(from_model_name.calc_obj.phase)
+            self.space_group_symbol = from_model_name.space_group_symbol
+        else:
+            raise ValueError(
+                "Either structure_file or from_model must be provided."
+            )
+        self.sgpar_names = []
         self._rebuild_graph()
-        self._hide_dependent_parameters()
 
-    def _hide_dependent_parameters(self):
-        dependent_par_names = [
-            r"\.U21$",
-            r"\.U31$",
-            r"\.U32$",  # U21=U12, U31=U13, U32=U23
-            r"\.Biso",
-            r"\.B\d{2}",  # Bij = Uij * 8 * pi^2
-            r"\.occupancy$",  # occupancy=oc
-        ]
+    def _hide_dependent_parameters(self, use_uiso=True):
+        if use_uiso:
+            dependent_par_names = [
+                r"\.U21$",
+                r"\.U31$",
+                r"\.U32$",  # U21=U12, U31=U13, U32=U23
+                r"\.Biso",
+                r"\.B\d{2}",  # Bij = Uij * 8 * pi^2
+                r"\.occupancy$",  # occupancy=oc
+            ]
+        else:
+            dependent_par_names = [
+                r"\.B21$",
+                r"\.B31$",
+                r"\.B32$",
+                r"\.Uiso",
+                r"\.U\d{2}",
+                r"\.occupancy$",  # occupancy=oc
+            ]
         regex = re.compile("|".join(dependent_par_names))
         for par_name in self.parameters.keys():
             if regex.search(par_name):
                 self._graph.nodes[par_name]["constrained_or_constant"] = True
 
-    def constrain_symmetry(self, spacegroup_symbol):
-        space_group_parset = constrain_as_space_group(
-            self.calc_obj.phase, spacegroup_symbol
-        )
+    def constrain_symmetry(self, spacegroup_symbol=None, use_uiso=True):
+        if spacegroup_symbol is None:
+            spacegroup_symbol = self.space_group_symbol
+        if isinstance(self.calc_obj.phase, DiffpyStructureParSet):
+            space_group_parset = constrain_as_space_group(
+                self.calc_obj.phase, spacegroup_symbol
+            )
+            self._hide_dependent_parameters(use_uiso=use_uiso)
+        elif isinstance(self.calc_obj.phase, ObjCrystCrystalParSet):
+            if use_uiso is not False:
+                use_uiso = False
+                logger.warning(
+                    "ObjCrystCrystalParSet prefers using the letter B instead "
+                    "U for ADP parameters."
+                )
+            self._hide_dependent_parameters(use_uiso=use_uiso)
+        else:
+            raise ValueError(
+                "Unsupported calculation object type."
+                "Currently supported types are "
+                "DiffpyStructureParSet and ObjCrystCrystalParSet."
+            )
         # hide constrained parameters in the graph
-        symmetry_par_names = [
-            r"\.a$",
-            r"\.b$",
-            r"\.c$",
-            r"\.alpha$",
-            r"\.beta$",
-            r"\.gamma$",
-            r"\.x$",
-            r"\.y$",
-            r"\.z$",
-            r"\.Uiso$",
-            r"\.U11$",
-            r"\.U22$",
-            r"\.U33$",
-            r"\.U12$",
-            r"\.U13$",
-            r"\.U23$",
-        ]
+        if use_uiso:
+            symmetry_par_names = [
+                r"\.a$",
+                r"\.b$",
+                r"\.c$",
+                r"\.alpha$",
+                r"\.beta$",
+                r"\.gamma$",
+                r"\.x$",
+                r"\.y$",
+                r"\.z$",
+                r"\.Uiso$",
+                r"\.U11$",
+                r"\.U22$",
+                r"\.U33$",
+                r"\.U12$",
+                r"\.U13$",
+                r"\.U23$",
+            ]
+        else:
+            symmetry_par_names = [
+                r"\.a$",
+                r"\.b$",
+                r"\.c$",
+                r"\.alpha$",
+                r"\.beta$",
+                r"\.gamma$",
+                r"\.x$",
+                r"\.y$",
+                r"\.z$",
+                r"\.Biso$",
+                r"\.B11$",
+                r"\.B22$",
+                r"\.B33$",
+                r"\.B12$",
+                r"\.B13$",
+                r"\.B23$",
+            ]
         free_variables = []
-        for latpar in space_group_parset.latpars:
-            free_variables.append(latpar)
-        for adpar in space_group_parset.adppars:
-            free_variables.append(adpar)
-        for xyzpar in space_group_parset.xyzpars:
-            free_variables.append(xyzpar)
+        if isinstance(self.calc_obj.phase, DiffpyStructureParSet):
+            for latpar in space_group_parset.latpars:
+                free_variables.append(latpar)
+            for adpar in space_group_parset.adppars:
+                free_variables.append(adpar)
+            for xyzpar in space_group_parset.xyzpars:
+                free_variables.append(xyzpar)
+        elif isinstance(self.calc_obj.phase, ObjCrystCrystalParSet):
+            for par in self.calc_obj.phase.sgpars:
+                free_variables.append(par)
         for i in range(len(free_variables)):
             while isinstance(free_variables[i], ParameterProxy):
                 free_variables[i] = free_variables[i].par
+            organized_name = ".".join(
+                [
+                    obj.name
+                    for obj in self.calc_obj._locate_managed_object(
+                        free_variables[i]
+                    )
+                ]
+            )
+            self.sgpar_names.append(organized_name)
         symmetry_par_regex = re.compile("|".join(symmetry_par_names))
         for par_name, par in self.parameters.items():
             if symmetry_par_regex.search(par_name):