Support for SAXS fitting (v2)

build_tools/requirements.txt

@@ -16,7 +16,7 @@ numpy
 packaging
 periodictable
 platformdirs
-pyausaxs==1.0.4
+pyausaxs
 pybind11
 pylint
 pyopencl

src/sas/qtgui/Calculators/GenericScatteringCalculator.py

@@ -27,6 +27,7 @@
 from sas.qtgui.Utilities.ModelEditors.TabbedEditor.TabbedModelEditor import TabbedModelEditor
 from sas.sascalc.calculator import sas_gen
 from sas.sascalc.calculator.geni import create_beta_plot, f_of_q, radius_of_gyration
+from sas.sascalc.calculator.sas_gen import ComputationType
 from sas.system.user import find_plugins_dir
 
 # Local UI
@@ -89,7 +90,7 @@ def __init__(self, parent=None):
         self.setup_display()
 
         # combox box
-        self.cbOptionsCalc.currentIndexChanged.connect(self.change_is_avg)
+        self.cbOptionsCalc.currentIndexChanged.connect(self.change_computation_type)
         # prevent layout shifting when widget hidden
         # TODO: Is there a way to lcoate this policy in the ui file?
         sizePolicy = self.cbOptionsCalc.sizePolicy()
@@ -621,7 +622,7 @@ def update_cbOptionsCalc_visibility(self):
         self.cbOptionsCalc.setVisible(allow)
         if (allow):
             # A helper function to set up the averaging system
-            self.change_is_avg()
+            self.change_computation_type()
         else:
             # If magnetic data present then no averaging is allowed
             self.is_avg = False
@@ -633,7 +634,7 @@ def update_cbOptionsCalc_visibility(self):
             self.checkboxLogSpace.setEnabled(not self.is_mag)
 
 
-    def change_is_avg(self):
+    def change_computation_type(self):
         """Adjusts the GUI for whether 1D averaging is enabled
 
         If the user has chosen to carry out Debye full averaging then the magnetic sld
@@ -658,6 +659,22 @@ def change_is_avg(self):
         self.checkboxLogSpace.setEnabled(self.is_avg)
         self.checkboxPluginModel.setEnabled(self.is_avg)
 
+        # set the type of calculation
+        self.model.set_computation_type(ComputationType(self.cbOptionsCalc.currentIndex()))
+        match self.cbOptionsCalc.currentIndex():
+            case 0 | 1 | 2:
+                pass
+            case 3:
+                self.checkboxPluginModel.setEnabled(False)
+                self.checkboxPluginModel.setChecked(True)
+                self.txtFileName.setText("saxs_fitting")
+                self.txtFileName.setEnabled(False)
+                self.cmdCompute.setText("Generate plugin model")
+                return
+            case _:
+                raise RuntimeError(f"Unknown computation type selected: {self.cbOptionsCalc.currentIndex()}")
+
+        self.cmdCompute.setText("Compute")
         if self.is_avg:
             self.txtMx.setText("0.0")
             self.txtMy.setText("0.0")
@@ -704,13 +721,20 @@ def loadFile(self):
             load_nuc = self.sender() == self.cmdNucLoad
             # request a file from the user
             if load_nuc:
-                f_type = """
-                    All supported files (*.SLD *.sld *.pdb *.PDB, *.vtk, *.VTK);;
-                        SLD files (*.SLD *.sld);;
-                        PDB files (*.pdb *.PDB);;
-                        VTK files (*.vtk *.VTK);;
-                        All files (*.*)
-                """
+                if self.model.type is ComputationType.SAXS:
+                    f_type = """
+                        All supported files (*.CIF *.cif *.pdb *.PDB);;
+                            CIF files (*.CIF *.cif);;
+                            PDB files (*.pdb *.PDB);;
+                    """
+                else:
+                    f_type = """
+                        All supported files (*.SLD *.sld *.pdb *.PDB, *.vtk, *.VTK);;
+                            SLD files (*.SLD *.sld);;
+                            PDB files (*.pdb *.PDB);;
+                            VTK files (*.vtk *.VTK);;
+                            All files (*.*)
+                    """
             else:
                 f_type = """
                     All supported files (*.OMF *.omf *.SLD *.sld, *.vtk, *.VTK);;
