Merge pull request #76 from rs-station/dials-v-3-22-1

Merge serial + latest DIALS-compatible branch into main

README.rst

@@ -71,6 +71,12 @@ Note that you need to import the image data using ``dials.import``. For informat
 
 If any issues occur either with installation or use of the software, please file an issue at `issue tracker`_. Any and all questions, concerns, or feature requests are welcome.
 
+============
+Publications
+============
+
+`Laue-DIALS: Open-source software for polychromatic x-ray diffraction data <https://doi.org/10.1063/4.0000265>`_. **Hewitt, Rick A.** and Kevin M. Dalton, Derek Mendez, Harrison K. Wang, Margaret A. Klureza, Dennis E. Brookner, Jack B. Greisman, David McDonagh, Vukica Šrajer, Nicholas K. Sauter, Aaron S. Brewster, Doeke R. Hekstra. *Structural Dynamics, 2024*
+
 .. _careless: https://github.com/rs-station/careless
 .. _DIALS: https://dials.github.io/index.html
 .. _issue tracker: https://github.com/rs-station/laue-dials/issues

docs/api/laue_dials.command_line.rst

@@ -77,3 +77,11 @@ laue\_dials.command\_line.sequence\_to\_stills module
    :members:
    :undoc-members:
    :show-inheritance:
+
+laue\_dials.command\_line.version module
+----------------------------------------
+
+.. automodule:: laue_dials.command_line.version
+   :members:
+   :undoc-members:
+   :show-inheritance:

docs/cli/functions.md

@@ -6,6 +6,7 @@
 .. toctree::
    :maxdepth: 1
 
+   version
    find_spots
    index
    sequence_to_stills

docs/cli/version.rst

@@ -0,0 +1,9 @@
+.. _version:
+
+laue_dials.version
+==================
+
+Introduction
+------------
+
+.. python_string:: laue_dials.command_line.version.help_message

docs/requirements.txt

@@ -1,7 +1,7 @@
 # Requirements file for ReadTheDocs, check .readthedocs.yml.
 # To build the module reference correctly, make sure every external package
 # under `install_requires` in `setup.cfg` is also listed here!
-dials
+dials>=3.22.1
 sphinx
 myst-parser
 sphinxcontrib-autoprogram

docs/tutorials/tutorial_README.rst

@@ -1,11 +1,12 @@
-==========
+==========================
 Jupyter Notebook Tutorials
-==========
+==========================
 
 These tutorials can help guide you through analyzing a data set with Laue-DIALS. The data for
 these tutorials are stored on SBGrid, and can be downloaded using the following in any bash terminal:
 
 .. code:: bash
+
    # HEWL
    rsync -av rsync://data.sbgrid.org/10.15785/SBGRID/1118 .
    cd 1118 ; shasum -c files.sha

examples/serial_pipeline.sh

