diff --git a/xfel/small_cell/command_line/powder_from_spots.py b/xfel/small_cell/command_line/powder_from_spots.py index 933fbe1c51a..0f1940366ba 100644 --- a/xfel/small_cell/command_line/powder_from_spots.py +++ b/xfel/small_cell/command_line/powder_from_spots.py @@ -95,6 +95,9 @@ .type = space_group .help = Show positions of miller indices from this unit_cell and space \ group. Not implemented. + n_max = None + .type = int + .help = Stop plotting after this many experiments. filter { enable = False .type = bool @@ -121,6 +124,10 @@ .type = str xy_file = None .type = str + xy_file_units = *d q + .type = choice + .help = X-axis units for xy_file output: d-spacing in Angstroms (default) \ + or inverse-d-spacing (q) in inverse Angstroms. peak_file = None .type = str .help = Optionally, specify an output file for interactive peak picking in \ @@ -139,9 +146,21 @@ .type = float d_max = None .type = float + q_min = None + .type = float + .help = Minimum q (inverse angstroms). Converted to d_max if provided. + q_max = None + .type = float + .help = Maximum q (inverse angstroms). Converted to d_min if provided. step_px = None .type = float .multiple = True + d_target = None + .type = float + .help = If set, enables detector distance refinement along z-axis + q_target = None + .type = float + .help = Target q value (inverse angstroms). Converted to d_target if provided. } plot { interactive = True @@ -185,11 +204,33 @@ def run(self): experiments = params.input.experiments[0].data reflections = params.input.reflections[0].data + # Convert q parameters to d parameters if needed (q = 1/d) + if params.center_scan.q_min is not None: + params.center_scan.d_max = 1.0 / params.center_scan.q_min + if params.center_scan.q_max is not None: + params.center_scan.d_min = 1.0 / params.center_scan.q_max + if params.center_scan.q_target is not None: + params.center_scan.d_target = 1.0 / params.center_scan.q_target + if params.center_scan.d_min: assert params.center_scan.d_max cscan = Center_scan(experiments, reflections, params) + + # First XY refinement sequence for step in params.center_scan.step_px: cscan.search_step(step) + + # Z-axis refinement if d_target is set + if params.center_scan.d_target is not None: + # Default Z steps in microns + z_steps_um = [4000, 2000, 1000, 500, 250, 125, 80, 40, 20, 10, 5] + for step_um in z_steps_um: + cscan.search_step_z(step_um, params.center_scan.d_target) + + # Second XY refinement sequence + for step in params.center_scan.step_px: + cscan.search_step(step) + if params.output.geom_file is not None: experiments.as_file(params.output.geom_file) diff --git a/xfel/small_cell/command_line/powder_refine_geometry.py b/xfel/small_cell/command_line/powder_refine_geometry.py new file mode 100644 index 00000000000..26338239c55 --- /dev/null +++ b/xfel/small_cell/command_line/powder_refine_geometry.py @@ -0,0 +1,158 @@ +# LIBTBX_SET_DISPATCHER_NAME cctbx.xfel.powder_refine_geometry +from __future__ import division +import logging + +from iotbx.phil import parse +from dials.util import log +from dials.util import show_mail_on_error +from dials.util.options import ArgumentParser +from xfel.small_cell.geometry_refiner import PowderGeometryRefiner + + +logger = logging.getLogger("dials.command_line.powder_refine_geometry") + +help_message = """ +Refine detector geometry using powder diffraction d-spacings. + +Examples of usage: + +# Basic usage with defaults (LaB6) +$ cctbx.xfel.powder_refine_geometry spots.expt spots.refl + +# Custom d-spacings (e.g., for Si standard) +$ cctbx.xfel.powder_refine_geometry spots.expt spots.refl \\ + reference_d_spacings=3.135,1.920,1.637,1.357 + +# Or use unit cell and space group instead +$ cctbx.xfel.powder_refine_geometry spots.expt spots.refl \\ + unit_cell=4.156,4.156,4.156,90,90,90 space_group=Pm-3m + +# Refine only XY shift (fix distance and tilt) +$ cctbx.xfel.powder_refine_geometry spots.expt spots.refl \\ + refine.distance=False refine.tilt=False + +# Specify output file +$ cctbx.xfel.powder_refine_geometry spots.expt spots.refl \\ + output.experiments=calibrated.expt + +This tool uses spotfinding output from a powder standard with known d-spacings +to refine detector geometry. It optimizes a 5-parameter detector model: +- distance: shift along detector normal +- shift1/shift2: XY shifts along fast/slow axes +- tau2/tau3: tilts around fast/slow axes + +The refinement minimizes the sum of squared differences between observed +d-spacings and the nearest reference d-spacing. + +Default reference d-spacings are for LaB6 (SRM 660): + 4.156 A (100), 2.939 A (110), 2.399 A (111), 2.078 A (200), 1.858 A (210) +""" + +phil_scope = parse( + """ +reference_d_spacings = 4.156 2.939 2.399 2.078 1.858 + .type = floats + .help = Reference d-spacings in Angstroms. Default: first 5 LaB6 peaks. \ + Either use this OR specify unit_cell and space_group. + +unit_cell = None + .type = unit_cell + .help = Unit cell to generate reference d-spacings (use with space_group). \ + If specified, this overrides reference_d_spacings. + +space_group = None + .type = space_group + .help = Space group to generate reference d-spacings (use with unit_cell). \ + If specified, this overrides reference_d_spacings. + +d_min = 1.5 + .type = float + .help = Minimum d-spacing to include in refinement + +d_max = 20 + .type = float + .help = Maximum d-spacing to include in refinement + +max_distance_inv_ang = 0.002 + .type = float + .help = Maximum distance from reference d-spacing in inverse Angstroms. \ + Reflections further than this from any reference will be excluded. + +refine { + distance = True + .type = bool + .help = Refine detector distance along normal + shift = True + .type = bool + .help = Refine XY shift (shift1 and shift2) + tilt = True + .type = bool + .help = Refine detector tilts (tau2 and tau3) +} + +output { + experiments = refined.expt + .type = path + .help = Output filename for refined experiments + log = powder_refine_geometry.log + .type = path + .help = Output log file +} +""" +) + + +class Script(object): + def __init__(self): + usage = "$ cctbx.xfel.powder_refine_geometry EXPERIMENTS REFLECTIONS [options]" + self.parser = ArgumentParser( + usage=usage, + phil=phil_scope, + epilog=help_message, + check_format=False, + read_reflections=True, + read_experiments=True, + ) + + def run(self): + params, options = self.parser.parse_args(show_diff_phil=True) + + # Validate input + if len(params.input.experiments) != 1 or len(params.input.reflections) != 1: + raise ValueError("Please provide exactly one experiments file and " + "one reflections file") + + experiments = params.input.experiments[0].data + reflections = params.input.reflections[0].data + + print(f"\nLoaded {len(experiments)} experiments with " + f"{len(reflections)} reflections") + + if params.unit_cell is not None and params.space_group is not None: + print(f"Using unit_cell: {params.unit_cell}") + print(f"Using space_group: {params.space_group}") + else: + print(f"Reference d-spacings: {params.reference_d_spacings}") + + # Create and run refiner + refiner = PowderGeometryRefiner(experiments, reflections, params) + result = refiner.run() + + if result is not None: + # Report final geometry + refiner.report_geometry_changes() + + # Save refined experiments + refined_experiments = refiner.get_refined_experiments() + print(f"\nSaving refined experiments to {params.output.experiments}") + refined_experiments.as_file(params.output.experiments) + + print("\nRefinement complete.") + else: + print("\nNo refinement performed.") + + +if __name__ == "__main__": + with show_mail_on_error(): + script = Script() + script.run() diff --git a/xfel/small_cell/geometry_refiner.py b/xfel/small_cell/geometry_refiner.py new file mode 100644 index 00000000000..03454dd14d7 --- /dev/null +++ b/xfel/small_cell/geometry_refiner.py @@ -0,0 +1,366 @@ +from __future__ import division +import numpy as np + +from dials.array_family import flex +from cctbx import uctbx, miller +from scitbx import matrix +from scipy.optimize import minimize + + +class PowderGeometryRefiner: + """ + Refines detector geometry using powder diffraction d-spacings. + + Uses a 5-parameter model: + - dist: distance from origin to detector plane along normal (mm) + - shift1: shift along detector fast axis (mm) + - shift2: shift along detector slow axis (mm) + - tau2: rotation around fast axis (mrad) - tilt + - tau3: rotation around slow axis (mrad) - tilt + + NOT refined: tau1 (rotation around normal) - indeterminate from powder + """ + + def __init__(self, experiments, reflections, params): + self.experiments = experiments + self.reflections = reflections + self.params = params + + # Compute reference d-spacings from unit_cell and space_group if provided + if params.unit_cell is not None and params.space_group is not None: + print(f"Computing d-spacings from unit_cell={params.unit_cell} " + f"and space_group={params.space_group}") + + # Get beam and detector for resolution calculation + beam = experiments[0].beam + detector = experiments[0].detector + d_min = params.d_min + d_max = params.d_max + + # Average unit cell over space group + unit_cell = params.unit_cell + + # Generate Miller indices within resolution range + generator = miller.index_generator(unit_cell, params.space_group.type(), False, d_min) + indices = generator.to_array() + + # Compute d-spacings and filter by resolution range + all_spacings = unit_cell.d(indices) + spacings_in_range = flex.sorted(flex.double([d for d in all_spacings if d_min <= d <= d_max])) + + self.reference_d = np.array(spacings_in_range) + print(f"Generated {len(self.reference_d)} reference d-spacings in range " + f"[{d_min:.3f}, {d_max:.3f}] A") + else: + self.reference_d = np.array(params.reference_d_spacings) + + # Store initial detector state + self.detector = experiments[0].detector + self.initial_state = self._get_detector_state() + + # Track which parameters to refine + self.refine_distance = params.refine.distance + self.refine_shift = params.refine.shift + self.refine_tilt = params.refine.tilt + + # Build parameter vector and bounds + self._setup_parameters() + + # Prepare reflections + self._prepare_reflections() + + def _get_detector_state(self): + """Extract fast, slow, origin from detector hierarchy.""" + hierarchy = self.detector.hierarchy() + fast = matrix.col(hierarchy.get_local_fast_axis()) + slow = matrix.col(hierarchy.get_local_slow_axis()) + origin = matrix.col(hierarchy.get_local_origin()) + normal = fast.cross(slow) + + # For multipanel detectors, compute center using first panel dimensions + # (this is just for rotation center - all panels move together) + panel = self.detector[0] + size = panel.get_image_size() + pixel_size = panel.get_pixel_size() + center = origin + (size[0]/2 * pixel_size[0]) * fast + (size[1]/2 * pixel_size[1]) * slow + + return { + 'fast': fast, + 'slow': slow, + 'origin': origin, + 'normal': normal, + 'center': center, + } + + def _setup_parameters(self): + """Set up parameter vector based on what is being refined.""" + # Initial values: [dist, shift1, shift2, tau2, tau3] + # All start at 0 (representing no change from initial geometry) + self.param_names = [] + self.initial_params = [] + self.bounds = [] + + if self.refine_distance: + self.param_names.append('dist') + self.initial_params.append(0.0) + self.bounds.append((-50.0, 50.0)) # +/- 50 mm + + if self.refine_shift: + self.param_names.append('shift1') + self.initial_params.append(0.0) + self.bounds.append((-10.0, 10.0)) # +/- 10 mm + + self.param_names.append('shift2') + self.initial_params.append(0.0) + self.bounds.append((-10.0, 10.0)) # +/- 10 mm + + if self.refine_tilt: + self.param_names.append('tau2') + self.initial_params.append(0.0) + self.bounds.append((-50.0, 50.0)) # +/- 50 mrad + + self.param_names.append('tau3') + self.initial_params.append(0.0) + self.bounds.append((-50.0, 50.0)) # +/- 50 mrad + + self.initial_params = np.array(self.initial_params) + + def _prepare_reflections(self): + """Filter reflections by d-range and proximity to reference d-spacings.""" + refls = self.reflections + + # Compute initial d-spacings to filter + refls.centroid_px_to_mm(self.experiments) + refls.map_centroids_to_reciprocal_space(self.experiments) + d_star_sq = flex.pow2(refls['rlp'].norms()) + refls['d'] = uctbx.d_star_sq_as_d(d_star_sq) + + # Filter by d-range + d_min = self.params.d_min + d_max = self.params.d_max + sel = (refls['d'] >= d_min) & (refls['d'] <= d_max) + refls_in_range = refls.select(sel) + + print(f"Found {len(refls_in_range)} reflections in d-range " + f"[{d_min:.3f}, {d_max:.3f}] A") + + # Filter by proximity to reference d-spacings (in inverse angstroms) + max_dist = self.params.max_distance_inv_ang + dvals = refls_in_range['d'].as_numpy_array() + dvals_inv = 1.0 / dvals + ref_d_inv = 1.0 / self.reference_d + + # For each reflection, find minimum distance to any reference + min_distances = np.min( + np.abs(dvals_inv[:, np.newaxis] - ref_d_inv[np.newaxis, :]), + axis=1 + ) + close_to_ref = min_distances <= max_dist + sel_close = flex.bool(close_to_ref.tolist()) + self.refls_filtered = refls_in_range.select(sel_close) + + print(f"Using {len(self.refls_filtered)} reflections within " + f"{max_dist:.4f} inv. A of reference d-spacings") + + # Store pixel positions for later use + self.xyzobs_px = self.refls_filtered['xyzobs.px.value'] + self.panels = self.refls_filtered['panel'] + self.ids = self.refls_filtered['id'] + + # Pre-compute unique experiment IDs and create subset for optimization + unique_exp_ids_list = sorted(set(self.refls_filtered['id'])) + print(f"Filtered reflections span {len(unique_exp_ids_list)} unique experiments " + f"out of {len(self.experiments)}") + + # Create a subset of experiments and remap reflection IDs + # This dramatically speeds up coordinate transformations + if len(unique_exp_ids_list) < len(self.experiments): + from dxtbx.model.experiment_list import ExperimentList + self.experiments_subset = ExperimentList() + old_to_new_id = {} + for new_id, old_id in enumerate(unique_exp_ids_list): + self.experiments_subset.append(self.experiments[old_id]) + old_to_new_id[old_id] = new_id + + # Remap reflection IDs to the subset + old_ids = self.refls_filtered['id'] + new_ids = flex.size_t([old_to_new_id[old_id] for old_id in old_ids]) + self.refls_filtered['id'] = new_ids + + # Use the subset for all operations + self.experiments = self.experiments_subset + print(f"Created experiment subset with {len(self.experiments)} experiments") + else: + print(f"Using all experiments (no subset needed)") + + def apply_params(self, x): + """Apply parameter vector to detector geometry.""" + # Parse parameter vector + params_dict = {} + for i, name in enumerate(self.param_names): + params_dict[name] = x[i] + + # Convert to Python floats for scitbx matrix compatibility + dist = float(params_dict.get('dist', 0.0)) + shift1 = float(params_dict.get('shift1', 0.0)) + shift2 = float(params_dict.get('shift2', 0.0)) + tau2 = float(params_dict.get('tau2', 0.0)) / 1000.0 # mrad to rad + tau3 = float(params_dict.get('tau3', 0.0)) / 1000.0 # mrad to rad + + # Get initial state + d1 = self.initial_state['fast'] + d2 = self.initial_state['slow'] + dn = self.initial_state['normal'] + origin = self.initial_state['origin'] + center = self.initial_state['center'] + + # Apply rotations around the panel CENTER (not origin) + # tau2: rotation around fast axis (d1) + # tau3: rotation around slow axis (d2) + if abs(tau2) > 1e-10 or abs(tau3) > 1e-10: + # Use axis-angle rotation around the actual detector axes + # This correctly handles the left-handed detector frame + R2 = d1.axis_and_angle_as_r3_rotation_matrix(tau2, deg=False) + R3 = d2.axis_and_angle_as_r3_rotation_matrix(tau3, deg=False) + + # Combined rotation: first tau2 around fast, then tau3 around slow + R = R3 * R2 + + # Apply rotation to axes + d1_new = R * d1 + d2_new = R * d2 + + # Rotate origin around center to preserve panel center position + origin_rotated = center + R * (origin - center) + else: + d1_new = d1 + d2_new = d2 + origin_rotated = origin + + # Apply shifts and distance change (relative to initial axes) + new_origin = (origin_rotated + + shift1 * d1 + + shift2 * d2 + + dist * dn) + + # Update detector hierarchy (moves all panels together as rigid body) + hierarchy = self.detector.hierarchy() + hierarchy.set_local_frame( + d1_new.elems, + d2_new.elems, + new_origin.elems + ) + + def compute_dvals(self): + """Compute d-spacings for all reflections with current geometry.""" + # Reset mm coordinates to force recalculation + if 'xyzobs.mm.value' in self.refls_filtered: + del self.refls_filtered['xyzobs.mm.value'] + if 'rlp' in self.refls_filtered: + del self.refls_filtered['rlp'] + + self.refls_filtered.centroid_px_to_mm(self.experiments) + self.refls_filtered.map_centroids_to_reciprocal_space(self.experiments) + d_star_sq = flex.pow2(self.refls_filtered['rlp'].norms()) + self.dvals = uctbx.d_star_sq_as_d(d_star_sq) + return self.dvals + + def find_nearest_reference(self, d_obs): + """Find closest reference d-spacing.""" + distances = np.abs(self.reference_d - d_obs) + return self.reference_d[np.argmin(distances)] + + def find_nearest_references_vectorized(self, d_obs_array): + """Find closest reference d-spacing for array of observations (vectorized).""" + # d_obs_array shape: (n_obs,) + # reference_d shape: (n_ref,) + # Compute distances matrix: (n_obs, n_ref) + distances = np.abs(d_obs_array[:, np.newaxis] - self.reference_d[np.newaxis, :]) + # Find index of minimum distance for each observation + nearest_indices = np.argmin(distances, axis=1) + return self.reference_d[nearest_indices] + + def objective(self, x): + """Compute sum of squared residuals.""" + self.apply_params(x) + dvals = self.compute_dvals() + + # Vectorized computation of residuals + dvals_np = dvals.as_numpy_array() if hasattr(dvals, 'as_numpy_array') else np.array(dvals) + d_ref = self.find_nearest_references_vectorized(dvals_np) + residuals = dvals_np - d_ref + + # Progress counter + if not hasattr(self, '_obj_eval_count'): + self._obj_eval_count = 0 + self._obj_eval_count += 1 + if self._obj_eval_count % 10 == 0: + print(f" Objective evaluation {self._obj_eval_count}, value: {np.sum(residuals ** 2):.6f}") + + return np.sum(residuals ** 2) + + def run(self): + """Run refinement and return results.""" + if len(self.param_names) == 0: + print("No parameters selected for refinement!") + return None + + # Compute initial objective + initial_obj = self.objective(self.initial_params) + print(f"\nInitial objective: {initial_obj:.6f}") + print(f"Initial RMS residual: {np.sqrt(initial_obj / len(self.refls_filtered)):.6f} A") + + # Run minimization using Powell method (derivative-free, more robust) + print(f"\nRefining {len(self.param_names)} parameters: {self.param_names}") + result = minimize( + self.objective, + x0=self.initial_params, + method='Powell', + options={'maxiter': 50, 'disp': True, 'xtol': 0.01, 'ftol': 0.001} + ) + + # Apply final parameters + self.apply_params(result.x) + final_dvals = self.compute_dvals() + + # Compute final statistics + final_obj = result.fun + print(f"\nFinal objective: {final_obj:.6f}") + print(f"Final RMS residual: {np.sqrt(final_obj / len(self.refls_filtered)):.6f} A") + print(f"Total objective evaluations: {self._obj_eval_count}") + + # Report parameter changes + print("\nRefined parameters:") + for i, name in enumerate(self.param_names): + value = result.x[i] + unit = "mm" if name in ['dist', 'shift1', 'shift2'] else "mrad" + print(f" {name}: {value:.4f} {unit}") + + return result + + def get_refined_experiments(self): + """Return experiments with refined detector.""" + return self.experiments + + def report_geometry_changes(self): + """Print summary of geometry changes.""" + hierarchy = self.detector.hierarchy() + new_origin = matrix.col(hierarchy.get_local_origin()) + new_fast = matrix.col(hierarchy.get_local_fast_axis()) + new_slow = matrix.col(hierarchy.get_local_slow_axis()) + + old_origin = self.initial_state['origin'] + old_fast = self.initial_state['fast'] + old_slow = self.initial_state['slow'] + + origin_delta = new_origin - old_origin + + # Compute tilt changes from axis differences + # tau2 (tilt around fast) shows in slow axis z-component + # tau3 (tilt around slow) shows in fast axis z-component + tau2_change = np.arcsin(new_slow[2]) - np.arcsin(old_slow[2]) + tau3_change = np.arcsin(-new_fast[2]) - np.arcsin(-old_fast[2]) + + print("\nGeometry changes:") + print(f" Origin shift: ({origin_delta[0]:.4f}, {origin_delta[1]:.4f}, {origin_delta[2]:.4f}) mm") + print(f" Tilt changes: tau2={tau2_change*1000:.3f} mrad, tau3={tau3_change*1000:.3f} mrad") diff --git a/xfel/small_cell/powder_util.py b/xfel/small_cell/powder_util.py index fcdc3d85326..57530b2f604 100644 --- a/xfel/small_cell/powder_util.py +++ b/xfel/small_cell/powder_util.py @@ -3,7 +3,6 @@ import numpy as np import matplotlib.pyplot as plt import matplotlib.ticker as tick -import numpy as np import copy from dials.array_family import flex from cctbx import uctbx @@ -44,8 +43,9 @@ def _process_pixel(self, i_panel, s0, panel, xy, value): i_bin = int( n_bins * (res_inv - d_max_inv ) / (d_min_inv - d_max_inv) ) - if i_bin < 0 or i_bin >= n_bins: return - self.current_panelsums[i_panel][i_bin] += value + if 0 <= i_bin < n_bins: + self.current_panelsums[i_panel][i_bin] += value + return res def _nearest_peak(self, x, xvalues, yvalues): i = np.searchsorted(xvalues, x, side="left") @@ -115,32 +115,28 @@ def calculate(self): else: self.use_current_expt = True if i % 1000 == 0: print("experiment ", i) + if self.params.n_max is not None and i >= self.params.n_max: + break s0 = expt.beam.get_s0() sel = refls['id'] == i refls_sel = refls.select(sel) xyzobses = refls_sel['xyzobs.px.value'] intensities = refls_sel['intensity.sum.value'] panels = refls_sel['panel'] - shoeboxes = refls_sel['shoebox'] for i_refl in range(len(refls_sel)): self.expt_count += 1 i_panel = panels[i_refl] panel = expt.detector[i_panel] - peak_height = intensities[i_refl] if params.peak_position=="xyzobs": xy = xyzobses[i_refl][0:2] if params.peak_weighting == "intensity": value = intensities[i_refl] else: value = 1 - self._process_pixel(i_panel, s0, panel, xy, value) - if params.peak_position=="shoebox": - sb = shoeboxes[i_refl] - sbpixels = zip(sb.coords(), sb.values()) - for (x,y,_), value in sbpixels: - self._process_pixel(i_panel, s0, panel, (x,y), value) + res = self._process_pixel(i_panel, s0, panel, xy, value) + for i in range(len(self.panelsums)): self.panelsums[i] = self.panelsums[i] + self.current_panelsums[i] if self.params.filter.enable and self.params.filter.select_mode=='any': @@ -167,7 +163,7 @@ def plot(self): d_max_inv = 1/params.d_max d_min_inv = 1/params.d_min xvalues = np.linspace(d_max_inv, d_min_inv, params.n_bins) - fig, ax = plt.subplots() + fig, ax = plt.subplots(figsize=(10, 3)) ps_maxes = [max(ps) for ps in self.panelsums] ps_max = max(ps_maxes) @@ -178,6 +174,8 @@ def plot(self): yvalues = np.array(sums) plt.plot(xvalues, yvalues+0.5*i_sums*offset) elif params.filter.enable and params.filter.plot_mode=="ratio": + print('filtered: ', self.filtered_expt_count) + print('antifiltered: ', self.antifiltered_expt_count) for x in self.filtered_panelsums: x /= self.filtered_expt_count for x in self.antifiltered_panelsums: @@ -200,8 +198,10 @@ def plot(self): if params.output.xy_file: with open(params.output.xy_file, 'w') as f: + use_q = getattr(params.output, 'xy_file_units', 'd') == 'q' for x,y in zip(xvalues, yvalues): - f.write("{:.6f}\t{}\n".format(1/x, y)) + xout = x if use_q else 1/x + f.write("{:.6f}\t{}\n".format(xout, y)) # Now plot the predicted peak positions if requested if params.unit_cell or params.space_group: @@ -211,18 +211,21 @@ def plot(self): ) hkl_list = cctbx.miller.build_set(sym, False, d_min=params.d_min) dspacings = params.unit_cell.d(hkl_list.indices()) +# for hkl, d in sorted(zip(hkl_list.indices(), dspacings), key=lambda x:x[1]): +# print('{:.3f}: {}'.format(d, hkl)) dspacings_inv = 1/dspacings - pplot_min = -.05*ps_max + pplot_min = -.05*max(yvalues) for d in dspacings_inv: plt.plot((d,d),(pplot_min,0), 'r-', linewidth=1) if params.output.plot_file: + fig.tight_layout() plt.savefig(params.output.plot_file) if params.plot.interactive and params.output.peak_file: backend_list = ["TkAgg","QtAgg"] assert (plt.get_backend() in backend_list), """Matplotlib backend not compatible with interactive peak picking. -You can set the MPLBACKEND environment varibale to change this. +You can set the MPLBACKEND environment variable to change this. Currently supported options: %s""" %backend_list #If a peak list output file is specified, do interactive peak picking: with open(params.output.peak_file, 'w') as f: @@ -277,6 +280,7 @@ def __init__(self, experiments, reflections, params): assert len(self.experiments.detectors())==1 self.net_origin_shift = np.array([0.,0.,0.]) + self.net_distance_shift = 0.0 self.centroid_px_mm_done = False self.px_size = self.experiments.detectors()[0][0].get_pixel_size() self.target_refl_count = 0 @@ -325,6 +329,16 @@ def width(self): print(f'width {result:.5f} from {sel.count(True)} dvals') return result + def mean_squared_error_from_target(self, d_target): + """Calculate mean squared difference between d-spacings and target d-spacing""" + if len(self.dvals) < 3: return 999 + dvals_array = flex.double(self.dvals) + differences = dvals_array - d_target + mse = flex.mean(differences * differences) + rmse = np.sqrt(mse) + print(f'RMSE from target {d_target:.5f}: {rmse:.5f} from {len(self.dvals)} dvals') + return mse + def search_step(self, step_px, nsteps=3, update=True): step_size_mm = np.array(self.px_size + (0.,)) * step_px assert nsteps%2 == 1, "nsteps should be odd" @@ -366,6 +380,74 @@ def search_step(self, step_px, nsteps=3, update=True): print(f'net shift: {self.net_origin_shift}') return width_start, width_end, self.net_origin_shift + def search_step_z(self, step_um, d_target, nsteps=3, update=True, min_refl_fraction=0.8): + """Search along detector z-axis (distance) to minimize MSE from d_target + + Args: + step_um: step size in microns + d_target: target d-spacing value + nsteps: number of steps to search (default 3: -1, 0, +1) + update: whether to update the geometry with the best result + min_refl_fraction: reject steps if reflection count drops below this fraction (default 0.8) + """ + step_size_mm = step_um / 1000.0 # convert microns to mm + assert nsteps % 2 == 1, "nsteps should be odd" + step_min = -1 * (nsteps // 2) + step_max = nsteps // 2 + 1e-6 # make the range inclusive + step_arange = np.arange(step_min, step_max) + + detector = self.experiments.detectors()[0] + hierarchy = detector.hierarchy() + fast = hierarchy.get_local_fast_axis() + slow = hierarchy.get_local_slow_axis() + origin = hierarchy.get_local_origin() + + # Compute detector normal (z-axis direction) + # Normal points from sample toward detector + normal = np.cross(fast, slow) + normal = normal / np.linalg.norm(normal) + + results = [] + self.update_dvals() + initial_count = len(self.dvals) + mse_start = self.mean_squared_error_from_target(d_target) + print(f'Z-axis search start (MSE): {mse_start:.5f}, n_refl: {initial_count}') + + for step_idx in step_arange: + step_mm = step_idx * step_size_mm + # Shift along normal direction + origin_shift = normal * step_mm + new_origin = origin + origin_shift + hierarchy.set_local_frame(fast, slow, new_origin) + self.update_dvals() + current_count = len(self.dvals) + + # Reject if reflection count dropped by more than (1 - min_refl_fraction) + if current_count < min_refl_fraction * initial_count: + result = 999 # penalty value + print(f' Step {step_mm:.6f} mm: rejected (n_refl={current_count}, {100*current_count/initial_count:.1f}%)') + else: + result = self.mean_squared_error_from_target(d_target) + results.append(result) + + mse_end = min(results) + i_best = results.index(mse_end) + distance_shift = step_arange[i_best] * step_size_mm + if update: + origin_shift = normal * distance_shift + new_origin = origin + origin_shift + self.net_origin_shift += origin_shift + self.net_distance_shift += distance_shift + else: + new_origin = origin + hierarchy.set_local_frame(fast, slow, new_origin) + self.update_dvals() + final_count = len(self.dvals) + print(f'Z step: {distance_shift:.6f} mm') + print(f'Z end (MSE): {mse_end:.5f}, n_refl: {final_count}') + print(f'net distance shift: {self.net_distance_shift:.6f} mm') + return mse_start, mse_end, self.net_distance_shift + def augment(expts, refls, d_min, d_max): """ Add pairwise 3D spot distances to the d-spacing histogram """ lab = flex.vec3_double() diff --git a/xfel/small_cell/small_cell.py b/xfel/small_cell/small_cell.py index cc218f2f730..cc38186ee41 100644 --- a/xfel/small_cell/small_cell.py +++ b/xfel/small_cell/small_cell.py @@ -1,6 +1,7 @@ from __future__ import absolute_import, division, print_function from six.moves import range from six.moves import zip +import copy #-*- Mode: Python; c-basic-offset: 2; indent-tabs-mode: nil; tab-width: 8 -*- # # LIBTBX_PRE_DISPATCHER_INCLUDE_SH export PHENIX_GUI_ENVIRONMENT=1 @@ -924,196 +925,223 @@ def small_cell_index_detail(experiments, reflections, horiz_phil, write_output = beam = imageset.get_beam() s0 = col(beam.get_s0()) - lattice_results = small_cell_index_lattice_detail(experiments, reflections, horiz_phil) - if not lattice_results: - return None - - max_clique_len, all_spots_len, ori, indexed = lattice_results - integrated_count = 0 - - if ori is not None and horiz_phil.small_cell.write_gnuplot_input: - write_cell(ori,beam,indexed,horiz_phil) - - indexed_hkls = flex.vec2_double() - indexed_intensities = flex.double() - indexed_sigmas = flex.double() - - if ori is not None: # ok to integrate - results = [] - buffers = [] - backgrounds = [] - indexed_hkls = flex.miller_index() - indexed_intensities = flex.double() - indexed_sigmas = flex.double() - mapped_predictions = flex.vec2_double() - mapped_panels = flex.size_t() - max_signal = flex.double() - xyzobs = flex.vec3_double() - xyzvar = flex.vec3_double() - shoeboxes = flex.shoebox() - s1 = flex.vec3_double() - bbox = flex.int6() - - raw_data = imageset[0] - if not isinstance(raw_data, tuple): - raw_data = (raw_data,) - rmsd = 0 - rmsd_n = 0 - for spot in indexed: - if spot.pred is None: continue - peakpix = [] - peakvals = [] - tmp = [] - is_bad = False - panel = detector[spot.pred_panel_id] - panel_raw_data = raw_data[spot.pred_panel_id] - for p in spot.peak_pixels: - #if is_bad_pixel(panel_raw_data,p): - # is_bad = True - # break - p = (p[0]+.5,p[1]+.5) - peakpix.append(p) - tmp.append(p) - peakvals.append(panel_raw_data[int(p[1]),int(p[0])]) - if is_bad: continue - - buffers.append(grow_by(peakpix,1)) - - tmp.extend(buffers[-1]) - backgrounds.append(grow_by(tmp,1)) - tmp.extend(backgrounds[-1]) - backgrounds[-1].extend(grow_by(tmp,1)) - - background = [] - bg_vals = [] - raw_bg_sum = 0 - for p in backgrounds[-1]: - try: - i = panel_raw_data[int(p[1]),int(p[0])] - except IndexError: - continue - if i is not None and i > 0: - background.append(p) - bg_vals.append(i) - raw_bg_sum += i - - ret = reject_background_outliers(background, bg_vals) - if ret is None: - print("Not enough background pixels to integrate spot %d"%spot.ID) - continue - background, bg_vals = ret - backgrounds[-1] = background - - bp_a,bp_b,bp_c = get_background_plane_parameters(bg_vals, background) - - intensity = 0 - bg_peak = 0 - for v,p in zip(peakvals,peakpix): - intensity += v - (bp_a*p[0] + bp_b*p[1] + bp_c) - bg_peak += bp_a*p[0] + bp_b*p[1] + bp_c + if horiz_phil.indexing.stills.reflection_subsampling.enable: + all_reflections = copy.deepcopy(reflections) + subsets = range( + horiz_phil.indexing.stills.reflection_subsampling.step_start, + horiz_phil.indexing.stills.reflection_subsampling.step_stop + - horiz_phil.indexing.stills.reflection_subsampling.step_size, + -horiz_phil.indexing.stills.reflection_subsampling.step_size, + ) + else: + subsets = [100] + for i_subset, pct in enumerate(subsets): + if pct != 100: + reflections = all_reflections.select( + flex.random_permutation(len(all_reflections)) + )[: int(len(all_reflections) * pct / 100)] - gain = panel.get_gain() - sigma = math.sqrt(gain * (intensity + bg_peak + ((len(peakvals)/len(bg_vals))**2) * raw_bg_sum)) - print("ID: %3d, ohkl: %s, ahkl: %s, I: %9.1f, sigI: %9.1f, RDiff: %9.6f"%( \ - spot.ID, spot.hkl.get_ohkl_str(), spot.hkl.get_ahkl_str(), intensity, sigma, - (sqr(ori.reciprocal_matrix())*spot.hkl.ohkl - spot.xyz).length())) + try: + lattice_results = small_cell_index_lattice_detail(experiments, reflections, horiz_phil) + if not lattice_results: + continue - max_sig = panel_raw_data[int(spot.spot_dict['xyzobs.px.value'][1]),int(spot.spot_dict['xyzobs.px.value'][0])] + max_clique_len, all_spots_len, ori, indexed = lattice_results + integrated_count = 0 + + if ori is not None and horiz_phil.small_cell.write_gnuplot_input: + write_cell(ori,beam,indexed,horiz_phil) + + indexed_hkls = flex.vec2_double() + indexed_intensities = flex.double() + indexed_sigmas = flex.double() + + if ori is not None: # ok to integrate + results = [] + buffers = [] + backgrounds = [] + indexed_hkls = flex.miller_index() + indexed_intensities = flex.double() + indexed_sigmas = flex.double() + mapped_predictions = flex.vec2_double() + mapped_panels = flex.size_t() + max_signal = flex.double() + xyzobs = flex.vec3_double() + xyzvar = flex.vec3_double() + shoeboxes = flex.shoebox() + s1 = flex.vec3_double() + bbox = flex.int6() + + raw_data = imageset[0] + if not isinstance(raw_data, tuple): + raw_data = (raw_data,) + rmsd = 0 + rmsd_n = 0 + for spot in indexed: + if spot.pred is None: continue + peakpix = [] + peakvals = [] + tmp = [] + is_bad = False + panel = detector[spot.pred_panel_id] + panel_raw_data = raw_data[spot.pred_panel_id] + for p in spot.peak_pixels: + #if is_bad_pixel(panel_raw_data,p): + # is_bad = True + # break + p = (p[0]+.5,p[1]+.5) + peakpix.append(p) + tmp.append(p) + peakvals.append(panel_raw_data[int(p[1]),int(p[0])]) + if is_bad: continue + + buffers.append(grow_by(peakpix,1)) + + tmp.extend(buffers[-1]) + backgrounds.append(grow_by(tmp,1)) + tmp.extend(backgrounds[-1]) + backgrounds[-1].extend(grow_by(tmp,1)) + + background = [] + bg_vals = [] + raw_bg_sum = 0 + for p in backgrounds[-1]: + try: + i = panel_raw_data[int(p[1]),int(p[0])] + except IndexError: + continue + if i is not None and i > 0: + background.append(p) + bg_vals.append(i) + raw_bg_sum += i + + ret = reject_background_outliers(background, bg_vals) + if ret is None: + print("Not enough background pixels to integrate spot %d"%spot.ID) + continue + background, bg_vals = ret + backgrounds[-1] = background + + bp_a,bp_b,bp_c = get_background_plane_parameters(bg_vals, background) + + intensity = 0 + bg_peak = 0 + for v,p in zip(peakvals,peakpix): + intensity += v - (bp_a*p[0] + bp_b*p[1] + bp_c) + bg_peak += bp_a*p[0] + bp_b*p[1] + bp_c + + gain = panel.get_gain() + sigma = math.sqrt(gain * (intensity + bg_peak + ((len(peakvals)/len(bg_vals))**2) * raw_bg_sum)) + + print("ID: %3d, ohkl: %s, ahkl: %s, I: %9.1f, sigI: %9.1f, RDiff: %9.6f"%( \ + spot.ID, spot.hkl.get_ohkl_str(), spot.hkl.get_ahkl_str(), intensity, sigma, + (sqr(ori.reciprocal_matrix())*spot.hkl.ohkl - spot.xyz).length())) + + max_sig = panel_raw_data[int(spot.spot_dict['xyzobs.px.value'][1]),int(spot.spot_dict['xyzobs.px.value'][0])] + + s = "Orig HKL: % 4d % 4d % 4d "%(spot.hkl.ohkl.elems) + s = s + "Asu HKL: % 4d % 4d % 4d "%(spot.hkl.ahkl.elems) + s = s + "I: % 10.1f sigI: % 8.1f I/sigI: % 8.1f "%(intensity, sigma, intensity/sigma) + s = s + "Size (pix): %3d Max pix val: %6d\n"%(len(spot.peak_pixels),max_sig) + results.append(s) + + if spot.pred is None: + mapped_predictions.append((spot.spot_dict['xyzobs.px.value'][0], spot.spot_dict['xyzobs.px.value'][1])) + mapped_panels.append(spot.spot_dict['panel']) + else: + mapped_predictions.append((spot.pred[0],spot.pred[1])) + mapped_panels.append(spot.pred_panel_id) + xyzobs.append(spot.spot_dict['xyzobs.px.value']) + xyzvar.append(spot.spot_dict['xyzobs.px.variance']) + shoeboxes.append(spot.spot_dict['shoebox']) + + indexed_hkls.append(spot.hkl.ohkl.elems) + indexed_intensities.append(intensity) + indexed_sigmas.append(sigma) + max_signal.append(max_sig) + s1.append(s0+spot.xyz) + bbox.append(spot.spot_dict['bbox']) + + if spot.pred is not None: + rmsd_n += 1 + rmsd += measure_distance(col((spot.spot_dict['xyzobs.px.value'][0],spot.spot_dict['xyzobs.px.value'][1])),col(spot.pred))**2 + + if len(results) >= horiz_phil.small_cell.min_spots_to_integrate: + # Uncomment to get a text version of the integration results + #f = open(os.path.splitext(os.path.basename(path))[0] + ".int","w") + #for line in results: + # f.write(line) + #f.close() + + if write_output: + info = dict( + xbeam = refined_bcx, + ybeam = refined_bcy, + distance = distance, + wavelength = wavelength, + pointgroup = horiz_phil.small_cell.spacegroup, + observations = [cctbx.miller.set(sym,indexed_hkls).array(indexed_intensities,indexed_sigmas)], + mapped_predictions = [mapped_predictions], + mapped_panels = [mapped_panels], + model_partialities = [None], + sa_parameters = [None], + max_signal = [max_signal], + current_orientation = [ori], + current_cb_op_to_primitive = [sgtbx.change_of_basis_op()], # identity. only support primitive lattices. + pixel_size = pixel_size, + ) + G = open("int-" + os.path.splitext(os.path.basename(path))[0] +".pickle","wb") + import pickle + pickle.dump(info,G,pickle.HIGHEST_PROTOCOL) + + crystal = ori_to_crystal(ori, horiz_phil.small_cell.spacegroup) + experiments = ExperimentListFactory.from_imageset_and_crystal(imageset, crystal) + if write_output: + experiments.as_file( + os.path.splitext(os.path.basename(path).strip())[0] + "_integrated.expt" + ) + + refls = flex.reflection_table() + refls['id'] = flex.int(len(indexed_hkls), 0) + refls['panel'] = mapped_panels + refls['intensity.sum.value'] = indexed_intensities + refls['intensity.sum.variance'] = indexed_sigmas**2 + refls['xyzobs.px.value'] = xyzobs + refls['xyzobs.px.variance'] = xyzvar + refls['miller_index'] = indexed_hkls + refls['xyzcal.px'] = flex.vec3_double(mapped_predictions.parts()[0], mapped_predictions.parts()[1], flex.double(len(mapped_predictions), 0)) + refls['shoebox'] = shoeboxes + refls['entering'] = flex.bool(len(refls), False) + refls['s1'] = s1 + refls['bbox'] = bbox + + refls.centroid_px_to_mm(experiments) + + refls.set_flags(flex.bool(len(refls), True), refls.flags.indexed) + if write_output: + refls.as_pickle(os.path.splitext(os.path.basename(path).strip())[0]+"_integrated.refl") + + print("cctbx.small_cell: integrated %d spots."%len(results), end=' ') + integrated_count = len(results) + else: + raise RuntimeError("cctbx.small_cell: not enough spots to integrate (%d)."%len(results)) - s = "Orig HKL: % 4d % 4d % 4d "%(spot.hkl.ohkl.elems) - s = s + "Asu HKL: % 4d % 4d % 4d "%(spot.hkl.ahkl.elems) - s = s + "I: % 10.1f sigI: % 8.1f I/sigI: % 8.1f "%(intensity, sigma, intensity/sigma) - s = s + "Size (pix): %3d Max pix val: %6d\n"%(len(spot.peak_pixels),max_sig) - results.append(s) + if rmsd_n > 0: + print(" RMSD: %f"%math.sqrt((1/rmsd_n)*rmsd)) + else: + print(" Cannot calculate RMSD. Not enough integrated spots or not enough clique spots near predictions.") - if spot.pred is None: - mapped_predictions.append((spot.spot_dict['xyzobs.px.value'][0], spot.spot_dict['xyzobs.px.value'][1])) - mapped_panels.append(spot.spot_dict['panel']) + print("IMAGE STATS %s: spots %5d, max clique: %5d, integrated %5d spots"%(path,all_spots_len,max_clique_len,integrated_count)) + except Exception: + if i_subset == len(subsets)-1: + raise else: - mapped_predictions.append((spot.pred[0],spot.pred[1])) - mapped_panels.append(spot.pred_panel_id) - xyzobs.append(spot.spot_dict['xyzobs.px.value']) - xyzvar.append(spot.spot_dict['xyzobs.px.variance']) - shoeboxes.append(spot.spot_dict['shoebox']) - - indexed_hkls.append(spot.hkl.ohkl.elems) - indexed_intensities.append(intensity) - indexed_sigmas.append(sigma) - max_signal.append(max_sig) - s1.append(s0+spot.xyz) - bbox.append(spot.spot_dict['bbox']) - - if spot.pred is not None: - rmsd_n += 1 - rmsd += measure_distance(col((spot.spot_dict['xyzobs.px.value'][0],spot.spot_dict['xyzobs.px.value'][1])),col(spot.pred))**2 - - if len(results) >= horiz_phil.small_cell.min_spots_to_integrate: - # Uncomment to get a text version of the integration results - #f = open(os.path.splitext(os.path.basename(path))[0] + ".int","w") - #for line in results: - # f.write(line) - #f.close() - - if write_output: - info = dict( - xbeam = refined_bcx, - ybeam = refined_bcy, - distance = distance, - wavelength = wavelength, - pointgroup = horiz_phil.small_cell.spacegroup, - observations = [cctbx.miller.set(sym,indexed_hkls).array(indexed_intensities,indexed_sigmas)], - mapped_predictions = [mapped_predictions], - mapped_panels = [mapped_panels], - model_partialities = [None], - sa_parameters = [None], - max_signal = [max_signal], - current_orientation = [ori], - current_cb_op_to_primitive = [sgtbx.change_of_basis_op()], # identity. only support primitive lattices. - pixel_size = pixel_size, - ) - G = open("int-" + os.path.splitext(os.path.basename(path))[0] +".pickle","wb") - import pickle - pickle.dump(info,G,pickle.HIGHEST_PROTOCOL) - - crystal = ori_to_crystal(ori, horiz_phil.small_cell.spacegroup) - experiments = ExperimentListFactory.from_imageset_and_crystal(imageset, crystal) - if write_output: - experiments.as_file( - os.path.splitext(os.path.basename(path).strip())[0] + "_integrated.expt" - ) - - refls = flex.reflection_table() - refls['id'] = flex.int(len(indexed_hkls), 0) - refls['panel'] = mapped_panels - refls['intensity.sum.value'] = indexed_intensities - refls['intensity.sum.variance'] = indexed_sigmas**2 - refls['xyzobs.px.value'] = xyzobs - refls['xyzobs.px.variance'] = xyzvar - refls['miller_index'] = indexed_hkls - refls['xyzcal.px'] = flex.vec3_double(mapped_predictions.parts()[0], mapped_predictions.parts()[1], flex.double(len(mapped_predictions), 0)) - refls['shoebox'] = shoeboxes - refls['entering'] = flex.bool(len(refls), False) - refls['s1'] = s1 - refls['bbox'] = bbox - - refls.centroid_px_to_mm(experiments) - - refls.set_flags(flex.bool(len(refls), True), refls.flags.indexed) - if write_output: - refls.as_pickle(os.path.splitext(os.path.basename(path).strip())[0]+"_integrated.refl") - - print("cctbx.small_cell: integrated %d spots."%len(results), end=' ') - integrated_count = len(results) - else: - raise RuntimeError("cctbx.small_cell: not enough spots to integrate (%d)."%len(results)) - - if rmsd_n > 0: - print(" RMSD: %f"%math.sqrt((1/rmsd_n)*rmsd)) + continue else: - print(" Cannot calculate RMSD. Not enough integrated spots or not enough clique spots near predictions.") + print('indexed on subset ', i_subset) + break - print("IMAGE STATS %s: spots %5d, max clique: %5d, integrated %5d spots"%(path,all_spots_len,max_clique_len,integrated_count)) return max_clique_len, experiments, refls def spots_rmsd(spots):