Source code for dials.algorithms.refinement.refiner

"""Refiner is the refinement module public interface. RefinerFactory is
what should usually be used to construct a Refiner."""


import copy
import logging
import math

import psutil

import libtbx
from dxtbx.model.experiment_list import ExperimentList
from libtbx.phil import parse

import dials.util
from dials.algorithms.refinement import DialsRefineConfigError
from dials.algorithms.refinement.constraints import ConstraintManagerFactory
from dials.algorithms.refinement.engine import AdaptLstbx, refinery_phil_str
from dials.algorithms.refinement.parameterisation import (
    build_prediction_parameterisation,
)
from dials.algorithms.refinement.parameterisation import (
    phil_str as parameterisation_phil_str,
)
from dials.algorithms.refinement.parameterisation.autoreduce import AutoReduce
from dials.algorithms.refinement.parameterisation.parameter_report import (
    ParameterReporter,
)
from dials.algorithms.refinement.prediction.managed_predictors import (
    ExperimentsPredictorFactory,
)
from dials.algorithms.refinement.refinement_helpers import ordinal_number, string_sel
from dials.algorithms.refinement.reflection_manager import ReflectionManagerFactory
from dials.algorithms.refinement.reflection_manager import (
    phil_str as reflections_phil_str,
)
from dials.algorithms.refinement.restraints import RestraintsParameterisation
from dials.algorithms.refinement.target import TargetFactory
from dials.algorithms.refinement.target import phil_str as target_phil_str
from dials.array_family import flex

logger = logging.getLogger(__name__)

# The include scope directive does not work here. For example:
#
#   include scope dials.algorithms.refinement.outlier_detection.phil_scope
#
# results in:
#
#   AttributeError: 'module' object has no attribute 'refinement'
#
# to work around this, just include external phil scopes as strings
format_data = {
    "reflections_phil": reflections_phil_str,
    "target_phil": target_phil_str,
    "parameterisation_phil": parameterisation_phil_str,
    "refinery_phil": refinery_phil_str,
}
phil_scope = parse(
    """
refinement
  .help = "Parameters to configure the refinement"
{

  mp
    .expert_level = 2
  {
    nproc = 1
      .type = int(value_min=1)
      .help = "The number of processes to use. Not all choices of refinement"
              "engine support nproc > 1. Where multiprocessing is possible,"
              "it is helpful only in certain circumstances, so this is not"
              "recommended for typical use."
  }

  parameterisation
    .help = "Parameters to control the parameterisation of experimental models"
  {
    %(parameterisation_phil)s
  }

  %(refinery_phil)s

  target
    .help = "Parameters to configure the target function"
    .expert_level = 1
  {
    %(target_phil)s
  }

  reflections
    .help = "Parameters used by the reflection manager"
  {
    %(reflections_phil)s
  }

}
"""
    % format_data,
    process_includes=True,
)

RAD2DEG = 180 / math.pi


def _copy_experiments_for_refining(experiments):
    """
    Make a partial copy of experiments, copying properties used in refinement.

    Any experiment property that can be altered by refinement e.g. Beam,
    Detector and Goniometer will be deep-copied, whereas anything that the
    refiner doesn't touch (e.g. Scan, ImageSet) will be left as original
    references.

    This makes it safe to pass an object into the refiner or get an object
    out of the refiner without having to worry about your copy being
    unexpectedly altered, but saves time by avoiding the copying of potentially
    expensive experiment properties (e.g. ImageSet and its attributes).

    Args
      experiments (Experiment or ExperimentList or Iterable[Experiment]):

    Returns:
      ExperimentList: The copied experiments in new ExperimentList
    """
    # Look for a non-list e.g. a single experiment and convert to a list
    if not hasattr(experiments, "__iter__"):
        experiments = [experiments]
    out_list = ExperimentList()
    # Save a map of object id to copies so shared objects remain shared
    id_memo = {}

    # Copy each experiment individually
    for experiment in experiments:
        # Be inclusive about the initial copy
        new_exp = copy.copy(experiment)

        # Ensure every 'refined' attribute is uniquely copied
        for model in ["beam", "goniometer", "detector", "crystal"]:
            original = getattr(experiment, model)
            if id(original) not in id_memo:
                id_memo[id(original)] = copy.deepcopy(original)
            # assign the new copy to the experiment
            setattr(new_exp, model, id_memo[id(original)])

        # Collect this together
        out_list.append(new_exp)

    return out_list