@@ -720,6 +744,11 @@ def loadFile(self):
                         All files (*.*)
                 """
             self.datafile = QtWidgets.QFileDialog.getOpenFileName(self, "Choose a file", "", f_type)[0]
+
+            if self.model.type is ComputationType.SAXS:
+                self.txtNucData.setText(os.path.basename(str(self.datafile)))
+                return
+
             # If a file has been sucessfully chosen
             if self.datafile:
                 # set basic data about the file
@@ -1412,6 +1441,49 @@ def onCompute(self):
 
         Copied from previous version
         """
+
+        if self.model.type is ComputationType.SAXS:
+            if self.datafile is None:
+                raise RuntimeError("No structure file is loaded! SAXS calculations require a structure file.")
+            from sas.qtgui.Calculators.SAXSPluginModelGenerator import get_base_plugin_name, write_plugin_model
+            write_plugin_model(self.datafile)
+            self.manager.communicator().customModelDirectoryChanged.emit() # notify that a new plugin model is available
+
+            # try to bring the fit panel into focus and select the newly generated plugin
+            try:
+                self.manager.actionFitting() # switch to fitting window
+                per = self.manager.perspective() # internal access into the fitting window's state
+                # currentFittingWidget is provided by the Fitting perspective
+                fw = getattr(per, 'currentFittingWidget', None)
+                if fw is not None:
+                    # select the plugin models category & our newly generated model
+                    idx = fw.cbCategory.findText("Plugin Models")
+                    if idx == -1: return
+
+                    # force population of model combobox
+                    fw.cbCategory.setCurrentIndex(idx)
+                    fw.onSelectCategory()
+
+                    # plugin name base is 'SAXS fit'
+                    # the actual model name includes a structure tag, e.g. 'SAXS fit (2epe)'
+                    model_name = get_base_plugin_name()
+                    midx = fw.cbModel.findText(model_name, QtCore.Qt.MatchStartsWith)
+                    if midx == -1: return
+
+                    # load the model into the parameter table
+                    fw.cbModel.setCurrentIndex(midx)
+                    fw.onSelectModel()
+
+                    # make sure the perspective window is visible and focused
+                    self.close() # close the calculator window to highlight the changes to the fitting window
+                    per.show()
+
+            except Exception:
+                logger.warning("Failed to bring the newly generated plugin model into focus. Please report this issue.")
+                pass
+
+            return
+
         try:
             # create the combined sld data and update from gui
             sld_data = self.create_full_sld_data()

src/sas/qtgui/Calculators/SAXSPluginModelGenerator.py

@@ -0,0 +1,74 @@
+from pathlib import Path
+
+from sas.system.user import find_plugins_dir
+
+
+def get_base_plugin_name() -> str:
+    """
+    Get the base name for the AUSAXS SAXS plugin model.
+
+    :return: The base name of the plugin model.
+    """
+
+    return "SAXS fit"
+
+def write_plugin_model(structure_path: str):
+    """
+    Write the AUSAXS SAXS plugin model to the plugins directory.
+    The current version will be overwritten if it exists.
+
+    :param structure_path: Path to the structure file to be used by the plugin.
+    """
+
+    path = Path(find_plugins_dir()) / "ausaxs_saxs_plugin.py"
+    text = get_model_text(structure_path)
+    with open(path, 'w') as f:
+        f.write(text)
+
+def get_model_text(structure_path: str) -> str:
+    """
+    Generate the text of the AUSAXS SAXS plugin model.
+
+    :param structure_path: Path to the structure file to be used by the plugin.
+    :return: The text of the plugin model.
+    """
+
+    return (
+f'''\
+r"""
+This file is auto-generated, and any changes will be overwritten. 
+
+This plugin model uses the AUSAXS library (https://doi.org/10.1107/S160057672500562X) to fit the provided SAXS data to the file:
+ * \"{structure_path}\"
+If this is not the intended structure file, please regenerate the plugin model from the generic scattering calculator.
+"""
+'''
+
+f'''\
+name = "{get_base_plugin_name()} ({Path(structure_path).name.split('.')[0]})"
+title = "AUSAXS"
+description = "Structural validation using AUSAXS"
+category = "plugin"
+parameters = [
+    # name, units, default, [min, max], type, description
+    ['c', '', 1, [0, 100], '', 'Solvent density'],
+    #['d', '', 1, [0, 2], '', 'Excluded volume parameter']
+]
+
+###
+import pyausaxs as ausaxs
+ausaxs.settings.set(\"allow_unknown_atoms\", \"false\")
+ausaxs.settings.set(\"allow_unknown_residues\", \"false\")
+
+structure_path = "{str(Path(structure_path).as_posix())}"
+
+def Iq(q, c):
+    # Initialize on first call to keep objects alive for function lifetime
+    if not hasattr(Iq, '_initialized'):
+        Iq._mol = ausaxs.create_molecule(structure_path)
+        Iq._mol.hydrate()
+        Iq._fitobj = ausaxs.manual_fit(Iq._mol)
+        Iq._initialized = True
+    return Iq._fitobj.evaluate([c], q)
+Iq.vectorized = True
+''')

