diff --git a/.azure-pipelines/azure-pipelines.yml b/.azure-pipelines/azure-pipelines.yml new file mode 100644 index 0000000..d2f048d --- /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-screen19 + + - bash: | + 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/.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..e671074 --- /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 + 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=screen19 --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)") diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 2b49cbd..c483a09 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 @@ -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/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/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 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/setup.cfg b/setup.cfg index 71e9849..83659df 100644 --- a/setup.cfg +++ b/setup.cfg @@ -3,15 +3,15 @@ current_version = 0.213 commit = True tag = True parse = (?P\d+)\.(?P\d+)(\.(?P\d+))? -serialize = +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}" @@ -22,23 +22,51 @@ exclude = docs test = pytest [metadata] -classifiers = +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 +long_description = file: README.rst +licence = BSD license +licence_file = LICENSE +classifiers = Development Status :: 4 - Beta 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 +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 +zip_safe = False + +[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 = screen19 + screen19.minimum_exposure = screen19.minimum_exposure +libtbx.precommit = + screen19 = screen19 [tool:pytest] collect_ignore = ['setup.py'] diff --git a/setup.py b/setup.py index d37e415..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,36 +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", - 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", - ], - "libtbx.precommit": ["screen19 = screen19"], - }, - install_requires=['typing;python_version<"3.5"', "procrunner"], - license="BSD license", - long_description="\n\n".join([readme, changelog_header, changelog]), - include_package_data=True, - name="screen19", - packages=find_packages(), - test_suite="tests", - tests_require=["mock>=2.0", "pytest>=4.5"], - version="0.213", - zip_safe=False, -) +setup() diff --git a/screen19/__init__.py b/src/screen/__init__.py similarity index 58% rename from screen19/__init__.py rename to src/screen/__init__.py index 2d957fc..ebef663 100644 --- a/screen19/__init__.py +++ b/src/screen/__init__.py @@ -1,48 +1,70 @@ -"""Common tools for the I19 module.""" +"""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(".")) + + +import argparse 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" +from dials.util import version 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. +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] -def terminal_size() -> Tuple[int, int]: - """ - Find the current size of the terminal window. +def terminal_size() -> tuple[int, int]: + """Find the current size of the terminal window. - :return: Number of columns; number of rows. + :return: Tuple[int, int]: Number of columns; number of rows """ - columns, rows = 80, 25 + 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, - ) + 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: # ignore any errors and use default size - pass # FIXME: Can we be more specific about the type of exception? + except Exception as err: + logger.exception(err) + logger.warning("Exception caught and ignored, using default size.") + pass columns = min(columns, 120) - rows = min(rows, int(columns / 3)) + rows = min(rows, columns // 3) return columns, rows @@ -68,9 +90,9 @@ def prettyprint_dictionary(d): ) -def prettyprint_procrunner(d): +def prettyprint_process(d): """ - Produce a nice string representation of a procrunner ReturnObject, for printing. + Produce a nice string representation of a subprocess CompletedProcess object, for printing. :param d: ReturnObject to be printed. :return: String representation of :param d:. @@ -86,37 +108,6 @@ def prettyprint_procrunner(d): ) -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, @@ -158,14 +149,12 @@ def plot_intensities( debug("running %s with:\n %s\n", " ".join(command), "\n ".join(plot_commands)) try: - result = procrunner.run( + result = subprocess.run( command, - stdin="\n".join(plot_commands).encode("utf-8") + b"\n", + input="\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": ""}, + capture_output=True, + env={"LD_LIBRARY_PATH": ""}, ) except (OSError, subprocess.TimeoutExpired): info(traceback.format_exc()) @@ -175,7 +164,7 @@ def plot_intensities( ) return else: - debug("result = %s", prettyprint_procrunner(result)) + debug("result = %s", prettyprint_process(result)) returncode = getattr(result, "returncode") if returncode: diff --git a/src/screen/inputs.py b/src/screen/inputs.py new file mode 100644 index 0000000..29bbae8 --- /dev/null +++ b/src/screen/inputs.py @@ -0,0 +1,75 @@ +import argparse + +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, +) + +options_parser = argparse.ArgumentParser(add_help=False) +options_parser.add_argument( + "--find-spots", + type=str, + default=[], + nargs="*", + help="Additional options for dials.find_spots.", +) +options_parser.add_argument( + "--index", + type=str, + default=[], + nargs="*", + help="Additional options for dials.index.", +) +options_parser.add_argument( + "--refine", + type=str, + default=[], + nargs="*", + help="Additional options for dials.refine.", +) +options_parser.add_argument( + "--integrate", + type=str, + default=[], + 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/screen19/minimum_exposure.py b/src/screen/minimum_exposure.py similarity index 66% rename from screen19/minimum_exposure.py rename to src/screen/minimum_exposure.py index 406d794..c27f71c 100644 --- a/screen19/minimum_exposure.py +++ b/src/screen/minimum_exposure.py @@ -1,114 +1,124 @@ """ 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 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 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 + 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, Sequence, Union import numpy as np 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.version import dials_version from dxtbx.model import ExperimentList -from screen19 import __version__, d_ticks, plot_intensities, terminal_size +from screen import ( # FIXME TODO change to relative import + config_parser, + d_ticks, + plot_intensities, + terminal_size, + version_parser, +) # 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) - .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." - } - """, + .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 + .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 { + 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' + 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, ) @@ -176,7 +186,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. @@ -217,8 +227,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: """ @@ -272,7 +282,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. @@ -297,24 +307,29 @@ 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( miller.set( - expts[0].crystal.get_crystal_symmetry(), + expts.crystal.get_crystal_symmetry(), refls["miller_index"], anomalous_flag=False, ), @@ -415,77 +430,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("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.") + - 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 = OptionParser( - 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)) - params, options = parser.parse_args(args=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 = 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, - ) + log.config(params.output.verbosity, params.output.log) - expts = params.input.experiments[0].data - refls = params.input.reflections[0].data + expt = ExperimentList.from_file(args.expt) + refl = flex.reflection_table.from_file(args.refl) - 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.expt}, 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 new file mode 100644 index 0000000..82fbbd5 --- /dev/null +++ b/src/screen/screen.py @@ -0,0 +1,303 @@ +"""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 -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 +import logging +import re +import subprocess +import sys +from pathlib import Path + +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, + import_scope, + index_scope, + integrate_scope, + options_parser, + refine_scope, +) + +# from screen.minimum_exposure import phil_scope as minimum_exposure_scope + +# Custom types +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]*\#+).(.*)") + + +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" + output { + log = screen19.log + .type = str + .help = The log file name. + } + dials_import { + include scope screen.inputs.import_scope + } + dials_find_spots{ + include scope screen.inputs.find_spots_scope + } + dials_index { + include scope screen.inputs.index_scope + } + dials_refine { + include scope screen.inputs.refine_scope + } + dials_integrate { + include scope screen.inputs.integrate_scope + } + include scope screen.minimum_exposure.phil_scope + """, + process_includes=True, +) + + +class _ImportImages(argparse.Action): + def __call__(self, parser, namespace, values, option_string=None): + 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, 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, None + elif match: + 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, None + + +usage = "%(prog)s /path/to/data [options]" + +parser = argparse.ArgumentParser( + usage=usage, + formatter_class=argparse.RawTextHelpFormatter, + description=__doc__, + epilog=help_message, + parents=[version_parser, config_parser, options_parser], +) +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: list | str | None, params: ScopeExtract): + # 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()]) + + +def run_find_spots(params: ScopeExtract, options: list = []): + find_spots_params = find_spots_scope.format(python_object=params) + + subprocess.run( + ["dials.find_spots", "imported.expt", find_spots_params.as_str(), *options] + ) + + +def run_indexing(params: ScopeExtract, options: list = []): + index_params = index_scope.format(python_object=params) + + subprocess.run( + ["dials.index", "imported.expt", "strong.refl", index_params.as_str(), *options] + ) + + +def run_refine(params: ScopeExtract, options: list = []): + refine_params = refine_scope.format(python_object=params) + + subprocess.run( + [ + "dials.refine", + "indexed.expt", + "indexed.refl", + refine_params.as_str(), + *options, + ] + ) + + +def run_integrate(params: ScopeExtract, options: list = []): + integrate_params = integrate_scope.format(python_object=params) + + subprocess.run( + [ + "dials.integrate", + "refined.expt", + "refined.refl", + integrate_params.as_str(), + *options, + ] + ) + + +def run_minimum_exposure(choice: str, options: list = []): + if choice == "indexed": + subprocess.run( + [ + "screen19.minimum_exposure", + "indexed.expt", + "indexed.refl", + f"minimum_exposure.data={choice}", + *options, + ] + ) + else: + subprocess.run( + [ + "screen19.minimum_exposure", + "integrated.expt", + "integrated.refl", + f"minimum_exposure.data={choice}", + *options, + ] + ) + + +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 [] + params.dials_import.geometry.scan.image_range = ( + args.image_range if args.image_range else None + ) + + # 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) + run_indexing(params.dials_index, indexing_options) + subprocess.run(["dev.dials.pixel_histogram", "indexed.refl"]) + if args.data == "integrated": + run_refine(params.dials_refine, refinement_options) + run_integrate(params.dials_integrate, integration_options) + run_minimum_exposure(args.data, min_exp_options) + else: + run_minimum_exposure(args.data, min_exp_options) + + +def main(args=None): + args = parser.parse_args(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() + + pipeline(args, working_phil) + + +# if __name__ == "__main__": +# main(None) diff --git a/tests/test_screen.py b/tests/test_screen.py new file mode 100644 index 0000000..c700dcc --- /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=[None]) + + +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 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