@@ -0,0 +1,80 @@
+# The values in this script are for the DHFR data at
+# https://zenodo.org/records/10199220. You may need
+# to change these values if analyzing a different data set.
+
+FILE_INPUT_TEMPLATE="/PATH/TO/DATA"
+MASK="/PATH/TO/MASK/FILE"
+N=8
+
+# Import data
+# You may need to change some of these values to match your data set
+dials.import geometry.scan.oscillation=0,1 \
+    geometry.goniometer.axes=0,1,0 \
+    convert_sequences_to_stills=True \
+    geometry.beam.wavelength=1.04 \
+    geometry.detector.panel.pixel=0.08854,0.08854 \
+    lookup.mask=$MASK \
+    input.template=$FILE_INPUT_TEMPLATE
+
+# Get a monochromatic geometry model
+dials.find_spots imported.expt \
+    output.shoeboxes=False \
+    spotfinder.mp.nproc=$N \
+    spotfinder.threshold.dispersion.gain=0.15 \
+    spotfinder.filter.max_separation=10 \
+    lookup.mask=$MASK
+
+# Get a monochromatic geometry model
+laue.index imported.expt strong.refl \
+    indexer.indexing.nproc=$N \
+    indexer.indexing.method="pink_indexer" \
+    indexer.indexing.pink_indexer.wavelength=1.1 \
+    indexer.indexing.pink_indexer.percent_bandwidth=15 \
+    indexer.indexing.pink_indexer.max_refls=50 \
+    indexer.indexing.pink_indexer.min_lattices=5 \
+    indexer.indexing.pink_indexer.rotogram_grid_points=360 \
+    indexer.indexing.pink_indexer.voxel_grid_points=250 \
+    indexer.indexing.known_symmetry.space_group=19 \
+    indexer.indexing.known_symmetry.unit_cell=34.297,45.552,99.035,90,90,90 \
+    indexer.indexing.refinement_protocol.mode=None \
+    indexer.indexing.stills.ewald_proximal_volume_max=0.03 \
+    indexer.indexing.stills.rmsd_min_px=20 \
+    indexer.indexing.stills.refine_all_candidates=True \
+    indexer.indexing.joint_indexing=False \
+    indexer.refinement.reflections.outlier.algorithm=None \
+    indexer.refinement.parameterisation.auto_reduction.action=fix \
+    indexer.refinement.parameterisation.scan_varying=False \
+    laue_output.index_only=True
+
+# Polychromatic pipeline
+laue.optimize_indexing monochromatic.expt monochromatic.refl \
+    output.experiments="optimized.expt" \
+    output.reflections="optimized.refl" \
+    output.log="laue.optimize_indexing.log" \
+    wavelengths.lam_min=1.00 \
+    wavelengths.lam_max=1.20 \
+    reciprocal_grid.d_min=2.0 \
+    geometry.unit_cell=34.297,45.552,99.035,90,90,90 \
+    n_macrocycles=5 \
+    keep_unindexed=False \
+    filter_spectrum=True \
+    nproc=$N
+
+laue.refine optimized.* \
+    output.experiments="poly_refined.expt" \
+    output.reflections="poly_refined.refl" \
+    output.log="laue.poly_refined.log" \
+    nproc=$N
+
+laue.predict poly_refined.* \
+    output.reflections="predicted.refl" \
+    output.log="laue.predict.log" \
+    wavelengths.lam_min=1.00 \
+    wavelengths.lam_max=1.20 \
+    reciprocal_grid.d_min=1.4 \
+    nproc=$N
+
+laue.integrate poly_refined.expt predicted.refl \
+    output.filename="integrated.mtz" \
+    output.log="laue.integrate.log" \
+    nproc=$N

examples/tutorial.ipynb