src/sas/qtgui/Calculators/UI/GenericScatteringCalculator.ui

@@ -2194,6 +2194,11 @@ NOTE: Currently not impacted by Solvent SLD </string>
        <string>Debye full avg. w/ β(Q)</string>
       </property>
      </item>
+     <item>
+      <property name="text">
+       <string>SAXS fitting</string>
+      </property>
+     </item>
     </widget>
    </item>
   </layout>

src/sas/sascalc/calculator/ausaxs/ausaxs_sans_debye.py

@@ -1,6 +1,6 @@
 import logging
 
-from pyausaxs import AUSAXS
+import pyausaxs as ausaxs
 
 from sas.sascalc.calculator.ausaxs.sasview_sans_debye import sasview_sans_debye
 
@@ -14,8 +14,7 @@ def evaluate_sans_debye(q, coords, w):
     *w* is the weight associated with each point.
     """
     try:
-        ausaxs = AUSAXS()
-        Iq = ausaxs.debye(q, coords[0,:], coords[1,:], coords[2,:], w)
+        Iq = ausaxs.sasview.debye_no_ff(q, coords[0,:], coords[1,:], coords[2,:], w)
         return Iq
     except Exception as e:
         logging.warning("AUSAXS Debye calculation failed: %s. Falling back to default implementation.", e)

src/sas/sascalc/calculator/sas_gen.py

@@ -9,6 +9,7 @@
 import logging
 import os
 import sys
+from enum import Enum
 
 import numpy as np
 from periodictable import formula, nsf
@@ -53,6 +54,12 @@ def transform_center(pos_x, pos_y, pos_z):
     posz = pos_z - (min(pos_z) + max(pos_z)) / 2.0
     return posx, posy, posz
 
+class ComputationType(Enum):
+    SANS_2D = 0
+    SANS_1D = 1
+    SANS_1D_BETA = 2
+    SAXS = 3
+
 class GenSAS:
     """
     Generic SAS computation Model based on sld (n & m) arrays
@@ -75,6 +82,7 @@ def __init__(self):
         self.data_vol = None # [A^3]
         self.is_avg = False
         self.is_elements = False
+        self.type = ComputationType.SANS_2D
         ## Name of the model
         self.name = "GenSAS"
         ## Define parameters
@@ -106,6 +114,12 @@ def __init__(self):
         # fixed parameters
         self.fixed = []
 