def _trim_scans_to_observations(experiments, reflections):
    """Check the range of each scan matches the range of observed data and
    trim the scan to match if it is too wide"""

    # Get observed image number (or at least observed phi)
    obs_phi = reflections["xyzobs.mm.value"].parts()[2]
    try:
        obs_z = reflections["xyzobs.px.value"].parts()[2]
    except KeyError:
        obs_z = None

    # Get z_min and z_max from shoeboxes if present
    try:
        shoebox = reflections["shoebox"]
        bb = shoebox.bounding_boxes()
        z_min, z_max = bb.parts()[4:]
        if z_min.all_eq(0):
            shoebox = None
    except KeyError:
        shoebox = None

    for iexp, exp in enumerate(experiments):

        sel = reflections["id"] == iexp
        isel = sel.iselection()
        if obs_z is not None:
            exp_z = obs_z.select(isel)
        else:
            exp_phi = obs_phi.select(isel)
            exp_z = exp.scan.get_array_index_from_angle(exp_phi, deg=False)

        start, stop = exp.scan.get_array_range()
        min_exp_z = flex.min(exp_z)
        max_exp_z = flex.max(exp_z)

        # If observed array range is correct, skip to next experiment
        if int(min_exp_z) == start and int(math.ceil(max_exp_z)) == stop:
            continue

        # Extend array range either by shoebox size, or 0.5 deg if shoebox not available
        if shoebox is not None:
            obs_start = flex.min(z_min.select(isel))
            obs_stop = flex.max(z_max.select(isel))
        else:
            obs_start = int(min_exp_z)
            obs_stop = int(math.ceil(max_exp_z))
            half_deg_in_images = int(math.ceil(0.5 / exp.scan.get_oscillation()[1]))
            obs_start -= half_deg_in_images
            obs_stop += half_deg_in_images

        # Convert obs_start, obs_stop from position in array range to integer image number
        if obs_start > start or obs_stop < stop:
            im_start = max(start, obs_start) + 1
            im_stop = min(obs_stop, stop)

            logger.warning(
                "The reflections for experiment {0} do not fill the scan range. The scan will be trimmed "
                "to images {{{1},{2}}} to match the range of observed data".format(
                    iexp, im_start, im_stop
                )
            )

            # Ensure the scan is unique to this experiment and set trimmed limits
            exp.scan = copy.deepcopy(exp.scan)
            new_oscillation = (
                exp.scan.get_angle_from_image_index(im_start),
                exp.scan.get_oscillation()[1],
            )
            exp.scan.set_image_range((im_start, im_stop))
            exp.scan.set_oscillation(new_oscillation)

    return experiments