@@ -41,7 +41,7 @@
    "source": [
     "# Importing Data\n",
     "\n",
-    "We can use `dials.import` as a way to import the data files written at experimental facilities into a format that is friendly to both `DIALS` and `laue-dials`. Feel free to use any data set you'd like below, but a [sample lysozyme data set](https://zenodo.org/record/6407157) has been uploaded to zenodo for your convenience, and this notebook has been tested using that dataset.\n",
+    "We can use `dials.import` as a way to import the data files written at experimental facilities into a format that is friendly to both `DIALS` and `laue-dials`. Feel free to use any data set you'd like below, but a [sample DHFR data set](https://zenodo.org/record/6407157) has been uploaded to zenodo for your convenience, and this notebook has been tested using that dataset.\n",
     "\n",
     "### Masks\n",
     "\n",

setup.cfg

@@ -72,6 +72,7 @@ testing =
 # console_scripts =
 #     script_name = laue_dials.module:function
 console_scripts =
+    laue.version = laue_dials.command_line.version:run
     laue.find_spots = laue_dials.command_line.find_spots:run
     laue.index = laue_dials.command_line.index:run
     laue.sequence_to_stills = laue_dials.command_line.sequence_to_stills:run

setup.py

@@ -1,6 +1,6 @@
 """
-    Setup file for laue_dials.
-    Use setup.cfg to configure.
+Setup file for laue_dials.
+Use setup.cfg to configure.
 """
 
 from setuptools import setup

src/laue_dials/algorithms/laue.py

@@ -186,13 +186,16 @@ class LaueAssigner(LaueBase):
     An object to assign miller indices to a laue still
     """
 
-    def __init__(self, s0, s1, cell, R, lam_min, lam_max, dmin, spacegroup="1"):
+    def __init__(
+        self, s0, s1, cell, R, lam_min, lam_max, lam_peak, dmin, spacegroup="1"
+    ):
         """
         Parameters:
             s1 (np.ndarray): n x 3 array indicating the direction of the scattered beam wavevector.
         """
         super().__init__(s0, cell, R, lam_min, lam_max, dmin, spacegroup)
 
+        self.lam_peak = lam_peak
         self._s1 = s1 / np.linalg.norm(s1, axis=-1)[:, None]
         self._qobs = self._s1 - self.s0
         self._qpred = np.zeros_like(self._s1)
@@ -349,11 +352,30 @@ def assign(self):
         _, idx, counts = np.unique(
             Raypred, return_index=True, return_counts=True, axis=0
         )
-        harmonics = counts > 1
-
-        # Remove harmonics from the feasible set
-        Hall = Hall[idx]
-        qall = qall[idx]
+        unique_rays, counts = np.unique(Raypred, return_counts=True, axis=0)
+        harmonic_rays = counts > 1
+
+        # Keep only harmonics closest to peak wavelength in feasible set
+        harmonics = np.zeros(len(Hall), dtype=bool)
+        to_keep = np.ones(len(Hall), dtype=bool)
+        for i, ray in enumerate(unique_rays[harmonic_rays]):
+            idh = np.where((Raypred == ray).all(axis=1))[0]
+            if len(idh) == 1:  # Not a harmonic
+                to_keep[idh] = True
+            else:
+                harmonics[idh] = True
+                qharm = qall[idh]
+                harmonic_wavelengths = (
+                    -2.0 * (self.s0 * qharm).sum(-1) / (qharm * qharm).sum(-1)
+                )
+                kept_harmonic = np.argmin(np.abs(harmonic_wavelengths - self.lam_peak))
+                to_keep[idh] = False
+                to_keep[idh[kept_harmonic]] = True
+
+        # Remove non-optimal harmonics from the feasible set
+        Hall = Hall[to_keep]
+        qall = qall[to_keep]
+        harmonics = harmonics[to_keep]
 
         dmat = rs.utils.angle_between(self.qobs[..., None, :], qall[None, ..., :])
         cost = dmat

src/laue_dials/command_line/compute_rmsds.py

@@ -19,9 +19,6 @@
 
 from laue_dials.utils.version import laue_version
 
-# Print laue-dials + DIALS versions
-laue_version()
-
 logger = logging.getLogger("laue-dials.command_line.compute_rmsds")
 
 help_message = """
@@ -47,6 +44,10 @@
     .type = str
     .help = "Save a CSV of the RMSDs per image with this filename."
 
+  image = None
+    .type = int
+    .help = "Show a histogram of residuals for this image."
+
   output = "residuals.png"
     .type = str
     .help = "The filename for the generated plot."
@@ -126,6 +127,9 @@ def run(args=None, *, phil=working_phil):
     xfel_logger.setLevel(loglevel)
     fh.setLevel(loglevel)
 
+    # Print version information
+    logger.info(laue_version())
+
     # Print help if no input
     if not params.input.experiments or not params.input.reflections:
         parser.print_help()
@@ -184,12 +188,20 @@ def run(args=None, *, phil=working_phil):
     y_resids = data["ycal"] - data["yobs"]
     sqr_resids = x_resids**2 + y_resids**2
     mean_resids = np.zeros(len(images))
+    n_refls = np.zeros(len(images))
     for i, img_num in enumerate(images):
         sel = data["image"] == img_num
         mean_resids[i] = np.mean(sqr_resids[sel])
+        n_refls[i] = sum(sel)
+        if params.image == img_num:
+            plt.hist(np.sqrt(sqr_resids[sel]))
+            plt.xlabel("Residuals (px)")
+            plt.xlabel("# Reflections")
+            plt.title(f"Residuals for Image {img_num}")
+
     rmsds = np.sqrt(mean_resids)
 
-    resid_data = pd.DataFrame({"Image": images, "RMSD (px)": rmsds})
+    resid_data = pd.DataFrame({"Image": images, "RMSD (px)": rmsds, "N Spots": n_refls})
 
     pd.set_option("display.max_rows", None)
     logger.info(f"RMSDs per image: \n{resid_data}")

src/laue_dials/command_line/find_spots.py

@@ -13,9 +13,6 @@
 from laue_dials.algorithms.monochromatic import find_spots
 from laue_dials.utils.version import laue_version
 
-# Print laue-dials + DIALS versions
-laue_version()
-
 logger = logging.getLogger("laue-dials.command_line.find_spots")
 
 help_message = """
@@ -118,6 +115,9 @@ def run(args=None, *, phil=working_phil):
     xfel_logger.setLevel(loglevel)
     fh.setLevel(loglevel)
 
+    # Print version information
+    logger.info(laue_version())
+
     # Log diff phil
     diff_phil = parser.diff_phil.as_str()
     if diff_phil != "":

src/laue_dials/command_line/index.py

@@ -15,9 +15,6 @@
                                                  scan_varying_refine)
 from laue_dials.utils.version import laue_version
 
-# Print laue-dials + DIALS versions
-laue_version()
-
 logger = logging.getLogger("laue-dials.command_line.index")
 
 help_message = """
@@ -218,6 +215,9 @@ def run(args=None, *, phil=working_phil):
     xfel_logger.setLevel(loglevel)
     fh.setLevel(loglevel)
 
+    # Print version information
+    logger.info(laue_version())
+
     # Log diff phil
     diff_phil = parser.diff_phil.as_str()
     if diff_phil != "":
@@ -234,6 +234,10 @@ def run(args=None, *, phil=working_phil):
         params.input.reflections, params.input.experiments
     )
 