+    def set_computation_type(self, computation_type : ComputationType):
+        """
+        Set the computation type. This will determine which calculation is performed.
+        """
+        self.type = computation_type
+
     def set_pixel_volumes(self, volume):
         """
         Set the volume of a pixel in (A^3) unit
@@ -180,33 +194,41 @@ def calculate_Iq(self, qx, qy=None):
         x, y, z = self.transform_positions()
         sld = self.data_sldn - self.params['solvent_SLD']
         vol = self.data_vol
-        if qy is not None and len(qy) > 0:
-            # 2-D calculation
-            qx, qy = _vec(qx), _vec(qy)
-            # MagSLD can have sld_m = None, although in practice usually a zero array
-            # if all are None can continue as normal, otherwise set None to array of zeroes to allow rotations
-            mx, my, mz = self.transform_magnetic_slds()
-            in_spin = self.params['Up_frac_in']
-            out_spin = self.params['Up_frac_out']
-            # transform angles from environment to beamline coords
-            s_theta, s_phi = self.transform_angles()
-
-            if self.is_elements:
-                I_out = Iqxy(
-                    qx, qy, x, y, z, sld, vol, mx, my, mz,
-                    in_spin, out_spin, s_theta, s_phi,
-                    self.data_elements, self.is_elements)
-            else:
-                I_out = Iqxy(
-                    qx, qy, x, y, z, sld, vol, mx, my, mz,
-                    in_spin, out_spin, s_theta, s_phi,
-                    )
-        else:
-            # 1-D calculation
-            q = _vec(qx)
-            if self.is_avg:
-                x, y, z = transform_center(x, y, z)
-            I_out = Iq(q, x, y, z, sld, vol, is_avg=self.is_avg)
+        match self.type:
+            case ComputationType.SANS_2D:
+                if not (qy is not None and len(qy) > 0):
+                    raise ValueError("For a SANS_2D computation, qy cannot be None or empty")
+
+                # 2-D calculation
+                qx, qy = _vec(qx), _vec(qy)
+                # MagSLD can have sld_m = None, although in practice usually a zero array
+                # if all are None can continue as normal, otherwise set None to array of zeroes to allow rotations
+                mx, my, mz = self.transform_magnetic_slds()
+                in_spin = self.params['Up_frac_in']
+                out_spin = self.params['Up_frac_out']
+                # transform angles from environment to beamline coords
+                s_theta, s_phi = self.transform_angles()
+
+                if self.is_elements:
+                    I_out = Iqxy(
+                        qx, qy, x, y, z, sld, vol, mx, my, mz,
+                        in_spin, out_spin, s_theta, s_phi,
+                        self.data_elements, self.is_elements)
+                else:
+                    I_out = Iqxy(
+                        qx, qy, x, y, z, sld, vol, mx, my, mz,
+                        in_spin, out_spin, s_theta, s_phi,
+                        )
+
+            case ComputationType.SANS_1D | ComputationType.SANS_1D_BETA:
+                # 1-D calculation
+                q = _vec(qx)
+                if self.is_avg:
+                    x, y, z = transform_center(x, y, z)
+                I_out = Iq(q, x, y, z, sld, vol, is_avg=self.is_avg)
+
+            case ComputationType.SAXS:
+                raise RuntimeError("SAXS calculations can only be performed through a plugin model! Please click the \"plugin model\" button instead.")
 
         vol_correction = self.data_total_volume / self.params['total_volume']
         result = ((self.params['scale'] * vol_correction) * I_out
@@ -266,6 +288,7 @@ def run(self, x=0.0):
             if len(x[1]) > 0:
                 raise ValueError("Not a 1D vector.")
             # 1D I is found at y=0 in the 2D pattern
+            self.set_computation_type(ComputationType.SANS_1D)
             out = self.calculate_Iq(x[0])
             return out
         else:

test/sascalculator/utest_sas_gen.py

@@ -226,15 +226,15 @@ def test_debye_impl(self):
         """
         Test that the Debye algorithm supplied by the external AUSAXS library agrees with the default implementation.
         """
-        from pyausaxs import AUSAXS
+        import pyausaxs as ausaxs
 
         from sas.sascalc.calculator.ausaxs import ausaxs_sans_debye, sasview_sans_debye
 
         rng = np.random.default_rng(1984)
-        ausaxs = AUSAXS()
 
         # ensure the library is available and ready to run on all CI systems
-        assert ausaxs.ready(), "AUSAXS library not available, test cannot be run."
+        # this awkward syntax will be improved in a future version of pyausaxs ...
+        assert ausaxs.wrapper.AUSAXS.AUSAXS().ready(), "AUSAXS library not available, test cannot be run."
 
         # get all pdb files in the data folder
         import glob