diff --git a/src/dials/algorithms/indexing/indexer.py b/src/dials/algorithms/indexing/indexer.py index e401213020..0b91506b1e 100644 --- a/src/dials/algorithms/indexing/indexer.py +++ b/src/dials/algorithms/indexing/indexer.py @@ -930,14 +930,24 @@ def show_experiments( if unindexed_reflections: sel = unindexed_reflections["imageset_id"] == i unindexed_count += sel.count(True) - rows.append( - [ - str(i), - str(indexed_count), - str(unindexed_count), - f"{indexed_count / (indexed_count + unindexed_count)*100:.1f}", - ] - ) + if indexed_count + unindexed_count > 0: + rows.append( + [ + str(i), + str(indexed_count), + str(unindexed_count), + f"{indexed_count / (indexed_count + unindexed_count)*100:.1f}", + ] + ) + else: + rows.append( + [ + str(i), + str(indexed_count), + str(unindexed_count), + f"{0:.1f}", + ] + ) logger.info(dials.util.tabulate(rows, headers="firstrow")) def find_max_cell(self): diff --git a/src/dials/algorithms/refinement/prediction/managed_predictors.py b/src/dials/algorithms/refinement/prediction/managed_predictors.py index a4f5ca6e35..1aab182f80 100644 --- a/src/dials/algorithms/refinement/prediction/managed_predictors.py +++ b/src/dials/algorithms/refinement/prediction/managed_predictors.py @@ -85,6 +85,8 @@ def __call__(self, reflections): # select the reflections for this experiment only sel = reflections["id"] == iexp refs = reflections.select(sel) + if len(refs) == 0: + continue self._predict_one_experiment(e, refs) refs = self._post_predict_one_experiment(e, refs) @@ -175,6 +177,8 @@ def _predict_one_experiment(self, experiment, reflections): class LaueExperimentsPredictor(ExperimentsPredictor): def _predict_one_experiment(self, experiment, reflections): + if len(reflections) == 0: + return min_s0_idx = min( range(len(reflections["wavelength"])), @@ -197,6 +201,9 @@ def _predict_one_experiment(self, experiment, reflections): class TOFExperimentsPredictor(LaueExperimentsPredictor): def _post_predict_one_experiment(self, experiment, reflections): + if len(reflections) == 0: + return + # Add ToF to xyzcal.mm wavelength_cal = reflections["wavelength_cal"] distance = experiment.beam.get_sample_to_source_distance() * 10**-3 diff --git a/src/dials/command_line/tof_integrate.py b/src/dials/command_line/tof_integrate.py index 364f122057..be3fca3eea 100644 --- a/src/dials/command_line/tof_integrate.py +++ b/src/dials/command_line/tof_integrate.py @@ -9,6 +9,7 @@ import cctbx.array_family.flex import libtbx +from dxtbx import flumpy import dials.util.log from dials.algorithms.integration.fit.tof_line_profile import ( @@ -61,6 +62,17 @@ .type = str .help = "The log filename" } +integration_type = *observed calculated + .type = choice + .help = "observed integrates only reflections observed during spot finding" + "calculated integrates all reflections out to calculated.dmin" + +calculated{ + dmin = none + .type = float (value_min=0.5) + .help = "The resolution spots are integrated to when using integration_type.calculated" + +} method{ line_profile_fitting = False .type = bool @@ -293,6 +305,182 @@ def join_reflections(list_of_reflections): return reflections +def get_predicted_observed_reflections(params, experiments, reflections): + min_s0_idx = min( + range(len(reflections["wavelength"])), key=reflections["wavelength"].__getitem__ + ) + min_s0 = reflections["s0"][min_s0_idx] + dmin = None + for experiment in experiments: + expt_dmin = experiment.detector.get_max_resolution(min_s0) + if dmin is None or expt_dmin < dmin: + dmin = expt_dmin + + predicted_reflections = None + miller_indices = reflections["miller_index"] + for idx, experiment in enumerate(experiments): + + predictor = TOFReflectionPredictor(experiment, dmin) + if predicted_reflections is None: + predicted_reflections = predictor.indices_for_ub(miller_indices) + predicted_reflections["id"] = cctbx.array_family.flex.int( + len(predicted_reflections), idx + ) + predicted_reflections["imageset_id"] = cctbx.array_family.flex.int( + len(predicted_reflections), idx + ) + else: + r = predictor.indices_for_ub(miller_indices) + r["id"] = cctbx.array_family.flex.int(len(r), idx) + r["imageset_id"] = cctbx.array_family.flex.int(len(r), idx) + predicted_reflections.extend(r) + predicted_reflections["s0"] = predicted_reflections["s0_cal"] + predicted_reflections.calculate_entering_flags(experiments) + + for i in range(len(experiments)): + predicted_reflections.experiment_identifiers()[i] = experiments[i].identifier + + # Updates flags to set which reflections to use in generating reference profiles + matched, reflections, unmatched = predicted_reflections.match_with_reference( + reflections + ) + # sel = predicted_reflections.get_flags(predicted_reflections.flags.reference_spot) + predicted_reflections = predicted_reflections.select(matched) + if "idx" in reflections: + predicted_reflections["idx"] = reflections["idx"] + + predicted_reflections["xyzobs.px.value"] = reflections["xyzobs.px.value"] + + tof_padding = params.bbox_tof_padding + xy_padding = params.bbox_xy_padding + image_size = experiments[0].detector[0].get_image_size() + tof_size = len(experiments[0].scan.get_property("time_of_flight")) + bboxes = flex.int6(len(predicted_reflections)) + partiality = flex.double(len(predicted_reflections)) + for i in range(len(predicted_reflections)): + + bboxes[i], partiality[i] = update_bounding_box( + reflections["bbox"][i], + reflections["xyzobs.px.value"][i], + predicted_reflections["xyzcal.px"][i], + (int(xy_padding), int(xy_padding), int(tof_padding)), + (0, image_size[0], 0, image_size[1], 0, tof_size), + ) + predicted_reflections["bbox"] = bboxes + predicted_reflections["partiality"] = partiality + + return predicted_reflections + + +def get_predicted_calculated_reflections(params, experiments, reflections): + + dmin = params.calculated.dmin + + predicted_reflections = None + for idx, experiment in enumerate(experiments): + + predictor = TOFReflectionPredictor(experiment, dmin) + if predicted_reflections is None: + predicted_reflections = predictor.for_ub(experiment.crystal.get_A()) + predicted_reflections["id"] = cctbx.array_family.flex.int( + len(predicted_reflections), idx + ) + predicted_reflections["imageset_id"] = cctbx.array_family.flex.int( + len(predicted_reflections), idx + ) + else: + r = predictor.for_ub(experiment.crystal.get_A()) + r["id"] = cctbx.array_family.flex.int(len(r), idx) + r["imageset_id"] = cctbx.array_family.flex.int(len(r), idx) + predicted_reflections.extend(r) + predicted_reflections["s0"] = predicted_reflections["s0_cal"] + predicted_reflections.calculate_entering_flags(experiments) + + for i in range(len(experiments)): + predicted_reflections.experiment_identifiers()[i] = experiments[i].identifier + + logger.info(f"Predicted {len(predicted_reflections)} using dmin {dmin}") + + used_in_ref = reflections.get_flags(reflections.flags.used_in_refinement) + model_reflections = reflections.select(used_in_ref) + logger.info( + f"Using {len(model_reflections)} observed reflections to calculate bounding boxes" + ) + + # Get the closest observed reflection for each calculated reflection + # Use the observed bounding box as a basis for the calculated reflection + + tof_padding = params.bbox_tof_padding + xy_padding = params.bbox_xy_padding + image_size = experiments[0].detector[0].get_image_size() + tof_size = len(experiments[0].scan.get_property("time_of_flight")) + predicted_reflections["bbox"] = flex.int6(len(predicted_reflections)) + predicted_reflections["partiality"] = flex.double(len(predicted_reflections)) + + # If no observed reflections on a given panel, use the average bbox size + x0, x1, y0, y1, z0, z1 = reflections["bbox"].parts() + num_reflections = len(reflections) + avg_bbox_size = ( + (sum(x1 - x0) / num_reflections) * 0.5, + (sum(y1 - y0) / num_reflections) * 0.5, + (sum(z1 - z0) / num_reflections) * 0.5, + ) + + for panel_idx in range(len(experiments[0].detector)): + for expt_idx in range(len(experiments)): + + p_sel = (predicted_reflections["id"] == expt_idx) & ( + predicted_reflections["panel"] == panel_idx + ) + expt_p_reflections = predicted_reflections.select(p_sel) + bboxes = flex.int6(len(expt_p_reflections)) + partiality = flex.double(len(expt_p_reflections)) + + m_sel = (model_reflections["id"] == expt_idx) & ( + model_reflections["panel"] == panel_idx + ) + expt_m_reflections = model_reflections.select(m_sel) + + for i in range(len(expt_p_reflections)): + if len(expt_m_reflections) == 0: + c = expt_p_reflections["xyzcal.px"][i] + bbox = ( + c[0] - avg_bbox_size[0], + c[0] + avg_bbox_size[0], + c[1] - avg_bbox_size[1], + c[1] + avg_bbox_size[1], + c[2] - avg_bbox_size[2], + c[2] + avg_bbox_size[2], + ) + + bboxes[i], partiality[i] = update_bounding_box( + bbox, + expt_p_reflections["xyzcal.px"][i], + expt_p_reflections["xyzcal.px"][i], + (int(xy_padding), int(xy_padding), int(tof_padding)), + (0, image_size[0], 0, image_size[1], 0, tof_size), + ) + + else: + distances = ( + expt_m_reflections["xyzobs.px.value"] + - expt_p_reflections["xyzcal.px"][i] + ).norms() + min_distance_idx = np.argmin(flumpy.to_numpy(distances)) + + bboxes[i], partiality[i] = update_bounding_box( + expt_m_reflections["bbox"][min_distance_idx], + expt_m_reflections["xyzobs.px.value"][min_distance_idx], + expt_p_reflections["xyzcal.px"][i], + (int(xy_padding), int(xy_padding), int(tof_padding)), + (0, image_size[0], 0, image_size[1], 0, tof_size), + ) + predicted_reflections["bbox"].set_selected(p_sel, bboxes) + predicted_reflections["partiality"].set_selected(p_sel, partiality) + + return predicted_reflections + + def run(): """ Input setup @@ -372,68 +560,14 @@ def run_integrate(params, experiments, reflections): Predict reflections using experiment crystal """ - min_s0_idx = min( - range(len(reflections["wavelength"])), key=reflections["wavelength"].__getitem__ - ) - min_s0 = reflections["s0"][min_s0_idx] - dmin = None - for experiment in experiments: - expt_dmin = experiment.detector.get_max_resolution(min_s0) - if dmin is None or expt_dmin < dmin: - dmin = expt_dmin - - predicted_reflections = None - miller_indices = reflections["miller_index"] - for idx, experiment in enumerate(experiments): - - predictor = TOFReflectionPredictor(experiment, dmin) - if predicted_reflections is None: - predicted_reflections = predictor.indices_for_ub(miller_indices) - predicted_reflections["id"] = cctbx.array_family.flex.int( - len(predicted_reflections), idx - ) - predicted_reflections["imageset_id"] = cctbx.array_family.flex.int( - len(predicted_reflections), idx - ) - else: - r = predictor.indices_for_ub(miller_indices) - r["id"] = cctbx.array_family.flex.int(len(r), idx) - r["imageset_id"] = cctbx.array_family.flex.int(len(r), idx) - predicted_reflections.extend(r) - predicted_reflections["s0"] = predicted_reflections["s0_cal"] - predicted_reflections.calculate_entering_flags(experiments) - - for i in range(len(experiments)): - predicted_reflections.experiment_identifiers()[i] = experiments[i].identifier - - # Updates flags to set which reflections to use in generating reference profiles - matched, reflections, unmatched = predicted_reflections.match_with_reference( - reflections - ) - # sel = predicted_reflections.get_flags(predicted_reflections.flags.reference_spot) - predicted_reflections = predicted_reflections.select(matched) - if "idx" in reflections: - predicted_reflections["idx"] = reflections["idx"] - - predicted_reflections["xyzobs.px.value"] = reflections["xyzobs.px.value"] - - tof_padding = params.bbox_tof_padding - xy_padding = params.bbox_xy_padding - image_size = experiments[0].detector[0].get_image_size() - tof_size = len(experiments[0].scan.get_property("time_of_flight")) - bboxes = flex.int6(len(predicted_reflections)) - partiality = flex.double(len(predicted_reflections)) - for i in range(len(predicted_reflections)): - - bboxes[i], partiality[i] = update_bounding_box( - reflections["bbox"][i], - reflections["xyzobs.px.value"][i], - predicted_reflections["xyzcal.px"][i], - (int(xy_padding), int(xy_padding), int(tof_padding)), - (0, image_size[0], 0, image_size[1], 0, tof_size), + if params.integration_type == "observed": + predicted_reflections = get_predicted_observed_reflections( + params=params, experiments=experiments, reflections=reflections + ) + elif params.integration_type == "calculated": + predicted_reflections = get_predicted_calculated_reflections( + params=params, experiments=experiments, reflections=reflections ) - predicted_reflections["bbox"] = bboxes - predicted_reflections["partiality"] = partiality predicted_reflections.compute_d(experiments) # predicted_reflections.compute_partiality(experiments) @@ -543,7 +677,6 @@ def run_integrate(params, experiments, reflections): tof_calculate_shoebox_mask(expt_reflections, expt) expt_reflections.is_overloaded(experiments) expt_reflections.contains_invalid_pixels() - expt_reflections["partiality"] = flex.double(len(expt_reflections), 1.0) # Background calculated explicitly to expose underlying algorithm background_algorithm = SimpleBackgroundExt( @@ -561,10 +694,8 @@ def run_integrate(params, experiments, reflections): expt_reflections = compute_line_profile_intensity(expt_reflections) predicted_reflections.set_selected(sel, expt_reflections) else: - if params.corrections.lorentz: - logger.info("Applying Lorentz correction to target run") - print("Calculating summed intensities") for idx, expt in enumerate(experiments): + print(f"Computing for experiment {idx}") sel = predicted_reflections["id"] == idx expt_reflections = predicted_reflections.select(sel) expt_data = expt.imageset @@ -578,7 +709,6 @@ def run_integrate(params, experiments, reflections): tof_calculate_shoebox_mask(expt_reflections, expt) expt_reflections.is_overloaded(experiments) expt_reflections.contains_invalid_pixels() - expt_reflections["partiality"] = flex.double(len(expt_reflections), 1.0) # Background calculated explicitly to expose underlying algorithm background_algorithm = SimpleBackgroundExt( @@ -588,11 +718,14 @@ def run_integrate(params, experiments, reflections): expt_reflections.set_flags( ~success, expt_reflections.flags.failed_during_background_modelling ) + if params.corrections.lorentz: + logger.info(" Applying Lorentz correction to target run") + print(" Calculating summed intensities") expt_reflections.compute_summed_intensity() if params.method.line_profile_fitting: - print(f"Calculating line profile fitted intensities for expt {idx}") + print(f" Calculating line profile fitted intensities for expt {idx}") expt_reflections = compute_line_profile_intensity(expt_reflections) predicted_reflections.set_selected(sel, expt_reflections) diff --git a/src/dials/util/batch_handling.py b/src/dials/util/batch_handling.py index a013ce0f0d..f31d8d0a0c 100644 --- a/src/dials/util/batch_handling.py +++ b/src/dials/util/batch_handling.py @@ -94,7 +94,10 @@ def batch_plot_shapes_and_annotations(self): def assign_batches_to_reflections(reflections, batch_offsets): """Assign a 'batch' column to the reflection table""" for batch_offset, refl in zip(batch_offsets, reflections): - xdet, ydet, zdet = [flex.double(x) for x in refl["xyzobs.px.value"].parts()] + if "xyzobs.px.value" in refl: + xdet, ydet, zdet = [flex.double(x) for x in refl["xyzobs.px.value"].parts()] + else: + xdet, ydet, zdet = [flex.double(x) for x in refl["xyzcal.px"].parts()] # compute BATCH values - floor() to get (fortran) image captured within # +1 because FORTRAN counting; zdet+1=image_index # +off because image_index+o=batch diff --git a/src/dials/util/export_mtz.py b/src/dials/util/export_mtz.py index 13469a98ae..0c1a200ab3 100644 --- a/src/dials/util/export_mtz.py +++ b/src/dials/util/export_mtz.py @@ -293,7 +293,7 @@ def add_batch_list( umat_array[i, j] = U_t_elements[j] # We ignore panels beyond the first one, at the moment - panel = experiment.detector[0] + panel = experiment.detector[6] panel_size = panel.get_image_size() panel_distance = panel.get_directed_distance() @@ -549,16 +549,29 @@ def write_columns_tof(mtz, reflection_table): nref = len(reflection_table["miller_index"]) assert nref - xdet, ydet, _ = [ - flex.double(x) for x in reflection_table["xyzobs.px.value"].parts() - ] + if "xyzcal.px" in reflection_table: + xdet, ydet, _ = [flex.double(x) for x in reflection_table["xyzcal.px"].parts()] + else: + xdet, ydet, _ = [ + flex.double(x) for x in reflection_table["xyzobs.px.value"].parts() + ] type_table = { "H": "H", "K": "H", "L": "H", + "PACK_ID": "B", + "PLATE": "I", + "XF": "R", + "YF": "R", + "LAMBDA": "R", "I": "J", "SIGI": "Q", + "MULT": "I", + "MINHARM": "I", + "MAXHARM": "I", + "NOVPIX": "R", + "FLAGS": "I", "IPR": "J", "SIGIPR": "Q", "BG": "R", @@ -575,42 +588,35 @@ def write_columns_tof(mtz, reflection_table): "FRACTIONCALC": "R", "ROT": "R", "QE": "R", - "LAMBDA": "R", } mtz_data = pd.DataFrame( flumpy.to_numpy(reflection_table["miller_index"]).astype("float32"), columns=["H", "K", "L"], ) - mtz_data.insert(3, "M/ISYM", np.zeros(nref, dtype="float32")) + mtz.add_column("PACK_ID", type_table["PACK_ID"]) + mtz_data.insert( + len(mtz_data.columns), "PACK_ID", flumpy.to_numpy(reflection_table["panel"]) + ) - # H, K, L are in the base dataset, but we have to add M/ISYM - mtz.add_column("M/ISYM", type_table["M_ISYM"]) - mtz.add_column("BATCH", type_table["BATCH"]) + mtz.add_column("PLATE", type_table["PLATE"]) + mtz_data.insert(len(mtz_data.columns), "PLATE", np.ones(nref).astype("int32")) + + mtz.add_column("XF", type_table["XF"]) mtz_data.insert( - 4, "BATCH", flumpy.to_numpy(reflection_table["batch"]).astype("float32") + len(mtz_data.columns), "XF", flumpy.to_numpy(xdet).astype("float32") + ) + mtz.add_column("YF", type_table["YF"]) + mtz_data.insert( + len(mtz_data.columns), "YF", flumpy.to_numpy(ydet).astype("float32") ) - if "intensity.prf.value" in reflection_table: - if "intensity.sum.value" in reflection_table: - col_names = ("IPR", "SIGIPR") - else: - col_names = ("I", "SIGI") - I_profile = reflection_table["intensity.prf.value"] - V_profile = reflection_table["intensity.prf.variance"] - assert V_profile.all_gt(0) # Trap negative variances - mtz.add_column(col_names[0], type_table["I"]) - mtz_data.insert( - len(mtz_data.columns), - col_names[0], - flumpy.to_numpy(I_profile.as_float()).astype("float32"), - ) - mtz.add_column(col_names[1], type_table["SIGI"]) - mtz_data.insert( - len(mtz_data.columns), - col_names[1], - flumpy.to_numpy(flex.sqrt(V_profile)).astype("float32"), - ) + mtz.add_column("LAMBDA", type_table["LAMBDA"]) + mtz_data.insert( + len(mtz_data.columns), + "LAMBDA", + flumpy.to_numpy(reflection_table["wavelength_cal"]).astype("float32"), + ) if "intensity.sum.value" in reflection_table: I_sum = reflection_table["intensity.sum.value"] @@ -627,6 +633,48 @@ def write_columns_tof(mtz, reflection_table): flumpy.to_numpy(flex.sqrt(V_sum)).astype("float32"), ) + mtz.add_column("MULT", type_table["MULT"]) + mult_dict = Counter(reflection_table["miller_index"]) + mult = [mult_dict[reflection_table["miller_index"][i]] for i in range(nref)] + mtz_data.insert(len(mtz_data.columns), "MULT", np.array(mult).astype("int32")) + + mtz.add_column("MINHARM", type_table["MINHARM"]) + mtz_data.insert(len(mtz_data.columns), "MINHARM", np.ones(nref).astype("int32")) + mtz.add_column("MAXHARM", type_table["MAXHARM"]) + mtz_data.insert(len(mtz_data.columns), "MAXHARM", np.ones(nref).astype("int32")) + + mtz.add_column("NOVPIX", type_table["NOVPIX"]) + x0, x1, y0, y1, z0, z1 = reflection_table["bbox"].parts() + novpix = (x1 - x0) * (y1 - y0) * (z1 - z0) + mtz_data.insert( + len(mtz_data.columns), "NOVPIX", flumpy.to_numpy(novpix).astype("int32") + ) + + mtz_data.insert(len(mtz_data.columns), "M/ISYM", np.zeros(nref, dtype="float32")) + + # H, K, L are in the base dataset, but we have to add M/ISYM + mtz.add_column("M/ISYM", type_table["M_ISYM"]) + mtz.add_column("BATCH", type_table["BATCH"]) + mtz_data.insert( + len(mtz_data.columns), + "BATCH", + flumpy.to_numpy(reflection_table["batch"]).astype("float32"), + ) + + if "intensity.prf.value" in reflection_table: + I_prf = reflection_table["intensity.prf.value"] + assert V_sum.all_gt(0) # Trap negative variances + mtz.add_column("IPR", type_table["IPR"]) + mtz_data.insert( + len(mtz_data.columns), "IPR", flumpy.to_numpy(I_prf).astype("float32") + ) + mtz.add_column("SIGIPR", type_table["SIGIPR"]) + mtz_data.insert( + len(mtz_data.columns), + "SIGIPR", + flumpy.to_numpy(flex.sqrt(V_sum)).astype("float32"), + ) + if ( "background.sum.value" in reflection_table and "background.sum.variance" in reflection_table @@ -652,13 +700,6 @@ def write_columns_tof(mtz, reflection_table): flumpy.to_numpy(reflection_table["fractioncalc"]).astype("float32"), ) - mtz.add_column("LAMBDA", type_table["LAMBDA"]) - mtz_data.insert( - len(mtz_data.columns), - "LAMBDA", - flumpy.to_numpy(reflection_table["wavelength_cal"]).astype("float32"), - ) - mtz.add_column("XDET", type_table["XDET"]) mtz_data.insert( len(mtz_data.columns), "XDET", flumpy.to_numpy(xdet).astype("float32") @@ -668,39 +709,6 @@ def write_columns_tof(mtz, reflection_table): len(mtz_data.columns), "YDET", flumpy.to_numpy(ydet).astype("float32") ) - if "ROT" in reflection_table: - mtz.add_column("ROT", type_table["ROT"]) - mtz_data.insert( - len(mtz_data.columns), - "ROT", - flumpy.to_numpy(reflection_table["ROT"]).astype("float32"), - ) - - if "lp" in reflection_table: - mtz.add_column("LP", type_table["LP"]) - mtz_data.insert( - len(mtz_data.columns), - "LP", - flumpy.to_numpy(reflection_table["lp"]).astype("float32"), - ) - if "qe" in reflection_table: - mtz.add_column("QE", type_table["QE"]) - mtz_data.insert( - len(mtz_data.columns), - "QE", - flumpy.to_numpy(reflection_table["qe"]).astype("float32"), - ) - elif "dqe" in reflection_table: - mtz.add_column("QE", type_table["QE"]) - mtz_data.insert( - len(mtz_data.columns), - "QE", - flumpy.to_numpy(reflection_table["dqe"]).astype("float32"), - ) - else: - mtz.add_column("QE", type_table["QE"]) - mtz_data.insert(len(mtz_data.columns), "QE", np.ones(nref).astype("float32")) - mtz.switch_to_original_hkl() mtz.set_data(mtz_data)