+    # Remove extraneous data from params
+    params.input.reflections = None
+    params.input.experiments = None
+
     # Get initial time for process
     start_time = time.time()
 

src/laue_dials/command_line/integrate.py

@@ -23,9 +23,6 @@
 from laue_dials.algorithms.integration import SegmentedImage
 from laue_dials.utils.version import laue_version
 
-# Print laue-dials + DIALS versions
-laue_version()
-
 logger = logging.getLogger("laue-dials.command_line.integrate")
 
 help_message = """
@@ -191,6 +188,9 @@ def run(args=None, *, phil=working_phil):
     xfel_logger.setLevel(loglevel)
     fh.setLevel(loglevel)
 
+    # Print version information
+    logger.info(laue_version())
+
     # Log diff phil
     diff_phil = parser.diff_phil.as_str()
     if diff_phil != "":
@@ -208,6 +208,10 @@ def run(args=None, *, phil=working_phil):
     )
     preds = reflections[0]  # Get predictions
 
+    # Remove duplicate expt + refl data
+    params.input.experiments = None
+    params.input.reflections = None
+
     # Sanity checks
     if len(expts) == 0:
         parser.print_help()
@@ -226,8 +230,11 @@ def run(args=None, *, phil=working_phil):
     # Multiprocess integration
     num_processes = params.nproc
     logger.info("Starting integration.")
-    with Pool(processes=num_processes) as pool:
-        refls_arr = pool.starmap(integrate_image, inputs, chunksize=1)
+    if num_processes == 1:
+        refls_arr = [integrate_image(*i) for i in inputs]
+    else:
+        with Pool(processes=num_processes) as pool:
+            refls_arr = pool.starmap(integrate_image, inputs, chunksize=1)
     logger.info("Integration finished.")
 
     # Construct an integrated reflection table