[docs]class RefinerFactory: """Factory class to create refiners""" @staticmethod def _filter_reflections(reflections): """Return a copy of the input reflections filtered to keep only those columns that are required by refinement""" cols = [ "id", "miller_index", "panel", "s1", "xyzobs.mm.value", "xyzobs.px.value", "xyzcal.px", "xyzobs.mm.variance", "flags", "shoebox", "delpsical.weights", ] # NB xyzobs.px.value & xyzcal.px required by SauterPoon outlier rejector # NB delpsical.weights is used by ExternalDelPsiWeightingStrategy rt = flex.reflection_table() # copy columns to the new table. Could use the select method # for this except that 's1' is optional in the input so would want # to copy that in like this if present anyway for k in cols: if k in reflections: rt[k] = reflections[k] return rt
[docs] @classmethod def from_parameters_data_experiments(cls, params, reflections, experiments): # TODO Checks on the input # E.g. does every experiment contain at least one overlapping model with at # least one other experiment? Are all the experiments either rotation series # or stills (the combination of both not yet supported)? # copy the experiments experiments = _copy_experiments_for_refining(experiments) # copy and filter the reflections reflections = cls._filter_reflections(reflections) return cls._build_components(params, reflections, experiments)
@classmethod def _build_components(cls, params, reflections, experiments): """low level build""" # Currently a refinement job can only have one parameterisation of the # prediction equation. This can either be of the XYDelPsi (stills) type, the # XYPhi (scans) type or the scan-varying XYPhi type with a varying crystal # model single_as_still = params.refinement.parameterisation.treat_single_image_as_still exps_are_stills = [] for exp in experiments: if exp.scan is None: exps_are_stills.append(True) elif exp.scan.get_num_images() == 1: if single_as_still: exps_are_stills.append(True) elif exp.scan.is_still(): exps_are_stills.append(True) else: exps_are_stills.append(False) else: if exp.scan.get_oscillation()[1] <= 0.0: raise DialsRefineConfigError("Cannot refine a zero-width scan") exps_are_stills.append(False) # check experiment types are consistent if not all(exps_are_stills[0] == e for e in exps_are_stills): raise DialsRefineConfigError("Cannot refine a mixture of stills and scans") do_stills = exps_are_stills[0] # If experiments are stills, ensure scan-varying refinement won't be attempted if do_stills: params.refinement.parameterisation.scan_varying = False # Refiner does not accept scan_varying=Auto. This is a special case for # doing macrocycles of refinement in dials.refine. if params.refinement.parameterisation.scan_varying is libtbx.Auto: params.refinement.parameterisation.scan_varying = False # Calculate reflection block_width for scan-varying refinement. Trim scans # to the extent of the observations, if requested. if params.refinement.parameterisation.scan_varying: if params.refinement.parameterisation.trim_scan_to_observations: experiments = _trim_scans_to_observations(experiments, reflections) from dials.algorithms.refinement.reflection_manager import BlockCalculator block_calculator = BlockCalculator(experiments, reflections) if params.refinement.parameterisation.compose_model_per == "block": reflections = block_calculator.per_width( params.refinement.parameterisation.block_width, deg=True ) elif params.refinement.parameterisation.compose_model_per == "image": reflections = block_calculator.per_image() logger.debug("\nBuilding reflection manager") logger.debug("Input reflection list size = %d observations", len(reflections)) # create reflection manager refman = ReflectionManagerFactory.from_parameters_reflections_experiments( params.refinement.reflections, reflections, experiments, do_stills ) logger.debug( "Number of observations that pass initial inclusion criteria = %d", refman.get_accepted_refs_size(), ) sample_size = refman.get_sample_size() if sample_size > 0: logger.debug("Working set size = %d observations", sample_size) logger.debug("Reflection manager built\n") # configure use of sparse data types params = cls.config_sparse(params, experiments) do_sparse = params.refinement.parameterisation.sparse # create managed reflection predictor ref_predictor = ExperimentsPredictorFactory.from_experiments( experiments, force_stills=do_stills, spherical_relp=params.refinement.parameterisation.spherical_relp_model, ) # Predict for the managed observations, set columns for residuals and set # the used_in_refinement flag to the predictions obs = refman.get_obs() ref_predictor(obs) x_obs, y_obs, phi_obs = obs["xyzobs.mm.value"].parts() x_calc, y_calc, phi_calc = obs["xyzcal.mm"].parts() obs["x_resid"] = x_calc - x_obs obs["y_resid"] = y_calc - y_obs obs["phi_resid"] = phi_calc - phi_obs # determine whether to do basic centroid analysis to automatically # determine outlier rejection block if params.refinement.reflections.outlier.block_width is libtbx.Auto: ca = refman.get_centroid_analyser() analysis = ca(calc_average_residuals=False, calc_periodograms=False) else: analysis = None # Now predictions and centroid analysis are available, so we can finalise # the reflection manager refman.finalise(analysis) # Create model parameterisations logger.debug("Building prediction equation parameterisation") pred_param = build_prediction_parameterisation( params.refinement.parameterisation, experiments, refman, do_stills ) # Build a constraints manager, if requested cmf = ConstraintManagerFactory(params, pred_param) constraints_manager = cmf() # Test for parameters that have too little data to refine and act accordingly autoreduce = AutoReduce( params.refinement.parameterisation.auto_reduction, pred_param, refman, constraints_manager, cmf, ) autoreduce() # if reduction was done, constraints_manager will have changed constraints_manager = autoreduce.constraints_manager # Build a restraints parameterisation (if requested). # Only unit cell restraints are supported at the moment. restraints_parameterisation = cls.config_restraints( params.refinement.parameterisation, pred_param ) # Parameter reporting logger.debug("Prediction equation parameterisation built") logger.debug("Parameter order : name mapping") for i, e in enumerate(pred_param.get_param_names()): logger.debug("Parameter %03d : %s", i + 1, e) param_reporter = ParameterReporter( pred_param.get_detector_parameterisations(), pred_param.get_beam_parameterisations(), pred_param.get_crystal_orientation_parameterisations(), pred_param.get_crystal_unit_cell_parameterisations(), pred_param.get_goniometer_parameterisations(), ) # Create target function logger.debug("Building target function") target = cls.config_target( params.refinement.target, experiments, refman, ref_predictor, pred_param, restraints_parameterisation, do_stills, do_sparse, ) logger.debug("Target function built") # create refinery logger.debug("Building refinement engine") refinery = cls.config_refinery(params, target, pred_param, constraints_manager) logger.debug("Refinement engine built") nparam = len(pred_param) ndim = target.dim nref = len(refman.get_matches()) logger.info( "There are %s parameters to refine against %s reflections in %s dimensions", nparam, nref, ndim, ) if not params.refinement.parameterisation.sparse and isinstance( refinery, AdaptLstbx ): dense_jacobian_gigabytes = ( nparam * nref * ndim * flex.double.element_size() ) / 1e9 avail_memory_gigabytes = psutil.virtual_memory().available / 1e9 # Report if the Jacobian requires a large amount of storage if ( dense_jacobian_gigabytes > 0.2 * avail_memory_gigabytes or dense_jacobian_gigabytes > 0.5 ): logger.info( "Storage of the Jacobian matrix requires %.1f GB", dense_jacobian_gigabytes, ) # build refiner interface and return if params.refinement.parameterisation.scan_varying: refiner = ScanVaryingRefiner else: refiner = Refiner return refiner( experiments, pred_param, param_reporter, refman, target, refinery )
[docs] @staticmethod def config_sparse(params, experiments): """Configure whether to use sparse datatypes""" # Automatic selection for sparse parameter if params.refinement.parameterisation.sparse == libtbx.Auto: if len(experiments) > 1: params.refinement.parameterisation.sparse = True else: params.refinement.parameterisation.sparse = False if params.refinement.refinery.engine == "SparseLevMar": params.refinement.parameterisation.sparse = True if params.refinement.mp.nproc > 1: if params.refinement.refinery.engine != "SparseLevMar": # sparse vectors cannot be pickled, so can't use easy_mp here params.refinement.parameterisation.sparse = False else: pass # but SparseLevMar requires sparse jacobian; does not implement mp # Check incompatible selection elif ( params.refinement.parameterisation.sparse and params.refinement.mp.nproc > 1 ): logger.warning( "Could not set sparse=True and nproc=%s", params.refinement.mp.nproc ) logger.warning("Resetting sparse=False") params.refinement.parameterisation.sparse = False return params
[docs] @staticmethod def config_restraints(params, pred_param): """Given a set of user parameters plus a model parameterisation, create restraints plus a parameterisation of these restraints Params: params: The input PHIL parameters pred_param: A PredictionParameters object Returns: A restraints parameterisation or None """ if not any( ( params.crystal.unit_cell.restraints.tie_to_target, params.crystal.unit_cell.restraints.tie_to_group, ) ): return None if params.scan_varying and not params.crystal.unit_cell.force_static: logger.warning("Restraints will be ignored for scan_varying=True") return None det_params = pred_param.get_detector_parameterisations() beam_params = pred_param.get_beam_parameterisations() xl_ori_params = pred_param.get_crystal_orientation_parameterisations() xl_uc_params = pred_param.get_crystal_unit_cell_parameterisations() gon_params = pred_param.get_goniometer_parameterisations() rp = RestraintsParameterisation( detector_parameterisations=det_params, beam_parameterisations=beam_params, xl_orientation_parameterisations=xl_ori_params, xl_unit_cell_parameterisations=xl_uc_params, goniometer_parameterisations=gon_params, ) # Shorten params path cell_r = params.crystal.unit_cell.restraints for tie in cell_r.tie_to_target: if len(tie.values) != 6: raise DialsRefineConfigError( "6 cell parameters must be provided as the tie_to_target.values." ) if len(tie.sigmas) != 6: raise DialsRefineConfigError( "6 sigmas must be provided as the tie_to_target.sigmas. " "Note that individual sigmas of 0.0 will remove " "the restraint for the corresponding cell parameter." ) if tie.id is None: # get one experiment id for each parameterisation to apply to all tie.id = [e.get_experiment_ids()[0] for e in xl_uc_params] for exp_id in tie.id: rp.add_restraints_to_target_xl_unit_cell(exp_id, tie.values, tie.sigmas) for tie in cell_r.tie_to_group: if len(tie.sigmas) != 6: raise DialsRefineConfigError( "6 sigmas must be provided as the tie_to_group.sigmas. " "Note that individual sigmas of 0.0 will remove " "the restraint for the corresponding cell parameter." ) if tie.id is None: rp.add_restraints_to_group_xl_unit_cell(tie.target, "all", tie.sigmas) else: rp.add_restraints_to_group_xl_unit_cell(tie.target, tie.id, tie.sigmas) return rp
[docs] @staticmethod def config_refinery(params, target, pred_param, constraints_manager): """Given a set of parameters, a target class, a prediction parameterisation class and a constraints_manager (which could be None), build a refinery Params: params The input parameters Returns: The refinery instance """ # Shorten parameter path options = params.refinement.refinery if options.engine == "SimpleLBFGS": from dials.algorithms.refinement.engine import SimpleLBFGS as refinery elif options.engine == "LBFGScurvs": from dials.algorithms.refinement.engine import LBFGScurvs as refinery elif options.engine == "GaussNewton": from dials.algorithms.refinement.engine import ( GaussNewtonIterations as refinery, ) elif options.engine == "LevMar": from dials.algorithms.refinement.engine import ( LevenbergMarquardtIterations as refinery, ) elif options.engine == "SparseLevMar": from dials.algorithms.refinement.sparse_engine import ( SparseLevenbergMarquardtIterations as refinery, ) else: raise RuntimeError( "Refinement engine " + options.engine + " not recognised" ) logger.debug("Selected refinement engine type: %s", options.engine) engine = refinery( target=target, prediction_parameterisation=pred_param, constraints_manager=constraints_manager, log=options.log, tracking=options.journal, max_iterations=options.max_iterations, ) if params.refinement.mp.nproc > 1: nproc = params.refinement.mp.nproc try: engine.set_nproc(nproc) except NotImplementedError: logger.warning( "Could not set nproc=%s for refinement engine of type %s", nproc, options.engine, ) return engine
# Overload to allow subclasses of RefinerFactory to use a different # TargetFactory
[docs] @staticmethod def config_target( params, experiments, reflection_manager, predictor, pred_param, restraints_param, do_stills, do_sparse, ): target = TargetFactory.from_parameters_and_experiments( params, experiments, reflection_manager, predictor, pred_param, restraints_param, do_stills, do_sparse, ) return target
[docs]class Refiner: """Public interface for performing DIALS refinement. Public methods: run rmsds get_experiments get_matches get_param_reporter parameter_correlation_plot selection_used_for_refinement predict_for_reflection_table predict_for_indexed Notes: * The return value of run is a recorded history of the refinement * The experiments accessor provides a copy of the experiments used by refinement * get_matches exposes the function of the same name from the privately stored reflection manager * The return value of selection_used_for_refinement is a flex.bool """
[docs] def __init__( self, experiments, pred_param, param_reporter, refman, target, refinery ): """ Mandatory arguments: experiments - a dxtbx ExperimentList object pred_param - An object derived from the PredictionParameterisation class param_reporter -A ParameterReporter object refman - A ReflectionManager object target - An object derived from the Target class refinery - An object derived from the Refinery class """ # the experimental models self._experiments = experiments # refinement module main objects self._pred_param = pred_param self._refman = refman self._target = target self._refinery = refinery # parameter reporter self._param_report = param_reporter # Keep track of whether this is stills or scans type refinement self.experiment_type = refman.experiment_type return
[docs] def get_experiments(self): """Return a copy of the current refiner experiments""" return _copy_experiments_for_refining(self._experiments)
[docs] def rmsds(self): """Return rmsds of the current model""" # ensure predictions for the matches are up to date self._refinery.prepare_for_step() return self._target.rmsds()
[docs] def rmsds_for_reflection_table(self, reflections): """Calculate unweighted RMSDs for the specified reflections""" # ensure predictions for these reflections are up to date preds = self.predict_for_reflection_table(reflections) return self._target.rmsds_for_reflection_table(preds)
[docs] def get_matches(self): """Delegated to the reflection manager""" return self._refman.get_matches()
[docs] def get_free_reflections(self): """Delegated to the reflection manager""" return self._refman.get_free_reflections()
[docs] def get_param_reporter(self): """Get the ParameterReport object linked to this Refiner""" return self._param_report
[docs] def get_parameter_correlation_matrix(self, step, col_select=None): """Return the correlation matrix between columns of the Jacobian at the specified refinement step. The parameter col_select can be used to select subsets of the full number of columns. The column labels are also returned as a list of strings""" corrmats = self._refinery.get_correlation_matrix_for_step(step) if corrmats is None: return None, None all_labels = self._pred_param.get_param_names() if col_select is None: col_select = list(range(len(all_labels))) sel = string_sel(col_select, all_labels) labels = [e for e, s in zip(all_labels, sel) if s] num_cols = len(labels) if num_cols == 0: return None, None for k, corrmat in corrmats.items(): assert corrmat.is_square_matrix() idx = flex.bool(sel).iselection() sub_corrmat = flex.double(flex.grid(num_cols, num_cols)) for (i, x) in enumerate(idx): for (j, y) in enumerate(idx): sub_corrmat[i, j] = corrmat[x, y] corrmats[k] = sub_corrmat return (corrmats, labels)
@property def history(self): """Get the refinement engine's step history""" return self._refinery.history
[docs] def print_step_table(self): """print useful output about refinement steps in the form of a simple table""" logger.info("\nRefinement steps:") rmsd_multipliers = [] header = ["Step", "Nref"] for (name, units) in zip(self._target.rmsd_names, self._target.rmsd_units): if units == "mm": header.append(name + "\n(mm)") rmsd_multipliers.append(1.0) elif units == "rad": # convert radians to degrees for reporting header.append(name + "\n(deg)") rmsd_multipliers.append(RAD2DEG) else: # leave unknown units alone header.append(name + "\n(" + units + ")") rows = [] for i in range(self._refinery.history.get_nrows()): rmsds = [ r * m for (r, m) in zip(self._refinery.history["rmsd"][i], rmsd_multipliers) ] rows.append( [str(i), str(self._refinery.history["num_reflections"][i])] + [f"{r:.5g}" for r in rmsds] ) logger.info(dials.util.tabulate(rows, header)) logger.info(self._refinery.history.reason_for_termination) return
[docs] def print_out_of_sample_rmsd_table(self): """print out-of-sample RSMDs per step, if these were tracked""" # check if it makes sense to proceed if "out_of_sample_rmsd" not in self._refinery.history: return nref = len(self.get_free_reflections()) if nref < 10: return # don't do anything if very few refs logger.info("\nRMSDs for out-of-sample (free) reflections:") rmsd_multipliers = [] header = ["Step", "Nref"] for (name, units) in zip(self._target.rmsd_names, self._target.rmsd_units): if units == "mm": header.append(name + "\n(mm)") rmsd_multipliers.append(1.0) elif units == "rad": # convert radians to degrees for reporting header.append(name + "\n(deg)") rmsd_multipliers.append(RAD2DEG) else: # leave unknown units alone header.append(name + "\n(" + units + ")") rows = [] for i in range(self._refinery.history.get_nrows()): rmsds = [ r * m for r, m in zip( self._refinery.history["out_of_sample_rmsd"][i], rmsd_multipliers ) ] rows.append([str(i), str(nref)] + [f"{e:.5g}" for e in rmsds]) logger.info(dials.util.tabulate(rows, header)) return
[docs] def print_exp_rmsd_table(self): """print useful output about refinement steps in the form of a simple table""" logger.info("\nRMSDs by experiment:") header = ["Exp\nid", "Nref"] for (name, units) in zip(self._target.rmsd_names, self._target.rmsd_units): if name == "RMSD_X" or name == "RMSD_Y" and units == "mm": header.append(name + "\n(px)") elif name == "RMSD_Phi" and units == "rad": # will convert radians to images for reporting of scans header.append("RMSD_Z" + "\n(images)") elif units == "rad": # will convert other angles in radians to degrees (e.g. for # RMSD_DeltaPsi and RMSD_2theta) header.append(name + "\n(deg)") else: # skip other/unknown RMSDs pass rows = [] for iexp, exp in enumerate(self._experiments): detector = exp.detector px_sizes = [p.get_pixel_size() for p in detector] it = iter(px_sizes) px_size = next(it) if not all(tst == px_size for tst in it): logger.info( "The detector in experiment %d does not have the same pixel " + "sizes on each panel. Skipping...", iexp, ) continue px_per_mm = [1.0 / e for e in px_size] scan = exp.scan try: images_per_rad = 1.0 / abs(scan.get_oscillation(deg=False)[1]) except (AttributeError, ZeroDivisionError): images_per_rad = None raw_rmsds = self._target.rmsds_for_experiment(iexp) if raw_rmsds is None: continue # skip experiments where rmsd cannot be calculated num = self._target.get_num_matches_for_experiment(iexp) rmsds = [] for (name, units, rmsd) in zip( self._target.rmsd_names, self._target.rmsd_units, raw_rmsds ): if name == "RMSD_X" and units == "mm": rmsds.append(rmsd * px_per_mm[0]) elif name == "RMSD_Y" and units == "mm": rmsds.append(rmsd * px_per_mm[1]) elif name == "RMSD_Phi" and units == "rad": rmsds.append(rmsd * images_per_rad) elif units == "rad": rmsds.append(rmsd * RAD2DEG) rows.append([str(iexp), str(num)] + [f"{r:.5g}" for r in rmsds]) if len(rows) > 0: logger.info(dials.util.tabulate(rows, header)) return
[docs] def print_panel_rmsd_table(self): """print useful output about refinement steps in the form of a simple table""" if len(self._experiments.scans()) > 1: logger.warning( "Multiple scans present. Only the first scan will be used " "to determine the image width for reporting RMSDs" ) scan = self._experiments.scans()[0] try: images_per_rad = 1.0 / abs(scan.get_oscillation(deg=False)[1]) except AttributeError: images_per_rad = None for idetector, detector in enumerate(self._experiments.detectors()): if len(detector) == 1: continue logger.info("\nDetector %s RMSDs by panel:", idetector + 1) header = ["Panel\nid", "Nref"] for (name, units) in zip(self._target.rmsd_names, self._target.rmsd_units): if name == "RMSD_X" or name == "RMSD_Y" and units == "mm": header.append(name + "\n(px)") elif ( name == "RMSD_Phi" and units == "rad" ): # convert radians to images for reporting of scans header.append("RMSD_Z" + "\n(images)") elif ( name == "RMSD_DeltaPsi" and units == "rad" ): # convert radians to degrees for reporting of stills header.append(name + "\n(deg)") else: # skip RMSDs that cannot be expressed in image/scan space pass rows = [] for ipanel, panel in enumerate(detector): px_size = panel.get_pixel_size() px_per_mm = [1.0 / e for e in px_size] num = self._target.get_num_matches_for_panel(ipanel) if num <= 0: continue raw_rmsds = self._target.rmsds_for_panel(ipanel) if raw_rmsds is None: continue # skip panels where rmsd cannot be calculated rmsds = [] for (name, units, rmsd) in zip( self._target.rmsd_names, self._target.rmsd_units, raw_rmsds ): if name == "RMSD_X" and units == "mm": rmsds.append(rmsd * px_per_mm[0]) elif name == "RMSD_Y" and units == "mm": rmsds.append(rmsd * px_per_mm[1]) elif name == "RMSD_Phi" and units == "rad": rmsds.append(rmsd * images_per_rad) elif name == "RMSD_DeltaPsi" and units == "rad": rmsds.append(rmsd * RAD2DEG) rows.append([str(ipanel), str(num)] + [f"{r:.5g}" for r in rmsds]) if len(rows) > 0: logger.info(dials.util.tabulate(rows, header)) return
[docs] def run(self): """Run refinement""" #################################### # Do refinement and return history # #################################### logger.debug("\nExperimental models before refinement:") for i, beam in enumerate(self._experiments.beams()): logger.debug(ordinal_number(i) + " " + str(beam)) for i, detector in enumerate(self._experiments.detectors()): logger.debug(ordinal_number(i) + " " + str(detector)) for i, goniometer in enumerate(self._experiments.goniometers()): if goniometer is None: continue logger.debug(ordinal_number(i) + " " + str(goniometer)) for i, scan in enumerate(self._experiments.scans()): if scan is None: continue logger.debug(ordinal_number(i) + " " + str(scan)) for i, crystal in enumerate(self._experiments.crystals()): logger.debug(ordinal_number(i) + " " + str(crystal)) self._refinery.run() # These involve calculation, so skip them when output is quiet if logger.getEffectiveLevel() < logging.ERROR: self.print_step_table() self.print_out_of_sample_rmsd_table() self.print_exp_rmsd_table() det_npanels = [len(d) for d in self._experiments.detectors()] if any(n > 1 for n in det_npanels): self.print_panel_rmsd_table() # Perform post-run tasks to write the refined states back to the models self._update_models() logger.debug("\nExperimental models after refinement:") for i, beam in enumerate(self._experiments.beams()): logger.debug(ordinal_number(i) + " " + str(beam)) for i, detector in enumerate(self._experiments.detectors()): logger.debug(ordinal_number(i) + " " + str(detector)) for i, goniometer in enumerate(self._experiments.goniometers()): if goniometer is None: continue logger.debug(ordinal_number(i) + " " + str(goniometer)) for i, scan in enumerate(self._experiments.scans()): if scan is None: continue logger.debug(ordinal_number(i) + " " + str(scan)) for i, crystal in enumerate(self._experiments.crystals()): logger.debug(ordinal_number(i) + " " + str(crystal)) # Report on the refined parameters logger.debug(str(self._param_report)) # Return the refinement history return self._refinery.history
def _update_models(self): """Perform any extra tasks required to update the models after refinement. Does nothing here, but used by subclasses""" pass
[docs] def selection_used_for_refinement(self): """Return a selection as a flex.bool in terms of the input reflection data of those reflections that were used in the final step of refinement.""" matches = self._refman.get_matches() selection = flex.bool(len(self._refman.get_indexed()), False) try: # new reflection table format for matches isel = matches["iobs"] selection.set_selected(isel, True) except TypeError: # old ObsPredMatch format for matches for m in matches: selection[m.iobs] = True return selection
[docs] def predict_for_indexed(self): """perform prediction for all the indexed reflections passed into refinement and additionally set the used_in_refinement flag. Do not compose the derivatives of states of the model as this is expensive and they are not needed outside of a refinement run""" reflections = self.predict_for_reflection_table( self._refman.get_indexed(), skip_derivatives=True ) reflections.sort("iobs") mask = self.selection_used_for_refinement() reflections.set_flags(mask, reflections.flags.used_in_refinement) return reflections
[docs] def predict_for_reflection_table(self, reflections, skip_derivatives=False): """perform prediction for all reflections in the supplied table""" # delegate to the target object, which has access to the predictor return self._target.predict_for_reflection_table(reflections, skip_derivatives)
[docs]class ScanVaryingRefiner(Refiner): """Includes functionality to update the models with their states at scan-points after scan-varying refinement""" def _update_models(self): for iexp, exp in enumerate(self._experiments): ar_range = exp.scan.get_array_range() obs_image_numbers = list(range(ar_range[0], ar_range[1] + 1)) # write scan-varying s0 vectors back to beam models s0_list = self._pred_param.get_varying_s0(obs_image_numbers, iexp) if s0_list is not None: exp.beam.set_s0_at_scan_points(s0_list) # write scan-varying setting rotation matrices back to goniometer models S_list = self._pred_param.get_varying_setting_rotation( obs_image_numbers, iexp ) if S_list is not None: exp.goniometer.set_setting_rotation_at_scan_points(S_list) # write scan-varying crystal setting matrices back to crystal models A_list = self._pred_param.get_varying_UB(obs_image_numbers, iexp) if A_list is not None: exp.crystal.set_A_at_scan_points(A_list) # Calculate scan-varying errors if requested if self._pred_param.set_scan_varying_errors: # get state covariance matrices the whole range of images. We select # the first element of this at each image because crystal scan-varying # parameterisations are not multi-state state_cov_list = [ self._pred_param.calculate_model_state_uncertainties( obs_image_number=t, experiment_id=iexp ) for t in range(ar_range[0], ar_range[1] + 1) ] if "U_cov" in state_cov_list[0]: u_cov_list = [e["U_cov"] for e in state_cov_list] else: u_cov_list = None if "B_cov" in state_cov_list[0]: b_cov_list = [e["B_cov"] for e in state_cov_list] else: b_cov_list = None # return these to the model parameterisations to be set in the models self._pred_param.set_model_state_uncertainties( u_cov_list, b_cov_list, iexp )