From d9639a883a56c1435a0cca9f532a374a7509e64c Mon Sep 17 00:00:00 2001 From: Ben Williams Date: Wed, 26 Oct 2022 14:06:12 +0100 Subject: [PATCH 01/33] Raze the old stuff --- screen19/__init__.py | 200 ------ screen19/minimum_exposure.py | 493 -------------- screen19/screen.py | 1192 ---------------------------------- tests/conftest.py | 9 - tests/test_screen19.py | 89 --- 5 files changed, 1983 deletions(-) delete mode 100644 screen19/__init__.py delete mode 100644 screen19/minimum_exposure.py delete mode 100644 screen19/screen.py delete mode 100644 tests/conftest.py delete mode 100644 tests/test_screen19.py diff --git a/screen19/__init__.py b/screen19/__init__.py deleted file mode 100644 index 2d957fc..0000000 --- a/screen19/__init__.py +++ /dev/null @@ -1,200 +0,0 @@ -"""Common tools for the I19 module.""" - -import logging -import os -import re -import subprocess -import sys -import traceback -from typing import Dict, Tuple # noqa: F401 - -import procrunner - -# Flake8 does not detect typing yet (https://gitlab.com/pycqa/flake8/issues/342) - -__version__ = "0.213" - -logger = logging.getLogger("dials.screen19") -debug, info, warn = logger.debug, logger.info, logger.warning - - -# Set axis tick positions manually. Accounts for reciprocal(-square) d-scaling. -d_ticks = [5, 3, 2, 1.5, 1, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4] - - -def terminal_size() -> Tuple[int, int]: - """ - Find the current size of the terminal window. - - :return: Number of columns; number of rows. - """ - columns, rows = 80, 25 - if sys.stdout.isatty(): - try: - result = procrunner.run( - ["stty", "size"], - timeout=1, - raise_timeout_exception=True, - print_stdout=False, - print_stderr=False, - ) - rows, columns = (int(i) for i in result.stdout.decode("utf-8").split()) - except Exception: # ignore any errors and use default size - pass # FIXME: Can we be more specific about the type of exception? - columns = min(columns, 120) - rows = min(rows, int(columns / 3)) - - return columns, rows - - -def prettyprint_dictionary(d): - """ - Produce a nice string representation of a dictionary, for printing. - - :param d: Dictionary to be printed. - :type d: Dict[Optional[Any]] - :return: String representation of :param d:. - :rtype: str - """ - return "{\n%s\n}" % "\n".join( - " %s: %s" - % ( - k, - str(v.decode("latin-1") if isinstance(v, bytes) else v).replace( - "\n", "\n%s" % (" " * (4 + len(k))) - ), - ) - for k, v in d.items() - ) - - -def prettyprint_procrunner(d): - """ - Produce a nice string representation of a procrunner ReturnObject, for printing. - - :param d: ReturnObject to be printed. - :return: String representation of :param d:. - :rtype: str - """ - return prettyprint_dictionary( - { - "command": d.args, - "exitcode": d.returncode, - "stdout": d.stdout, - "stderr": d.stderr, - } - ) - - -def make_template(f): - """ - Generate a xia2-style filename template. - - From a given filename, generate a template filename by substituting a hash - character (#) for each numeral in the last contiguous group of numerals - before the file extension. - For example, the filename example_01_0001.cbf becomes example_01_####.cbf. - - :param f: Filename, with extension. - :type f: str - :return: Filename template, with extension; image number. - :rtype: Tuple(str, int) - """ - # Split the file from its path - directory, f = os.path.split(f) - # Split off the file extension, assuming it begins at the first full stop, - # also split the last contiguous group of digits off the filename root - parts = re.split(r"([0-9#]+)(?=\.\w)", f, 1) - # Get the number of digits in the group we just isolated and their value - try: - # Combine the root, a hash for each digit and the extension - length = len(parts[1]) - template = parts[0] + "#" * length + parts[2] - image = int(parts[1].replace("#", "0")) - except IndexError: - template = parts[0] - image = None - return os.path.join(directory, template), image - - -def plot_intensities( - bins, - hist_value_factor, - title="'Pixel intensity distribution'", - xlabel="'% of maximum'", - ylabel="'Number of pixels'", - xticks="", - style="with boxes", -): - """ - Create an ASCII art histogram of intensities. - - :param bins: - :param hist_value_factor: - :param title: - :param xlabel: - :param ylabel: - :param xticks: - :param style: - """ - columns, rows = terminal_size() - - command = ["gnuplot"] - plot_commands = [ - "set term dumb %d %d" % (columns, rows - 2), - "set title %s" % title, - "set xlabel %s" % xlabel, - "set ylabel %s offset character %d,0" % (ylabel, len(ylabel) // 2), - "set logscale y", - "set boxwidth %f" % hist_value_factor, - "set xtics %s out nomirror" % xticks, - "set ytics out", - "plot '-' using 1:2 title '' %s" % style, - ] - for x in sorted(bins.keys()): - plot_commands.append("%f %d" % (x * hist_value_factor, bins[x])) - plot_commands.append("e") - - debug("running %s with:\n %s\n", " ".join(command), "\n ".join(plot_commands)) - - try: - result = procrunner.run( - command, - stdin="\n".join(plot_commands).encode("utf-8") + b"\n", - timeout=120, - raise_timeout_exception=True, - print_stdout=False, - print_stderr=False, - environment_override={"LD_LIBRARY_PATH": ""}, - ) - except (OSError, subprocess.TimeoutExpired): - info(traceback.format_exc()) - warn( - "Error running gnuplot. Cannot plot intensity distribution. " - "No exit code." - ) - return - else: - debug("result = %s", prettyprint_procrunner(result)) - - returncode = getattr(result, "returncode") - if returncode: - warn( - "Error running gnuplot. Cannot plot intensity distribution. " - "Exit code %d", - returncode, - ) - else: - star = re.compile(r"\*") - state = set() - for line in result.stdout.decode("utf-8").split("\n"): - if line.strip() != "": - stars = {m.start(0) for m in re.finditer(star, line)} - if not stars: - state = set() - else: - state |= stars - line = list(line) - for s in state: - line[s] = "*" - info("".join(line)) diff --git a/screen19/minimum_exposure.py b/screen19/minimum_exposure.py deleted file mode 100644 index 406d794..0000000 --- a/screen19/minimum_exposure.py +++ /dev/null @@ -1,493 +0,0 @@ -""" -Perform straight-line Wilson plot fit. Draw the Wilson plot. - -Reflection d-spacings are determined from the crystal symmetry (from -indexing) and the Miller indices of the indexed reflections. The -atomic displacement parameter is assumed isotropic. Its value is -determined from a fit to the reflection data: - I = A * exp(-B / (2 * d^2)), -where I is the intensity and the scale factor, A, and isotropic -displacement parameter, B, are the fitted parameters. - -An I/σ condition for 'good' diffraction statistics is set by the -instance variable min_i_over_sigma, and the user's desired -resolution is set by the instance variable desired_d. A crude -error model is assumed, whereby σ² = I, and so the I/σ condition -translates trivially to a threshold I. - -The value of the fitted intensity function at the desired -resolution is compared with the threshold I. The ratio of these -values is used to determine a recommended exposure (flux × exposure time) -for the full data collection. - -The Wilson plot of I as a function of d is drawn as the file -'wilson_plot.png'. The plot can optionally be saved in other formats. - -Examples: - - screen19.minimum_exposure integrated.expt integrated.refl - - screen19.minimum_exposure indexed.expt indexed.refl - - screen19.minimum_exposure min_i_over_sigma=2 desired_d=0.84 wilson_fit_max_d=4 \ - integrated.expt integrated.refl - -""" - -import logging -import time -from io import StringIO -from typing import Iterable, List, Optional, Sequence, Union - -import numpy as np -from scipy.optimize import curve_fit -from tabulate import tabulate - -import iotbx.phil -from cctbx import miller -from libtbx.phil import scope, scope_extract - -from dials.array_family import flex -from dials.util import log -from dials.util.options import OptionParser -from dials.util.version import dials_version -from dxtbx.model import ExperimentList - -from screen19 import __version__, d_ticks, plot_intensities, terminal_size - -# Custom types -FloatSequence = Sequence[float] -Fit = Union[np.ndarray, Iterable, int, float] - -help_message = __doc__ - - -phil_scope = iotbx.phil.parse( - """ - verbosity = 0 - .type = int(value_min=0) - .caption = 'Verbosity level of log output' - .help = "Possible values:\n" - "\t• 0: Info log output to stdout/logfile\n" - "\t• 1: Info & debug log output to stdout/logfile" - minimum_exposure - .caption = 'Parameters for the calculation of the lower exposure bound' - { - desired_d = None - .multiple = True - .type = float - .caption = u'Desired resolution limit, in Ångströms, of diffraction data' - .help = 'This is the resolution target for the lower-bound exposure ' \ - 'recommendation.' - min_i_over_sigma = 2 - .type = float - .caption = u'Target I/σ value for lower-bound exposure recommendation' - .help = 'The lower-bound exposure recommendation provides an estimate of ' \ - u'the exposure (flux × exposure time) required to ensure that the' \ - 'majority of expected reflections at the desired resolution limit' \ - u'have I/σ greater than or equal to this value.' - wilson_fit_max_d = 4 # Å - .type = float - .caption = u'Maximum d-value (in Ångströms) for displacement parameter fit' - .help = 'Reflections with lower resolution than this value will be ' \ - 'ignored for the purposes of the Wilson plot.' - } - output - .caption = 'Parameters to control the output' - { - log = 'screen19.minimum_exposure.log' - .type = str - .caption = 'Location for the info log' - debug_log = 'screen19.minimum_exposure.debug.log' - .type = str - .caption = 'Location for the debug log' - wilson_plot = 'wilson_plot' - .type = str - .caption = 'Filename for the Wilson plot image' - .help = "By default, the extension '.png' is appended. If you include " \ - "a different extension, either '.pdf', '.ps', '.eps' or '.svg', " \ - "a file of that format will be created instead." - } - """, - process_includes=True, -) - -logger_name = "dials.screen19.minimum_exposure" -logger = logging.getLogger(logger_name) -debug, info, warn = logger.debug, logger.info, logger.warning - - -def scaled_debye_waller(x: float, b: float, a: float) -> float: - """ - Calculate a scaled isotropic Debye-Waller factor. - - By assuming a single isotropic disorder parameter, :param:`b`, this factor - approximates the decay of diffracted X-ray intensity increasing resolution - (decreasing d, increasing sin(θ)). - - Args: - x: Equivalent to 1/d². - b: Isotropic displacement parameter. - a: A scale factor. - - Returns: - Estimated value of scaled isotropic Debye-Waller factor. - """ - return a * np.exp(-b / 2 * x) - - -def wilson_fit( - d_star_sq: FloatSequence, - intensity: FloatSequence, - sigma: FloatSequence, - wilson_fit_max_d: float, -) -> Fit: - """ - Fit a simple Debye-Waller factor, assume isotropic disorder parameter. - - Reflections with d ≥ :param:`wilson_fit_max_d` are ignored. - - Args: - d_star_sq: 1/d² (equivalently d*²), sequence of values for the observed - reflections (units of Å⁻² assumed). - intensity: Sequence of reflection intensities. - sigma: Sequence of uncertainties in reflection intensity. - wilson_fit_max_d: The minimum resolution for reflections against which to - fit. - - Returns: - - The fitted isotropic displacement parameter (units of Ų assumed); - - The fitted scale factor. - - """ - # Eliminate reflections with d > wilson_fit_max_d from the fit - sel = d_star_sq > 1 / wilson_fit_max_d**2 - - # Perform a weighted Wilson plot fit to the reflection intensities - fit, cov = curve_fit( - scaled_debye_waller, - d_star_sq.select(sel), - intensity.select(sel), - sigma=sigma.select(sel), - bounds=(0, np.inf), - ) - - return fit - - -def wilson_plot_ascii( - miller_array: miller.array, d_ticks: Optional[Sequence] = None -) -> None: - """ - Print an ASCII-art Wilson plot of reflection intensities. - - Equivalent reflections will be merged according to the crystal symmetry. - - Args: - miller_array: An array of integrated intensities, bundled with appropriate - crystal symmetry and unit cell info. - d_ticks: d location of ticks on 1/d² axis. - """ - # Draw the Wilson plot, using existing functionality in cctbx.miller - columns, rows = terminal_size() - n_bins = min(columns, miller_array.data().size()) - miller_array.setup_binner_counting_sorted(n_bins=n_bins, reflections_per_bin=1) - wilson = miller_array.wilson_plot(use_binning=True) - # Get the relevant plot data from the miller_array: - binned_intensity = [x if x else 0 for x in wilson.data[1:-1]] - bins = dict(zip(wilson.binner.bin_centers(1), binned_intensity)) - if d_ticks: - tick_positions = ", ".join([f'"{d:g}" {1 / d ** 2}' for d in d_ticks]) - tick_positions = tick_positions.join(["(", ")"]) - else: - tick_positions = "" - # Draw the plot: - plot_intensities( - bins, - 1, - title="'Wilson plot'", - xlabel="'d (Angstrom) (inverse-square scale)'", - ylabel="'I (counts)'", - xticks=tick_positions, - style="with lines", - ) - - -def wilson_plot_image( - d_star_sq: FloatSequence, - intensity: FloatSequence, - fit: Fit, - max_d: Optional[float] = None, - ticks: Optional[FloatSequence] = None, - output: str = "wilson_plot", -) -> None: - """ - Generate the Wilson plot as a PNG image. - - :param:`max_d` allows greying out of the reflections not included in the - isotropic Debye-Waller fit. - - Args: - d_star_sq: 1/d² values of reflections. - intensity: Intensities of reflections. - fit: Fitted parameters (tuple of fitted isotropic displacement parameter and - fitted scale factor). - max_d: The minimum resolution for reflections used in the Debye-Waller fit. - ticks: d location of ticks on 1/d² axis. - output: Output filename. The extension `.png` will be added automatically. - """ - import matplotlib - - matplotlib.use("Agg") - from matplotlib import pyplot as plt - - plt.xlabel("d (Å) (inverse-square scale)") - plt.ylabel("Intensity (counts)") - if ticks: - plt.xticks([1 / d**2 for d in ticks], [f"{d:g}" for d in ticks]) - - # Matplotlib v3.3.0 includes API change 'nonposy' → 'nonpositive' - # https://matplotlib.org/api/api_changes.html#log-symlog-scale-base-ticks-and-nonpos-specification - try: - plt.yscale("log", nonpositive="clip") - except ValueError: - plt.yscale("log", nonposy="clip") - - plt.plot(d_star_sq, intensity, "b.", label=None) - plt.plot( - d_star_sq, scaled_debye_waller(d_star_sq, *fit), "r-", label="Debye-Waller fit" - ) - if max_d: - plt.fill_betweenx( - plt.ylim(), - 1 / np.square(max_d), - color="k", - alpha=0.5, - zorder=2.1, - label="Excluded from fit", - ) - plt.legend(loc=0) - plt.savefig(output) - plt.close() - - -def suggest_minimum_exposure( - expts: ExperimentList, refls: flex.reflection_table, params: scope_extract -) -> None: - """ - Suggest an estimated minimum sufficient exposure to achieve a certain resolution. - - The estimate is based on a fit of a Debye-Waller factor under the assumption that a - single isotropic displacement parameter can be used to adequately describe the - decay of intensities with increasing sin(θ). - - An ASCII-art Wilson plot is printed, along with minimum exposure recommendations for - a number of different resolution targets. The Wilson plot, including the fitted - isotropic Debye-Waller factor, is saved as a PNG image. - - Args: - expts: Experiment list containing a single experiment, from which the crystal - symmetry will be extracted. - refls: Reflection table of observed reflections. - params: Parameters for calculation of minimum exposure estimate. - """ - # Ignore reflections without an index, since uctbx.unit_cell.d returns spurious - # d == -1 values, rather than None, for unindexed reflections. - refls.del_selected(refls["id"] == -1) - # Ignore all spots flagged as overloaded - refls.del_selected(refls.get_flags(refls.flags.overloaded).iselection()) - - # Work from profile-fitted intensities where possible but if the number of - # profile-fitted intensities is less than 75% of the number of summed - # intensities, use summed intensities instead. This is a very arbitrary heuristic. - sel_prf = refls.get_flags(refls.flags.integrated_prf).iselection() - sel_sum = refls.get_flags(refls.flags.integrated_sum).iselection() - if sel_prf.size() < 0.75 * sel_sum.size(): - refls = refls.select(sel_sum) - intensity = refls["intensity.sum.value"] - sigma = flex.sqrt(refls["intensity.sum.variance"]) - else: - refls = refls.select(sel_prf) - intensity = refls["intensity.prf.value"] - sigma = flex.sqrt(refls["intensity.prf.variance"]) - - # Apply French-Wilson scaling to ensure positive intensities. - miller_array = miller.array( - miller.set( - expts[0].crystal.get_crystal_symmetry(), - refls["miller_index"], - anomalous_flag=False, - ), - data=intensity, - sigmas=sigma, - ) - miller_array.set_observation_type_xray_intensity() - miller_array = miller_array.merge_equivalents().array() - cctbx_log = StringIO() # Prevent idiosyncratic CCTBX logging from polluting stdout. - miller_array = miller_array.french_wilson(log=cctbx_log).as_intensity_array() - logger.debug(cctbx_log.getvalue()) - - d_star_sq = miller_array.d_star_sq().data() - intensity = miller_array.data() - sigma = miller_array.sigmas() - - # Parameters for the lower-bound exposure estimate: - min_i_over_sigma = params.minimum_exposure.min_i_over_sigma - wilson_fit_max_d = params.minimum_exposure.wilson_fit_max_d - desired_d = params.minimum_exposure.desired_d - # If no target resolution is given, use the following defaults: - if not params.minimum_exposure.desired_d: - desired_d = [ - 1, # Å - 0.84, # Å (IUCr publication requirement) - 0.6, # Å - 0.4, # Å - ] - desired_d.sort(reverse=True) - - # Perform the Wilson plot fit - fit = wilson_fit(d_star_sq, intensity, sigma, wilson_fit_max_d) - - # Get recommended exposure factors - # Use the fact that σ² = I for indexed data, so I/σ = √̅I - desired_d_star_sq = [1 / d**2 for d in desired_d] - target_i = min_i_over_sigma**2 - recommended_factor = [ - (target_i / scaled_debye_waller(target_d, *fit)) - for target_d in desired_d_star_sq - ] - - # Get the achievable resolution at the current exposure - desired_d += [np.sqrt(fit[0] / (2 * np.log(fit[1] / target_i)))] - recommended_factor += [1] - - # Draw the ASCII art Wilson plot - wilson_plot_ascii(miller_array, d_ticks) - - recommendations = zip(desired_d, recommended_factor) - recommendations = sorted(recommendations, key=lambda rec: rec[0], reverse=True) - - # Print a recommendation to the user. - info("\nFitted isotropic displacement parameter, B = %.3g Ų", fit[0]) - for target, recommendation in recommendations: - if recommendation < 1: - debug( - "\nIt is likely that you can achieve a resolution of %g Å using a " - "lower exposure (lower transmission and/or shorter exposure time).", - target, - ) - elif recommendation > 1: - debug( - "\nIt is likely that you need a higher exposure (higher transmission " - "and/or longer exposure time to achieve a resolution of %g Å.", - target, - ) - debug( - "The estimated minimal sufficient exposure (flux × exposure time) to " - "achievea resolution of %.2g Å is %.3g times the exposure used for this " - "data collection.", - target, - recommendation, - ) - - summary = "\nRecommendations summarised:\n" - summary += tabulate( - recommendations, - ["Resolution (Å)", "Suggested\nexposure factor"], - floatfmt=(".3g", ".3g"), - tablefmt="rst", - ) - summary += ( - "\nExposure is flux × exposure time." - "\nYou can achieve your desired exposure factor by modifying " - "transmission and/or exposure time." - ) - info(summary) - - # Draw the Wilson plot image and save to file - wilson_plot_image( - d_star_sq, - intensity, - fit, - max_d=params.minimum_exposure.wilson_fit_max_d, - ticks=d_ticks, - output=params.output.wilson_plot, - ) - - -def run( - phil: scope = phil_scope, - args: Optional[List[str]] = None, - set_up_logging: bool = False, -) -> None: - """ - Parse command-line arguments, run the script. - - Uses the DIALS option parser to extract an experiment list, reflection table and - parameters, then passes them to :func:`suggest_minimum_exposure`. - Optionally, sets up the logger. - - Args: - phil: PHIL scope for option parser. - args: Arguments to parse. If None, :data:`sys.argv[1:]` will be used. - set_up_logging: Choose whether to configure :module:`screen19` logging. - """ - usage = "%(prog)s [options] integrated.expt integrated.refl" - - parser = OptionParser( - usage=usage, - phil=phil, - read_experiments=True, - read_reflections=True, - check_format=False, - epilog=help_message, - ) - - params, options = parser.parse_args(args=args) - - if set_up_logging: - # Configure the logging - log.config(params.verbosity, params.output.log) - - if not (params.input.experiments and params.input.reflections): - version_information = ( - f"screen19.minimum_exposure v{__version__} using {dials_version()} " - f"({time.strftime('%Y-%m-%d %H:%M:%S')})" - ) - - print(help_message) - print(version_information) - return - - if len(params.input.experiments) > 1: - warn( - "You provided more than one experiment list (%s). Only the " - "first, %s, will be used.", - ", ".join([expt.filename for expt in params.input.experiments]), - params.input.experiments[0].filename, - ) - if len(params.input.reflections) > 1: - warn( - "You provided more than one reflection table (%s). Only the " - "first, %s, will be used.", - ", ".join([refls.filename for refls in params.input.reflections]), - params.input.reflections[0].filename, - ) - - expts = params.input.experiments[0].data - refls = params.input.reflections[0].data - - if len(expts) > 1: - warn( - "The experiment list you provided, %s, contains more than one " - "experiment object (perhaps multiple indexing solutions). Only " - "the first will be used, all others will be ignored.", - params.input.experiments[0].filename, - ) - - suggest_minimum_exposure(expts, refls, params) - - -def main() -> None: - """Dispatcher for command-line call.""" - run(set_up_logging=True) diff --git a/screen19/screen.py b/screen19/screen.py deleted file mode 100644 index 0ae71a6..0000000 --- a/screen19/screen.py +++ /dev/null @@ -1,1192 +0,0 @@ -""" -Process screening data obtained at Diamond Light Source Beamline I19. - -This program presents the user with recommendations for adjustments to beam -flux, based on a single-sweep screening data collection. It presents an -upper- and lower-bound estimate of suitable flux. - • The upper-bound estimate is based on a comparison of a histogram of - measured pixel intensities with the trusted intensity range of the detector. - The user is warned when the measured pixel intensities indicate that the - detector would have a significant number of overloaded or untrustworthy - pixels. - • The lower-bound estimate is based on a linear fit of isotropic disorder - parameter, B, to a Wilson plot of reflection intensities. From this, - an estimate is made of the minimum exposure (flux × exposure time) required - to achieve a target I/σ ratio (by default, target I/σ = 2) at one or more values - of desired resolution, d, (by default, desired d = 1 Å, 0.84 Å, 0.6 Å & 0.4 Å). - -Target I/σ and target d (in Ångström) can be set using the parameters -'min_i_over_sigma' and 'desired_d'. One can set multiple values of the latter. - -By default the disorder parameter fit is conducted on the -integrated data. This ought to provide a reasonably true fit, but requires -an integration step, which can take some time. You can achieve a quicker, -dirtier answer by fitting to the indexed data (i.e. only the stronger -spots), using 'minimum_exposure.data=indexed'. - -Examples: - - screen19 imported_experiments.json - - screen19 *.cbf - - screen19 /path/to/data/ - - screen19 /path/to/data/image0001.cbf:1:100 - - screen19 min_i_over_sigma=2 desired_d=0.84 - - screen19 minimum_exposure.data=indexed - -""" - -import json -import logging -import math -import os -import re -import sys -import time -import timeit -from glob import glob -from pickle import PickleError -from typing import Dict, List, Optional, Sequence, Tuple - -import procrunner - -import iotbx.phil -from libtbx import Auto -from libtbx.introspection import number_of_processors -from libtbx.phil import scope - -import dials.command_line.integrate -import dials.util.version -from dials.algorithms.indexing import DialsIndexError -from dials.algorithms.indexing.bravais_settings import ( - refined_settings_from_refined_triclinic, -) -from dials.algorithms.shoebox import MaskCode -from dials.array_family import flex -from dials.command_line.dials_import import MetaDataUpdater -from dials.command_line.index import index -from dials.command_line.refine import run_dials_refine -from dials.command_line.refine_bravais_settings import ( - bravais_lattice_to_space_group_table, - eliminate_sys_absent, - map_to_primitive, -) -from dials.util import Sorry, log, version -from dials.util.ascii_art import spot_counts_per_image_plot -from dials.util.options import OptionParser -from dxtbx.model import ExperimentList -from dxtbx.model.experiment_list import ( - BeamComparison, - DetectorComparison, - ExperimentListFactory, - GoniometerComparison, -) - -import screen19 -from screen19.minimum_exposure import suggest_minimum_exposure - -Templates = List[Tuple[str, Tuple[int, int]]] - -phil_scope = iotbx.phil.parse( - """ - verbosity = 0 - .type = int(value_min=0) - .caption = 'Verbosity level of log output' - .help = "Possible values:\n" - "\t• 0: Info log output to stdout/logfile\n" - "\t• 1: Info & debug log output to stdout/logfile" - - output - .caption = 'Options to control the output files' - { - log = 'screen19.log' - .type = str - .caption = "The log filename" - } - nproc = Auto - .type = int - .caption = 'Number of processors to use' - .help = "The chosen value will apply to all the DIALS utilities with a " - "multi-processing option. If 'False' or 'Auto', all available " - "processors will be used." - - minimum_exposure - .caption = 'Options for screen19.minimum_exposure' - { - include scope screen19.minimum_exposure.phil_scope - data = indexed *integrated - .type = choice - .caption = 'Choice of data for the displacement parameter fit' - .help = 'For the lower-bound exposure estimate, choose whether to use ' - 'indexed (quicker) or integrated (better) data in fitting ' - 'the isotropic displacement parameter.' - } - - maximum_flux - .caption = 'Options for avoiding detector paralysation' - { - trusted_range_correction = 0.25 - .type = float(value_min=0, value_max=1) - .caption = 'Factor by which to multiply the maximum trusted flux.' - .help = "The detector manufacturer's photon count correction, to correct " - "for pixel paralysation, is often found to be unreliable at photon " - "counts in the upper part of the nominal trusted range. In such " - "cases, this factor can be used to adjust the upper limit of the " - "trusted range. Pilatus detectors, for example, have been found " - "not to give reliable correction for photon counts greater than " - "0.25 × the manufacturer's trusted range. It is therefore " - "sensible to present the user with a correspondingly reduced upper-" - "limit flux recommendation." - } - - dials_import - .caption = 'Options for dials.import' - { - include scope dials.command_line.dials_import.phil_scope - - input - { - include scope dials.util.options.tolerance_phil_scope - - experiments = None - .help = "The experiment list file path" - .type = str - .multiple = True - .optional = True - } - } - - dials_find_spots - .caption = 'Options for dials.find_spots' - { - include scope dials.command_line.find_spots.phil_scope - } - - dials_index - .caption = 'Options for dials.index' - { - include scope dials.command_line.index.phil_scope - } - - dials_refine - .caption = 'Options for dials.refine' - { - include scope dials.command_line.refine.phil_scope - } - - dials_refine_bravais - .caption = 'Options for dials.refine_bravais_settings' - { - include scope dials.command_line.refine_bravais_settings.phil_scope - } - - dials_create_profile - .caption = 'Options for dials.create_profile_model' - { - include scope dials.command_line.create_profile_model.phil_scope - } - - dials_integrate - .caption = 'Options for dials.integrate' - { - include scope dials.command_line.integrate.phil_scope - } - - dials_report - .caption = 'Options for dials.report' - { - include scope dials.command_line.report.phil_scope - } - """, - process_includes=True, -) - -logger = logging.getLogger("dials.screen19") -debug, info, warning = logger.debug, logger.info, logger.warning - - -def _run_integration( - phil_scope: scope, experiments_file: str, reflections_file: str -) -> Tuple[ExperimentList, flex.reflection_table]: - """Run integration programmatically, compatible with multiple DIALS versions. - - Args: - phil_scope: The dials.integrate phil scope - experiments_file: Path to the experiment list file - reflections_file: Path to the reflection table file - """ - - if hasattr(dials.command_line.integrate, "run_integration"): - # DIALS 3.1+ interface - expts, refls, _ = dials.command_line.integrate.run_integration( - phil_scope.extract(), - ExperimentList.from_file(experiments_file), - flex.reflection_table.from_file(reflections_file), - ) - elif hasattr(dials.command_line.integrate, "Script"): - # Pre-3.1-style programmatic interface - expts, refls = dials.command_line.integrate.Script(phil=phil_scope).run( - [experiments_file, reflections_file] - ) - else: - raise RuntimeError( - "Could not find dials.integrate programmatic interface 'run_integration' or 'Script'" - ) - - return expts, refls - - -def overloads_histogram( - d_spacings: Sequence[float], - ticks: Optional[Sequence[float]] = None, - output: Optional[str] = "overloads", -) -> None: - """ - Generate a histogram of reflection d-spacings as an image, default is .png. - - Args: - d_spacings: d-spacings of the reflections. - ticks (optional): d-values for the tick positions on the 1/d axis. - output (optional): Output filename root, to which the extension `.png` will - be appended. Default is `overloads`. - """ - import matplotlib - - matplotlib.use("Agg") - from matplotlib import pyplot as plt - - plt.xlabel("d (Å) (inverse scale)") - plt.ylabel("Number of overloaded reflections") - if ticks: - plt.xticks([1 / d for d in ticks], [f"{d:g}" for d in ticks]) - - # Matplotlib v3.3.0 includes API change 'nonposy' → 'nonpositive' - # https://matplotlib.org/api/api_changes.html#log-symlog-scale-base-ticks-and-nonpos-specification - try: - plt.yscale("log", nonpositive="clip") - except ValueError: - plt.yscale("log", nonposy="clip") - - plt.hist(d_spacings, min(100, d_spacings.size())) - plt.savefig(output) - plt.close() - - -class Screen19: - """Encapsulates the screening script.""" - - def __init__(self): - # Throughout the pipeline, retain the state of the processing. - self.expts = ExperimentList([]) - self.refls = flex.reflection_table() - # Get some default parameters. These must be extracted from the 'fetched' - # PHIL scope, rather than the 'definition' phil scope returned by - # iotbx.phil.parse. Confused? Blame PHIL. - self.params = phil_scope.fetch(iotbx.phil.parse("")).extract() - - def _quick_import(self, files: List[str]) -> bool: - """ - Generate xia2-style templates from file names and attempt a quick import. - - From each given filename, generate a filename template by substituting a hash - character (#) for each numeral in the last contiguous group of numerals - before the file extension. For example, the filename `example_01_0001.cbf` - becomes `example_01_####.cbf`. - - Contiguous image ranges are recorded by associating the start and end image - number of the range with the relevant filename template. - - dials.import is then run with options to extrapolate header information from - the first image file, thereby running more quickly than reading each image - header individually. - - Args: - files: List of image filenames. - - Returns: - Boolean flag indicating whether the quick import has succeeded. - """ - if len(files) == 1: - # No point in quick-importing a single file - return False - debug("Attempting quick import...") - files.sort() - templates: Dict[str, List[Optional[List[int]]]] = {} - for f in files: - template, image = screen19.make_template(f) - if template not in templates: - image_range = [image, image] if image else [] - templates.update({template: [image_range]}) - elif image == templates[template][-1][-1] + 1: - templates[template][-1][-1] = image - elif image == templates[template][-1][-1]: - # We have a duplicate input file name. Do nothing. - pass - else: - templates[template].append([image, image]) - # Return tuple of template and image range for each unique image range - templates: Templates = [ - (t, tuple(r)) for t, ranges in templates.items() for r in ranges - ] - return self._quick_import_templates(templates) - - def _quick_import_templates(self, templates: Templates) -> bool: - """ - Take image file templates and frame number ranges and try to run dials.import. - - dials.import is run with options to extrapolate header information from - the first image file, thereby running more quickly than reading each image - header individually. - - Args: - templates: A list of tuples, each tuple containing a xia2-style filename - template and the start and end image numbers of the associated - sweep. - - Returns: - Boolean flag indicating whether the quick import has succeeded. - """ - debug("Quick import template summary:\n\t%s", templates) - if len(templates) > 1: - debug("Cannot currently run quick import on multiple templates.") - return False - - try: - scan_range: Tuple[int, int] = templates[0][1] - if not scan_range: - raise IndexError - except IndexError: - debug("Cannot run quick import: could not determine image naming template.") - return False - - info("Running quick import.") - self.params.dials_import.input.template = [templates[0][0]] - self.params.dials_import.geometry.scan.image_range = scan_range - self.params.dials_import.geometry.scan.extrapolate_scan = True - self._run_dials_import() - - return True - - def _import(self, files: List[str]) -> None: - """ - Try to run a quick call of dials.import. Failing that, run a slow call. - - Try initially to construct file name templates contiguous groups of files. - Failing that, pass a full list of the files to the importer (slower). - - Args: - files: List of image filenames. - """ - info("\nImporting data...") - if len(files) == 1: - if os.path.isdir(files[0]): - debug( - "You specified a directory. Importing all CBF files in " - "that directory." - ) - # TODO Support HDF5. - files = [ - os.path.join(files[0], f) - for f in os.listdir(files[0]) - if f.endswith(".cbf") - or f.endswith(".cbf.gz") - or f.endswith(".cbf.bz2") - ] - elif len(files[0].split(":")) == 3: - debug( - "You specified an image range in the xia2 format. " - "Importing all specified files." - ) - template, start, end = files[0].split(":") - template = screen19.make_template(template)[0] - start, end = int(start), int(end) - if not self._quick_import_templates([(template, (start, end))]): - warning("Could not import specified image range.") - sys.exit(1) - info("Quick import successful.") - return - elif files[0].endswith(".expt"): - debug( - "You specified an existing experiment list file. " - "No import necessary." - ) - try: - self.expts = ExperimentList.from_file(files[0]) - except (OSError, PickleError, ValueError): - pass - else: - self.params.dials_import.output.experiments = files[0] - if self.expts: - return - - if not files: - warning("No images found matching input.") - sys.exit(1) - - # Can the files be quick-imported? - if self._quick_import(files): - info("Quick import successful.") - return - - self.params.dials_import.input.experiments = files - self._run_dials_import() - - def _run_dials_import(self): - """ - Perform a minimal version of dials.import to get an experiment list. - - Use some filleted bits of dials.import and dials.util.options.Importer. - """ - # Get some key data format arguments. - try: - format_kwargs = { - "dynamic_shadowing": self.params.dials_import.format.dynamic_shadowing, - "multi_panel": self.params.dials_import.format.multi_panel, - } - except AttributeError: - format_kwargs = {} - - # If filenames contain wildcards, expand - args = [] - for arg in self.params.dials_import.input.experiments: - if "*" in arg: - args.extend(glob(arg)) - else: - args.append(arg) - - if args: - # Are compare{beam,detector,goniometer} and scan_tolerance necessary? - # They are cargo-culted from the DIALS option parser. - tol_params = self.params.dials_import.input.tolerance - compare_beam = BeamComparison( - wavelength_tolerance=tol_params.beam.wavelength, - direction_tolerance=tol_params.beam.direction, - polarization_normal_tolerance=tol_params.beam.polarization_normal, - polarization_fraction_tolerance=tol_params.beam.polarization_fraction, - ) - compare_detector = DetectorComparison( - fast_axis_tolerance=tol_params.detector.fast_axis, - slow_axis_tolerance=tol_params.detector.slow_axis, - origin_tolerance=tol_params.detector.origin, - ) - compare_goniometer = GoniometerComparison( - rotation_axis_tolerance=tol_params.goniometer.rotation_axis, - fixed_rotation_tolerance=tol_params.goniometer.fixed_rotation, - setting_rotation_tolerance=tol_params.goniometer.setting_rotation, - ) - scan_tolerance = tol_params.scan.oscillation - - # Import an experiment list from image data. - try: - experiments = ExperimentListFactory.from_filenames( - args, - compare_beam=compare_beam, - compare_detector=compare_detector, - compare_goniometer=compare_goniometer, - scan_tolerance=scan_tolerance, - format_kwargs=format_kwargs, - ) - except OSError as e: - warning("%s '%s'", e.strerror, e.filename) - sys.exit(1) - - # Record the imported experiments for use elsewhere. - # Quit if there aren't any. - self.expts.extend(experiments) - if not self.expts: - warning("No images found.") - sys.exit(1) - - else: - # Use the template importer. - if len(self.params.dials_import.input.template) > 0: - experiments = ExperimentList.from_templates( - self.params.dials_import.input.template, - format_kwargs=format_kwargs, - ) - # Record the imported experiments for use elsewhere. - # Quit if there aren't any. - self.expts.extend(experiments) - if not self.expts: - warning( - "No images found matching template %s", - self.params.dials_import.input.template[0], - ) - sys.exit(1) - - # Setup the metadata updater - metadata_updater = MetaDataUpdater(self.params.dials_import) - - # Extract the experiments and loop through - self.expts = metadata_updater(self.expts.imagesets()) - - def _count_processors(self, nproc: Optional[int] = None) -> None: - """ - Determine the number of processors and save it as an instance variable. - - The user may specify the number of processors to use. If no value is - given, the number of available processors is returned. - - Args: - nproc: Number of processors. - """ - if nproc and nproc is not Auto: - self.nproc = nproc - return - - # if environmental variable NSLOTS is set to a number then use that - try: - self.nproc = int(os.environ.get("NSLOTS")) - return - except (ValueError, TypeError): - pass - - self.nproc = number_of_processors(return_value_if_unknown=-1) - - if self.nproc <= 0: - warning( - "Could not determine number of available processors. Error code %d", - self.nproc, - ) - sys.exit(1) - - def _count_images(self) -> int: - """ - Attempt to determine the number of diffraction images. - - The number of diffraction images is determined from the imported_experiments - JSON file. - - Returns: - Number of images. - """ - # FIXME: This exception handling should be redundant. Empty experiment - # lists should get caught at the import stage. Is this so? - try: - return self.expts[0].imageset.size() - except IndexError: - warning("Could not determine number of images in dataset.") - sys.exit(1) - - def _check_intensities(self, mosaicity_correction: bool = True) -> None: - """ - Run xia2.overload and plot a histogram of pixel intensities. - - If `mosaicity_correction` is true, the pixel intensities are approximately - adjusted to take account of a systematic defect in the detector count rate - correction. See https://github.com/xia2/screen19/wiki#mosaicity-correction - - Args: - mosaicity_correction: default is `True`. - """ - info("\nTesting pixel intensities...") - command = ["xia2.overload", f"nproc={self.nproc}", "indexed.expt"] - - debug("running %s", command) - start = timeit.default_timer() - result = procrunner.run(command, print_stdout=False) - - debug("result = %s", screen19.prettyprint_procrunner(result)) - info("Successfully completed (%.1f sec)", timeit.default_timer() - start) - - if result.returncode: - warning("Failed with exit code %d", result.returncode) - sys.exit(1) - - with open("overload.json") as fh: - overload_data = json.load(fh) - - info("Pixel intensity distribution:") - count_sum = 0 - hist = {} - if "bins" in overload_data: - for b in range(overload_data["bin_count"]): - if overload_data["bins"][b] > 0: - hist[b] = overload_data["bins"][b] - count_sum += b * overload_data["bins"][b] - else: - hist = {int(k): v for k, v in overload_data["counts"].items() if int(k) > 0} - count_sum = sum([k * v for k, v in hist.items()]) - - average_to_peak = 1 - if mosaicity_correction: - # Adjust for the detector count rate correction - if self._sigma_m: - delta_z = self._oscillation / self._sigma_m / math.sqrt(2) - average_to_peak = ( - math.sqrt(math.pi) * delta_z * math.erf(delta_z) - + math.exp(-(delta_z**2)) - - 1 - ) / delta_z**2 - info("Average-to-peak intensity ratio: %f", average_to_peak) - - scale = 100 * overload_data["scale_factor"] / average_to_peak - info("Determined scale factor for intensities as %f", scale) - - debug( - "intensity histogram: { %s }", - ", ".join([f"{k:d}:{hist[k]:d}" for k in sorted(hist)]), - ) - max_count = max(hist.keys()) - hist_max = max_count * scale - hist_granularity, hist_format = 1, ".0f" - if hist_max < 50: - hist_granularity, hist_format = 2, ".1f" - if hist_max < 15: - hist_granularity, hist_format = 10, ".1f" - rescaled_hist = {} - for x in hist.keys(): - rescaled = round(x * scale * hist_granularity) - if rescaled > 0: - rescaled_hist[rescaled] = hist[x] + rescaled_hist.get(rescaled, 0) - hist = rescaled_hist - debug( - "rescaled histogram: { %s }", - ", ".join( - [ - f"{k / hist_granularity:{hist_format}}:{hist[k]:d}" - for k in sorted(hist) - ] - ), - ) - - screen19.plot_intensities(hist, 1 / hist_granularity) - - linear_response_limit = 100 * self.params.maximum_flux.trusted_range_correction - marginal_limit = max(70, linear_response_limit) - - text = ( - f"Strongest pixel ({max_count:d} counts) reaches {hist_max:.1f}% of the detector count rate limit", - ) - if hist_max > 100: - warning("Warning: %s!", text) - else: - info(text) - if ( - "overload_limit" in overload_data - and max_count >= overload_data["overload_limit"] - ): - warning( - "Warning: THE DATA CONTAIN REGULAR OVERLOADS!\n" - " The photon incidence rate is outside the specified " - "limits of the detector.\n" - " The built-in detector count rate correction cannot " - "adjust for this.\n" - " You should aim for count rates below " - f"{self.params.maximum_flux.trusted_range_correction:.0%} of the " - "detector limit." - ) - elif hist_max > marginal_limit: - warning( - "Warning: The photon incidence rate is well outside the " - "linear response region of the detector " - f"(<{self.params.maximum_flux.trusted_range_correction:.0%}).\n" - " The built-in detector count rate correction may not be " - "able to adjust for this." - ) - elif hist_max > linear_response_limit: - info( - "The photon incidence rate is outside the linear response " - "region of the detector " - f"(<{self.params.maximum_flux.trusted_range_correction:.0%}).\n" - " The built-in detector count rate correction may be able " - "to adjust for this." - ) - if not mosaicity_correction: - warning( - "Warning: Not enough data for proper profile estimation." - " The spot intensities are not corrected for mosaicity.\n" - " The true photon incidence rate will be higher than the " - "given estimate." - ) - - info("Total sum of counts in dataset: %d", count_sum) - - def _find_spots(self, args: Optional[List[str]] = None) -> None: - """ - Call `dials.find_spots` on the imported experiment list. - - Args: - args: List of any additional PHIL parameters to be used by dials.import. - """ - info("\nFinding spots...") - - dials_start = timeit.default_timer() - - # Use some choice fillets from dials.find_spots - # Ignore `args`, use `self.params` - - # Loop through all the imagesets and find the strong spots - - self.refls = flex.reflection_table.from_observations( - self.expts, self.params.dials_find_spots - ) - - # Add n_signal column - before deleting shoeboxes - - good = MaskCode.Foreground | MaskCode.Valid - self.refls["n_signal"] = self.refls["shoebox"].count_mask_values(good) - - # Delete the shoeboxes - if not self.params.dials_find_spots.output.shoeboxes: - del self.refls["shoebox"] - - info( - 60 * "-" + "\n%s\n" + 60 * "-" + "\nSuccessfully completed (%.1f sec)", - spot_counts_per_image_plot(self.refls), - timeit.default_timer() - dials_start, - ) - - def _index(self) -> bool: - """ - Call `dials.index` on the output of spot finding. - - Returns: - Boolean value indicating whether indexing was successful. - """ - dials_start = timeit.default_timer() - - # Prepare max_cell constraint strategies. - max_cell = self.params.dials_index.indexing.max_cell - # By default, try unconstrained max_cell followed by max_cell=20. - # If the user has already specified a max_cell < 20, do not relax to 20Å. - cell_constraints = [([], max_cell)] - if not max_cell or max_cell is Auto or max_cell > 20: - cell_constraints += [(["max_cell constraint"], 20)] - - # Prepare indexing methods, preferring the real_space_grid_search if a - # known unit cell has been specified, otherwise using 3D FFT, then 1D FFT. - methods = ( - [(["real space grid search"], "real_space_grid_search")] - if self.params.dials_index.indexing.known_symmetry.unit_cell - else [] - ) - methods += [(["3D FFT"], "fft3d"), (["1D FFT"], "fft1d")] - - # Cycle through the indexing methods for each of the max_cell constraint - # strategies until an indexing solution is found. - for i, (max_cell_msg, max_cell) in enumerate(cell_constraints): - # Set the max_cell constraint strategy. - self.params.dials_index.indexing.max_cell = max_cell - for j, (method_msg, method) in enumerate(methods): - # Set the indexing method. - self.params.dials_index.indexing.method = method - # Log a handy message to the user. - msg = ( - "Retrying with " + " and ".join(method_msg + max_cell_msg) - if i + j - else "Indexing" - ) - info("\n%s...", msg) - try: - # If indexing is successful, break out of the inner loop. - self.expts, self.refls = index( - self.expts, [self.refls], self.params.dials_index - ) - break - except (DialsIndexError, ValueError) as e: - # If indexing is unsuccessful, try again with the next - # strategy. - warning("Failed: %s", str(e)) - continue - else: - # When all the indexing methods are unsuccessful, move onto - # the next max_cell constraint strategy and try again. - continue - # We should only get here if successfully indexed. Break out of the loop - break - else: - # Indexing completely unsuccessful. - return False - - sg_type = self.expts[0].crystal.get_crystal_symmetry().space_group().type() - symb = sg_type.universal_hermann_mauguin_symbol() - unit_cell = self.expts[0].crystal.get_unit_cell() - - self.refls.as_file(self.params.dials_index.output.reflections) - self.expts.as_file(self.params.dials_index.output.experiments) - self.refls.as_file(self.params.dials_index.output.reflections) - info( - "Found primitive solution: %s %s using %s reflections\n" - "Indexed experiments and reflections saved as %s, %s\n" - "Successfully completed (%.1f sec)", - symb, - unit_cell, - self.refls["id"].count(0), - self.params.dials_index.output.experiments, - self.params.dials_index.output.reflections, - timeit.default_timer() - dials_start, - ) - - # Report the indexing successful. - return True - - def _wilson_calculation(self) -> None: - """ - Run `screen19.minimum_exposure` on an experiment list and reflection table. - - For best results, the reflections and experiment list should contain the - results of integration or scaling. If only strong spots are used, the Wilson - plot fit may be poor. - """ - dials_start = timeit.default_timer() - info("\nEstimating lower exposure bound...") - - suggest_minimum_exposure(self.expts, self.refls, self.params.minimum_exposure) - - info("Successfully completed (%.1f sec)", timeit.default_timer() - dials_start) - - def _refine(self) -> None: - """ - Run `dials.refine` on the results of indexing. - """ - dials_start = timeit.default_timer() - info("\nRefining...") - - try: - self.expts, self.refls, _, _ = run_dials_refine( - self.expts, self.refls, self.params.dials_refine - ) - except Sorry as e: - warning("dials.refine failed: %d\nGiving up.\n", e) - sys.exit(1) - - info("Successfully refined (%.1f sec)", timeit.default_timer() - dials_start) - - def _create_profile_model(self) -> bool: - """ - Run `dials.create_profile_model` on indexed reflections. - - The indexed experiment list will be overwritten with a copy that includes - the profile model but is otherwise identical. - - Returns: - Boolean value indicating whether it was possible to determine a profile - model from the data. - """ - info("\nCreating profile model...") - command = [ - "dials.create_profile_model", - self.params.dials_index.output.experiments, - self.params.dials_index.output.reflections, - f"output = {self.params.dials_index.output.experiments}", - ] - - start = timeit.default_timer() - result = procrunner.run(command, print_stdout=False) - - debug("result = %s", screen19.prettyprint_procrunner(result)) - self._sigma_m = None - if result.returncode == 0: - db = ExperimentList.from_file(self.params.dials_index.output.experiments)[0] - self._oscillation = db.imageset.get_scan().get_oscillation()[1] - self._sigma_m = db.profile.sigma_m() - info( - "%d images, %s° oscillation, σ_m=%.3f°", - db.imageset.get_scan().get_num_images(), - str(self._oscillation), - self._sigma_m, - ) - info("Successfully completed (%.1f sec)", timeit.default_timer() - start) - return True - warning("Failed with exit code %d", result.returncode) - return False - - def _integrate(self) -> None: - """Run `dials.integrate` to integrate reflection intensities.""" - dials_start = timeit.default_timer() - info("\nIntegrating...") - - # Don't waste time recreating the profile model - self.params.dials_integrate.create_profile_model = False - # Get the dials.integrate PHIL scope, populated with parsed input parameters - integrate_scope = phil_scope.get("dials_integrate").objects[0] - integrate_scope.name = "" - integrate_scope = integrate_scope.format(self.params.dials_integrate) - - try: - integrated_experiments, integrated_reflections = _run_integration( - integrate_scope, - self.params.dials_index.output.experiments, - self.params.dials_index.output.reflections, - ) - # Save the output to files - integrated_reflections.as_file( - self.params.dials_integrate.output.reflections - ) - integrated_experiments.as_file( - self.params.dials_integrate.output.experiments - ) - # ... and also store the output internally - self.expts, self.refls = integrated_experiments, integrated_reflections - info( - "Successfully completed (%.1f sec)", - timeit.default_timer() - dials_start, - ) - except SystemExit as e: - if e.code: - warning("dials.integrate failed with exit code %d\nGiving up.", e.code) - sys.exit(1) - - # This is a hacky check but should work for as long as DIALS 2.0 is supported. - if version.dials_version() < "DIALS 2.1": - - def _refine_bravais( - self, experiments: ExperimentList, reflections: flex.reflection_table - ) -> None: - """ - Run `dials.refine_bravais_settings` on an experiments and reflections. - - Args: - experiments: An experiment list.. - reflections: The corresponding reflection table. - """ - info("\nRefining Bravais settings...") - command = ["dials.refine_bravais_settings", experiments, reflections] - - start = timeit.default_timer() - result = procrunner.run(command, print_stdout=False) - - debug("result = %s", screen19.prettyprint_procrunner(result)) - if result.returncode == 0: - m = re.search( - r"[-+]{3,}\n[^\n]*\n[-+|]{3,}\n(.*\n)*[-+]{3,}", - result.stdout.decode("utf-8"), - ) - if m: - info(m.group(0)) - else: - info( - "Could not interpret dials.refine_bravais_settings output, " - "please check dials.refine_bravais_settings.log" - ) - info( - "Successfully completed (%.1f sec)", timeit.default_timer() - start - ) - else: - warning("Failed with exit code %d", result.returncode) - sys.exit(1) - - else: - - def _refine_bravais(self) -> None: - """Run `dials.refine_bravais_settings` to determine the space group.""" - dials_start = timeit.default_timer() - info("\nRefining Bravais settings...") - - self.refls = eliminate_sys_absent(self.expts, self.refls) - map_to_primitive(self.expts, self.refls) - - # We can only refine against reflections with non-zero variance in - # observed centroid position. - x_valid, y_valid, z_valid = ( - part > 0 for part in self.refls["xyzobs.mm.variance"].parts() - ) - nonzero_variance = x_valid | y_valid | z_valid - - try: - refined_settings = refined_settings_from_refined_triclinic( - self.expts, - self.refls.select(nonzero_variance), - self.params.dials_refine_bravais, - ) - except RuntimeError as e: - warning("dials.refine_bravais_settings failed.\nGiving up.") - sys.exit(e) - - possible_bravais_settings = { - solution["bravais"] for solution in refined_settings - } - bravais_lattice_to_space_group_table(possible_bravais_settings) - try: - # Old version of dials with as_str() method - logger.info(refined_settings.as_str()) - except AttributeError: - # Newer versions of dials (>= 2.2.2) has proper __str__ method - logger.info(refined_settings) - - info( - "Successfully completed (%.1f sec)", - timeit.default_timer() - dials_start, - ) - - def _report( - self, experiments: ExperimentList, reflections: flex.reflection_table - ) -> None: - """ - Run `dials.report` on an experiment list and reflection table. - - Args: - experiments: An experiment list. - reflections: The corresponding reflection table. - """ - info("\nCreating report...") - command = ["dials.report", experiments, reflections] - - start = timeit.default_timer() - result = procrunner.run(command, print_stdout=False) - - debug("result = %s", screen19.prettyprint_procrunner(result)) - if result.returncode == 0: - info("Successfully completed (%.1f sec)", timeit.default_timer() - start) - # if sys.stdout.isatty(): - # info("Trying to start browser") - # try: - # import subprocess - # d = dict(os.environ) - # d["LD_LIBRARY_PATH"] = "" - # subprocess.Popen(["xdg-open", "dials-report.html"], env=d) - # except Exception as e: - # debug("Could not open browser\n%s", str(e)) - else: - warning("Failed with exit code %d", result.returncode) - sys.exit(1) - - def run( - self, - args: Optional[List[str]] = None, - phil: scope = phil_scope, - set_up_logging: bool = False, - ) -> None: - """ - TODO: Docstring. - - Args: - args: - phil: - set_up_logging: - - Returns: - - """ - usage = "%(prog)s [options] image_directory | image_files.cbf | imported.expt" - - parser = OptionParser( - usage=usage, epilog=__doc__, phil=phil, check_format=False - ) - - self.params, options, unhandled = parser.parse_args( - args=args, show_diff_phil=True, return_unhandled=True, quick_parse=True - ) - - version_information = ( - f"screen19 v{screen19.__version__} " - f"using {dials.util.version.dials_version()} " - f"({time.strftime('%Y-%m-%d %H:%M:%S')})" - ) - - start = timeit.default_timer() - - if len(unhandled) == 0: - print(__doc__) - print(version_information) - return - - if set_up_logging: - # Configure the logging - log.config(verbosity=self.params.verbosity, logfile=self.params.output.log) - # Unless verbose output has been requested, suppress generation of - # debug and info log records from any child DIALS command, retaining - # those from screen19 itself. - if not self.params.verbosity: - logging.getLogger("dials").setLevel(logging.WARNING) - logging.getLogger("dials.screen19").setLevel(logging.INFO) - - info(version_information) - debug("Run with:\n%s\n%s", " ".join(unhandled), parser.diff_phil.as_str()) - - self._count_processors(nproc=self.params.nproc) - debug("Using %s processors", self.nproc) - # Set multiprocessing settings for spot-finding, indexing and - # integration to match the top-level specified number of processors - self.params.dials_find_spots.spotfinder.mp.nproc = self.nproc - self.params.dials_index.indexing.nproc = self.nproc - # Setting self.params.dials_refine.refinement.mp.nproc is not helpful - self.params.dials_integrate.integration.mp.nproc = self.nproc - - # Set the input and output parameters for the DIALS components - # TODO: Compare to diff_phil and start from later in the pipeline if - # appropriate - self._import(unhandled) - imported_name = self.params.dials_import.output.experiments - - self._find_spots() - - if not self._index(): - info("\nRetrying for stronger spots only...") - strong_refls = self.refls - self.params.dials_find_spots.spotfinder.threshold.dispersion.sigma_strong = ( - 15 - ) - self._find_spots() - - if not self._index(): - warning("Giving up.") - self.expts.as_file(imported_name) - strong_refls.as_file("strong.refl") - self.refls.as_file("stronger.refl") - info( - "Could not find an indexing solution. You may want to " - "have a look at the reciprocal space by running:\n\n" - " dials.reciprocal_lattice_viewer %s %s\n\n" - "or, to only include stronger spots:\n\n" - " dials.reciprocal_lattice_viewer %s %s\n", - imported_name, - "strong.refl", - imported_name, - "stronger.refl", - ) - sys.exit(1) - - if not self._create_profile_model(): - info("\nRefining model to attempt to increase number of valid spots...") - self._refine() - if not self._create_profile_model(): - warning("Giving up.") - info( - "The identified indexing solution may not be correct. " - "You may want to have a look at the reciprocal space by " - "running:\n\n" - " dials.reciprocal_lattice_viewer indexed.expt indexed.refl\n" - ) - sys.exit(1) - - self._check_intensities() - - if self.params.minimum_exposure.data == "integrated": - self._integrate() - - self._wilson_calculation() - - experiments = self.params.dials_integrate.output.experiments - reflections = self.params.dials_integrate.output.reflections - else: - self._wilson_calculation() - - experiments = self.params.dials_create_profile.output - reflections = self.params.dials_index.output.reflections - - # This is a hacky check but should work for as long as DIALS 2.0 is supported. - if version.dials_version() < "DIALS 2.1": - self._refine_bravais(experiments, reflections) - else: - self._refine_bravais() - - self._report(experiments, reflections) - - runtime = timeit.default_timer() - start - debug( - "Finished at %s, total runtime: %.1f", - time.strftime("%Y-%m-%d %H:%M:%S"), - runtime, - ) - info("screen19 successfully completed (%.1f sec).", runtime) - - -def main() -> None: - """Dispatcher for command-line call.""" - Screen19().run(set_up_logging=True) diff --git a/tests/conftest.py b/tests/conftest.py deleted file mode 100644 index a4bda99..0000000 --- a/tests/conftest.py +++ /dev/null @@ -1,9 +0,0 @@ -import pytest - -try: - import dials_data as _ # noqa: F401 -except ImportError: - - @pytest.fixture(scope="session") - def dials_data(): - pytest.skip("Test requires python package dials_data") diff --git a/tests/test_screen19.py b/tests/test_screen19.py deleted file mode 100644 index f19c874..0000000 --- a/tests/test_screen19.py +++ /dev/null @@ -1,89 +0,0 @@ -import pytest - -from screen19 import minimum_exposure -from screen19.screen import Screen19 - -# A list of tuples of example sys.argv[1:] cases and associated image count. -import_checks = [ - ([""], 900), - (["/"], 900), - (["x3_1_00*.cbf.gz"], 99), - (["x3_1_####.cbf.gz:1:99"], 99), - (["x3_1_00##.cbf.gz:1:99"], 99), - (["x3_1_0001.cbf.gz:1:99"], 99), - ( - [ - "x3_1_0001.cbf.gz", - "x3_1_0002.cbf.gz", - "x3_1_0003.cbf.gz", - "x3_1_0004.cbf.gz", - "x3_1_0005.cbf.gz", - ], - 5, - ), -] - - -def test_screen19_command_line_help_does_not_crash(): - Screen19().run([]) - - -def test_minimum_exposure_help_does_not_crash(): - minimum_exposure.run(args=[]) - - -@pytest.mark.parametrize("import_checks", import_checks) -def test_screen19_inputs(dials_data, tmpdir, import_checks): - """Test various valid input argument styles""" - data_files, image_count = import_checks - data = [ - dials_data("small_molecule_example").join(filename).strpath - for filename in data_files - ] - - foo = Screen19() - # The tmpdir should only be necessary for DIALS v1 — no output expected for DIALS v2 - foo._import(data) - - # Check that the import has resulted in the creation of a single experiment. - assert len(foo.expts) == 1 - # Check that the associated imageset has the expected number of images. - assert foo.expts[0].imageset.size() == image_count - - -def test_screen19(dials_data, tmpdir): - """An integration test. Check the full functionality of screen19.""" - data_dir = dials_data("x4wide").join("X4_wide_M1S4_2_####.cbf:1:30").strpath - - # Test screen19 first. - with tmpdir.as_cwd(): - Screen19().run([data_dir], set_up_logging=True) - - logfile = tmpdir.join("screen19.log").read() - - assert "screen19 successfully completed" in logfile - assert "photon incidence rate is outside the linear response region" in logfile - - # Then check screen19.minimum_exposure. - with tmpdir.as_cwd(): - minimum_exposure.run( - args=["integrated.expt", "integrated.refl"], set_up_logging=True - ) - - m_e_logfile = tmpdir.join("screen19.minimum_exposure.log").read() - - assert ( - "You can achieve your desired exposure factor by modifying transmission " - "and/or exposure time." in m_e_logfile - ) - - -def test_screen19_single_frame(dials_data, tmpdir): - image = dials_data("x4wide").join("X4_wide_M1S4_2_0001.cbf").strpath - - with tmpdir.as_cwd(): - Screen19().run([image], set_up_logging=True) - - logfile = tmpdir.join("screen19.log").read() - - assert "screen19 successfully completed" in logfile From f78a02d7615a7a108341600d0924cddcf7e18fac Mon Sep 17 00:00:00 2001 From: Ben Williams Date: Wed, 26 Oct 2022 14:06:58 +0100 Subject: [PATCH 02/33] Raze the old entry points --- setup.py | 16 ++-------------- 1 file changed, 2 insertions(+), 14 deletions(-) diff --git a/setup.py b/setup.py index d37e415..27dec15 100644 --- a/setup.py +++ b/setup.py @@ -20,20 +20,8 @@ description="Screening program for small-molecule single-crystal X-ray diffraction " "data", entry_points={ - "console_scripts": [ - "i19.screen = screen19.screen:main", - "screen19 = screen19.screen:main", - "i19.stability_fft = screen19.stability_fft:main", - "i19.minimum_exposure = screen19.minimum_exposure:main", - "screen19.minimum_exposure = screen19.minimum_exposure:main", - ], - "libtbx.dispatcher.script": [ - "i19.screen = i19.screen", - "screen19 = screen19", - "i19.stability_fft = i19.stability_fft", - "i19.minimum_exposure = i19.minimum_exposure", - "screen19.minimum_exposure = screen19.minimum_exposure", - ], + "console_scripts": [], + "libtbx.dispatcher.script": [], "libtbx.precommit": ["screen19 = screen19"], }, install_requires=['typing;python_version<"3.5"', "procrunner"], From 8038d10a86287be3d4c1801e5194d2b6d4616bae Mon Sep 17 00:00:00 2001 From: Ben Williams Date: Wed, 26 Oct 2022 14:07:15 +0100 Subject: [PATCH 03/33] Set up Azure Pipelines CI (maybe?) --- .azure-pipelines/azure-pipelines.yml | 130 ++++++++++++++++++ .../check-working-directory-is-clean | 0 .azure-pipelines/ci.yml | 41 ++++++ .azure-pipelines/flake8-validation.py | 43 ++++++ {.travis => .azure-pipelines}/run-flake8 | 4 +- {.travis => .azure-pipelines}/setup-base | 6 - .azure-pipelines/syntax-validation.py | 32 +++++ 7 files changed, 248 insertions(+), 8 deletions(-) create mode 100644 .azure-pipelines/azure-pipelines.yml rename {.travis => .azure-pipelines}/check-working-directory-is-clean (100%) create mode 100644 .azure-pipelines/ci.yml create mode 100644 .azure-pipelines/flake8-validation.py rename {.travis => .azure-pipelines}/run-flake8 (63%) rename {.travis => .azure-pipelines}/setup-base (80%) create mode 100644 .azure-pipelines/syntax-validation.py diff --git a/.azure-pipelines/azure-pipelines.yml b/.azure-pipelines/azure-pipelines.yml new file mode 100644 index 0000000..f8e79eb --- /dev/null +++ b/.azure-pipelines/azure-pipelines.yml @@ -0,0 +1,130 @@ +trigger: + branches: + include: + - '*' + tags: + include: + - '*' + +stages: +- stage: static + displayName: Static Analysis + jobs: + - job: checks + displayName: static code analysis + pool: + vmImage: ubuntu-22.04 + steps: + # Use Python >=3.8 for syntax validation + - task: UsePythonVersion@0 + displayName: Set up python + inputs: + versionSpec: 3.8 + + # Run syntax validation on a shallow clone + - bash: | + libtbx.python "$(Pipeline.Workspace)/.azure-pipelines/syntax-validation.py" + displayName: Syntax validation + + # Run flake8 validation on a shallow clone + - bash: | + libtbx.pip install flake8 + libtbx.python "$(Pipeline.Workspace)/.azure-pipelines/flake8-validation.py" + displayName: Flake8 validation + +- stage: build + displayName: Build + dependsOn: + jobs: + - job: build + displayName: build package + pool: + vmImage: ubuntu-22.04 + steps: + - task: UsePythonVersion@0 + displayName: Set up python + inputs: + versionSpec: 3.10 + + - bash: | + $(Pipeline.Workspace)/.azure-pipelines/setup-base + displayName: Set up DIALS + + - bash: | + set -ex + libtbx.python setup.py sdist bdist_wheel + mkdir -p dist/pypi + shopt -s extglob + mv -v dist/!(pypi) dist/pypi + git archive HEAD | gzip > dist/repo-source.tar.gz + ls -laR dist + displayName: Build python package + + - task: PublishBuildArtifacts@1 + displayName: Store artifact + inputs: + pathToPublish: dist/ + artifactName: package + + - bash: libtbx.python "$(Pipeline.Workspace)/setup.py" checkdocs + displayName: Check package description + +- stage: tests + displayName: Run unit tests + dependsOn: + - static + - build + jobs: + - job: linux + pool: + vmImage: ubuntu-22.04 + strategy: + matrix: + python38: + PYTHON_VERSION: 3.8 + python39: + PYTHON_VERSION: 3.9 + python310: + PYTHON_VERSION: 3.10 + steps: + - template: ci.yml + +- stage: deploy + displayName: Publish release + dependsOn: + - tests + condition: and(succeeded(), startsWith(variables['Build.SourceBranch'], 'refs/tags/')) + jobs: + - job: pypi + displayName: Publish pypi release + pool: + vmImage: ubuntu-22.04 + steps: + - checkout: none + + - task: UsePythonVersion@0 + displayName: Set up python + inputs: + versionSpec: 3.8 + + - task: DownloadBuildArtifacts@0 + displayName: Get pre-built package + inputs: + buildType: current + downloadType: single + artifactName: package + downloadPath: $(System.ArtifactsDirectory) + + - script: | + libtbx.pip install -U pip + libtbx.pip install twine + displayName: Install twine + + - task: TwineAuthenticate@1 + displayName: Set up credentials + inputs: + pythonUploadServiceConnection: pypi-tristan + + - bash: | + libtbx.python -m twine upload -r pypi-tristan --config-file $(PYPIRC_PATH) $(System.ArtifactsDirectory)/package/pypi/*.tar.gz $(System.ArtifactsDirectory)/package/pypi/*.whl + displayName: Publish package diff --git a/.travis/check-working-directory-is-clean b/.azure-pipelines/check-working-directory-is-clean similarity index 100% rename from .travis/check-working-directory-is-clean rename to .azure-pipelines/check-working-directory-is-clean diff --git a/.azure-pipelines/ci.yml b/.azure-pipelines/ci.yml new file mode 100644 index 0000000..456b2dd --- /dev/null +++ b/.azure-pipelines/ci.yml @@ -0,0 +1,41 @@ +steps: +- checkout: none + +- task: UsePythonVersion@0 + inputs: + versionSpec: $(PYTHON_VERSION) + displayName: Use Python $(PYTHON_VERSION) + +- task: DownloadBuildArtifacts@0 + displayName: Get pre-built package + inputs: + buildType: current + downloadType: single + artifactName: package + downloadPath: $(System.ArtifactsDirectory) + +- task: ExtractFiles@1 + displayName: Checkout sources + inputs: + archiveFilePatterns: $(System.ArtifactsDirectory)/package/repo-source.tar.gz + destinationFolder: $(Pipeline.Workspace)/src + +- script: | + set -eux + pip install -r "$(Pipeline.Workspace)/src/requirements_dev.txt" + pip install "$(Pipeline.Workspace)/src" + displayName: Install package + +- script: | + PYTHONDEVMODE=1 pytest -ra --cov=tristan --cov-report=xml --cov-branch + displayName: Run tests + workingDirectory: $(Pipeline.Workspace)/src + +- bash: | + curl -Os https://uploader.codecov.io/latest/linux/codecov + chmod +x codecov + ./codecov -B $(Build.SourceBranchName) -C $(Build.SourceVersion) -t $(CODECOV_TOKEN) -n "Python $(PYTHON_VERSION) $(Agent.OS)" + displayName: Publish coverage stats + continueOnError: true + workingDirectory: $(Pipeline.Workspace)/src + timeoutInMinutes: 2 diff --git a/.azure-pipelines/flake8-validation.py b/.azure-pipelines/flake8-validation.py new file mode 100644 index 0000000..da2ff7e --- /dev/null +++ b/.azure-pipelines/flake8-validation.py @@ -0,0 +1,43 @@ +import os +import subprocess + +# Flake8 validation +failures = 0 +try: + flake8 = subprocess.run( + [ + "flake8", + "--exit-zero", + "--max-line-length=88", + "--select=E401,E711,E712,E713,E714,E721,E722,E901,F401,F402,F403,F405,F631,F632,F633,F811,F812,F821,F822,F841,F901,W191,W291,W292,W293,W602,W603,W604,W605,W606", + ], + capture_output=True, + check=True, + encoding="latin-1", + timeout=300, + ) +except (subprocess.CalledProcessError, subprocess.TimeoutExpired) as e: + print( + "##vso[task.logissue type=error;]flake8 validation failed with", + str(e.__class__.__name__), + ) + print(e.stdout) + print(e.stderr) + print("##vso[task.complete result=Failed;]flake8 validation failed") + exit() +else: + for line in flake8.stdout.split("\n"): + if ":" not in line: + continue + filename, lineno, column, error = line.split(":", maxsplit=3) + errcode, error = error.strip().split(" ", maxsplit=1) + filename = os.path.normpath(filename) + failures += 1 + print( + f"##vso[task.logissue type=error;sourcepath={filename};" + f"linenumber={lineno};columnnumber={column};code={errcode};]" + error + ) + +if failures: + print(f"##vso[task.logissue type=warning]Found {failures} flake8 violation(s)") + print(f"##vso[task.complete result=Failed;]Found {failures} flake8 violation(s)") diff --git a/.travis/run-flake8 b/.azure-pipelines/run-flake8 similarity index 63% rename from .travis/run-flake8 rename to .azure-pipelines/run-flake8 index 195ec6e..0113f67 100755 --- a/.travis/run-flake8 +++ b/.azure-pipelines/run-flake8 @@ -6,7 +6,7 @@ CHANGED_FILES=$(git diff --name-only --diff-filter=AM $TRAVIS_BRANCH...HEAD -- " # List is empty if this is not a pull request run or if no python files have changed [ -z "$CHANGED_FILES" ] && echo There are no relevant changes. Skipping test. && exit 0 -pip install -q flake8 +libtbx.pip install -q flake8 echo echo Running flake8 on changed files -flake8 --select=E401,E711,E712,E713,E714,E721,E722,E901,F401,F402,F403,F405,F631,F632,F633,F811,F812,F821,F822,F841,F901,W191,W291,W292,W293,W602,W603,W604,W605,W606 $CHANGED_FILES +libtbx.flake8 --select=E401,E711,E712,E713,E714,E721,E722,E901,F401,F402,F403,F405,F631,F632,F633,F811,F812,F821,F822,F841,F901,W191,W291,W292,W293,W602,W603,W604,W605,W606 $CHANGED_FILES diff --git a/.travis/setup-base b/.azure-pipelines/setup-base similarity index 80% rename from .travis/setup-base rename to .azure-pipelines/setup-base index 0441be2..c239870 100755 --- a/.travis/setup-base +++ b/.azure-pipelines/setup-base @@ -2,12 +2,6 @@ EXPECTED_CACHE_REVISION=20191115-python3 -if [ "$TRAVIS_EVENT_TYPE" == "cron" ]; then - echo -e "\e[31;1mCron job build detected. Invalidating cache.\e[0m" - rm -f $HOME/build_dials/.cache_valid - rm -f $HOME/build_dials/.build_complete -fi - if [ -f $HOME/build_dials/.cache_valid ] && [ "$EXPECTED_CACHE_REVISION" == "$(cat $HOME/build_dials/.cache_valid)" ]; then echo -e "\e[1mCache probably valid\e[0m" else diff --git a/.azure-pipelines/syntax-validation.py b/.azure-pipelines/syntax-validation.py new file mode 100644 index 0000000..2d74948 --- /dev/null +++ b/.azure-pipelines/syntax-validation.py @@ -0,0 +1,32 @@ +from __future__ import annotations + +import ast +import os +import sys + +print("Python", sys.version, "\n") + +failures = 0 + +for base, _, files in os.walk("."): + for f in files: + if not f.endswith(".py"): + continue + filename = os.path.normpath(os.path.join(base, f)) + try: + with open(filename) as fh: + ast.parse(fh.read()) + except SyntaxError as se: + failures += 1 + print( + f"##vso[task.logissue type=error;sourcepath={filename};" + f"linenumber={se.lineno};columnnumber={se.offset};]" + f"SyntaxError: {se.msg}" + ) + print(" " + se.text + " " * se.offset + "^") + print(f"SyntaxError: {se.msg} in {filename} line {se.lineno}") + print() + +if failures: + print(f"##vso[task.logissue type=warning]Found {failures} syntax error(s)") + print(f"##vso[task.complete result=Failed;]Found {failures} syntax error(s)") From cd951e84faab22dbe180661b7d43dfbb1944a96f Mon Sep 17 00:00:00 2001 From: Ben Williams Date: Wed, 26 Oct 2022 14:53:11 +0100 Subject: [PATCH 04/33] Set up basic scripts and tests --- src/screen/__init__.py | 0 src/screen/minimum_exposure.py | 493 +++++++++++++++++++++++++++++++++ src/screen/screen.py | 26 ++ tests/conftest.py | 9 + tests/test_screen.py | 52 ++++ 5 files changed, 580 insertions(+) create mode 100644 src/screen/__init__.py create mode 100644 src/screen/minimum_exposure.py create mode 100644 src/screen/screen.py create mode 100644 tests/conftest.py create mode 100644 tests/test_screen.py diff --git a/src/screen/__init__.py b/src/screen/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/screen/minimum_exposure.py b/src/screen/minimum_exposure.py new file mode 100644 index 0000000..406d794 --- /dev/null +++ b/src/screen/minimum_exposure.py @@ -0,0 +1,493 @@ +""" +Perform straight-line Wilson plot fit. Draw the Wilson plot. + +Reflection d-spacings are determined from the crystal symmetry (from +indexing) and the Miller indices of the indexed reflections. The +atomic displacement parameter is assumed isotropic. Its value is +determined from a fit to the reflection data: + I = A * exp(-B / (2 * d^2)), +where I is the intensity and the scale factor, A, and isotropic +displacement parameter, B, are the fitted parameters. + +An I/σ condition for 'good' diffraction statistics is set by the +instance variable min_i_over_sigma, and the user's desired +resolution is set by the instance variable desired_d. A crude +error model is assumed, whereby σ² = I, and so the I/σ condition +translates trivially to a threshold I. + +The value of the fitted intensity function at the desired +resolution is compared with the threshold I. The ratio of these +values is used to determine a recommended exposure (flux × exposure time) +for the full data collection. + +The Wilson plot of I as a function of d is drawn as the file +'wilson_plot.png'. The plot can optionally be saved in other formats. + +Examples: + + screen19.minimum_exposure integrated.expt integrated.refl + + screen19.minimum_exposure indexed.expt indexed.refl + + screen19.minimum_exposure min_i_over_sigma=2 desired_d=0.84 wilson_fit_max_d=4 \ + integrated.expt integrated.refl + +""" + +import logging +import time +from io import StringIO +from typing import Iterable, List, Optional, Sequence, Union + +import numpy as np +from scipy.optimize import curve_fit +from tabulate import tabulate + +import iotbx.phil +from cctbx import miller +from libtbx.phil import scope, scope_extract + +from dials.array_family import flex +from dials.util import log +from dials.util.options import OptionParser +from dials.util.version import dials_version +from dxtbx.model import ExperimentList + +from screen19 import __version__, d_ticks, plot_intensities, terminal_size + +# Custom types +FloatSequence = Sequence[float] +Fit = Union[np.ndarray, Iterable, int, float] + +help_message = __doc__ + + +phil_scope = iotbx.phil.parse( + """ + verbosity = 0 + .type = int(value_min=0) + .caption = 'Verbosity level of log output' + .help = "Possible values:\n" + "\t• 0: Info log output to stdout/logfile\n" + "\t• 1: Info & debug log output to stdout/logfile" + minimum_exposure + .caption = 'Parameters for the calculation of the lower exposure bound' + { + desired_d = None + .multiple = True + .type = float + .caption = u'Desired resolution limit, in Ångströms, of diffraction data' + .help = 'This is the resolution target for the lower-bound exposure ' \ + 'recommendation.' + min_i_over_sigma = 2 + .type = float + .caption = u'Target I/σ value for lower-bound exposure recommendation' + .help = 'The lower-bound exposure recommendation provides an estimate of ' \ + u'the exposure (flux × exposure time) required to ensure that the' \ + 'majority of expected reflections at the desired resolution limit' \ + u'have I/σ greater than or equal to this value.' + wilson_fit_max_d = 4 # Å + .type = float + .caption = u'Maximum d-value (in Ångströms) for displacement parameter fit' + .help = 'Reflections with lower resolution than this value will be ' \ + 'ignored for the purposes of the Wilson plot.' + } + output + .caption = 'Parameters to control the output' + { + log = 'screen19.minimum_exposure.log' + .type = str + .caption = 'Location for the info log' + debug_log = 'screen19.minimum_exposure.debug.log' + .type = str + .caption = 'Location for the debug log' + wilson_plot = 'wilson_plot' + .type = str + .caption = 'Filename for the Wilson plot image' + .help = "By default, the extension '.png' is appended. If you include " \ + "a different extension, either '.pdf', '.ps', '.eps' or '.svg', " \ + "a file of that format will be created instead." + } + """, + process_includes=True, +) + +logger_name = "dials.screen19.minimum_exposure" +logger = logging.getLogger(logger_name) +debug, info, warn = logger.debug, logger.info, logger.warning + + +def scaled_debye_waller(x: float, b: float, a: float) -> float: + """ + Calculate a scaled isotropic Debye-Waller factor. + + By assuming a single isotropic disorder parameter, :param:`b`, this factor + approximates the decay of diffracted X-ray intensity increasing resolution + (decreasing d, increasing sin(θ)). + + Args: + x: Equivalent to 1/d². + b: Isotropic displacement parameter. + a: A scale factor. + + Returns: + Estimated value of scaled isotropic Debye-Waller factor. + """ + return a * np.exp(-b / 2 * x) + + +def wilson_fit( + d_star_sq: FloatSequence, + intensity: FloatSequence, + sigma: FloatSequence, + wilson_fit_max_d: float, +) -> Fit: + """ + Fit a simple Debye-Waller factor, assume isotropic disorder parameter. + + Reflections with d ≥ :param:`wilson_fit_max_d` are ignored. + + Args: + d_star_sq: 1/d² (equivalently d*²), sequence of values for the observed + reflections (units of Å⁻² assumed). + intensity: Sequence of reflection intensities. + sigma: Sequence of uncertainties in reflection intensity. + wilson_fit_max_d: The minimum resolution for reflections against which to + fit. + + Returns: + - The fitted isotropic displacement parameter (units of Ų assumed); + - The fitted scale factor. + + """ + # Eliminate reflections with d > wilson_fit_max_d from the fit + sel = d_star_sq > 1 / wilson_fit_max_d**2 + + # Perform a weighted Wilson plot fit to the reflection intensities + fit, cov = curve_fit( + scaled_debye_waller, + d_star_sq.select(sel), + intensity.select(sel), + sigma=sigma.select(sel), + bounds=(0, np.inf), + ) + + return fit + + +def wilson_plot_ascii( + miller_array: miller.array, d_ticks: Optional[Sequence] = None +) -> None: + """ + Print an ASCII-art Wilson plot of reflection intensities. + + Equivalent reflections will be merged according to the crystal symmetry. + + Args: + miller_array: An array of integrated intensities, bundled with appropriate + crystal symmetry and unit cell info. + d_ticks: d location of ticks on 1/d² axis. + """ + # Draw the Wilson plot, using existing functionality in cctbx.miller + columns, rows = terminal_size() + n_bins = min(columns, miller_array.data().size()) + miller_array.setup_binner_counting_sorted(n_bins=n_bins, reflections_per_bin=1) + wilson = miller_array.wilson_plot(use_binning=True) + # Get the relevant plot data from the miller_array: + binned_intensity = [x if x else 0 for x in wilson.data[1:-1]] + bins = dict(zip(wilson.binner.bin_centers(1), binned_intensity)) + if d_ticks: + tick_positions = ", ".join([f'"{d:g}" {1 / d ** 2}' for d in d_ticks]) + tick_positions = tick_positions.join(["(", ")"]) + else: + tick_positions = "" + # Draw the plot: + plot_intensities( + bins, + 1, + title="'Wilson plot'", + xlabel="'d (Angstrom) (inverse-square scale)'", + ylabel="'I (counts)'", + xticks=tick_positions, + style="with lines", + ) + + +def wilson_plot_image( + d_star_sq: FloatSequence, + intensity: FloatSequence, + fit: Fit, + max_d: Optional[float] = None, + ticks: Optional[FloatSequence] = None, + output: str = "wilson_plot", +) -> None: + """ + Generate the Wilson plot as a PNG image. + + :param:`max_d` allows greying out of the reflections not included in the + isotropic Debye-Waller fit. + + Args: + d_star_sq: 1/d² values of reflections. + intensity: Intensities of reflections. + fit: Fitted parameters (tuple of fitted isotropic displacement parameter and + fitted scale factor). + max_d: The minimum resolution for reflections used in the Debye-Waller fit. + ticks: d location of ticks on 1/d² axis. + output: Output filename. The extension `.png` will be added automatically. + """ + import matplotlib + + matplotlib.use("Agg") + from matplotlib import pyplot as plt + + plt.xlabel("d (Å) (inverse-square scale)") + plt.ylabel("Intensity (counts)") + if ticks: + plt.xticks([1 / d**2 for d in ticks], [f"{d:g}" for d in ticks]) + + # Matplotlib v3.3.0 includes API change 'nonposy' → 'nonpositive' + # https://matplotlib.org/api/api_changes.html#log-symlog-scale-base-ticks-and-nonpos-specification + try: + plt.yscale("log", nonpositive="clip") + except ValueError: + plt.yscale("log", nonposy="clip") + + plt.plot(d_star_sq, intensity, "b.", label=None) + plt.plot( + d_star_sq, scaled_debye_waller(d_star_sq, *fit), "r-", label="Debye-Waller fit" + ) + if max_d: + plt.fill_betweenx( + plt.ylim(), + 1 / np.square(max_d), + color="k", + alpha=0.5, + zorder=2.1, + label="Excluded from fit", + ) + plt.legend(loc=0) + plt.savefig(output) + plt.close() + + +def suggest_minimum_exposure( + expts: ExperimentList, refls: flex.reflection_table, params: scope_extract +) -> None: + """ + Suggest an estimated minimum sufficient exposure to achieve a certain resolution. + + The estimate is based on a fit of a Debye-Waller factor under the assumption that a + single isotropic displacement parameter can be used to adequately describe the + decay of intensities with increasing sin(θ). + + An ASCII-art Wilson plot is printed, along with minimum exposure recommendations for + a number of different resolution targets. The Wilson plot, including the fitted + isotropic Debye-Waller factor, is saved as a PNG image. + + Args: + expts: Experiment list containing a single experiment, from which the crystal + symmetry will be extracted. + refls: Reflection table of observed reflections. + params: Parameters for calculation of minimum exposure estimate. + """ + # Ignore reflections without an index, since uctbx.unit_cell.d returns spurious + # d == -1 values, rather than None, for unindexed reflections. + refls.del_selected(refls["id"] == -1) + # Ignore all spots flagged as overloaded + refls.del_selected(refls.get_flags(refls.flags.overloaded).iselection()) + + # Work from profile-fitted intensities where possible but if the number of + # profile-fitted intensities is less than 75% of the number of summed + # intensities, use summed intensities instead. This is a very arbitrary heuristic. + sel_prf = refls.get_flags(refls.flags.integrated_prf).iselection() + sel_sum = refls.get_flags(refls.flags.integrated_sum).iselection() + if sel_prf.size() < 0.75 * sel_sum.size(): + refls = refls.select(sel_sum) + intensity = refls["intensity.sum.value"] + sigma = flex.sqrt(refls["intensity.sum.variance"]) + else: + refls = refls.select(sel_prf) + intensity = refls["intensity.prf.value"] + sigma = flex.sqrt(refls["intensity.prf.variance"]) + + # Apply French-Wilson scaling to ensure positive intensities. + miller_array = miller.array( + miller.set( + expts[0].crystal.get_crystal_symmetry(), + refls["miller_index"], + anomalous_flag=False, + ), + data=intensity, + sigmas=sigma, + ) + miller_array.set_observation_type_xray_intensity() + miller_array = miller_array.merge_equivalents().array() + cctbx_log = StringIO() # Prevent idiosyncratic CCTBX logging from polluting stdout. + miller_array = miller_array.french_wilson(log=cctbx_log).as_intensity_array() + logger.debug(cctbx_log.getvalue()) + + d_star_sq = miller_array.d_star_sq().data() + intensity = miller_array.data() + sigma = miller_array.sigmas() + + # Parameters for the lower-bound exposure estimate: + min_i_over_sigma = params.minimum_exposure.min_i_over_sigma + wilson_fit_max_d = params.minimum_exposure.wilson_fit_max_d + desired_d = params.minimum_exposure.desired_d + # If no target resolution is given, use the following defaults: + if not params.minimum_exposure.desired_d: + desired_d = [ + 1, # Å + 0.84, # Å (IUCr publication requirement) + 0.6, # Å + 0.4, # Å + ] + desired_d.sort(reverse=True) + + # Perform the Wilson plot fit + fit = wilson_fit(d_star_sq, intensity, sigma, wilson_fit_max_d) + + # Get recommended exposure factors + # Use the fact that σ² = I for indexed data, so I/σ = √̅I + desired_d_star_sq = [1 / d**2 for d in desired_d] + target_i = min_i_over_sigma**2 + recommended_factor = [ + (target_i / scaled_debye_waller(target_d, *fit)) + for target_d in desired_d_star_sq + ] + + # Get the achievable resolution at the current exposure + desired_d += [np.sqrt(fit[0] / (2 * np.log(fit[1] / target_i)))] + recommended_factor += [1] + + # Draw the ASCII art Wilson plot + wilson_plot_ascii(miller_array, d_ticks) + + recommendations = zip(desired_d, recommended_factor) + recommendations = sorted(recommendations, key=lambda rec: rec[0], reverse=True) + + # Print a recommendation to the user. + info("\nFitted isotropic displacement parameter, B = %.3g Ų", fit[0]) + for target, recommendation in recommendations: + if recommendation < 1: + debug( + "\nIt is likely that you can achieve a resolution of %g Å using a " + "lower exposure (lower transmission and/or shorter exposure time).", + target, + ) + elif recommendation > 1: + debug( + "\nIt is likely that you need a higher exposure (higher transmission " + "and/or longer exposure time to achieve a resolution of %g Å.", + target, + ) + debug( + "The estimated minimal sufficient exposure (flux × exposure time) to " + "achievea resolution of %.2g Å is %.3g times the exposure used for this " + "data collection.", + target, + recommendation, + ) + + summary = "\nRecommendations summarised:\n" + summary += tabulate( + recommendations, + ["Resolution (Å)", "Suggested\nexposure factor"], + floatfmt=(".3g", ".3g"), + tablefmt="rst", + ) + summary += ( + "\nExposure is flux × exposure time." + "\nYou can achieve your desired exposure factor by modifying " + "transmission and/or exposure time." + ) + info(summary) + + # Draw the Wilson plot image and save to file + wilson_plot_image( + d_star_sq, + intensity, + fit, + max_d=params.minimum_exposure.wilson_fit_max_d, + ticks=d_ticks, + output=params.output.wilson_plot, + ) + + +def run( + phil: scope = phil_scope, + args: Optional[List[str]] = None, + set_up_logging: bool = False, +) -> None: + """ + Parse command-line arguments, run the script. + + Uses the DIALS option parser to extract an experiment list, reflection table and + parameters, then passes them to :func:`suggest_minimum_exposure`. + Optionally, sets up the logger. + + Args: + phil: PHIL scope for option parser. + args: Arguments to parse. If None, :data:`sys.argv[1:]` will be used. + set_up_logging: Choose whether to configure :module:`screen19` logging. + """ + usage = "%(prog)s [options] integrated.expt integrated.refl" + + parser = OptionParser( + usage=usage, + phil=phil, + read_experiments=True, + read_reflections=True, + check_format=False, + epilog=help_message, + ) + + params, options = parser.parse_args(args=args) + + if set_up_logging: + # Configure the logging + log.config(params.verbosity, params.output.log) + + if not (params.input.experiments and params.input.reflections): + version_information = ( + f"screen19.minimum_exposure v{__version__} using {dials_version()} " + f"({time.strftime('%Y-%m-%d %H:%M:%S')})" + ) + + print(help_message) + print(version_information) + return + + if len(params.input.experiments) > 1: + warn( + "You provided more than one experiment list (%s). Only the " + "first, %s, will be used.", + ", ".join([expt.filename for expt in params.input.experiments]), + params.input.experiments[0].filename, + ) + if len(params.input.reflections) > 1: + warn( + "You provided more than one reflection table (%s). Only the " + "first, %s, will be used.", + ", ".join([refls.filename for refls in params.input.reflections]), + params.input.reflections[0].filename, + ) + + expts = params.input.experiments[0].data + refls = params.input.reflections[0].data + + if len(expts) > 1: + warn( + "The experiment list you provided, %s, contains more than one " + "experiment object (perhaps multiple indexing solutions). Only " + "the first will be used, all others will be ignored.", + params.input.experiments[0].filename, + ) + + suggest_minimum_exposure(expts, refls, params) + + +def main() -> None: + """Dispatcher for command-line call.""" + run(set_up_logging=True) diff --git a/src/screen/screen.py b/src/screen/screen.py new file mode 100644 index 0000000..537175f --- /dev/null +++ b/src/screen/screen.py @@ -0,0 +1,26 @@ +"""The main screening script.""" + +import argparse +import subprocess + +parser = argparse.ArgumentParser(description=__doc__) +parser.add_argument("images", help="Path to the images to process.") + + +def pipeline(images): + subprocess.run(["dials.import", images]) + subprocess.run(["dials.find_spots", "imported.expt"]) + subprocess.run(["dials.pixel_histogram", "strong.refl"]) + subprocess.run(["dials.index", "imported.expt", "strong.refl"]) + subprocess.run(["dials.refine", "indexed.expt", "indexed.refl"]) + subprocess.run(["dials.integrate", "refined.expt", "refined.refl"]) + subprocess.run(["screen19.minimum_exposure", "integrated.expt", "integrated.refl"]) + + +def main(args=None): + args = parser.parse_args(args) + pipeline(args.images) + + +if __name__ == "__main__": + main() diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 0000000..a4bda99 --- /dev/null +++ b/tests/conftest.py @@ -0,0 +1,9 @@ +import pytest + +try: + import dials_data as _ # noqa: F401 +except ImportError: + + @pytest.fixture(scope="session") + def dials_data(): + pytest.skip("Test requires python package dials_data") diff --git a/tests/test_screen.py b/tests/test_screen.py new file mode 100644 index 0000000..485ae7c --- /dev/null +++ b/tests/test_screen.py @@ -0,0 +1,52 @@ +import pytest + +import dxtbx + +from screen import minimum_exposure, screen + +# A list of tuples of example sys.argv[1:] cases and associated image count. +import_checks = [ + ([""], 900), + (["/"], 900), + (["x3_1_00*.cbf.gz"], 99), + (["x3_1_####.cbf.gz:1:99"], 99), + (["x3_1_00##.cbf.gz:1:99"], 99), + (["x3_1_0001.cbf.gz:1:99"], 99), + ( + [ + "x3_1_0001.cbf.gz", + "x3_1_0002.cbf.gz", + "x3_1_0003.cbf.gz", + "x3_1_0004.cbf.gz", + "x3_1_0005.cbf.gz", + ], + 5, + ), +] + + +def test_screen19_help_does_not_crash(): + screen.main(args=[]) + + +def test_minimum_exposure_help_does_not_crash(): + minimum_exposure.run(args=[]) + + +@pytest.mark.parametrize("import_checks", import_checks) +def test_screen19_inputs(dials_data, tmp_path, import_checks): + """Test various valid input argument styles""" + data_files, image_count = import_checks + data = [ + dials_data("small_molecule_example").join(filename).strpath + for filename in data_files + ] + + # Run screen19 + screen.main(data) + + # Check that the import has resulted in the creation of a single experiment. + expts = dxtbx.model.ExperimentList.from_files(tmp_path / "imported.expt") + assert len(expts) == 1 + # Check that the associated imageset has the expected number of images. + assert expts[0].imageset.size() == image_count From 6a7db0bf84f7e9ffa6f96c13fcd68640a658cf7d Mon Sep 17 00:00:00 2001 From: Ben Williams Date: Wed, 26 Oct 2022 15:12:50 +0100 Subject: [PATCH 05/33] Setuptools nonsense... --- setup.cfg | 7 +++---- setup.py | 16 +++++++++++----- 2 files changed, 14 insertions(+), 9 deletions(-) diff --git a/setup.cfg b/setup.cfg index 71e9849..c07885e 100644 --- a/setup.cfg +++ b/setup.cfg @@ -27,18 +27,17 @@ classifiers = License :: OSI Approved :: BSD License Natural Language :: English Programming Language :: Python :: 3 - Programming Language :: Python :: 3.6 - Programming Language :: Python :: 3.7 Programming Language :: Python :: 3.8 Programming Language :: Python :: 3.9 + Programming Language :: Python :: 3.10 Operating System :: OS Independent -project-urls = +project_urls = Download = https://github.com/xia2/screen19/releases GitHub = https://github.com/xia2/screen19 Bug-Tracker = https://github.com/xia2/screen19/issues [options] -python_requires = >=3.6 +python_requires = >=3.8 [tool:pytest] collect_ignore = ['setup.py'] diff --git a/setup.py b/setup.py index 27dec15..fd4eed6 100644 --- a/setup.py +++ b/setup.py @@ -19,19 +19,25 @@ author_email="scientificsoftware@diamond.ac.uk", description="Screening program for small-molecule single-crystal X-ray diffraction " "data", + package_dir={"": "src"}, entry_points={ - "console_scripts": [], - "libtbx.dispatcher.script": [], + "console_scripts": [ + "screen19 = screen.screen:main" + "screen19.minimum_exposure = screen.minimum_exposure:main" + ], + "libtbx.dispatcher.script": [ + "screen19.minimum_exposure = screen19.minimum_exposure" + ], "libtbx.precommit": ["screen19 = screen19"], }, - install_requires=['typing;python_version<"3.5"', "procrunner"], + install_requires=['python_version>="3.8"'], license="BSD license", long_description="\n\n".join([readme, changelog_header, changelog]), include_package_data=True, name="screen19", - packages=find_packages(), + packages=find_packages(where="src/screen"), test_suite="tests", - tests_require=["mock>=2.0", "pytest>=4.5"], + tests_require=["pytest"], version="0.213", zip_safe=False, ) From b86ba39ef0ef36f3cb1892499be855780b20d88a Mon Sep 17 00:00:00 2001 From: Ben Williams Date: Wed, 26 Oct 2022 15:42:11 +0100 Subject: [PATCH 06/33] More setuptools nonsense... --- MANIFEST.in | 7 ------- setup.cfg | 29 ++++++++++++++++++++++++++++- setup.py | 30 ++---------------------------- 3 files changed, 30 insertions(+), 36 deletions(-) delete mode 100644 MANIFEST.in diff --git a/MANIFEST.in b/MANIFEST.in deleted file mode 100644 index fc4ad3c..0000000 --- a/MANIFEST.in +++ /dev/null @@ -1,7 +0,0 @@ -include CHANGELOG.rst -include LICENSE -include README.rst - -recursive-include tests * -recursive-exclude * __pycache__ -recursive-exclude * *.py[co] diff --git a/setup.cfg b/setup.cfg index c07885e..c5f589b 100644 --- a/setup.cfg +++ b/setup.cfg @@ -22,7 +22,13 @@ exclude = docs test = pytest [metadata] -classifiers = +name = screen19 +author = Diamond Light Source +author_email = data_analysis@diamond.ac.uk +description = Screening program for small-molecule single-crystal X-ray diffraction data +licence = BSD license +licence_file = LICENSE +classifiers = Development Status :: 4 - Beta License :: OSI Approved :: BSD License Natural Language :: English @@ -37,7 +43,28 @@ project_urls = Bug-Tracker = https://github.com/xia2/screen19/issues [options] +package_dir = + =src +packages = find: python_requires = >=3.8 +include_package_data = True +install_requires = python_version>="3.8" +test_suite = tests +tests_require = pytest + +[options.packages.find] +where = src + +[options.entry_points] +console_scripts = + screen19 = screen.screen:main + screen19.minimum_exposure = screen.minimum_exposure:main +libtbx.dispatcher.script = screen19.minimum_exposure = screen19.minimum_exposure" +libtbx.precommit = screen19 = screen19 [tool:pytest] collect_ignore = ['setup.py'] + + + +# long_description="\n\n".join([readme, changelog_header, changelog]), diff --git a/setup.py b/setup.py index fd4eed6..a822c2e 100644 --- a/setup.py +++ b/setup.py @@ -1,6 +1,6 @@ #!/usr/bin/env python -from setuptools import find_packages, setup +from setuptools import setup with open("README.rst") as readme_file: readme = readme_file.read() @@ -14,30 +14,4 @@ with open("CHANGELOG.rst") as changelog_file: changelog = changelog_file.read() -setup( - author="Diamond Light Source", - author_email="scientificsoftware@diamond.ac.uk", - description="Screening program for small-molecule single-crystal X-ray diffraction " - "data", - package_dir={"": "src"}, - entry_points={ - "console_scripts": [ - "screen19 = screen.screen:main" - "screen19.minimum_exposure = screen.minimum_exposure:main" - ], - "libtbx.dispatcher.script": [ - "screen19.minimum_exposure = screen19.minimum_exposure" - ], - "libtbx.precommit": ["screen19 = screen19"], - }, - install_requires=['python_version>="3.8"'], - license="BSD license", - long_description="\n\n".join([readme, changelog_header, changelog]), - include_package_data=True, - name="screen19", - packages=find_packages(where="src/screen"), - test_suite="tests", - tests_require=["pytest"], - version="0.213", - zip_safe=False, -) +setup() From 20d8a90cf069808b2ff4c1fb983bd9c117a53eb7 Mon Sep 17 00:00:00 2001 From: Noemi Frisina Date: Wed, 26 Oct 2022 16:09:21 +0100 Subject: [PATCH 07/33] Even more setuptools stuff ... --- setup.cfg | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/setup.cfg b/setup.cfg index c5f589b..2a4af54 100644 --- a/setup.cfg +++ b/setup.cfg @@ -26,6 +26,7 @@ name = screen19 author = Diamond Light Source author_email = data_analysis@diamond.ac.uk description = Screening program for small-molecule single-crystal X-ray diffraction data +long_description = file: README.rst licence = BSD license licence_file = LICENSE classifiers = @@ -48,9 +49,10 @@ package_dir = packages = find: python_requires = >=3.8 include_package_data = True -install_requires = python_version>="3.8" +install_requires = python_version>='3.8' test_suite = tests tests_require = pytest +zip_safe = False [options.packages.find] where = src @@ -59,12 +61,10 @@ where = src console_scripts = screen19 = screen.screen:main screen19.minimum_exposure = screen.minimum_exposure:main -libtbx.dispatcher.script = screen19.minimum_exposure = screen19.minimum_exposure" -libtbx.precommit = screen19 = screen19 +libtbx.dispatcher.script = + screen19.minimum_exposure = screen19.minimum_exposure +libtbx.precommit = + screen19 = screen19 [tool:pytest] collect_ignore = ['setup.py'] - - - -# long_description="\n\n".join([readme, changelog_header, changelog]), From 9f2f58df0827079d57115ee338fb650afdd9f9f9 Mon Sep 17 00:00:00 2001 From: Noemi Frisina Date: Wed, 26 Oct 2022 16:43:52 +0100 Subject: [PATCH 08/33] Get correct versioning from setup.cfg --- setup.cfg | 5 +++-- src/screen/__init__.py | 6 ++++++ 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/setup.cfg b/setup.cfg index 2a4af54..5866489 100644 --- a/setup.cfg +++ b/setup.cfg @@ -7,11 +7,11 @@ serialize = {major}.{minor}.{patch} {major}.{minor} -[bumpversion:file:setup.py] +[bumpversion:file:setup.cfg] search = version="{current_version}", replace = version="{new_version}", -[bumpversion:file:screen19/__init__.py] +[bumpversion:file:screen/__init__.py] search = __version__ = "{current_version}" replace = __version__ = "{new_version}" @@ -23,6 +23,7 @@ test = pytest [metadata] name = screen19 +version = 0.213 author = Diamond Light Source author_email = data_analysis@diamond.ac.uk description = Screening program for small-molecule single-crystal X-ray diffraction data diff --git a/src/screen/__init__.py b/src/screen/__init__.py index e69de29..00a5ea0 100644 --- a/src/screen/__init__.py +++ b/src/screen/__init__.py @@ -0,0 +1,6 @@ +"""Screening program for small-molecule single-crystal X-ray diffraction data""" + +__author__ = "Diamond Light Source - Scientific Software" +__email__ = "data_analysis@diamond.ac.uk" +__version__ = "0.213" +__version_tuple__ = tuple(int(x) for x in __version__.split(".")) \ No newline at end of file From 08b612fa03ed8e7994b87c6260675ded3993ccba Mon Sep 17 00:00:00 2001 From: Noemi Frisina Date: Mon, 31 Oct 2022 14:32:48 +0000 Subject: [PATCH 09/33] update pre-commit to 3.8 --- .pre-commit-config.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 2b49cbd..a852f41 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -6,7 +6,7 @@ repos: rev: v2.38.2 hooks: - id: pyupgrade - args: ["--py36-plus"] + args: ["--py38-plus"] # Automatically sort imports with isort - repo: https://github.com/PyCQA/isort From 4e1d345303b29629ce65962421ab8f2236dd139d Mon Sep 17 00:00:00 2001 From: Noemi Frisina Date: Mon, 31 Oct 2022 14:48:44 +0000 Subject: [PATCH 10/33] Add missing bits for basic script to run. Exchange procrunner with subprocess --- src/screen/__init__.py | 157 ++++++++++++++++++++++++++++++++- src/screen/minimum_exposure.py | 19 ++-- src/screen/screen.py | 8 +- 3 files changed, 168 insertions(+), 16 deletions(-) diff --git a/src/screen/__init__.py b/src/screen/__init__.py index 00a5ea0..b19c3a1 100644 --- a/src/screen/__init__.py +++ b/src/screen/__init__.py @@ -1,6 +1,161 @@ """Screening program for small-molecule single-crystal X-ray diffraction data""" +from __future__ import annotations __author__ = "Diamond Light Source - Scientific Software" __email__ = "data_analysis@diamond.ac.uk" __version__ = "0.213" -__version_tuple__ = tuple(int(x) for x in __version__.split(".")) \ No newline at end of file +__version_tuple__ = tuple(int(x) for x in __version__.split(".")) + + +import logging +import re +import subprocess +import sys +import traceback + +logger = logging.getLogger("dials.screen19") +debug, info, warn = logger.debug, logger.info, logger.warning +# FIXME setting verbosity to 1 doesn't seem to print any debugging. + +# Set axis tick positions manually. Accounts for reciprocal(-square) d-scaling +d_ticks = [5, 3, 2, 1.5, 1, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4] + + +def terminal_size() -> tuple[int, int]: + """Find the current size of the terminal window. + + :return: Tuple[int, int]: Number of columns; number of rows + """ + columns, rows = (80, 25) + if sys.stdout.isatty(): + try: + result = subprocess.run(["stty", "size"], timeout=1, capture_output=True) + rows, columns = (int(i) for i in result.stdout.decode("utf-8").split()) + except Exception as err: + logger.exception(err) + logger.warning("Exception caught and ignored, using default size.") + pass + columns = min(columns, 120) + rows = min(rows, columns // 3) + + return columns, rows + + +def prettyprint_dictionary(d): + """ + Produce a nice string representation of a dictionary, for printing. + + :param d: Dictionary to be printed. + :type d: Dict[Optional[Any]] + :return: String representation of :param d:. + :rtype: str + """ + return "{\n%s\n}" % "\n".join( + " %s: %s" + % ( + k, + str(v.decode("latin-1") if isinstance(v, bytes) else v).replace( + "\n", "\n%s" % (" " * (4 + len(k))) + ), + ) + for k, v in d.items() + ) + + +def prettyprint_process(d): + """ + Produce a nice string representation of a subprocess CompletedProcess object, for printing. + + :param d: ReturnObject to be printed. + :return: String representation of :param d:. + :rtype: str + """ + return prettyprint_dictionary( + { + "command": d.args, + "exitcode": d.returncode, + "stdout": d.stdout, + "stderr": d.stderr, + } + ) + + +def plot_intensities( + bins, + hist_value_factor, + title="'Pixel intensity distribution'", + xlabel="'% of maximum'", + ylabel="'Number of pixels'", + xticks="", + style="with boxes", +): + """ + Create an ASCII art histogram of intensities. + + :param bins: + :param hist_value_factor: + :param title: + :param xlabel: + :param ylabel: + :param xticks: + :param style: + """ + columns, rows = terminal_size() + + command = ["gnuplot"] + plot_commands = [ + "set term dumb %d %d" % (columns, rows - 2), + "set title %s" % title, + "set xlabel %s" % xlabel, + "set ylabel %s offset character %d,0" % (ylabel, len(ylabel) // 2), + "set logscale y", + "set boxwidth %f" % hist_value_factor, + "set xtics %s out nomirror" % xticks, + "set ytics out", + "plot '-' using 1:2 title '' %s" % style, + ] + for x in sorted(bins.keys()): + plot_commands.append("%f %d" % (x * hist_value_factor, bins[x])) + plot_commands.append("e") + + debug("running %s with:\n %s\n", " ".join(command), "\n ".join(plot_commands)) + + try: + result = subprocess.run( + command, + input="\n".join(plot_commands).encode("utf-8") + b"\n", + timeout=120, + capture_output=True, + env={"LD_LIBRARY_PATH": ""}, + ) + except (OSError, subprocess.TimeoutExpired): + info(traceback.format_exc()) + warn( + "Error running gnuplot. Cannot plot intensity distribution. " + "No exit code." + ) + return + else: + debug("result = %s", prettyprint_process(result)) + + returncode = getattr(result, "returncode") + if returncode: + warn( + "Error running gnuplot. Cannot plot intensity distribution. " + "Exit code %d", + returncode, + ) + else: + star = re.compile(r"\*") + state = set() + for line in result.stdout.decode("utf-8").split("\n"): + if line.strip() != "": + stars = {m.start(0) for m in re.finditer(star, line)} + if not stars: + state = set() + else: + state |= stars + line = list(line) + for s in state: + line[s] = "*" + info("".join(line)) diff --git a/src/screen/minimum_exposure.py b/src/screen/minimum_exposure.py index 406d794..ba540f6 100644 --- a/src/screen/minimum_exposure.py +++ b/src/screen/minimum_exposure.py @@ -43,26 +43,27 @@ from scipy.optimize import curve_fit from tabulate import tabulate -import iotbx.phil +import libtbx.phil from cctbx import miller -from libtbx.phil import scope, scope_extract from dials.array_family import flex from dials.util import log -from dials.util.options import OptionParser +from dials.util.options import ArgumentParser # OptionParser deprecated anyway from dials.util.version import dials_version from dxtbx.model import ExperimentList -from screen19 import __version__, d_ticks, plot_intensities, terminal_size +from screen import __version__, d_ticks, plot_intensities, terminal_size # Custom types FloatSequence = Sequence[float] Fit = Union[np.ndarray, Iterable, int, float] +Scope = libtbx.phil.scope +ScopeExtract = libtbx.phil.scope_extract help_message = __doc__ -phil_scope = iotbx.phil.parse( +phil_scope = libtbx.phil.parse( """ verbosity = 0 .type = int(value_min=0) @@ -272,7 +273,7 @@ def wilson_plot_image( def suggest_minimum_exposure( - expts: ExperimentList, refls: flex.reflection_table, params: scope_extract + expts: ExperimentList, refls: flex.reflection_table, params: ScopeExtract ) -> None: """ Suggest an estimated minimum sufficient exposure to achieve a certain resolution. @@ -416,7 +417,7 @@ def suggest_minimum_exposure( def run( - phil: scope = phil_scope, + phil: Scope = phil_scope, args: Optional[List[str]] = None, set_up_logging: bool = False, ) -> None: @@ -434,7 +435,7 @@ def run( """ usage = "%(prog)s [options] integrated.expt integrated.refl" - parser = OptionParser( + parser = ArgumentParser( usage=usage, phil=phil, read_experiments=True, @@ -443,7 +444,7 @@ def run( epilog=help_message, ) - params, options = parser.parse_args(args=args) + params, _ = parser.parse_args(args=args) if set_up_logging: # Configure the logging diff --git a/src/screen/screen.py b/src/screen/screen.py index 537175f..b8ed9d4 100644 --- a/src/screen/screen.py +++ b/src/screen/screen.py @@ -4,14 +4,14 @@ import subprocess parser = argparse.ArgumentParser(description=__doc__) -parser.add_argument("images", help="Path to the images to process.") +parser.add_argument("images", type=str, help="Path to the images to process.") def pipeline(images): subprocess.run(["dials.import", images]) subprocess.run(["dials.find_spots", "imported.expt"]) - subprocess.run(["dials.pixel_histogram", "strong.refl"]) subprocess.run(["dials.index", "imported.expt", "strong.refl"]) + subprocess.run(["dev.dials.pixel_histogram", "indexed.refl"]) subprocess.run(["dials.refine", "indexed.expt", "indexed.refl"]) subprocess.run(["dials.integrate", "refined.expt", "refined.refl"]) subprocess.run(["screen19.minimum_exposure", "integrated.expt", "integrated.refl"]) @@ -20,7 +20,3 @@ def pipeline(images): def main(args=None): args = parser.parse_args(args) pipeline(args.images) - - -if __name__ == "__main__": - main() From 0ef17be5bc4a89f999d15168aae2cd3e7d12197b Mon Sep 17 00:00:00 2001 From: Ben Williams Date: Mon, 31 Oct 2022 15:48:44 +0000 Subject: [PATCH 11/33] Correct logging levels Also remove redundant debug log file specification. --- src/screen/minimum_exposure.py | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/src/screen/minimum_exposure.py b/src/screen/minimum_exposure.py index ba540f6..b6fd5f4 100644 --- a/src/screen/minimum_exposure.py +++ b/src/screen/minimum_exposure.py @@ -65,12 +65,12 @@ phil_scope = libtbx.phil.parse( """ - verbosity = 0 - .type = int(value_min=0) + verbosity = 1 + .type = int(value_min=1) .caption = 'Verbosity level of log output' .help = "Possible values:\n" - "\t• 0: Info log output to stdout/logfile\n" - "\t• 1: Info & debug log output to stdout/logfile" + "\t• 1: Info log output to stdout/logfile\n" + "\t• 2: Info & debug log output to stdout/logfile" minimum_exposure .caption = 'Parameters for the calculation of the lower exposure bound' { @@ -99,9 +99,6 @@ log = 'screen19.minimum_exposure.log' .type = str .caption = 'Location for the info log' - debug_log = 'screen19.minimum_exposure.debug.log' - .type = str - .caption = 'Location for the debug log' wilson_plot = 'wilson_plot' .type = str .caption = 'Filename for the Wilson plot image' From 961a61b49c2847cd0ac0607c33024ede4eea35e5 Mon Sep 17 00:00:00 2001 From: Noemi Frisina Date: Mon, 31 Oct 2022 16:08:44 +0000 Subject: [PATCH 12/33] Avoid pytest error with argparse --- tests/test_screen.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_screen.py b/tests/test_screen.py index 485ae7c..c700dcc 100644 --- a/tests/test_screen.py +++ b/tests/test_screen.py @@ -26,7 +26,7 @@ def test_screen19_help_does_not_crash(): - screen.main(args=[]) + screen.main(args=[None]) def test_minimum_exposure_help_does_not_crash(): From 6194b8a308b6f418eb97b82330bb5b8f76dde501 Mon Sep 17 00:00:00 2001 From: Noemi Frisina Date: Mon, 31 Oct 2022 16:24:00 +0000 Subject: [PATCH 13/33] Remove the last tristan bits from packaging --- .azure-pipelines/azure-pipelines.yml | 4 ++-- .azure-pipelines/ci.yml | 6 +++--- setup.cfg | 8 ++++---- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/.azure-pipelines/azure-pipelines.yml b/.azure-pipelines/azure-pipelines.yml index f8e79eb..d2f048d 100644 --- a/.azure-pipelines/azure-pipelines.yml +++ b/.azure-pipelines/azure-pipelines.yml @@ -123,8 +123,8 @@ stages: - task: TwineAuthenticate@1 displayName: Set up credentials inputs: - pythonUploadServiceConnection: pypi-tristan + pythonUploadServiceConnection: pypi-screen19 - bash: | - libtbx.python -m twine upload -r pypi-tristan --config-file $(PYPIRC_PATH) $(System.ArtifactsDirectory)/package/pypi/*.tar.gz $(System.ArtifactsDirectory)/package/pypi/*.whl + libtbx.python -m twine upload -r pypi-screen19 --config-file $(PYPIRC_PATH) $(System.ArtifactsDirectory)/package/pypi/*.tar.gz $(System.ArtifactsDirectory)/package/pypi/*.whl displayName: Publish package diff --git a/.azure-pipelines/ci.yml b/.azure-pipelines/ci.yml index 456b2dd..e671074 100644 --- a/.azure-pipelines/ci.yml +++ b/.azure-pipelines/ci.yml @@ -22,12 +22,12 @@ steps: - script: | set -eux - pip install -r "$(Pipeline.Workspace)/src/requirements_dev.txt" - pip install "$(Pipeline.Workspace)/src" + libtbx.pip install -r "$(Pipeline.Workspace)/src/requirements_dev.txt" + libtbx.pip install "$(Pipeline.Workspace)/src" displayName: Install package - script: | - PYTHONDEVMODE=1 pytest -ra --cov=tristan --cov-report=xml --cov-branch + PYTHONDEVMODE=1 pytest -ra --cov=screen19 --cov-report=xml --cov-branch displayName: Run tests workingDirectory: $(Pipeline.Workspace)/src diff --git a/setup.cfg b/setup.cfg index 5866489..87c9e82 100644 --- a/setup.cfg +++ b/setup.cfg @@ -3,7 +3,7 @@ current_version = 0.213 commit = True tag = True parse = (?P\d+)\.(?P\d+)(\.(?P\d+))? -serialize = +serialize = {major}.{minor}.{patch} {major}.{minor} @@ -62,9 +62,9 @@ where = src console_scripts = screen19 = screen.screen:main screen19.minimum_exposure = screen.minimum_exposure:main -libtbx.dispatcher.script = - screen19.minimum_exposure = screen19.minimum_exposure -libtbx.precommit = +libtbx.dispatcher.script = + screen19.minimum_exposure = screen.minimum_exposure +libtbx.precommit = screen19 = screen19 [tool:pytest] From e0a3664712b022160bbfec102b26e2f52913bf80 Mon Sep 17 00:00:00 2001 From: Noemi Frisina Date: Mon, 31 Oct 2022 16:24:26 +0000 Subject: [PATCH 14/33] Add requirements file --- requirements_dev.txt | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 requirements_dev.txt diff --git a/requirements_dev.txt b/requirements_dev.txt new file mode 100644 index 0000000..cffeec6 --- /dev/null +++ b/requirements_dev.txt @@ -0,0 +1,2 @@ +pytest +pytest-cov \ No newline at end of file From de4baf7bf728c107a41f62777b6da8eae8d3c769 Mon Sep 17 00:00:00 2001 From: Ben Williams Date: Mon, 31 Oct 2022 18:47:57 +0000 Subject: [PATCH 15/33] Correct logging levels, properly this time --- src/screen/minimum_exposure.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/screen/minimum_exposure.py b/src/screen/minimum_exposure.py index b6fd5f4..7b4b665 100644 --- a/src/screen/minimum_exposure.py +++ b/src/screen/minimum_exposure.py @@ -65,12 +65,15 @@ phil_scope = libtbx.phil.parse( """ - verbosity = 1 - .type = int(value_min=1) + verbosity = 0 + .type = int(value_min=0) .caption = 'Verbosity level of log output' .help = "Possible values:\n" - "\t• 1: Info log output to stdout/logfile\n" - "\t• 2: Info & debug log output to stdout/logfile" + "\t• 0: Info log output to stdout/logfile\n" + "\t• 1: Info log output to stdout/logfile, logfile contains timing" + "information\n" + "\t• 2: Info & debug log output to stdout/logfile, logfile contains" + "timing information" minimum_exposure .caption = 'Parameters for the calculation of the lower exposure bound' { From 2d9f569d8d27b1bce49b8b279c9c283a0b16e74f Mon Sep 17 00:00:00 2001 From: Noemi Frisina Date: Fri, 18 Nov 2022 11:19:11 +0000 Subject: [PATCH 16/33] Add libtbx dispatcher --- setup.cfg | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/setup.cfg b/setup.cfg index 87c9e82..83659df 100644 --- a/setup.cfg +++ b/setup.cfg @@ -63,7 +63,8 @@ console_scripts = screen19 = screen.screen:main screen19.minimum_exposure = screen.minimum_exposure:main libtbx.dispatcher.script = - screen19.minimum_exposure = screen.minimum_exposure + screen19 = screen19 + screen19.minimum_exposure = screen19.minimum_exposure libtbx.precommit = screen19 = screen19 From 97877fb8c58a8482c7d34777eb5c6bb709147240 Mon Sep 17 00:00:00 2001 From: Noemi Frisina Date: Mon, 21 Nov 2022 16:25:40 +0000 Subject: [PATCH 17/33] Get dials.import to work --- src/screen/screen.py | 100 +++++++++++++++++++++++++++++++++++++++---- 1 file changed, 91 insertions(+), 9 deletions(-) diff --git a/src/screen/screen.py b/src/screen/screen.py index b8ed9d4..87aacbf 100644 --- a/src/screen/screen.py +++ b/src/screen/screen.py @@ -1,22 +1,104 @@ """The main screening script.""" +from __future__ import annotations import argparse +import re import subprocess +from pathlib import Path + +import libtbx.phil + +template_pattern = re.compile(r"(.*)_(?:[0-9]*\#+).(.*)") + +import_scope = libtbx.phil.parse( + """ + include scope dials.command_line.dials_import.phil_scope + """, + process_includes=True, +) + +phil_scope = libtbx.phil.parse( + """ + log = False + .type = bool + dials_import { + include scope screen.import_scope + } + dials_find_spots{ + include scope dials.command_line.find_spots.phil_scope + } + """, + process_includes=True, +) + + +class _ImportImages(argparse.Action): + def __call__(self, parser, namespace, values, option_string=None): + if len(values) > 1: # ie. already a list + expt, directory, template = (values, None, None) + else: + expt, directory, template = self.find_import_arguments(values) + setattr(namespace, self.dest, expt) + setattr(namespace, "directory", directory) + setattr(namespace, "template", template) + + @staticmethod + def find_import_arguments(val) -> tuple[str | None]: # expts, dir, template + in_value = Path(val[0]).expanduser().resolve() + match = template_pattern.fullmatch(in_value.stem) + if in_value.is_dir() is True: + return None, in_value.as_posix(), None + elif match: + return None, None, in_value.as_posix() + else: + return in_value.as_posix(), None, None + parser = argparse.ArgumentParser(description=__doc__) -parser.add_argument("images", type=str, help="Path to the images to process.") +parser.add_argument( + "experiments", type=str, nargs="+", action=_ImportImages, help="" +) # FIXME TODO add file.cbf:1:100 +parser.add_argument("phil_args", nargs="*") + +def run_import(images, import_params): + # Ugly, but works + if images and type(images) is str: + subprocess.run(["dials.import", images, import_params.as_str()]) + elif images and type(images) is list: + subprocess.run(["dials.import", *images, import_params.as_str()]) + else: + subprocess.run(["dials.import", import_params.as_str()]) # FIXME -def pipeline(images): - subprocess.run(["dials.import", images]) + +def pipeline(args, working_phil): + params = working_phil.extract() + + # Set directory/template if that's what's been parsed. + params.dials_import.input.directory = [args.directory] if args.directory else [] + params.dials_import.input.template = [args.template] if args.template else [] + + import_params = import_scope.format(python_object=params.dials_import) + + run_import(args.experiments, import_params) subprocess.run(["dials.find_spots", "imported.expt"]) - subprocess.run(["dials.index", "imported.expt", "strong.refl"]) - subprocess.run(["dev.dials.pixel_histogram", "indexed.refl"]) - subprocess.run(["dials.refine", "indexed.expt", "indexed.refl"]) - subprocess.run(["dials.integrate", "refined.expt", "refined.refl"]) - subprocess.run(["screen19.minimum_exposure", "integrated.expt", "integrated.refl"]) + + +# def pipeline(images): +# subprocess.run(["dials.index", "imported.expt", "strong.refl"]) +# subprocess.run(["dev.dials.pixel_histogram", "indexed.refl"]) +# subprocess.run(["dials.refine", "indexed.expt", "indexed.refl"]) +# subprocess.run(["dials.integrate", "refined.expt", "refined.refl"]) +# subprocess.run(["screen19.minimum_exposure", "integrated.expt", "integrated.refl"]) def main(args=None): args = parser.parse_args(args) - pipeline(args.images) + cl = phil_scope.command_line_argument_interpreter() + working_phil = phil_scope.fetch(cl.process_and_fetch(args.phil_args)) + + pipeline(args, working_phil) + + +if __name__ == "__main__": + main(None) From e45809e028ace694788e5ce69df946f417a72c97 Mon Sep 17 00:00:00 2001 From: Noemi Frisina Date: Mon, 21 Nov 2022 16:44:56 +0000 Subject: [PATCH 18/33] Add basic find_spots and indexing --- src/screen/screen.py | 49 ++++++++++++++++++++++++++++++++++++-------- 1 file changed, 40 insertions(+), 9 deletions(-) diff --git a/src/screen/screen.py b/src/screen/screen.py index 87aacbf..f9f8bd9 100644 --- a/src/screen/screen.py +++ b/src/screen/screen.py @@ -17,6 +17,20 @@ process_includes=True, ) +find_spots_scope = libtbx.phil.parse( + """ + include scope dials.command_line.find_spots.phil_scope + """, + process_includes=True, +) + +index_scope = libtbx.phil.parse( + """ + include scope dials.command_line.index.phil_scope + """, + process_includes=True, +) + phil_scope = libtbx.phil.parse( """ log = False @@ -25,7 +39,10 @@ include scope screen.import_scope } dials_find_spots{ - include scope dials.command_line.find_spots.phil_scope + include scope screen.find_spots_scope + } + dials_index { + include scope screen.index_scope } """, process_includes=True, @@ -61,14 +78,30 @@ def find_import_arguments(val) -> tuple[str | None]: # expts, dir, template parser.add_argument("phil_args", nargs="*") -def run_import(images, import_params): +def run_import(images, params): # Ugly, but works + import_params = import_scope.format(python_object=params) + if images and type(images) is str: subprocess.run(["dials.import", images, import_params.as_str()]) elif images and type(images) is list: subprocess.run(["dials.import", *images, import_params.as_str()]) else: - subprocess.run(["dials.import", import_params.as_str()]) # FIXME + subprocess.run(["dials.import", import_params.as_str()]) + + +def run_find_spots(params): + find_spots_params = import_scope.format(python_object=params) + + subprocess.run(["dials.find_spots", "imported.expt", find_spots_params.as_str()]) + + +def run_indexing(params): + index_params = import_scope.format(python_object=params) + + subprocess.run( + ["dials.index", "imported.expt", "strong.refl", index_params.as_str()] + ) def pipeline(args, working_phil): @@ -78,15 +111,13 @@ def pipeline(args, working_phil): params.dials_import.input.directory = [args.directory] if args.directory else [] params.dials_import.input.template = [args.template] if args.template else [] - import_params = import_scope.format(python_object=params.dials_import) - - run_import(args.experiments, import_params) - subprocess.run(["dials.find_spots", "imported.expt"]) + run_import(args.experiments, params.dials_import) + run_find_spots(params.dials_find_spots) + run_indexing(params.dials_index) + subprocess.run(["dev.dials.pixel_histogram", "indexed.refl"]) # def pipeline(images): -# subprocess.run(["dials.index", "imported.expt", "strong.refl"]) -# subprocess.run(["dev.dials.pixel_histogram", "indexed.refl"]) # subprocess.run(["dials.refine", "indexed.expt", "indexed.refl"]) # subprocess.run(["dials.integrate", "refined.expt", "refined.refl"]) # subprocess.run(["screen19.minimum_exposure", "integrated.expt", "integrated.refl"]) From f26d5c14be1b4f673ebee9e70359e3c4ff8782ef Mon Sep 17 00:00:00 2001 From: Noemi Frisina Date: Mon, 21 Nov 2022 16:54:55 +0000 Subject: [PATCH 19/33] Fill in pipeline - TBC --- src/screen/screen.py | 51 ++++++++++++++++++++++++++++++++++++++------ 1 file changed, 45 insertions(+), 6 deletions(-) diff --git a/src/screen/screen.py b/src/screen/screen.py index f9f8bd9..c95832b 100644 --- a/src/screen/screen.py +++ b/src/screen/screen.py @@ -31,6 +31,20 @@ process_includes=True, ) +refine_scope = libtbx.phil.parse( + """ + include scope dials.command_line.refine.phil_scope + """, + process_includes=True, +) + +integrate_scope = libtbx.phil.parse( + """ + include scope dials.command_line.integrate.phil_scope + """, + process_includes=True, +) + phil_scope = libtbx.phil.parse( """ log = False @@ -44,6 +58,12 @@ dials_index { include scope screen.index_scope } + dials_refine { + include scope screen.refine_scope + } + dials_integrate { + include scope screen.integrate_scope + } """, process_includes=True, ) @@ -104,6 +124,27 @@ def run_indexing(params): ) +def run_refine(params): + refine_params = import_scope.format(python_object=params) + + subprocess.run( + ["dials.refine", "indexed.expt", "indexed.refl", refine_params.as_str()] + ) + + +def run_integrate(params): + integrate_params = import_scope.format(python_object=params) + + subprocess.run( + ["dials.integrate", "refined.expt", "refined.refl", integrate_params.as_str()] + ) + + +def run_minimum_exposure(): + # subprocess.run(["screen19.minimum_exposure", "integrated.expt", "integrated.refl"]) + pass + + def pipeline(args, working_phil): params = working_phil.extract() @@ -112,15 +153,13 @@ def pipeline(args, working_phil): params.dials_import.input.template = [args.template] if args.template else [] run_import(args.experiments, params.dials_import) + # FIXME at the moment after dials_import it doesn't seem to correctly read input phil parameters ... run_find_spots(params.dials_find_spots) run_indexing(params.dials_index) subprocess.run(["dev.dials.pixel_histogram", "indexed.refl"]) - - -# def pipeline(images): -# subprocess.run(["dials.refine", "indexed.expt", "indexed.refl"]) -# subprocess.run(["dials.integrate", "refined.expt", "refined.refl"]) -# subprocess.run(["screen19.minimum_exposure", "integrated.expt", "integrated.refl"]) + run_refine(params.dials_refine) + run_integrate(params.dials_integrate) + run_minimum_exposure() def main(args=None): From 16f867a3fabbf5474a0b1ab2b7b0b56e64b8b130 Mon Sep 17 00:00:00 2001 From: Noemi Frisina Date: Mon, 21 Nov 2022 17:22:24 +0000 Subject: [PATCH 20/33] Add version parser, and a slightly working show_config --- src/screen/__init__.py | 30 +++++++++++++++++++++++++++++- src/screen/minimum_exposure.py | 7 ++++++- src/screen/screen.py | 15 ++++++++++++++- 3 files changed, 49 insertions(+), 3 deletions(-) diff --git a/src/screen/__init__.py b/src/screen/__init__.py index b19c3a1..ebef663 100644 --- a/src/screen/__init__.py +++ b/src/screen/__init__.py @@ -7,15 +7,43 @@ __version_tuple__ = tuple(int(x) for x in __version__.split(".")) +import argparse import logging import re import subprocess import sys import traceback +from dials.util import version + logger = logging.getLogger("dials.screen19") debug, info, warn = logger.debug, logger.info, logger.warning -# FIXME setting verbosity to 1 doesn't seem to print any debugging. + +version_parser = argparse.ArgumentParser(add_help=False) +version_parser.add_argument( + "-V", + "--version", + action="version", + version=f"Screen19 version {__version__}, using DIALS {version.dials_version()}.", +) + +config_parser = argparse.ArgumentParser(add_help=False) +config_parser.add_argument( + "-c", + "--show-config", + action="store_true", + default=False, + dest="show_config", + help="Show the configuration parameters.", +) +config_parser.add_argument( + "-a", + "--attributes-level", + default=0, + type=int, + dest="attributes_level", + help="Set the attributes level for showing the configuration parameters.", +) # Set axis tick positions manually. Accounts for reciprocal(-square) d-scaling d_ticks = [5, 3, 2, 1.5, 1, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4] diff --git a/src/screen/minimum_exposure.py b/src/screen/minimum_exposure.py index 7b4b665..1fc52dd 100644 --- a/src/screen/minimum_exposure.py +++ b/src/screen/minimum_exposure.py @@ -52,7 +52,12 @@ from dials.util.version import dials_version from dxtbx.model import ExperimentList -from screen import __version__, d_ticks, plot_intensities, terminal_size +from screen import ( # FIXME TODO change to relative import + __version__, + d_ticks, + plot_intensities, + terminal_size, +) # Custom types FloatSequence = Sequence[float] diff --git a/src/screen/screen.py b/src/screen/screen.py index c95832b..508b320 100644 --- a/src/screen/screen.py +++ b/src/screen/screen.py @@ -4,8 +4,14 @@ import argparse import re import subprocess +import sys from pathlib import Path +from __init__ import ( # FIXME TODO change to relative import + config_parser, + version_parser, +) + import libtbx.phil template_pattern = re.compile(r"(.*)_(?:[0-9]*\#+).(.*)") @@ -91,7 +97,9 @@ def find_import_arguments(val) -> tuple[str | None]: # expts, dir, template return in_value.as_posix(), None, None -parser = argparse.ArgumentParser(description=__doc__) +parser = argparse.ArgumentParser( + description=__doc__, parents=[version_parser, config_parser] +) parser.add_argument( "experiments", type=str, nargs="+", action=_ImportImages, help="" ) # FIXME TODO add file.cbf:1:100 @@ -167,6 +175,11 @@ def main(args=None): cl = phil_scope.command_line_argument_interpreter() working_phil = phil_scope.fetch(cl.process_and_fetch(args.phil_args)) + if args.show_config: + # FIXME doesn't work unless some words are passed as positional arguments + working_phil.show(attributes_level=args.attributes_level) + sys.exit() + pipeline(args, working_phil) From aa341e469a1d9883eba752315f5e86b0381392c8 Mon Sep 17 00:00:00 2001 From: Noemi Frisina Date: Fri, 9 Dec 2022 11:49:20 +0000 Subject: [PATCH 21/33] Hopefully fix imports --- .pre-commit-config.yaml | 4 ++-- src/screen/screen.py | 17 +++++++---------- 2 files changed, 9 insertions(+), 12 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index a852f41..c483a09 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -23,8 +23,8 @@ repos: args: [--safe, --quiet] # Enforce style with Flake8 -- repo: https://gitlab.com/pycqa/flake8 - rev: 3.9.2 +- repo: https://github.com/pycqa/flake8 + rev: 4.0.1 hooks: - id: flake8 args: ['--max-line-length=88', '--select=E401,E711,E712,E713,E714,E721,E722,E901,F401,F402,F403,F405,F631,F632,F633,F811,F812,F821,F822,F841,F901,W191,W291,W292,W293,W602,W603,W604,W605,W606'] diff --git a/src/screen/screen.py b/src/screen/screen.py index 508b320..4fa1a02 100644 --- a/src/screen/screen.py +++ b/src/screen/screen.py @@ -7,13 +7,10 @@ import sys from pathlib import Path -from __init__ import ( # FIXME TODO change to relative import - config_parser, - version_parser, -) - import libtbx.phil +from screen import config_parser, version_parser # FIXME TODO change to relative import + template_pattern = re.compile(r"(.*)_(?:[0-9]*\#+).(.*)") import_scope = libtbx.phil.parse( @@ -56,19 +53,19 @@ log = False .type = bool dials_import { - include scope screen.import_scope + include scope screen.screen.import_scope } dials_find_spots{ - include scope screen.find_spots_scope + include scope screen.screen.find_spots_scope } dials_index { - include scope screen.index_scope + include scope screen.screen.index_scope } dials_refine { - include scope screen.refine_scope + include scope screen.screen.refine_scope } dials_integrate { - include scope screen.integrate_scope + include scope screen.screen.integrate_scope } """, process_includes=True, From 5a2086294836a44ebccc1bae887104f687f75d65 Mon Sep 17 00:00:00 2001 From: Noemi Frisina Date: Fri, 9 Dec 2022 13:57:03 +0000 Subject: [PATCH 22/33] Get correct scopes --- src/screen/screen.py | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/src/screen/screen.py b/src/screen/screen.py index 4fa1a02..93945b0 100644 --- a/src/screen/screen.py +++ b/src/screen/screen.py @@ -48,24 +48,25 @@ process_includes=True, ) + phil_scope = libtbx.phil.parse( """ log = False .type = bool dials_import { - include scope screen.screen.import_scope + include scope screen.import_scope } dials_find_spots{ - include scope screen.screen.find_spots_scope + include scope screen.find_spots_scope } dials_index { - include scope screen.screen.index_scope + include scope screen.index_scope } dials_refine { - include scope screen.screen.refine_scope + include scope screen.refine_scope } dials_integrate { - include scope screen.screen.integrate_scope + include scope screen.integrate_scope } """, process_includes=True, @@ -116,13 +117,13 @@ def run_import(images, params): def run_find_spots(params): - find_spots_params = import_scope.format(python_object=params) + find_spots_params = find_spots_scope.format(python_object=params) subprocess.run(["dials.find_spots", "imported.expt", find_spots_params.as_str()]) def run_indexing(params): - index_params = import_scope.format(python_object=params) + index_params = index_scope.format(python_object=params) subprocess.run( ["dials.index", "imported.expt", "strong.refl", index_params.as_str()] @@ -130,7 +131,7 @@ def run_indexing(params): def run_refine(params): - refine_params = import_scope.format(python_object=params) + refine_params = refine_scope.format(python_object=params) subprocess.run( ["dials.refine", "indexed.expt", "indexed.refl", refine_params.as_str()] @@ -138,7 +139,7 @@ def run_refine(params): def run_integrate(params): - integrate_params = import_scope.format(python_object=params) + integrate_params = integrate_scope.format(python_object=params) subprocess.run( ["dials.integrate", "refined.expt", "refined.refl", integrate_params.as_str()] From ed2e3f674b095d6a46ca9bc12556d0834663cb01 Mon Sep 17 00:00:00 2001 From: Noemi Frisina Date: Fri, 9 Dec 2022 14:00:26 +0000 Subject: [PATCH 23/33] Try different scope path --- src/screen/screen.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/screen/screen.py b/src/screen/screen.py index 93945b0..f7d72a0 100644 --- a/src/screen/screen.py +++ b/src/screen/screen.py @@ -54,19 +54,19 @@ log = False .type = bool dials_import { - include scope screen.import_scope + include scope screen.screen.import_scope } dials_find_spots{ - include scope screen.find_spots_scope + include scope screen.screen.find_spots_scope } dials_index { - include scope screen.index_scope + include scope screen.screen.index_scope } dials_refine { - include scope screen.refine_scope + include scope screen.screen.refine_scope } dials_integrate { - include scope screen.integrate_scope + include scope screen.screen.integrate_scope } """, process_includes=True, From 0187d8acbf8d74c432bc21e72ec7bc651a114f3f Mon Sep 17 00:00:00 2001 From: Noemi Frisina Date: Fri, 9 Dec 2022 14:04:16 +0000 Subject: [PATCH 24/33] Ugly scope workaround --- src/screen/scopes.py | 36 ++++++++++++++++++++++++++++++ src/screen/screen.py | 52 ++++++++++---------------------------------- 2 files changed, 48 insertions(+), 40 deletions(-) create mode 100644 src/screen/scopes.py diff --git a/src/screen/scopes.py b/src/screen/scopes.py new file mode 100644 index 0000000..2469d97 --- /dev/null +++ b/src/screen/scopes.py @@ -0,0 +1,36 @@ +import libtbx.phil + +import_scope = libtbx.phil.parse( + """ + include scope dials.command_line.dials_import.phil_scope + """, + process_includes=True, +) + +find_spots_scope = libtbx.phil.parse( + """ + include scope dials.command_line.find_spots.phil_scope + """, + process_includes=True, +) + +index_scope = libtbx.phil.parse( + """ + include scope dials.command_line.index.phil_scope + """, + process_includes=True, +) + +refine_scope = libtbx.phil.parse( + """ + include scope dials.command_line.refine.phil_scope + """, + process_includes=True, +) + +integrate_scope = libtbx.phil.parse( + """ + include scope dials.command_line.integrate.phil_scope + """, + process_includes=True, +) diff --git a/src/screen/screen.py b/src/screen/screen.py index f7d72a0..3feea5f 100644 --- a/src/screen/screen.py +++ b/src/screen/screen.py @@ -10,43 +10,15 @@ import libtbx.phil from screen import config_parser, version_parser # FIXME TODO change to relative import - -template_pattern = re.compile(r"(.*)_(?:[0-9]*\#+).(.*)") - -import_scope = libtbx.phil.parse( - """ - include scope dials.command_line.dials_import.phil_scope - """, - process_includes=True, -) - -find_spots_scope = libtbx.phil.parse( - """ - include scope dials.command_line.find_spots.phil_scope - """, - process_includes=True, -) - -index_scope = libtbx.phil.parse( - """ - include scope dials.command_line.index.phil_scope - """, - process_includes=True, +from screen.scopes import ( + find_spots_scope, + import_scope, + index_scope, + integrate_scope, + refine_scope, ) -refine_scope = libtbx.phil.parse( - """ - include scope dials.command_line.refine.phil_scope - """, - process_includes=True, -) - -integrate_scope = libtbx.phil.parse( - """ - include scope dials.command_line.integrate.phil_scope - """, - process_includes=True, -) +template_pattern = re.compile(r"(.*)_(?:[0-9]*\#+).(.*)") phil_scope = libtbx.phil.parse( @@ -54,19 +26,19 @@ log = False .type = bool dials_import { - include scope screen.screen.import_scope + include scope screen.scopes.import_scope } dials_find_spots{ - include scope screen.screen.find_spots_scope + include scope screen.scopes.find_spots_scope } dials_index { - include scope screen.screen.index_scope + include scope screen.scopes.index_scope } dials_refine { - include scope screen.screen.refine_scope + include scope screen.scopes.refine_scope } dials_integrate { - include scope screen.screen.integrate_scope + include scope screen.scopes.integrate_scope } """, process_includes=True, From aa6b62d4582d2feb83eae0e3a592e4d9ae5665fe Mon Sep 17 00:00:00 2001 From: Noemi Frisina Date: Fri, 9 Dec 2022 18:39:54 +0000 Subject: [PATCH 25/33] Update parser to accept also image range as input eg. file_00.cbf:0:10 --- src/screen/screen.py | 32 +++++++++++++++++++++++--------- 1 file changed, 23 insertions(+), 9 deletions(-) diff --git a/src/screen/screen.py b/src/screen/screen.py index 3feea5f..9e8f4c4 100644 --- a/src/screen/screen.py +++ b/src/screen/screen.py @@ -47,24 +47,35 @@ class _ImportImages(argparse.Action): def __call__(self, parser, namespace, values, option_string=None): - if len(values) > 1: # ie. already a list - expt, directory, template = (values, None, None) - else: - expt, directory, template = self.find_import_arguments(values) + expt, directory, template, image_range = self.find_import_arguments(values) setattr(namespace, self.dest, expt) setattr(namespace, "directory", directory) setattr(namespace, "template", template) + setattr(namespace, "image_range", image_range) @staticmethod - def find_import_arguments(val) -> tuple[str | None]: # expts, dir, template + def find_import_arguments( + val, + ) -> tuple[str | None]: # expts, dir, template, image_range + if len(val) > 1: + return val, None, None, None + in_value = Path(val[0]).expanduser().resolve() match = template_pattern.fullmatch(in_value.stem) if in_value.is_dir() is True: - return None, in_value.as_posix(), None + return None, in_value.as_posix(), None, None elif match: - return None, None, in_value.as_posix() + return None, None, in_value.as_posix(), None + elif ":" in in_value.as_posix(): + if len(in_value.name.split(":")) != 3: + raise OSError(5, "Please specify both start and end of image range.") + filename, start, end = in_value.name.split(":") + temp = re.split(r"([0-9#]+)(?=\.\w)", filename)[1] + template = in_value.parent / filename.replace(temp, "#" * len(temp)) + image_range = (int(start), int(end)) + return None, None, template.as_posix(), image_range else: - return in_value.as_posix(), None, None + return in_value.as_posix(), None, None, None parser = argparse.ArgumentParser( @@ -74,6 +85,7 @@ def find_import_arguments(val) -> tuple[str | None]: # expts, dir, template "experiments", type=str, nargs="+", action=_ImportImages, help="" ) # FIXME TODO add file.cbf:1:100 parser.add_argument("phil_args", nargs="*") +options_group = parser.add_argument_group("") def run_import(images, params): @@ -129,9 +141,11 @@ def pipeline(args, working_phil): # Set directory/template if that's what's been parsed. params.dials_import.input.directory = [args.directory] if args.directory else [] params.dials_import.input.template = [args.template] if args.template else [] + params.dials_import.geometry.scan.image_range = ( + args.image_range if args.image_range else None + ) run_import(args.experiments, params.dials_import) - # FIXME at the moment after dials_import it doesn't seem to correctly read input phil parameters ... run_find_spots(params.dials_find_spots) run_indexing(params.dials_index) subprocess.run(["dev.dials.pixel_histogram", "indexed.refl"]) From b935dfcb2326302f7d8ec2db3a3fcc5558f503ee Mon Sep 17 00:00:00 2001 From: Noemi Frisina Date: Mon, 12 Dec 2022 11:59:24 +0000 Subject: [PATCH 26/33] Workaround for additional options in dials pipeline --- src/screen/{scopes.py => inputs.py} | 28 +++++++++++++++ src/screen/screen.py | 55 +++++++++++++++++++---------- 2 files changed, 64 insertions(+), 19 deletions(-) rename src/screen/{scopes.py => inputs.py} (53%) diff --git a/src/screen/scopes.py b/src/screen/inputs.py similarity index 53% rename from src/screen/scopes.py rename to src/screen/inputs.py index 2469d97..fda57f3 100644 --- a/src/screen/scopes.py +++ b/src/screen/inputs.py @@ -1,3 +1,5 @@ +import argparse + import libtbx.phil import_scope = libtbx.phil.parse( @@ -34,3 +36,29 @@ """, process_includes=True, ) + +options_parser = argparse.ArgumentParser(add_help=False) +options_parser.add_argument( + "--find-spots", + type=str, + default=[], + nargs="*", + help="Additional options for spot finding.", +) +options_parser.add_argument( + "--index", type=str, default=[], nargs="*", help="Additional options for indexing." +) +options_parser.add_argument( + "--refine", + type=str, + default=[], + nargs="*", + help="Additional options for refinement.", +) +options_parser.add_argument( + "--integrate", + type=str, + default=[], + nargs="*", + help="Additional options for integration.", +) diff --git a/src/screen/screen.py b/src/screen/screen.py index 9e8f4c4..eae0fd8 100644 --- a/src/screen/screen.py +++ b/src/screen/screen.py @@ -10,11 +10,12 @@ import libtbx.phil from screen import config_parser, version_parser # FIXME TODO change to relative import -from screen.scopes import ( +from screen.inputs import ( find_spots_scope, import_scope, index_scope, integrate_scope, + options_parser, refine_scope, ) @@ -79,13 +80,10 @@ def find_import_arguments( parser = argparse.ArgumentParser( - description=__doc__, parents=[version_parser, config_parser] + description=__doc__, parents=[version_parser, config_parser, options_parser] ) -parser.add_argument( - "experiments", type=str, nargs="+", action=_ImportImages, help="" -) # FIXME TODO add file.cbf:1:100 +parser.add_argument("experiments", type=str, nargs="+", action=_ImportImages, help="") parser.add_argument("phil_args", nargs="*") -options_group = parser.add_argument_group("") def run_import(images, params): @@ -100,33 +98,47 @@ def run_import(images, params): subprocess.run(["dials.import", import_params.as_str()]) -def run_find_spots(params): +def run_find_spots(params, options=[]): find_spots_params = find_spots_scope.format(python_object=params) - subprocess.run(["dials.find_spots", "imported.expt", find_spots_params.as_str()]) + subprocess.run( + ["dials.find_spots", "imported.expt", find_spots_params.as_str(), *options] + ) -def run_indexing(params): +def run_indexing(params, options=[]): index_params = index_scope.format(python_object=params) subprocess.run( - ["dials.index", "imported.expt", "strong.refl", index_params.as_str()] + ["dials.index", "imported.expt", "strong.refl", index_params.as_str(), *options] ) -def run_refine(params): +def run_refine(params, options=[]): refine_params = refine_scope.format(python_object=params) subprocess.run( - ["dials.refine", "indexed.expt", "indexed.refl", refine_params.as_str()] + [ + "dials.refine", + "indexed.expt", + "indexed.refl", + refine_params.as_str(), + *options, + ] ) -def run_integrate(params): +def run_integrate(params, options=[]): integrate_params = integrate_scope.format(python_object=params) subprocess.run( - ["dials.integrate", "refined.expt", "refined.refl", integrate_params.as_str()] + [ + "dials.integrate", + "refined.expt", + "refined.refl", + integrate_params.as_str(), + *options, + ] ) @@ -145,13 +157,18 @@ def pipeline(args, working_phil): args.image_range if args.image_range else None ) + spot_finding_options = args.find_spots + indexing_options = args.index + refinement_options = args.refine + integration_options = args.integrate + run_import(args.experiments, params.dials_import) - run_find_spots(params.dials_find_spots) - run_indexing(params.dials_index) + run_find_spots(params.dials_find_spots, spot_finding_options) + run_indexing(params.dials_index, indexing_options) subprocess.run(["dev.dials.pixel_histogram", "indexed.refl"]) - run_refine(params.dials_refine) - run_integrate(params.dials_integrate) - run_minimum_exposure() + run_refine(params.dials_refine, refinement_options) + run_integrate(params.dials_integrate, integration_options) + run_minimum_exposure() # If minimum exposure at indexing, stope there, else go on to integrate def main(args=None): From 6175996bf11ed2eb2e035a9159dbd15eecfa711b Mon Sep 17 00:00:00 2001 From: Noemi Frisina Date: Mon, 12 Dec 2022 12:04:11 +0000 Subject: [PATCH 27/33] Fix the names --- src/screen/screen.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/screen/screen.py b/src/screen/screen.py index eae0fd8..dc6a0db 100644 --- a/src/screen/screen.py +++ b/src/screen/screen.py @@ -27,19 +27,19 @@ log = False .type = bool dials_import { - include scope screen.scopes.import_scope + include scope screen.inputs.import_scope } dials_find_spots{ - include scope screen.scopes.find_spots_scope + include scope screen.inputs.find_spots_scope } dials_index { - include scope screen.scopes.index_scope + include scope screen.inputs.index_scope } dials_refine { - include scope screen.scopes.refine_scope + include scope screen.inputs.refine_scope } dials_integrate { - include scope screen.scopes.integrate_scope + include scope screen.inputs.integrate_scope } """, process_includes=True, From 01a81a283536d2b1888f99777d4e60e4d1f0e87a Mon Sep 17 00:00:00 2001 From: Noemi Frisina Date: Mon, 12 Dec 2022 16:53:26 +0000 Subject: [PATCH 28/33] Get minimum exposure to run - at least for integration --- src/screen/minimum_exposure.py | 215 ++++++++++++++------------------- src/screen/screen.py | 81 ++++++++++--- 2 files changed, 156 insertions(+), 140 deletions(-) diff --git a/src/screen/minimum_exposure.py b/src/screen/minimum_exposure.py index 1fc52dd..ed4323f 100644 --- a/src/screen/minimum_exposure.py +++ b/src/screen/minimum_exposure.py @@ -1,43 +1,46 @@ """ Perform straight-line Wilson plot fit. Draw the Wilson plot. +""" +help_message = """ Reflection d-spacings are determined from the crystal symmetry (from indexing) and the Miller indices of the indexed reflections. The atomic displacement parameter is assumed isotropic. Its value is -determined from a fit to the reflection data: - I = A * exp(-B / (2 * d^2)), +determined from a fit to the reflection data: \n + I = A * exp(-B / (2 * d^2)), \n where I is the intensity and the scale factor, A, and isotropic -displacement parameter, B, are the fitted parameters. - +displacement parameter, B, are the fitted parameters. \n +\n An I/σ condition for 'good' diffraction statistics is set by the instance variable min_i_over_sigma, and the user's desired resolution is set by the instance variable desired_d. A crude error model is assumed, whereby σ² = I, and so the I/σ condition -translates trivially to a threshold I. - +translates trivially to a threshold I. \n +\n The value of the fitted intensity function at the desired resolution is compared with the threshold I. The ratio of these values is used to determine a recommended exposure (flux × exposure time) -for the full data collection. - +for the full data collection. \n +\n The Wilson plot of I as a function of d is drawn as the file 'wilson_plot.png'. The plot can optionally be saved in other formats. +\n +Examples:\n -Examples: - - screen19.minimum_exposure integrated.expt integrated.refl + screen19.minimum_exposure integrated.expt integrated.refl \n - screen19.minimum_exposure indexed.expt indexed.refl + screen19.minimum_exposure indexed.expt indexed.refl \n screen19.minimum_exposure min_i_over_sigma=2 desired_d=0.84 wilson_fit_max_d=4 \ - integrated.expt integrated.refl + integrated.expt integrated.refl \n """ +import argparse import logging -import time +import sys from io import StringIO -from typing import Iterable, List, Optional, Sequence, Union +from typing import Iterable, Optional, Sequence, Union import numpy as np from scipy.optimize import curve_fit @@ -48,15 +51,17 @@ from dials.array_family import flex from dials.util import log -from dials.util.options import ArgumentParser # OptionParser deprecated anyway -from dials.util.version import dials_version + +# from dials.util.options import ArgumentParser # OptionParser deprecated anyway +# from dials.util.version import dials_version from dxtbx.model import ExperimentList -from screen import ( # FIXME TODO change to relative import - __version__, +from screen import ( # FIXME TODO change to relative import; __version__, + config_parser, d_ticks, plot_intensities, terminal_size, + version_parser, ) # Custom types @@ -65,56 +70,49 @@ Scope = libtbx.phil.scope ScopeExtract = libtbx.phil.scope_extract -help_message = __doc__ - phil_scope = libtbx.phil.parse( """ verbosity = 0 - .type = int(value_min=0) - .caption = 'Verbosity level of log output' - .help = "Possible values:\n" - "\t• 0: Info log output to stdout/logfile\n" - "\t• 1: Info log output to stdout/logfile, logfile contains timing" - "information\n" - "\t• 2: Info & debug log output to stdout/logfile, logfile contains" - "timing information" + .type = int(value_min=0) + .help = "Verbosity level of log output. Possible values:\n" + "\t• 0: Info log output to stdout/logfile\n" + "\t• 1: Info log output to stdout/logfile, logfile contains timing" + "information\n" + "\t• 2: Info & debug log output to stdout/logfile, logfile contains" + "timing information" minimum_exposure - .caption = 'Parameters for the calculation of the lower exposure bound' - { - desired_d = None - .multiple = True - .type = float - .caption = u'Desired resolution limit, in Ångströms, of diffraction data' - .help = 'This is the resolution target for the lower-bound exposure ' \ - 'recommendation.' - min_i_over_sigma = 2 - .type = float - .caption = u'Target I/σ value for lower-bound exposure recommendation' - .help = 'The lower-bound exposure recommendation provides an estimate of ' \ - u'the exposure (flux × exposure time) required to ensure that the' \ - 'majority of expected reflections at the desired resolution limit' \ - u'have I/σ greater than or equal to this value.' - wilson_fit_max_d = 4 # Å - .type = float - .caption = u'Maximum d-value (in Ångströms) for displacement parameter fit' - .help = 'Reflections with lower resolution than this value will be ' \ - 'ignored for the purposes of the Wilson plot.' - } - output - .caption = 'Parameters to control the output' - { - log = 'screen19.minimum_exposure.log' - .type = str - .caption = 'Location for the info log' - wilson_plot = 'wilson_plot' - .type = str - .caption = 'Filename for the Wilson plot image' - .help = "By default, the extension '.png' is appended. If you include " \ - "a different extension, either '.pdf', '.ps', '.eps' or '.svg', " \ - "a file of that format will be created instead." - } - """, + .caption = 'Parameters for the calculation of the lower exposure bound' + { + desired_d = None + .multiple = True + .type = float + .help = u'Desired resolution limit, in Ångströms, of diffraction data' \ + 'This is the resolution target for the lower-bound exposure recommendation.' + min_i_over_sigma = 2 + .type = float + .help = u'Target I/σ value for lower-bound exposure recommendation' \ + 'The lower-bound exposure recommendation provides an estimate of ' \ + u'the exposure (flux × exposure time) required to ensure that the' \ + 'majority of expected reflections at the desired resolution limit' \ + u'have I/σ greater than or equal to this value.' + wilson_fit_max_d = 4 # Å + .type = float + .help = u'Maximum d-value (in Ångströms) for displacement parameter fit' \ + 'Reflections with lower resolution than this value will be ' \ + 'ignored for the purposes of the Wilson plot.' + } + output { + log = 'screen19.minimum_exposure.log' + .type = str + .help = 'Location for the info log' + wilson_plot = 'wilson_plot' + .type = str + .help = "Filename for the Wilson plot image. By default, the extension '.png' is appended. " \ + "If you include a different extension, either '.pdf', '.ps', '.eps' or '.svg', " \ + "a file of that format will be created instead." + } + """, process_includes=True, ) @@ -320,7 +318,7 @@ def suggest_minimum_exposure( # Apply French-Wilson scaling to ensure positive intensities. miller_array = miller.array( miller.set( - expts[0].crystal.get_crystal_symmetry(), + expts.crystal.get_crystal_symmetry(), refls["miller_index"], anomalous_flag=False, ), @@ -421,77 +419,52 @@ def suggest_minimum_exposure( ) -def run( - phil: Scope = phil_scope, - args: Optional[List[str]] = None, - set_up_logging: bool = False, -) -> None: - """ - Parse command-line arguments, run the script. +usage = "%(prog)s [options] integrated.expt integrated.refl" +parser = argparse.ArgumentParser( + usage=usage, + formatter_class=argparse.RawTextHelpFormatter, + description=__doc__, + epilog=help_message, + parents=[version_parser, config_parser], +) +parser.add_argument("experiments", type=str, help="") +parser.add_argument("reflections", type=str, help="") +parser.add_argument("phil_args", nargs="*", help="") + - Uses the DIALS option parser to extract an experiment list, reflection table and - parameters, then passes them to :func:`suggest_minimum_exposure`. - Optionally, sets up the logger. +def run(set_up_logging: bool = False): + """ + Parse command-line arguments, run the script. Optionally, sets up the logger. Args: - phil: PHIL scope for option parser. - args: Arguments to parse. If None, :data:`sys.argv[1:]` will be used. - set_up_logging: Choose whether to configure :module:`screen19` logging. + set_up_logging (bool, optional): Choose whether to configure `screen19` logging. Defaults to False. """ - usage = "%(prog)s [options] integrated.expt integrated.refl" - - parser = ArgumentParser( - usage=usage, - phil=phil, - read_experiments=True, - read_reflections=True, - check_format=False, - epilog=help_message, - ) + args = parser.parse_args() + cl = phil_scope.command_line_argument_interpreter() + working_phil = phil_scope.fetch(cl.process_and_fetch(args.phil_args)) + + if args.show_config: + # FIXME doesn't work unless some words are passed as positional argument (experiments) + working_phil.show(attributes_level=args.attributes_level) + sys.exit() - params, _ = parser.parse_args(args=args) + params = working_phil.extract() if set_up_logging: # Configure the logging log.config(params.verbosity, params.output.log) - if not (params.input.experiments and params.input.reflections): - version_information = ( - f"screen19.minimum_exposure v{__version__} using {dials_version()} " - f"({time.strftime('%Y-%m-%d %H:%M:%S')})" - ) - - print(help_message) - print(version_information) - return - - if len(params.input.experiments) > 1: - warn( - "You provided more than one experiment list (%s). Only the " - "first, %s, will be used.", - ", ".join([expt.filename for expt in params.input.experiments]), - params.input.experiments[0].filename, - ) - if len(params.input.reflections) > 1: - warn( - "You provided more than one reflection table (%s). Only the " - "first, %s, will be used.", - ", ".join([refls.filename for refls in params.input.reflections]), - params.input.reflections[0].filename, - ) - - expts = params.input.experiments[0].data - refls = params.input.reflections[0].data + expt = ExperimentList.from_file(args.experiments) + refl = flex.reflection_table.from_file(args.reflections) - if len(expts) > 1: + if len(expt) > 1: warn( - "The experiment list you provided, %s, contains more than one " + f"The experiment list you provided, {args.experiments}, contains more than one " "experiment object (perhaps multiple indexing solutions). Only " - "the first will be used, all others will be ignored.", - params.input.experiments[0].filename, + "the first will be used, all others will be ignored." ) - suggest_minimum_exposure(expts, refls, params) + suggest_minimum_exposure(expt[0], refl, params) def main() -> None: diff --git a/src/screen/screen.py b/src/screen/screen.py index dc6a0db..ffa479f 100644 --- a/src/screen/screen.py +++ b/src/screen/screen.py @@ -19,13 +19,29 @@ refine_scope, ) +# Custom types +Scope = libtbx.phil.scope +ScopeExtract = libtbx.phil.scope_extract + template_pattern = re.compile(r"(.*)_(?:[0-9]*\#+).(.*)") phil_scope = libtbx.phil.parse( """ - log = False - .type = bool + verbosity = 0 + .type = int(value_min=0) + .multiple = True + .help = "Verbosity level of log output. Possible values:\n" + "\t• 0: Info log output to stdout/logfile\n" + "\t• 1: Info log output to stdout/logfile, logfile contains timing" + "information\n" + "\t• 2: Info & debug log output to stdout/logfile, logfile contains" + "timing information" + output { + log = screen19.log + .type = str + .help = The log file name. + } dials_import { include scope screen.inputs.import_scope } @@ -41,6 +57,7 @@ dials_integrate { include scope screen.inputs.integrate_scope } + include scope screen.minimum_exposure.phil_scope """, process_includes=True, ) @@ -82,11 +99,27 @@ def find_import_arguments( parser = argparse.ArgumentParser( description=__doc__, parents=[version_parser, config_parser, options_parser] ) -parser.add_argument("experiments", type=str, nargs="+", action=_ImportImages, help="") -parser.add_argument("phil_args", nargs="*") +parser.add_argument( + "experiments", + type=str, + nargs="+", + action=_ImportImages, + help="The experiment path - either a directory or a list of files.", +) +parser.add_argument( + "phil_args", nargs="*", help="Phil parameters for pipeline." +) # I think at this point this might only be needed for import and minimum_exposure??? +parser.add_argument( + "-d", + "--data", + type=str, + choices=["indexed", "integrated"], + default="integrated", + help="At which point the disorder parameter fit should be conducted. Defaults to 'integrated'.", +) -def run_import(images, params): +def run_import(images, params: ScopeExtract): # Ugly, but works import_params = import_scope.format(python_object=params) @@ -98,7 +131,7 @@ def run_import(images, params): subprocess.run(["dials.import", import_params.as_str()]) -def run_find_spots(params, options=[]): +def run_find_spots(params: ScopeExtract, options: list = []): find_spots_params = find_spots_scope.format(python_object=params) subprocess.run( @@ -106,7 +139,7 @@ def run_find_spots(params, options=[]): ) -def run_indexing(params, options=[]): +def run_indexing(params: ScopeExtract, options: list = []): index_params = index_scope.format(python_object=params) subprocess.run( @@ -114,7 +147,7 @@ def run_indexing(params, options=[]): ) -def run_refine(params, options=[]): +def run_refine(params: ScopeExtract, options: list = []): refine_params = refine_scope.format(python_object=params) subprocess.run( @@ -128,7 +161,7 @@ def run_refine(params, options=[]): ) -def run_integrate(params, options=[]): +def run_integrate(params: ScopeExtract, options: list = []): integrate_params = integrate_scope.format(python_object=params) subprocess.run( @@ -142,12 +175,17 @@ def run_integrate(params, options=[]): ) -def run_minimum_exposure(): - # subprocess.run(["screen19.minimum_exposure", "integrated.expt", "integrated.refl"]) - pass +def run_minimum_exposure(choice): + if choice == "indexed": + # subprocess.run + pass + else: + subprocess.run( + ["screen19.minimum_exposure", "integrated.expt", "integrated.refl"] + ) -def pipeline(args, working_phil): +def pipeline(args: argparse.Namespace, working_phil: Scope): params = working_phil.extract() # Set directory/template if that's what's been parsed. @@ -157,6 +195,8 @@ def pipeline(args, working_phil): args.image_range if args.image_range else None ) + print(params.minimum_exposure.desired_d) # SIGH + spot_finding_options = args.find_spots indexing_options = args.index refinement_options = args.refine @@ -166,9 +206,12 @@ def pipeline(args, working_phil): run_find_spots(params.dials_find_spots, spot_finding_options) run_indexing(params.dials_index, indexing_options) subprocess.run(["dev.dials.pixel_histogram", "indexed.refl"]) - run_refine(params.dials_refine, refinement_options) - run_integrate(params.dials_integrate, integration_options) - run_minimum_exposure() # If minimum exposure at indexing, stope there, else go on to integrate + if args.data == "integrated": + run_refine(params.dials_refine, refinement_options) + run_integrate(params.dials_integrate, integration_options) + run_minimum_exposure(args.data) + else: + run_minimum_exposure(args.data) def main(args=None): @@ -177,12 +220,12 @@ def main(args=None): working_phil = phil_scope.fetch(cl.process_and_fetch(args.phil_args)) if args.show_config: - # FIXME doesn't work unless some words are passed as positional arguments + # FIXME doesn't work unless some words are passed as positional argument (experiments) working_phil.show(attributes_level=args.attributes_level) sys.exit() pipeline(args, working_phil) -if __name__ == "__main__": - main(None) +# if __name__ == "__main__": +# main(None) From ef9d52d485c6e8a9bf4fcbdadfaefe421b785bb0 Mon Sep 17 00:00:00 2001 From: Noemi Frisina Date: Mon, 12 Dec 2022 17:17:19 +0000 Subject: [PATCH 29/33] Add some HOWTO info --- src/screen/inputs.py | 12 ++++--- src/screen/minimum_exposure.py | 31 +++++++++--------- src/screen/screen.py | 60 +++++++++++++++++++++++++++++----- 3 files changed, 75 insertions(+), 28 deletions(-) diff --git a/src/screen/inputs.py b/src/screen/inputs.py index fda57f3..7dd42e2 100644 --- a/src/screen/inputs.py +++ b/src/screen/inputs.py @@ -43,22 +43,26 @@ type=str, default=[], nargs="*", - help="Additional options for spot finding.", + help="Additional options for dials.find_spots.", ) options_parser.add_argument( - "--index", type=str, default=[], nargs="*", help="Additional options for indexing." + "--index", + type=str, + default=[], + nargs="*", + help="Additional options for dials.index.", ) options_parser.add_argument( "--refine", type=str, default=[], nargs="*", - help="Additional options for refinement.", + help="Additional options for dials.refine.", ) options_parser.add_argument( "--integrate", type=str, default=[], nargs="*", - help="Additional options for integration.", + help="Additional options for dials.integrate.", ) diff --git a/src/screen/minimum_exposure.py b/src/screen/minimum_exposure.py index ed4323f..cbf1e2f 100644 --- a/src/screen/minimum_exposure.py +++ b/src/screen/minimum_exposure.py @@ -2,6 +2,8 @@ Perform straight-line Wilson plot fit. Draw the Wilson plot. """ +from __future__ import annotations + help_message = """ Reflection d-spacings are determined from the crystal symmetry (from indexing) and the Miller indices of the indexed reflections. The @@ -40,7 +42,7 @@ import logging import sys from io import StringIO -from typing import Iterable, Optional, Sequence, Union +from typing import Iterable, Sequence, Union import numpy as np from scipy.optimize import curve_fit @@ -51,12 +53,9 @@ from dials.array_family import flex from dials.util import log - -# from dials.util.options import ArgumentParser # OptionParser deprecated anyway -# from dials.util.version import dials_version from dxtbx.model import ExperimentList -from screen import ( # FIXME TODO change to relative import; __version__, +from screen import ( # FIXME TODO change to relative import config_parser, d_ticks, plot_intensities, @@ -73,14 +72,6 @@ phil_scope = libtbx.phil.parse( """ - verbosity = 0 - .type = int(value_min=0) - .help = "Verbosity level of log output. Possible values:\n" - "\t• 0: Info log output to stdout/logfile\n" - "\t• 1: Info log output to stdout/logfile, logfile contains timing" - "information\n" - "\t• 2: Info & debug log output to stdout/logfile, logfile contains" - "timing information" minimum_exposure .caption = 'Parameters for the calculation of the lower exposure bound' { @@ -103,6 +94,14 @@ 'ignored for the purposes of the Wilson plot.' } output { + verbosity = 0 + .type = int(value_min=0) + .help = "Verbosity level of log output. Possible values:\n" + "\t• 0: Info log output to stdout/logfile\n" + "\t• 1: Info log output to stdout/logfile, logfile contains timing" + "information\n" + "\t• 2: Info & debug log output to stdout/logfile, logfile contains" + "timing information" log = 'screen19.minimum_exposure.log' .type = str .help = 'Location for the info log' @@ -180,7 +179,7 @@ def wilson_fit( def wilson_plot_ascii( - miller_array: miller.array, d_ticks: Optional[Sequence] = None + miller_array: miller.array, d_ticks: Sequence | None = None ) -> None: """ Print an ASCII-art Wilson plot of reflection intensities. @@ -221,8 +220,8 @@ def wilson_plot_image( d_star_sq: FloatSequence, intensity: FloatSequence, fit: Fit, - max_d: Optional[float] = None, - ticks: Optional[FloatSequence] = None, + max_d: float | None = None, + ticks: FloatSequence | None = None, output: str = "wilson_plot", ) -> None: """ diff --git a/src/screen/screen.py b/src/screen/screen.py index ffa479f..d2cd346 100644 --- a/src/screen/screen.py +++ b/src/screen/screen.py @@ -1,6 +1,43 @@ -"""The main screening script.""" +"""Pipeline to process screening data obtained at Diamond Light Source Beamline I19.""" + from __future__ import annotations +help_message = """ +This program presents the user with recommendations for adjustments to beam +flux, based on a single-sweep screening data collection. It presents an +upper- and lower-bound estimate of suitable flux.\n + • The upper-bound estimate is based on a comparison of a histogram of + measured pixel intensities with the trusted intensity range of the detector. + The user is warned when the measured pixel intensities indicate that the + detector would have a significant number of overloaded or untrustworthy + pixels.\n + • The lower-bound estimate is based on a linear fit of isotropic disorder + parameter, B, to a Wilson plot of reflection intensities. From this, an estimate + is made of the minimum exposure (flux × exposure time) required to achieve a target + I/σ ratio (by default, target I/σ = 2) at one or more values of desired resolution, d, + (by default, desired d = 1 Å, 0.84 Å, 0.6 Å & 0.4 Å).\n +\n +Target I/σ and target d (in Ångström) can be set using the parameters +'min_i_over_sigma' and 'desired_d'. One can set multiple values of the latter.\n +\n +By default the disorder parameter fit is conducted on the integrated data. +This ought to provide a reasonably true fit, but requires an integration step, +which can take some time. You can achieve a quicker, dirtier answer by fitting +to the indexed data (i.e. only the stronger spots), using 'minimum_exposure.data=indexed'.\n +\n +Examples:\n + + screen19 *.cbf\n + + screen19 /path/to/data/\n + + screen19 /path/to/data/image0001.cbf:1:100\n + + screen19 min_i_over_sigma=2 desired_d=0.84 \n + + screen19 minimum_exposure.data=indexed \n +""" + import argparse import re import subprocess @@ -30,7 +67,6 @@ """ verbosity = 0 .type = int(value_min=0) - .multiple = True .help = "Verbosity level of log output. Possible values:\n" "\t• 0: Info log output to stdout/logfile\n" "\t• 1: Info log output to stdout/logfile, logfile contains timing" @@ -96,8 +132,14 @@ def find_import_arguments( return in_value.as_posix(), None, None, None +usage = "%(prog)s [options] /path/to/data" + parser = argparse.ArgumentParser( - description=__doc__, parents=[version_parser, config_parser, options_parser] + usage=usage, + formatter_class=argparse.RawTextHelpFormatter, + description=__doc__, + epilog=help_message, + parents=[version_parser, config_parser, options_parser], ) parser.add_argument( "experiments", @@ -119,7 +161,7 @@ def find_import_arguments( ) -def run_import(images, params: ScopeExtract): +def run_import(images: list | str | None, params: ScopeExtract): # Ugly, but works import_params = import_scope.format(python_object=params) @@ -175,9 +217,10 @@ def run_integrate(params: ScopeExtract, options: list = []): ) -def run_minimum_exposure(choice): +def run_minimum_exposure(params, choice): if choice == "indexed": # subprocess.run + print(params.__dict__) pass else: subprocess.run( @@ -195,8 +238,9 @@ def pipeline(args: argparse.Namespace, working_phil: Scope): args.image_range if args.image_range else None ) - print(params.minimum_exposure.desired_d) # SIGH + # print(params.minimum_exposure.desired_d) # SIGH + # This probably makes the include scope in phil_scope a bit redundant... 2 choices to do this I guess? spot_finding_options = args.find_spots indexing_options = args.index refinement_options = args.refine @@ -209,9 +253,9 @@ def pipeline(args: argparse.Namespace, working_phil: Scope): if args.data == "integrated": run_refine(params.dials_refine, refinement_options) run_integrate(params.dials_integrate, integration_options) - run_minimum_exposure(args.data) + run_minimum_exposure(params.minimum_exposure, args.data) else: - run_minimum_exposure(args.data) + run_minimum_exposure(params.minimum_exposure, args.data) def main(args=None): From 1269ba73ef06509d9b4a73252dce1bf4122136b0 Mon Sep 17 00:00:00 2001 From: Noemi Frisina Date: Mon, 12 Dec 2022 17:34:26 +0000 Subject: [PATCH 30/33] Get minimum exposure to work correctly for indexing --- src/screen/minimum_exposure.py | 36 ++++++++++++++++++++++------------ src/screen/screen.py | 25 ++++++++++++++++++----- 2 files changed, 44 insertions(+), 17 deletions(-) diff --git a/src/screen/minimum_exposure.py b/src/screen/minimum_exposure.py index cbf1e2f..43a7603 100644 --- a/src/screen/minimum_exposure.py +++ b/src/screen/minimum_exposure.py @@ -75,6 +75,13 @@ minimum_exposure .caption = 'Parameters for the calculation of the lower exposure bound' { + data = indexed *integrated + .type = choice + .caption = '' + .help = 'Choice of data for the displacement parameter fit' + 'For the lower-bound exposure estimate, choose whether to use ' + 'indexed (quicker) or integrated (better) data in fitting ' + 'the isotropic displacement parameter.' desired_d = None .multiple = True .type = float @@ -300,19 +307,24 @@ def suggest_minimum_exposure( # Ignore all spots flagged as overloaded refls.del_selected(refls.get_flags(refls.flags.overloaded).iselection()) - # Work from profile-fitted intensities where possible but if the number of - # profile-fitted intensities is less than 75% of the number of summed - # intensities, use summed intensities instead. This is a very arbitrary heuristic. - sel_prf = refls.get_flags(refls.flags.integrated_prf).iselection() - sel_sum = refls.get_flags(refls.flags.integrated_sum).iselection() - if sel_prf.size() < 0.75 * sel_sum.size(): - refls = refls.select(sel_sum) + if params.minimum_exposure.data == "integrated": + # Work from profile-fitted intensities where possible but if the number of + # profile-fitted intensities is less than 75% of the number of summed + # intensities, use summed intensities instead. This is a very arbitrary heuristic. + sel_prf = refls.get_flags(refls.flags.integrated_prf).iselection() + sel_sum = refls.get_flags(refls.flags.integrated_sum).iselection() + if sel_prf.size() < 0.75 * sel_sum.size(): + refls = refls.select(sel_sum) + intensity = refls["intensity.sum.value"] + sigma = flex.sqrt(refls["intensity.sum.variance"]) + else: + refls = refls.select(sel_prf) + intensity = refls["intensity.prf.value"] + sigma = flex.sqrt(refls["intensity.prf.variance"]) + else: + # This is still a bit rough intensity = refls["intensity.sum.value"] sigma = flex.sqrt(refls["intensity.sum.variance"]) - else: - refls = refls.select(sel_prf) - intensity = refls["intensity.prf.value"] - sigma = flex.sqrt(refls["intensity.prf.variance"]) # Apply French-Wilson scaling to ensure positive intensities. miller_array = miller.array( @@ -451,7 +463,7 @@ def run(set_up_logging: bool = False): if set_up_logging: # Configure the logging - log.config(params.verbosity, params.output.log) + log.config(params.output.verbosity, params.output.log) expt = ExperimentList.from_file(args.experiments) refl = flex.reflection_table.from_file(args.reflections) diff --git a/src/screen/screen.py b/src/screen/screen.py index d2cd346..451eae7 100644 --- a/src/screen/screen.py +++ b/src/screen/screen.py @@ -55,6 +55,7 @@ options_parser, refine_scope, ) +from screen.minimum_exposure import phil_scope as minimum_exposure_scope # Custom types Scope = libtbx.phil.scope @@ -217,14 +218,28 @@ def run_integrate(params: ScopeExtract, options: list = []): ) -def run_minimum_exposure(params, choice): +def run_minimum_exposure(params: ScopeExtract, choice: str): + min_exp_params = minimum_exposure_scope.format(python_object=params) + if choice == "indexed": - # subprocess.run - print(params.__dict__) - pass + subprocess.run( + [ + "screen19.minimum_exposure", + "indexed.expt", + "indexed.refl", + f"minimum_exposure.data={choice}", + min_exp_params.as_str(), + ] + ) else: subprocess.run( - ["screen19.minimum_exposure", "integrated.expt", "integrated.refl"] + [ + "screen19.minimum_exposure", + "integrated.expt", + "integrated.refl", + f"minimum_exposure.data={choice}", + min_exp_params.as_str(), + ] ) From ead2b2de52376acb162f2084204727988da33816 Mon Sep 17 00:00:00 2001 From: Noemi Frisina Date: Mon, 12 Dec 2022 17:57:28 +0000 Subject: [PATCH 31/33] Parser workaround for minimum exposure --- src/screen/inputs.py | 7 +++++++ src/screen/minimum_exposure.py | 12 ++++++------ src/screen/screen.py | 22 ++++++++++------------ 3 files changed, 23 insertions(+), 18 deletions(-) diff --git a/src/screen/inputs.py b/src/screen/inputs.py index 7dd42e2..29bbae8 100644 --- a/src/screen/inputs.py +++ b/src/screen/inputs.py @@ -66,3 +66,10 @@ nargs="*", help="Additional options for dials.integrate.", ) +options_parser.add_argument( + "--min-exposure", + type=str, + default=[], + nargs="*", + help="Additional options for screen19.minimum_exposure.", +) diff --git a/src/screen/minimum_exposure.py b/src/screen/minimum_exposure.py index 43a7603..cf18e5d 100644 --- a/src/screen/minimum_exposure.py +++ b/src/screen/minimum_exposure.py @@ -438,9 +438,9 @@ def suggest_minimum_exposure( epilog=help_message, parents=[version_parser, config_parser], ) -parser.add_argument("experiments", type=str, help="") -parser.add_argument("reflections", type=str, help="") -parser.add_argument("phil_args", nargs="*", help="") +parser.add_argument("expt", type=str, help="Experiment file") +parser.add_argument("refl", type=str, help="Reflection file") +parser.add_argument("phil_args", nargs="*", help="Phil parameters for fit.") def run(set_up_logging: bool = False): @@ -465,12 +465,12 @@ def run(set_up_logging: bool = False): # Configure the logging log.config(params.output.verbosity, params.output.log) - expt = ExperimentList.from_file(args.experiments) - refl = flex.reflection_table.from_file(args.reflections) + expt = ExperimentList.from_file(args.expt) + refl = flex.reflection_table.from_file(args.refl) if len(expt) > 1: warn( - f"The experiment list you provided, {args.experiments}, contains more than one " + f"The experiment list you provided, {args.expt}, contains more than one " "experiment object (perhaps multiple indexing solutions). Only " "the first will be used, all others will be ignored." ) diff --git a/src/screen/screen.py b/src/screen/screen.py index 451eae7..f42b0ea 100644 --- a/src/screen/screen.py +++ b/src/screen/screen.py @@ -33,9 +33,9 @@ screen19 /path/to/data/image0001.cbf:1:100\n - screen19 min_i_over_sigma=2 desired_d=0.84 \n + screen19 -d indexed\n - screen19 minimum_exposure.data=indexed \n + screen19 --find-spots nproc=4 --index space_group=P222 """ import argparse @@ -55,7 +55,8 @@ options_parser, refine_scope, ) -from screen.minimum_exposure import phil_scope as minimum_exposure_scope + +# from screen.minimum_exposure import phil_scope as minimum_exposure_scope # Custom types Scope = libtbx.phil.scope @@ -218,9 +219,7 @@ def run_integrate(params: ScopeExtract, options: list = []): ) -def run_minimum_exposure(params: ScopeExtract, choice: str): - min_exp_params = minimum_exposure_scope.format(python_object=params) - +def run_minimum_exposure(choice: str, options: list = []): if choice == "indexed": subprocess.run( [ @@ -228,7 +227,7 @@ def run_minimum_exposure(params: ScopeExtract, choice: str): "indexed.expt", "indexed.refl", f"minimum_exposure.data={choice}", - min_exp_params.as_str(), + *options, ] ) else: @@ -238,7 +237,7 @@ def run_minimum_exposure(params: ScopeExtract, choice: str): "integrated.expt", "integrated.refl", f"minimum_exposure.data={choice}", - min_exp_params.as_str(), + *options, ] ) @@ -253,13 +252,12 @@ def pipeline(args: argparse.Namespace, working_phil: Scope): args.image_range if args.image_range else None ) - # print(params.minimum_exposure.desired_d) # SIGH - # This probably makes the include scope in phil_scope a bit redundant... 2 choices to do this I guess? spot_finding_options = args.find_spots indexing_options = args.index refinement_options = args.refine integration_options = args.integrate + min_exp_options = args.min_exposure run_import(args.experiments, params.dials_import) run_find_spots(params.dials_find_spots, spot_finding_options) @@ -268,9 +266,9 @@ def pipeline(args: argparse.Namespace, working_phil: Scope): if args.data == "integrated": run_refine(params.dials_refine, refinement_options) run_integrate(params.dials_integrate, integration_options) - run_minimum_exposure(params.minimum_exposure, args.data) + run_minimum_exposure(args.data, min_exp_options) else: - run_minimum_exposure(params.minimum_exposure, args.data) + run_minimum_exposure(args.data, min_exp_options) def main(args=None): From c9cd7140944b2864d7b64adf13a136711f7965aa Mon Sep 17 00:00:00 2001 From: Noemi Frisina Date: Thu, 15 Dec 2022 13:02:21 +0000 Subject: [PATCH 32/33] Tidy up help essage --- src/screen/minimum_exposure.py | 2 +- src/screen/screen.py | 4 +++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/screen/minimum_exposure.py b/src/screen/minimum_exposure.py index cf18e5d..c27f71c 100644 --- a/src/screen/minimum_exposure.py +++ b/src/screen/minimum_exposure.py @@ -31,7 +31,7 @@ screen19.minimum_exposure integrated.expt integrated.refl \n - screen19.minimum_exposure indexed.expt indexed.refl \n + screen19.minimum_exposure data=indexed indexed.expt indexed.refl \n screen19.minimum_exposure min_i_over_sigma=2 desired_d=0.84 wilson_fit_max_d=4 \ integrated.expt integrated.refl \n diff --git a/src/screen/screen.py b/src/screen/screen.py index f42b0ea..f41e7f5 100644 --- a/src/screen/screen.py +++ b/src/screen/screen.py @@ -36,6 +36,8 @@ screen19 -d indexed\n screen19 --find-spots nproc=4 --index space_group=P222 + + screen19 --min-exposure desired_d=0.84 min_i_over_sigma=2 """ import argparse @@ -134,7 +136,7 @@ def find_import_arguments( return in_value.as_posix(), None, None, None -usage = "%(prog)s [options] /path/to/data" +usage = "%(prog)s /path/to/data [options]" parser = argparse.ArgumentParser( usage=usage, From 00d9b9ffa867687267384daa9d37d82c51d64297 Mon Sep 17 00:00:00 2001 From: Noemi Frisina Date: Thu, 15 Dec 2022 18:54:47 +0000 Subject: [PATCH 33/33] Start setting up logging in screen script --- src/screen/screen.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/src/screen/screen.py b/src/screen/screen.py index f41e7f5..82fbbd5 100644 --- a/src/screen/screen.py +++ b/src/screen/screen.py @@ -41,6 +41,7 @@ """ import argparse +import logging import re import subprocess import sys @@ -48,6 +49,8 @@ import libtbx.phil +from dials.util import log + from screen import config_parser, version_parser # FIXME TODO change to relative import from screen.inputs import ( find_spots_scope, @@ -64,6 +67,11 @@ Scope = libtbx.phil.scope ScopeExtract = libtbx.phil.scope_extract +# Logging set up +logger_name = "dials.screen19" +logger = logging.getLogger(logger_name) +debug, info, warn = logger.debug, logger.info, logger.warning + template_pattern = re.compile(r"(.*)_(?:[0-9]*\#+).(.*)") @@ -247,6 +255,11 @@ def run_minimum_exposure(choice: str, options: list = []): def pipeline(args: argparse.Namespace, working_phil: Scope): params = working_phil.extract() + log.config(params.verbosity, logfile=params.output.log) + if not params.verbosity: + logging.getLogger("dials").setLevel(logging.WARNING) + logging.getLogger("dials.screen19").setLevel(logging.INFO) + # Set directory/template if that's what's been parsed. params.dials_import.input.directory = [args.directory] if args.directory else [] params.dials_import.input.template = [args.template] if args.template else []